To celebrate the release of GATK 4.0, we are giving away free credits for running the GATK4 Best Practices pipelines in FireCloud, our secure online analysis portal. It’s first come first serve, so sign up now to claim your free credits worth $250. Sponsored by Google Cloud. Learn more at https://software.broadinstitute.org/firecloud/documentation/freecredits

A 50 bp insertion (INDEL) varaint calling

Hi All,

I am using both UnifiedGenotyper and HaplotypeCaller to genotype given alleles from amplicon based sequencing data.
It seems that both caller did not performe well on variant rs1799752, a 50 bp insertion.
Sample was verified by agarose gel method on PCR products, which was considered as gold standard for this variant.

I looked into the sequences from bam files, the insertion sequence is indeed in the bam file, but can not be called by both caller.

Any suggestion?

Command line I used

primer sequences were hard clipped after bwa mapping, calling was performed under given alleles mode

  1. UG:
    java -jar GenomeAnalysisTK.jar -T UnifiedGenotyper -R /data/iGenomes/Homo_sapiens/NCBI/build37.2/Sequence/WholeGenomeFasta/genome.fa -stand_call_conf 10 -D /data/NGS/analysis/snps.vcf --genotyping_mode GENOTYPE_GIVEN_ALLELES --alleles /data/NGS/analysis/snps.vcf -I /data/NGS/analysis/GWI127_combined.primerclipped.bam -L 17:61565590-61566290 -o GWI127.rs1799752.ug.vcf
  2. HC:
    java -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R /data/iGenomes/Homo_sapiens/NCBI/build37.2/Sequence/WholeGenomeFasta/genome.fa -stand_call_conf 10 -D /data/NGS/analysis/snps.vcf --genotyping_mode GENOTYPE_GIVEN_ALLELES --alleles /data/NGS/analysis/snps.vcf -I /data/NGS/analysis/GWI127_combined.primerclipped.bam -L 17:61565590-61566290 -o GWI127.rs1799752.hc.vcf
  3. HC without given alleles:
    java -jar GenomeAnalysisTK.jar -T UnifiedGenotyper -R /data/iGenomes/Homo_sapiens/NCBI/build37.2/Sequence/WholeGenomeFasta/genome.fa -stand_call_conf 10 -D /data/NGS/analysis/snps.vcf -I /data/NGS/analysis/GWI127_combined.primerclipped.bam -L 17:61565590-61566290 -o GWI127.rs1799752.hc.v2.vcf

I also added -activeRegionMaxSize 1000/3000 as suggested in the forum, but it did not work as well.

