cannin/biostars_qa huggingface dataset:
Question: When I analyze somatic mutations, I usually judge them according to many annotation filters such as conserved genomic regions, segmental_duplications and integrated databases (cosmic). Especially, I didn't believe any filter-based annotation in annovar which are SIFT or mutation_assessor because some cases have lower score but they are likely to be relevant to the cancer-disease or have dominant effects for disease.

From this points, I am wondering what useful annotation filter you use for mutation in exome seq?

Does any definite and reliable annotation filter exist?

It is always hard for bioinformaticians to differentiate artifacts or false_positives from real variants.

Answer: I do not believe that there is a definite and reliable filter but I do like [KGGSeq][1] very much because:

- it combines several prediction programs into 1 variable
- it uses many different databases
- you can add your own databases

So strategy depends on the data you have, we usually have trio's and look for denovo variants, so:

1. only keep denovo: WT in parents and HET in affected
2. remove all that are present in ExAc database
3. remove synonymous, non-frameshifts and non-coding
4. prioritize based on the prediction programs score (so no filtering here)

[1]: http://statgenpro.psychiatry.hku.hk/limx/kggseq/
Question: Alright, that didn't take long. So I made my script to make my bed file which worked.

```
Track name=G_anatisUMN179 description=G_anatisUMN179 genes
G_anatisUMN179_dna 0 2687335 UMN179_RS12665 +
G_anatisUMN179_dna 210 1068 UMN179_RS00010 +
G_anatisUMN179_dna 1167 1980 UMN179_RS00015 -
G_anatisUMN179_dna 1972 2215 UMN179_RS00020 -
G_anatisUMN179_dna 2285 2618 UMN179_RS00025 +
G_anatisUMN179_dna 2671 4069 UMN179_RS00030 +
G_anatisUMN179_dna 4188 5592 UMN179_RS00035 -
G_anatisUMN179_dna 5861 7073 UMN179_RS00040 +
G_anatisUMN179_dna 7241 8294 UMN179_RS00045 -
G_anatisUMN179_dna 8309 9716 UMN179_RS00050 -
G_anatisUMN179_dna 9768 10383 UMN179_RS00055 +
G_anatisUMN179_dna 10460 10676 UMN179_RS00060 -
G_anatisUMN179_dna 10647 11184 UMN179_RS00065 -
G_anatisUMN179_dna 11294 11600 UMN179_RS00070 -
```

4 mandatory columns ripped right from the genbank file, however when I try to load the file into IGV I get the following error.

Error loading /home/sbrimer/genomes/G_anatisUMN179.bed: Error parsing line at byte position: htsjdk.tribble.readers.LineIteratorImpl@42f867b5, for input source: /home/sbrimer/genomes/G_anatisUMN179.bed

I've tried googling the error and looking into the igv-help group on google but I'm not really finding an answer, did I format my file incorrectly or what step did I miss?

Answer: So besides having the header off, it turns out that the bed file I generated from the genbank file had a different id then the fasta version I download as originally to do the alignements. So, in the future I have learned to use the script from the start and I have have become very, very familiar with the BED formate.
Question: Hi all,

Mammalian whole Genome/Exome Sequencing data analysis in particular de-novo assembly based projects usually require >512 Gbs or even >1TB RAM. I get such dataset twice or thrice a year. Instead of buying a computational machine with this much of RAM will cost >10,000 Euros. Are there any institutions/companies based in Europe which provide paid access to such massive RAM computational clusters by taking care of data privacy and security?

Thanks for your feedback!
Rahul

Answer: Thanks guys for the helpful comments. Here are the services I could short-list along with the pricing information:

Google Compute:
https://cloud.google.com/compute/pricing?authuser=1

Microsoft Azure:
https://azure.microsoft.com/en-us/pricing/details/virtual-machines/linux/

Amazon AWS:
https://aws.amazon.com/ec2/pricing/on-demand/
, https://aws.amazon.com/ec2/pricing/reserved-instances/

For me Google Compute fits well for computation intensive jobs.

