GATK Mutect2类介绍

GATK 的 Mutect2 是一个用于检测体细胞变异（somatic mutations）的工具，主要用于癌症研究。体细胞变异是指在个体发育过程中体细胞中发生的 DNA 变化，这些变异不遗传给下一代，但在癌症发展过程中起到重要作用。Mutect2 是 GATK 工具集的一部分，专门设计用于识别这些变异。

Mutect2功能与用途

主要功能包括：

1. 体细胞突变检测：Mutect2 可以检测癌症样本（如肿瘤样本）中的体细胞突变，并可以处理纯度不高的样本，甚至是低频突变。

2. 支持多种数据类型：支持不同的测序数据类型，包括全外显子组测序（WES）、全基因组测序（WGS）、以及其他靶向测序数据。

3. 支持单样本和多样本：Mutect2 可以对单一肿瘤样本或成对的肿瘤-正常样本进行分析。

4. 肿瘤正常配对分析：当提供了匹配的正常样本时，Mutect2 可以更准确地过滤掉常见的胚系变异，减少假阳性率。

5. 低频变异检测：即使在杂合率较低或样本纯度较差的情况下，Mutect2 仍然可以检测出罕见的突变。

Mutect2类继承关系

Mutect2 （实现apply虚方法：调用Mutect2Engine的callRegion方法）<- AssemblyRegionWalker（实现 traverse()，其中processReadShard方法调用apply虚方法） <- WalkerBase <- GATKTool（实现doWork（）含traverse虚方法） <- CommandLineProgram

GATK命令行程序的执行逻辑

每个最终的类如 Mutect2、HaplotypeCaller 会重写 onTraversalStart方法，onTraversalStart执行方法初始化engine等操作，包含各自的engine，通过apply方法 调用engine的callRegion方法
如 Mutect2 apply方法 <-  Mutect2Engine callRegion 方法

Mutect2源码

package org.broadinstitute.hellbender.tools.walkers.mutect;

import htsjdk.variant.variantcontext.writer.VariantContextWriter;
import org.broadinstitute.barclay.argparser.Argument;
import org.broadinstitute.barclay.argparser.ArgumentCollection;
import org.broadinstitute.barclay.argparser.CommandLineException;
import org.broadinstitute.barclay.argparser.CommandLineProgramProperties;
import org.broadinstitute.barclay.help.DocumentedFeature;
import org.broadinstitute.hellbender.cmdline.StandardArgumentDefinitions;
import org.broadinstitute.hellbender.cmdline.programgroups.ShortVariantDiscoveryProgramGroup;
import org.broadinstitute.hellbender.engine.*;
import org.broadinstitute.hellbender.engine.filters.ReadFilter;
import org.broadinstitute.hellbender.exceptions.UserException;
import org.broadinstitute.hellbender.tools.walkers.annotator.*;
import org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCallerArgumentCollection;
import org.broadinstitute.hellbender.tools.walkers.haplotypecaller.ReferenceConfidenceMode;
import org.broadinstitute.hellbender.transformers.ReadTransformer;
import org.broadinstitute.hellbender.utils.downsampling.MutectDownsampler;
import org.broadinstitute.hellbender.utils.downsampling.ReadsDownsampler;
import org.broadinstitute.hellbender.cmdline.ModeArgumentUtils;
import org.broadinstitute.hellbender.utils.variant.writers.SomaticGVCFWriter;

import java.io.File;
import java.util.*;

