cannin/biostars_qa huggingface dataset:
Question: Hi guys,

I am not sure, in terms of sequencing what is referred to as color space and nucleotide space. What is a difference between a normal reference genome and a color space genome. Does it has to do something with base calling, that the original data is retained without the conversions to textual nucleotide data!!!! What we use in normal chip-seq and rna-seq analysis.

There's a kind of similar question here but the answer is missing the link.

Thanks

Answer: Short answer: color space refers to the native format of ABI SoLID technology. Color space is translated to nucleotide, or base space (same thing) so that it can be understood.

That technology is not growing in market share, so in the next few years it will become less common. ABI is putting most of their effort behind the Ion Torrent now.
Question: Hi,

Is it possibly to only output dosages from a generated vcf dosage file from minimac?

Input format is

##fileformat=VCFv4.1
##filedate=2017.7.5
##source=Minimac3
##contig=
##FORMAT=
##FORMAT=
##INFO=
##INFO=
##INFO=
##INFO=
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	1242658141	1364665948	1242658615
29	11	Chr29:11	A	G	.	PASS	.	GT:DS	0|0:0.193	0|0:0.193	0|0:0.193

Is there an option in vcftools to output the dosages for each genotype seperatly, without the genotype call?

Answer: bcftools query -f "%CHROM\t%POS[\t%DS]\n" minimac3.dose.vcf.gz
Question: How to download all genes inside all gene sets in C6 (oncogenic signature gene sets, 189 gene sets) in Molecular signature database (Msigdb) I want to download the all genes inside 189 gene sets in R in a data frame?

Thanks in advance

Answer: ```r
library(msigdbr)

genesets = msigdbr(species = "Homo sapiens", category = "C6",subcategory = NULL)
View(genesets)

unique_genes <- genesets%>% distinct(gene_symbol)
View(unique_genes)
```
Question: Hello everyone,
Is it possible to get the telomeric locations in a file for every chromosome?

Answer: [UCSC Table Browser][1]

Select
group: All Tables
table: gap

Then filter "create" button and type match "telomere". Then submit and get output.

[1]: https://genome.ucsc.edu/cgi-bin/hgTables
Question: Hello everyone,

I have a question regarding my single-cell RNA-seq data. I have the following pair-end data in `fastq.gz` format.

Read1 (contains 6bp UMI, followed by 6bp cell barcode info and the rest is a polyT stretch):

@J00182:79:HV2WWBBXX:6:1101:11160:38873 1:N:0:ACAGTG
GAGAAGACAGTGTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCTT
+
AAAFFJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJ--A-----7AJJFAF-AJJJJJJJJJJJ

Answer: Do yourself a favour and simply use `umi_tools extract`. The headers won't be formatted how you want, but will instead be formatted in a way more compatible with other tools.
Question: Hello,

I want to ask a basic question. So, I read from the readme in the Ensembl ftp for Hg38 reference and it seems there are several type of file of dna, which is only dna, dna_sm, and dna_rm. If I want to align some fastq file to Ensembl, to which kind of file I should choose? I notice that I should merge all of the fasta file which is per chromosome based fasta into 1 single fasta file and then create index and then do the alignment process. Should I merge all of the file that I download from Ensembl or I should choose only what I need? If I need to choose, how to choose what I need based on my fastq. Thank you for your help and explanation. By the way, I already merge all of the fasta file into 1 single file (including the one with PATCH name in it) and currently I'm aligning my data.

Answer: The answer is to use the primary assembly. ftp://ftp.ensembl.org/pub/release-79/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
Question: I have generated a gVCF (genomic VCF) with GATK using this command line:

gatk HaplotypeCaller -R GRCh38.p13.genome.fa.gz -I bam.bam -D dbsnp_146.hg38.vcf.gz -O genome.g.vcf.gz -ERC GVCF --output-mode EMIT_ALL_ACTIVE_SITES

I'm using GRCh38 patch 13, my BAM file is also based on it.

I want to analyze my genotype for rs4633: https://www.snpedia.com/index.php/Rs4633

As I understand from SNPedia, the reference allele is "C". As I have two variants from my mother and my father, I can have either (C;C), (C;T) or (T;T)

In my gVCF I have the corresponding line:

CHROM POS ID REF ALT QUAL FILTER INFO FORMAT bar
chr22 19962712 rs4633 C T, 1117.03 . DB;DP=29;ExcessHet=3.0103;MLEAC=2,0;MLEAF=1.00,0.00;RAW_MQandDP=95148,29

The reference allele in chromosome 22 at position 19962712 is "C". But what exactly does "T," mean in the ALT column? How should I interpret this? How can I extract the two alleles (father/mother) from this line and map them to the ones from SNPedia: (C;C), (C;T) or (T;T)?

To ask the question in a different way: how can I obtain both my alleles from the ALT column?

To add to this confusion. Another example. I want to know my genotype for rs1801131. SNPedia says: https://www.snpedia.com/index.php/rs1801131
Possible alleles: (A;A), (A;C), (C;C)

Why the "letters" are different than rs4633? Is it because of A=T and G=C?

The line in the gVCF is:

CHROM POS ID REF ALT QUAL FILTER INFO FORMAT bar
chr1 11794419 rs1801131 T G, 1193.60 . BaseQRankSum=0.306;DB;DP=54;ExcessHet=3.0103;MLEAC=1,0;MLEAF=0.500,0.00;MQRankSum=-4.833;RAW_MQandDP=151742,54;ReadPosRankSum=1.170 GT:AD:DP:GQ:PGT:PID:PL:PS:SB 0|1:19,35,0:54:99:0|1:11794400_G_A:1201,0,583,1258,688,1945:11794400:12,7,20,15

So reference allele "T" and "alternative" is "G,". So what are my two alleles which I can map to the possible alleles from SNPedia: (A;A), (A;C), (C;C)?

My assumption is that REF column is the allele from GRCh38 and ALT is somehow my own, but there should be two - one from mom and one from dad. What I'm doing wrong?

Answer: Two things to help you:

> But what exactly does "T,<\non_ref>" mean in the ALT column?

You have only made a `g.vcf` file with `HaplotypeCaller`. This is a file that is 'poised' to be genotyped, but that has not yet been performed. You need to use the `g.vcf` file as input in the genotyping step of GATK using the [GenotypeGVCFs][1] tool. This will give you the genotypes that you are expecting to see.

> Why the "letters" are different than rs4633? Is it because of A=T and G=C?

Yes. the MTHFR gene is transcribed on the reverse strand of DNA and SNPedia is showing the alleles corresponding to the mRNA. In the reference sequence these alleles will be the complimentary alleles. from SNPedia:

> rs1801131 is a SNP in the MTHFR gene, representing an A>C mutation at mRNA position 1298, resulting in a glu429-to-ala (E429A) substitution (hence this SNP is also known as A1298C or E429A).

[1]: https://gatk.broadinstitute.org/hc/en-us/articles/360036351612-GenotypeGVCFs
Question: I'm exploring the results from Cuffdiff:

1. `.gene_differential_expression_testing`
2. `.transcript_differential_expression_testing`

On some genes I see multiple entries in the 2) and with opposite log2 direction as:

```
TMEM51 chr1:15479027-15546974 A B OK 0.000151584 0.563561 11.8602
TMEM51 chr1:15479027-15546974 A B OK 0.741354 3.92E-05 -14.2062
TMEM51 chr1:15479027-15546974 A B OK 2.39194 0.460979 -2.37541
```

Moreover the gene TMEM51 is missing in 1)......is that normal?

In the 2) are these multiple gene names are different isoforms, if yes how to know if the isoform is new or known one?

Answer: On your transcripts file, everything should have a TCONS ID, if there are two entries with the same gene name, or gene ID, then chances are they'll have a different TCONS ID, implying they're different isoforms of the gene.

If you're seeing the same gene name in the Genes differential expression file from cuffdiff, then that's more concerning, as everything should be collapsed down to an XLOC ID, or essentially a locus that encapsulates the gene.

If you want to know if it's 'novel' under the tuxedo method of detection, you should look at the class code that it's been given, better yet, load up some of the alignments in IGV and judge for yourself. In my experience, I'd take what cufflinks is calling 'novel' with a bucket of salt.
Question: I am still new to bioinformatics and I have not yet fully understood the definition of contig. I have read a few explanations and what I understand is that contigs are fragments of the genome for which we are certain that the order of the bases is correct. Then, we make scaffolds out of the contigs and the goal is to get one scaffold to represent the entire genome.

Right now, I am trying to obtain the full reference genome in FASTA format of Streptococcus pneumoniae BR1064. I found this at ENA and in the top right category under "Send Feedback" it appears "Genome Representation: full". From there, one can get over to the assembly contig and there are 245 contigs. Can I just put all this contigs together and obtain the full genome of the organism? If so, is there a particular way to do it? Should it just be in increasing numerical order?