Hope this will be helpful.

Best
Rahul
Question: Hi Community,

I have a multi-column file, but here I am mentioning just two columns:

Contig0.100010-snap.2	93.4065934066
Contig0.100014-snap.1	67.3611111111
Contig0.100025-snap.1	72.9411764706
Contig0.100025-snap.1	74.1176470588
Contig0.100025-snap.1	82.3529411765
Contig0.100051-snap.1 95

Column one has three duplicated values, that is, Contig0.100025-snap.1. So, I want that if column one has duplicated values, remove those line in which the corresponding values in column 2 are smaller. My desired output would be:

Contig0.100010-snap.2	93.4065934066
Contig0.100014-snap.1	67.3611111111
Contig0.100025-snap.1	82.35229411765
Contig0.100051-snap.1 95

Any unix sed awk etc or help from excel???

Thanks,

Waqas.

Answer: Simple bash sort solution:

sort -k1,1 -k2,2nr test_file | sort -u -k1,1
Question: I have reads from a bacterial sequence and I want to identify variants in them. I have performed an alignment using `bwa mem` and the reads aligned well. Is there a good pipeline for identifying SNPs, small indels etc in microbes? I could use mpileup but I don't find many people recommending that.

Answer: Well, dunno about other people, but I like `samtools mpileup`. It has always worked well for me in the microbial context:

samtools mpileup -uf [genome] [bam] | bcftools view -vcg - > vcfs > [filename].vcf
Question: Dear all,

Do you have any idea how to easy extract contigs from fasta file wich contains specific sequence?

For example:

My sequence:

ACCGTACCC

My FASTA:

```
>c1042
ACCGTACCC
>c1043
GCTACAGTTGAAAGGGGACCGTACCC
>c1044
ATGAATAAAATAATTTTGTATCATAAATCGAGCTGTTAATTATT
>c1044
TTCATATTTGTAGCTAAGCAGAGGCGAAGCGTTCTTGTATCG
```

My output:

```
>c1042
ACCGTACCC
>c1043
GCTACAGTTGAAAGGGGACCGTACCC
```

Thank you so much for any ideas and help.

Answer: cat fasta.fa | grep -B1 "ACCGTACCC" > out.fa
Question: So firstly, I'm completely aware that this type of question has been asked multiple times (I know this since I've been scrolling over these type of questions for the past 2 days), but I'm actually more interested in knowing the reasons as to why some people prefer

I've performed differential expression analysis using DESeq2 and I want to see which Gene ontology terms, KEGG pathway terms etc are enriched in my data set. I've initially tried using clusterProfileR, but I keep getting 3 enriched terms for all my differentially expressed genes using ``enrichGO()``. I also know that some input in clusterProfileR requires you to put logFC values, so I wasn't sure if that was for ALL the genes analysed, or just the differentially expressed genes.

I've also used goseq but my main issue with that is the GO terms are too broad.

I also only have about 300 DEGs, so I'm not really sure if this sort of analysis is best performed when you have a myriad of DEGs, or can be done with a small number.

Anyway, looking forward to hearing people's responses :)

Answer: I'm a big user of `GOseq`, which allows to control for the fact that longer and more highly expressed genes are more likely to be found to be differential, if long or highly expressed genes are not evenly distributed between pathways/categories, then this can bias your enriched pathways.

I also like the `GSEA`-like algorithms, because you do not have to set an artificial limit on what you consider significant. The version of this where you rank on significance suffers from the same gene-length bias that traditional GO tests suffer, but this may be lessened by ranking on some suitably strunken logFC metric. `cameraPR` from limma is a good example of this.

Finally we've used `SPIA` before, which is a pathway enrichment tool that takes the topology of the network into account. Its a great idea, let down by the quality of the pathway annotations it runs on.
Question: I am new to bioinformatics and just for the sake of this question, suppose that I am looking for the full genome sequence of [*Escherichia coli*][1] and [*Saccharomyces cerevisiae*][2] in FASTA format. I thought that these would be fairly easy to find online.

