Genome Analysis Toolkit

Variant Discovery in High-Throughput Sequencing Data

GATK process banner

Need Help?

Search our documentation

Community Forum

Hi, How can we help?

Developed in the Data Sciences Platform at the Broad Institute, the toolkit offers a wide variety of tools with a primary focus on variant discovery and genotyping. Its powerful processing engine and high-performance computing features make it capable of taking on projects of any size. Learn more

How to recover a mutation from secondary alignment

0

8 comments

  • Avatar
    David Benjamin

    Laurent MANCHON The command line given above filters secondary and supplementary alignments (it turns on the corresponding read filters, which is also the default behavior).  The program log most likely reports no reads filtered in this way because they fail the mapping quality filter first.

    However, even if you disable these filters you will end up with a lot of mapping errors and double-counted variants.  Simply running Mutect2 would not produce sensible output.

    0
    Comment actions Permalink
  • Avatar
    David Benjamin

    By the way, to turn off a read filter, use the --disable-read-filter argument.

    0
    Comment actions Permalink
  • Avatar
    SkyWarrior

    Hi Laurent MANCHON

    Same region in my bam files do not contain any MAPQ0 reads. Could you be using a reference genome without unmasked alternate and unmapped contigs ? Can you try the very same analysis using GIAB hg38v2 reference genome?

    https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/references/GRCh38/GCA_000001405.15_GRCh38_no_alt_analysis_set_maskedGRC_exclusions_v2.fasta.gz 

    0
    Comment actions Permalink
  • Avatar
    Laurent MANCHON

    that's what I did in my analysis, I used the genome GRCh38_p13 (08.2022) from gencode consortium. I have removed all the alts contigs and use only the chromosome sequences.

     

    0
    Comment actions Permalink
  • Avatar
    SkyWarrior

    Interesting. Did you use bowtie2 with end-to-end mapping settings to align your reads ? 

    0
    Comment actions Permalink
  • Avatar
    Laurent MANCHON

    no.
    this is the command i use:

    bwa-mem2 mem -M -w 150 -d 80 -c 2000 -t 24 -R "@RG\tID:${lib}\tSM:${lib}\tLB:${lib}\tPL:ILLUMINA" $GENOME/GRCh38_p13.fa $f1 $f2 | samtools fixmate -O bam -@ 24 - - | sambamba-0.8 sort -t 24 -m 128G -o $OUTPUT_BAM/$lib".bam" /dev/stdin &>
     $BWA_LOG/$lib.log

    $f1 and $f2 are the input fastq file from my paired-end library.

     

    0
    Comment actions Permalink
  • Avatar
    SkyWarrior

    Even more interesting because those reads seem to have aligned also to another location on chromosome 21 with MAPQ0. I use bwa mem 0.7.17 with default settings using GIAB hg38v2 genome and the very same locus does not have any reads with XA tags. I was wondering if you have a chance to test it using default settings and GIAB hg38v2 genome (Gencode genomes are not masked like GIAB reference genomes. Gencode genomes are more suitable for RNASeq applications which require unmasked regions sometimes). I also see that you are using bwa-mem2 which is a vector accelerated version. github page says outputs are the same as regular bwa however I am concerned that there may be some discrepancies. 

    So here are a couple of things to test

    1- Using GIAB hg38v2 

    2- Using bwa mem with regular settings 

    3- All of the above. 

    If any of the parameters could save those reads to to go MAPQ0 then that is the most likely culprit here. I do not recommend removing read filters or accepting MAPQ0 reads into variant calling which will open up a whole bunch of can of worms later. 

    0
    Comment actions Permalink
  • Avatar
    Laurent MANCHON

    for now this problem is over for me because i have added those parameters to Mutect2:

    --read-filter NotSecondaryAlignmentReadFilter \
    --read-filter NotSupplementaryAlignmentReadFilter \
    --interval-padding 10 \

    also i use VarDict in collaboration with Mutect2 and i combine the results of these two callers.
    sometimes VarDict rescue some SNV not identified by Mutect2.

    best

    0
    Comment actions Permalink

Please sign in to leave a comment.

Powered by Zendesk