SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
VarScan output to VCF OlgaMikh Bioinformatics 4 06-07-2012 11:46 PM
varscan output jgSoton Bioinformatics 1 12-20-2011 08:39 AM
QC in GATK VCF output jorge Bioinformatics 2 07-26-2011 10:34 PM
Cufflinks and CuffDiff bugs? lewewoo Bioinformatics 7 07-14-2011 11:33 PM

Reply
 
Thread Tools
Old 05-12-2012, 01:47 AM   #1
skilpinen
Junior Member
 
Location: Helsinki, Finland

Join Date: Sep 2011
Posts: 3
Default Varscan VCF output bugs(?)

Hello,

This is both observation and question does anybody know (Varscan developers?) is this really a bug or something else.

By using latest Varscan 2.2.11 and mpileup2snp command. Mpileup file from samtools with phred+33 scale quality values produces following VCF file:

##fileformat=VCFv4.0
##source=VarScan2
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##FILTER=<ID=str10,Description="Less than 10% or more than 90% of variant supporting reads on one strand">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample1
1 565286 C T . PASS DP=13 GT:GQP 1/1:7:13
1 569492 T C . PASS DP=45 GT:GQP 1/1:26:45


There are two problems:

i) First of all ID column is empty, and I mean that even though header line has "ID" the following lines do not contain that at all => messes up mandatory 8 fields of VCF file => does not work in combination with Annovar (convert2annovar.pl). I managed to get convert2annovar.pl to accept the vcf file from Varscan only by running it through following (basicly add "null" to each row of ID column):

perl -nle 'if ($_ !~ /^#/){@tmp=split(/\t/,$_);$,="\t";print @tmp[0..1],"null",@tmp[2..(scalar(@tmp)-1)]} else {print}' result.vcf > result_fixed.vcf

Does anybody know some other solution or reason for this output?

ii) QUAL column contains dot for every variant. Shouldn't there be a number? (This doesn't seem bother Annovar, but it is potentially problematic bug(?)

Last edited by skilpinen; 05-12-2012 at 01:51 AM.
skilpinen is offline   Reply With Quote
Old 05-12-2012, 01:58 AM   #2
skilpinen
Junior Member
 
Location: Helsinki, Finland

Join Date: Sep 2011
Posts: 3
Default Update on the issue

Interestingly, apparently data was not phred+33 even though Tophat claimed so... conversion to phred+33 and redo everything:


##fileformat=VCFv4.0
##source=VarScan2
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##FILTER=<ID=str10,Description="Less than 10% or more than 90% of variant supporting reads on one strand">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample1
1 565286 . C T PASS DP=13 GT:GQP 1/1:7:13
1 569492 . T C PASS DP=45 GT:GQP 1/1:26:45
1 22974297 . C A PASS DP=19 GT:GQP 1/1:10:19


Now Varscan produces dot in the ID column (convert2annovar.pl is happy with this), but dot disappeared from QUAL column (which is empty)...?
skilpinen is offline   Reply With Quote
Old 05-15-2012, 01:26 PM   #3
dkoboldt
Member
 
Location: St. Louis

Join Date: Mar 2009
Posts: 62
Default

Hello, I'm the developer of VarScan and this appears to be a bug affecting mpileup2snp, mpileup2indel, and mpileup2cns but not the somatic functions. The QUAL field should contain an empty character (".") but is skipped altogether.

I've already put in a fix and will push that out as soon as possible.

At present, the QUAL value will remain an empty value for multi-sample calling in VarScan. This meets VCF specification, and it seems most appropriate since it's difficult to assess the overall quality of a variant called across multiple samples.

If you have further questions, feel free to e-mail me directly: dkoboldt (at) genome (dot) wustl (dot) edu.
dkoboldt is offline   Reply With Quote
Old 09-04-2012, 06:36 AM   #4
tommivat
Junior Member
 
Location: Finland

Join Date: Aug 2012
Posts: 4
Default Further output problems

Hi,

I was processing varscan vcf output with GATK and got the following error:

Code:
##### ERROR MESSAGE: The provided VCF file is malformed at approximately line number 1479: Unparsable vcf record with allele T/C
Apparently, T/C is not valid vcf. Any plans of correcting this in the near future?

br,
Tommi
tommivat is offline   Reply With Quote
Old 09-13-2012, 08:04 AM   #5
oliviajm
Member
 
Location: france

Join Date: Apr 2012
Posts: 13
Default

Hi,

I'm wondering if it's valid vcf too.
Actually, the vcf format in the varscan output looks like

chr1 14454194 . C G . str10 ADP=11;WT=0;HET=0;HOM=1;NC=0 GT:GQ:SDPP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR 1/1:5:11:11:0:11:100%:1,4176E-6:0:25:0:0:11:0
chr1 16300817 . A -T . PASS ADP=14;WT=0;HET=1;HOM=0;NC=0 GT:GQ:SDPP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR 0/1:4:14:14:4:10:71,43%:7,6278E-5:40:40:4:0:10:0
chr1 16459021 . G +CA . PASS ADP=13;WT=0;HET=1;HOM=0;NC=0 GT:GQ:SDPP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR 0/1:3:16:13:3:8:53,33%:5,16E-4:40:30:3:0:0:8