For Escherichia coli I was only able to find some links to the ncbi website but I could not find anywhere in the website a link to download the genome in FASTA format. For Saccharomyces cerevisiae, I found [this FTP site][3] linked by ensembl with the sequence for each chromosome. I suppose that I could concatenate all the files in order to obtain the full genome in one file but isn't there already one file with the entire sequence? I saw the top level files but those only appear for masked and soft masked DNA, I think it would be doable to work with those if modern sequencers masked the DNA automatically, do they?

So, in summary, how do experienced people find full reference genomes for those two organisms and possibly other common species? I have only been using Google and the only page that appears to provide FASTA files for some species is Ensembl.

[1]: http://en.wikipedia.org/wiki/Escherichia_coli
[2]: http://en.wikipedia.org/wiki/Saccharomyces_cerevisiae
[3]: ftp://ftp.ensembl.org/pub/current_fasta/saccharomyces_cerevisiae/dna/

Answer: The short answer is:

this is a rare case where simply "Googling" will not help you very much; Google works using "popularity" and most people are not interested in "genome sequences fasta format" :)

"experienced people" know that there are a few key sites where genomic data can be found, as covered in previous answers - NCBI, Ensembl, UCSC, JGI

and then it is simply a case of exploring those sites and bookmarking the most useful areas for future reference - voilà, you are now experienced
Question: I want to produce a list of all the somatic mutations that are known using the 1000 genome project and using ensembl. How do I start this project?

Answer: While you could use BioMart or the API to retrieve these data, you might be better off with VCF file dumps:

ftp://ftp.ensembl.org/pub/release-79/variation/vcf/homo_sapiens/

`Homo_sapiens_somatic.vcf.gz` (all somatic mutations from COSMIC and dbSNP) and Homo_sapiens.vcf.gz (all non-somatic variants from various sources, see http://www.ensembl.org/info/genome/variation/sources_documentation.html) from this directory should be a good place to start.

If you specifically only want those from the 1000 genomes project, you could grep for the 1000 genomes evidence flag:

curl ftp://ftp.ensembl.org/pub/release-79/variation/vcf/homo_sapiens/Homo_sapiens.vcf.gz | gzip -dc | grep E_1000G | head

or download the population-specific files from the above directory, or the VCF files from the original project site (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/).
Question: Hi all,

Our lab has sequenced 20 different strains of the same bacterial species. I now want to analyse these data, i.e., link genomic data to virulence (some of the strains are virulent, others aren't) and other phenotype differences.

I started with de novo assemblies and mapping to an annotated reference genome. However, it is unclear to me how to proceed now in a good way. Most of the information for comparative genomics I find on the web applies only to the human genome or one-to-one comparisons of bacterial genomes.

Is there a standard workflow for this kind of analysis? Does anybody know of any good tutorials/references/articles? How would you proceed?

Answer: Part of my thesis was on comparative genomics to get closer to uncovering the genomic basis of virulence in *Neisseria meningitidis*. I think that two good papers to look at are [Schoen et al 2008][1] and my own paper [Katz et al 2011][2].

Schoen compared genomes on the gene level; my group looked for SNP-level differences. For the genome project part of it, ie assembly-to-annotation, we used [CG-Pipeline][3]. However, a more interactive service like [RAST][4] or [BASys][5] might be easier (although less comprehensive in annotations). NCBI's PGAAP annotations however would make it easiest to submit the genomes to NCBI at the end of the study.

Anyway, I think that the methods in Schoen et al and Katz et al would be very relevant to your questions.

[1]: http://www.pnas.org/content/105/9/3473.full.pdf
[2]: http://jb.asm.org/content/193/14/3633.short
[3]: http://sourceforge.net/projects/cg-pipeline/
[4]: http://rast.nmpdr.org/
[5]: http://basys.ca/basys/cgi/submit.pl
Question: Is there a way to download individual mitochondrial sequences from NCBI without having to scrape from their webpages or to download thousands of other sequences alongside. Is there some programmatic way to access these fasta data, for example:

https://www.ncbi.nlm.nih.gov/nuccore/NC_012920.1?report=fasta

Answer: Following this [post][1]:

curl "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=NC_012920.1&rettype=fasta&retmode=text" >mtDNA.fa

[1]: http://angus.readthedocs.io/en/2016/howe-ncbi.html
Question: Excuse me:

The SureSelect protocol said that we should add 'A' bases to the 3' end of the DNA fragments. Could someone tell me the purpose of this procedure?

Many thanks in advance!

Answer: As far as I know, this is not specific to SureSelect, it's all Illumina ligation-based protocols.

When you ligate the sequencing adapters to the target insert, it helps get adapter to insert ligation, because the adapters have a T overhang on their 3' ends and can pair. The inserts or adapters ligate to themselves much less efficiently because they will have A-A or T-T mismatches.
Question: Dear all,

I have indexed bam file, and I want to print all alignments on chromosome "2" using pysam.

My code is

import pysam
bam = pysam.AlignmentFile("Aligned.sortedByCoord_rep1.out.bam", "rb")
for line in bam.fetch("2"):
print line

What I get is:

D00733:162:CADM2ANXX:2:1311:7466:63232	163	18	3328	255	1S47M2S	18	3419	47	TTGCTTAGTGTCCGAAATACCATCCTCAAGGCTAAGAACTAAATCGATTA	array('B', [33, 33, 33, 33, 33, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38])	[('NH', 1), ('HI', 1), ('AS', 70), ('nM', 10)]
D00733:162:CADM2ANXX:2:1215:11810:89291	163	18	3356	255	44M6S	18	3441	44	GGCTAAGAACTAAATCGATTATTCTGGCTCGTAACGCATATAATATGGCA	array('B', [33, 33, 33, 32, 33, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 35, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 37, 38, 38, 38, 38, 38, 38, 33])	[('NH', 1), ('HI', 1), ('AS', 70), ('nM', 10)]

So there are 2 questions:
1. Why pysam does not preserve format of bam file?
2. Why does it print out chromosome 18 instead of 2?

I would appreciate if somebody points out on what am I doing wrong.
Thanks

Answer: In order to print an AlignedSegment out as a formatted SAM entry you need to use the `tostring` method and supply the AlignmentFile that the AlignedSegment came from. This contains the mapping of the numeric contig ids to the string contig names and also the encoding for the quality scores.

Thus your code should be:

import pysam
bam = pysam.AlignmentFile("Aligned.sortedByCoord_rep1.out.bam", "rb")
for line in bam.fetch("2"):
print line.tostring(bam)

following on from your comment, to filter on an NH tag you would do the following:

import pysam
bam = pysam.AlignmentFile("Aligned.sortedByCoord_rep1.out.bam", "rb")
for line in bam.fetch("2"):
if not (line.has_tag("NH") and line.get_tag("NH") == 1):
print line.tostring(bam)

you only need the `has_tag` if only some of your lines have the NH tag.
Question: Hi

I have a question about how to make cluster for relative expression (on read counts) using k-means clustering, I couldn't find any good manual for solving this problem.

I need to make a plot like below: to have different cluster on genes count:

Does anybody know how to make it?

![enter image description here][1]

[1]: https://i.postimg.cc/qRmPVccS/Screen-Shot-1397-07-20-at-20-13-12.png

Answer: I have solved the problem, the code should be something like below:

You need to run DESeq2 before processing rowMeans:

library(ggplot2)
library(DESeq2)
#suppose you have one column per each time/tissue condition

cds <- DESeq(cds)
vst <- getVarianceStabilizedData(cds)
df <- vst
df <- df - rowMeans(df)
df <- as.data.frame(df)
colnames(df)
res_km <- kmeans(df, 5, nstart = 10)
data_plot <- data.table(melt(data.table(class = as.factor(res_km$cluster), df)))
data_plot[, Time := rep(1:ncol(df), each = nrow(df))]
data_plot[, ID := rep(1:nrow(df), ncol(df))]
head(data_plot)
# prepare centroids
centers <- data.table(melt(res_km$centers))
setnames(centers, c("Var1", "Var2"), c("class", "Time"))
centers[, ID := class]
centers[, gr := as.numeric(as.factor(Time))]
head(centers)
head(data_plot)
# plot the results
ggplot(data_plot, aes(variable, value, group = ID)) +
facet_wrap(~class, ncol = 2, scales = "free_y") +
geom_line(color = "grey10", alpha = 0.65) +

geom_line(data = centers, aes(gr, value),
color = "firebrick1", alpha = 0.80, size = 1.2) +
labs(x = "Time", y = "Load (normalised)") +
theme_bw()
Question: I have a list of gene symbols.

APOB, TTC39B, ATF3, RGS1, LIPG,

I am trying to get the genomic coordinates (as in bp) with +/-5KB position of these genes via UCSC TableBrowser / MySQL server using NCBI 36/hg18 build. I have tried to get this information via TableBrowser, but I can't find the chromStart and chromEnd field in the given table. Am I missing something ?

Also, please share your favorite tutorial / docs that explain the schema/tables in UCSC MySQL server.

Thanks in advance.

Answer: We used to love the annotation database overview here: http://genome.ucsc.edu/goldenPath/gbdDescriptionsOld.html And we used to show it to people in trainings all the time. They loved it. It was much easier to scan through and get a sense of what was available.

But it went away a few years back the upgrade that gave use the "describe table schema" button. We've been asking for a page like it again, but haven't come across it yet.
Question: I have been trying to get GroopM to run on our local servers for quite some time now.
GroopM has been recommended n some fairly recent reviews, but i simply can't figure out how they got it ro run...

Installation with pip consistently fails with error messages mentioning someting about the "oldest supported numpy version" not being available

installation with conda seems to work at first, BUT is non-functiontal. Everytime you try to run that GroopM installation on actual data fails with an error about decarpicated code being used.

compiling from source does not work, despite all dependencies being installed, as i get tons of "syntax errors" while running setup.py...

Can anybody give me any tips please?

Answer: so I finally managed to get it installed by adapting v.berriosfarias procedure to:

a.) first create env with most dependencies:

conda create -n groomp2 python=2.7 matplotlib scipy bamm cython pillow pytables tempdir diamond hmmer pplacer mfqe orfm

b.) obtain pip3 and use it to install singlem globally (while currently in the base environment OR any python3 environment)

