SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Varscan v2.3.5 headers missing wdemos Bioinformatics 11 04-28-2015 09:29 PM
How to generate a figure like this? gene_x Bioinformatics 2 03-07-2013 02:23 PM
VarScan mpileup2snp --p-value affects genotype call Graham Etherington Bioinformatics 2 10-22-2012 01:41 PM
how to get quality information from VarScan mpileup2snp alg Bioinformatics 0 07-10-2012 06:18 PM
First Call Missing on SOLiD johnadam33 Bioinformatics 9 09-20-2011 08:44 PM

Reply
 
Thread Tools
Old 08-29-2013, 12:19 PM   #1
eeyun
Junior Member
 
Location: 90049

Join Date: Jan 2013
Posts: 9
Default VarScan mpileup2snp missing a call - can't figure out why

VarScan mpileup to snp seems to be missing a call

the mpileup file looks like this...

Code:
chr5	13894894	T		40	.a.AAAAAaAA,aaaaA,,aAA,AaaAaAaAaAAa,.,.,	B02/8/0///0G/3/10GH///C/5//0511//10HHH3A	55	a,....Aa,.AAa,aAAA.a,,,,...,,a..,,A,.,.,AaAAaA.aaA.,,,a	!DHHHH0/HH1//H////G/HHGHHGFHH/GGHH/HHHFC////1/G/00FHHH)
Then, I'm using the following options in VarScan
Code:
mpileup2snp --min-coverage 20 --min-var-freq 0.2 --min-reads2 4 --strand-filter 1 --output-vcf 1 >whatever.vcf
And this is the line in the VCF where the problem is...

Code:
chr5	13894894	.	T	A	.	PASS	ADP=32;WT=1;HET=2;HOM=0;NC=1	GT:GQ:SDP:DP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR		0/1:53:40:25:11:14:56%:4.7528E-6:33:16:4:7:7:7		0/0:32:55:37:32:5:13.51%:2.706E-2:38:15:15:17:3:2
Attached is an image from IVG of the two alignments (front to back replicates of the same specimen, variant is clearly in both, but not in VCF)

Any thoughts? Apologies is the answer is obvious, but I can't seem to figure it out.

Thanks in advance
Attached Images
File Type: png MissedCall.PNG (13.2 KB, 14 views)
eeyun is offline   Reply With Quote
Old 09-03-2013, 09:18 AM   #2
dkoboldt
Member
 
Location: St. Louis

Join Date: Mar 2009
Posts: 62
Default

Thank you for posting. From the VarScan VCF entry for sample 2:
0/0:32:55:37:32:5:13.51%:2.706E-2:38:15:15:17:3:2

You can see that the VAF as computed by VarScan is 13.51%, which is below your minimum threshold of 20%.

Notably, VarScan depth (DP, field 4) is 37 whereas the raw SAMtools depth (SDP, field 3) is 55. This suggests that the raw pileup and IGV are showing about 18 reads whose base qualities are below VarScan's threshold.
dkoboldt is offline   Reply With Quote
Old 09-03-2013, 09:36 AM   #3
eeyun
Junior Member
 
Location: 90049

Join Date: Jan 2013
Posts: 9
Default

Quote:
Originally Posted by dkoboldt View Post
Thank you for posting. From the VarScan VCF entry for sample 2:
0/0:32:55:37:32:5:13.51%:2.706E-2:38:15:15:17:3:2

You can see that the VAF as computed by VarScan is 13.51%, which is below your minimum threshold of 20%.

Notably, VarScan depth (DP, field 4) is 37 whereas the raw SAMtools depth (SDP, field 3) is 55. This suggests that the raw pileup and IGV are showing about 18 reads whose base qualities are below VarScan's threshold.
Thanks very much for the reply.

I had wondered whether this was a result of the base quality being downgraded by samtools, and therefore not read by varscan. I've tried mpileup with -B and -E, but neither seems to make much difference.

It's not a huge issue as overall concordance between the two specimens is very good, but I like to chase down the "misses". I might try a different quality threshold in varscan and see if that helps at all.
eeyun is offline   Reply With Quote
Old 10-01-2013, 02:59 PM   #4
eeyun
Junior Member
 
Location: 90049

Join Date: Jan 2013
Posts: 9
Default

I'm sorry to re-ask the same question again, but I've noticed another missed variant call using the same above pipeline, and this one I can't figure out based on the base qualities. If anyone can point out what I'm missing, I'd really appreciate it.

mpileup at position
Code:
chr11	89017961	G	26	..A.A,a....,a....A,,.a,.,,	GG5?5H6GEGGG52GG06FHC6BAHH
As far as I can tell, the variant bases have an average quality of around 20 and the reference bases have an average quality of around 30 (am I wrong?). There should be no strand bias issues.

Yet, using the following VarScan parameters does not yield the variant in the vcf.

Code:
--min-coverage 20 --min-var-freq 0.2 --min-reads2 4 --output-vcf 1 --strand-filter 1
Any help greatly appreciated, I'm really lost on this one.
eeyun is offline   Reply With Quote
Reply

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 01:48 AM.


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