Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We encourage our fourm members to be more involved, jump in and help out your fellow researchers with their questions. GATK forum is a community forum and helping each other with using GATK tools and research is the cornerstone of our success as a genomics research community.We appreciate your help!

Test-drive the GATK tools and Best Practices pipelines on Terra


Check out this blog post to learn how you can get started with GATK and try out the pipelines in preconfigured workspaces (with a user-friendly interface!) without having to install anything.

LeftAlignAndTrimVariants --splitMultiallelics changes GT from known to unknown

tommycarstensentommycarstensen United KingdomMember ✭✭✭

I have a VCF file with this line (i.e. GT=0/1=G/T):

20  10120854    .   G   T,A 32175.56    .   AC=399,18;AF=0.111,5.006e-03;AN=3596;BaseQRankSum=1

.03;DP=6710;FS=2.485;GQ_MEAN=15.45;GQ_STDDEV=20.21;InbreedingCoeff=0.1235;MLEAC=416,17;MLEAF=0.116,4.727e-03;MQ=60.00;MQ0=0
;MQRankSum=0.358;NCC=189;QD=18.08;ReadPosRankSum=0.358 GT:AD:DP:GQ:PL 0/1:1,3,0:.:34:123,0,34,126,43,169

When I run it through version 3.2 of LeftAlignAndTrimVariants with the --splitMultiallelics flag, then the genotype information is lost; i.e. GT=./. and the output is:

20  10120854    .   G   T   32175.56    .   BaseQRankSum=1.03;DP=6710;FS=2.485;GQ_MEAN=15.45;GQ_STDDEV=20.21;InbreedingCoeff=0.1235;MLEAC=416,17;MLEAF=0.116,4.727e-03;MQ=60.00;MQ0=0;MQRankSum=0.358;NCC=189;QD=18.08;ReadPosRankSum=0.358 GT  ./.
20  10120854    .   G   A   32175.56    .   BaseQRankSum=1.03;DP=6710;FS=2.485;GQ_MEAN=15.45;GQ_STDDEV=20.21;InbreedingCoeff=0.1235;MLEAC=416,17;MLEAF=0.116,4.727e-03;MQ=60.00;MQ0=0;MQRankSum=0.358;NCC=189;QD=18.08;ReadPosRankSum=0.358 GT  ./.

I have attached the input VCF. I also got rid of the PL information, but still got the unexpected output.

Maybe I should just write some code for normalizing variants myself in order to save time :) Thank you very much as always!

Issue · Github
by Geraldine_VdAuwera

Issue Number
834
State
closed
Last Updated
Assignee
Array
Closed By
ronlevine

Best Answer