sudo apt update
sudo apt install python3-pip
pip3 install singlem

(singlem will automatically be in PATH for all environments, since it was installed globally)

c.) switch to groopm2 conda environment, download groopm2 and install

conda activate groopm2
git clone https://github.com/timbalam/GroopM.git
cd GroopM
python setup.py install

Thanks so much for your help!
Question: Hi everyone,

I am trying to find allele specific expression from RNA-seq data. I have the results from `GATK ASEReadCounter` and the output looks like this:

contig	position	variantID	refAllele	altAllele	refCount	altCount	totalCount	lowMAPQDepth	lowBaseQDepth	rawDepth	otherBases	improperPairs
chr2	38021	rs113895774	C	T	0	1	1	0	0	1	0	0
chr2	217334	rs6709534	G	C	0	3	3	0	0	3	0	0
chr2	218386	rs9213	G	A	84	101	185	0	0	185	0	0
chr2	220889	rs3828165	G	A	9	7	16	0	0	16	0	0
chr2	221560	rs60484953	G	A	11	11	22	0	0	22	0	0
chr2	221981	rs3791224	C	T	3	4	7	0	0	7	0	0
chr2	222336	rs3791223	T	C	3	8	11	0	0	11	0	0
chr2	224086	rs1474053	T	A	3	2	5	0	0	5	0	0
chr2	224919	rs2290911	A	G	55	16	71	0	0	71	0	0

I had an impression that this tool implements some test(s) to find statistically significant sites. But if I am not wrong, it only calculates the counts of ref and alt alleles based on RNA-seq - more or less like `bam-readcount` or `samtools mpileup` followed by [counting bases][1].

