# Question
Title: How to process multiple VCF files present in a 2 folders and create afd plots
*This question has also been asked on Biostars*

I have 2 directories that contains multiple VCF files.

I want to write code in python that will process all vcf files and create an AFD (allele frequency difference) plot for each VCF file present in a folder

I have started with this code, which attempts to convert the VCF files into a pandas dataframe:

```
# Import Module
import os
import pandas as pd
# Folder Path
path1 = "C://Users//USER//Desktop//isiah/VCFs_1/"
path2 = "C://Users//USER//Desktop//isiah/VCFs_2/"
#os.chdir(path1)

# Read text File

def read_text_file(file_path1,file_path2):
    with open(file_path, 'r') as f:
        print(f.read())

# iterate through all file
for file in os.listdir():
    # Check whether file is in text format or not
    if file.endswith(".vcf"):
        file_path1 = f"{path1}\{file}"
        file_path2 = f"{path2}\{file}"
        print(file_path1,"\n\n",file_path2)

# call read text file function
#data = read_text_file(path1,path2)
#print(read_text_file(path1,path2))
Final= pd.DataFrame(read_text_file(path1,path2))
print(Final)

```

and its giving me output without calling function and also gives me empty dataframe. Also i am not sure if result is showing all vcf files record or not.

Kindly help me to modify this code. I need my all result in dataframe so i can make plot easily

Folder 1 contains 3 test vcfs whereas folder 2 contains 4 test vcfs:

This is test 1 vcf file:

```
##fileformat=VCFv4.2
##ALT=<ID=NON_REF,Description="Represents any possible alternative allele at this location">
##FILTER=<ID=FS,Description="FS>60.0">
##FILTER=<ID=HS,Description="HaplotypeScore<13.0">
##FILTER=<ID=LowQual,Description="Low quality">
##FILTER=<ID=MQ,Description="MQ<40.0">
##FILTER=<ID=MQRS,Description="MQRankSum < -12.5">
##FILTER=<ID=Mask,Description="Doesn't overlap a user-input mask">
##FILTER=<ID=PASS,Description="All filters passed">
##FILTER=<ID=QD,Description="QD<2.0">
##FILTER=<ID=RPRS,Description="ReadPosRankSum < -8.0">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=PGT,Number=1,Type=String,Description="Physical phasing haplotype information, describing how the alternate alleles are phased in relation to one another">
##FORMAT=<ID=PID,Number=1,Type=String,Description="Physical phasing ID information, where each unique ID within a given sample (but not across samples) connects records within a phasing group">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
##FORMAT=<ID=RGQ,Number=1,Type=Integer,Description="Unconditional reference genotype confidence, encoded as a phred quality -10*log10 p(genotype call is wrong)">
##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias.">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##INFO=<ID=BaseQRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt Vs. Ref base qualities">
##INFO=<ID=ClippingRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref number of hard clipped bases">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth; some reads may have been filtered">
##INFO=<ID=DS,Number=0,Type=Flag,Description="Were any of the samples downsampled?">
##INFO=<ID=ExcessHet,Number=1,Type=Float,Description="Phred-scaled p-value for exact test of excess heterozygosity">
##INFO=<ID=FS,Number=1,Type=Float,Description="Phred-scaled p-value using Fisher's exact test to detect strand bias">
##INFO=<ID=HaplotypeScore,Number=1,Type=Float,Description="Consistency of the site with at most two segregating haplotypes">
##INFO=<ID=InbreedingCoeff,Number=1,Type=Float,Description="Inbreeding coefficient as estimated from the genotype likelihoods per-sample when compared against the Hardy-Weinberg expectation">
##INFO=<ID=MLEAC,Number=A,Type=Integer,Description="Maximum likelihood expectation (MLE) for the allele counts (not necessarily the same as the AC), for each ALT allele, in the same order as listed">
##INFO=<ID=MLEAF,Number=A,Type=Float,Description="Maximum likelihood expectation (MLE) for the allele frequency (not necessarily the same as the AF), for each ALT allele, in the same order as listed">
##INFO=<ID=MQ,Number=1,Type=Float,Description="RMS Mapping Quality">
##INFO=<ID=MQRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities">
##INFO=<ID=QD,Number=1,Type=Float,Description="Variant Confidence/Quality by Depth">
##INFO=<ID=RAW_MQ,Number=1,Type=Float,Description="Raw data for RMS Mapping Quality">
##INFO=<ID=ReadPosRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias">
##INFO=<ID=SOR,Number=1,Type=Float,Description="Symmetric Odds Ratio of 2x2 contingency table to detect strand bias">
##reference=file:///Users/sian_bray/Dropbox/Yant/Cochlearia/02_Data/References/pseudohap_cochlearia_135M_3kb.fasta
##source=SelectVariants
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  ALO_006 ALO_007 ALO_013 ALO_017 BEA_002 BEA_004 BEA_010 BRI_002 BRI_005 BRI_006 BRI_009 ELI_001 ELI_002 ELI_003 ELI_004 FRE_013 JON_001 KVA_002 KVA_003 KVA_009 KVA_010 LAB_004 LAB_300 LAB_400 LAB_500 LNL_001 LNL_002 LNL_003 LNL_008 MEL_001 MEL_002 NEN_001 NEN_003 NEN_200 NEN_300 ROT_004 ROT_006 ROT_007 ROT_013 RUZ_001 SAL_003 SAL_004 SAL_006 SKF_002 SKF_003 SKF_005 SKF_009 SKI_004 SKI_005 SKI_008 SKN_001 SKN_002 SKN_005 SKN_008 SPU_006 SPU_008 SPU_009 SPU_010 TET_002 TET_004 TET_006 TET_008 TRO_001 TRO_003 TRO_005 TRO_009 VAG_003 VAG_004 VAG_007 VAG_102 VEG_003 VEG_004 WOL_002 WOL_006 WOL_009 WOL_010
237 50041   .   G   T   4099.26 PASS    AC=35;AF=0.133;AN=264;BaseQRankSum=0.00;ClippingRankSum=0.00;DP=1228;ExcessHet=4.3562;FS=0.644;InbreedingCoeff=-0.0953;MQ=57.84;MQRankSum=0.00;QD=25.94;ReadPosRankSum=-1.570e-01;SOR=0.768 GT:AD:DP:GQ:PL  0/0/0/0:11,0:11:14:0,14,33,66,425   0/0/0/0:11,1:12:0:0,0,0,27,369  0/0/0/0:13,3:16:0:0,0,0,0,355   0/0/0/0:41,0:41:50:0,50,120,241,1800    0/0/1/1:5,8:13:1:293,15,0,1,170 0/0/0/1:5,3:8:0:104,0,0,10,191  0/0/0/0:8,0:8:10:0,10,24,48,272 0/1/1/1:1,8:9:11:318,33,11,0,26 ./././.:0,0 1/1/1/1:0,2:2:2:87,12,6,2,0 0/0/0/0:5,2:7:0:0,0,0,0,120 0/0/0/0:16,2:18:0:0,0,0,21,515  0/0/0/0:8,1:9:0:0,0,0,13,270    0/0/0/0:15,3:18:0:0,0,0,0,459   0/0/0/0:10,1:11:0:0,0,0,21,363  0/0/0/0/0/0:15,0:15:12:0,12,26,45,71,116,589    0/0/0/0/0/0:7,2:9:0:0,0,0,0,0,0,151 ./././.:0,0 0/0/0/0:14,0:14:14:0,14,33,66,495   0/0/0/0:1,0:1:1:0,1,3,6,42  0/0/0/0:9,0:9:10:0,10,24,48,360 0/0:21,2:23:14:0,14,726 0/0:14,2:16:11:0,11,452 0/0:14,2:16:0:0,0,467   0/0:23,9:32:0:0,0,527   0/0/0/0:10,0:10:12:0,12,30,60,359   0/0/0/0:11,1:12:0:0,0,0,27,404  0/0/0/0:13,0:13:16:0,16,39,78,526   0/0/0/0:20,3:23:0:0,0,0,2,655   1/1/1/1:0,4:4:5:155,24,12,5,0   0/0/0/0:2,3:5:0:0,0,0,0,8   0/0:23,2:25:0:0,0,798   0/0:12,2:14:0:0,0,381   0/0:14,0:14:36:0,36,540 0/0:24,0:24:63:0,63,945 0/0/0/0:8,1:9:0:0,0,0,12,245    0/0/0/0:27,0:27:32:0,32,78,157,1170 0/0/0/0:18,5:23:0:0,0,0,0,474   0/0/0/0:22,1:23:15:0,15,50,114,736  0/0:10,0:10:27:0,27,405 0/0:8,0:8:24:0,24,322   0/0:13,0:13:36:0,36,540 0/0:4,0:4:12:0,12,155   0/0/0/0/0/0:5,2:7:0:0,0,0,0,0,0,125 0/0/0/0/0/0:19,1:20:0:0,0,0,18,50,106,731   0/0/0/0/0/0:28,3:31:0:0,0,0,0,1,82,978  0/0/0/0/0/0:53,10:63:0:0,0,0,0,0,0,1843 0/0/0/0:1,1:2:0:0,0,0,0,4   1/1/1/1:0,22:22:27:917,132,66,27,0  ./././.:0,0 0/0/0/0:15,0:15:17:0,17,42,84,630   0/0/0/0:14,2:16:0:0,0,0,27,421  0/0/0/0:17,0:17:21:0,21,51,102,670  0/0/0/0:17,7:24:0:0,0,0,0,323   0/0:11,2:13:0:0,0,380

```

# Answer
> 2 votes
*Answer from @user324810, converted from comment:*

There is an alternative than doing it with Python unless it is absolutely necessary. If you want to merge multiple samples into one, use bcftools merge. If you want simply to concatenate multiple VCFs, use bcftools concatenate. Once you obtain the output file, you can load the dataframe and plot it as you wish.

---
Tags: python, vcf
---


# Question
Title: How should I imagine the pre-indexing of suffix arrays in the STAR RNA Seq-aligner?
I'm using STAR for an internship (fusion-genes in ALL cancer) and I'd like to understand what exactly it is doing. The people at the internship don't know either how it works, given their biology backgrounds. They told me "its enough to know that it works". That was not an satisfactory answer for me which is why I am fighting my way trough the suplementary materials trying to wrap my head around the algo. My background is bio as well but I've spent reading a book about algos and algos in BI this year.

I am trying to understand the section *Pre-indexing of suffix arrays* from supplementary materials of the *Dobin et al.* STAR RNA seq-aligner Publication.

> While suffix array search is theoretically fast owing to its binary nature, in practice it may suffer from non-locality resulting in persistent cache misses which deteriorate the performance. To alleviate this problem we developed a pre-indexing strategy. After the SA is generated, we find the locations of all possible $L-mers$ in the SA, $L\<=L\_{max}$, where L\_{max} is user defined and is typically 12-15. Since the nucleotide alphabet contains only four letters, there are NL=22L different L-mers for which the SA locations have to be stored.

The authors also state:

> Suffix Array (SA) of the whole genome is utilized to find the Maximum Mappable Prefixes (MMP).

1. I am confused by the use of the term $L-mer$. When googleing I only keep finding $k-mers$ in combination with suffix arrays. I did however find this GitHub page mentioning $L-mers$ in the following way:

> This program computes either the suffix array or the L-mer array of a string. The suffix array gives the lexicographical order of all the suffixes of a string. The L-mer array gives the lexicographical order of all the L-mers in the string (sorts the suffixes with respect to only the first L characters).

