1 用途

mutect2是GATK（The Genome Analysis Toolkit，是Broad Institute开发的用于二代重测序数据分析的一款软件，里面包含了很多有用的工具）中的一个工具，主要功能是检测体细胞变异，包含SNV和INDEL。GATK4中的mutect2是基于一种贝叶斯模型来检测变异的，这与 MuTect （ Cibulskis et al., 2013）不同。Mutect2 v4.1.0.0 之后的版本支持多个样本的joint calling。

GATK官网有关于mutect2 的详细介绍，也提供了最佳实践共大家参考。但对于新手来说更希望有拿来就能用的流程，而不需要研究官方给的资料。下面小编为大家整理了一套mutect2检测体细胞变异的流程（核心还是最佳实践。这里仅针对配对样本），共初学者们参考。默认上游已经完成了比对（比对推荐BWA），我们的分析将从对bam文件的BQSR校正开始。

2 实战

本文提供的代码经修改后（路径，样本名等）写到shell文件中，就可以直接运行。

2.1 BQSR校正

#!/bin/bash
SID=样本名
INPUT=/全路径/${SID}.bam
##数据库路径
REF_GENOME=/全路径/ucsc.hg19.fasta
DBSNP=/全路径/dbsnp_138.hg19.vcf
MILLS=/全路径/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf
INDEL=/全路径/1000G_phase1.indels.hg19.sites.vcf
###软件路径GATK_HOME=GATK软件路径
PICARD_HOME=Picard软件路径
TMP_DIR=Tmp路径
###修改bam文件中RG标签
time java -Djava.io.tmpdir=${TMP_DIR} -XX:ParallelGCThreads=1 -jar \
     ${PICARD_HOME}/picard.jar AddOrReplaceReadGroups \ I=${INPUT} \             
     O=${INPUT%.bam}.new.bam \
     RGID=${SID} \
     RGLB=${SID} \
     RGPL=ILLUMINA \
     RGPU=flowcell-barcode.lane \
     RGSM=${SID}
 rm -f ${INPUT}
 ##BQSR校正
time java -Djava.io.tmpdir=${TMP_DIR} -XX:ParallelGCThreads=1 -jar \
    ${GATK_HOME}/gatk-package-4.1.2.0-local.jar BaseRecalibrator \
    -R ${REF_GENOME} \
     --known-sites $DBSNP \
    --known-sites $MILLS \
     --known-sites $INDEL \
     -I ${INPUT%.bam}.new.bam \
     -O ${INPUT%.bam}.recal.table
time java -Djava.io.tmpdir=${TMP_DIR} -XX:ParallelGCThreads=1 -jar \
    ${GATK_HOME}/gatk-package-4.1.2.0-local.jar ApplyBQSR \
     -R ${REF_GENOME} \
     -I ${INPUT%.bam}.new.bam \
     --bqsr-recal-file ${INPUT%.bam}.recal.table \
     -O ${INPUT}
 rm -f ${INPUT%.bam}.new.bam ${INPUT%.bam}.recal.table

BQSR的校正对象是所有样本，也就是说对于成对（pair）样本中的tumor和normal都要进行校正。

2.2 生成PON(panel of normals)文件

先对所有的对照样本（normal）进行预处理（过滤掉配对reads不在同一条染色体上的reads）：

#!/bin/bash
SID=样本名#这里仅为normal样本
PANEL=捕获区间bed文件#全基因组测序不需要该参数
 INPUT=${SID}.bam#经过BQSR校正后的bam
OUTPUT=${SID}.vcf.gz
 #数据库
REF_GENOME=/全路径/ucsc.hg19.fasta
GNOMAD=/全路径/af-only-gnomad.raw.sites.b37.vcf.gz
 ###软件路径
GATK_HOME=GATK软件路径
 PICARD_HOME=Picard软件路径
 TMP_DIR=Tmp路径
 time java -Djava.io.tmpdir=${TMP_DIR} -XX:ParallelGCThreads=1 -jar \
    ${GATK_HOME}/gatk-package-4.1.2.0-local.jar Mutect2 \
     -R ${REF_GENOME} \ -I ${INPUT} \ -O ${OUTPUT} \
     -L ${PANEL} \#全基因组测序不需要改参数
     --germline-resource ${GNOMAD} \
     --max-mnp-distance 0 \
     --disable-read-filter MateOnSameContigOrNoMappedMateReadFilter

接着生成PON文件：

#!/bin/bash
SAMPLE_MAP=对照样本经过预处理后所有的${SID}.vcf.gz文件的全路径列表，每个样本一行
PANEL=捕获区间bed文件#全基因组测序不需要该参数
OUTPUT=pon.${PANEL}.vcf.gz
#数据库
REF_GENOME=/全路径/ucsc.hg19.fasta
GNOMAD=/全路径/af-only-gnomad.raw.sites.b37.vcf.gz
 ###软件路径
GATK_HOME=GATK软件路径
 PICARD_HOME=Picard软件路径
 TMP_DIR=Tmp路径

