Detecting called Indels with low read support on both sides?

Hi,

I am mapping chimpanzee samples to the human reference hg19. I mappend the samples using the standard protocol (BWA mem, remove duplicates, indel realigner) and called them with GATK 3.7 Haplotype Caller. After all variant filtering (hard filter + remove duplicated and low mappability regions from external bed files), I found an interesting insertion in one of my samples:

In chr2:48033272 there is a deletion of this sequence TTTTTGTTTTAATTCCT . The human reference has GCA|TTTTTGTTTTAATTCCTT|TTTTGTTTTAATTCCTT|TG this sequence duplicated. This sample is called homozygous for the deletion.

A few bp after this, GATK calls an insertion:
chr2 48033352 . C CAACCGATGTTGCTTTTCTGTCCTAGCATTTTTGTTTTAATTCCTT 108.02 PASS

Long story short:

  • There are only 6 reads supporting this insertion.
  • Of them, only 3 have the full "GATK-ALT-insertion". All of these 3 have, at least, 1 bp more.
  • None of these reads have the 3' side of the reference. It should be: CTTTAACAGGAAGAGGTAC ins TGCAACATTTGATGGG
  • I lied. One of these reads does have the full sequence:
    TAACAGGAAGAGGTAC | AACCGATGTTGCTTTTCTGTCCTAGCATTTTTGTTTTAATTCCTTTGAGTTACTTCCTTATGCATATTTTACTTTAACAGGAAGAGGTAC | TGCAACATTTGATGGGACAGCAATAGCAAATGCAGTTGTTAAAGA

It is a duplication of the whole previous sequence, including the deletion 80bp upstream. I want to run functional analysis of the variants detected, and I am changing from a frameshift insertion to a non-frameshift insertion.

Ok, I have detected this wrong indel, but I am calling 1.9M Indels in this dataset, and 1.7M more in another one and I am worried about reporting strong functional annotation to erroneous variants.

Is there any method I can use to detect this kind of indel (not enough reads supporting both tips of the insertion)? Or I can only filter by QD?
My hard filter removes QUAL<50 and QD<2 . This particular variant has QUAL=108.02, QD=4.91 and Genotyping Quality=99

Any advice?

Thanks in advance,