/**
 * <p>Call somatic short mutations via local assembly of haplotypes.
 * Short mutations include single nucleotide (SNA) and insertion and deletion (indel) alterations.
 * The caller uses a Bayesian somatic genotyping model that differs from the original MuTect by
 * <a href='http://www.nature.com/nbt/journal/v31/n3/full/nbt.2514.html'>Cibulskis et al., 2013</a>
 * and uses the assembly-based machinery of
 * <a href="https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_haplotypecaller_HaplotypeCaller.php">HaplotypeCaller</a>.
 * Of note, Mutect2 v4.1.0.0 onwards enables joint analysis of multiple samples.</p>
 *
 * <p>This tool is featured in the <i>Somatic Short Mutation calling Best Practice Workflow</i>.
 * See <a href="https://software.broadinstitute.org/gatk/documentation/article?id=11136">Tutorial#11136</a> for a
 * step-by-step description of the workflow and <a href="https://software.broadinstitute.org/gatk/documentation/article?id=11127">Article#11127</a>
 * for an overview of what traditional somatic calling entails. For the latest pipeline scripts, see the
 * <a href="https://github.com/broadinstitute/gatk/tree/master/scripts/mutect2_wdl">Mutect2 WDL scripts directory</a>.
 * For pipelines with example data, see the <a href="https://github.com/gatk-workflows/gatk4-somatic-snvs-indels">gatk-workflows repository</a>.
 * Although we present the tool for somatic calling, it may apply to other contexts, such as mitochondrial variant calling and detection of somatic mosaicism.</p>
 *
 * <p>Starting with v4.1.0.0 Mutect2 accomodates extreme high depths, e.g. 20,000X. See the following articles for details on this and additional applications.
 * <ul>
 *   <li><a href="https://software.broadinstitute.org/gatk/blog?id=23400">Blog#23400</a> details general improvements to
 *   Mutect2 v4.1.0.0.</li>
 *   <li><a href="https://software.broadinstitute.org/gatk/blog?id=23598">Blog#23598</a> details Mutect2 mitochondrial
 *   mode.</li>
 *   <li><a href="https://software.broadinstitute.org/gatk/blog?id=XXX">Blog#XXX</a> (link to come) details use of Mutect2 in extremely
 *   low allele fraction variant detection.</li>
 * </ul></p>
 *
 * <h3>Usage examples</h3>
 * <p>Example commands show how to run Mutect2 for typical scenarios. The three modes are (i) tumor-normal mode where a
 * tumor sample is matched with a normal sample in analysis, (ii) tumor-only mode where a single sample's alignment
 * data undergoes analysis, and (iii) mitochondrial mode where sensitive calling at high depths is desirable.
 * <ul>
 *     <li>As of v4.1, there is no longer a need to specify the tumor sample name with -tumor. You need only specify
 *     the normal sample name with <nobr>-normal,</nobr> if you include a normal.</li>
 *     <li>Starting with v4.0.4.0, GATK recommends the default setting of --af-of-alleles-not-in-resource, which the
 *     tool dynamically adjusts for different modes. tumor-only calling sets the default to 5e-8, tumor-normal calling
 *     sets it to 1e-6 and mitochondrial mode sets it to 4e-3. For previous versions, the default was 0.001, the
 *     average heterozygosity of humans. For other organisms, change --af-of-alleles-not-in-resource to
 *     1/(ploidy*samples in resource).</li>
 * </ul></p>
 *
 * <h4>(i) Tumor with matched normal</h4>
 * <p>Given a matched normal, Mutect2 is designed to call somatic variants only. The tool includes logic to skip
 * emitting variants that are clearly present in the germline based on provided evidence, e.g. in the matched normal.
 * This is done at an early stage to avoid spending computational resources on germline events. If the variant's
 * germline status is borderline, then Mutect2 will emit the variant to the callset for subsequent filtering by
 * <a href="https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_mutect_FilterMutectCalls.php">FilterMutectCalls</a>
 * and review.</p>
 *
 * <pre>
 *     gatk Mutect2 \
 *     -R reference.fa \
 *     -I tumor.bam \
 *     -I normal.bam \
 *     -normal normal_sample_name \
 *     --germline-resource af-only-gnomad.vcf.gz \
 *     --panel-of-normals pon.vcf.gz \
 *     -O somatic.vcf.gz
 * </pre>
 *
 * <p>
 *     Mutect2 also generates a stats file names [output vcf].stats.  That is, in the above example the stats file would be named somatic.vcf.gz.stats
 *     and would be in the same folder as somatic.vcf.gz.  As of GATK 4.1.1 this file is a required input to FilterMutectCalls.
 * </p>
 *
 * <p> As of v4.1 Mutect2 supports joint calling of multiple tumor and normal samples from the same individual. The
 * only difference is that -I and <nobr>-normal</nobr> must be specified for the extra samples.</p>
 *
 * <pre>
 *     gatk Mutect2 \
 *     -R reference.fa \
 *     -I tumor1.bam \
 *     -I tumor2.bam \
 *     -I normal1.bam \
 *     -I normal2.bam \
 *     -normal normal1_sample_name \
 *     -normal normal2_sample_name \
 *     --germline-resource af-only-gnomad.vcf.gz \
 *     --panel-of-normals pon.vcf.gz \
 *     -O somatic.vcf.gz
 * </pre>
 *
 * <h4>(ii) Tumor-only mode</h4>
 * <p>This mode runs on a single type of sample, e.g. the tumor or the normal.
 * To create a PoN, call on each normal sample in this mode, then use
 * <a href="https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_mutect_CreateSomaticPanelOfNormals.php">CreateSomaticPanelOfNormals</a>
 * to generate the PoN.</p>
 *
 * <pre>
 *  gatk Mutect2 \
 *   -R reference.fa \
 *   -I sample.bam \
 *   -O single_sample.vcf.gz
 * </pre>
 *
 * <p>To call mutations on a tumor sample, call in this mode using a PoN and germline resource. After FilterMutectCalls
 * filtering, consider additional filtering by functional significance with
 * <a href="https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_funcotator_Funcotator.php">Funcotator</a>.</p>
 *
 * <pre>
 *  gatk Mutect2 \
 *  -R reference.fa \
 *  -I sample.bam \
 *  --germline-resource af-only-gnomad.vcf.gz \
 *  --panel-of-normals pon.vcf.gz \
 *  -O single_sample.vcf.gz
 * </pre>
 *
 * <h4>(iii) Mitochondrial mode</h4>
 * <p>Mutect2 automatically sets parameters appropriately for calling on mitochondria with the <nobr>--mitochondria</nobr> flag.
 * Specifically, the mode sets <nobr>--initial-tumor-lod</nobr> to 0, <nobr>--tumor-lod-to-emit</nobr> to 0, <nobr>--af-of-alleles-not-in-resource</nobr> to
 * 4e-3, and the advanced parameter <nobr>--pruning-lod-threshold</nobr> to -4.</p>
 *
 * <pre>
 *  gatk Mutect2 \
 *  -R reference.fa \
 *  -L chrM \
 *  --mitochondria-mode \
 *  -I mitochondria.bam \
 *  -O mitochondria.vcf.gz
 * </pre>
 *
 * <p>The mode accepts only a single sample, which can be provided in multiple files.</p>
 *
 * <h4>(iv) Force-calling mode</h4>
 * <p>This mode force-calls all alleles in force-call-alleles.vcf in addition to any other variants Mutect2 discovers.
 *
 * <pre>
 *  gatk Mutect2 \
 *   -R reference.fa \
 *   -I sample.bam \
 *   -alleles force-call-alleles.vcf
 *   -O single_sample.vcf.gz
 * </pre>
 *
 * <p>
 *     If the sample is suspected to exhibit orientation bias artifacts (such as in the case of FFPE tumor samples) one should also
 *     collect F1R2 metrics by adding an --f1r2-tar-gz argument as shown below.  This file contains information that can then be passed
 *     to LearnReadOrientationModel, which generate an artifact prior table for each tumor sample for FilterMutectCalls to use.
 *
 * </p>
 *
 * <pre>
 * gatk Mutect2 \
 *  -R reference.fa \
 *  -I sample.bam \
 *  --f1r2-tar-gz f1r2.tar.gz \
 *  -O single_sample.vcf.gz
 *  </pre>
 *
 * <h3>Notes</h3>
 * <ol>
 *     <li>Mutect2 does not require a germline resource nor a panel of normals (PoN) to run, although both are recommended. The tool prefilters sites
 *     for the matched normal and the PoN.
 *     <li>If a variant is absent from a given germline
 *     resource, then the value for --af-of-alleles-not-in-resource is used as an imputed allele frequency. Below is an excerpt of a known variants
 *     resource with population allele frequencies.
 *
 *     <pre>
 *     #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO
 *      1       10067   .       T       TAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCC      30.35   PASS    AC=3;AF=7.384E-5
 *      1       10108   .       CAACCCT C       46514.32        PASS    AC=6;AF=1.525E-4
 *      1       10109   .       AACCCTAACCCT    AAACCCT,*       89837.27        PASS    AC=48,5;AF=0.001223,1.273E-4
 *      1       10114   .       TAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAACCCTAACCCTAACCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAAACCCTA  *,CAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAACCCTAACCCTAACCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAAACCCTA,T      36728.97        PASS    AC=55,9,1;AF=0.001373,2.246E-4,2.496E-5
 *      1       10119   .       CT      C,*     251.23  PASS    AC=5,1;AF=1.249E-4,2.498E-5
 *      1       10120   .       TA      CA,*    14928.74        PASS    AC=10,6;AF=2.5E-4,1.5E-4
 *      1       10128   .       ACCCTAACCCTAACCCTAAC    A,*     285.71  PASS    AC=3,1;AF=7.58E-5,2.527E-5
 *      1       10131   .       CT      C,*     378.93  PASS    AC=7,5;AF=1.765E-4,1.261E-4
 *      1       10132   .       TAACCC  *,T     18025.11        PASS    AC=12,2;AF=3.03E-4,5.049E-5
 *     </pre>
 *     </li>
 *
 * <li>Additional parameters that factor towards filtering are available in
 * <a href="https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_mutect_FilterMutectCalls.php">FilterMutectCalls</a>.
 * While the tool calculates normal-lod assuming a diploid genotype, it calculates
 * normal-artifact-lod with the same approach it uses for tumor-lod, i.e. with a variable ploidy assumption.
 *
 * <ul>
 *     <li>If the normal artifact log odds becomes large, then FilterMutectCalls applies the artifact-in-normal filter..</li>
 *     <li>The tumor log odds, which is calculated independently of any matched normal, determines whether to filter a tumor
 *     variant. Variants with tumor LODs exceeding the threshold pass filtering.</li>
 * </ul></p>
 *</li>
 * </ol>
 */
 @CommandLineProgramProperties(
         summary = "Call somatic SNVs and indels via local assembly of haplotypes",
         oneLineSummary = "Call somatic SNVs and indels via local assembly of haplotypes",
         programGroup = ShortVariantDiscoveryProgramGroup.class
 )