#数据导入
 time java -Djava.io.tmpdir=${TMP_DIR} -XX:ParallelGCThreads=4 -jar \
    ${GATK_HOME}/gatk-package-4.1.2.0-local.jar GenomicsDBImport \
     -R ${REF_GENOME} \
     --sample-name-map ${SAMPLE_MAP} \
     --genomicsdb-workspace-path ${SAMPLE_MAP%.txt} \
     -L chr1 -L chr2 -L chr3 -L chr4 -L chr5 -L chr6 -L chr7 -L chr8 -L chr9 \
     -L chr10 -L chr11 -L chr12 -L chr13 -L chr14 -L chr15 -L chr16 -L chr17 \
     -L chr18 -L chr19 -L chr20 -L chr21 -L chr22 -L chrX -L chrY
    -L chrM \
     --batch-size 50 \
     --reader-threads 5

#PON生成
 time java -Djava.io.tmpdir=${TMP_DIR} -XX:ParallelGCThreads=4 -jar \
    ${GATK_HOME}/gatk-package-4.1.2.0-local.jar CreateSomaticPanelOfNormals \
     -R ${REF_GENOME} \
     --germline-resource ${GNOMAD} \
     -V gendb://${SAMPLE_MAP%.txt} \
     -O ${OUTPUT}

2.3 somatic变异检测与变异过滤

#!/bin/bash
TUMOR=$1 #tumor样本名
NORMAL=$2 #normal样本名
PANEL=捕获区间bed文件#全基因组测序不需要该参数
PON=/全路径/pon.$(PANEL).vcf.gz
TUMOR_BAM=/全路径/${TUMOR}.bam #经过BQSR校正后的bam
NORMAL_BAM=/全路径/${NORMAL}.bam #经过BQSR校正后的bam
RAW_VCF=/全路径/raw/${TUMOR}.vcf.gz
OUTPUT=/全路径/${TUMOR}.vcf.gz

#数据库,涉及到的数据库自己去相应官网下载，漫漫科研路会频繁使用的
REF_GENOME=/全路径/ucsc.hg19.fasta
GNOMAD=/全路径/af-only-gnomad.raw.sites.b37.vcf.gz
EXAC=/全路径/small_exac_common_3_b37.vcf.gz
BWA_IMAGE=/全路径/GRCh38.primary.index.img

###软件路径
GATK_HOME=GATK软件路径
 PICARD_HOME=Picard软件路径
 TMP_DIR=Tmp路径

 # Mutect2
time java -Djava.io.tmpdir=${TMP_DIR} -XX:ParallelGCThreads=1 -jar \
    ${GATK_HOME}/gatk-package-4.1.2.0-local.jar Mutect2 \
     -R ${REF_GENOME} \
     -I ${TUMOR_BAM} \
     -I ${NORMAL_BAM} \
     -normal ${NORMAL} \
     -L ${PANEL} \
     --panel-of-normals ${PON} \
     --germline-resource ${GNOMAD} \
     --disable-read-filter MateOnSameContigOrNoMappedMateReadFilter \
     --f1r2-tar-gz ${TUMOR_BAM%.bam}.f1r2.tar.gz \
     --bam-output ${TUMOR_BAM%.bam}.bamout.bam \
     -O ${RAW_VCF}

 # Estimate cross-sample contamination
    time java -Djava.io.tmpdir=${TMP_DIR} -XX:ParallelGCThreads=1 -jar \
    ${GATK_HOME}/gatk-package-4.1.2.0-local.jar GetPileupSummaries \
     -R ${REF_GENOME} \
     -I ${TUMOR_BAM} \
     -V ${EXAC} \
     -L ${EXAC} \
     -L ${PANEL} \
     --interval-set-rule INTERSECTION \
     -O ${TUMOR_BAM%.bam}.getpileupsum.table

time java -Djava.io.tmpdir=${TMP_DIR} -XX:ParallelGCThreads=1 -jar \
    ${GATK_HOME}/gatk-package-4.1.2.0-local.jar GetPileupSummaries \
     -R ${REF_GENOME} \
     -I ${NORMAL_BAM} \
     -V ${EXAC} \
     -L ${EXAC} \
     -L ${PANEL} \
     --interval-set-rule INTERSECTION \
     -O ${NORMAL_BAM%.bam}.getpileupsum.table

 time java -Djava.io.tmpdir=${TMP_DIR} -XX:ParallelGCThreads=1 -jar \
    ${GATK_HOME}/gatk-package-4.1.2.0-local.jar CalculateContamination \
     -I ${TUMOR_BAM%.bam}.getpileupsum.table \
     -matched ${NORMAL_BAM%.bam}.getpileupsum.table \
     -O ${TUMOR_BAM%.bam}.contamination.table \
     --tumor-segmentation ${TUMOR_BAM%.bam}.segmentation.table
rm  -f ${TUMOR_BAM%.bam}.getpileupsum.table ${NORMAL_BAM%.bam}.getpileupsum.table
rm -f ${TUMOR_BAM} ${NORMAL_BAM} ${TUMOR_BAM%.bam}.bai ${NORMAL_BAM%.bam}.bai
 
