Filtering variants by depth of coverage using SelectVariants

I am trying to manually filter the results from HaplotypeCaller because my organism is anything but model. I sequenced 60 samples and had an average coverage of ~60X per sample across the genome.

As shown in an excerpt from my HaplotypeCaller output below, "DP=3201". Is this the total depth? i.e. across 60 samples the average is 53.35? This number makes sense. If this is the case, is there any way to exclude SNPs in which ANY INDIVIDUAL sample has a depth of less than 20? As opposed to an average depth of less than 20.

I've been using SelectVariants as such:

java -jar /home/berrya/bin/GenomeAnalysisTK.jar -T SelectVariants -R /path/to/reference.fasta -V HaplotypeCaller_output.vcf -select "DP >=20" -select "QUAL > 30.0" -select "QD < 2.0 || FS > 60.0 || SOR > 4.0 || MQ < 40.0 || MQRankSum < -12.5 ||ReadPosRankSum < -8.0" -o output.vcf

When I raise the DP variable to "DP >=1000", no additional SNPs are cut, which is how I know it's not doing what I'm expecting it to be doing.

I appreciate any guidance.

Alex

KB223232.1 200 . T A 112397.32 .
AC=120;AF=1.00;AN=120;BaseQRankSum=1.50;ClippingRankSum=0.945;DP=3201;ExcessHet=3.0103;FS=0.000;InbreedingCoeff=-0.0000;MLEAC=120;MLEAF=1.00;MQ=36.90;MQRankSum=0.394;QD=34.24;ReadPosRankSum=-1.339e+00;SOR=0.223 GT:AD:DP:GQ:PL 1/1:0,53:53:99:1865,159,0 1/1:0,47:47:99:1744,141,0 1/1:0,179:179:99:6518,536,0 1/1:0,40:40:99:1432,119,0 1/1:0,15:15:45:497,45,0 1/1:0,20:20:59:679,59,0 1/1:0,39:39:99:1317,116,0 1/1:0,63:63:99:2167,187,0 1/1:0,34:34:99:1274,102,0 1/1:0,30:30:90:1105,90,0 1/1:1,43:44:99:1507,121,0 1/1:0,140:140:99:5024,418,0 1/1:0,17:17:50:558,50,0 1/1:0,39:39:99:1315,116,0 1/1:0,24:24:72:872,72,0 1/1:0,17:17:50:534,50,0 1/1:0,34:34:99:1158,101,0 1/1:0,33:33:97:1053,97,0 1/1:0,36:36:99:1133,106,0 1/1:0,71:71:99:2572,212,0 1/1:0,35:35:99:1220,104,0 1/1:0,42:42:99:1512,126,0 1/1:0,39:39:99:1402,116,0 1/1:0,44:44:99:1602,131,0 1/1:0,42:42:99:1419,124,0 1/1:0,126:126:99:4695,378,0 1/1:0,48:48:99:1708,143,0 1/1:0,41:41:99:1487,123,0 1/1:0,132:132:99:4789,396,0 1/1:0,39:39:99:1340,116,0 1/1:0,131:131:99:4785,393,0 1/1:0,34:34:99:1161,101,0 1/1:0,35:35:99:1304,105,0 1/1:0,27:27:81:1001,81,0 1/1:0,32:32:96:1190,96,0 1/1:0,51:51:99:1840,153,0 1/1:0,38:38:99:1347,114,0 1/1:0,64:64:99:2151,190,0 1/1:0,47:47:99:1684,140,0 1/1:0,66:66:99:2346,197,0 1/1:0,64:64:99:2210,190,0 1/1:0,51:51:99:1846,153,0 1/1:0,48:48:99:1651,142,0 1/1:0,20:20:60:744,60,0 1/1:0,38:38:99:1369,114,0 1/1:0,20:20:59:691,59,0 1/1:0,65:65:99:2307,193,0 1/1:0,181:181:99:6516,542,0 1/1:0,60:60:99:2022,177,0 1/1:0,44:44:99:1586,131,0 1/1:0,133:133:99:4473,394,0 1/1:0,44:44:99:1527,131,0 1/1:0,42:42:99:1369,124,0 1/1:0,51:51:99:1794,152,0 1/1:0,38:38:99:1321,113,0 1/1:0,85:85:99:2941,255,0 1/1:0,42:42:99:1299,122,0 1/1:0,28:28:83:957,83,0 1/1:0,26:26:76:856,76,0

Best Answer

Answers

  • I guess my question, put more simply is: Can you filter by FORMAT (variant-level) instead of INFO (site-level)?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Variant-level is the same thing as site-level; I think you mean sample-level or genotype-level.

    And yes, you can filter at the FORMAT level. Have a look at the VariantFiltration tool doc, or the doc on filtering using JEXL.

  • Hi Geraldine, thanks for the quick reply. I have looked at the docs for VariantFiltration, SelectVariants, and JEXL. (I only chose SV over VF because it was easier to broadly test filtering strictness by looking at file sizes.) I'm comfortable filtering by INFO (Qual, DP, etc. as used above) but I have not had success using FORMAT commands such as:

    java -jar GenomeAnalysisTK.jar -T SelectVariants -R reference.fasta -V input.vcf -select "GQ > 50" -env -o output.vcf

    or

    java -jar GenomeAnalysisTK.jar -T SelectVariants -R reference.fasta -V input.vcf -select "AD >= 20" -env -o output.vcf

    My assumption was that, for the latter example, SelectVariants would use any site for which all 60 samples have a depth (AD) above 20. However, that commands yields no sites whatsoever.

    Is there a document on how to filter by GT, AD, DP, GQ, or PL? It's likely that I'm making a silly mistake.

    Thanks again,

    Alex

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @aberry814
    Hi Alex,

    To do what you want, I think you need to use VariantFiltration. You can use -G_filter and -G_filterName.

    Also, have a look at this article under "Using JEXL to evaluate arrays" for an example using AD.

    -Sheila

Sign In or Register to comment.