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

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

Sign In or Register to comment.