cannin/biostars_qa huggingface dataset:
Question: Today I had a lunchtime discussion with some colleagues about the existence of poly-allelic Single Nucleotide Polymorphisms (SNPs). Quoting Wikipedia to explain what I mean:

>For example, two sequenced DNA fragments from different individuals, AAGCCTA to AAGCTTA, contain a difference in a single nucleotide. In this case we say that there are two alleles : C and T. Almost all common SNPs have only two alleles.

My point was that there could be some but I don't know any SNPs that have eg. 3 alleles such that there is either e.g. A/G/T, or all 4 bases in a population. If I wish to query a database, say BioMart, UCSC, HapMap for this, how would I do this?

If wanted to detect poly-allelic SNPs from sequencing data (*de novo*) using a reference sequence, which method/program could be used?

Edit: So far this seems to be real great stuff, got at least two different methods of getting the query right!
But what about *de novo* (sorry, of course a reference sequence can be used, that's not really de novo) detection? Would for example [Maq](http://maq.sourceforge.net/) detect tri-allelic SNPs in short-reads?

Answer: You can query UCSC snp130 table. The fields of interest would be observed and class.
You can use their [public mysql server](http://genome.ucsc.edu/FAQ/FAQdownloads.html#download29) and run this sql query:

SELECT name,chrom,chromStart,chromEnd,observed,class,avHet
FROM snp130
WHERE length(observed)>3 and class="single" and avHet>0

I've added a condition to find only SNPs with average heterozygosity over zero. It finds 4949 SNPs in hg19, 490 are 4-allelic.
Question: Hello Forum,

What pipelines are available for exome sequencing analysis? Thanks for information.

Answer: Pipelines are built from existing tools. Pipelines are usually customized. You will benefit from reading up on GATK's best practices, and exploring Galaxy documentation on if they have an exome sequencing tutorial.
Question: I ran FASTQC on a human gut metagenome sample and found that I have a duplication rate of about 80% Does this seem too high? I checked out some environmental samples and saw approximately the same rates.

I've read papers that recommend de-duplicating reads before analysis because they're most likely PCR artefacts. But I've read papers that recommend keeping all reads since some high abundance species will be sequenced deeply and some reads may be seen more than once. Any thoughts on the matter?

Answer: High read depth on small templates or in your case low-diversity metagenomes does indeed lead to high duplicate rates (see your answer to Josh Herr above). We see this on amplicon datasets all the time.

For your analysis I would attempt analysis - if mapping reads - with and without the duplicates. Remember to add a Mapping Quality filter afterwards too. There will be information in both analysis methods. Remember to check which parts of the genome duplicates are coming from.

Typically, in our pipeline we remove duplicate and low-quality reads before analysis with several tools (mapping to human + bacterial genomes) and after mapping (eg. with Picard).
Question: From literature I know that a few TCGA samples have also information about treatment. It isn't clear to me how to obtain that information. After downloading the clinical data sets with TCGAbiolinks package I see one column "treatment_id". Each TCGA sample has a code in this column, looking like this: dbf26982-97b4-59f7-b165-d3889623bc21.

I can't figure out what that means, or how I can translate the code into drug name.
Any ideas?
Or is there an other way, how to obtain the treatment information outside the "clinical data information"?

Answer: The treatment/drug information can be downloaded using the TCGAbiolinks package.
Here is the command from the [vignette](https://bioconductor.org/packages/devel/bioc/vignettes/TCGAbiolinks/inst/doc/clinical.html#parse_xml_clinical_data):

query <- GDCquery(project = "TCGA-COAD",
data.category = "Clinical",
file.type = "xml",
barcode = c("TCGA-RU-A8FL","TCGA-AA-3972"))
GDCdownload(query)
clinical <- GDCprepare_clinic(query, clinical.info = "drug")

You can extract other kinds of clinical data such as "radiation" therapy using the same command. Just set `clinical.info` argument with `radiation`.
Question: Hi all,

I have these two .dat files (I only show the first 20 lines for both):

GO:0005509	PDCD6
GO:0004672	CDK1
GO:0005524	CDK1
GO:0005634	CDK1
GO:0005737	CDK1
GO:0006468	CDK1
GO:0005615	SERPINB6
GO:0006629	APOC2
GO:0006869	APOC2
GO:0008047	APOC2
GO:0042627	APOC2
GO:0043085	APOC2
GO:0001932	TADA2L
GO:0003677	TADA2L
GO:0005671	TADA2L
GO:0006357	TADA2L
GO:0007067	TADA2L
GO:0008270	TADA2L
GO:0016573	TADA2L

And

GO:0000001	mitochondrion inheritance
GO:0000002	mitochondrial genome maintenance
GO:0000003	reproduction
GO:0000005	ribosomal chaperone activity
GO:0000006	high affinity zinc uptake transmembrane transporter activity
GO:0000007	low-affinity zinc ion transmembrane transporter activity
GO:0000008	thioredoxin
GO:0000009	alpha-1,6-mannosyltransferase activity
GO:0000010	trans-hexaprenyltranstransferase activity
GO:0000011	vacuole inheritance
GO:0000012	single strand break repair
GO:0000014	single-stranded DNA specific endodeoxyribonuclease activity
GO:0000015	phosphopyruvate hydratase complex
GO:0000016	lactase activity
GO:0000017	alpha-glucoside transport
GO:0000018	regulation of DNA recombination
GO:0000019	regulation of mitotic recombination
GO:0000020	negative regulation of recombination within rDNA repeats
(...)

When I try to make a join for both files, I only get a few results (exactly 10). The complete code is:

ls *gene_association* | while read file;
do
echo;
echo @@@ File: $file;
echo;

# New file "assoc_specie.txt"
IFS='_' read -r -a array <<< "$file"
SPECIE=${array[2]}

#Filtering comments (!comment...)
cat $file | grep -v '!' > assoc_$SPECIE.txt;
gawk 'BEGIN{OFS="\t";FS="\t"}{print $5, $3}' assoc_$ESPECIE.txt > goTerms_$ESPECIE.dat;
join goTerms_$SPECIE.dat gene_ontology.dat > join.dat

echo
done;

I don't know what I am doing wrong, but it's obvious that join is not showing all the results.

Thanks in advance

PS: assoc_specie.txt file has this format (only showing first line):

UniProtKB	A0A024QZ42	PDCD6	GO:0005509	GO_REF:0000002	IEA	InterPro:IPR002048	F	HCG1985580, isoform CRA_c	A0A024QZ42_HUMAN|PDCD6|hCG_1985580	protein	taxon:9606	20160312	InterPro
(...)

Answer: Hi all,

I just saw the mistake. I forgot to sort the file:

(...)
gawk 'BEGIN{OFS="\t";FS="\t"}{print $5, $3}' assoc_$ESPECIE.txt | sort | uniq > goTerms_$ESPECIE.dat;
(...)
Question: I want to get all the GO terms under the biological process GO domain for any organism -right now, I need all the biological process associated to all gene products of E. coli (non-IEA ideally).

I'm new in the field of bioinformatics (I'm a software engineer by education) and right now I'm struggling with the great number of interconnected databases I'm finding.

At first, I tried with the R library biomaRt[1], but I discovered they don't support the Ensembl Bacteria database anymore.

Then I tried to do it by hand using AmiGO (actually GOlr[2]), but then I noticed that GO Central focused on human health-related annotations[3], so I supposed I would loose a lot of useful annotations for non-human species. After that I wanted to try QuickGO, but then I thought that with so many databases, repositories and tools around, it's a bit hard to know how exhaustive my request will be.

Should I stick with Ensemble databases and ditch biomaRt as a library?
Or using a combination of GOlr(AmiGO) and GAnnotation(QuickGO) enough?

**References**

- [1] Following this question https://www.biostars.org/p/1226/
- [2] http://wiki.geneontology.org/index.php/GOlr
- [3] http://geneontology.org/faq/what-are-differences-between-data-available-amigo-and-those-quickgo

Answer: Check out https://www.biostars.org/p/108669/

**Please use recent version of GeneSCF for successful results.**

You can extract all GO terms from gene ontology for your desired organism (Example below, Escherichia coli [ecocyc]) using simple command from GeneSCF,

(Complete GO for E.coli)

./prepare_database -db=GO_all -org=ecocyc
OR (Biological Process)

./prepare_database -db=GO_BP -org=ecocyc

OR (Molecular Functions)

./prepare_database -db=GO_MF -org=ecocyc

OR (Cellular Components)

./prepare_database -db=GO_CC -org=ecocyc

**Advantages**

- Real-time analysis, do not have to depend on enrichment tools to get updated.

- Easy for computational biologists to integrate this simple tool with their NGS pipeline.

- GeneSCF supports more organisms.

- Enrichment analysis for Multiple gene list in single run.

- Enrichment analysis for Multiple gene list using Multiple source database (GO,KEGG, REACTOME and NCG) in single run.

- **Download complete GO terms/Pathways/Functions with associated genes as simple table format in a plain text file (Check "Two step process" below in "GeneSCF USAGE" section)**.
Question: Hi all,

I have a fasta file of an assembled genome of a bacterial strain. I would like to identify the species of this strain (or the closest one possible) using the 16s rRNA sequence.

I tried downloading multiple 16s rRNA sequences from the NCBI database and to perform blastn analysis but I don't think I am using the most efficient method.

Do you have any recommendations on what tool to use to identify the species of a particular strain using 16s rRNA ?

Thank you in advance for your help!

Audrey

Answer: This is a good 16S reference file:
https://www.arb-silva.de/fileadmin/silva_databases/current/Exports/SILVA_138.1_SSURef_tax_silva.fasta.gz

And the map that goes with it:
https://www.arb-silva.de/fileadmin/silva_databases/current/Exports/taxonomy/taxmap_slv_ssu_ref_138.1.txt.gz

Although you may want to read the SILVA site first to determine if you want these files or perhaps e.g. a non-redundant higher quality version

https://www.arb-silva.de/

Blasting against the fasta file is fine as long as you place proper thresholds for hit reporting. Also, it will obviously be way faster if you're just blasting the 16S sequence and not the entire genome assembly (for ideas how to extract see Michael's reply)

That being said, I don't think 16S is the way to go for species identification, especially if you want strain/serotype level resolution
Question: Hi guys :D

I'm working with distance matrices produced by clustal omega for moderately large fasta files combining sequences of two different plant species in each .

When I was about to finish the script and code the final pipeline step ; which is retrieving the actual sequences corresponding to ID's given in the distance matrices using the biopython function `SeqIO.index()` ... I realized that the original fasta files have duplicate ID's for different sequences resulting from different positions of SSR's on the same sequence , in which I extracted the left and right flanking regions for each SSR .

```
Traceback (most recent call last):
File "C:\Users\Al-Hammad\Desktop\Test Sample\dictionary.py", line 9, in
dictionary=SeqIO.index("Left(Brachypodium_Brachypodium).fasta","fasta",IUPAC.unambiguous_dna)
File "C:\Python34\lib\site-packages\Bio\SeqIO\__init__.py", line 856, in index
key_function, repr, "SeqRecord")
File "C:\Python34\lib\site-packages\Bio\File.py", line 275, in __init__
raise ValueError("Duplicate key '%s'" % key)
ValueError: Duplicate key 'BRADI5G06067.1'

Tool completed with exit code 1
```

Here's a sample of one of my files :

```
>BRADI5G06067.1 cdna:novel chromosome:v1.0:5:7642747:7642899:-1 gene:BRADI5G06067 transcript:BRADI5G06067.1 description:"" Startpos_in_parent=24 Startpos_here=24 Length=26
ATGTATCTCCAACAACAACAACA
>BRADI5G06067.1 cdna:novel chromosome:v1.0:5:7642747:7642899:-1 gene:BRADI5G06067 transcript:BRADI5G06067.1 description:"" Startpos_in_parent=54 Startpos_here=54 Length=34
ATGTATCTCCAACAACAACAACAACGACGACGACGACGACGACGACGACAACG
>BRADI5G06067.1 cdna:novel chromosome:v1.0:5:7642747:7642899:-1 gene:BRADI5G06067 transcript:BRADI5G06067.1 description:"" Startpos_in_parent=102 Startpos_here=102 Length=26 ATGTATCTCCAACAACAACAACAACGACGACGACGACGACGACGACGACAACGACAACAACAACAACAACAACAACAACAACAACAAGAACGACGACGACG
```

**My question is**: What is the best, safest and most efficient way to rename the duplicate ID's for different sequences ?! and do I have to recompute the distance matrices again with the unique ID's after renaming or can I simply map the duplicates with their corresponding new unique values on the surface ?!

I'm really confused about that , and a little worried about the recomputing if considered since it's time consuming and takes nearly 4 days to produce the matrices .

I found this: http://stackoverflow.com/questions/7815553/iterate-through-fasta-entries-and-rename-duplicates/7836747#7836747 but it wasn't useful in my case, I'm working on a windows 7 64bit platform and python 3.4

Also I found this: https://www.biostars.org/p/10625/ but I believe it was the opposite of my case , I tried it though and ran on my files infinitely !! It wasn't that clear to me , for my bad luck :\

I desperately need this :( 😔

Any help would be appreciated, thanks in advance.

Answer: For the moment, I'd suggest replacing blank spaces in the header with double underscore characters (just a placeholder, you can use anything but underscores are among the most docile of characters)

You can install GnuWin32 tools for windows. Then, use this command:

sed -re "/^>/s/ /__/g" input.fa | sed -e "s/\_\_$//" >trInput.fa
Question: I have a dataframe which contains value of log2fold change but it contains inf and NA values i searched all over stack-exchange tried their solution but most of them seems not working few of them works but its not giving me the desired output ,need help it should't be that difficult i suppose .

my data frame

GENE_NAME `HSC_VS_CMP` `HSC_VS_GMP``HSC_VS_Monocytes`

ACTL6A -0.20084399 0.297430 -0.350876000
ACTR8 -0.2925280 -0.158551 1.10747
AICDA inf NA inf
ANP32B -0.6549 -0.615725 -0.35858

I get error like this " default method not implemented for type 'list'
"

So please suggest me how to remove all the inf and NA from the data frame

Answer: Alternative (borrowed from [this SO answer][1]):

# test data
d <- mtcars

# add NAs and Inf.
d[1,1] <- NA
d[2,2] <- NA
d[2,1] <- Inf
d[1,2] <- Inf

head(d)
mpg cyl disp hp drat wt qsec vs am gear carb
Mazda RX4 NA Inf 160 110 3.90 2.620 16.46 0 1 4 4
Mazda RX4 Wag Inf NA 160 110 3.90 2.875 17.02 0 1 4 4
Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1

# the magic:
d <- do.call(data.frame, lapply(d, function(x) {
replace(x, is.infinite(x) | is.na(x), 0)
})
)

# note that you lose the row names.
head(d)
mpg cyl disp hp drat wt qsec vs am gear carb
1 0.0 0 160 110 3.90 2.620 16.46 0 1 4 4
2 0.0 0 160 110 3.90 2.875 17.02 0 1 4 4
3 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
4 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
5 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
6 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1

[1]: http://stackoverflow.com/a/12188551/3742244
Question: We sequenced two years ago two RNA samples using RNA-seq with 200 bp Ion Torrent kits. Now we need biological replicated for these two samples. The sequencing service offers now the possibility to get 200 or 400 bp.

​Do you think I should upgrade to 400 bp or stay with 200 bp in order to avoid technical variability?

Answer: We have a similar experience. We integrated 50 bp reads produced by Illumina HiSeq2000 and 75 bp reads produced by Illumina NextSeq500 in a chicken project. Fortunately, it didn’t make a problem in our case. However, a paper showed that read length could be a confounding factor as below for your reference.

Danielsson F, James T, Gomez-Cabrero D, Huss M (2015) Assessing the consistency of public human tissue RNA-seq data sets. Brief Bioinform.
Question: Dear all, I have easy script in bash using awk language. I filtrate bam file by GC content (see at code), but I need to get on output bam file too with header. Could you help me with this? thank you My code:

#!/bin/bash

for i in 54321*.bam
do
samtools view -h /home/filip/Desktop/Analyza\ NIFTY\ and\ CNVs/NIFTY\ pooling009/$i | awk '{ n=length($10); print gsub(/[GCCgcs]/,"",$10)/n,"\t",$0 }' | awk '($1 <= 0.6 && $1 >= 0.3){print $0}' | awk '!($1="")' | samtools -bS - > z$i;
done

Answer: As long as you don't have any reads beginning with "@" (and you shouldn't), just add an if(substr($1,1,1) == "@") { print $0} else {...your current code...}. You just want to avoid all of the piping and have a single awk command (unless you want to modify each of the pipes, which would quickly become annoying).

BTW, this might be easier with pysam. Awk is great, but it can get a bit unwieldy with longer commands.
Question: Hi all,

I'm trying to run CNVkit (from source, in linux) after following the [installation instructions][1] on the [github page][2] but running into import errors which seem to be rare as I couldn't find any references online:

Traceback (most recent call last):
File "cnvkit.py", line 8, in
from cnvlib import commands
File "/home/alon/Desktop/BioinformaticsPrograms/cnvkit/cnvlib/__init__.py", line 4, in
from .cmdutil import read_cna as read
File "/home/alon/Desktop/BioinformaticsPrograms/cnvkit/cnvlib/cmdutil.py", line 9, in
from .cnary import CopyNumArray as CNA
File "/home/alon/Desktop/BioinformaticsPrograms/cnvkit/cnvlib/cnary.py", line 14, in
from .segmetrics import segment_mean
ImportError: No module named segmetrics

Will greatly appreciate any help!

Alon

[1]: https://github.com/etal/cnvkit#from-source
[2]: https://github.com/etal/cnvkit

Answer: This is fixed now. Sorry for the trouble, I forgot to check in a new file while refactoring.
Question: Hi all,

I want to do survival analysis of clinical data but I am not sure about censoring TCGA clinical data? I have ''daysto_lastfollowup", "daysto_death", "vital status" and "overall_survival_months". What I did is

**Method 1:**

daysToEvent <- rep(NA, nrow(tcgaClinical))
daysToEvent[vitalStatus == "Alive"] <- daysToLastFollowup[vitalStatus == "Alive"]
daysToEvent[vitalStatus == "Dead"] <- daysToDeath[vitalStatus == "Dead"]
eventStatus <- rep(NA, nrow(tcgaClinical))
eventStatus[vitalStatus == "Alive"] <- 1
eventStatus[vitalStatus == "Dead"] <- 0

tcgaOS <- Surv(daysToEvent/30, eventStatus == 0)
rownames(tcgaOS) <- rownames(tcgaClinical)

Or

**Method 2:**

eventStatus <- rep(NA, nrow(tcgaClinical))
eventStatus[vitalStatus == "Alive"] <- 1
eventStatus[vitalStatus == "Dead"] <- 0

tcgaOSM <- Surv(tcgaClinical$OS_MONTHS, eventStatus == 0)
rownames(tcgaOSM) <- rownames(tcgaClinical)

Which method makes more sense ? Thank you

Answer: Hello, I recommend you download clinical data from [cBioPortal](http://www.cbioportal.org/) instead of TCGA/GDC. If you download clinical data from cBioPortal you will see fields `Overall Survival (Months)` and `Overall Survival Status` thats what you need for OS(Overall survival) analysis.

As OS the event is dead, so `alive` is censored data, you should give value `1` to `dead` and `0` to `alive`. However I recommend use `1` and `2`, give value `2` to `dead` and `1` to `alive`.

Here, some tutorial may be useful for you:
https://www.emilyzabor.com/tutorials/survival_analysis_in_r_tutorial.html#introduction
https://www.datacamp.com/community/tutorials/survival-analysis-R
http://www.sthda.com/english/wiki/survival-analysis-basics
Question: Hi Everyone,

I am trying to run a mysql query and save output into a text file. The command works fine for 10 SNPs and the output file has entries for all 10 SNPs, however when I try to run the same for large list (100 K) of SNPs the output file only contains the entry for the last SNP only.

My command -

for i in $(cat /path/to/directory/listofSNPs.txt); do echo "select chrom,chromStart+1,name,alleles from snp150 where (name= \"$i\")" | mysql -N --user=genome --host=genome-mysql.cse.ucsc.edu -A -P 3306 -D hg38 done > '/path/to/directory/output.txt'

I would appreciate any advice on how can I do this?

Thanks,
Kanika

Answer: faster:

awk '{printf("select chrom,chromStart+1,name,alleles from snp150 where name= \"%s\";\n",$1);}' listofSNPs.txt | \
mysql -N --user=genome --host=genome-mysql.cse.ucsc.edu -A -P 3306 -D hg38 > /path/to/directory/output.txt
Question: Hi everyone,

I've been looking at gene expression patterns in one dataset (dataset A) for which I have Ensembl gene identifiers obtained using the annotation file corresponding to Ensembl version 86 (GRCh38.p16). I now want to look at the patterns in another dataset (dataset B) for which the genes are identified by gene symbols (annotated using Gencode release 31 GRCh38.p12). I would like to match the identifiers from both datasets and considered several strategies, however, I'm struggling to find which solution would be the best:

1/ Convert Ensembl gene identifiers from dataset A to gene symbols using biomaRt and the version corresponding to Gencode release 31

Possible issues: I'm unsure which version of Ensembl would be the best to do that, will it be the first Ensembl version with patch 12 (that would be version 92 from April 2018)? Or is this more complicated than that?

2/ Convert gene symbols from dataset B to Ensembl gene identifiers using biomaRt and the Ensembl version corresponding to Gencode release 31 (used for the annotation of dataset B)

Possible issues: As previously, I'm not sure which Ensembl version to use and also I'm not entirely sure Ensembl gene identifiers will refer to the exact same things between the version?

3/ Convert gene symbols from dataset B to Ensembl gene identifiers using biomaRt and the version corresponding to Ensembl version 86 (used for the annotation of dataset A)

Possible issues: Pretty much the same ones as 2/

4/ Re-annotate one of the two sets so that they are both annotated with the same file

Possible issues: I don't have access to the raw data...

I'm also thinking that I could be misunderstanding the 'stability' part of Ensembl gene identifiers or possibly overcomplicating the whole thing...

In any case, I'd appreciate some suggestions/explanations to help clear the fog away from my mind...

Thank you!

Answer: The stable gene IDs (ENSGs) should refer to the same entity between releases, and are far more stable than the gene names, which change more frequently. I would recommend using the current BioMart to convert your list of gene names from the current release (GENCODE 31 = Ensembl 97 = latest data) to ENSGs, then you can reasonably compare these to your data from Ensembl 86.
Question: I have multiple vcf of CASES and CONTROLS variations annotated by VEP, SNPEff, SnpSift.

first pair vcf -> only variations| CASES and CONTROLS
second pair vcf -> variations + SnpEff | CASES and CONTROLS
third pair vcf-> variations + SnpEff+VEP+SnpSIFT| CASES and CONTROLS | at so on

Because it have non-overlapping samples, i tried to merge corresponding vcf of CASE and CONTROLS into one file with bcftools merge.

bgzip -c file.vcf > file.vcf.gz
tabix -p vcf file.vcf.gz
bcftools merge -o merged.vcf.gz -Oz file.vcf.gz file_2.vcf.gz

Unfortunatelly in files with SnpSift annotations i got this error:

Error at chr1:20724637: wrong number of fields in dbNSFP_ada_score?

Please anybody knows how to solved this error or what does it mean?

Answer: your vcf is badly annotated.

there is a VCF with 6 ALT alleles :

chr1 20724637 rs1280663306;rs1445195248 TGGGAGGGAGGGAGAGAGGT TGGGAGGGAGAGAGGT,TGGGAGGGAGGGAGGGAGGT,TGGGAGGGAGGGAGAGAGGG,TGGGGGGGAGGGAGAGAGGT,TGGGTGGGAGGGAGAGAGGT,TGGGAGAGAGGGAGAGAGGT

and the annotation for dbNSFP_ada_score is `INFO=