@DocumentedFeature
public final class Mutect2 extends AssemblyRegionWalker {
     public static final String MUTECT_STATS_SHORT_NAME = "stats";
     public static final String DEFAULT_STATS_EXTENSION = ".stats";

    @ArgumentCollection
    protected M2ArgumentCollection MTAC = new M2ArgumentCollection();

    @Argument(fullName = StandardArgumentDefinitions.OUTPUT_LONG_NAME, shortName = StandardArgumentDefinitions.OUTPUT_SHORT_NAME, doc = "File to which variants should be written")
    public GATKPath outputVCF;

    private VariantContextWriter vcfWriter;

    private Mutect2Engine m2Engine;

    @Override
    public boolean useVariantAnnotations() { return true;}

    @Override
    public List<ReadFilter> getDefaultReadFilters() {
        return Mutect2Engine.makeStandardMutect2ReadFilters();
    }

    @Override
    public ReadTransformer makePostReadFilterTransformer() {
        return super.makePostReadFilterTransformer().andThen(Mutect2Engine.makeStandardMutect2PostFilterReadTransformer(referenceArguments.getReferencePath(), !MTAC.dontClipITRArtifacts));
    }

    @Override
    public List<Class<? extends Annotation>> getDefaultVariantAnnotationGroups() {
        return Mutect2Engine.getStandardMutect2AnnotationGroups();
    }

