SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
samtools: BAQ differs between mpileup and calmd adeirossi Bioinformatics 1 04-12-2018 03:57 PM
VarScan2 somatic output (genotype call) desmo Bioinformatics 1 12-27-2012 08:49 AM
Mpileup output bfantinatti Bioinformatics 7 11-07-2012 03:42 AM
too many indels output by varscan2 mrfox Bioinformatics 2 09-25-2012 09:17 AM

Reply
 
Thread Tools
Old 02-28-2014, 02:47 AM   #1
JohanF
Junior Member
 
Location: Sweden

Join Date: Feb 2014
Posts: 1
Default varscan2 mpileup2snp output differs from mpileup output

I have a problem understanding the output of VarScan mpileup2snp.

I use samtools mpileup to generate the pileup (samtools 0.1.18).

samtools mpileup –q0 –Q0 -r 5:175346000-175347000 -f reference.fasta inputfile.bam > outputfile.pileup

The output for the position I am interested in is:

5 175346528 G 94 .$,$Tt.t,,tTtTT+1T,,T.T+1Tt+1t.t,t..ttt,.t.T+1T,Tttttt,t.,.Ttt+1t,.,,,Tt+1t,t...,T.,..tTt..,,.,,T..T,Tt,,.t+1tt,.TT,,
@>!!@!BB!!!!(@A!A(+B!A!>B!!!A@!A&A!!!!!!A!@AB!!'@@AA)!(@!@?@>!@@>?!!!?@3A@@?!@@!=!!A?A(!?<!!>>


(mpileup output format - http://samtools.sourceforge.net/pileup.shtml)

There are many reads with inserts. Since the BAQ calculations in mpileup are applied by default (downweighting the scores close to indels), the base quality scores are low (! = 0) for most variant supporting reads.
I count 51 reads supporting the reference and 43 reads supporting a variant.
However, there are only 7 reads supporting a variant with a base quality score > 0.

I then use the pileup as input to VarScan mpileup2snp (VarScan 2):

VarScan.jar mpileup2snp outputfile.pileup > mpileup2snp_output.snpcounts

The output for the position I am interested in is now:

5 175346528 G T K:50:29:23:38,98%:5,1348E-9 Pass:12:17:7:16:1E0 0 1 0 0 K:50:29:23:38,98%:5,1348E-9

(mpileup2snp output format - http://varscan.sourceforge.net/using-varscan.html)

There are two things I don’t understand in this output:
1) Why is the total depth (50) lower than the number of reads that are counted (52)?
Total depth of coverage is 50, but the sum of reads supporting reference (29) and reads supporting variant (23) is 52. The sum of the reads counted with reference to strand is also 52 (12+17+7+16=52).
2) Which reads are counted by mpileup2snp?
The default threshold for base qualities (the option –min-avg-qual, Minimum base quality at a position to count a read) is 15. If there are only 7 reads supporting the variant with a base quality >0 in the pileup, why does mpileup2snp count 23 variant supporting reads?

To try to understand this discrepancy I tried the analysis with different base quality thresholds.
With –min-avg-qual 0 I get the following output:
5 175346528 G T K:94:51:43:42,57%:2,0929E-16 Pass:24:27:17:26:1E0 0 1 0 0 K:94:51:43:42,57%:2,0929E-16

I can understand this output since all read counts (total, reference and variant supporting) are the same as the ones I count in the mpileup output.

However, with higher base quality thresholds (as with the default value in the example above) the read counts from mpileup2snp do not seem to coincide with what I am seeing in the mpileup output.

How does varscan mpileup2snp translate the quality scores?
I translate the scores as the ASCII of the character minus 33.
Example: ! = 0, @ = 31, > = 29.

I am thankful for any comments that can help me understand the mpileup2snp output.
JohanF is offline   Reply With Quote
Reply

Tags
mpileup, varscan2

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 10:06 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO