cannin/biostars_qa huggingface dataset:
Question: Biopython comes with methods to interface with the PAML package for phylogenetic analysis.

In particular I am using Bio.Phylo.PAML to run analyses using PAML's codeml.exe program which in my case does Ka/Ks (dN/dS) ratio analysis on pairs of orthologous gene sequences.

After running the analysis using results = cml.run() I can see it has successfully generated result.out files that look about right. Most importantly is the final line of the file that I need to parse into Python:

t= 0.2173 S= 703.9 N= 1489.1 dN/dS= 0.2247 dN = 0.0344 dS = 0.1529

What I am most interested in is dN/dS = 0.2247

According to Biopython's PAML wiki this value can be obtained from Python by doing results = cml.run() which generates a dictionary with a set of values I am interested in after running the analysis. The wiki claims I can find the values I need in a key called 'parameters'. But this only returns one of the values I need t= 0.2173, look:

>>> results.values()
['Fcodon', 'One dN/dS ratio for branches, ', '4.7b', {0: {'description': 'one-ratio', 'parameters': {'t': 0.1982}}}, {'htlv': {}, 'stlv': {}}]

Notice, that my parameters key only contains the t= 0.2173 and has omitted S= 703.9 N= 1489.1 dN/dS= 0.2247 dN = 0.0344 dS = 0.1529

Could anybody with codeml experience explain to me why the parser fails to yield most of the parameters (values) I am interested in?

Thank you

Answer: Hi,

You can obtain omega values for each branch using the following commands:

from Bio.Phylo.PAML import codeml

results = codeml.read(paml_outputfile)

print results["NSsites"][0]["parameters"]["omega"]

This gives you a list of omega (dn/ds) for each branch

I am using version 4.7b and biopython-1.6.3

Hope that helps!

A
Question: Hi Everyone,

I'm trying to find orthologs and lineage specific paralogs between two species. I tried using the ensembl homology pipeline but both my species are not on the database. Therefore, I tried to write my own similar pipeline. So far I've accomplished the following:

1. Blast all for every gene in both genomes
2. Filtering of blast results based of evalue and alignment length
3. Single Linkage clustering with MCL to form gene families
4. For each gene family, I did a protein alignment with PRANK and built a tree with Treebest, which also takes in the species trees and tries to build a gene tree accordingly.

My question deals with parsing these gene trees. I want to use these gene trees to find paralogs and orthologs between my two species but I'm not sure how to parse all of the topologies of these trees and how to determine paralogous or orthologous relationships.

Are there any programs that can take in gene trees and output a list of paralogs and orthologs?

Thanks.

LC

Answer: The ETE toolkit is indeed capable of doing that. You can use a species overlap algorithm to detect duplication and speciation events, or reconcile your gene tree with the expected species trees (this paper includes a comparison of both methods).

Briefly, you will need to load your gene tree as a PhyloTree object, and then call any of the tree.get_descendant_evol_events or tree.reconcile methods. Both methods will process your tree, label the nodes as speciation or duplication and return a list of speciation and duplication events. Then you can visualize your tree with `tree.show()` to see the predictions or process the list of events. If you used `tree.reconcile`, a reconciled-tree PhyloTree object will also be returned, including the inferred lost branches (also visible con tree.show). There are some examples in the ETE tutorial showing how to get orthology/paralogy prediction based on gene trees: http://pythonhosted.org/ete2/tutorial/tutorial_phylogeny.html#detecting-evolutionary-events
Question: Hi,

Has anyone determined how signal-to-noise changes in MNase-Seq profiles with 100 vs 200 vs 400 million reads? Does it improve? Or does it stay the same regardless?

Thanks!

Answer: This [method paper][1] for human MNase-seq recommends >10x coverage.
In [this other paper][2] on yeast, they reach between 13x and 177x fold coverage and consider it sufficient.

So I would say that at least 10x coverage is enough for proper peak definition.