    @Override
    protected ReadsDownsampler createDownsampler() {
        return new MutectDownsampler(assemblyRegionArgs.maxReadsPerAlignmentStart, MTAC.maxSuspiciousReadsPerAlignmentStart, MTAC.downsamplingStride);
    }

    @Override
    public AssemblyRegionEvaluator assemblyRegionEvaluator() { return m2Engine; }

    @Override
    public boolean shouldTrackPileupsForAssemblyRegions() {
        return MTAC.pileupDetectionArgs.usePileupDetection;
    }

    @Override
    public void onTraversalStart() {
        VariantAnnotatorEngine annotatorEngine = new VariantAnnotatorEngine(makeVariantAnnotations(), null, Collections.emptyList(), false, false);
        m2Engine = new Mutect2Engine(MTAC, assemblyRegionArgs, createOutputBamIndex, createOutputBamMD5, getHeaderForReads(), referenceArguments.getReferenceSpecifier(), annotatorEngine);
        vcfWriter = createVCFWriter(outputVCF);
        if (m2Engine.emitReferenceConfidence()) {
            logger.warn("Note that the Mutect2 reference confidence mode is in BETA -- the likelihoods model and output format are subject to change in subsequent versions.");
            if ( MTAC.emitReferenceConfidence == ReferenceConfidenceMode.GVCF ) {
                try {
                    vcfWriter = new SomaticGVCFWriter(vcfWriter, new ArrayList<Number>(MTAC.GVCFGQBands));
                } catch ( IllegalArgumentException e ) {
                    throw new CommandLineException.BadArgumentValue("GQBands", "are malformed: " + e.getMessage());
                }
            }
        }

        if (MTAC.f1r2TarGz != null && !MTAC.f1r2TarGz.getAbsolutePath().endsWith(".tar.gz")) {
            throw new UserException.CouldNotCreateOutputFile(MTAC.f1r2TarGz, M2ArgumentCollection.F1R2_TAR_GZ_NAME + " file must end in .tar.gz");
        }
        m2Engine.writeHeader(vcfWriter, getDefaultToolVCFHeaderLines());
    }

    @Override
    public Collection<Annotation> makeVariantAnnotations(){
        final Collection<Annotation> annotations = super.makeVariantAnnotations();

        if (MTAC.mitochondria) {
            annotations.add(new OriginalAlignment());
        }

        return annotations;
    }

