SEQanswers

SEQanswers (http://seqanswers.com/forums/index.php)
-   Bioinformatics (http://seqanswers.com/forums/forumdisplay.php?f=18)
-   -   vcf file - filter based on DP4 (not total) and removal alternative alleles (http://seqanswers.com/forums/showthread.php?t=50702)

clarissaboschi 03-02-2015 04:14 AM

vcf file - filter based on DP4 (not total) and removal alternative alleles
 
Hi all,

I'm using samtools/mpileup to find variants in my resequencing data of multiple samples together. I and got one vcf file with SNPs and Indels.
I am applying many different filtrations but I would like to filter based on strand bias, for example

Chr1_66769 G A 135 PASS DP=259;VDB=9.53273e-06;SGB=-37.5721;RPB=0.988496;MQB=2.5526e-0
6;MQSB=0.0499642;BQB=0.235831;MQ0F=0;ICB=0.0671329;HOB=0.0246914;AC=6;AN=54;DP4=87,156,0,10 ...

I would like to remove variants with 0 reads in the forward or reverse strand for the alternative allele. In other words variant supported by both forward and reverse strands.

I did not find any command to do it, only manual. Any suggestion?

Other question is how to remove from vcf file the alternative alleles, I do not want this info in my filtered file.
Example

Chr11_6676982 C G,A .....

I would like to have
Chr11_6676982 C G .....

Thanks very much :D:D
Clarissa

sarvidsson 03-02-2015 04:38 AM

Try vcf-annotate in VCFtools http://vcftools.sourceforge.net/perl...l#vcf-annotate - look under "read even more" to get examples on how to create custom filters for DP4.

clarissaboschi 03-02-2015 10:19 AM

Thanks sarvidsson

I saw this option. In my case I think the one that I need is:

#Only loci with enough reads supporting the variant will pass the filter
{
tag => 'INFO/DP4',
name => 'FewAlts',
desc => 'Too few reads supporting the variant',
apply_to => 'SNPs',
test => sub {
if ( !($MATCH =~ /^([^,]+),([^,]+),([^,]+),(.+)$/) )
{
error("Could not parse INFO/DP4: $CHROM:$POS [$MATCH]");
}
if ( 0.1*($1+$2) > $3+$4 ) { return $PASS; }
return $FAIL;
},
},

But I tried to edit this script to my need and did not worked... The way that the script is, is not what I want :(

lindenb 03-02-2015 02:12 PM

using my tool vcffilterjs: https://github.com/lindenb/jvarkit/wiki/VCFFilterJS

Code:

curl -Ls "https://raw.githubusercontent.com/cs418-familysnps/familysnps/3ab8d34c7a7afd7336a40c8f79a32787dcb01389/SNPPhasingProject/vcf/test/bcftools.vcf" |\
java -jar dist-1.128/vcffilterjs.jar -e 'function accept(v) { if(!v.hasAttribute("DP4")) return false; var tokens=v.getAttribute("DP4"); return tokens.get(2)>0 && tokens.get(3)>0;} accept(variant);'

output:

Code:

#CHROM        POS        ID        REF        ALT        QUAL        FILTER        INFO        FORMAT        sample
MT        321        .        T        .        283        .        AC1=0;AF1=0;DP=289;DP4=270,3,1,1;FQ=-282;MQ=56;PV4=0.029,0.00028,0.12,0.0089;VDB=0.0302        PL        0
MT        344        .        T        .        283        .        AC1=0;AF1=0;DP=392;DP4=319,49,1,1;FQ=-282;MQ=55;PV4=0.25,6.8e-05,0.16,1;VDB=0.0394        PL        0
MT        354        .        C        .        283        .        AC1=0;AF1=0;DP=481;DP4=328,104,1,1;FQ=-282;MQ=53;PV4=0.43,0.00091,0.012,1;VDB=0.0447        PL        0
MT        360        .        A        .        283        .        AC1=0;AF1=0;DP=515;DP4=355,125,1,1;FQ=-282;MQ=52;PV4=0.45,2.4e-06,0.29,0.0068;VDB=0.0090        PL        0
MT        366        .        G        .        283        .        AC1=0;AF1=0;DP=553;DP4=385,146,1,2;FQ=-282;MQ=52;PV4=0.19,5.9e-09,0.12,1;VDB=0.0429        PL        0
MT        375        .        C        .        283        .        AC1=0;AF1=0;DP=620;DP4=406,186,7,4;FQ=-282;MQ=51;PV4=0.75,1,0.49,1;VDB=0.0394        PL        0
MT        405        .        T        .        283        .        AC1=0;AF1=0;DP=766;DP4=479,272,1,1;FQ=-282;MQ=51;PV4=1,2.2e-21,1,0.19;VDB=0.0227        PL        0
MT        416        .        T        .        283        .        AC1=0;AF1=0;DP=788;DP4=473,306,1,1;FQ=-282;MQ=52;PV4=1,5.1e-12,0.31,1;VDB=0.0205        PL        0
MT        434        .        C        .        283        .        AC1=0;AF1=0;DP=779;DP4=442,285,3,1;FQ=-282;MQ=53;PV4=1,8.7e-16,1,0.0067;VDB=0.0367        PL        0
MT        450        .        T        .        283        .        AC1=0;AF1=0;DP=758;DP4=446,248,1,1;FQ=-282;MQ=56;PV4=1,0.00099,0.031,0.3;VDB=0.0242        PL        0


clarissaboschi 03-03-2015 03:15 AM

Thanks lindenb, I will try your script

Your script does maintain the header of the vcf? Because I need.

Actually I want the alternative allele (SNP allele), but sometimes there are different alleles for it, and the first one is the more frequent, so I would like to have only the more frequent alt allele, for example

Chr11_667698 ....C G,A .....
I would like to have
Chr11_667698 ....C G .....

thanks
Clarissa


All times are GMT -8. The time now is 08:49 AM.

Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.