Parse VCF file

Question

[aj123]

Hello,

I'm trying to parse a VCF file and filter it based on 1000Gp1_AF < 0.05. Any ideas on how to go about this? Thank you!

[Ram]

You may wanna be more careful with the terminology. Parsing is pretty specific process and you're only looking to filter VCF. You've elicited what you wish to do very clearly, so you're good at stating requirements, but the tech community can be a bit touchy about jargon.

Answer 1

[Biomonika (Noolean)]

It's almost always better to use existing tools than to write parser on your own. Take a look here at section ALLELE FILTERING.

[Ram]

I think OP didn't actually meant parse as in semantic parsing. OP might have used it as a synonym for processing.

Answer 2

[Tariq Daouda]

Hi,

I wrote a python parser that you can find here.

Here is the doc on how to use it.

Cheers

Answer 3

[Pierre Lindenbaum]

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

$ curl -s "ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20130502/ALL.wgs.phase3_shapeit2_mvncall_integrated_v5.20130502.sites.vcf.gz" | gunzip -c | java -jar ~/src/jvarkit-git/dist/vcffilterjs.jar -e 'variant.hasAttribute("AF") && variant.getAlternateAlleles().size()==1 && variant.getAttributeAsDouble("AF",1.0)<0.05' 
(...)
1    10235    .    T    TA    100    PASS    AA=|||unknown(NO_COVERAGE);AC=6;AF=0.00119808;AFR_AF=0;AMR_AF=0.0014;AN=5008;DP=78015;EAS_AF=0;EUR_AF=0;NS=2504;SAS_AF=0.0051
1    10505    .    A    T    100    PASS    AA=.|||;AC=1;AF=0.000199681;AFR_AF=0.0008;AMR_AF=0;AN=5008;DP=9632;EAS_AF=0;EUR_AF=0;NS=2504;SAS_AF=0
1    10506    .    C    G    100    PASS    AA=.|||;AC=1;AF=0.000199681;AFR_AF=0.0008;AMR_AF=0;AN=5008;DP=9676;EAS_AF=0;EUR_AF=0;NS=2504;SAS_AF=0
1    10511    .    G    A    100    PASS    AA=.|||;AC=1;AF=0.000199681;AFR_AF=0;AMR_AF=0.0014;AN=5008;DP=9869;EAS_AF=0;EUR_AF=0;NS=2504;SAS_AF=0
1    10539    .    C    A    100    PASS    AA=.|||;AC=3;AF=0.000599042;AFR_AF=0;AMR_AF=0.0014;AN=5008;DP=9203;EAS_AF=0;EUR_AF=0.001;NS=2504;SAS_AF=0.001
1    10542    .    C    T    100    PASS    AA=.|||;AC=1;AF=0.000199681;AFR_AF=0;AMR_AF=0;AN=5008;DP=9007;EAS_AF=0.001;EUR_AF=0;NS=2504;SAS_AF=0
1    10579    .    C    A    100    PASS    AA=.|||;AC=1;AF=0.000199681;AFR_AF=0.0008;AMR_AF=0;AN=5008;DP=5502;EAS_AF=0;EUR_AF=0;NS=2504;SAS_AF=0
1    10642    .    G    A    100    PASS    AA=.|||;AC=21;AF=0.00419329;AFR_AF=0.0129;AMR_AF=0.0014;AN=5008;DP=1360;EAS_AF=0.003;EUR_AF=0;NS=2504;SAS_AF=0
1    11063    .    T    G    100    PASS    AA=.|||;AC=15;AF=0.00299521;AFR_AF=0.0106;AMR_AF=0.0014;AN=5008;DP=2834;EAS_AF=0;EUR_AF=0;NS=2504;SAS_AF=0
1    13011    .    T    G    100    PASS    AA=t|||;AC=3;AF=0.000599042;AFR_AF=0.0023;AMR_AF=0;AN=5008;DP=35822;EAS_AF=0;EUR_AF=0;NS=2504;SAS_AF=0

[Ram]

Why biallelic only?

[Pierre Lindenbaum]

Yes, I'm lazy. The real script would look like:

function accept(v)
{
if(!v.hasAttribute("AF")) return false;
if(v.getAlternateAlleles().size()==1) return v.getAttributeAsDouble("AF",1.0)<0.05
for(var af: v.getAttribute("AF")) if(af < 0.05) return true;
return false;
}

accept(variant)

[Ram]

Well, you're working for free so you can afford to be lazy :)

[aj123]

Thank you very much! Could you please explain this so that I can understand it and use this as a template to write my own script? also, is this in python?

[Pierre Lindenbaum]

As described on the doc page, my program is a java program. It takes as input a javascript program and inject in this script a variable called 'variant' which is an instance of https://samtools.github.io/htsjdk/javadoc/htsjdk/htsjdk/variant/variantcontext/VariantContext.html. The small problem here is that there can be more than one AF attribute if there is more than one ALT allele.

[aj123]

Thanks Pierre! Will take a look and post my script up sometime to check.

[Ram]

Why re-invent the wheel/write your own script? And why attempt to, when you have difficulties differentiating between JS and Python? Believe me, you're better off using existing tools.