cannin/biostars_qa huggingface dataset:
Question: Does anyone have a simple solution to downloading all the refseq genomes for a particular taxon?

Using `ncbi-genome-download` its possible to specify the species or genus TaxIDs and download them, but apparently you can't go higher up the taxonomic ranks (even though enterobacteria has a TaxID of 543 for instance).

If anyone knows of a way to download all the Enterobacteria I'm all ears.

Alternatively, if there is a method of extracting the species TaxIDs from the Enterobacterial taxid in NCBI such that I can pass them all directly to `ncbi-genome-download` that would work too.

Answer: From speaking with a few other pros, this was the solution in the end (though only very rough at the mo):

Use the `ete3` toolkit to get a list of IDs:

from ete3 import NCBITaxa
import sys

taxon_name = sys.argv[1]

ncbi = NCBITaxa()
ncbi.update_taxonomy_database()
ebact = ncbi.get_descendant_taxa(taxon_name)

with open('./taxids', 'w') as ofh:
for i in ebact:
ofh.write("%s\n" % i)

# At this point, one could import ncbi-genome-download as a python method and continue

Which gave me a list of IDs (though this includes ALL descendent taxa, even ones without complete genomes etc).

I passed these to the latest version of `ncbi-genome-download` which accepts a `--taxid 12345,65890` format for specifiying the IDs.

So I just ran:

for file in * ;
do python ~/bin/ncbi-genome-download/ncbi-genome-download-runner.py -l complete -v -p 10 --taxid $(paste -s -d ',' "$file") bacteria ;
done

I had to run this iteratively on many files after I split my `taxids` file up as there is a limit to how many args can be passed to `--taxid` at once.

EDIT Sept 2018:

I contributed a script to the `ncbi-genome-download` repo to make getting the TaxIDs nice and easy. It uses the approach above, but there’s no need to rewrite it for oneself now.
Question: Hi all,
I'm looking for a simple solution for renaming fasta headers.

I have this fasta header

>trpE___AA_HMM___6fa05435949258489b608db9e58e5ba38821f2f26fffe5755daff43abin_id:MALBOS1|source:AA_HMM|e_value:5.2e99|contig:MALBOS1_000000117228|gene_callers_id:113772|start:215745|stop:217260|length:1515

And I would like to rename it only like this

>MALBOS1_000000117228

That means, remove everything before the pattern "contig:" and after "|gene_callers_id"

Any ideas?

thanks

Answer: perl -pe 's/>.+contig:(.+?)\|.+/>$1/' < FASTA_IN > FASTA_OUT
Question: Hi, all.

Can anybody tell help me explain which parent comes first when writing the genotype?

if phased,

mom 0/1 dad 0/0
child 1|0 or 0|1?

many thanks.

Answer: I am not working with trios myself - so I am just thinking out loud.

Since the child may have novel variants and the parents may also share variants, could we always correctly assign a variant to either parent? Many times would possible, but other times not.

I believe that the order `0|1` vs `1|0` has meaning only in terms of phasing. Every variant in a block lists the genotypes to be in sync with other variants in that block. That is the sole purpose of the order. For unphased variants, the order is always low to high `0/1`

To figure out which parent shares the variant look at the parent genotype.
Question: I have gene lists that I want to compare based on how well they were conserved across human and mice. Is there an established way of doing this? My first impulse was to compare identity directly but since genes have differing lengths and differing functional sizes of differing importance, a measure just based on sequence identity seems like a bad way of comparing any two genes to each other.

Answer: Hi oganm, you can do this the complicated way or the easy way. The easy way is with identities (from Ensembl BioMart for instance) and the complicated way could be with a BLASTP similarity network.

Yes orthologs have differing length but in a pairwise comparison (like human-mouse), you have limited info.
Question: I am trying to do exactly the same pipeline for RNA-seq data process as the TCGA does. Usually when we ask a sequencing service, we can get a fastQ file. It contains sequence and read quality information. The alignment step comes next.
However in the case of TCGA, as their pipeline suggested ([https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/][1]). It seems that they used BAM files as one of the inputs. I was wondering why they used BAM files as inputs and how can I repeat what they did?
In addition, why isn't there seem to be a adaptor trimming process?

[1]: ![TCGA mRNA-seq pipeline schematic][1]https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/