Answers

  • ravichavravichav United StatesMember

    Have the features been updated? I tried GATK3.3 and while the VCF is left aligned, the genotypes become no calls.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @ravichav
    Hi,

    The fix has not been made yet. One of our developers is working on it now. I will let you know when the fix is in.

    -Sheila

  • jsreddy82jsreddy82 Mayo Clinic, Jacksonville, FL.Member ✭✭

    Hi @Sheila,

    Has the tool been updated to retain appropriate genotype calls (and all other related annotations/fields specific to the alternate allele across samples) after splitting multi-allelic sites?

    I wish to split multi-allelic sites (containing SNPs and/or INDELs) in a multi-sample VCF (generated after joint genotyping) into biallelic form and then left normalize (trim) alleles. Subsequently, I wish to perform VQSR and obtain a FILTER score for all variants. Is this currently possible?

    Thanks,
    Joseph.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @jsreddy82
    Hi Joseph,

    Yes, the fix went in last year.

    For VQSR, yes you will obtain a value in the FILTER column for all your variants. It will say PASS for passing variants and the tranche the failing variants fall into.

    -Sheila

  • jsreddy82jsreddy82 Mayo Clinic, Jacksonville, FL.Member ✭✭
    edited July 2016

    @Sheila Thank you for your response. In relation to VQSR, earlier we used to split alleles and their corresponding genotypes after VQSR using in-house scripts, thereby "copying" VQSR FILTER score across all split variants. If LeftAlignAndTrimVariants splits genotypes (and all related information) across alleles, [A->G,C 0/2 to a A->G ./. and A->C 0/1] it will now help us assign independent VQSR scores to each. Which version of GATK should I be using? We are currently using v3.3

    Thanks,
    Joseph.

    Post edited by jsreddy82 on
  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin
    edited July 2016

    @ravichav
    Hi,

    I did not let you know the fix went in like I told you I would! Sorry about that. The fix has been made now. If you upgrade to the latest version, you will have the fix.

    @jsreddy82
    Hi Joseph,

    It looks like the fix went in last December. However, it is best to use the latest version of GATK (3.6) so you will have all the latest features. We recommend using the same version of GATK for your entire analysis. https://www.broadinstitute.org/gatk/documentation/article?id=3536 So, if your final VCF was made with 3.3, you should really re-do your entire analysis with 3.6 for best results.

    -Sheila

  • jsreddy82jsreddy82 Mayo Clinic, Jacksonville, FL.Member ✭✭

    Thank you @Sheila for your recommendations. Will keep in mind.

    Thanks,
    Joseph.

  • ravichavravichav United StatesMember

    I tried using this walker with GATK3.6. Some variants give out INFO and some splits doesnt have an INFO. IS there a reasoning behind that?

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @ravichav
    Hi,

    It looks like the sites where the INFO is missing are sites where the genotypes are ./.
    For example, I have this site in my VCF:
    1 1 . G A,C 299.59 . AC=1,1;AF=0.500,0.500;AN=2;BaseQRankSum=-2.335e+00;DP=46;FS=8.450;GQ_MEAN=272.00;GQ_STDDEV=324.74;MQ=59.30;MQ0=0;MQRankSum=0.00;NCC=0;QD=27.52;ReadPosRankSum=-1.090e+00;SOR=0.108 GT:AD:DP:GQ:PL ./. 1/2:0,28,18:46:99:1334,720,636,679,0,1026

    In my VCF from LeftAlignAndTrimVariants, that site looks like:
    1 1 . G A 299.59 . . GT:DP ./. ./.:46
    1 1 . G C 299.59 . . GT:DP ./. ./.:46

    Notice the ./. genotypes are there because one sample was already ./. in the original VCF. But, the second sample was 1/2 genotype. Because you split the multiallelic site, the 1 and 2 alleles in the original genotype are no longer present and cannot be called properly.

    -Sheila

  • sheilaztsheilazt ValenciaMember

    Hi Sheila,
    I'm currently running the latest version of GATK (v3.7-0-gcfedb67) and the LeftAlignAndTrimVariants tool on a GVCF file I obtained with the GenotypeGVCF tool this way:
    gatk -T LeftAlignAndTrimVariants -R GRCh38.fasta --variant all_files.g.vcf -o output_no_multi.vcf --splitMultiallelics

    Original file:
    chr1 22511093 . CTT CT,C 1209.37 . AC=4,2;AF=0.333,0.167;AN=12;BaseQRankSum=1.76;ClippingRankSum=0.00;DP=195;ExcessHet=7.2016;FS=21.329;MLEAC=4,2;MLEAF=0.333,0.167;MQ=57.80;MQRankSum=0.00;QD=11.63;ReadPosRankSum=-3.520e-01;SOR=0.076 GT:AD:DP:GQ:PL ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. 1/2:2,28,11:41:99:816,217,145,496,0,572 ./. ./. ./. ./. ./. ./.:11,0,0:11:.:0,0,0,0,0,0 0/2:9,0,5:14:99:131,159,416,0,257,241 0/1:14,7,0:21:38:38,0,184,77,202,279 0/1:4,5,0:9:36:39,0,36,49,49,98 0/0:42,0,0:42:0:0,0,545,0,545,545 ./. ./. ./. ./. ./. ./. ./. ./. ./. 0/1:6,13,0:19:56:218,0,56,235,94,329 ./. ./. ./. ./. ./. ./. ./. ./. ./. ./.

    After applying the LeftAlignAndTrimVariants tool:
    chr1 22511093 . CT C 1209.37 . . GT:DP ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./.:41 ./. ./. ./. ./. ./. ./.:11 ./.:14 ./.:21 ./.:9 ./.:42 ./. ./. ./. ./. ./. ./. ./. ./. ./. ./.:19 ./. ./. ./. ./. ./. ./. ./. ./. ./. ./.
    chr1 22511093 . CTT C 1209.37 . . GT:DP ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./.:41 ./. ./. ./. ./. ./. ./.:11 ./.:14 ./.:21 ./.:9 ./.:42 ./. ./. ./. ./. ./. ./. ./. ./. ./. ./.:19 ./. ./. ./. ./. ./. ./. ./. ./. ./. ./.

    The old multi-allelic variants loose the correct genotype, the AD, the PL and other flags as first reported in this thread. Is this the right behavior for the tool?

    Any help would be much appreciated.

    Thanks very much in advance.

    Regards,

    Sheila

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin
    That is indeed the expected behavior right now, but we have a ticket to change it because it's not actually necessay in most records.
  • chaozhongyinxiangchaozhongyinxiang Philadelphia, PAMember

    @Geraldine_VdAuwera said:
    That is indeed the expected behavior right now, but we have a ticket to change it because it's not actually necessay in most records.

    Hi Geraldine, do we have an ETA on when this feature will be ready? Thank you.

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Hi @chaozhongyinxiang,

    You can track the issue with our bug tracker. Our engineer is currently actively working on this and I'll let you know when I find out more.

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Alright, it appears this fix is just waiting for review. I'll let you know when you can download the latest nightly build to test it.

  • mattimatti FinlandMember

    Dear Shlee, what is the status of the fix. Tried to check that out already from your bug tracker, but without a great success.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin
    The fix was merged two weeks ago so you should be able to use the nightly build successfully.
  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    I must have missed the notice about the fix going through. Apologies and thanks for the reminder.

  • mattimatti FinlandMember

    Hi,
    downloaded the nightly build. To best of my understanding, I still suffer from the same error.

    java -Xmx5g -Djava.io.tmpdir=./tmp -jar GenomeAnalysisTK.jar -T LeftAlignAndTrimVariants -R homo_sapiens_82.fa --variant vcf/31-03-15-EN.final.vcf -o test.1.norm.vcf --splitMultiallelics

    INFO 11:49:31,612 HelpFormatter - The Genome Analysis Toolkit (GATK) vnightly-2017-03-15-ga4f7203, Compiled 2017/03/15 00:01:14
    INFO 11:49:31,612 HelpFormatter - Copyright (c) 2010-2016 The Broad Institute
    INFO 11:49:31,613 HelpFormatter - For support and documentation go to https://software.broadinstitute.org/gatk
    INFO 11:49:31,613 HelpFormatter - [Wed Mar 15 11:49:31 EET 2017] Executing on Linux 3.10.0-229.7.2.el7.x86_64 amd64
    INFO 11:49:31,614 HelpFormatter - OpenJDK 64-Bit Server VM 1.8.0_65-b17
    INFO 11:49:31,620 HelpFormatter - Program Args: -T LeftAlignAndTrimVariants -R homo_sapiens_82.fa --variant vcf/31-03-15-EN.final.vcf -o test.1.norm.vcf --splitMultiallelics --reference_window_stop 300
    INFO 11:49:31,632 HelpFormatter - Executing as [email protected] on Linux 3.10.0-229.7.2.el7.x86_64 amd64; OpenJDK 64-Bit Server VM 1.8.0_65-b17.
    INFO 11:49:31,632 HelpFormatter - Date/Time: 2017/03/15 11:49:31

    My input VCF has the following line:

    16 78764083 . G A,C 413.77 . AC=1,1;AF=0.500,0.500;AN=2;DP=11;FS=0.000;MLEAC=1,1;MLEAF=0.500,0.500;MQ=60.00;QD=35.75;SOR=0.859 GT:AD:DP:GQ:PL 1/2:0,5,6:11:99:442,235,237,210,0,192

    that becomes:

    16 78764083 . G A 413.77 . . GT:DP ./.:11
    16 78764083 . G C 413.77 . . GT:DP ./.:11

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin
    I believe the fix was to retain genotypes like 0/1 or 0/2 since no information changes for them. The example you show (1/2) is more difficult since one of the genotype alleles is no longer present in each resulting record -- we would have to emit a half-genotype for those. We haven't yet decided how that should be handled. In a case like that ideally you wouldn't split the call.
  • dmarkiedmarkie Dunedin, New ZealandMember

    This still seems to be an issue with GATK4.1.0.0 when using LeftAlignAndTrimVariants to splitMultiallelics. It works well with many multiallelic variants (including trisomic regions it seems!), but when there are one or more calls that are heterozygous for two ALT alleles, then all genotypes are converted to uncalled (./.) and the INFO field is deleted. With increasing cohort size this situation is becoming more frequent. The normalise function from BCFtools handles this in a mostly acceptable way - it will convert a 1/2 genotype into a 1/0 and a 0/1 genotype over two lines - so although the reference allele is somewhat misrepresented, there is not the dramatic (and sometimes unnoticed) loss of information that can occur. Perhaps another option would be to represent the genotypes as 1/. and ./1 to make it explicit that there is another allele that is unmentioned (instead of mislabelling it as reference).
    BCFtools would be a solution for our needs except it currently doen't function beyond a ploidy of 2, and we would like to keep the occasional sex chromosome trisomies accurately represented in our cohorts.
    Is there any resolution on the way for this problem - even if it is only to fail with an appropriate warning that data is about to be lost?
    Many thanks

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin
    edited March 7

    HI @dmarkie

    Would you please post an example of this issue. This will help our dev team understand better why this might be happening.

  • dmarkiedmarkie Dunedin, New ZealandMember

    Hi, I have chopped down the vcf to just 4 illustrative examples in Test.vcf.gz and then run LeftAlignAndTrimVariants --split-multi-allelics --dont-trim-alleles --keep-original-ac using GATK4.1.0.0 to produce Test_Split.vcf.gz (see attached files). The biallelic variant is unaffected, one of the multiallelic variants is split as expected but the first two examples are split but the INFO fields are empty and all the genotypes are ./. This seems to be associated with having a heterozgous ALT call in the original vcf eg 1/2, but I'm not completely sure that is always the case.

Sign In or Register to comment.