$samtools view /data/NGS/analysis/bamclipper/GWI127_combined.primerclipped.bam | grep ATACAGTCACTTTTTTTTTTTTTTTGAGACGGAGTCTCGCTCTGTCGCCC
E00491:102:H23N2CCXY:7:1218:18791:48107 97 1 100426255 0 56S94M 17 61565891 0 GAGAGCCACTCCCATCCTTTCTCCCATTTCTCTAGACCTGCTGCCTATACAGTCACTTTTTTTTTTTTTTTGAGACGGAGTCTCGCTCTGTCGCCCCGGCTGGAGTGCTGTGGCGGGATCTCGGCTCCCTGCAAGCTCCGCCTCCCGGGT AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJ<JJJJJJJJJJJJJJJJJJFJFFJJJJJJJJJJJJJJJ<<-<FAAFFJJJF<JA7<FJ<<FFF-77<J7
-77<A<-7F--F7<77FAAF-AA-7---<7-7-AFFAFFAF<AA7<7) NM:i:2 MD:Z:52A18A22 AS:i:84 XS:i:84 RG:Z:GWI127_combined SA:Z:17,61565845,+,60M90S,0,0;
E00491:102:H23N2CCXY:7:2116:24150:44697 99 17 61565813 60 19S92M39S = 61565891 181 AGGTGTCTGCAGCATGTGGCCCCAGGCCGGGGACTCTGTAAGCCACTGCTGGAGAGCCACTCCCATCCTTTCTCCCATTTCTCTAGACCTGCTGCCTATACAGTCACTTTTTTTTTTTTTTTGAGACGGAGTCTCGCTCTGTCGCCCAGG AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ
JJJJJFJJJFJJJJJJJJJJJJJJJJJJ-7--FJJFFF)7-7AJJJ<J-7JJJ<FF AS:i:111 XS:i:24 RG:Z:GWI127_combined SA:Z:2,149045941,-,47M103S,0,0;
E00491:102:H23N2CCXY:7:2116:19136:33902 99 17 61565813 60 19S92M39S = 61565891 181 AGGTGTCTGCAGCATGTGGCCCCAGGCCGGGGACTCTGTAAGCCACTGCTGGAGAGCCACTCCCATCCTTTCTCCCATTTCTCTAGACCTGCTGCCTATACAGTCACTTTTTTTTTTTTTTTGAGACGGAGTCTCGCTCTGTCGCCCAGG AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJFJJFJJJJJJJJJJJJJJJ
JJJJJJJJJJJJJJFJJJJJJJJJJJJJ-77<FFJ7FJJFFAJJ<)7F<AJJJ-<- AS:i:111 XS:i:24 RG:Z:GWI127_combined SA:Z:6,17717274,+,103S47M,0,0;
E00491:102:H23N2CCXY:7:2124:29011:16727 99 17 61565813 60 19S92M39S = 61565891 181 AGGTGTCTGCAGCATGTGGCCCCAGGCCGGGGACTCTGTAAGCCACTGCTGGAGAGCCACTCCCATCCTTTCTCCCATTTCTCTAGACCTGCTGCCTATACAGTCACTTTTTTTTTTTTTTTGAGACGGAGTCTCGCTCTGTCGCCCAGT AAFFFAJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJJFJF<JJJJJJ<JJJJJJFJJFFJAJJJJJJJF7-FJJJJJJFJJJJJJJJJJJFJJJ
JJFA7FJAJFJJJJAJJJJAF<<A<-<FA-A-777<F-F7-)<)7A7JF)7)<))) AS:i:111 XS:i:24 RG:Z:GWI127_combined SA:Z:6,17717274,+,103S47M,0,1;
E00491:102:H23N2CCXY:7:2122:19928:51869 99 17 61565813 60 19S92M39S = 61565891 181 AGGTGTCTGCAGCATGTGGCCCCAGGCCGGGGACTCTGTAAGCCACTGCTGGAGAGCCACTCCCATCCTTTCTCCCATTTTTCTAGACCTGCTGCCTATACAGTCACTTTTTTTTTTTTTTTGAGACGGAGTCTCGCTCTGTCGCCCAGG AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJAJJJJJJJJJJJJJJJJFJFJJFJJJJJJJJJJJJJJJJFJJJJJFJ
JJJFJFJFJJFJJJJJFJJJJJJJJJJJ--7-AAFJJJ<FAJJJ7FAFJAF<<)7F AS:i:106 XS:i:24 RG:Z:GWI127_combined SA:Z:6,17717274,+,103S47M,0,0;
E00491:102:H23N2CCXY:7:2114:4827:50568 99 17 61565813 60 19S92M39S = 61565891 181 AGGTGTCTGCAGCATGTGGCCCCAGGCCGGGGACTCTGTAAGCCACTGCTGGAGAGCCACTCCCATCCTTTCTCCCATTTCTCTAGACCTGCTGCCTATACAGTCACTTTTTTTTTTTTTTTGAGACGGAGTCTCGCTCTGTCGCCCAGG AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFFJJFJJJJJFJJJJJJJJJJJJJJJJFFAJJJJJJJ
JJJFJFJJJFJFJJFJJJJJJJJJJJJF--7-7<-7A7<-<7-)77F7AFAFF)A< AS:i:111 XS:i:24 RG:Z:GWI127_combined SA:Z:6,17717274,+,103S47M,0,0;
...

$samtools view /data/NGS/analysis/bamclipper/GWI127_combined.primerclipped.bam | grep ATACAGTCACTTTTTTTTTTTTTTTGAGACGGAGTCTCGCTCTGTCGCCC| wc -l
21

VCF file:

CHROM POS ID REF ALT QUAL FILTER INFO FORMAT GWI127_combined

17 61565890 rs1799752 T TATACAGTCACTTTTTTTTTTTTTTTGAGACGGAGTCTCGCTCTGTCGCCC 0 LowQual AC=0;AF=0.00;AN=2;DB;DP=71;ExcessHet=3.0103;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=59.99;SOR=0.368 GT:AD:DP:GQ:PL:SB 0/0:3,0:3:10:0,10,137:2,1,0,0