[1]: https://www-ncbi-nlm-nih-gov.insb.bib.cnrs.fr/pmc/articles/PMC3541821/
[2]: https://academic.oup.com/nar/article/42/8/4934/1076121/Fuzziness-and-noise-in-nucleosomal-architecture#55493469
Question: Hi Guys,

I am using pileup in Rsamtool from R package. I need to get the total "count" for each position. For example, for chr11, position 1643082, there are a total of 1 count of A, 3 counts of C and 18 counts of T. I want them in this order:

chr11:1643082 A:1, C:3, T:18

and so forth for position chr11, 47663948 to be shown in the result. There are thousands of different positions and I need to calculate the counts for each nucleotide.Thanks!

`table1`

```r
seqnames pos strand nucleotide count which_label
1 chr11 1643082 + A 1 chr11:1643082-1643082
2 chr11 1643082 - C 3 chr11:1643082-1643082
3 chr11 1643082 + T 15 chr11:1643082-1643082
4 chr11 1643082 - T 3 chr11:1643082-1643082
5 chr11 47663948 + C 16 chr11:47663948-47663948
6 chr11 47663948 - C 11 chr11:47663948-47663948
7 chr11 47663948 + T 2 chr11:47663948-47663948
```

Answer: The alternative using dplyr:

```r
library(dplyr)
df %>%
group_by(seqnames,pos,nucleotide) %>%
summarise(count=sum(count))
```

Edit: to get the exact output you described, you have to:

```r
> biost %>%
group_by(seqnames,pos,nucleotide) %>%
summarise(count=sum(count))%>%
mutate(mypos=paste(seqnames, pos, sep=':')) %>%
group_by(mypos) %>%
arrange(nucleotide) %>% # Make sure nucleotides are always sorted the same way
summarise(counts=paste(nucleotide, count, sep=':',collapse=', '))

Source: local data frame [2 x 2]

mypos counts
1 chr11:1643082 A:1, C:3, T:18
2 chr11:47663948 C:27, T:2
```
Question: Hi Biostars Users,

I am rookie in bio-informatics interested in undergoing NGS training to learn how to build NGS pipeline for whole exome sequencing analyses.

May you please guide me to find a suitable training program?

Best,
Kizuna

Answer: Some links.

- https://www.biostars.org/p/1268/
- http://www.cureffi.org/2012/09/19/exome-sequencing-pipeline-using-gatk/
- http://seqanswers.com/wiki/How-to/exome_analysis
- http://www.gqinnovationcenter.com/services/bioinformatics/tutorials/tutorials.aspx?l=e
Question: Hello, I have RNAseq data with a major perturbation where we expect a large number of genes (>50% of the dataset) to be differentially expressed. I'm aware that this can be a problem for TMM normalization which assumes that the majority of genes are not DE between samples. However, I haven't seen any other approaches suggested for when this assumption isn't met.

I suspect it may be necessary in our data because, in addition to the large perturbation, we observe systematic differences in the normalization factors between conditions that are independent of library size. This leads to sign-changes where the cpm data show genes as down-regulated but the normalized data show the the genes as up-regulated. I understand this could be produced by outlier genes (what TMM is trying to correct for), but given the genome-wide differences I'm unsure if applying the correction here is appropriate.

Any advice or suggestions regarding a better way to normalize datasets with a major perturbation?

Answer: You can always inspect the MA-plot to judge how normalization goes. If you have a set of control genes that you can confidently identify as probably non-DE then you can just run the TMM calculation on these and then feed this back to the full dataset. After all the edgeR procedure is to calculate CPM based on library size and then correct this for composition with these TMMs. Or you can try to use only genes with large logCPMs as these usually tend to be rather non-DE. Can you show an MA-plot?
Question: Hi,

I have to plot allele frequencies of two different SNP chip datasets. I have two VCF files and would like to make a scatterplot in which these 2 datasets are plotted one against each other. What is the easiest way to do this?
I apologize if this question was already posted.

