cannin/biostars_qa huggingface dataset:
Question: When I use BWA(version 0.7.15-r1140), I select the BWA-MEM algorithm and here is my command:

`bwa mem -t 2 reference_file fastq1.fq.gz fastq2.fq.gz > result1.sam`

This works well and then I use this command to improve the threads:

`bwa mem -t 10 reference_file fastq1.fq.gz fastq2.fq.gz > result2.sam`

This also works well but when I compare the result1.sam with result2.sam, they are different! I also test with `-t 6, -t 16`, and all the results are different. However, when I run with the same threads twice, the results are identical. So I found that BWA-MEM will get different results with different threads.

Then I read the source code and found this:

kt_for(opt->n_threads, worker1, &w, (opt->flag&MEM_F_PE)? n>>1 : n); // find mapping positions
for (i = 0; i < opt->n_threads; ++i) smem_aux_destroy(w.aux[i]);
free(w.aux);
if (opt->flag&MEM_F_PE) { // infer insert sizes if not provided
if (pes0) memcpy(pes, pes0, 4 * sizeof(mem_pestat_t));
else mem_pestat(opt, bns->l_pac, n, w.regs, pes);
}
kt_for(opt->n_threads, worker2, &w, (opt->flag&MEM_F_PE)? n>>1 : n);

That is, BWA-MEM use n_threads(such as -t 6, n = 6) to find mapping positions, but use only 1 thread to execute the function `mem_pestat` to calculate `avg`(average of the insert size) and `std`(standard deviation of the insert size), which are important to find pair information. According to BWA, every thread will process around 10000000bp, so:

If I use -t 2, BWA will calculate avg and std with 2 x 10000000bp

If I use -t 10, BWA will calculate avg and std with 10 x 10000000bp

If I use -t 16, BWA will calculate avg and std with 16 x 10000000bp

So I know why the results are different with different threads.

I wonder to know if there is anything wrong with my opinion? If it's correct, I want to know how to evaluate the difference? The difference will change which filed of the SAM record( such as RNAME or POS)? If it's wrong, I want to know the real reason to make the difference.

Any reply will be much appreciated!

Answer: It's almost guaranteed that the differences are in alignments with low (0 or so) MAPQs. I have [a little C program](https://github.com/dpryan79/Misc/blob/master/CompareBAMs.c) that uses htslib to compare two name sorted BAM files (or at the very least they need to be in the same order). It'll print out differences due to chromosome/position/cigar. If you then looked at a couple of those my guess would be that the differences are only occurring at low confidence alignments, which is wholly unsurprising.
Question: Dear All,

As per my understanding from STAR manual, I am about to run a STAR 2.7.0f mapping pipeline with 2-pass mode for multiple samples of patiets of diseases and healthy peoples as follows:

Could you please help me to validate all the commands I am running correctly or do you have any suggestions?

**1) Indexing genome with annotations**

STAR --runMode genomeGenerate --genomeDir ~/db/hg38/ --genomeFastaFiles ~/db/hg38/hg38.fa --sjdbGTFfile ~/db/hg38/hg38.gtf --runThreadN 30 --sjdbOverhang 89

**Note:**

- Indexing for maximum read length **90 bp**.

**2) 1-pass mapping with indexed genome**

STAR --genomeDir ~/db/hg38/ --readFilesIn sample1.R1.fastq.gz sample1.R2.fastq.gz --readFilesCommand zcat --outSAMunmapped Within --outFileNamePrefix sample1. --runThreadN 30

**Notes:**

- The same command has been run for multiple samples in the for loop, therefore, it will generate **SJ.out.tab** file for each sample.

- Next, I have copied SJ.out.tab files of all the samples into a single folder "**SJ_out**"

**3) Indexing genome with annotations and SJ.out.tab files**