Issue · Github
by Sheila

Issue Number
2175
State
closed
Last Updated
Assignee
Array
Milestone
Array
Closed By
chandrans

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @asaki
    Hi,

    Can you please post IGV screenshots of the original BAM file and bamout file in the region?

    Thanks,
    Sheila

  • asakiasaki chinaMember

    imageimage

    Hi Shelia @Sheila
    It's weird that I can not see sequences containing insert alleles, which can be observed in bam file.
    Plz take a look at the screenshot of the original and bamout file.

    Thanks

  • asakiasaki chinaMember

    I also submit the zip file as "GWI127.rs1799752.bam.tar.gz", plz take a look at your convenience.

    Thanks

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @asaki
    Hi,

    I am going to check with the team and get back to you.

    -Sheila

  • asakiasaki chinaMember

    @Sheila
    Hi Sheila, any updates?

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @asaki
    Hi,

    I am just about to put in a bug report so one of the developers has a look. It looks like something else may be going on there, because of the odd coverage in the region (perhaps a copy number event). I will let you know when the developers have had a look.

    -Sheila

    Issue · Github
    by Sheila

    Issue Number
    3149
    State
    open
    Last Updated
  • asakiasaki chinaMember

    @Sheila

    Thanks Sheila.
    We happened to run the gel to get the genotypes for this variant.

    Two methods have been tested these days, since I found the BadMate reads in this region.

    • I tried to add the -rf BadMate option in command line, the variants can be called, however, 18 of 36 samples showed disconcordant genotypes.
    java -jar $gatk -T UnifiedGenotyper -R $ref -stand_call_conf 10 -D $vcf --genotyping_mode GENOTYPE_GIVEN_ALLELES --alleles $vcf --output_mode EMIT_ALL_SITES -A StrandBiasBySample -A AlleleBalanceBySample -dt NONE -I $dir/$bam -o $accession.rs1799752.ug.vcf -L 17:61565790-61565990 -rf BadMate -rf BadCigar
    
    • I tried to extract forward sequences with matched primer sequences only, 11 of 36 samples showed disconcordant genotypes.
    samtools view -h $dir/$bam 17:61565813-61565904 | awk '{if(\$10~/^AGGTGTCTGCAGCATGTGG/||\$0~/@/){print}}' > $accession.rs1799752.sam
    samtools view -bS $accession.rs1799752.sam -o $accession.rs1799752.bam
    samtools index $accession.rs1799752.bam
    java -jar /bioinfo/software/bin/GenomeAnalysisTK.jar -T UnifiedGenotyper -R /bioinfo/data/iGenomes/Homo_sapiens/NCBI/build37.2/Sequence/WholeGenomeFasta/genome.fa -stand_call_conf 10 -D $vcf --genotyping_mode GENOTYPE_GIVEN_ALLELES --alleles $vcf --output_mode EMIT_ALL_SITES -dt NONE -I $accession.rs1799752.bam -o $accession.rs1799752.ug.vcf -L 17:61565790-61565990
    

    I found that many genotype mistakes are due to ABHet ratio. Then I did the genotype adjust, by converting 0-0.1 as hom_ref, 0.4-0.6 as Het, 0.9-1.0 as hom_alt.

    • option 1, 15 out of 36 samples showed disconcordance
    • option 2, only 1 of 36 samples showed disconcordance

    Although option 2 can largely increase the concordance, I am still wondering if there is any other better option that I can use GATK directly to call the variants?

    Thanks

    Junfeng

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @asaki
    Hi Junfeng,

    Are you saying that the exact indel you are looking for is called when you add -rf BadMate to your command? Did you try it with both HaplotypeCaller and UnifiedGenotyper?

    Thanks
    Sheila

  • asakiasaki chinaMember

    @Sheila

    Yes, adding -rf BadMate can make calls, however, most of the results are disconcordant with gel results.

    Both HC and UG are applied to directly call the variants, which result in similar genotype callings.

    The only method I can find now is first extract forward sequences only and make callings.
    Hope there is a convenient method to directly call the variants from original bam.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @asaki
    Hi,

    I am waiting for the developers to add some input. Let me see what they say.

    -Sheila

Sign In or Register to comment.