CombineVariants prioritize option

meharmehar Member ✭✭
edited September 2014 in Ask the GATK team

Hi,

I have used CombineVariants to combine variants from GATK and samtools as shown below:

java -jar GenomeAnalysisTK.jar -T CombineVariants -R ref.fa --variant:GatkSNP GATKsnp.vcf --variant:GatkINDEL GATKind.vcf --variant:SamSNP Samsnp.vcf --variant:SamINDEL Samind.vcf -o allvar.vcf -genotypeMergeOptions PRIORITIZE -priority GatkSNP,GatkINDEL,SamSNP,SamINDEL --filteredrecordsmergetype KEEP_UNCONDITIONAL

This merges all the variants. However, with the above command, i do get the variants present in both GATK and samtools emitted from samtools.

I would like to get all the variants such that:

  • variants present in both GATK and samtools emitted from GATK vcf files
  • variants in only GATK
  • variants in only samtools

could someone suggest any ideas or of there is something to be fixed in the command.

Thanks

Post edited by Geraldine_VdAuwera on

Best Answer

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @mehar‌

    Hi,

    Can you please clarify which variants exactly you want?

    Another option is after running CombineVariants, you can use SelectVariants and JEXL expressions to retrieve exactly what you want.

    Please read about SelectVariants here: http://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.html
    and how to use JEXL expressions here: http://www.broadinstitute.org/gatk/guide/article?id=1255

    -Sheila

  • meharmehar Member ✭✭

    Thanks for your suggestion. I thought i clarified in the initial post what i need. But anyway here is what i need: 1. variants present in both GATK and samtools emitted from GATK vcf files 2.variants in only GATK 3.variants in only samtools.

    The command i have used is doing this, but the variants which are present in both GATK and samtools are emitted from samtools vcf whereas i want them to be emitted from GATK vcf files.

    And SelectVariants is something to make a subset of the variants based on different criteria and that is not what i need here. I hope i clarified to some extent atleast!!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @mehar It's still not clear what you want to end up with. Do you want one file with the various combinations tagged with an identifying annotation, or do you want separate files containing the different combinations?

  • meharmehar Member ✭✭

    I want all the variants into a single file. i.e. variants that are common in all the four vcfs and also the unique variants in each vcf file.
    For the variants that are common in all the vcfs, i want them to be emitted from gatk vcf files.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    What do you mean when you say "i want them to be emitted from gatk vcf files"?

  • meharmehar Member ✭✭

    May be i will put it in a different way:

    For instance,

    **[[email protected] VCF]$ grep -w "37842119" GATKsnp.vcf
    chr17   37842119    .   G   A   2838.77 .   AC=2;AF=1.00;AN=2;BaseQRankSum=1.615;DP=94;Dels=0.00;FS=0.000
    HaplotypeScore=0.9332;MLEAC=2;MLEAF=1.00;MQ=57.32;MQ0=3;MQRankSum=1.693;QD=30.20;ReadPosRankSum=1.304         
    GT:AD:DP:GQ:PL    1/1:5,89:94:99:2867,262,0
    
    [[email protected] VCF]$ grep -w "37842119" SAMsnp.vcf
    chr17   37842119    .   G   A   124 .   DP=89;VDB=3.271741e-02;AF1=1;AC1=2;DP4=0,0,20,0;MQ=59;FQ=-87    GT:PL:GQ     
    1/1:157,60,0:99
    

    When CombineVariants is used to combine variants from GATK and Samtools, the above variant is reported as shown below:

    chr17   37842119    .   G   A   2838.77 .   AC=2;AC1=2;AF=1.00;AF1=1;AN=2;BaseQRankSum=1.615;DP=94;DP4=0,0,20,0;    
    Dels=0.00;FQ=-87;FS=0.000;HaplotypeScore=0.9332;MLEAC=2;MLEAF=1.00;MQ0=3;MQRankSum=1.693;QD=30.20;
    ReadPosRankSum=1.304;VDB=3.271741e-02;set=GatkSNP-SamSNP    GT:AD:DP:GQ:PL  1/1:5,89:94:99:2867,262,0
    

    Incase of this single alternate allele SNP, the resultant variant in CombineVariants has the INFO and FORMAT fields that are present in GATK vcf file, which means this SNP is emitted from GATK vcf file. I have no issues in single allelic SNPs.**

    _But in case of the below biallelic SNP,

    [[email protected] VCF]$ grep -w "20608169" GATKsnp.vcf
    chr27   20608169    .   A   G   2571.77 .   AC=2;AF=1.00;AN=2;DP=84;Dels=0.00;FS=0.000;HaplotypeScore=5.2101;MLEAC=2;
    MLEAF=1.00;MQ=59.70;MQ0=0;QD=30.62  GT:AD:DP:GQ:PL  1/1:0,82:84:99:2600,238,0
    
    [[email protected] VCF]$ grep -w "20608169" SAMsnp.vcf
    chr27   20608169    .   A   G,T 135 .   DP=83;VDB=3.012910e-01;AF1=1;AC1=2;DP4=0,0,0,53;MQ=60;FQ=-184   GT:PL:GQ    
    1/1:168,157,0,158,129,155:99
    

    When CombineVariants is used to combine variants from GATK and Samtools, the above variant is reported as shown below:

    chr27   20608169    .   A   G,T 2571.77 .   AC=2,0;AC1=2;AF=1.00,0.00;AF1=1;AN=2;DP=84;DP4=0,0,0,53;Dels=0.00;FQ=-184;
    FS=0.000;HaplotypeScore=5.2101;MLEAC=2;MLEAF=1.00;MQ0=0;QD=30.62;VDB=3.012910e-01;set=GatkSNP-SamSNP    
    GT:DP:GQ    1/1:84:99
    

    In this case the INFO and FORMAT fields completely differ in the content from both GATK and Samtools VCF files.

    As per my understanding, The INFO field is modified such that it carries all the information from INFO fields in both GATK and Samtools. Inother words, it adds the extra fields in INFO column of samtools VCF to GATK INFO column. This is fine so far.

    But, the FORMAT column differs completely. Is it not possible to retain the GATK style of FORMAT column?_

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Oh I see now. Thanks for the example, that is much clearer. What's happening here is that for annotations like AD and PL, the tool is not able to produce meaningful merged values when some alleles are absent in one of the callsets, so it simply drops those annotations. There is no workaround for this aside from regenotyping the site with the given alleles.

  • meharmehar Member ✭✭
    edited September 2014

    Thanks for your time.

    I have another issue. After combining the variants, i used VariantsToTable to convert into a tabular format. But, if find difference in number of variants between the CombineVariants VCF files and VariantsToTable file.

     Command used for VariantsToTable:
      java -jar GenomeAnalysisTK.jar -T VariantsToTable -R %s -V %s --allowMissingData -F CHROM -F POS 
     -F ID -F REF -F ALT -F QUAL -F MQ -F DP4 -F AD -GF GT -GF AD -F FILTER -F VQSLOD -F set -F resource.EFF  -o %s
    
     $grep -v "#"SNP_INDELS_Annotated.vcf | wc -l
     35927
     $ wc -l DDT076SNP_INDELS_Annotated.table 
     33533 DDT076SNP_INDELS_Annotated.table
    

    Could you comment on what has been going on here!!

    Post edited by mehar on
  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @mehar‌

    Hi,

    Do you know which records exactly are missing? Are they mixed records of SNPs and Indels? It will help us greatly if you can narrow down what types of records are missing.

    Thanks,
    Sheila

  • meharmehar Member ✭✭

    Hi,

    The missing records has only SNPs. I have attached the file of missing SNPs if that can help.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @mehar‌

    Hi,

    I see that all of the records that are not included have a LOW_VQSLOD. Are there any of the included records that have a LOW_VQSLOD?

    Maybe if you leave out -F VQSLOD, the records that have LOW_VQSLOD will be included in the table.

    -Sheila

  • meharmehar Member ✭✭

    Thanks. Included records do not have LOW_VQSLOD. But, if i leave -F VQSLOD, VQSR annotation will be missing in the final variant table and it makes no sense after performing VQSR. I wonder why these variants are filtered when VariantsToTable is used, which should only convert the vcf to tabular format.

  • meharmehar Member ✭✭
    edited September 2014

    Thanks sheila!! Now i have another minor concern.

    After using CombineVariants tools to combine variants from GATK and Samtools, MQ annotation is missing for some of the variants specifically for variants which are present in both GATK and Samtools. Below is an example:

    $grep "16203159" VCF/GatkSNP.vcf
    chr8    16203159    .   G   T   4609.77 .   AC=1;AF=0.500;AN=2;BaseQRankSum=-7.105;ClippingRankSum=0.652;DP=333;     
    FS=1.794;MLEAC=1;MLEAF=0.500;MQ=57.94;MQ0=0;MQRankSum=1.205;QD=13.84;ReadPosRankSum=0.157   
    GT:AD:DP:GQ:PL  0/1:163,170:333:99:4638,0,4651
    
    $ grep "16203159" VCF/SamSNP.vcf
    chr8    16203159    .   G   T   225 .   DP=346;VDB=4.969420e-01;RPB=-1.470555e+00;AF1=0.5;AC1=1;DP4=82,84,81,92;  
    MQ=58;FQ=225;PV4=0.66,6.2e-09,1,1   GT:PL:GQ    0/1:255,0,255:99
    

    MQ annotation is present in both the VCF files. But after combining MQ is missing as shown below:

    $ grep "16203159" TMP_files/SNP_INDEL.vcf
    chr8    16203159    .   G   T   4609.77 .   AC=1;AC1=1;AF=0.500;AF1=0.5;AN=2;BaseQRankSum=-7.105;  
    ClippingRankSum=0.652;DP=679;DP4=82,84,81,92;FQ=225;FS=1.794;MLEAC=1;MLEAF=0.500;MQ0=0;MQRankSum=1.205;  
    PV4=0.66,6.2e-09,1,1;QD=13.84;RPB=-1.470555e+00;ReadPosRankSum=0.157;VDB=4.969420e-01;set=GatkSNP-SamSNP    
    GT:AD:DP:GQ:PL    0/1:163,170:333:99:4638,0,4651
    

    Variants which are present in only of the VCF files do have MQ annotation after using CombineVariants. This happens only for variants which are present in both VCF files. Could you help to get MQ annotations in this case?

    Here is the command used:

     java -jar GenomeAnalysisTK.jar -T CombineVariants -R ref.fa --variant:GatkSNP GatkSNP.vcf --variant:GatkINDEl GatkINDEL.vcf  
    --variant:SamSNP SamSNP.vcf --variant:SamINDEL SamINDEL.vcf -o SNP_INDEL.vcf -genotypeMergeOptions PRIORITIZE -priority    
    GatkSNP,GatkINDEL,SamSNP,SamINDEL --filteredrecordsmergetype KEEP_UNCONDITIONAL
    
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @mehar,

    I'm not sure but I think it's because CombineVariants does not have a rule to combine the MQ value. I will check if this is correct.

  • meharmehar Member ✭✭

    Hi Geraldine,

    I have noticed few other annotations are also missing from CombineVariants output.

    As mentioned earlier, MQ is missing for variants (SNPs+ INDELs) present in both GATK and samtools. In addition,
    1. QD and FS annotations are missing for SNPs+ INDELs identified by samtools, this could be bcoz samtools doesn't have these
    annotations
    2. MQRankSum and ReadPosRankSum are missing for all the variants

    Lack of these annotation values doesn't allow to use hard filtering. Could you help to get this information. Thanks.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @mehar, you shouldn't be trying to filter variants together if they come from different programs' output. Even if the annotations are present in both, they are usually not calculated exactly the same way by the different programs. If you have variants from different programs, you should filter them before combining them.

  • meharmehar Member ✭✭

    Thanks. I agree with you that annotations will be different in different programs in terms of calculation. However, annoations like MQ is missing in the CombineVariants output. And It is obvious that MQRankSum and ReadPosRankSum are absent for samtools variants however i would like to bring to your notice that MQRankSum and ReadPosRankSum are also missing for GATK variants in CombineVariants output.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @mehar Sorry for the late response. CombineVariants doesn't know how to process those annotations so they are dropped.

Sign In or Register to comment.