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.

The variant which is seen single sample gvcf is missing after combing multiple samples

NeethuNeethu IndiaMember

@Sheila, I generally use the following command to combine gvcf of my samples. when new sample comes I keep adding them to already combined gvcfs
java -Xmx20g -jar /mnt/exome/Softwares/GenomeAnalysisTK.jar -T CombineGVCFs -R /mnt/exome/ReferenceFiles/human_g1k_v37.fasta -V Combined_200.gvcf -V samplea.gvcf -V sampleb.gvcf -V samplec.gvcf -o $outname.gvcf
But one of the variant which is seen in gvcf of a single sample is not seen in combined gvcf. Please help me to rectify the issue.

Thanks and Regards
Neethu

Tagged:

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @Neethu, did you check that your samples all have distinct sample names? Sometimes we see this happen when two samples had the same SM tag.

  • NeethuNeethu IndiaMember

    Hi @Geraldine_VdAuwera , All of our samples have unique id. The one variant which we are missed in combined GVCF resides on X-chromosome. Based on gvcf record of the sample its a hemizygous variant and many other samples are non variant for this location. But still the location is completely absent in combined gvcf . Then I recombined the gvcf using small set (30 samples) and combined gvcf shows the record for the particular location with zygosity './.' for all the 30 samples.

  • NeethuNeethu IndiaMember

    @Geraldine_VdAuwera , There are some gross mistakes happening for us in combined gvcf/or in joint genotyping. A good quality heterozygous variant '0/1' (validated using Sanger sequencing ) which is in single sample gvcf record, shows homozygous for wildtype '0/0' in genotyped gvcf. Could you please help us to rectify this issue as soon as possible.

  • NeethuNeethu IndiaMember

    @Geraldine_VdAuwera, to troubleshoot the issue I took only 100 gvcfs (includes the samples which has given the errors). When I combined these 100 and done the joint genotyping, I could find the missed variants in the in multi sample VCF file. Also the zygosity problem got resolved. Do you think, this has to do something with sample size and our system configuration. My sample size is 800. Please let me know.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Neethu
    Hi Neethu,

    Can you post the GVCF record and VCF record for the missing site?

    Thanks,
    Sheila

  • NeethuNeethu IndiaMember

    GVCF record:
    X 70341452 . G A, 1103.77 . DP=42;ExcessHet=3.0103;MLEAC=2,0;MLEAF=1.00,0.00;RAW_MQ=151200.00 GT:AD:DP:GQ:PL:SB 1/1:0,42,0:42:99:1132,126,0,1132,126,1132:0,0,18,24

    This record is absent in Genotyped VCF of 800 people. But when I have done the joint genotyping for 100 samples, I could see this record.

    X 70341452 . G A 1066.73 . AC=2;AF=0.010;AN=200;DP=5334;ExcessHet=0.0109;FS=0.000;InbreedingCoeff=1.0000;MLEAC=2;MLEAF=0.010;MQ=60.00;QD=25.40;SOR=1.005GT:AD:DP:GQ:PL 0/0:47,0:47:99:0,107,1221 0/0:52,0:52:99:0,120,1800 0/0:53,0:53:99:0,105,1575 0/0:37,0:37:99:0,105,1575 0/0:44,0:44:99:0,113,1356 0/0:35,0:35:99:0,99,1485 0/0:40,0:40:99:0,107,1351 0/0:42,0:42:99:0,111,1281 0/0:26,0:26:75:0,75,708 0/0:90,0:90:99:0,120,1800 0/0:41,0:41:99:0,108,1620 0/0:30,0:30:90:0,90,801 0/0:36,0:36:99:0,105,1387 0/0:38,0:38:99:0,99,1340 1/1:0,42:42:99:1132,126,0 0/0:38,0:38:99:0,108,985

  • NeethuNeethu IndiaMember

    @Sheila There is also an another case were, a particular variant is shown as '0/1' after combining and joint genotyping '0/1' has changed to '0/0'. Again in my 100 samples this has come correctly.
    4 187539753 . C T, 374.77 . BaseQRankSum=0.974;ClippingRankSum=0.000;DP=47;ExcessHet=3.0103;MLEAC=1,0;MLEAF=0.500,0.00;MQRankSum=0.000;RAW_MQ=169200.00;ReadPosRankSum=-0.697 GT:AD:DP:GQ:PL:SB 0/1:30,17,0:47:99:403,0,756,493,807,1300:15,15,8,9

    Generally I combine all the GVCF files and make a single combine GVCF and give to GenotypeGVCFs.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Neethu
    Hi,

    Which version of GATK are you using? Can you check if this happens with the very latest version? Also, can you try GenomicsDBImport?

    Thanks,
    Sheila

  • NeethuNeethu IndiaMember

    @Sheila
    We created gvcf file for each of the samples using 3.6 . Now I am using GATK4 and tried to import the gvcfs using GenomicsDBImport . After running for a while it got terminated with the following error

    terminate called after throwing an instance of 'RunConfigException'
    what(): RunConfigException : Syntax error in JSON file /tmp/root/loader_3026847754078477265.json

    I have used following command to run this

    /mnt/exome/Softwares/gatk-4.0.7.0/gatk --java-options "-Xmx4g -Xms4g" GenomicsDBImport -L /mnt/exome/Softwares/HG19/nexterarapidcapture_exome_target_region.bed --genomicsdb-workspace-path /mnt/exome/PROCESSING/GNOMIC_DB/DB \
    -V a.gvcf \
    -V b.gvcf \
    -V c.gvcf \
    -V d.gvcf \
    -V e.gvcf \
    -V f.gvcf \
    -V g.gvcf

    I have given 100 sample gvcf file . It would be really great if you can help us to rectify this issue.
    Thank you

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin
    edited August 2018

    @Neethu
    Hi,

    I see. Then I guess there may be an incompatibility between 3.6 and GATK4. Sorry for the hassle, but I need to make sure this is an issue in GATK4 before I ask for test data, as the team is not working on GATK3 anymore.

    Can you try to narrow down the issue to a few samples and re-run HaplotypeCaller in GATK4 and then GenotypeGVCFs? It would be great if you could narrow this down to a small subset of samples if we need to get test data as well.

    -Sheila

  • NeethuNeethu IndiaMember

    @sheila,
    The above mentioned error has been rectified by clearing the tmp folder. I have given only 6 sample gvcf. But now I am getting an error as follows. Looking forward to hearing from you. Thank you.

    terminate called after throwing an instance of 'LoadOperatorException'
    what(): LoadOperatorException : Could not open TileDB array for loading
    TileDB error message : [TileDB::StorageManager] Error: Cannot list TileDB directory; Directory buffer overflow

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Neethu
    Hi,

    Can you confirm this happens with the very latest version? I think a fix went in for this recently.

    Thanks,
    Sheila

  • NeethuNeethu IndiaMember

    @Sheila
    I have used the version : gatk-4.0.7.0

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Neethu
    Hi,

    Just to confirm, can you make sure this happens in 4.0.8.1?

    Thanks,
    Sheila

Sign In or Register to comment.