But on the vcf format reference page on the 1000genome website (http://www.1000genomes.org/wiki/Anal...mat-version-41) there is no "+CA" or "-T" in the fourth column, but only letters (sometimes separated by a comma, but never with a + or a -).

I'm using annovar after varscan and I get some weird results (no insertion or deletion) and I'm wondering if this strange behaviour could come from this format differences. Does anyone have met the same kind of problem?


Olivia
oliviajm is offline   Reply With Quote
Old 01-29-2013, 04:07 PM   #6
thondeboer
Member
 
Location: Bay Area

Join Date: Jan 2009
Posts: 24
Default Problem with VCF file produced by VarScan

I also see the problem with the + and - values for the ALT allele in the VCF produced by VarScan...

I also see "str10' appear in some of the Qual fields as well (See Below).
This does not always appear with the InDels....

i think the +A is to try to signify an insertion, but that needs to be defined as G GA if the A is an insertion after the G...

Conversely, a deletion always includes the base that is left

So: GATC G is a deletion of the 3 bases AFTER the G...

Code:
1	228405118	.	G	+A	.	str10	ADP=29;WT=0;HET=0;HOM=1;NC=0	GT:GQ:SDP:DP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR	1/1:0:29:29:1:28:96.55%:9.8E-1:28:30:0:1:0:28
1	237730169	.	C	T	.	str10	ADP=58;WT=0;HET=0;HOM=1;NC=0	GT:GQ:SDP:DP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR	1/1:0:67:58:0:58:100%:9.8E-1:0:17:0:0:3:55
Thanks,

Thon
__________________
Thon
__________________________________
Thon de Boer, Ph.D.
Director of Product Management, Software
Strand Life Sciences
548 Market Street, Suite 82804
San Francisco, CA 94104, USA
[email protected]
www.strandls.com
Pioneers in Discovery Research Informatics
_______________________________________
thondeboer is offline   Reply With Quote
Old 01-29-2013, 11:08 PM   #7
oliviajm
Member
 
Location: france

Join Date: Apr 2012
Posts: 13
Default

Hi Thon,

I've solved the problem by writing a small script to convert the files obtained with VarScan in the "right" vcf format. By this way, it's OK with Annovar.

Concerning the str10 values that you see, they're not in the Qual field, but in the Filter one and mean that this variation would not pass the filter "Less than 10% or more than 90% of variant supporting reads on one strand".
And the Qual field is missing. There should be a "." in this column.

Hope it will help you.

Olivia
oliviajm is offline   Reply With Quote
Old 01-30-2013, 09:30 AM   #8
dkoboldt
Member
 
Location: St. Louis

Join Date: Mar 2009
Posts: 62
Default

Hello all,

Thank you for following up on this issue -- we certainly do want VarScan's VCF output to meet the specification. I believe that I've already addressed the missing QUAL field, and the way that indels are represented, but I'll need to verify that the alleles-separated-by-slash issue was addressed so that they're separated by commas in VCF output. There will almost certainly be a new release soon that addresses these issues.

Olivia, would you mind sending me your script to dkoboldt (at) genome [dot] wustl [dot] edu ? I could use that to determine if more changes are necessary to meet VCF compliance.
dkoboldt is offline   Reply With Quote
Old 01-31-2013, 09:01 AM   #9
tommivat
Junior Member
 
Location: Finland

Join Date: Aug 2012
Posts: 4
Default

Hello Dan and others,

First, thanks for the great piece of software! It would space us some work if somaticFilter supported .vcf files as well. I don't know if it is tedious to implement.

Another thing I wanted to ask, not related to vcf, concerns false-positive filtering (fpfilter.pl). I'm using bam-readcount to produce input for the script, but even if I do it chromosome by chromosome, the files are too big (>50G) and my computer (with 8Gb memory) just gets jammed when running the fpfilter.pl. Is there a way to do modify the script to support pipeing? And please tell me if it already does.

br,
Tommi
tommivat is offline   Reply With Quote
Old 02-01-2013, 04:00 AM   #10
oliviajm
Member
 
Location: france

Join Date: Apr 2012
Posts: 13
Default

Hello all,

I realised that the missing qual field had been added in one of the last versions of VarScan. As I use it in a pipeline, I did not update it recently to avoid compatibility problems.

But after a few tests, it seems to me that the insertion and deletion are still coded with + and - in the ref and alt column, which don't match with the vcf specifications. I think an insertion of a T after a C should be written C in the ref field and CT in the alt field (and not by +T) for example.

Regards,

Olivia
oliviajm is offline   Reply With Quote
Old 02-01-2013, 04:19 AM   #11
tommivat
Junior Member
 
Location: Finland

Join Date: Aug 2012
Posts: 4
Default

Olivia, +/- issue is fixed in the latest version 2.3.4.

However, the way variant alleles are coded is still unconventional. vcf format uses comma to separate alleles whereas varscan uses slash so I hope this can be fixed in future releases:

Code:
A/C -> A,C
ACG/CG -> ACG,CG
br,
Tommi
tommivat is offline   Reply With Quote
Old 02-01-2013, 06:11 AM   #12
oliviajm
Member
 
Location: france

Join Date: Apr 2012
Posts: 13
Default

Is it fixed in all VarScan tools?

Which one do you use?

I use VarScan.v2.3.4.jar mpileup2indel and I still get some "C +AAAG" or "G -AT" in my vcf output file.


Olivia
oliviajm is offline   Reply With Quote
Old 02-01-2013, 06:17 AM   #13
tommivat
Junior Member
 
Location: Finland

Join Date: Aug 2012
Posts: 4
Default

Quote:
Originally Posted by oliviajm View Post
Is it fixed in all VarScan tools?
Which one do you use?

I use VarScan.v2.3.4.jar mpileup2indel and I still get some "C +AAAG" or "G -AT" in my vcf output file.
That explains.. I use somatic for tumor-normal pairs.

Tommi
tommivat is offline   Reply With Quote
Old 02-05-2013, 02:38 AM   #14
NestorNotabilis
Member
 
Location: Cardiff

Join Date: Dec 2011
Posts: 19
Default

Quote:
Originally Posted by tommivat View Post
Olivia, +/- issue is fixed in the latest version 2.3.4.
To add to Olivia's comment, I'm also using 2.3.4 and still getting the +/- issue when using mpileup2snp.
NestorNotabilis is offline   Reply With Quote
Old 05-15-2013, 02:55 AM   #15
sophiespo
Member
 
Location: australia

Join Date: Apr 2013
Posts: 15
Default

Quote:
Originally Posted by NestorNotabilis View Post
To add to Olivia's comment, I'm also using 2.3.4 and still getting the +/- issue when using mpileup2snp.
I am as well when using the somatic/processSomatic functions.

Annovar doesn't like this.. can anyone help?
sophiespo is offline   Reply With Quote
Old 05-15-2013, 10:32 AM   #16
eeyun
Junior Member
 
Location: 90049

Join Date: Jan 2013
Posts: 9
Default

We are having the same problem with 2.3.5

<pre>#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample1
chr1 6529182 . TTCC TCC . PASS ADP=314;WT=0;HET=1;HOM=0;NC=0 GT:GQ:SDPP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR 0/1:255:322:314:178:138:43.67%:1.1101E-50:34:31:88:90:70:68</pre>
eeyun is offline   Reply With Quote
Old 05-15-2013, 10:38 AM   #17
eeyun
Junior Member
 
Location: 90049

Join Date: Jan 2013
Posts: 9
Default

Quote:
Originally Posted by eeyun View Post
We are having the same problem with 2.3.5

<pre>#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample1
chr1 6529182 . TTCC TCC . PASS ADP=314;WT=0;HET=1;HOM=0;NC=0 GT:GQ:SDPP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR 0/1:255:322:314:178:138:43.67%:1.1101E-50:34:31:88:90:70:68</pre>
As far as I can tell, it should be ref = TTCC and alt = T
eeyun is offline   Reply With Quote
Old 05-15-2013, 10:39 AM   #18
eeyun
Junior Member
 
Location: 90049

Join Date: Jan 2013
Posts: 9
Default

Quote:
Originally Posted by eeyun View Post
As far as I can tell, it should be ref = TTCC and alt = T
Attachment included here to show the variant in question.
Attached Images
File Type: png varscan problem.PNG (1.8 KB, 19 views)
eeyun is offline   Reply With Quote
Old 07-25-2013, 02:17 PM   #19
IsmailM
Junior Member
 
Location: Los Angeles

Join Date: Apr 2013
Posts: 7
Default solved indel vcf format with awk command

Here is an awk command that can change your indel vcf format into the correct format.

cat Original_VCF | awk 'BEGIN {OFS="\t"} NR <= 24' > FINAL_VCF && cat Original_VCF | awk 'BEGIN {OFS="\t"} NR >= 25 { if (length($4)>length($5)) {$5 = substr($4, 0, 1)}; print }' >> FINAL_VCF


It uses two awk commands because the second command changes the header of the file if you run it on the whole file. So the first awk command transfers the header(assumed to be 24 lines) and then from the 25th line down is the vcf indels that are changed to the correct indel format using the second awk command.
IsmailM is offline   Reply With Quote
Old 12-23-2013, 07:33 PM   #20
rnahar
Junior Member
 
Location: Singapore

Join Date: Aug 2013
Posts: 6
Default

I am also facing the +/- issue in the varscan indel notations - however i do not use the vcf output but prefer the regular tabular output of Varscan. Is there a way that this indel notation can be changed so as to be compatible with annovar ? I use Varscan 2.3.6
rnahar is offline   Reply With Quote
Reply

Tags
annovar, varscan, vcf

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 11:50 PM.


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