→ Can I assume that $L-mer$ was used instead of $k-mer$ to stress the lexicographical ordering? But as far as I understood *all* suffix arrays contain the suffixes' index based on its lexicographic ordering as illustrated by wikipedia's page on suffix arrays:

2. Do I have to imagine the entire genome being a single string that is converted into one single suffix array? Could you provide sources fit for a beginner with biology background and very little idea about algorithms?

# Answer
You can take a look at slide 40 of the lectures slides I use when I teach suffix arrays (reproduced here)

An L-mer *is* a K-mer — STAR builds an explicit lookup table from all length 14 (by default) strings that maps the length 14 string to interval of the suffix array that contains **all suffixes that begin with this length 14 prefix**. This means that, to find a MMP, instead of having to do a binary search on the entire suffix array, you only have to do a binary search on the sub-interval that begins with the length 14 prefix that your query string starts with. If all substrings suffixes began with length 14 prefixes with equal frequency, then this would reduce the size of the interval you have to look at by a factor of $268435456 = 4^{14}$. Of course, in biological sequence, all prefixes are not equally likely, so the speedup is more variable in practice. Nonetheless, it makes lookup in the suffix array much, much faster.

Of course, in addition to the core algorithm, STAR contains many heuristics not all of which are likely described in the paper (especially since the tool has been and continues to be under active development and improvement since the paper's publication). Nonetheless, this is the idea behind the suffix array pre-indexing. Since STAR still has the whole suffix array, if there is no length 14 exact prefix match, it is still possible to find a smaller match using a more standard binary search in the suffix array. However, I do not know (and am not sure where it is documented) exactly when/if STAR will fall back to a more naive binary search if a prefix L-mer is not found. If you are truly curious about this, I'd recommend asking Alex Dobin (the STAR author and maintainer) over on the STAR GitHub repository. He is generally very responsive to user queries there.
> 2 votes

---
Tags: algorithms, star, indexing
---


# Question
Title: How to consume non-CRAN packages in custom package?
In a custom package which imports a package from Bioconductor, what is the right DESCRIPTION file?

I added in my DESCRIPTION a line with Remote like so

```
Remotes:
    bioc::3.14/dada2

```

as suggested by https://cran.r-project.org/web/packages/remotes/index.html vignette, but it still failed the R `check()` with error "Packages required not available". I am very confused about what is the right way to solve this problem, since multiple sources propose a different solution:

What is the appropriate way to declare dependencies on packages on Bioconductor? The remotes support helped me satisfy a dependency on a package from Github, but not for bioconductor

# Answer
> 2 votes
As I answered on the linked question, the right way to declare packages from repositories recognized by R is to just declare normally on the Imports, Depends, Enhances and LinkingTo.

Bioconductor is an accepted repository so it doesn't need to be added on Additional\_repository field. BiocViews is a required field from Bioconductor. As CRAN accepts other fields on the DESCRIPTION it could be used by CRAN packages.

You don't mention it but if you found this problem and it wasn't on your computer were you are developing the package it means you are using a continuous integration test. The problem with CI like Github actions, travis, appveyor or other is that dependencies must be installed before (as they should on your computer).

To install these dependencies before that depends on the CI you use and the strategy used there (install from ppa, use RSPM, use `bspm::install`, ...) . However, I recommend to do it via `BiocManager::install` because this way you can install all packages from both CRAN and Bioconductor and their dependencies.

Even if the package is on CRAN, github or just Bioconductor with `BiocManager::install` you'll be able to install all the dependencies and continue to check your package normally on the CI.

---

Additional\_repositories is for other non standard repositories like a personal one or other a specific of a project. It is used by only ~100 packages of the 18765 packages on CRAN. But there are many more dependencies from CRAN to Bioconductor (and vice versa) .

---
Tags: bioconductor
---


title: Polygenic Risk Scores

text:
# Question
Title: Polygenic Risk Scores
I have run a GWAS and used PRSice (https://www.prsice.info/) to calculate polygenic risk scores on the data collected from the GWAS. I need to check whether alcohol dependence can predict gambling disorder, through the PRS. However I am having trouble with how exactly to do that. I have a GWAS on alcohol disorder from the one population and a GWAS on gambling and alcohol disorder from the other population, then I calculated PRS for the three GWAS datasets, but how do I use the one to predict in the other? Through the literature it was suggested to use linear mixed models, which I have done but I still feel as if I am doing it wrong as I only did the model on the one dataset and not the two combined.

# Answer
> 1 votes
Before doing anything else, you should split your dataset into training and test groups, and \[this is important\] don't touch the test dataset until you're fully confident that your approach is correct: no post-test modification.

A split that's considered acceptable in Computer Science is a 10%/90% split, i.e. use 10% of the data for training / algorithm generation and the remaining 90% for testing. This will give you good confidence that the algorithm works.

I personally think it's a good idea to do population sub-sampling (i.e. further subsample that 10%) to determine an informative marker set prior to any risk calculations. What you should be looking for are markers that are consistently linked with the trait of interest (i.e. in your case alcohol dependence) in multiple sub-sampled populations. I go into a lot of detail about how to do this here \[bearing in mind it's an old pre-print, and doesn't do the 10/90 split I've recommended above\].

After you're happy with the marker set, calculate the risk scores *only using that set* (hopefully PRSice can do this), then carry that risk score \[only\] through to test against the dependent trait of interest in the training dataset (i.e. in your case gambling disorder). If your dependent trait is a binary trait, this could be done using ROC analysis (as in my preprint) to calculate area under the curve. Or use linear mixed models, whatever works for you.

Make any changes you think you need to improve the test result (e.g. changing the composition of the marker set).

After convincing yourself that the approach is correct, do the same risk score calculation (using the same set of markers) on the test set, find out how well it predicts the dependent trait, and report on the results. Don't change anything.

---
Tags: python, gwas, polygenic-risk-score
---


title: Organelle genome genes

text:
# Question
Title: Organelle genome genes
How can I get all proteins coded by mitochondrial genome or plastid genome from UniProt. I see UniProt has multiple filtering and search parameters but I could not distinguish between proteins/genes coded by nuclear DNA and mitochondrial or plastid genomes.

I specifically want UniProt (rather than e.g. Organelle Genome Resources) because usually annotation of proteins uses UniProt ids.

Any suggestions?

# Answer
> 2 votes
You can use the UniProt Advanced search for this.

If you look at this help page https://www.uniprot.org/help/advanced\_search, in particular the first example screenshot, you can type "encoded" where the screenshot has "structure". Then you can choose your organelle. For Plastid, for example, this will translate into the query syntax organelle:plastid :

https://www.uniprot.org/uniprot/?query=organelle%3APlastid&sort=score

Several other organelles are available, of course, e.g. Mitochondrion :

https://www.uniprot.org/uniprot/?query=organelle%3AMitochondrion&sort=score

Please don't hesitate to contact the UniProt helpdesk. We can unfortunately not guarantee to monitor all available forums as closely as we would like to, and can answer much more promptly when contacted directly.

---
Tags: proteins, genome, uniprot
---


# Question
Title: It is possible to have a UniProt protein entry that does not have any genomic coordinates, because lacks of an Ensembl link?
I have a list of UniProt codes and I am trying to fetch the genomic coordinates through the API from programmatic access to UniProtKB (https://www.ebi.ac.uk/proteins/api/doc/).

Moreover, I read that UniProt in collaboration with Ensembl have mapped the protein sequences on the GRCh38 genome (doi: 10.1093/nar/gkw1099)

However, it is true that the number of protein codes in UniProt is updated more frequently than the GRCh38 assembly.

Considering only the protein genomic coordinates data: It is possible that a new submission in UniProt (in TrEMBL for instance) does not have any correspondence with Ensembl, due to the possibility that the genomic coordinates of a new protein do not match any gene model stored in GRCh38?

If not: It is possible to say that all the protein genomic coordinates in UniProt have been mapped on the GRCh38 (in connection with Ensembl) and it is impossible that a new submission in UniProt lacks of this data?

Thanks

# Answer
If you look at this query https://www.uniprot.org/uniprot/?query=NOT%20database:(type:ensembl)&fil=organism:%22Homo%20sapiens%20(Human)%20\[9606\]%22&sort=score

and restrict it to the reference proteome (by using the link on the left hand side): https://www.uniprot.org/uniprot/?query=NOT%20database:(type:ensembl)&fil=proteome%3AUP000005640+AND+organism%3A%22Homo+sapiens+%28Human%29+%5B9606%5D%22&sort=score

you will see that there are "only" 1100 entries, and that they are all in UniProtKB/Swiss-Prot, but none in UniProtKB/TrEMBL, the unreviewed section.

The manual review of sequences for human and other proteomes is one of our most important and crucial tasks but as you can imagine this is an ongoing process. First because it represents an immense work, second because knowledge is still evolving, the human reference genome being for instance constantly improved by the Genome Reference Consortium. For the human proteome, there are still a small but non negligible number of UniProtKB/Swiss-Prot entries for which the sequence is not identical or is not represented at all by Ensembl sequences. We are constantly collaborating with Ensembl and other genome annotation resources trying to agree on identical sequences to fill the gaps. UniProtKB human entries are automatically mapped to the latest Ensembl release but you have to take into account that UniProtKB and Ensembl releases are not synchronized.

Starting with release 2022\_01, we will have cross-references from relevant UniProtKB/Swiss-Prot entries to MANE-Select (Matched Annotation from NCBI and EMBL-EBI, see https://www.ensembl.org/info/genome/genebuild/mane.html and https://www.ncbi.nlm.nih.gov/refseq/MANE/)
> 1 votes

---
Tags: ensembl, uniprot
---


# Question
Title: connect same ids with values in id rows (result of effectiveELD)
I can't cope with reworking the result. If id is the same I want to stick those items.

```
with open(x,'r') as rf:
    l = fr.readlines()
    previous_line_id = ''
    corrected_lines = []

    for i in l:
       i_split = i.split('\t')
       id_ = i_split[0]
       if id_ == previous_line
           xxxxxxxxxx
       else:
           corrected_lines .append(i)

```

result:

```
    WP_XXXXXXXX2.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5
    WP_192XXXXXX.1  PF13489 Methyltransferase domain    4
    XXXXXXX45.1 PF02310 B12 binding domain  12  <-----------  DOUBLE IDS
    XXXXXXX45.1 PF04055 Radical SAM superfamily 5 <-----------
    WP_XXXXXXX29.1  PF13489 Methyltransferase domain    4
    WP_0XXXXXXXX.1  PF13489 Methyltransferase domain    4
    WP_XXXXXXX90.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5  <--------------   DOUBLE ID
    WP_XXXXXXX90.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  7  <--------------
    WP_0XXXXXX62.1  PF13489 Methyltransferase domain    4

```

result I need:

```
WP_XXXXXXXX2.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5
WP_192XXXXXX.1  PF13489 Methyltransferase domain    4
XXXXXXX45.1 PF02310 B12 binding domain  12; PF04055 Radical SAM superfamily 5 <----------- STICKED SAME ID,OTHER VALUES
WP_XXXXXXX29.1  PF13489 Methyltransferase domain    4
WP_0XXXXXXXX.1  PF13489 Methyltransferase domain    4
WP_XXXXXXX90.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5; PF06325  Ribosomal protein L11 methyltransferase (PrmA)  7  <--------- STICKED STICKED SAME ID,OTHER VALUES
WP_0XXXXXX62.1  PF13489 Methyltransferase domain    4

```

# Answer
Assuming that the file `test.txt` is tab-separated values with these fields:

```
WP_XXXXXXXX2.1  PF06325     Ribosomal protein L11 methyltransferase (PrmA)      5
WP_192XXXXXX.1  PF13489     Methyltransferase domain                            4
XXXXXXX45.1     PF02310     B12 binding domain                                  12
XXXXXXX45.1     PF04055     Radical SAM superfamily                             5
WP_XXXXXXX29.1  PF13489     Methyltransferase domain                            4
WP_0XXXXXXXX.1  PF13489     Methyltransferase domain                            4
WP_XXXXXXX90.1  PF06325     Ribosomal protein L11 methyltransferase (PrmA)      5
WP_XXXXXXX90.1  PF06325     Ribosomal protein L11 methyltransferase (PrmA)      7
WP_0XXXXXX62.1  PF13489     Methyltransferase domain                            4

```

You can proceed first by reading your file:

```
with open("test.txt",'r') as rf:
    l = rf.readlines() # all the lines are stored in variable 'l'

```

You should create a dictionary with your ID as keys and the remaining line as values:

```
# create empty dictionary
dict_id = {}

for i in l:                                                               # for each line
    key_id = i.split("\t")[0]                                             # get the id
    if key_id not in dict_id.keys():                                      # if id as key is not found in dictionary
        value = "\t".join(i.split("\t")[1:]).rstrip()                     # parse the value
        dict_id[key_id] = [value]                                         # add the id as key to dictionary with its value
    else:                                                                 # else if id as key is already in dictionary
        dict_id[key_id].append("\t".join(i.split("\t")[1:]).rstrip())     # append to id as key

```

And now, you can print the result:

```
for key, value in dict_id.items():         # for each key, values in dictionary
    print(key + ": " + "; ".join(value))   # print the key and their values separated by '; '

```

And you should obtain:

```
WP_XXXXXXXX2.1: PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5
WP_192XXXXXX.1: PF13489 Methyltransferase domain        4
XXXXXXX45.1: PF02310    B12 binding domain      12; PF04055     Radical SAM superfamily 5
WP_XXXXXXX29.1: PF13489 Methyltransferase domain        4
WP_0XXXXXXXX.1: PF13489 Methyltransferase domain        4
WP_XXXXXXX90.1: PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5; PF06325      Ribosomal protein L11 methyltransferase (PrmA)  7
WP_0XXXXXX62.1: PF13489 Methyltransferase domain        4

```

The final script:

```
with open("test.txt",'r') as rf:
    l = rf.readlines()

dict_id = {}
for i in l:
    key_id = i.split("\t")[0]
    if key_id not in dict_id.keys():
        value = "\t".join(i.split("\t")[1:]).rstrip()
        dict_id[key_id] = [value]
    else:
        dict_id[key_id].append("\t".join(i.split("\t")[1:]).rstrip())

for key, value in dict_id.items():
    print(key + ": " + "; ".join(value))

```
> 2 votes

---
Tags: python
---


title: Topmed File Requirements

text:
# Question
Title: Topmed File Requirements
I am trying to run my VCF's through Topmed Imputation (https://imputation.biodatacatalyst.nhlbi.nih.gov/#!pages/home)

I am trying to make sure that I have the correct requirements for my file to be able to run through the software. These requirements are:

Create a separate vcf.gz file for each chromosome.

Variations must be sorted by genomic position.

GRCh37 or GRCh38 coordinates are required.

I have managed to split my VCF into each chromosome but I am unsure how to check which coordinates are being used and if my file is sorted by genomic position. Are there unix commands I can use to check, or do I just look in the file, Ive looked and not noticed anything.

# Answer
There are tools that allows you to sort your VCF file by genomic position like bcftools sort or vcf-sort. But if you don't want to install these software, you could always do with bash tools for one vcf:

```
# get the header of your VCF and write to file 'sorted.vcf'
# if your file has the extension .gz, use zgrep instead
grep "^#" your.vcf > sorted.vcf

# get the rest of the lines not starting with '#', sort them first by chromosome, then position and add them to file 'sorted.vcf'
grep -v "^#" your.vcf | sort -V -k1,1 -k2,2n >> sorted.vcf

```

But since you have one file for each chromosome and assuming that your files look like `sample_chrZ.vcf.gz`, where `Z` is the chromosome number, you will have to do a loop:

```
mkdir sorted_vcf
for chr in *.vcf.gz; do
   filename="${chr%.*}"
   zgrep "^#" "$chr" > sorted_vcf/"$filename"
   zgrep -v "^#" "$chr" | sort -V -k1,1 -k2,2n  >> sorted_vcf/"$filename"
   bgzip sorted_vcf/"$filename"
done

```

And if you want to know whether your VCF has been aligned against GRCh38 or GRCh37, check the header in the VCF file with a line starting with `##reference=`. A command-line like `grep "##reference=" your.vcf` should give you the result. However, if you don't have the reference mentioned in your VCF file, this is going to be harder to guess. This post should help you figure it out.
> 3 votes

---
Tags: vcf, imputation
---


title: Extracting a column from a vcf file

text:
# Question
Title: Extracting a column from a vcf file
I have a zipped vcf file of dbSNP hg38 version

No space left on my device to unzip that

I want to extract a column from that

I have tried this

```
cut -f 3 00-All.vcf.gz | sed 's/[\t]/,/g' > cols.csv

```

But that needs my vcf to be unzipped

Do you know any solution please?

# Answer
Either the `zcat` command or `gunzip -c` will decompress the file in a streaming fashion to the standard output without writing the decompressed data to disk. Pipe the output to `cut`, `sed`, etc.

Try

```
zcat 00-All.vcf.gz | cut -f 3 | sed 's/[\t]/,/g' > cols.csv

```

or

```
gunzip -c 00-All.vcf.gz | cut -f 3 | sed 's/[\t]/,/g' > cols.csv

```

But if you just want a single 'column' (I presume you mean field), then you are better off just doing `bcftools query`.
> 3 votes

---
Tags: vcf, snp, wes
---


# Question
Title: How to execute several processes sequentially, without channels?
Using DSL2.

I have processes A, B, C which I want to run in sequence - B waits for A to finish, C waits for B to finish.

They are not communicating via channels so I can't do A.out.collect().

I have tried making dummy channels but this is problematic and a hack. The onComplete subscription doesn't seem to let me trigger another process.

What is the correct way to run a sequence of processes which are not communicating or triggered by channels?

P.S. I am aware that this is not really the intended use case of Nextflow (which is oriented around channel communication/triggers... but... I just need to get this done.)

# Answer
> 1 votes
I would guess it is not possible without channels unless you are using `beforeScript` and `afterScript`. I came up with this solution:

```
#!/usr/bin/env nextflow

nextflow.enable.dsl=2

workflow {
  // C(B(A("")))

  next = ""
  next = A(next)
  next = B(next)
  next = C(next)
}

process A{
  echo true
  input: val _
  output: val ""

  """
  echo $task.process
  """
}

process B{
  echo true
  input: val _
  output: val ""

  """
  echo $task.process
  """
}

process C{
  echo true
  input: val _
  output: val ""

  """
  echo $task.process
  """
}

```

By using workflows you can change the execution order. But it is still a hack and I would prefer something like a `next:` keyword to specify the next process.

---

And for those not interested in DSL2:

```
process A {
  echo true

  output:
  val "" into B

  """
  echo "$task.process"
  """
}

process B {
  echo true

  input:
  val _ from B

  output:
  val "" into C

  """
  echo "$task.process"
  """
}

process C {
  echo true

  input:
  val _ from C

  """
  echo "$task.process"
  """
}

```

Also, I am quite curious to see other solutions.

---
Tags: nextflow, workflow-management
---


# Question
Title: Assign multiple taxids to a sequence when constructing a local BLAST database
I recently had a script fail due to poor handling of BLAST output. The BLAST `-outfmt staxids` field usually returns a single taxid, but occasionally it returns two or more taxids separated by a semicolon, such as `556514;701533`. Fixing the script to handle this should be fairly straightforward. But the original failure happened when I was searching a local copy of the NCBI nt database, and I'd like to create a test case with a small query and small database that reproduces this issue more rapidly.

I am familiar with the `-taxid` and `-taxid_map` flags of the `makeblastdb` command, but I have been unable to reproduce this behavior with BLAST; that is, make BLAST return multiple taxids for a sequence.

* If I provide `556514;701533` or `556514,701533` as a sequence's taxid with `-taxid` or `-taxid_map`, the `makeblastdb` command fails citing an integer conversion error.
* If I provide `556514 701533` as the taxid in the `taxid_map` file, the database will be indexed successfully but BLAST will return the invalid value `0` for `staxids`.
* If I duplicate each entry in the `-taxid_map` file so that there is an entry for each taxid (see below), only a single taxid is reported by BLAST.
```
seq1 556514
seq1 701533
seq2 556514
seq2 701533
seq3 556514
seq3 701533
...

```

When constructing a custom local BLAST database, how do I associated a database sequence with more than one taxid?

---

**UPDATE**: The offending entry (or one such entry) appears to be a duplicate or merged entry in the nt database. When I use `blastdbcmd` to retrieve the subject sequence by accession, nt returns a sequence with two deflines on single line, like so.

```
>accession1 plain text description >accession2 other description
ACGT...

```

Using either of the two accessions, `blastdbcmd` will return the exact same entry, and the two staxids shown above correspond to the two distinct accessions.

**So part of the issue seems to be how to construct a BLAST database such that multiple accessions can point to a single entry/subject.**

# Answer
> 4 votes
I inquired with the NCBI BLAST help desk today, and within an hour they provided a very helpful response.

To summarize, when redundant records in the database are collapsed, the deflines (accessions + descriptions) of the original entries are placed on a single line separated by an invisible control character (Ctrl-A or `\x01`). The `makeblastdb` program will recognize this and associate multiple accessions (from the combined defline) and taxids (from the taxid\_map file) to a single entry. If you are retrieving entries from e.g. nt or nr with `blastdbcmd`, you can use the `-ctrl_a` flag to format the data with these invisible control characters—by default the `>` character is used.

---

**tl;dr** I can reproduce the behavior described in the original question by:

1. creating a FASTA file containing an entry with multiple accessions separated by the invisible `\x01` character
2. create a taxid\_map file mapping the accessions to different taxids
3. format the FASTA file with `makeblastdb -taxid_map ...`
4. perform a BLAST search with a query that matches the entry in question; `-outfmt "6 staxids"` will return a string with multiple taxids.

---
Tags: blast, taxonomy, makeblastdb
---


title: if, else statment gives back else value

text:
# Question
Title: if, else statment gives back else value
I have a list of protein sequence ids (names). I want to sort the lines of a table, checking if the id is equal to the first column of the table. I split the table on lines and search for the id via regex:

code:

```
def map_line_to_ids(id):
    id= id.strip('\n')
    for s in lines:
        if re.search(id,s):
            return s #<------------ why always return else not if(s) when takes empty 'def'
        else:
            return id + '\n'

names = ['XXXXXXX45.1', 'YP_343434.1', 'WP_XXXXXXXX2.1', 'WP_XXXXXXX90.1',
         'WP_0XXXXXX62.1', 'XP_222222222.2', 'WP_192XXXXXX.1',
         'WP_XXXXXXX29.1', 'WP_0XXXXXXXX.1']

sorted_lines = list(map(map_line_to_ids, names))

```

names (order):

```
XXXXXXX45.1
YP_343434.1
WP_XXXXXXXX2.1
WP_XXXXXXX90.1
WP_0XXXXXX62.1
XP_222222222.2
WP_192XXXXXX.1
WP_XXXXXXX29.1
WP_0XXXXXXXX.1

```

Input table:

```
WP_XXXXXXXX2.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5
WP_192XXXXXX.1  PF13489 Methyltransferase domain    4
XXXXXXX45.1 PF02310 B12 binding domain  12; PF04055 Radical SAM superfamily 5
WP_XXXXXXX29.1  PF13489 Methyltransferase domain    4
WP_0XXXXXXXX.1  PF13489 Methyltransferase domain    4
WP_XXXXXXX90.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5; PF06325  Ribosomal protein L11 methyltransferase (PrmA)  7
WP_0XXXXXX62.1  PF13489 Methyltransferase domain    4

```

Expected output:

```
XXXXXXX45.1 PF02310 B12 binding domain  12; PF04055 Radical SAM superfamily 5
YP_343434.1
WP_XXXXXXXX2.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5
WP_XXXXXXX90.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5; PF06325  Ribosomal protein L11 methyltransferase (PrmA)  7
XP_222222222.2
WP_0XXXXXX62.1  PF13489 Methyltransferase domain    4
WP_192XXXXXX.1  PF13489 Methyltransferase domain    4
WP_XXXXXXX29.1  PF13489 Methyltransferase domain    4
WP_0XXXXXXXX.1  PF13489 Methyltransferase domain    4

```

# Answer
> 3 votes
Following on your last post, your initial file `test.txt` looks like that:

```
WP_XXXXXXXX2.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5
WP_192XXXXXX.1  PF13489 Methyltransferase domain    4
XXXXXXX45.1 PF02310 B12 binding domain  12
XXXXXXX45.1 PF04055 Radical SAM superfamily 5
WP_XXXXXXX29.1  PF13489 Methyltransferase domain    4
WP_0XXXXXXXX.1  PF13489 Methyltransferase domain    4
WP_XXXXXXX90.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5
WP_XXXXXXX90.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  7
WP_0XXXXXX62.1  PF13489 Methyltransferase domain    4

```

And you have another file `order.txt` that contains your ordered IDs:

```
XXXXXXX45.1
YP_343434.1
WP_XXXXXXXX2.1
WP_XXXXXXX90.1
WP_0XXXXXX62.1
XP_222222222.2
WP_192XXXXXX.1
WP_XXXXXXX29.1
WP_0XXXXXXXX.1

```

Since you want to merge by the first field in `test.txt` and sort them according to `order.txt`, you have to read first both of the files:

```
# read 'test.txt'
with open("test.txt",'r') as rf:
    l = rf.readlines()

# read sorted IDs 'order.txt'
with open("order.txt",'r') as rf:
    l_order = rf.readlines()

# remove newlines from each element in list
l_order = [i.strip() for i in l_order]

```

My previous answers describe how we obtain the merged fields in the file, so I will not go over it again. Since we need to sort the dictionary according to the list, we simply have to do:

```
# for each element in ordered list
for i in l_order:
    if i in dict_id.keys():                        # if element is in my dictionary's keys
        print(i + "\t" + "; ".join(dict_id[i]))    # print the key and the value
    else:
        print(i)                                   # else prints the element

```

This should print:

```
XXXXXXX45.1     PF02310 B12 binding domain      12; PF04055     Radical SAM superfamily 5
YP_343434.1
WP_XXXXXXXX2.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5
WP_XXXXXXX90.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5; PF06325      Ribosomal protein L11 methyltransferase (PrmA)  7
WP_0XXXXXX62.1  PF13489 Methyltransferase domain        4
XP_222222222.2
WP_192XXXXXX.1  PF13489 Methyltransferase domain        4
WP_XXXXXXX29.1  PF13489 Methyltransferase domain        4
WP_0XXXXXXXX.1  PF13489 Methyltransferase domain        4

```

Full code:

```
with open("test.txt",'r') as rf:
    l = rf.readlines()

with open("order.txt",'r') as rf:
    l_order = rf.readlines()

l_order = [i.strip() for i in l_order]

dict_id = {}
for i in l:
    key_id = i.split("\t")[0]
    if key_id not in dict_id.keys():
        value = "\t".join(i.split("\t")[1:]).rstrip()
        dict_id[key_id] = [value]
    else:
        dict_id[key_id].append("\t".join(i.split("\t")[1:]).rstrip())

for i in l_order:
    if i in dict_id.keys():
        print(i + "\t" + "; ".join(dict_id[i]))
    else:
        print(i)

```

PS: please for further pure programming questions, do post your questions in StackOverflow instead of Bioinformatics.

---
Tags: python
---


# Question
Title: Integrating scRNA-seq data using raw data
I am trying to reporduce the the following article Chen, M. B., Jiang, X., Quake, S. R., & Südhof, T. C. (2020). Persistent transcriptional programmes are associated with remote memory. Nature, 587(7834), 437-442.. In the methods section **Bioinformatics and data analysis - Mapping to the genome** they say "Next, the ‘canonical correlation analysis’ function from the Seurat package27 was used to align raw data from multiple experiments. " But they had previously aligned the data and generated gene counts. So by this do they mean anchor? I'm confused why would some one run CCA on the raw unaligned data. They only provided gene count data not raw data, but they use the CCA to restrict the gene's used from ~23000 down to just ~2500 of the more positively and negatively correlated (to the resulting CCA generated principle components) without providing the exact details of which genes they included or excluded before running their analysis.

Can I just do the CCA on the count data and hope for the same results? Have I misunderstood their filtering and alignment?

**EDIT** For clarity for future readers, my question boiled down to the the term alignment. I was confused because I thought alignment specifically referred to aligning sequenced fragments to a reference genome, here alignment basically refers to correcting for batch effects in scRNA-seq.

# Answer
> 1 votes
I believe when you say alignment, you mean aligning reads to a genome (sometimes to transcriptome) and count these to get count matrices. In the aforementioned paper, however, what is meant is "bringing different data sets to a level where they can be compared/integrated/...". Basically scRNA-seq data are heavily prone to batch effetcs and if you try to integrate different data sets (let's say data generated with 10x 3' and 10x 5' chemistries), there is a good chance is most of the signal would come from the technical differences rather than biology, for example resulting in clustering by scRNA-seq technology and not cell type/stat, hence there are many batch-effect-removal/data-integration/data-alignment tools for scRNA-seq.

Seurat has a nice CCA tutorial, you would not need the raw data, count data would be fine but you will need to be able to specify what data are coming from what source, the authors must have provided this information.

When it comes to what an "anchor" is, you would need to read about the CCA technique but in a rather naïve way, you can think of these as commonalities between data sets or some kind of a reference and are used to "align" these different data sets.

Moreover, reducing the number of genes from tens of thousands to thousands is just fine and is not really specific to CCA: not all genes are expected to be expressed at all times in all cell types and hence are not that informative and also scRNA-seq is not able to pick up lowly-expressed genes, causing these to show up with 0 counts in the count matrix, so again not that informative/interesting for the downstream analysis steps.

---
Tags: scrnaseq, seurat
---


# Question
Title: Filtering a vcf with a text file of SNP rsIDs
I have a vcf file containing an ID field with a snps and a list of human SNPs from dbSNP.I want to remove whatever in dbSNP from vcf file. The text file is like this:

```
rs978760828
rs1008829651
rs1052373574
rs1326880612
rs768019142
rs892501864
rs1010989343
rs1489251879
rs1022805358
rs1322538365

```

How to remove the common snps from dbSNP and my vcf?

# Answer
> 1 votes
I'd advise against using `grep`, `sed`, `awk` for processing vcfs unless you really know what you're doing - there's almost always an easier way using `bcftools`.

`bcftools view -i'ID=@snplist.txt' in.vcf`

---
Tags: vcf, snp, bcftools
---


# Question
Title: Visualization of mutations in (part of) viral genome - existing tools?
I'm trying to visualize mutations in SARS-CoV-2 sequences, somewhat like it's done on outbreak.info:

but for individual sequences, probably on the base of the vcf files generated by the ARTIC pipeline. A pipeline focused on virus nanopore output.

The ready-made tools I've found for visualization of mutations so far were specifically tailored to human sequences and/or for visualizing aggregated data for multiple genomes, or were web tools, which can't be used for medical data due to confidentiality.

Is there any existing offline tool that

* is Linux-based
* can be integrated into an automated pipeline
* is tailored towards viral genomes or is "organism-agnostic"
* is based on comparing an individual sequence to a reference sequence
* shows the sequence as "blocks" of genes with the mutations highlighted
* (ideally plays well with MultiQC for integration into an existing report, but that's a nice-to-have)

?

# Answer
> 3 votes
I ended up using `lolliplots` from the trackViewer R package; it can do VCF-based visualization *and* is flexible enough I can feed it data from other sources (I wound up parsing Nextclade mutation calls for the amino acids since I had these anyway).

---
Tags: vcf, visualization, covid-19, software-recommendation, sars-cov-2
---


# Question
Title: Help to create a dataframe in Python from a FASTA file
I want to create a dataframe in Python starting from a FASTA format file. Given the toy FASTA file that I am attaching, I built this program in Python that returns four colums corresponding to id, sequence length, sequence, animal name and rows corresponding to all the data available.

However, I am trying to understand how to modify this code in order to create a dataframe in which classes Human and Dog have the same number of data. For example, I want to say to Python: "Append to record (that is the empty list) id, sequence length, sequence and animal for Human, but do it a number of times that is equal to the number of data of the class with minimum number of data (that is Dog)".

I think that a `while` loop is needed but I have a bit troubles to understand how to do it. Any suggestion ?

Below the Python code I wrote and the FASTA format file I used.

```
import pandas as pd
import re
def read_fasta(file_path, columns) :
    from Bio.SeqIO.FastaIO import SimpleFastaParser
    with open("Proof.txt") as fasta_file :
        records = [] # create empty list
        for title, sequence in SimpleFastaParser(fasta_file): #SimpleFastaParser Iterate over Fasta records as string tuples. For each record a tuple of two strings is returned, the FASTA title line (without the leading ‘>’ character), and the sequence (with any whitespace removed).
            record = []
            title_splits=re.findall(r"[\w']+", title) # Data cleaning is needed

            record.append(title_splits[0])  #First values are ID (Append adds element to a list)
            record.append(len(sequence)) #Second values are sequences lengths
            sequence = " ".join(sequence) #It converts into one line
            record.append(sequence)#Third values are sequences

            #Fourth column will contain the species
            if "Human" in title_splits:
                    record.append("Human")
            else:
                    record.append("Dog")

            records.append(record)
    return pd.DataFrame(records, columns = columns) #We have created a function that returns a dataframe

#Now let's use this function by inserting in the first argument the file name (or file path if your working directory is different from where the fasta file is)
#And in the second one the names of columns
data = read_fasta("Proof.txt", columns=["id","sequence_length", "sequence", "animal"])
data

```

The FASTA format file is this:

```
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Dog|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Dog|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Dog|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Dog|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK

```

My code prints a dataframe like:

```
       id  sequence_length                               sequence animal
0   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
1   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
2   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
3   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
4   Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
5   Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
6   Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
7   Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
8   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
9   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
10  Numer               19  H S S F I E I V N I E H V I E H I V K  Human
11  Numer               19  H S S F I E I V N I E H V I E H I V K  Human

```

But I would like that the number of rows for Human is the same for Dog (because, in other words, I would like the same number of data for each class).

Hoping to have been clear, I thank you in advance.

P.S.: My program is an adapted version of this notebook .

# Answer
Assuming what I said in the comments is true. Here's my attempt, there may be an easier solution but this should do it:

```
import pandas as pd
import re

# a function to return the count number of animals in the fasta in a dictionary
def count_animals(file_path, separator, animal_index):
  animal_dictionary={}
  with open("Proof.txt") as fasta_file:
    fasta_contents = fasta_file.readlines()
    for line in fasta_contents:
      if line.startswith(">"):
        animal = line.split(separator)[animal_index]
        if animal not in animal_dictionary:
          animal_dictionary[animal] = 1
        else:
          animal_dictionary[animal] = animal_dictionary[animal] + 1
  return animal_dictionary

def read_fasta(file_path, columns) :
  counter = 0  # counter used to print no more than number of times 'Dog'
  animal_d = count_animals(file_path, '|', 1) # execute function
  max_animal =  max(animal_d.values()) # get max number of animals
  min_animal =  min(animal_d.values()) # get min number of animals

  from Bio.SeqIO.FastaIO import SimpleFastaParser
  with open("Proof.txt") as fasta_file :
    records = [] # create empty list
    for title, sequence in SimpleFastaParser(fasta_file): #SimpleFastaParser Iterate over Fasta records as string tuples. For each record a tuple of two strings is returned, the FASTA title line (without the leading ‘>’ character), and the sequence (with any whitespace removed).
      record = []
      title_splits=re.findall(r"[\w']+", title) # Data cleaning is needed

      record.append(title_splits[0])  #First values are ID (Append adds element to a list)
      record.append(len(sequence)) #Second values are sequences lengths
      sequence = " ".join(sequence) #It converts into one line
      record.append(sequence)#Third values are sequences
      # if 'Human' exists in title and the counter is less than min_animal('Dog'), add the record
      if 'Human' in title and counter < min_animal:
        counter += 1 # add +1 to counter

        #Fourth column will contain the species
        record.append("Human")

        records.append(record)
      elif 'Dog' in title:
        record.append("Dog")

        records.append(record)

    # print(records)
  return pd.DataFrame(records, columns = columns) #We have created a function that returns a dataframe

#Now let's use this function by inserting in the first argument the file name (or file path if your working directory is different from where the fasta file is)
#And in the second one the names of columns
data = read_fasta("Proof.txt", columns=["id","sequence_length", "sequence", "animal"])

print(data)

```

Should give:

```
      id  sequence_length                               sequence animal
0  Numer               19  H S S F I E I V N I E H V I E H I V K  Human
1  Numer               19  H S S F I E I V N I E H V I E H I V K  Human
2  Numer               19  H S S F I E I V N I E H V I E H I V K  Human
3  Numer               19  H S S F I E I V N I E H V I E H I V K  Human
4  Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
5  Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
6  Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
7  Numer               19  H S S F I E I V N I E H V I E H I V K    Dog

```
> 2 votes

# Answer
An example of contructing a pandas dataframe for fasta on Stack Exchange is **here** . This is clean code deploying `zip` outside Biopython. Constructing a pandas dataframe directly via dictionary.

Are you intending to use `split` method directly on the sequence, such that each residue is a separate column in a dataframe?

Anyway,

```
title_splits=re.findall(r"[\w']+", title) # Data cleaning is needed

```

This is a bad idea because you will remove the Genbank unique identifier (its an integer) and this should not leave the sequence.

You will need to switch out of the Fasta object. Your approacha and @user324810's use of

```
sequence = " ".join(sequence)

```

is okayish. However, I think it should be '', or "" (no space character). I spotted manually removing fasta, well the method should do this automatically. I am not familiar with Fastaparser, but AlignIO and SeqIO will both automatically process fasta without the need deal with "\>". Biopython will have a length method and whilst you have the Biopython object it is preferable to use Biopython (OOP) methods.

```
from Bio.SeqIO.FastaIO import SimpleFastaParser

```

This should be at the top of the code.

Overall, you manually open the data, pass it through the Bioparser then immediately dump the object. I do agree to write a pandas dataframe from a Biopython object these conversion to string, but Biopython will provide a method to print the object (pass it to a string) without all the parsing. Some of the parsing is removing key data (Genbank ID).

Personally I would look at the Stackoverflow example and base your code around this.

You clearly have an understanding of pandas, but I get the impression Biopython is something your not familiar with. My suspicion is that Biopython would infact answer the entirety of your question (whatever it might be) without passing to pandas.

pandas is cool and powerful for non-molecular data, but unless your using machine/deep learning (where a dataframe is necessary input) I think Biopython will provide most of the answers.

**Important Terminology** The class is FastaIO and the method is SimpleFastaParser. Humans and dogs are instances of this method and are therefore object(s).

Please don't write 'class' of dogs and humans in a formal report because it is incorrect OOP and will cause alarm if your reviewer is a coder. Your code is functional programming rather than OOP. You could say the human/dog fasta object, or just refer to it as human/dog dataframes and avoid any OOP terminology.
> 1 votes

# Answer
Just to add, usually the way to move a sequence object to a neat, decluttered string in Biopython is not to do this manually as done here:

```
title_splits=re.findall(r"[\w']+", title)
sequence = " ".join(sequence)

```

both being from `SimpleFastaParser(fastafile)`

It should be possible and, if so, better to state,

```
 str(title)
 str(sequence)

```

If you can issue a `print` command and the string appears nice and neat on the screen `str()` will work.

This becomes important for complex sequence/alignment objects, such as SeqIO and AlignIO. The reason it works is `__str__` has been declared or available through OOP inheritance.
> 0 votes

---
Tags: python, pandas, dataframe
---


# Question
Title: Translating Simpson and Shannon diversity indices?
So I am trying to compute diversity indices for a dataset I have but the math doesn't work out and I don't really understand how to translate Shannon and/or Simpson diversity indices into something that is understandable. Maybe I got everything wrong, this is certainly a possibility. However, this is what I have up to now:

```
# load libraries
library("phyloseq")
data("GlobalPatterns")

# load data
globp = GlobalPatterns

# compute transformation according to link below
simpson_transformed <- 1/(1-estimate_richness(globp, measures="Simpson"))
shannon_transformed <- exp(estimate_richness(globp, measures="Shannon"))

# print side by side
print(cbind.data.frame(simpson_transformed, shannon_transformed))

```

The **comparing\_indices** columns are computed according to what I found at this link (italic is mine):

> They *\[i.e. the values\]* should be converted to effective number of species, which are the true diversities. By referring to Table 1 we see that species richness is already a true diversity, so the true diversity according to this index is 5.000 species. Table 1 shows that the Shannon entropy or Shannon-Wiener index is converted to effective number of species or true diversity by taking the exponential: exp(1.609) is 5.000, so the true diversity according to the Shannon entropy is also 5.000 species. The Gini-Simpson index is converted to a true diversity by subtracting it from unity and inverting: 1/(1-0.8) = 5.000 species also. So in fact all these indices agree that the diversity of this community is 5.000 species.

I also verified how **phyloseq** is computing the indices. I believe the transformation for Shannon should be the same mentioned above. I have doubts about the Simpson, since they refer to the Gini-Simpson index and not the Simpson itself. However, the documentation for **vegan**'s `diversity()` function, which is the function that (I believe) is called by **phyloseq**'s `estimate_richness()` (see docs), mention that the Simpson index is returned as `1-D` which seems to be the Gini-Simpson index (see the `1-lambda`). In other words, it looks like the transformation I linked above should be exactly the one that are used in the code I provided. But, I was expecting the same values in the two columns. But they are not there and I wonder why this happens. What did I do wrong?

# Answer
> 2 votes
You are trying to obtain Hill numbers, sometimes called true diversity, and you have different results because you have uneven species frequencies.

`Vegan` calculates Shannon index like:

$$ H' = -\sum\_{i=1}^R p\_i \ln p\_i $$

And Gini-Simpson like:

$$ 1 - \lambda = 1 - \sum\_{i=1}^R p\_i^2 $$

So the Hill numbers should be:

$$ H1 = exp(H') $$

$$ H2 = 1 / 1 - \lambda $$

With `Phyloseq`:

```
H1 <- 1 / estimate_richness(globp, measures="Simpson")
H2 <- exp(estimate_richness(globp, measures="Shannon")

```

---
Tags: phyloseq
---


title: Getting error in running GATK

text:
# Question
Title: Getting error in running GATK
I am trying to making a fasta doc using GATK tools but I don't know what does this error say

Server2:/data2/RNASeq/Angel/gatk$ ./gatk CreateSequenceDictionary -R reference.fa Using GATK jar /data2/RNASeq/Angel/gatk/gatk-package-4.2.4.1-local.jar Running: java -Dsamjdk.use\_async\_io\_read\_samtools=false -Dsamjdk.use\_async\_io\_write\_samtools=true -Dsamjdk.use\_async\_io\_write\_tribble=false -Dsamjdk.compression\_level=2 -jar /data2/RNASeq/Angel/gatk/gatk-package-4.2.4.1-local.jar CreateSequenceDictionary -R reference.fa INFO 2022-01-28 15:04:19 CreateSequenceDictionary Output dictionary will be written in reference.dict 15:04:19.401 INFO NativeLibraryLoader - Loading libgkl\_compression.so from jar:file:/data2/RNASeq/Angel/gatk/gatk-package-4.2.4.1-local.jar!/com/intel/gkl/native/libgkl\_compression.so \[Fri Jan 28 15:04:19 GMT 2022\] CreateSequenceDictionary --REFERENCE reference.fa --TRUNCATE\_NAMES\_AT\_WHITESPACE true --NUM\_SEQUENCES 2147483647 --VERBOSITY INFO --QUIET false --VALIDATION\_STRINGENCY STRICT --COMPRESSION\_LEVEL 2 --MAX\_RECORDS\_IN\_RAM 500000 --CREATE\_INDEX false --CREATE\_MD5\_FILE false --GA4GH\_CLIENT\_SECRETS client\_secrets.json --help false --version false --showHidden false --USE\_JDK\_DEFLATER false --USE\_JDK\_INFLATER false Jan 28, 2022 3:04:19 PM shaded.cloud\_nio.com.google.auth.oauth2.ComputeEngineCredentials runningOnComputeEngine INFO: Failed to detect whether we are running on Google Compute Engine. \[Fri Jan 28 15:04:19 GMT 2022\] Executing as Server2 on Linux 4.15.0-166-generic amd64; OpenJDK 64-Bit Server VM 11.0.13+8-Ubuntu-0ubuntu1.18.04; Deflater: Intel; Inflater: Intel; Provider GCS is available; Picard version: Version:4.2.4.1 \[Fri Jan 28 15:04:19 GMT 2022\] picard.sam.CreateSequenceDictionary done. Elapsed time: 0.00 minutes. Runtime.totalMemory()=2155872256 To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp picard.PicardException: /data2/RNASeq/Angel/gatk/reference.dict already exists. Delete this file and try again, or specify a different output file. at picard.sam.CreateSequenceDictionary.doWork(CreateSequenceDictionary.java:220) at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:308) at org.broadinstitute.hellbender.cmdline.PicardCommandLineProgramExecutor.instanceMain(PicardCommandLineProgramExecutor.java:37) at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:160) at org.broadinstitute.hellbender.Main.mainEntry(Main.java:203) at org.broadinstitute.hellbender.Main.main(Main.java:289) Server2:/data2/RNASeq/Angel/gatk$

# Answer
```
picard.PicardException: /data2/RNASeq/Angel/gatk/reference.dict already exists.  Delete this file and try again, or specify a different output file.

```

I think that requires no further explanation.
> 2 votes

---
Tags: gatk, picard, wes
---


# Question
Title: tensorflow nn_model for DNA sequences: Matrix size-incompatible: In[0]: [2,1], In[1]: [784,300]
Hope anyone can help a beginner here. I'm building a proof-of concept tensorflow classifier for DNA sequences. However, the NN model does not let through train and test vectors saying the matrix size is incompatible. Where could be the mistake in my implementation?

```
import numpy as np
import pandas as pd

df=pd.DataFrame(columns=['sequences', 'label'])
a='ttttccagaattctcttagttt gtgatgtctttattgcttctattt'
b='ctcctgcttgctttttttcttg ggtttctgatattctttaaaggat'
c='tcctgcttgctttttttcttgg gtttctgatattctttaaaggatt'
df.sequences=[a,b,c]
df.label=[1,1,1]

from sklearn.preprocessing import LabelEncoder, OneHotEncoder
integer_encoder = LabelEncoder()
one_hot_encoder = OneHotEncoder(categories='auto')
input_features = []
sequences=df.sequences
for sequence in sequences:
     integer_encoded = integer_encoder.fit_transform(list(sequence))
     integer_encoded = np.array(integer_encoded).reshape(-1, 1)
     one_hot_encoded = one_hot_encoder.fit_transform(integer_encoded)
     input_features.append(one_hot_encoded.toarray())

np.set_printoptions(threshold=40)
input_features = np.stack(input_features)
print("Example sequence\n-----------------------")
print('DNA Sequence #1:\n',sequences[0][:10],'...',sequences[0][-10:])
print('One hot encoding of Sequence #1:\n',input_features[0].T)

labels=df.label
one_hot_encoder = OneHotEncoder(categories='auto')
labels = np.array(labels).reshape(-1, 1)
input_labels = one_hot_encoder.fit_transform(labels).toarray()
print('Labels:\n',labels.T)
print('One-hot encoded labels:\n',input_labels.T)

from sklearn.model_selection import train_test_split
train_features, test_features, train_labels, test_labels = train_test_split(input_features, input_labels, test_size=0.25, random_state=42)

def get_batch(x_data, y_data, batch_size):
     idxs = np.random.randint(0, len(y_data), batch_size)
     return x_data[idxs,:,:], y_data[idxs]

epochs = 10
batch_size = 100

x_train = train_features
x_test = train_labels
x_train = x_train / 255.0
x_test = x_test / 255.0
x_test = tf.Variable(x_test)

W1 = tf.Variable(tf.random.normal([784, 300], stddev=0.03), name='W1')
b1 = tf.Variable(tf.random.normal([300]), name='b1')
W2 = tf.Variable(tf.random.normal([300, 10], stddev=0.03), name='W2')
b2 = tf.Variable(tf.random.normal([10]), name='b2')

def nn_model(x_input, W1, b1, W2, b2):
     x_input = tf.reshape(x_input, (x_input.shape[0], -1))
     x = tf.add(tf.matmul(tf.cast(x_input, tf.float32), W1), b1) #ERROR FROM HERE?
     x = tf.nn.relu(x)
     logits = tf.add(tf.matmul(x, W2), b2)
     return logits

def loss_fn(logits, labels):
     cross_entropy = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(labels=labels,logits=logits))
     return cross_entropy

optimizer = tf.keras.optimizers.Adam()
total_batch = int(len(y_train) / batch_size)

for epoch in range(epochs):
    avg_loss = 0
    for i in range(total_batch):
        batch_x, batch_y = get_batch(x_train, y_train, batch_size=batch_size)
        # create tensors
        batch_x = tf.Variable(batch_x)
        batch_y = tf.Variable(batch_y)
        # create a one hot vector
        batch_y = tf.one_hot(batch_y, 10)
        with tf.GradientTape() as tape:
            logits = nn_model(batch_x, W1, b1, W2, b2)
            loss = loss_fn(logits, batch_y)
        gradients = tape.gradient(loss, [W1, b1, W2, b2])
        optimizer.apply_gradients(zip(gradients, [W1, b1, W2, b2]))
        avg_loss += loss / total_batch
    test_logits = nn_model(x_test, W1, b1, W2, b2) ##ERROR LOG APPEARS HERE
    max_idxs = tf.argmax(test_logits, axis=1)
    test_acc = np.sum(max_idxs.numpy() == y_test) / len(y_test)
    print(f"Epoch: {epoch + 1}, loss={avg_loss:.3f}, test set accuracy = {test_acc*100:.3f}%")
print("\nTraining complete!")

```

> Traceback (most recent call last): File "", line 16, in File "", line 4, in nn\_model File "/usr/local/lib/python3.8/dist-packages/tensorflow/python/util/traceback\_utils.py", line 153, in error\_handler raise e.with\_traceback(filtered\_tb) from None File "/usr/local/lib/python3.8/dist-packages/tensorflow/python/framework/ops.py", line 7107, in raise\_from\_not\_ok\_status raise core.\_status\_to\_exception(e) from None # pylint: disable=protected-access tensorflow.python.framework.errors\_impl.InvalidArgumentError: Matrix size-incompatible: In\[0\]: \[2,1\], In\[1\]: \[784,300\] \[Op:MatMul\]

# Answer
> 0 votes
The sci-kit learn looks fine, for a proof of concept. Furthermore, the deep learning looks good, e.g. Adams, loss etc .. all good. Of course it at ain't RNN or CNN but its proof of concept.

However, the error you encountered suggests it is about resizing the numpy array and that is better answered on Stackoverflow in the absence of an answer here. Dealing with `pandas` is one thing, munging a `numpy array` is pretty hardcore. Its a single line of code, but you've got to know it exactly.

An example is where I d focus,

```
tf.reshape(x_input, (x_input.shape[0], -1))

```

Before or whether you need to post there please check:

1. Where is `import tensorflow as tf` ?
2. Using keras without e.g. `from tensorflow.keras import optimizers` (well okay you can do that given point 1 but generally the keras method is specifically imported)

If I haven't misread the code it will fall over, as soon as it gets past `sci-kit learn`. The advantage of placing all imports at the top of the code is that its easy to read whats going on. There is another issue, but if a big deal Stackoverflow (or someone else here) will flag it.

---

Just to we explain point 2, I suspect the tensor flow on question is Tensorflow 2 if this is true then Keras is built in and thus tf.keras is fine

---
Tags: python, genome, sequence-analysis, encode, tensorflow
---


# Question
Title: Are there different PAM matrices for different protein families? - How is the inferred common ancestor known to start with?
I am confused about the sequences one starts with to generate the PAM matrices.

The sequences to start with are supposed to diverge by no more than 1% (in sequence identity) wrt a common "inferred" ancestor. - if I understood correctly.

How is the common ancestor sequence inferred without having the means to do any sequence alignment to start with?

Does this mean that a PAM matrix is specific for an inferred ancestor, i.e. a protein family? - I got the impression there was only one in practice.

Thanks in advance,

Gordon

# Answer
Ideally, but I wouldn't worry because the chance of meeting the criteria in practice are slim. Most importantly you need a tree which is 'independently' derived from the protein sequence family in question. The PAM matrix is absolutely dependent on this. Generating a "true" tree without reference to the sequences you are studying isn't trivial, so it's got to generated with a different data set and the question is whether that data set gives a true tree, if you're stuck the DNA sequence of is a rough approximation (not ideal but you could justify it). If you can adequately do this, I wouldn't worry about the remaining assumptions because it's simply an optimisation which can't proven.

The PAM matrix generation will estimate the sequence of a common ancestor, based on a given tree. What this criteria is saying is that if you want a 'perfect matrix' the sequence reconstruction of the ancestral node needs 'x' because that is pivotal to calculating the matrix. If that's a bit out your matrix is a bit by a margin of error that isn't quantified (maybe it has been, see below)

Whether you can get get something close to the 'hub' sequence within your data set is one question and whether that sequence was included in the first place is all trial and error. I would assume the criteria are ascertained has been set by simulation studies (its fairly easy to do if you've got the time), which only affects the margin of confidence within the confines of that given simulation. PAM matrix generation has been around since the dawn of time and founding principle is an independent tree, everything else you can calculate and arm-wave, ie admit the calculation may not be so precise. In other it doesn't stop you doing it and presenting the results, that is perfectly valid, as long as simulated criteria are understood. The problem would be is if you ignore a past study and present your results 'without fault'. The biggest problem would be if you generated a tree and used the same data.to generate the matrix, because that is whatever matrix you used to generate the tree and that's a no no.
> 0 votes

---
Tags: sequence-alignment, point-mutation
---


# Question
Title: No bootstrap value on single-copy gene tree created by OrthoFinder
I'm running analysis with OrthoFinder and it produces a single-copy gene tree. When I visualize the tree with iTOL, there's no bootstrap value on the branches or nodes. There's a published paper in BMC Genomics that also didn't show the bootstrap value using OrthoFinder.

https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-019-6065-7#Tab1

1. Does anybody know why there's no bootstrap values present on the phylogeny?
2. How can this be corrected?

Thanks

# Answer
> 1 votes
Please check the file

```
 `$PWD/SpeciesTree_unrooted.txt" `

```

You should find the bootstrap support\* values per branch there. If not make sure is the latest version is installed from https://github.com/davidemms/OrthoFinder

The papers are: genomebiology.biomedcentral.com/articles/10.1186/… and genomebiology.biomedcentral.com/articles/10.1186/… Genome Biology 16:157 ... It officially offers tree stats: "OrthoFinder also provides comprehensive statistics for comparative genomic analyses"

The documentation is here

---

The question was raised how to increase replicates by one, two or three logs. This is not really orthodox technique. The way the issue should be first addressed is the user performs replicates manually, e.g. 10 new runs, and then looks across the consistency of bootstrap values for each manual replicate. If all values of interest between manual replicates are \>80% or \<70% it really doesn't matter. However if a key branch ranges from 70 - 75% that is strictly unresolved and would generally be resolved using an alternative test. If the variation is \>= 10% that would need addressing. Generally variation between replicates is very low =\< 2%

---

* note, the original answer just contained the word 'support' and has now been more correctly described as "bootstrap support". Thanks to the OP for pointing this out in the comments below.

---
Tags: phylogenetics, phylogeny, orthofinder
---


title: Error in running Mutect2

text:
# Question
Title: Error in running Mutect2
I am trying to run Mutect2 on a .bam file but I get an error by Googling can not be tackled

Have you seen this error?

```
fi1d18@RBGO-Server2:/data2/RNASeq/Angel/gatk$ ./gatk Mutect2 -R reference.fa -I sample.bam --germline-resource af-only-gnomad.hg38.vcf.gz --panel-of-normals pon.vcf.gz -O single_sample.vcf.gz
Using GATK jar /data2/RNASeq/Angel/gatk/gatk-package-4.2.4.1-local.jar
Running:
    java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /data2/RNASeq/Angel/gatk/gatk-package-4.2.4.1-local.jar Mutect2 -R reference.fa -I sample.bam --germline-resource af-only-gnomad.hg38.vcf.gz --panel-of-normals pon.vcf.gz -O single_sample.vcf.gz
13:16:43.616 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/data2/RNASeq/Angel/gatk/gatk-package-4.2.4.1-local.jar!/com/intel/gkl/native/libgkl_compression.so
Jan 31, 2022 1:16:43 PM shaded.cloud_nio.com.google.auth.oauth2.ComputeEngineCredentials runningOnComputeEngine
INFO: Failed to detect whether we are running on Google Compute Engine.
13:16:43.733 INFO  Mutect2 - ------------------------------------------------------------
13:16:43.734 INFO  Mutect2 - The Genome Analysis Toolkit (GATK) v4.2.4.1
13:16:43.734 INFO  Mutect2 - For support and documentation go to https://software.broadinstitute.org/gatk/
13:16:43.734 INFO  Mutect2 - Executing as fi1d18@RBGO-Server2 on Linux v4.15.0-166-generic amd64
13:16:43.734 INFO  Mutect2 - Java runtime: OpenJDK 64-Bit Server VM v11.0.13+8-Ubuntu-0ubuntu1.18.04
13:16:43.734 INFO  Mutect2 - Start Date/Time: 31 January 2022 at 13:16:43 GMT
13:16:43.734 INFO  Mutect2 - ------------------------------------------------------------
13:16:43.734 INFO  Mutect2 - ------------------------------------------------------------
13:16:43.735 INFO  Mutect2 - HTSJDK Version: 2.24.1
13:16:43.735 INFO  Mutect2 - Picard Version: 2.25.4
13:16:43.735 INFO  Mutect2 - Built for Spark Version: 2.4.5
13:16:43.735 INFO  Mutect2 - HTSJDK Defaults.COMPRESSION_LEVEL : 2
13:16:43.735 INFO  Mutect2 - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
13:16:43.735 INFO  Mutect2 - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
13:16:43.735 INFO  Mutect2 - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
13:16:43.735 INFO  Mutect2 - Deflater: IntelDeflater
13:16:43.735 INFO  Mutect2 - Inflater: IntelInflater
13:16:43.735 INFO  Mutect2 - GCS max retries/reopens: 20
13:16:43.735 INFO  Mutect2 - Requester pays: disabled
13:16:43.735 INFO  Mutect2 - Initializing engine
13:16:44.015 INFO  Mutect2 - Shutting down engine
[31 January 2022 at 13:16:44 GMT] org.broadinstitute.hellbender.tools.walkers.mutect.Mutect2 done. Elapsed time: 0.01 minutes.
Runtime.totalMemory()=2155872256
***********************************************************************

A USER ERROR has occurred: Cannot read file:///data2/RNASeq/Angel/gatk/pon.vcf.gz because no suitable codecs found

***********************************************************************
Set the system property GATK_STACKTRACE_ON_USER_EXCEPTION (--java-options '-DGATK_STACKTRACE_ON_USER_EXCEPTION=true') to print the stack trace.

```

I also get this

```
fi1d18@RBGO-Server2:/data2/RNASeq/Angel/gatk$ file pon.vcf.gz
pon.vcf.gz: HTML document, ASCII text, with very long lines

```

# Answer
Your input file `pon.vcf.gz` should be a gzip file, but from your `file` command it looks like it is an HTML file instead. Perhaps you downloaded the VCF file incorrectly?

Here is what I get on a bona fide gzip file in Ubuntu:

```
$ echo hi | gzip -c > hi.gz
$ file hi.gz
hi.gz: gzip compressed data, from Unix, original size modulo 2^32 3

```
> 2 votes

# Answer
Downloading that vcf file from this link fixes the problem (but generates a new one)

```
wget https://storage.googleapis.com/gatk-best-practices/somatic-hg38/1000g_pon.hg38.vcf.gz

```
> 1 votes

---
Tags: gatk, mutations, wes
---


# Question
Title: (scRNA-seq) I have a dataset with 4 conditions (experiments) should I be anchoring aka integrating data for each comparison when looking for DEGs?
I have an experiment with 4 conditions in it two of them are controls. When I am correcting for batch effects (either with PCA or CCA) should I restrict my comparisons between particular groups? So first I should do control 1 with control 2 then hopefully there will be no DEGs between them. Next I could do Control 2 vs Condition 1 then condition 1 vs condition2 and control 2 vs condtion 2. Sorry if this is confusing but I am trying to avoid going into the experimental details in the interest of keeping this post brief.

As a sort of follow up question can anyone suggest the "standard" or best seurat method for correcting batch effects. I've found some tutorials that recommend using CCA and others suggest PCA.

# Answer
> 1 votes
Batch correction/data integration is a tricky subject for scRNA-seq, see a few recent papers for benchmarks: https://www.nature.com/articles/s41592-021-01336-8 and https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1850-9

Batch correction methods may not even impact differential expression. For example, if you correct the PCs, via a popular method like harmony, then this will impact the t-SNE/UMAP plots but will not impact the actual expression values in the counts matrix used to detect DEGs.

For detecting DEGs, the most common methods are simple stats tests like Wilcoxan Rank test that largely ignore batch effects. However, these are likely not the best approaches. I suggest this recent paper (https://www.nature.com/articles/s41467-021-25960-2). Here they state that the best methods for detecting DEGs consider biological replicates, like is done with bulk RNA-seq. By considering replicates, this is effectively like batch correction

---
Tags: scrnaseq, seurat, single-cell, batch-effects
---


# Question
Title: Bio.codonalign.codonseq module (cal_dn_ds) for SARS-COV-2
I have calculated the dN,dS ratio for a SARS-COV-2 data set ('NC\_045512.2' (Wuhan strain) and 'LC666924.1'), and results in a zero value (for most of the genes).

Is it possible and what is the interpretation, for example negative selection?

# Answer
Yes that is the expected result, because there are no amino acid mutations in the data sets you are using. Thus the numerator is 0 and anything divided by zero is always zero, even if every codon had a synonymous mutation.

The most precise interpretation is simply there are no amino acid mutations. It doesn't matter how you reword that, all divergence is resulting from synonymous mutations.

The only thing to watch out for is a rounding error, i.e. a single amino acid change is dropping off detection because dN/dS is incorrectly rounded.

---

A post was raised (now deleted) about what does that say about selection. dN/dS theory states that if it is significantly \>1 thats positive selection .. etc. However, a truer representation of the data is there is no amino acid mutation, rather than their is a significant excess of synonymous to non-synonymous mutations. A zero value makes things tricky, and is exemplified by two hypothetical data sets one where every codon had one synoymous mutation (strong negative selection) compared against one that had a single synonymous mutation (little mutation) and both data sets had zero amino acid mutations - both data sets would have a dN/dS of zero but that is not a reflection of negative selection.

---

To answer the point below a single mutation is not a great basis draw statistical significance and the dN/dS ratio is not a good test if either the numerator or denominator are zero.

Most data is under negative selection in strict accordance with dN/dS.

Neutral selection is extremely rare.

**Summary** Sometimes the statistical analysis of a data set does not say very much ... and sometimes we have to accept that as an answer.
> 1 votes

# Answer
So if that had a little mutation, Does that mean it's a neutral selection? And if the sequences (codons and amino) are exactly the same? We can say something about the interpretation?

See examples below-

**ORF8** Seq1: MKFLVFLGIITTVAAFHQECSLQSCTQHQPYVVDDPCPIHFYSKWYIRVGARKSAPLIELCVDEAGSKSPIQYIDIGNYTVSCLPFTINCQEPKLGSLVVRCSFYEDFLEYHDVRVVLDFI Seq2: MKFLVFLGIITTVAAFHQECSLQSCTQHQPYVVDDPCPIHFYSKWYIRVGARKSAPLIELCVDEAGSKSPIQYIDIGNYTVSCLPFTINCQEPKLGSLVVRCSFYEDFLEYHDVRVVL--I

**ORF10** Seq1 MGYINVFAFPFTIYSLLLCRMNSRNYIAQVDVVNFNLT Seq2 MGYINVFAFPFTIYSLLLCRMNSRNYIAQVDVVNFNLT

**E** (The same) Seq1 MYSFVSEETGTLIVNSVLLFLAFVVFLLVTLAILTALRLCAYCCNIVNVSLVKPSFYVYSRVKNLNSSRVPDLLV Seq2 MYSFVSEETGTLIVNSVLLFLAFVVFLLVTLAILTALRLCAYCCNIVNVSLVKPSFYVYSRVKNLNSSRVPDLLV
> 0 votes

---
Tags: python, phylogenetics, phylogeny, sars-cov-2, covid-19
---


# Question
Title: Is use_dimred deprecated or is this a runTSNE() bug?
NOTE: I posted the same question on Biostars but so far no replies there so I figured I would ask here, too.

So, I was running a pipeline which I developed some time ago and a very weird behaviour was found, which had quite an impact over our analyses.

Indeed, the scater's runTSNE() function was subjected to a change in the parameters names and, in particular, the old use\_dimred parameter now has a new name, which is dimred.

Interestingly, I was running the pipeline using the runTSNE() function and passing the dimensionality reduction through use\_dimred which is apparently deprecated and obtaining very weird results without being aware of that since no warning and/or error was returned.

A quick example:

```
# load some data
library("scRNAseq")
library("scater")
sce.416b <- LunSpikeInData(which="416b")

```

check the data and notice no dimensionality reduction is available:

```
> sce.416b
class: SingleCellExperiment
dim: 46604 192
metadata(0):
assays(1): counts
rownames(46604): ENSMUSG00000102693 ENSMUSG00000064842 ...
  ENSMUSG00000095742 CBFB-MYH11-mcherry
rowData names(1): Length
colnames(192): SLX-9555.N701_S502.C89V9ANXX.s_1.r_1
  SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 ...
  SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1
  SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1
colData names(9): Source Name cell line ... spike-in addition block
reducedDimNames(0):            <-- HERE IT IS EMPTY
mainExpName: endogenous
altExpNames(2): ERCC SIRV

```

at this point I simply run the runTSNE() function using a non-existent dimred:

```
# set seed
set.seed(001)
# add logcounts
sce.416b <- logNormCounts(sce.416b)
# add random colours to cells
colLabels(sce.416b) <- rep(c(1,2,3,4), 48)
# runTSNE, first try
sce.416b <- runTSNE(sce.416b, use_dimred="made_up_dimred", perplexity=30)
# and plot
dev.new(); plotTSNE(sce.416b, colour_by="label")
# runTSNE, second try
sce.416b <- runTSNE(sce.416b, use_dimred="another_made_up_dimred", perplexity=30)
# and plot
dev.new(); plotTSNE(sce.416b, colour_by="label")
# runTSNE, third try
sce.416b <- runTSNE(sce.416b, use_dimred="a_third_made_up_dimred", perplexity=30)
# and plot
dev.new(); plotTSNE(sce.416b, colour_by="label")

```

it, unfortunately, works. Indeed, I now have three TSNE available that I did plot and that are slightly different with each other.

```
> sce.416b
class: SingleCellExperiment
dim: 46604 192
metadata(0):
assays(2): counts logcounts
rownames(46604): ENSMUSG00000102693 ENSMUSG00000064842 ...
  ENSMUSG00000095742 CBFB-MYH11-mcherry
rowData names(1): Length
colnames(192): SLX-9555.N701_S502.C89V9ANXX.s_1.r_1
  SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 ...
  SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1
  SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1
colData names(10): Source Name cell line ... block sizeFactor
reducedDimNames(1): TSNE      <-- THERE YOU GO: A NICE TSNE COMPUTED OVER A NON-EXISTENT DIMRED
mainExpName: endogenous
altExpNames(2): ERCC SIRV

```

What is even more interesting is that if I run the function with the newly renamed parameter, i.e. dimred, I do get an error:

```
> sce.416b <- runTSNE(sce.416b, dimred="made_up_dimred", perplexity=30)
Error in value[[3L]](cond) :
  invalid subscript 'type' in 'reducedDim(<SingleCellExperiment>, type="character", ...)':
  'made_up_dimred' not in 'reducedDimNames(<SingleCellExperiment>)'

```

But then, if you have an old pipeline which used the runTSNE function with the use\_dimred and you're now running an updated version of scater, I have no idea what you may get.

So, as per my question...is this a bug? Has anyone else noticed this weird behaviour? Am I wrong?

# Answer
I believe it would be easier to answer the question by dividing it into related but more specific questions. Having said that, even if the answers below would be satisfactory, I would still post this on the package's Github repo, the author being informed about this can help other users of this package.

Question 1: Why does `runTSNE()` work when I use the deprecated `use_dimred` argument?

Let's create function and call it:

```
> my_func <- function(dimred) {
+   if(missing(dimred)) {
+     print("dimred argument is not provided so something else is being done...")
+   } else {
+     print(dimred)
+   }
+ }
> my_func(42)
[1] 42
> my_func(dimred = 42)
[1] 42
> my_func()
[1] "dimred argument is not provided so something else is being done..."
> my_func(use_dimred = 42)
Error in my_func(use_dimred = 42) : unused argument (use_dimred = 42)

```

As expected, we get an error if we call a function with an argument that does not belong to our function: `use_dimread`. If we define our function with the special argument `...` however, it works:

```
> my_func <- function(dimred, ...) {
+   if(missing(dimred)) {
+     print("dimred argument is not provided so something else is being done...")
+   } else {
+     print(dimred)
+   }
+ }
> my_func(42)
[1] 42
> my_func(dimred = 42)
[1] 42
> my_func()
[1] "dimred argument is not provided so something else is being done..."
> my_func(use_dimred = 42)
[1] "dimred argument is not provided so something else is being done..."

```

Basically when you call a function, the function you are calling can make calls to other functions and the `...` is basically allowing the user to feed arguments for these let's say "secondary" functions. Apparently `R` does not care if what you feed as an argument to a function is used or not if the function is defined with the special `...` argument.

Question 2: Why is `runTSNE()` working if the supplied `use_dimred` argument was not used (as shown above)?

If you take a look at the source code, `calculateTSNE()`, which `runTSNE()` runs under the hood, is using PCA if the `dimred` argument is not provided, its default is NULL anyway. It is been a while since I have used `Scater` but I think this would mean that the tSNE is being computed on (some number of) PCs.

And as a conclusion, I don't think we can really call this a bug but a rather unexpected behavior because of use of the special `...` argument in the function definition.
> 1 votes

---
Tags: single-cell, tsne
---


title: Some definitions about RNA-seq

text:
# Question
Title: Some definitions about RNA-seq
I want to select a 10x single cell kit

What does `2x 50 75 100 150 250` mean in paired end sequencing?

# Answer
> 2 votes
I am pretty sure that is read length. Both reads in a pair (2x) can be 50, 75, ... 250 bases long. And I am even more sure you will find this information on the webpage of the provider and if not, you should email them.

---
Tags: rna-seq, scrnaseq, 10x-genomics
---


# Question
Title: Is there a publication database or search engine offering geneID or UniProtID correlation?
Do we have a publication database or search engine that offers geneID or UniProtID correlation?

For example, I search "IPS cells" -\> I want to know what genes or proteins are under active research in IPS cells.

# Answer
I think this could be close to what you're looking for: http://glad4u.zhang-lab.org/index.php

Paper: GLAD4U: deriving and prioritizing gene lists from PubMed literature

It works with gene names though, so I am not sure it's a 100% fit.
> 4 votes

# Answer
Take a look at **GeneWalk**, https://churchman.med.harvard.edu/genewalk. It's a very interesting solution to this rather difficult problem. It is well documented and packaged. But It's open source software -- not a service. You'll have to run it on your own system.
> 1 votes

---
Tags: proteins, gene, database, phylogenetics, uniprot
---


# Question
Title: How to solve correlation problems between different samples in scRNA-seq?
I am trying to align and merge different samples from NCBI. I end up having correlation problem with these sample. The picture below shows an heatmap of the  by doing a linear regression between 2 samples' count matrices.

Samples **A**, **B** and **C** are from the same project and sample **D** from another. A, B and C are sequenced using the *10X Genomics for Single Cell V(D)J + 5' Gene Expression*, and *Chromium Next GEM Single Cell 3ʹ v3.1* for sample D.

All samples have been recovered from NCBI using `fastq-dump` — internaly threaded to be more efficient — and aligned using `cellranger` with the default human genome given by 10x.

QC shows some GC content small shifts along side some huge duplication level — even for scRNA-seq.

Am I missing something important when it comes to merge datasets?

Thank you very much

# Answer
> 1 votes
> The picture below shows an heatmap of the R² by doing a linear regression between 2 samples' count matrices.

In addition to doing a heatmap, I'd do a conventional UMAP/t-SNE plot and see how much seperation of the samples there are across all the different clusters/cell types.

> Am I missing something important when it comes to merge datasets?

Certainly there are 'technology batch effects' when comparing 3' and 5' data from the same tissue type. You can try to correct for these when merging the dataset by using a batch correction/integration tool (see reviews  and ), but will likely be imperfect.

> QC shows some GC content small shifts along side some huge duplication level — even for scRNA-seq.

If you're using a tool like fastQC then note that high duplication levels are expected for (single cell) RNAseq from highly expressed genes.

---
Tags: rna-seq, scrnaseq, single-cell, phylogenetics, cellranger
---


# Question
Title: How to submit a canu job on LSF high-performance computing cluster farm?
I am currently running canu on an LSF Linux server using the following script called assemble.sh:

```
canu -p species -d output genomeSize=100m -pacbio-raw subreads.fastq

```

I then submit this script with the following in a computer node:

```
bsub -q basement -E 'test -e /nfs/users/nfs_c/cr' -R "select[mem>25000] rusage[mem=25000]" -M25000 -o canu.o -e canu.e -J canu 'bash assemble.sh' -n8 -R"span[hosts=1]"

```

On my laptop, using 8 cores and 16 GB memory, assemble.sh took about 48h. On the cluster, it's been running for 72h and it is still running. I wonder if I am allocating resources correctly or if I can do something better. On my laptop, I can see in my activity monitor that canu is using many cores at the same time. Do you know how I can check this in an LSF server?

This page talks about submitting jobs on LSF. In particular, it mentions this: `gridOptions="your options list"` but I don't know how to use this option.

Any help would be appreciated.

checking the queue with `bjobs -u all` shows there are 125259 jobs running. My own job is running, it's not pending. Does the number of other jobs affect the running of my own job? I will ask IT if there are preferred writing locations. I did not consider that could change the speed of the job. Do you know a command to check how many resources my job is using?

# Answer
I haven't used canu for several years but looking back at my labbook I ran it with two additional flags:

`-minThreads=8 -maxThreads=8`

My memory on the subject is a bit hazy but think I had to play around with these settings to get it to work. I also set the flag:

`-maxMemory=6`
> 1 votes

---
Tags: assembly, linux, canu
---


title: Definition of "proband" in CNV analysis

text:
# Question
Title: Definition of "proband" in CNV analysis
I am reading this paper regarding *de novo* CNV analysis in Schizophrenia:

> Accordingly, such CNVs must frequently occur as recurrent de novo mutations. In a sample of 662 schizophrenia **proband–parent** trios, we found that rare de novo CNV mutations were significantly more frequent in cases (5.1% all cases, 5.5% family history negative) compared with 2.2% among 2623 controls, confirming the involvement of de novo CNVs in the pathogenesis of schizophrenia.

I assumed that **proband** stands for the sample being tested (child of the parents), however, I couldn't find an explicit definition for it. Is it just the sample being tested, or does it imply additional details?

# Answer
> 1 votes
*answer from @maximilian-press, converted from comments:*

> A **proband** is an individual affected with a disorder who is the first subject in a study (as of a genetic character in a family lineage) who brings a genetic disorder to the attention of the medical community.

---
Tags: cnv
---


# Question
Title: How to Sort and Index a SAM file without converting it to BAM?
Generally, I use `samtools` to `sort` and `index` BAM files. Samtools works also with BGZF compressed SAM but not with plain ones.

I have terabytes of plain SAMs at the moment and I would love to avoid compression/conversion. I have all the data in several external HD and I am not allowed to delete the original SAM files. I would need much more space than I have now to also store the BAM files.

So, I was wondering if there is a way to avoid it. Is there some tool/trick to sort and index uncompressed SAM files, without converting them to BAM or compressing them?

# Answer
> 7 votes
Not only will you save disk space by converting to BAM, but BAM files are faster to manipulate than SAM. Source: Dave Tang's SAMTools wiki.

`sort` supports uncompressed SAM format from a file or `stdin`, though `index` requires BGZIP-compressed SAM or BAM. I don't think you can get around this.

If you're dead-set on reading from and writing to uncompressed SAM format, you can use pipes to avoid writing BAM files. Specify the output of `view` to be an uncompressed bam, then pipe it to your next command, then pipe it to `view` again to output your sorted SAM file.

```
samtools view -u in.sam | samtools sort | samtools view -h > sort.sam
samtools view -u sort.sam | samtools index

```

---
Tags: samtools, sam, indexing
---


# Question
Title: fastqCleaner failing to launch in RStudio
When I execute:

`> FastqCleaner:::launch_fqc()()`

I get the following output

```
Listening on http://127.0.0.1:4863
Warning: Navigation containers expect a collection of `bslib::nav()`/`shiny::tabPanel()`s and/or `bslib::nav_menu()`/`shiny::navbarMenu()`s. Consider using `header` or `footer` if you wish to place content above (or below) every panel's contents.
Warning: Error in : Navigation containers expect a collection of `bslib::nav()`/`shiny::tabPanel()`s and/or `bslib::nav_menu()`/`shiny::navbarMenu()`s. Consider using `header` or `footer` if you wish to place content above (or below) every panel's contents.
  85: stop
  84: FUN
  83: lapply
  82: buildTabset
  81: tabsetPanel_
  80: func
  78: tabsetPanel
Warning: Navigation containers expect a collection of `bslib::nav()`/`shiny::tabPanel()`s and/or `bslib::nav_menu()`/`shiny::navbarMenu()`s. Consider using `header` or `footer` if you wish to place content above (or below) every panel's contents.
Warning: Error in : Navigation containers expect a collection of `bslib::nav()`/`shiny::tabPanel()`s and/or `bslib::nav_menu()`/`shiny::navbarMenu()`s. Consider using `header` or `footer` if you wish to place content above (or below) every panel's contents.
  85: stop
  84: FUN
  83: lapply
  82: buildTabset
  81: tabsetPanel_
  80: func
  78: tabsetPanel

```

Have I done something wrong or is the function broken in some way?

# Answer
> 2 votes
The issue is with the navbar page, downgrading shiny to 1.6 or in my case 1.4 fixed it for me, the detailed link to find to the solution is here: https://github.com/rstudio/shiny/issues/3556

Cheers.

---
Tags: r, bioconductor, fastq, fastqc, rstudio
---


title: What is DSSP algorithm for?

text:
# Question
Title: What is DSSP algorithm for?
The above article says that DSSP is a hydrogen bond estimation algorithm. Then the article also says that this algorithm is "assigns secondary structure to proteins".

If both are true, how does hydrogen bond estimation helps us to assign secondary structures of proteins?

# Answer
**"Secondary structure is formally defined by the pattern of hydrogen bonds between the amino hydrogen and carboxyl oxygen atoms in the peptide backbone"** (source)

To be more specific: secondary structures are formed explicitly because of hydrogen bonds. So when DSSP says that it estimates hydrogen bonds as a proxy for secondary structure classification, this is likely why.
> 3 votes

---
Tags: protein-structure
---