    @Override
    public Object onTraversalSuccess() {
        m2Engine.writeExtraOutputs(new File(outputVCF + DEFAULT_STATS_EXTENSION));

        return "SUCCESS";
    }

    @Override
    public void apply(final AssemblyRegion region, final ReferenceContext referenceContext, final FeatureContext featureContext ) {
        m2Engine.callRegion(region, referenceContext, featureContext).forEach(vcfWriter::add);
    }

    @Override
    public void closeTool() {
        if (vcfWriter != null) {
            vcfWriter.close();
        }
        if (m2Engine != null) {
            m2Engine.close();
        }
    }

    /**
     * mode adjustments
     * @return
     */
    @Override
    protected String[] customCommandLineValidation() {
        if (MTAC.flowMode != M2ArgumentCollection.FlowMode.NONE) {
            ModeArgumentUtils.setArgValues(
                    getCommandLineParser(),
                    MTAC.flowMode.getNameValuePairs(),
                    M2ArgumentCollection.FLOW_M2_MODE_LONG_NAME
                    );
        }
        return null;
    }
}

AssemblyRegionWalker源码

package org.broadinstitute.hellbender.engine;

import org.broadinstitute.barclay.argparser.Argument;
import org.broadinstitute.barclay.argparser.ArgumentCollection;
import org.broadinstitute.hellbender.engine.filters.CountingReadFilter;
import org.broadinstitute.hellbender.engine.filters.ReadFilter;
import org.broadinstitute.hellbender.engine.filters.ReadFilterLibrary;
import org.broadinstitute.hellbender.engine.filters.WellformedReadFilter;
import org.broadinstitute.hellbender.engine.spark.AssemblyRegionArgumentCollection;
import org.broadinstitute.hellbender.exceptions.UserException;
import org.broadinstitute.hellbender.utils.IGVUtils;
import org.broadinstitute.hellbender.utils.IntervalUtils;
import org.broadinstitute.hellbender.utils.SimpleInterval;
import org.broadinstitute.hellbender.utils.downsampling.PositionalDownsampler;
import org.broadinstitute.hellbender.utils.downsampling.ReadsDownsampler;

import java.io.IOException;
import java.io.PrintStream;
import java.util.ArrayList;
import java.util.Iterator;
import java.util.List;

/**
 * An AssemblyRegionWalker is a tool that processes an entire region of reads at a time, each marked as either "active"
 * (containing possible variation) or "inactive" (not likely to contain actual variation). Tool authors must implement
 * {@link #assemblyRegionEvaluator} to provide a means of determining whether a site is active or not, as well as
 * {@link #apply} to process each region. Authors must also implement methods providing default values for the
 * various traversal parameters.
 *
 * Each region passed to {@link #apply} will come pre-marked as either "active" or "inactive" using the results of the
 * configured {@link #assemblyRegionEvaluator}. The {@link #assemblyRegionEvaluator} is used to evaluate the probability
 * that each individual locus is active, and these probabilities are used to determine the bounds of each active and
 * inactive region.
 *
 * {@link #apply} will be called once for each active AND inactive region, and it is up to the implementation how to
 * handle/process active vs. inactive regions.
 *
 * Internally, the reads are loaded in chunks called read shards, which are then subdivided into active/inactive regions
 * for processing by the tool implementation. One read shard is created per contig.
 */
public abstract class AssemblyRegionWalker extends WalkerBase {

    @ArgumentCollection
    public final AssemblyRegionArgumentCollection assemblyRegionArgs = new AssemblyRegionArgumentCollection();

    /**
     * If provided, this walker will write out its assembly regions
     * to this file in the IGV formatted TAB-delimited output:
     *
     * http://www.broadinstitute.org/software/igv/IGV
     *
     * Intended to make debugging the active region calculations easier
     */
    @Argument(fullName = AssemblyRegionArgumentCollection.ASSEMBLY_REGION_OUT_LONG_NAME, doc="Output the assembly region to this IGV formatted file", optional = true)
    protected String assemblyRegionOut = null;

    private PrintStream assemblyRegionOutStream;

    @Override
    public final boolean requiresReads() { return true; }

    @Override
    public final boolean requiresReference() { return true; }

    @Override
    public String getProgressMeterRecordLabel() { return "regions"; }

    private List<MultiIntervalLocalReadShard> readShards;

    private boolean nonRandomDownsamplingMode = nonRandomDownsamplingMode();

    /**
     * Initialize data sources for traversal.
     *
     * Marked final so that tool authors don't override it. Tool authors should override onTraversalStart() instead.
     */
    @Override
    protected final void onStartup() {
        super.onStartup();

        assemblyRegionArgs.validate();

        final