Tweet
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.
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.