Answer: You're right that contigs are just fragments of the genome, and that scaffolding is the next step in assembly. Usually, this is done using genetic maps and SNPs (or other markers) so that the contigs can be anchored along that genetic map. Here's one recent example with the *Brassica oleracea* genome: [The *Brassica oleracea* genome reveals the asymmetrical evolution of polyploid genomes][1]

Looking at the [publication][2] for your genome, it doesn't look like they performed scaffolding, maybe because there is nothing to use as a reference to scaffold against. Therefore, it is likely that the 245 contigs are just numbered by the order they fell out of the assembly program, an order which doesn't reflect the 'real' genome. In that case, I wouldn't concatenate them.

What do you want to do with the contigs? If it's just SNP calling or something like that, I would leave the sequences as contigs.

[1]: http://www.nature.com/ncomms/2014/140523/ncomms4930/full/ncomms4930.html#methods
[2]: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3133277/
Question: In Cancer Gene Census (CGC) dataset, ([cgc link][1], there is a column named Role in cancer. The values of this attributes are 'fusion', 'oncogene', 'oncogene, fusion', 'oncogene, TSG', 'TSG, fusion', 'oncogene, TSG'.
What does the 'fusion' mean? Does this mean a mixture? Or a gene that sometimes acts as OG and sometimes as TSG?

[1]: http://cancer.sanger.ac.uk/census

Answer: I don't usually point to wikipedia, but this is not a bad description: https://en.wikipedia.org/wiki/Fusion_gene
Question: I am looking for *Pisum sativum* whole proteome FASTA file. Where can I download it in Uniprot?

Is there any way I can download it manually and also any R package?

Answer: From [**NCBI**][1]

From [**Ensembl**][2]

UniProt does not seem to have this proteome at the moment

[1]: https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/024/323/335/GCF_024323335.1_CAAS_Psat_ZW6_1.0/GCF_024323335.1_CAAS_Psat_ZW6_1.0_protein.faa.gz
[2]: http://%20https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-55/fasta/pisum_sativum/pep/Pisum_sativum.Pisum_sativum_v1a.pep.all.fa.gz
Question: Is there a way to convert from genome build to Ens.db version? (ex. grch38 -> EnsDb.Hsapiens.v87)

Answer: [A table from Ensembl][1] showing assemblies present in each release.

[1]: https://www.ensembl.org/info/website/archives/assembly.html
Question: I am confuse if padj and FDR is same? I am in the mid of analysis and confused if I have FDR or not. I have seen some explanation but not able to get the answer. please share your views.
Thanks

Answer: The "false-discovery rate" is the fraction of positives that are false positives at a given p-value threshold. It is a property of the threshold, not a property of a gene. So, if you do a 1000 tests, and get 100 positives at p<0.05, your FDR is 50% (as 50 false positives would be expected in 1000 tests at p<0.05).

Thus, technically speaking, it doesn't make sense to say a gene has an FDR. This is why many tools will use the term "adjusted p-value" or "q-value". The adjusted p-value is the FDR your experiment would have if you set the threshold at the p-value for this gene. Thus, in our example above, a gene with a p-value of 0.05 would have a padj/qvalue of 50% because if you set the threshold at 0.05, you would have a 50% FDR.
Question: Hi all,

I downloaded hg19.p13(GCF_000001405.25) fasta files from [here][1]. Release date is June 2013.

I needed gene annotation files for my downstream analysis. When I checked the UCSC website, they have the first release hg19 (GCA_000001405.1). Release date Feb 2009 and all the associated files. If I use UCSC annotation files for the hg19.p13, will that be alright? I think it is good to use the latest release of genomes. Since hg19 has all the associated file we decided to use hg19 instead of hg38. Is an annotation file available for hg19.p13?

My downstream pipeline includes both CNA(bowtie) and RNA (tophat,cufflinks,cuffdiff)analysis.

Any suggestions?

Thank you,
Deeps

[1]: http://www.ncbi.nlm.nih.gov/projects/genome/assembly/grc/human/data/index.shtml

Answer: Those patch numbers (p.13) don't affect the primary assembly. The sequence you actually use does not differ from hg19.p1 to hg19.p13.

That is why it's called a "freeze". Can you imagine the chaos that would ensue if the genomic sequence changed every few months?

https://www.biostars.org/p/45788/
Question: Hello,

I was reading through https://www.biorxiv.org/content/10.1101/2020.04.06.025635v1.full.pdf

which happens to be a pooled sequencing method using barcodes to test about 10000 covid samples at one go.

I came across compressed DNA barcoding space and DNA barcoding deconvolution for the first time and need some help understanding it.

Although, I searched about it on internet i couldn't find any information on it.

I will highly appreciate if someone can help me with explanation.

Thanks.

Answer: The general idea behind pooled sequencing is that we sequence N samples with X barcodes where X << N.

A major cost in NGS is uniquely barcoding each sample. For each unique barcode, a primer of ~ 60 - 90 bases (depending on design) needs to be synthesized and purified, and roughly adds a cost of ~ $1 - 2 per sample. For sequencing a large number of samples at a time to make full use of the sequencing capacity, say 10,000 samples per day like in Covid-19 testing, we need to uniquely barcode each sample so that we can identify each sample post sequencing. So now you can see the problem in terms of cost. Ordering 10,000 barcoded primers is going to cost several hundred thousand dollars and managing the workflow is going to be a non-trivial.

However, if you have add multiple barcodes to each sample, now you can uniquely tag each sample with a smaller set of barcodes. For example, with 10 barcodes and uniquely adding 5 barcodes to each sample, you can individually barcode ~ 30,000 samples (use permutation formula as barcode order also matters n! / (n-r)! ; n = 10, r = 5). Now you have drastically reduced the cost of ordering barcode primers. Sure, you are using more of each barcode primer but ordering a few barcode oligos in bulk is cheaper and makes managing the workflow easier.

Sample 1 gets barcodes B1,B2,B3,B4,B5

Sample 2 gets barcodes B1,B2,B3,B4,B6 and so on.

There are other ways to find the identities of N samples using X barcodes/NGS libraries where X << N, but here, we make pools of samples so that each sample is distributed over a unique set of pools and after sequencing the pool, we solve the sample IDs based on the pools in which the samples occurred. See the following papers for some simple examples.

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6134198/

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5109470/
Question: Hello,

I am trying to convert a batch of BAM files to FASTQs. I started out testing SAMTOOLS (collate/bam2fq) and PICARD (SAMTOFATQ). On the outset the numbers seemed OK but the statistics suggests that the SAMTOOLS out has twice the amount of duplicates as the Picard out.

Has anyone experienced this? I am not sure if it is a samtools problem or I am not comprehending the QC stats.

Any advice/recommendation/comments are welcome.

Thanks!

PS: In both cases I am outputting both first end of the pair and the second end of the pair as separate files.

**UPDATED:**
The commands used were:
##Samtools

samtools collate -o name-collate.bam sample.bam
samtools fastq -1 sample_1.fastq.gz -2 sample_2.fastq.gz -0 sample_0.fastq.gz name-collate.bam

##Picard

java -Xmx2g -jar picard.jar SamToFastq I=sample.bam FASTQ=sample_1p.fastq.gz SECOND_END_FASTQ=sample_2p.fastq.gz UNPAIRED_FASTQ=sample_0p.fastq.gz

##Fasqc check

fastqc -o fastqc_out/ sample_1p.fastq.gz
###Picard QC
![Picard][1]
###Samtools QC
![Samtools][2]

[1]: https://s15.postimg.cc/w6lv1rd8b/picard.png
[2]: https://s15.postimg.cc/i0646jhsr/samtools.png

Answer: By default, picard don't output non-primary alignments, and samtools does. These secondary alignments which `samtools fastq` outputs should have two effects: an increase in duplication rate, as you noticed, and a larger number of reads - can you confirm this?

Probably Picard behavior is what you want. If you read the samtools manual carefully, you will see how to avoid outputting non-primary alignments.
Question: Can anyone point me to a program [C/Python] that runs a Smith-Waterman pairwise local alignment but (unlike EMBOSS/water) has options for one or both of these:

Return the best N hits of query sequence.

Extend the/each local alignment to produce the best full local alignment for the (much smaller) query sequence.

Answer: I like exonerate, which can be run with the `affine:local` and `-n` options to get what you want. It is a nice tool that runs from the command line and does not need to format databases or further configurations

$: exonerate -q querysequence.txt -t targetsequences.txt -m affine:local -n 3 # n=3 to get the 3 best alignments.
Question: Briefly, I am trying to automate the process of aligning a query sequence (from FASTA) to an assembly graph (FASTG from Spades assembly). As output I need the sequences of the paths in the assembly graph corresponding to the alignment(s).

More detail: I have used Spades to assemble the genome from a diploid yeast starting with short reads (WGS sequencing). Using the wonderful program [Bandage](https://rrwick.github.io/Bandage/) I am then able to BLAST a certain query sequence against the assembly graph (FASTG file). Because of the diploid nature of the genome, the result looks like this:

There are two paths corresponding to this BLAST alignment. In Bandage I can select the nodes corresponding to a path, and then export that path's sequence to a FASTA file.

Doing this manually gives me exactly what I want (essentially, haplotypes derived from the assembled genome). However, I would very much like to automate this process. What tools should I be looking into?

Answer: I didn't realize it when I posted the question, but Bandage has a command line mode, including a command `querypaths` that accomplishes precisely the task in question.

First, find out where the Bandage executable lives. On a Mac, you select the Bandage application, do "Show Package Contents", and the executable is in `Contents/MacOS/Bandage`. I will just call this executable `Bandage`.

Then you can run

`Bandage assembly_graph.fastq query_sequence.fasta output_prefix`

and it will produce `output_prefix.tsv` with the exactly the paths desired! What a wonderful program. Note that it works with both `fastg` and `gfa` formats as input.
Question: Dear all,

I have several FASTA nucleotide file, containing millions identifiers and their sequence. Among them, there are some identifier like this, ">AT1G01340|AT1G01340.2 Sequence unavailable". Would you please let me know how I can remove them? looking forward to hearing your helpful commands. thanks

Answer: I have modified the perl script which should work for your cause

```
#!/usr/bin/perl
use strict;
use warnings;

$/="\n>";

while (<>) {
s/>//g;
my ($id, $seq) = split (/\n/, $_);
print ">$_" if ((length $seq) > 10 && $seq !~ "Sequence unavailable");
}
```

hth
Question: **We have encountered strange pattern in bam file that is generated from amplicon sequencing. (Nextera XT, illumina MiSeq).**

As you can see at the middle of bam file, rectangular-shaped coverage is formed. Half of the reads finished at the right side of the rectangular; while other half of reads finished at the left side of rectangular.

What can cause such abnormal coverage distribution? Any structural variant?

![enter image description here][1]
![enter image description here][2]
![enter image description here][3]
![enter image description here][4]

[1]: http://i.imgur.com/Z2zFnj5.png
[2]: http://i.imgur.com/uMrYmHO.png
[3]: http://i.imgur.com/uAqsVEj.png
[4]: http://i.imgur.com/ZVarV4M.png

Answer: **We solved the real cause of this pattern. It's duplication event of 200 bp area in that rectangular area.** We BLAST the unmapped parts of reads at the ends of rectangular area. And we found out that they are perfectly matched to region inside of this area.

![enter image description here][1]

[1]: http://i.imgur.com/u6ZjM9E.png
Question: Hello,
If I execute this in bcftools 1.10.2

clinvar.vcf.gz is https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh37/clinvar.vcf.gz

bcftools query -i "CLNREVSTAT='criteria_provided,_conflicting_interpretations'" -f '%CLNREVSTAT\n' clinvar.vcf.gz | sort -u

I get this results:

criteria_provided,_conflicting_interpretations
criteria_provided,_multiple_submitters,_no_conflicts
criteria_provided,_single_submitter

I was expecting to get only criteria_provided,_conflicting_interpretations. Is there anything that I don't understand.

Many thanks

Answer: I supect it's invoked as an OR operator. Can you please try

"CLNREVSTAT='criteria_provided' && CLNREVSTAT='_conflicting_interpretations' "
Question: Hi, I'm new in this area, so thanks a lot for any help in advance.

I have some fastq files, in which in some lines there are additional quotes " " added to the quality score in the beginning and the end sometimes and I want to remove them now.

For example:

@NGSNJ-086:647:GW2112051649th:1:1101:6506:1016 1:N:0:CTGAAGCT+ATAGCCTT
AAACTAAGTCAATTCTAATACGACTCACTATAGGAGCTCAGCCTTCACTGCTTCTTAAAGATGCGCACACAACACTCTTTACGTATGTACCGGCACCACGGTCGGATCCTAGATCGGAAGAGCACACGTCTGAACTCCAGTCACCTGAAG
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@NGSNJ-086:647:GW2112051649th:1:1101:7428:1078 1:N:0:CTGAAGCT+ATAGCCTT
AAACTAAGTCAATTCTAATACGACTCACTATAGGAGCTCAGCCTTCACTGCGACAAAATTGGCCATCTTTCCGACAAACAACATGCCCCACGGCACCACGGTCGGATCCTAGATCGGAAGAGCACACGTCTGAACTCCAGTCACCTGAAG
+
"FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF"

So I want to remove the " " in the last line, is there any efficient way to do this, thanks a lot

Answer: $ sed '0~4 s/"//g' test.fq