[1]: https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/images/gene-expression-quantification-pipeline-v3.png

Answer: https://gdc.cancer.gov/about-gdc/gdc-faqs - Read the answer to "How can I access GDC sequencing data in FASTQ format?"

Level 1 (aka the raw fastq) data is restricted; to request access to it, see instructions here: https://gdc.cancer.gov/access-data/obtaining-access-controlled-data

Adaptor trimming is unnecessary in RNA-seq read mapping; many papers have been written about it (e.g. https://academic.oup.com/nargab/article/2/3/lqaa068/5901066 ).
Question: Dear All,

I have to deal some small molecules screening for my protein-drug docking studies. Please suggest me few small molecule databases ..also the free software can use to screen them...thank you

Answer: What are you actually trying to achieve? Your question is extremely vague.

Do you just want a big collection of drug-like molecule? Grab GDB-17

Do you want a curated source of biologically interesting molecules? Go to ChEBI

Do you want small molecules with annotated activity? Try ChEMBL

Do you want small molecules bound to a protein? Try PDB

For docking, just go to any literature search engine and search for small molecule docking. You will get reviews, algorithms and papers describing software. There is plenty of literature, so start reading.

Edit: After several comments, it seems the interest lies in identifying any compound which may inhibit the activity of the given receptor. For this purpose I would suggest using ZINC, which also has vendor information.
Question: I am trying to convert BAM files to Bigwig format. I was originally using a combination of genomeCoverageBed and BedgraphToBigWig to go from BAM to Bigwig in 2 steps:

samtools sort file.bam -o file.sorted.bam;
genomeCoverageBed -bg -split -ibam file.sorted.bam -g file.chrom.sizes > file.bedgraph;
sort -k1,1 -k2,2n file.bedgraph >file.sorted.bedgraph;
bedGraphToBigWig file.sorted.bedgraph file.chrom.sizes file.bw;

After finding out that deeptools (function called bamCoverage) could convert BAM to Bigwig AND normalize all in one go, I decided to use it. I first decided to not normalize (just for testing purposes)

bamCoverage -b file.bam -o file.bw

This command works but I can't help but notice the output file.bw in both scenarios are drastically different. My starting BAM file is 1.1GB. When using bedGraphToBigWig, my output file is 122MB. When using bamCoverage, it's 23MB.

I know that you can change the --binSize in bamCoverage which will lead to a larger file. I tried a bin of 10 and my file was 45MB.

My question is: what is the relationship between file sizes between the 2 software? I am confused as to why the size is so different. From what I understand, in bamCoverage, the coverage is calculated as the number of reads per bin, where bins are short consecutive counting windows of a defined size. But what is the equivalent setting in BedGraphtoBigwig? Is it -blockSize or maybe itemsPerSlot? Their defaults are 256 and 1024 respectively so I am a bit confused.

I am not sure whether this is something to worry about? Presumably, all it means is that my 2 files (the one produced with BedgraphToBigwig and the one produced with bamCoverage) will have different 'resolutions', when visualized on a something like IGV or UCSC Genome Browser. Is this correct?

Answer: The reason is that `genomeCoverageBed` by default spills out bedGraphs at the base-pair level, so it piles up the depth for every base. In contrast, `bamCoverage` has a default 50bp window (option `-bs`, see documentation) that it aggregates reads over. I always use `-bs 1` because visually these tracks looks much nicer/smoother than the default `-bs 50`.
Question: Hi. I just have a quick question in regards to calculating the N50 for a list of contigs. Anybody know how to do it in excel, by any chance? I'm not much of a programmer just yet, and I was just wondering if there was an easier way than just constructing a python program or something.

Thanks!

Answer: Get the lengths of your contigs into a column, and sort that column. Sum over all lengths to get the total number bases N, and then just traverse your list of lengths from top to bottom until the cumulative sum reaches N/2. (You could use an additional column to keep track of the running sum if that makes it easier.)

The length of the contig where your running sum passes N/2 is the N50. (In the case where the running sum exactly equals N/2, then the N50 is the average of this length and the subsequent one in the list.)

For example, below, the cumulative sum exceeds 44/2=22 on a contig with length 6, so N50 = 6.

contigs length cumulative sum
AAAAAAA 7 7
AAAAAA 6 13
AAAAAA 6 19
AAAAAA 6 25
AAAA 4 29
AAAA 4 33
AAA 3 36
AAA 3 39
AA 2 41
AA 2 43
A 1 44
Question: Dear all,

referring to the batch correction methods for scRNA-seq, would you have any preference and/or comments ? among possible choices :

-- MNNCorrect, as outlined in SimpleSingleCell workflows :

https://bioconductor.org/packages/release/workflows/html/simpleSingleCell.html

-- ZINB-WAVE :

https://bioconductor.org/packages/release/bioc/html/zinbwave.html

-- HARMONY :

https://www.biorxiv.org/content/10.1101/461954v2

-- SCTransform :

https://satijalab.org/seurat/v3.0/integration.html

thanks a lot,

bogdan

Answer: The results seem to be very experiment-specific. For example, in today's [SCRIBE pre-print][1], all the methods (except the one introduced) perform poorly:

![enter image description here][2]

One thing to notice is that they all fail in different ways, so the problems don't seem to be due to some artifact in the data itself. For example, MNN mixes NF and TH, but Seurat splits PEP.

[1]: https://www.biorxiv.org/content/10.1101/793463v2
[2]: https://user-images.githubusercontent.com/6363505/66263486-70484680-e7c1-11e9-99ee-a0cb9bcd05db.png
Question: I have a pipeline script called `pipeline.sh` . I usually execute this for a single sample like so:

$ pipeline.sh sample1 sample1.R1.fastq.gz sample1.R2.fastq.gz

Where $1 is sample ID, $2 and $3 are the read files.

I use GNU parallel with a parameters file that specifies the paths to each file.

$ nohup parallel -j 4 -a params.pipeline.txt --colsep '\s+\ ./pipeline.sh

I want to automate my pipeline so that I don't need a parameters file and it looks through the folders for the reads (same file structure as they come out of the sequencer).

I have:

parallel -j 4 ./pipeline.sh {/1.} ::: *.R1.fastq.gz :::+ *.R2.fastqgz

However this assumes the fastqs are in the same folder, how can I change my parallel so that it searches through the folder structure for the files.

Answer: using nextflow (not tested, but it should look like this):

nextflow.enable.dsl = 1
params.directories=""

process scanDirectories {
output:
path("paths.txt") into paths
script:
"""
find ${params.directories} -type f -name "*.R1.fq.gz" \
awk -F '/' '{S=\$NF;gsub("\\.R1\\.fq\\.gz\$","",S);F2=\$0;gsub("\\.R1\\.fq\\.gz\$",".R2.fq.gz",F2);printf("%s,%s,%s\\n",S,\$0,F2);}' > paths.txt

"""
}

paths.splitCsv(header: false,sep:',',strip:true).set{pipe_in}

process runPipeline {
tag "${sample}"
input:
tuple val(sample),val(R1),val(R2) from pipe_in
output:
path("result.txt") into result_ch
script:
"""
echo "DO Something ${sample} ${R1} ${R2}" > result.txt
"""
}

and then something like

nextflow run -resume script.nf --directories "/path/to/dir1 /path/to/dir2"
Question: I designed a bait set to caputure genes but i would like to know whether my baits bind on other places other than the designed gene sequences. My idea is to blast it back to the reference genome and identify the "background noise" of the non gene sequences. An easy solution i think is to remove all gene sequences from my reference genome. Leaving only the intergenic and repetitive regions.

Does someone have a easy solution for this?
Or other suggestion i could try?

Answer: get a BED for your genes and then use bedtools mask fasta http://bedtools.readthedocs.io/en/latest/content/tools/maskfasta.html to replace the bases with N.
Question: Does Pindel generate genotypes or does it only report on sites? Adding extra characters to get to 80.

Answer: There is a `pindel2vcf` binary inside Pindel package to convert Pindel calls to genotyped vcf.
Question: I am trying to use the Bingo Plugin in Cytoscape but the output are empty and I can not figure out why. Can you guys give me some enlightment? Here is a screenshot of my input.

More one question: is 'Reference set' the background of genes?

Thank you guys.

![< image not found >][1]

[1]: http://imageshack.com/a/img661/9255/aHtdnw.png

Answer: I've not used it for a while, but be sure your Identifiers(Ensemble I guess in this case) can be mapped to the ontology files. It is better to update the ontology and annotation files in Bingo before using it.
Question: Hi all,

I'm looking to analyze mouse ERVs from bulk-RNAseq but have mostly found data resources and pipelines for human ([ERVmap][1], [HERVd][2]). I've also found a mouse [ERE database][3] that appears to no longer be supported.

The only reasonable resource I've found for mouse ERVs is [gEVE][4] but wanted to ask if anyone is aware of additional resources?

[1]: https://www.pnas.org/content/115/50/12565
[2]: https://herv.img.cas.cz/
[3]: https://www.ncbi.nlm.nih.gov/pubmed/22691267
[4]: https://academic.oup.com/database/article/doi/10.1093/database/baw087/2630466

Answer: Can you be more specific regarding what resources you're after?

You might have better luck with gEVE:
https://academic.oup.com/database/article/doi/10.1093/database/baw087/2630466
Question: This post is related to, but different (a continuation) from my earlier post at https://www.biostars.org/p/268060/ seeking help with running the genome assembler software - ABySS. I got my ABySS runs going.

But in order to assess the optimal k-mer, I was trying values 21-97 (increments of 2).

Related to this attempt, at k-mer value>64 (i.e. 65), the STDERR log terminated with the line - "ABYSS: ../Common/Kmer.h:48: static void Kmer::setLength(unsigned int): Assertion `length <= 64' failed."

Obviously my university's compute cluster has ABySS compiled and configured to allow k-mers <= 64.

Some context for my questions below - I have 2*150 PE Illumina HiSeq400 reads from fungal spore DNA (haploid genome). Read length trimming has resulted in size distribution between 41 (min) and 150 (max).

1a. How do I execute ABySS with k-mer value > 64, but WITHOUT re-compiling? Is that possible? I ask in the context of what is written at https://insidedna.me/tool_page_assets/pdf_manual/abyss.pdf, which reads:

> ./configure --enable-maxk=96
>
> "The default maximum k-mer size is 64 and may be decreased to reduce
> memory usage or increased at compile time. This value must be a
> multiple of 32 (i.e. 32, 64, 96, 128, etc):"

1b. If running k-mer > 64 is NOT possible without recompiling, is there a work-around that allows re-setting local variable(s) on just my compute cluster account?

2. Is there a theoretical / practical limit on k-mer value > default 64? I know SPAdes, for example, permits as high as 127. What is it for ABySS? Is it 128?

3. Is checking k-mer value > 64 **necessary**? Which is a different question from whether it is **possible**?

Thank you!

Answer: Hi Anand,

Glad to hear you were able to get some assemblies running!

To answer your questions:

* Unfortunately, there is no way to run an assembly with k-mer size > 64 without re-compiling the software.
* There is no upper bound for `--enable-maxk`. However: (i) increasing `--enable-maxk` increases the memory requirements of ABySS, and (ii) you will have to adjust some MPI settings for large k-mer sizes, as described here: https://github.com/bcgsc/abyss/wiki/ABySS-Users-FAQ#2-my-abyss-assembly-jobs-hang-when-i-run-them-with-high-k-values-eg-k250
* Yes, it is very likely that you will need to use k-mer sizes > 64 in order to achieve the best possible assembly of your data

If you are uncomfortable with compiling software from source, you can perhaps get your IT team to do it for you (if you have an IT team). That said, knowing how to compile software from source is a very useful skill in the bioinformatics domain. If you are interested in learning, here is an introductory tutorial that may be helpful: http://www.thegeekstuff.com/2012/06/install-from-source/
Question: I've been working on editing some code on the server at my work, I've come into a problem where we need a snp quality to go through our filters downstream in the pipeline that we use. Does anyone know if Snp Quality exists in the VCF file? If so where?

Answer: SNP Quality can be represented in several places in a VCF file.

1) the QUAL column, which is the phred-scaled quality score for the assertion made in ALT. In other words, it's: 10log_10 prob(call in ALT is wrong).