#计算Orientation bias模型参数
time java -Djava.io.tmpdir=${TMP_DIR} -XX:ParallelGCThreads=1 -jar \     ${GATK_HOME}/gatk-package-4.1.2.0-local.jar LearnReadOrientationModel \
     -I ${TUMOR_BAM%.bam}.f1r2.tar.gz \
     -O ${TUMOR_BAM%.bam}.artifact.tar.gz
 rm -f ${TUMOR_BAM%.bam}.f1r2.tar.gz
 # Filter for confident somatic calls

time java -Djava.io.tmpdir=${TMP_DIR} -XX:ParallelGCThreads=1 -jar \
    ${GATK_HOME}/gatk-package-4.1.2.0-local.jar FilterMutectCalls \
     -R ${REF_GENOME} \ -V ${RAW_VCF} \
     --contamination-table ${TUMOR_BAM%.bam}.contamination.table \
     --tumor-segmentation ${TUMOR_BAM%.bam}.segmentation.table \
     --ob-priors ${TUMOR_BAM%.bam}.artifact.tar.gz \
     -stats ${RAW_VCF}.stats \
     -O ${RAW_VCF%.vcf.gz}.filtered.vcf.gz
 rm -f ${TUMOR_BAM%.bam}.segmentation.table ${TUMOR_BAM%.bam}.artifact.tar.gz
rm -f ${TUMOR_BAM%.bam}.contamination.table
rm -f ${RAW_VCF} ${RAW_VCF}.tbi ${RAW_VCF}.stats
 #过滤由于比对引入的Artifacts

time java -Djava.io.tmpdir=${TMP_DIR} -XX:ParallelGCThreads=1 -jar \
     ${GATK_HOME}/gatk-package-4.1.2.0-local.jar FilterAlignmentArtifacts \
     -R ${REF_GENOME} \ -V ${RAW_VCF%.vcf.gz}.filtered.vcf.gz \
     -I ${TUMOR_BAM%.bam}.bamout.bam \
     --bwa-mem-index-image ${BWA_IMAGE} \
     -O ${OUTPUT}
 rm -f ${RAW_VCF%.vcf.gz}.filtered.vcf.gz ${RAW_VCF%.vcf.gz}.filtered.vcf.gz.tbi
 rm -f ${TUMOR_BAM%.bam}.bamout.bam ${TUMOR_BAM%.bam}.bamout.bai 

2.4 变异LeftAlign、indel碱基Trim、拆分复等位变异,提取

#!/bin/bash
 INPUT=上一步生成的变异文件（*.vcf.gz）
OUTPUT=/全路径/pass/${INPUT##*/}
#数据库,涉及到的数据库自己去相应官网下载，漫漫科研路会频繁使用的REF_GENOME=/全路径/ucsc.hg19.fasta
###软件路径
GATK_HOME=GATK软件路径 PICARD_HOME=Picard软件路径 TMP_DIR=Tmp路径
 
time java -Djava.io.tmpdir=${TMP_DIR} -XX:ParallelGCThreads=4 -jar \
    ${GATK_HOME}/gatk-package-4.1.2.0-local.jar LeftAlignAndTrimVariants \
     -R ${REF_GENOME} \
    -V ${INPUT} \
     -O ${OUTPUT}.biallele.vcf \
     --split-multi-allelics
 cat ${OUTPUT}.biallele.vcf | egrep '^##fileformat|^##FORMAT|^##INFO|^##contig' > ${OUTPUT}.vcf
cat ${OUTPUT}.biallele.vcf | grep -v '^##' | cut -f1-9,11 >> ${OUTPUT}.vcf
 rm -f ${OUTPUT}.biallele.vcf ${OUTPUT}.biallele.vcf.idx

 time java -Djava.io.tmpdir=${TMP_DIR} -XX:ParallelGCThreads=4 -jar \
    ${GATK_HOME}/gatk-package-4.1.2.0-local.jar SelectVariants \
     -V ${OUTPUT}.vcf \
     -O ${OUTPUT} \
     --exclude-non-variants \
     --exclude-filtered
 rm -f ${OUTPUT}.vcf

经过上述过程已经将高质量的somatic变异检测及提取出来啦~~

感谢流程原始撰写者WM大侠~

---

对上述流程进行简单总结：

最后为初学者提供一些关于mutec2 的介绍，虽然都是来自GATK官网，但毕竟官网辣么大，慢慢遨游也需要些时间~

GATK4 Mutect2官网介绍：https://gatk.broadinstitute.org/hc/en-us/articles/360037593851-Mutect2

最佳实践之预处理：https://gatk.broadinstitute.org/hc/en-us/articles/360035535912

最佳实践之somatic SNV/indel检测：https://gatk.broadinstitute.org/hc/en-us/articles/360035894731-Somatic-short-variant-discovery-SNVs-Indels-

讨论社区：https://gatk.broadinstitute.org/hc/en-us/community/topics

FilterMutectCalls：https://gatk.broadinstitute.org/hc/en-us/articles/360037225412-

当然，实现一次分析不是目的，有时间还是要好好研究一下官网~