Mutect2 — GATK4肿瘤体细胞变异检测工具¶

# Mutect2是GATK4的内置工具，安装GATK4即可
conda create -n mutect2 -c bioconda -c conda-forge gatk4=4.6.0.0
conda activate mutect2

# 验证
gatk Mutect2 --help   # 查看Mutect2帮助

# 下载配套资源数据库（需要GATK资源包）
# gnomAD生殖系资源（过滤germline变异）
# Panel of Normals（PON，正常样本面板，去背景噪音）
# GATK资源下载地址：https://gatk.broadinstitute.org/hc/en-us/articles/360035890811

# 第一步：调用体细胞变异（tumor+normal配对）
gatk Mutect2 \
    -R reference.fasta \                # 参考基因组
    -I tumor.bam \                      # 肿瘤样本BAM
    -I normal.bam \                     # 配对正常样本BAM
    -normal NORMAL_SAMPLE_NAME \        # 正常样本名（BAM中SM标签）
    --germline-resource gnomad.vcf.gz \ # gnomAD生殖系变异资源
    --panel-of-normals pon.vcf.gz \     # 正常样本面板（可选但推荐）
    -O unfiltered.vcf.gz                # 未过滤的原始结果

# 第二步：学习读段方向性偏差（可选，提升精度）
gatk LearnReadOrientationModel \
    -I f1r2.tar.gz \                    # Mutect2输出的read orientation模型数据
    -O read-orientation-model.tar.gz    # 输出模型文件

# 第三步：计算正常污染水平
gatk GetPileupSummaries \
    -I tumor.bam \
    -V known_variants.vcf.gz \          # 已知变异位点
    -L known_variants.vcf.gz \
    -O tumor_pileups.table

gatk CalculateContamination \
    -I tumor_pileups.table \
    -matched normal_pileups.table \
    -O contamination.table              # 污染率估算表

# 第四步：过滤体细胞变异
gatk FilterMutectCalls \
    -R reference.fasta \
    -V unfiltered.vcf.gz \
    --contamination-table contamination.table \
    --orientation-bias-artifact-priors read-orientation-model.tar.gz \
    -O filtered.vcf.gz                  # 过滤后VCF

gatk Mutect2 \
    -R reference.fasta \
    -I tumor.bam \
    --germline-resource gnomad.vcf.gz \
    --panel-of-normals pon.vcf.gz \     # PON非常重要！替代配对正常样本
    -O tumor_only_unfiltered.vcf.gz

# 同样需要过滤步骤
gatk FilterMutectCalls \
    -R reference.fasta \
    -V tumor_only_unfiltered.vcf.gz \
    -O tumor_only_filtered.vcf.gz

# 第一步：对每个正常样本跑Mutect2（单样本模式）
for normal in normal1.bam normal2.bam normal3.bam; do
    sample=$(samtools view -H $normal | grep "@RG" | grep -oP "SM:\K[^\t]+")
    gatk Mutect2 \
        -R reference.fasta -I $normal \
        --max-mnp-distance 0 \         # PON构建时设0
        -O ${sample}_for_pon.vcf.gz
done

# 第二步：合并正常样本GVCF
gatk GenomicsDBImport \
    -V normal1_for_pon.vcf.gz \
    -V normal2_for_pon.vcf.gz \
    -V normal3_for_pon.vcf.gz \
    --genomicsdb-workspace-path pon_db \
    -L chr1 -L chr2               # 按染色体分批

# 第三步：创建PON VCF
gatk CreateSomaticPanelOfNormals \
    -R reference.fasta \
    -V gendb://pon_db \
    -O pon.vcf.gz

# 完整Mutect2流程（配对模式，hg38）

REF="hg38.fasta"
GNOMAD="af-only-gnomad.hg38.vcf.gz"
PON="somatic-hg38_1000g_pon.hg38.vcf.gz"

# 1. 运行Mutect2（同时输出read orientation数据）
gatk Mutect2 \
    -R $REF \
    -I tumor.bam -I normal.bam \
    -normal NORMAL \
    --germline-resource $GNOMAD \
    --panel-of-normals $PON \
    -O unfiltered.vcf.gz \
    --f1r2-tar-gz f1r2.tar.gz \   # 输出read方向数据
    -L capture_targets.bed         # WES限定目标区间

# 2. 学习read orientation模型（对FFPE样本特别重要）
gatk LearnReadOrientationModel -I f1r2.tar.gz -O rom.tar.gz

# 3. 计算污染率
gatk GetPileupSummaries -I tumor.bam \
    -V small_exac_common_3.hg38.vcf.gz \
    -L small_exac_common_3.hg38.vcf.gz \
    -O tumor_pileups.table
gatk GetPileupSummaries -I normal.bam \
    -V small_exac_common_3.hg38.vcf.gz \
    -L small_exac_common_3.hg38.vcf.gz \
    -O normal_pileups.table
gatk CalculateContamination -I tumor_pileups.table \
    -matched normal_pileups.table -O contamination.table

# 4. 过滤
gatk FilterMutectCalls -R $REF \
    -V unfiltered.vcf.gz \
    --contamination-table contamination.table \
    --orientation-bias-artifact-priors rom.tar.gz \
    -O filtered.vcf.gz

# 5. 提取PASS变异
bcftools view -f PASS filtered.vcf.gz -o final_somatic.vcf.gz -O z
bcftools stats final_somatic.vcf.gz | grep "number of SNPs"