2) GQ, encoded in the FORMAT column is genotype quality, encoded as a phred score: -10log_10p(genotype call is wrong).

3) Especially if you're looking at tumor/normal pairs, you may see that it's represented as VAQ (variant quality).

The header of the VCF should give a description of exactly which fields are present in your files and help you determine which ones contain the quality scores that you're looking for.

For more info, check out the VCF format description
Question: I just want to make sure my interpretation is correct: a MAPQ value of 255 indicates that the mapping quality is "not available" because it's a unique alignment, and a unique alignment means that it has a 0 probability of mapping to the wrong place, therefore: -10*log10(0 probability of error) = inf (ie. "not available")

Answer: tldr: MAPQ of 255 has little actual meaning.

Whether 255 actually means that an alignment is unique or "I don't want to bother returning a meaningful value" is aligner dependent. In practice, tools that produce MAPQ of 255 are doing so because they aren't calculating a real MAPQ and are just indicating, "this is a more likely alignment than anything else I could find". How much more likely this particular alignment is is purely a matter of the settings you use, which are pretty arbitrary. As a rule of thumb, assume that MAPQ values from splice-aware aligners lack any standard phred-scale interpretation.
Question: I'm very new to bioinformatics so any suggestion how to start ?

my field is microbiology & molecular biology