In the manual of ASEReadCounter, they say that you can use the output format as input to [mamba][2] but I am not sure looking at the input format for mamba which requires an additional field: EXON_INFO - variant annotation label. Alternatively, I am thinking of using the refCount, altCount and totalCount, perform a `chi square test` to determine if there is allelic imbalance or not.

I would like to get suggestions on how to analyze this output or what downstream methods/statistical tests to use. Any help would be much appreciated.

Thanks!

[1]: https://github.com/riverlee/pileup2base/blob/master/pileup2baseindel.pl
[2]: http://www.well.ox.ac.uk/~rivas/mamba/

Answer: ASEReadCounter is doing what it says its doing..just counting the ref and alt reads per allele..

A chi square test or a two sided binomial test is fine for determining locus based ASE..for gene based there are other methods..you can get more information from this paper- http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=7297750
Question: Dear all,

I need help with parallel command. I have one script in shell. And I would like to run it 12time in same time. But for each script I need different name of output. output is .tsv and the name is same like name of input, could you help how to do that?

Thanks a lot

Answer: an example with echo:

seq 1 12 | parallel echo "Hello" '>' 'result.{}'
Question: Hi all,

I have made a script where I can extract protein sequences in a fasta file from protein IDs.

However, I have this issue that I can't solve. Here is an example of my text file containing protein IDs:

Lactococcus_104
Lactococcus_105

Here is an example of my fasta file containing protein sequences:

[...]
>Lactococcus_104
MXXXXX
>Lactococcus_105
MXXXXX
[...]
>Lactococcus_1050
MXXXXX
>Lactococcus_1051
MXXXXX
>Lactococcus_1052
MXXXXX
[...]

I extract the protein sequences with this command that works perfectly:

IFS=$'\n'; for i in $(cat IDs.txt);do line=$(grep -nr "$i" file.fasta); if [[ ! -z $line ]];then for j in $line;do lineNumber=$(echo $j | cut -d':' -f1); sed -n "$lineNumber p" file.fasta; awk -v nb=$lineNumber 'NR > nb {if ($0 ~ ">") exit; else print $0 }' file.fasta; done;fi;done > output.fasta

But the problem is that the script "extracts" all protein IDs that contains, in my example, "104" or "105" so it also "extracts" protein sequences that I don't want like Lactococcus_1051, Lactococcus_1052, etc. instead of only Lactococcus_105.

**Is there a way that I can modify my command line to specify that I strictly want the extraction of the protein IDs in my text file?**

If you have any suggestions, please let me know!
Thank you in advance for your help.

Answer: Have you tried `grep` with `-w` like below:

IFS=$'\n'; for i in $(cat IDs.txt);do line=$(grep -nw "$i" file.fasta); if [[ ! -z $line ]];then for j in $line;do lineNumber=$(echo $j | cut -d':' -f1); sed -n "$lineNumber p" file.fasta; awk -v nb=$lineNumber 'NR > nb {if ($0 ~ ">") exit; else print $0 }' file.fasta; done;fi;done > output.fasta
Question: I want to get all the mRNA or protein sequences in a specific pathway, e.g. glycolysis or oxidative phosphorylation, from fruitfly. Is there a convenient way to attain this goal?

I first thought of KEGG database; however, KEGG pathway does not possess pathways specific for fruitfly. (It seems that kegg pathway only has specific pathways for vertebrates.)

Then I turned to FlyBase and found QueryBuilder help a bit.But I am not sure whether I got all genes and no extra in oxidative phosphorylation for fruitfly by setting the controlled vocabularies to be "oxidative phosphorylation".

I also tried BioMart and failed, because I am not familiar with it. I hope that someone shows me how to do it if BioMart is the right choice.

If there are other tools that can do this job, I would like to try.

Many thanks!

Answer: http://www.genome.jp/kegg-bin/show_pathway?org_name=dme&mapno=00190&mapscale=&show_description=hide