Txema

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @txemaheredia
    Hi Txema,

    Interesting. Can you post IGV screenshots of the BAM and bamout files here? Please include ~300 bases before and after the indel.

    Thanks,
    Sheila

  • Hi Sheila,

    Many of the reads involved in the insertion have also the deletion at 5'. The images attached do not fully show the insertion because it is marked as soft clipped bp in the bam file.

    samtools tview

    IGV +- 300 bp

    IGV close-up

    CIGAR strings, with end of reference & beginning of insertion highlighted

    And the reads + cigar in plain text for your convenience:

    29M17D63M8S CTAACCGATGTTGCTTTTCTGTCCTAGCATTTTTGTTTTAATTCCTTTGAGTTACTTCCTTATGCATATTTTACTTTAACAGGAAGAGGTACAACCGATG
    27M17D73M AACCGATGTTGCTTTTCTGTCCTAGCATTTTTGTTTTAATTCCTTTGAGTTACTTCCTTATGCATATTTTACTTTAACAGGAAGAGGTACTGCAACATTT
    16S27M17D108M TAACAGGAAGAGGTACAACCGATGTTGCTTTTCTGTCCTAGCATTTTTGTTTTAATTCCTTTGAGTTACTTCCTTATGCATATTTTACTTTAACAGGAAGAGGTACTGCAACATTTGATGGGACAGCAATAGCAAATGCAGTTGTTAAAGA
    22M17D129M ATGTTGCTTTTCTGTCCTAGCATTTTTGTTTTAATTCCTTTGAGTTACTTCCTTATGCATATTTTACTTTAACAGGAAGAGGTACTGCAACATTTGATGGGACAGCAATAGCAAATGCAGTTGTTAAAGAACTTGCTGAGACTATAAAATG
    19M17D63M18S TTGCTTTTCTGTCCTAGCATTTTTGTTTTAATTCCTTTGAGTTACTTCCTTATGCATATTTTACTTTAACAGGAAGAGGTACAACCGATGTTGCTTTTCT
    11S63M26S CTGTCCTAGCATTTTTGTTTTAATTCCTTTGAGTTACTTCCTTATGCATATTTTACTTTAACAGGAAGAGGTACAACCGATGTTGCTTTTCTGTCCTAGC
    14S63M23S TTTCTGTCCTAGCATTTTTGTTTTAATTCCTTTGAGTTACTTCCTTATGCATATTTTACTTTAACAGGAAGAGGTACAACCGATGTTGATTTTCTGTCCT
    6S94M CTAGCATTTTTGTTTTAATTCCTTTGAGTTACTTCCTTATGCATATTTTACTTTAACAGGAAGAGGTACTGCAACATTTGATGGGACAGCAATAGCAAAT
    100M TGTTTTAATTCCTTTGAGTTACTTCCTTATGCATATTTTACTTTAACAGGAAGAGGTACTGCAACATTTGATGGGACAGCAATAGCAAATGCAGTTGTTA
    54M46S TAATTCCTTTGAGTTACTTCCTTATGCATATTTTACTTTAACAGGAAGAGGTACAACCGATGTTGCTTTTCTGTCCTAGCATTTTTGTTTTAATTCCTTT
    100M AATTCCTTTGAGTTACTTCCTTATGCATATTTTACTTTAACAGGAAGAGGTACTGCAACATTTGATGGGACAGCAATAGCAAATGCAGTTGTTAAAGAAC
    100M CTTTGAGTTACTTCCTTATGCATATTTTACTTTAACAGGAAGAGGTACTGCAACATTTGATGGGACAGCAATAGCAAATGCAGTTGTTAAAGAACTTGCT
    48M52S CTTTGAGTTACTTCCTTATGCATATTTTACTTTAACAGGAAGAGGTACAACCGATGTTGCTTTTCTGTCCTAGCATTTTTGTTTTAATTCCTTTGAGTTA
    100M TTACTTCCTTATGCATATTTTACTTTAACAGGAAGAGGTACTGCAACATTTGATGGGACAGCAATAGCAAATGCAGTTGTTAAAGAACTTGCTGAGACTA
    100M CATATTTTACTTTAACAGGAAGAGGTACTGCAACATTTGATGGGACAGCAATAGCAAATGCAGTTGTTAAAGAACTTGCTGAGACTATAAAATGTCGTAC
    100M ATTTTACTTTAACAGGAAGAGGTACTGCAACATTTGATGGGACAGCAATAGCAAATGCAGTTGTTAAAGAACTTGCTGAGACTATAAAATGTCGTACATT
    100M TTTTACTTTAACAGGAAGAGGTACTGCAACATTTGATGGGACAGCAATAGCAAATGCAGTTGTTAAAGAACTTGCTGAGACTATAAAATGTCGTACATTA
    100M CAGGAAGAGGTACTGCAACATTTGATGGGACAGCAATAGCAAATGCAGTTGTTAAAGAACTTGCTGAGACTATAAAATGTCGTACATTATTTTCAACTCA
    100M GGAAGAGGTACTGCAACATTTGATGGGACAGCAATAGCAAATGCAGTTGTTAAAGAACTTGCTGAGACTATAAAATGTCGTACATTATTTTCAACTCACT
    100M AGAGGTACTGCAACATTTGATGGGACAGCAATAGCAAATGCAGTTGTTAAAGAACTTGCTGAGACTATAAAATGTCGTACATTATTTTCAACTCACTACC
    100M GGTACTGCAACATTTGATGGGACAGCAATAGCAAATGCAGTTGTTAAAGAACTTGCTGAGACTATAAAATGTCGTACATTATTTTCAACTCACTACCATT
    100M GTACTGCAACATTTGATGGGACAGCAATAGCAAATGCAGTTGTTAAAGAACTTGCTGAGACTATAAAATGTCGTACATTATTTTCAACTCACTACCATTC

    Thanks for your help

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @txemaheredia
    Hi Txema,

    Are these the original BAM files? It looks like there might be a structural variant present which HaplotypeCaller is not able to detect. You may consider using a structural variant caller like GenomeSTRiP.

    Also, can you post the bamout?

    Thanks,
    Sheila

  • Thanks Sheila,

    I just wanted to call snps and small indels in my samples, I am not planning on resolving complex structural variants. It is not that I am interested in resolving this indel but in detecting/filtering when an indel call has gone wrong/little support to both sides.

    I attached the bamouts for the calling of ±300bp. I've added .txt to the binary files so I could upload them, but they are still binaries. The samout file is the text version.

    Cheers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @txemaheredia
    Hi Txema,

    I cannot view the bamouts because they are in txt format. But, can you post the VCF records for those sites? I may need you to submit a bug report.

    Thanks,
    Sheila

  • Hi Sheila,

    The bamouts I uploaded are the original bamout binary files, but with ".txt" added to the filename (so the forum lets me upload them). I've just downloaded them and renamed to .bamout / .bamout.bai and I was able to open them fine with samtools

    Do you still want me to send the VCF?

    Txema

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @txemaheredia
    Hi Txema,

    No, I just need to see the VCF records for those indel sites. I am curious what some of the other annotations look like.

    -Sheila

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @txemaheredia
    Hi Txema,

    I was hoping to see the final VCF with the final annotations, but this is getting a bit tough to follow. It might just be best if you submit a bug report so I can take a look locally. Instructions are here.

    Thanks,
    Sheila

Sign In or Register to comment.