实战之 mutect2检测somatic变异流程

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-

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

最后编辑于：2020.04.10 10:16:47

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 158,736评论 4赞 362
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 67,167评论 1赞 291
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 108,442评论 0赞 243
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 43,902评论 0赞 204
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 52,302评论 3赞 287
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 40,573评论 1赞 216
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 31,847评论 2赞 312
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 30,562评论 0赞 197
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 34,260评论 1赞 241
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 30,531评论 2赞 245
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 32,021评论 1赞 258
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 28,367评论 2赞 253
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 33,016评论 3赞 235
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 26,068评论 0赞 8
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 26,827评论 0赞 194
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 35,610评论 2赞 274
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 35,514评论 2赞 269