Answer: Hi!,

I have been there! Do you know any computer programming or scripting language! If not, start with a Linux or awk tutorial. There are several very user friendly bioinformatics tools. Install them into your computer, choose a project and try to analyze some sample data. Get some NGS data try to analyze it. If you do not have it yourself, try to download from publications or sample data from Illumina website. Start with simple things such as check sequence base quality, genome coverage, estimation of insert size, etc Following three bioinformatics tools will give you a good start: bedtools, samtools, bamtools. You can check official documentation for installation, if you are using a mac computer, here is a page explaining how to install them: http://genescripts.com/

Just try to use each command and figure out what they doing. You will learn a lot. Visit forums like this and read them. Try to take free courses online from Coursera, Mitedx, Udacity, etc. A great public tool is [Galaxy project](http://galaxyproject.org/). There are many tutorial video there. You will learn a lot. Have fun during your studies!
Question: Hi, currently, I am trying to calculate the coverage for my genes of interest (i.e. number of reads falling in or overlapping by the gene of interest) from DNA-seq data in a BAM file. Genes of interest are annotated in a BED file.

There are several [answers][1] on Biostar suggesting to use
`bedtools coverage` e.g., `bedtools coverage -abam sample.bam -b exons.bed -counts`. However, while checking the Samtools' manual, I stumbled upon `samtools bedcov` with the following description:

> read depth per BED region. Reports read depth per genomic region, as specified in the supplied BED file.

which can be used like this: `samtools bedcov gene.bed sample.bam`.

I have tested these two and I am getting different results! For instance, for `OR4F5`gene I get the following results (using the same BAM file for both):

chr1 69091 70008 "OR4F5" 61 #from bedtools coverage

and

chr1 69091 70008 "OR4F5" 4714 #from samtools bedcov

where `61` is very different than `4714`!

Another observation is that the result from `Samtools` is instantaneous for the above test while `Bedtools` takes a lot of time to produce the result. The speed itself has led me to think of using `samtools bedcov`. However, I was wondering whether there is catch that I am overlooking! I tried to find more information about `samtools bedcov` but I was not able to find anything more than that one-liner description.

I would be happy to hear about your feedback on this.

**[EDIT]:** This is how it looks like in IGV. It seems that `Bedtools` results is closer to reality.
![visualization in the IGV][2]

[1]: https://www.biostars.org/p/61748/
[2]: http://oi66.tinypic.com/2uztpue.jpg

Answer: The description that goes with `samtools bedcov` is wrong (I'll create an issue on github). It returns the sum of per-base coverage in each region. I'm not sure how that would be useful, but that's what it is.

**Edit:** [Here's the issue](https://github.com/samtools/samtools/issues/588) to get the description fixed.
Question: Hello guys,

It may seem like a basic question, but this is causing confusion.

What is the difference between fold change and Log fold change?

Regards,

Answer: If you think that fold-change is the 'expression level' in one set of samples (set A) divided by the 'expression level' in another set (set B), then log-fold-change is the log of that value (typically to base 2).

That is, if FC=, A/B, then log_FC = log(A / B) = log(A) - log(B)

and if log_FC=x, then FC=2^x