STAR --runMode genomeGenerate --genomeDir ~/db/hg38/SJ_Index/ --genomeFastaFiles ~/db/hg38/SJ_Index/hg38.fa --sjdbGTFfile ~/db/hg38/SJ_Index/hg38.gtf --runThreadN 30 --sjdbOverhang 89 --sjdbFileChrStartEnd SJ_out/*.SJ.out.tab

**Note:**

- Again indexing for maximum read length **90 bp**.

**4) 2-pass mapping with new indexed genome with annotations and SJ.out.tab files**

STAR --genomeDir ~/db/hg38/SJ_Index/ --readFilesIn sample1.R1.fastq.gz sample1.R2.fastq.gz --readFilesCommand zcat --outSAMunmapped Within --outFileNamePrefix sample1. --runThreadN 30

**Notes:**

- Again, the same command has been run for multiple samples in the for loop, therefore, it will generate mapping files for each sample.

Answer: **SOLVED!** Please look at my [recent conversations][1] with Alex regarding the issue.

[1]: https://groups.google.com/forum/#!topic/rna-star/Cpsf-_rLK9I
Question: Please let me know any tools or R packages that can create a mutation diagram showing mutations in protein domains like this figure from the MSKCC cBio Cancer Genomic Portal? Thanks in advance

Answer: Bioconductor package trackViewer gives beautiful lolliplot.
Question: I'm trying to get some `entrez` ids from some gene names using [biomaRt][1]

In my example below I have 2 genes `Igha` and `Mlc1` for mus musculus

My version of biomaRt is `biomaRt_2.38.0`

library("biomaRt")
mart <- useMart(biomart = "ensembl", dataset = "mmusculus_gene_ensembl")
getBM(attributes=c("ensembl_gene_id","mgi_symbol","entrezgene"), values = c("Igha", "Mlc1"), bmHeader = T, filters = "mgi_symbol", mart = mart)
# Gene stable ID MGI symbol NCBI gene ID
#1 ENSMUSG00000095079 Igha NA
#2 ENSMUSG00000035805 Mlc1 170790

I know that some Ensembl gene ID do not have Entrez gene ID https://www.biostars.org/p/275619/

But in this case, I'm able to find the Entrez gene ID from Ensembl

http://www.ensembl.org/Mus_musculus/Gene/Summary?db=core;g=ENSMUSG00000095079;r=12:113254830-113260236

Which is `Refseq Gene ID 238447`

Maybe https://www.biostars.org/u/7019/ will have an answer :)

[1]: https://bioconductor.org/packages/release/bioc/html/biomaRt.html

Answer: The gene listed on the gene page is not actually a proper external reference and doesn't come through BioMart. That comes through a different pipeline and it's not reliable and we're getting rid of it.
Question: Hi all,

I have EPIC array data analysed, but colleagues are concerned about annotation which is massively incomplete: 56% of probes do not have UCSC 'relation to CpG island' data, and 25% no gene annotated. This seems crazy to me, but is confirmed by Illumina. They say a lot of users go and annotate from public resources.

I presume others have come across this? Any chance anyone has a full(er) annotation? Or even ideas on where to source/create?

Thanks, Bruce.

Answer: The annotation that comes with [Minfi][1]'s workflow is generally pretty complete...but as with anything, it has gaps. The manifest that Minfi uses is [here][2]. Check out the annotation there and see if it's as incomplete as Illumina have suggested to you.

[1]: http://bioconductor.org/packages/release/bioc/html/minfi.html
[2]: https://bioconductor.org/packages/release/data/annotation/html/IlluminaHumanMethylationEPICmanifest.html
Question: Hi!

I'm parsing a series of gbff taken from NCBI. However, some of them don't have the DNA sequence inside (e.g.: see the .gbff.gz [here][1]).

What I want to do is to skip these files, but I dont know what to call on Bio::SeqIO. Could anyone just tell me what should I call?

A snippet of my code is:

```
my $in=Bio::SeqIO->new(-file=>"$mygbkff",-format=>'genbank',-alphabet=>"dna");

while (my $seq=$in->next_seq){
for my $feat ($seq -> get_SeqFeatures){
next unless (($feat->primary_tag eq "CDS") and (defined $feat->spliced_seq->seq)); #this was my attempt at finding only cds that had sequences annotated
print "\>";
print $feat->spliced_seq->seq,"\n";
}
}
```

[1]: ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF_000149425.1_ASM14942v1

Answer: Adding this as an answer as requested.

The solution was, rather than verifying if `$feat` is defined, to verify if `$seq` was defined:

```
my $in=Bio::SeqIO->new(-file=>"$mygbkff",-format=>'genbank',-alphabet=>"dna");

while (my $seq=$in->next_seq){
next unless defined $seq->seq;#This ignores those sequences whose DNA was missing
for my $feat ($seq -> get_SeqFeatures){
next unless (($feat->primary_tag eq "CDS") and (defined $feat->spliced_seq->seq)); #this was my attempt at finding only cds that had sequences annotated
print "\>";
print $feat->spliced_seq->seq,"\n";
}
}
```
Question: Hi, I am learning about vcf file formatting and sorry for asking a dumb question.

I have a multi-sample (>300) vcf file and I want to **group** (take) **samples with shared mutations/GT**.

Can someone suggest a tool/command/script solve this problem?

Thanks

Answer: I dont't know if this is exactly what you want.

For each genotype you are interested run this:

$ bcftools query -i 'GT="0/1"' -f '%CHROM\t%POS\t0/1\t%REF\t%ALT\t[%SAMPLE,]\n' input.vcf|sed 's/,$//'|awk -v OFS="\t" '{counts = split($6, samples, ","); print $1, $2, $3, counts, $4, $5, $6}' > stats_01.csv

In this example `bcftools query` looks for samples with genotype `0/1` and print out CHROM, POS, `0/1`, REF, ALT and a comma seperated list of sample names.

`sed` removes the additional `,` at the end. `awk` counts the number of sample names in each row and add this column.

fin swimmer
Question: Hi there,
I've constructed the tree with bootstrap supports in `Ape` package. But i'd prefer to get bootstrap values in percentage not in absolute number of replications. How can i do that in `Ape`?

Answer: I solved the issue by the `R` code:

plot(stree)
nodelabels(round(myBoots/1000*100),adj = 0.7)

where myBoots is numeric vector containing my bootstrap values
Question: What is the difference between gene FPKM and transcript FPKM values? Or rather how is FPKM value for gene obtained from FPKM values of transcript. Average is not suitable as only one transcripts may have expression out of several. So what method is used to do so.

Answer: Thak you Seidel for showing the way. I have now been able to find the calculation methods.

The general method for calculation of transcript FPKM is explained [here][1].

As I am using Ballgown tool, I found the answer to my question in [supplementary][2] article in one of their papers [Ballgown bridges the gap between transcriptome assembly and expression analysis][3]. The gene FPKM vale calulation is explained under the heading Expression data.

[1]: http://www.rna-seqblog.com/rpkm-fpkm-and-tpm-clearly-explained/
[2]: https://media.nature.com/original/nature-assets/nbt/journal/v33/n3/extref/nbt.3172-S1.pdf
[3]: https://www.nature.com/articles/nbt.3172
Question: Hi all

Recently, I deal with some human gene data and I want to know if the gene I handled is conserved gene.

Is any database provide human genes' conservation scores?

I need all human genes' conservation scores.

I know UCSC Genome Table Browser provide sequence conservation scores like phastCons and phyloP. But I can't find gene's conservation score.

Does anyone can teach me how to get it?

The format which I want is

```
Gene ID 1 (Ensembl ID or Refseq ID ) | conservation score
Gene ID 2 (Ensembl ID or Refseq ID ) | conservation score
```

Thanks

Answer: I think it depends on exactly what information you want. Is it conservation between human and a particular species? The UCSC conservation scores give you a measure of conservation based on an alignment between human and 45 other species (here is the information along with download links), which is good for a general impression of conservation levels, but less good for a specific question like "how different is this gene in the marmoset?" (for which you'd want a specific pairwise alignment). They also have a subset of scores for primates, and for placental mammals. So, one way to do it would be to download the data, then process it to extract the answer you're looking for (e.g. get the average phastCons score using the exons, or just coding sequence, of each human gene, then use that average as your gene conservation score). What you mean by "gene" would also depend on the question you were trying to answer - for example, introns may or may not be of interest to you.
Question: Samtools can be used to select reads above certain mapping quality.

samtools view -h -b -q 30 aligned.bam -o above.mapQ30.bam

**But, how to select a read below certain mapping quality - all aligned reads below mapQ 30?**

I know it can be done using awk. But, the pipeline gets lengthy and time consuming when first need to convert bam to sam - separate header - use awk for mapQ below 30 - add header - sam file - convert to bam.
Really, its taking a lots of time.

Thanks,

Answer: # Method B - requires samtools version 1.3

**first select the mapped reads only**
`samtools view -b -h -F 4 aligned.bam > mapped_only.bam` # this selects for mapped reads only.

`samtools index mapped_only.bam` # index the bam file

**Now separate the mapped reads above mapQ 40 and below**
`samtools -b -h -U below_mapQ40.bam -q 40 mapped_only.bam -o above_mapQ40.bam` # make sure to use right samtools version (1.3)

`samtools index above_mapQ40.bam` # this is bam file with reads with mapQ 40 and above

`samtools index below_mapQ40.bam` # this is bam file with reads with mapQ below 40.
Question: Hello,

I haven't worked with miRNA mapping..( the specie is mouse) It is the first time to do it.

I have read a lot of posting and web-site to map my microRNA data.

Here is what I follow.

**1) trim the adapter using cutadapt**

When I looked at my fastq file (single-end), read length was 36 bp.. by inspecting the eye-ball test, I could not find my common sequences. (So, I am not sure my fastq file has been already trimmed or not . So I kept both versions... )This is my command... (comes from http://www.ark-genomics.org/events-online-training-eu-training-course/adapter-and-quality-trimming-illumina-data)

```
cutadapt -a CGACAGGTTCAGAGTTCTACAGTCCGACGATC \
-a TACAGTCCGACGATC \
-a ATCTCGTATGCCGTCTTCTGCTTG \
-e 0.1 -O 5 -m 15 \
-o xxx_microRNA_adaprm.fastq xxx_microRNA.fastq
```

So, I have 2 version fastq files. (`xxx_microRNA.fastq`, `xxx_microRNA_adaprm.fastq`)

**2) Local bowtie2 alignment of miRNA data**

(I generated the reference with bowtie-build directly from the hairpin FASTA file downloaded from miRBase.)

```
Bowtie2-build hairpin_mms.fa hairpin_mms.fa
Bowtie2-build mature_mms.fa mature_mms.fa
Bowtie2 -local -N 1 -L 16 -x hairpin_mms.fa -U fastq/xxx_microRNA(_adaprm).fastq -S xxxx.sam
```

The bowtie2 mapping returns the following result.. (for both versions)

I should have done something wrong...

```
2849144 reads; of these:
2849144 (100.00%) were unpaired; of these:
2847350 (99.94%) aligned 0 times
933 (0.03%) aligned exactly 1 time
861 (0.03%) aligned >1 times
0.06% overall alignment rate
```

Could you please help me with this?

Oops, I forgot to convert U to T in hairpin.fa file.. (so stupid about it..) Once I run it, I will confirm the result again.

Answer: Try to clip the adapter of the Illumina HiSeq 2000 miRNA protocol (TGGAATTCTCGGGTGCCAAGGAACTCCAGTCAC)

And always check your read length distribution after clipping. There should be a peak at around 24nt.
Question: I often see the word depth in the manuals of the tools for NGS, what is its meaning ?

thanks.

Answer: I believe it is the same concept as coverage, it might come from a shorthand of saying depth of coverage.
Question: Hello,

My question is - How many mismatches are allowed (on average/or hard limit) in a bowtie 2 alignment before it moves on to find a new home for a read. I know the we can control the number of mismatches allowed in the seed with -N option. And I know that with the -D option, I control the seed extension attempts that "fail".

So I guess my question now is, when does an extension "fail"? The manual says, when an extension does not yield a new best or a new second-best alignment, it fails.

I'm not able to map this definition to a number (of mismatches). I've re-read the method to calculate the alignment score many times now. Would two consecutive SNPs be allowed, or do two SNPs need to have a certain number of matched bases between them so that the overall alignment doesn't "fail"?

I think, in my naive way, I've explained the core issues of my question. I'll edit if it's not entirely clear.

Thanks,
Ayush

Answer: -D and -N aren't directly related to the number of mismatches. An alignment fails, if it does not meet the --score-min criteria.

For global alignment costs are --ma 0, --mm 2-6 - depending on base quality. Score is given by --score-min L,-0.6,-0.6. Given a 100 bp read the --score-min computes to -60. Given only high quality SNPs (--mm 6), an alignment can have up to 10 mismatches and still meet the score threshold. However, with 10 mismatches randomly distributed along the read, you will need sensistive seeding parameters to initialize the alignment in the first place.

To limit bowtie2 to 2 mismatches per 100 bp in global mode, you can set --score-min L,0,-13. To limit reads of arbitrary length to 2 mismatches, you can set --score-min C,-13,0.

Given default gap costs (5,3), this would also allow 2 non-adjacent 1 bp gaps or one 3 bp gap in your alignment.
Question: I have a VCF file containing genotype date for a few thousand SNPs across a few thousand samples. I would like to firstly convert this to a matrix (possibly using the VariantAnnotation package) and then perform a PCA analysis on the samples followed by some sort of clustering algorithm. I have very little experience with any of SNP matrix packages, PCA or clustering algorithms so I was wondering if anyone knew of any good tutorials which may be able to help me.

It is also worth noting that due to the nature of the analysis I am running, the SNP matrix will be extremely sparse. I would therefore also like to get information on the fraction of missing genotypes for each sample and the fraction of missing samples for each SNP - is this possible?

Answer: To pursue PCA you can use [SNPRelate][1] or [TASSEL][2] that has a user-frendly platform.

[1]: http://corearray.sourceforge.net/tutorials/SNPRelate/
[2]: http://www.maizegenetics.net/tassel
Question: I'd be interested in retrieving via Biomart the protein length associated to a transcript, as per screenshot below

![scrennshot][1]

I have retrieved the list of attributes available from biomart but can't seem to find the right field.
Can someone confirm if this info is accessible programatically ? how ?

Thanks,

[1]: /media/images/f933ce1c-b0bb-453f-9f40-48e0117a

Answer: Hi iatz,

You can't retrieve protein length directly from BioMart. However, you can retrieve the CDS length and divide by 3.

![enter image description here][1]

[1]: /media/images/c887b1d6-5ce8-4e04-8274-37098c0b
Question: Hello, everyone

Someone can help-me at to access the tranche repository (Univ. of Michigan) for MS data acquisition?

The URL (https://www.proteomecommons.org/tranche/) does not work.
When I try to access, I'm redirected to another page.

Best regards,

Answer: The Proteome Commons Tranche sadly shut down around 2013 due to lack of funding. It is a good example why important bioinformatics infrastructure needs dedicated funding and should not have to operate out of research grants for individual groups.
Question: Hi All,

I have been using DESeq2 to perform differential gene expression analysis between two groups but I do not have house keeping genes in the data set. Is it still possible to do DESeq calculations without the house keeping genes? and I do not have all the genes expressed across all subjects/samples.

I have the following error message calculating size factors:

> dds <- DESeq(dds)

>estimating size factors

>Error in estimateSizeFactorsForMatrix(counts(object), locfunc = locfunc, :

>every gene contains at least one zero, cannot compute log geometric means

Any help would be appreciated.

Thank you.

Gajender

Answer: DESeq2 by default doesn't use housekeeping genes. It makes the assumption that most genes are not diferentially expressed, and corrects based on the gene with the median expression

But your problem is exactly as it says; the default normalization method calculates the geometric mean of each gene, that is not possible when a gene has a sample with zero counts. Some genes with some zeroes are fine, the normalization function just ignores them, but in your data, every single gene is being excluded for that reason. The default method of normalization will not work with your dataset as is. If you have a few samples with very few reads, omitting them might fix the problem.
Question: Hi All,

I am working with a viral sequence. Upon investigating the freebayes output vcf file I found that there are two positions where variants are reverse complement of the reference.

Reference	527	.	CCCGGGCGTCGGGCGAC	GTCGCCCGACGCCCGGG	0	.	AB=0;ABP=0;AC=0;AF=0;AN=1;AO=1837;CIGAR=2X2M2X2M1X2M2X2M2X;DP=8277;DPB=9767.53;DPRA=0;EPP=3415.19;EPPR=9687.22;GTI=0;LEN=17;MEANALT=29;MQM=40.3887;MQMR=39.4471;NS=1;NUMALT=1;ODDS=65195.1;PAIRED=0.997278;PAIREDR=0.968095;PAO=0;PQA=0;PQR=137846;PRO=4415;QA=0;QR=230712;RO=6394;RPL=63;RPP=3463.56;RPPR=11044.7;RPR=1774;RUN=1;SAF=1831;SAP=3940.06;SAR=6;SRF=6101;SRP=11459.1;SRR=293;TYPE=complex	GT:DP:DPR:RO:QR:AO:QA:GL	0:8277:8277,1837:6394:230712:1837:0:0,-30803.9

and

Reference	5586	.	GTCGCCCGACGCCCGGG	CCCGGGCGTCGGGCGAC	1.19962e-12	.	AB=0;ABP=0;AC=0;AF=0;AN=1;AO=1439;CIGAR=2X2M2X2M1X2M2X2M2X;DP=6715;DPB=8094.82;DPRA=0;EPP=2732.86;EPPR=7945.83;GTI=0;LEN=17;MEANALT=27;MQM=39.3489;MQMR=38.6015;NS=1;NUMALT=1;ODDS=54869.5;PAIRED=0.994441;PAIREDR=0.963359;PAO=0;PQA=0;PQR=125375;PRO=3860;QA=0;QR=188745;RO=5240;RPL=1401;RPP=2806.41;RPPR=9203.98;RPR=38;RUN=1;SAF=9;SAP=3050.08;SAR=1430;SRF=271;SRP=9149.39;SRR=4969;TYPE=complex	GT:DP:DPR:RO:QR:AO:QA:GL	0:6715:6715,1439:5240:188745:1439:0:0,-25889.3

When I visualized the alignment in IGV the variant was present at the specified location. This is visualization with IGV.

![enter image description here][1]

How can I know for sure that,

1. Is this variant a sequencing or mapping artifact or an inversion
variant?
2. Is visualization in IGV a correct method of validating a
variant?

[1]: https://s6.postimg.cc/6jwcjicg1/position527.png

Answer: Hello,

> Is this variant a sequencing or mapping artifact or an inversion
> variant?

this is very likely some kind of artifact. Have a look at the 6th column of your vcf. This is a quality value for the variant site. Your value is 0 or something very close to it. Normally this value should be at least something around 20. Depending on the read depth this value is much higher.

Also freebayes gives you the genotype `0` which means `REF`

> Is visualization in IGV a correct method of validating a variant?

For a quick check this is absolutely fine. In this case you can also see that the bases which doesn't match the ref are faded. This normally means they have bad quality values.

For a real validation, each variant that has an impact on your goal, has to be confirmed by another method like sanger sequencing.

fin swimmer
Question: Hello, everyone.

I want to download paired samples from the TCGA database of RNA-Seq experiments. I've been looking for information on how to download this data and it looks like it's from https://gdac.broadinstitute.org/. Specifically, I'm looking for breast cancer samples, so I looked in the mRNASeq section.

Inside this section there are several files to download (I want the raw counts) but I don't know what the differences are between the files:

illuminahiseq_rnaseqv2-RSEM_genes (MD5)
illuminahiseq_rnaseq-gene_expression (MD5)

I would also like to know how to filter these files to keep the paired samples. I found something about the sample codes in a previous post, but I couldn't access the link to the explanation (https://wiki.nci.nih.gov/display/TCGA/TCGA+barcode).

Could you help me with this problem please?

PS: Suggestions about other databases are welcome!

Answer: In your file manifest that you used to download the data, you will have a UUID. Use that to look up the TCGA barcode for each file via this function: https://www.biostars.org/p/306400/#306517

You can then infer Tumour-Normal pairings by matching on the TCGA barcode. Yet more information on barcodes:

- https://www.biostars.org/p/313063/#313066
- https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables
- https://github.com/kevinblighe/TCGAbarcode/blob/master/TCGAbarcode.pdf

I have already analysed the TCGA BRCA data many times and there are ~111 Tumour-Normal pairs for the RNA-seq data.

Kevin