Answer: Hi, your question motivated me to add a new Python API method called `pyvcf.plot_af_correlation` ([documentation](https://sbslee-fuc.readthedocs.io/en/0.28.0-dev/api.html?highlight=plot_af_correlation#fuc.api.pyvcf.plot_af_correlation)) to the `fuc` [package](https://github.com/sbslee/fuc) I wrote. Take a look at the example below and see if this achieves your goal.

![enter image description here][1]

```
>>> from fuc import pyvcf, common
>>> import matplotlib.pyplot as plt
>>> data1 = {
... 'CHROM': ['chr1', 'chr1', 'chr1', 'chr1', 'chr1', 'chr1'],
... 'POS': [100, 101, 102, 103, 104, 105],
... 'ID': ['.', '.', '.', '.', '.', '.'],
... 'REF': ['G', 'T', 'G', 'T', 'A', 'C'],
... 'ALT': ['A', 'C', 'C', 'G,A', 'C', 'T'],
... 'QUAL': ['.', '.', '.', '.', '.', '.'],
... 'FILTER': ['.', '.', '.', '.', '.', '.'],
... 'INFO': ['.', '.', '.', '.', '.', '.'],
... 'FORMAT': ['GT:DP', 'GT', 'GT', 'GT', 'GT', 'GT'],
... 'A': ['0/1:30', '0/0', '1/1', '0/1', '1/1', '0/1'],
... 'B': ['0/0:30', '0/0', '0/1', '0/1', '1/1', '0/1'],
... 'C': ['1/1:30', '0/0', '1/1', '0/1', '1/1', '0/1'],
... 'D': ['0/0:30', '0/0', '0/0', '0/0', '1/1', '0/1'],
... 'E': ['0/0:30', '0/0', '0/0', '1/2', '1/1', '0/1'],
... }
>>> vf1 = pyvcf.VcfFrame.from_dict([], data1)
>>> data2 = {
... 'CHROM': ['chr1', 'chr1', 'chr1', 'chr1', 'chr1'],
... 'POS': [101, 102, 103, 104, 105],
... 'ID': ['.', '.', '.', '.', '.'],
... 'REF': ['T', 'G', 'T', 'A', 'C'],
... 'ALT': ['C', 'C', 'G,A', 'C', 'T'],
... 'QUAL': ['.', '.', '.', '.', '.'],
... 'FILTER': ['.', '.', '.', '.', '.'],
... 'INFO': ['.', '.', '.', '.', '.'],
... 'FORMAT': ['GT', 'GT', 'GT', 'GT', 'GT'],
... 'F': ['0/0', '0/1', '0/1', '1/1', '0/0'],
... 'G': ['0/0', '0/1', '0/1', '1/1', './.'],
... 'H': ['0/0', '0/1', '0/1', '1/1', '1/1'],
... 'I': ['0/0', '0/1', '0/0', '1/1', '1/1'],
... 'J': ['0/0', '0/1', '1/2', '1/1', '0/1'],
... }
>>> vf2 = pyvcf.VcfFrame.from_dict([], data2)
>>> pyvcf.plot_af_correlation(vf1, vf2)
>>> plt.tight_layout()
```

Please note that since this method is still being implemented in a development branch (`0.28.0-dev`) of fuc, in order to use it you need to install the `fuc` package locally:

```
$ git clone https://github.com/sbslee/fuc
$ cd fuc
$ git checkout 0.28.0-dev
$ pip install .
```

Please let me know if you have any questions and/or suggestions regarding this.

[1]: /media/images/f1fc4877-479f-45ac-9309-c905c5ee
Question: I have a few hundred yeast sequences (20-80bp long) and I want to find common motifs (conserved bases at certain indices) in them. I am using a Mac

Answer: I would recommend MEME as others, however Weeder is also very easy to install and run on a Mac. It is word based and runs a lot faster than Meme, with a far fewer options.

If you would like something to run specifically on your Mac you could try iMotifs, which incorporates the NestedMICA algorithm for discovering over-represented motifs.

For a small set of sequences, such as yours, you could also use Meme, Weeder or others online without installing them locally. Finally, the Regulatory Sequence Analysis Tools (RSAT) ( http://rsat.ulb.ac.be/rsat/ or http://rsat.ccb.sickkids.ca/) are a great place to start for pattern matching and discovery.

In case it is of any use I have put up a short presentation about DNA motif finding from a course I gave earlier this year. You can find it here.
Question: Hi all,

I have RNAseq data of non-model organisms, which I'm planning to analyze by de novo assembly with Trinity and write a paper on.

The RNAseq paper contains basic statistics(for example, Total raw reads, Number of gene and Assembly N50 etc).
Are there any rules for these parameters?
I'm a beginner in bioinformatics, so I don't know what parameters are required.
Can you please advise me?

Thank you very much for your much appreciated help on this!

Answer: I think there are no fixed rules for these numbers that would determine if you can or cannot publish. Transcriptomes are easier to assemble than genomes due to less repetitive regions, so, N50 is limited by the actual lengths of transcripts. If your organism has different developmental stages and tissues you should rather devise a good sampling strategy to distribute your funds in order to detect tissue and stage-specific transcripts, or make a good pool of everything. I would say, this is more important than just how many paired-end reads you can afford. I think you can already get something decent from >50M fragments per library which is at the lower end of what you get today.

Wrt. the number of transcripts after assembly, how many transcripts including isoforms do you expect from related species? If you get an absurd number of unigenes (100k+ genes) then you might have to sequence deeper.

Here is an example of such a paper in Parasites & Vectors: https://parasitesandvectors.biomedcentral.com/articles/10.1186/s13071-020-04442-2#Sec9 (I noticed the authors report rather low % of database hits)

They have ~460M reads in total, that would mean a cost/effect-ratio of 1 Illumina 4000 lane per publication, but I think one can go with less in some journals, and of course you need to analyze the data properly and write something interesting after that.
Question: I am using 1000 Genomes data with my new project.
When I am inspecting the reference assembly they have been using, I found it contains a "decoy" contig.

The 1000 Genomes FAQ says:

For the final round of alignments the sequence data will be mapped to
a set of sequences derived from the GRCh37 assembly. This
GRCh37-derived alignment set includes chromosomal plus unlocalized and
unplaced contigs, the rCRS mitochondrial sequence (AC:NC_012920),
Human herpesvirus 4 type 1 (AC:NC_007605) and decoy sequence derived
from HuRef, Human Bac and Fosmid clones and NA12878. These files are
available in phase2_reference_assembly_sequence on the ftp site. All
human variant coordinates reported by the 1000 Genomes project are in
GRCh37 coordinates.

Here, I have no idea what the decoy sequence is and why it is included. Maybe to detect sample contamination?

Answer: The reference genome is incomplete, particularly around the centromeres, so often reads which truly belong elsewhere are wrongly mapped to a particular place in the genome because the true match is missing from the reference. These cause false positive calls, which were bothering us in the 1000 Genomes Project. The decoy is a pragmatic solution to this - it contains known true human genome sequence that is not in the reference genome, and will "suck up" reads that would otherwise map with low quality in the reference. The decoy was built by Heng Li, at the Broad, working with Richard Durbin (Sanger) and Deanna Church of the Genome Reference Consortium (who maintain the reference genome).
Question: Hi everybody, I want to use GCTA to estimate variance components. My data is for horse species.

This program doesn't recognize the command:

> gcta64 --reml --grm r --pheno remel.phen.txt --reml-pred-rand
> –qcovar tt.qcovar.txt --out test1

*The following error occurs: Error (invalid option "--qcovr.txt")*

my format of tt.qcovar.txt : (no header line; columns are family ID, individual ID, sex, EBV)

How can I solve my problem?

I appreciate your guiding me regarding this.
Thanks

Answer: Hi, it is likely due to the fact that you are using an em-dash, as follows:

–qcovar

It should, however, be two hyphens, like this:

--qcovar

Kevin
Question: I have tried to count mapped reads on the human miRNA regions from the miRBAse human GFF

the output was just:

__no_feature 70216353
__ambiguous 0
__too_low_aQual 0
__not_aligned 1383340
__alignment_not_unique 0

So there should be something I do wrong?

I have run it - htseq-count SAMfile Gfffile (SAMfile is from bowtie with option -S)

The GFF is from MirBase:

chr1 . miRNA_primary_transcript 17369 17436 . - . ID=MI0022705;Alias=MI0022705;Name=hsa-mir-6859-1
chr1 . miRNA 17409 17431 . - . ID=MIMAT0027618;Alias=MIMAT0027618;Name=hsa-miR-6859-5p;Derives_from=MI0022705
chr1 . miRNA 17369 17391 . - . ID=MIMAT0027619;Alias=MIMAT0027619;Name=hsa-miR-6859-3p;Derives_from=MI0022705
chr1 . miRNA_primary_transcript 30366 30503 . + . ID=MI0006363;Alias=MI0006363;Name=hsa-mir-1302-2
chr1 . miRNA 30438 30458 . + . ID=MIMAT0005890;Alias=MIMAT0005890;Name=hsa-miR-1302;Derives_from=MI0006363
chr1 . miRNA_primary_transcript 187891 187958 . - . ID=MI0026420;Alias=MI0026420;Name=hsa-mir-6859-2
chr1 . miRNA 187931 187953 . - . ID=MIMAT0027618_1;Alias=MIMAT0027618;Name=hsa-miR-6859-5p;Derives_from=MI0026420
chr1 . miRNA 187891 187913 . - . ID=MIMAT0027619_1;Alias=MIMAT0027619;Name=hsa-miR-6859-3p;Derives_from=MI0026420
chr1 . miRNA_primary_transcript 632325 632413 . - . ID=MI0022558;Alias=MI0022558;Name=hsa-mir-6723
chr1 . miRNA 632382 632403 . - . ID=MIMAT0025855;Alias=MIMAT0025855;Name=hsa-miR-6723-5p;Derives_from=MI0022558

Answer: Did you use the `-t` flag to change the feature type to miRNA_primary_transcript or miRNA? As far as HTSeq would be concerned (without using this flag) you have given it no features to look for.

The default looks for exon in the third column

You should also use the `-I` flag to use 'ID' as the id attribute (default looks in the last column for 'gene_id')

Ensure you read the manuals before using these tools.

http://www-huber.embl.de/users/anders/HTSeq/doc/count.html
Question: Efficient method(s) for comparison of data sets that are generated from different microarray experiments (meta-analysis) have been asked [previously][1]. Based on the discussion, I found [RankProd][2] R package. However, after using the method, I wanted to compare the results with an another method since RankProd was developed years back. I wold like to know any alternative method available for integrating multiple data set(RNASeq not included).

Out of 10 data sets, 8 are from Affymetrix and 2 are Agilent. Expression profiles are based on different developmental stages of an organ in Arabidopsis.

[1]: https://www.biostars.org/p/4418/
[2]: https://bioconductor.org/packages/release/bioc/html/RankProd.html

Answer: **Updated February 29, 2020**

**Updated October 9, 2018**

-----------------------------------------
---------------------------

#If the array versions are the same
If the Affymetrix arrays are all of the same version/type, then just process them all together at the same time. Same for the Agilent arrays. Afterwards, use one as the training dataset and the other as the validation dataset, i.e., obtain results from one, and then corroborate these [results] in the other.

You could attempt to perform statistical analyses on them combined, too (i.e. not as training and validation), but, in this case, include `ArrayVersion` (i.e. batch) as a covariate in downstream analyses - there will likely be some batch effect.

#If all data-set arrays are different

- Process each *independently*
- Summarise expression values across whole genes *independently*
- Create a list of common features/genes across all arrays
- Convert the data to Z-scores *independently*
- Merge (bind) the data together
- Build regression models predicting for your outcome and include
'`ArrayVersion`' as a covariate in the model

My advice is to avoid the use of programs that aim to directly correct for batch. It is preferable to include batch as a covarate in your statistical models.

---------------
-----

**NB - a word of caution about these approaches:** these are not 'fool proof'. If you are trying to analyse, e.g., brain tissue and skin tissue together, then the expression profiles will be so different such that you could not analyse them together in any way. Proceed with caution with these approaches that I mention, and take time to analyse the output via histograms and box-and-whisker plots.
Question: Hi,

I am wondering if any one can recommend a simple yet applicable tutorial for learning bash to be used for doing bioinformatics. The reason I am asking here is that there are so many tutorials out there and to my opinion some of them are too much and some are not targeting my needs, since many are written for developers and advanced users. So please let me know if you know good ones!

Cheers,
Parham

Answer: As a [Software Carpentry][1] (and Data Carpentry) instructor, I would recommend the [Introduction to the Shell tutorial materials][2] in the [lesson material][3]. Once you are comfortable with those tutorials, [Data Carpentry is developing Shell tutorials specifically for bioinformatics][4].

There are a ton of other resources out there for more specific questions... Best of luck!

[1]: http://software-carpentry.org/index.html
[2]: http://swcarpentry.github.io/shell-novice/
[3]: http://software-carpentry.org/lessons.html
[4]: https://github.com/datacarpentry
Question: Greetings all,

Let me start by saying I know how to add read groups to a BAM file. What I am looking for is the easiest way to do it. This problem has been around for a while so I though I would see what neat tricks people are using.

Thanks

Answer: In my opinion the easy way to do is using *AddOrReplaceReadGroups* in Picard:

- http://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups
- https://gatk.broadinstitute.org/hc/en-us/articles/360037226472-AddOrReplaceReadGroups
Question: Hi all,

I was recently given a data set (run through Illumina in 2014) to merge with a set run through Illumina this year. The set I just received was not demultiplexed and was dual barcoded. I know that Illumina's bcl2fastq can handle dual barcoded sets and demultiplex, but I don't have any of the BaseCall data it typically uses to demultiplex. I only have .fastq.gz files and a mapping file to work with. Does anyone know if I can still use the bcl2fastq without the BaseCall data, or if there is something that would work better?

tldr; I need to demultiplex dual barcoded .fastq.gz files

Answer: You can't use bcl2fastq to demux standalone fastq.gz files. You can take a look at this thread (https://www.biostars.org/p/270737/ ) to get inspiration on how to do the demultiplexing of your dataset.

Edit: If you have in-line index sequences then a different solution would be needed. If you have standard Illumina datafiles (post what you get from @h.mon's command) then `demuxbyname.sh` from BBMap will work as indicated in the thread in my answer above and the original solution referenced there.
Question: Hi,
I'm trying to get NCBI vecscreen (http://www.ncbi.nlm.nih.gov/VecScreen/VecScreen.html) to work as a command-line for examining many sequences. I've tried compiling NCBI toolbox, but NCBI toolbox is very large and didn't work for me very well. I tried compiling with only blast, but the dependencies for vecscreen don't seem to be included when I tried to run vecscreen. I tried using qblast (with blastn) in BioPython with the parameters for vecscreen on the website, but I got different results. I tried compiling NCBI toolbox on linux and windows. I couldn't find how to compile blast, etc. in Visual C++ on Windows, and the linux attempt didn't have the dependencies or something for vecscreen.

Does anyone know:

Can vecscreen be used in biopython somehow? Maybe with qblast?

What configuration options, etc. do I need to run for the NCBI toolbox to get vecscreen to compile on windows or linux?

The BioPython qblast parameters that I used are:
result_handle = NCBIWWW.qblast("blastn", "UniVec", record.seq, expect=700, filter="m D", nucl_penalty=-5, gapcosts="2 2", searchsp_eff=1.75e12)
VecScreen says that gapcosts should be "3 3", but qblast doesn't accept this. I got different results than running vecscreen from the website.

The VecScreen parameters are pre-set using blastn options:
-q -5 -G 3 -E 3 -F "m D" -e 700 -Y 1.75e12

Any help is appreciated.

Answer: You can download stand-alone-blast and the UniVec databases and then blast against UniVec locally. This will allow you to blast many sequences at once.