【2.6.3】使用GATK4 Mutect2 检测体细胞突变

二、常规流程

2.1 Mutect2检测体细胞突变

输入文件：

1. 肿瘤样本名：HCC1143_tumor
2. 配对癌旁样本名：HCC1143_normal
3. PoN正常样本集VCF文件：后面给出生成PoN文件的方法
4. germline资源：germline变异数据库VCF
5. 目标区域文件：interval

输出文件：

1. 原始的体细胞突变的VCF文件：1_somatic_m2.vcf.gz(.tbi)
2. 重组装reads的BAM文件：2_tumor_normal_m2.bam(.bai)

命令行格式：

gatk --java-options "-Xmx2g" Mutect2 \
-R hg38/Homo_sapiens_assembly38.fasta \
-I tumor.bam \
-I normal.bam \
-tumor HCC1143_tumor \
-normal HCC1143_normal \
-pon resources/chr17_pon.vcf.gz \
--germline-resource resources/chr17_af-only-gnomad_grch38.vcf.gz \
--af-of-alleles-not-in-resource 0.0000025 \
--disable-read-filter MateOnSameContigOrNoMappedMateReadFilter \
-L chr17plus.interval_list \
-O 1_somatic_m2.vcf.gz \
-bamout 2_tumor_normal_m2.bam 

参数的解释：

• -I提供case的BAM文件，-tumor提供BAM文件对应的read group id（BAM头文件中的@RG->SM值）

• -I提供control的BAM文件，-normal提供相应的read group id.如果tumor有配对的normal，Mutect2可以据此排除germline变异和个体特异性人工产物；如果tumor不能提供配对normal样本，得到的体细胞变异文件会产生更多假阳性

• PoN(正常样本callset)定义了一个预过滤变异位点集。②展示了如何创建一个PoN文件。PoN不仅呈现了通常的germline变异位点，也呈现了测序数据引入的噪音，如测序bias或比对bias。

• 默认情况下，出现在PoN中的变异位点不认为是tumor样本的somatic变异也不会进行重组装。

• 参数–genotype-pon-sites可以得到PoN位点的基因型信息

• 如果tumor位点和PoN位点变异不是精确匹配，Mutect2会重新组装此区域，在变异文件中输出此位点，并在INFO列标记为“IN_PON”

• 通过使用–germline-resource指定人群的germline信息来注释变异等位基因。

此文件必须包含等位基因特异性频率，即它必须在INFO字段中包含AF注释，这是因为Mutect2需要用群体等位基因频率注释变异。

使用这个参数，要考虑调整参数–af-of-alleles-not-in-resource，默认值是0.001。

示例中的gnomAD资源af-only-gnomad_grch38.vcf.gz代表~200k个外显子组和~16k个基因组，示例中的数据是外显子组数据，因此我们将–af-of-alleles-not-in-resource调整为0.0000025=1 /（2 外显子组样本数）=1/(2200k)

默认值0.001是基于人平均杂合率给出的，适用于没有任何群体信息的人样本分析。

在评估一个变异为体细胞突变的概率时，人群等位基因频率和参数af-of-alleles-not-in-resource都将考虑在内。

1. 不过滤比对到不同contig上的PE reads，大概是因为BWA的PE比对策略（alt-aware mapping,unclear），如果滤掉不配对reads,会丢失一部分与备用单体型相关的SNV和Indel信息。
2. -L指定目标区域，-bamout生成重组装比对文件用以人工核查重组装区域和人造单体型，可选项

2.2 使用Mutect2的tumor-only模式创建panel of normal（下面一律简称PoN）

1.对每个case样本，运行tumor-only命令,以样本HG00190示例：

gatk Mutect2 \
-R ~/Documents/ref/hg38/Homo_sapiens_assembly38.fasta \
-I HG00190.bam \
-tumor HG00190 \
--disable-read-filter MateOnSameContigOrNoMappedMateReadFilter \
-L chr17plus.interval_list \
-O 3_HG00190.vcf.gz

整理上一步生成的VCFs，运行命令，举例用三个样本，实际上，推荐PoN至少需要40个样本

gatk CreateSomaticPanelOfNormals \
-vcfs 3_HG00190.vcf.gz \
-vcfs 4_NA19771.vcf.gz \
-vcfs 5_HG02759.vcf.gz \
-O 6_threesamplepon.vcf.gz

这一步生成文件包含在>=2个样本中出现的变异allele和位点

理想情况下，用以创建PoN的样本最好是有代表性的tumor case样本，即样本来自同样的测序平台，同样的实验处理，使用同样的分析工具，但实际上，即使是不匹配的PoN样本也可以显著过滤掉大龙的测序引入的假阳。这是因为比对bias和PCR错误在短片段测序平台上，一般更倾向于发生在同样的基因区域。

除了PoN创建，tumor-only方法也可以灵活应用于其他场景，比如一群人的变异检测，或者线粒体DNA，或者简单比较两个个体间差异。

2.3 如何评估样本间交叉污染

对肿瘤BAM运行GetPileupSummaries以总结tumor样本在已知变异位点集上的reads支持情况

gatk GetPileupSummaries \
-I tumor.bam \
-V resources/chr17_small_exac_common_3_grch38.vcf.gz \
-O 7_tumor_getpileupsummaries.table

输出文件是一个6列表

此工具只考虑样本中纯和备用位点：等位基因频率范围在0.01-0.2（相关参数是

–minimum-population-allele-frequency和

–maximum-population-allele-frequency）,这样设计的理论基础是：如果某个纯和备用位点的次佳allele频率较低，当发生样本交叉污染时，我们更容易观测到ref allele（或更常见allele）的出现，这样我们可以更准确地定量污染比例。

• 一个可以加快速度的参数，是用-L参数

• GATK4.0.8.0需要同时指定-L和-V

• 用CalculateContamination来估计污染比例

  gatk CalculateContamination
  -I 7_tumor_getpileupsummaries.table
  -O 8_tumor_calculatecontamination.table

从生成文件中，我们可以得到污染比例和测序错误率

有两种运行模式，另一种是提供normal样本，以更精确估计污染比例，第一步对normal样本进行GetPileupSummaries，然后使用CalculateContamination时指定-matched参数

1. 注意：样本交叉污染不同于tumor的normal污染和normal样本的tumor污染。目前，这个工具还不能用来计算算肿瘤的纯度。
2. 如果发现了高度污染，那么第一件事情就是核查样本是否在read group水平上发生了交换，可以用Picard的CrosscheckFingerprints来实现。

DNA测序需多路复用一个跨lane样本，并对同一样本的多个read group重分组，这有可能 导致不相关样本包含在一个read group中，从而在体细胞分析中引入错误。

Picard工具（没有更多细节）可以（1）核对肿瘤样本和配对癌旁是否来自同一个个体 （2）在read group水平上核对每个read group数据是否来自同一个个体

如果正常样本有16个read group，而癌症样本有22个read group，我们设置参数

EXPECT_ALL_GROUPS_TO_MATCH=true，

CrosscheckFingerprints会告诉我们所有read groups的相关性

2.4 用FilterMutectCalls过滤检测到突变集

这个工具用来过滤Mutect2产生的callsets，采用预先调好的阈值参数以适用于人体细胞 突变分析

gatk FilterMutectCalls \
-V somatic_m2.vcf.gz \
--contamination-table tumor_calculatecontamination.table \
-O 9_somatic_oncefiltered.vcf.gz

生成文件VCF.gz: 如果候选calls被判为真阳变异，FILTER列标记为PASS，如果候选calls被标记为假阳变异，FILTER列给出原因，可以在VCF头文件中查看(grep ##FILTER)

如果一个filter事件在输入文件中没有出现，这个工具会忽略这个filter事件，而在输出文件的头文件中还写入这个filter注释。

在未来GATK版本,FilterMutectCalls将使用一个统计模型来过滤污染

2.5 用CollectSequencingArtifactMetrics估计artifacts并用FilterByOrientationBias过滤

这一步是可选项，使用时应始终先执行FilterMutectCalls.

FilterByOrientationBias可以过滤序列背景的artifacts, 例如OxoG和FFPE.

该工具需要Picard CollectSequencingArtifactMetrics工具中产生的pre_adapter_detail_metrics.

(1)用picard的CollectSequencingArtifactsMetric统计序列背景的artifacts.此工具可以区分杂交筛选（hybrid selection[2]）前的artifacts和杂交筛选时的诱饵bias.这一步工作可以从基因组全局的角度来决策，这是基于位点的分析做不到的。

gatk CollectSequencingArtifactMetrics \
-I tumor.bam \
-O 10_tumor_artifact \
–-FILE_EXTENSION ".txt" \
-R ~/Documents/ref/hg38/Homo_sapiens_assembly38.fasta

生成结果包括5个metrics(性能指标)文件，

Picard Metrics Definitions
​broadinstitute.github.io/picard/picard-metric-definitions.html

(2)用FilterByOrientationBias过滤方向bias

输入：FilterMutectCalls的输出VCF，(1)步生成文件pre_adapter_detail_metrics和指定的FFPE(C->T)、OxoG(G->T)序列背景，工具默认已经包含了反向互补链的情况

gatk FilterByOrientationBias \
-A G/T \
-A C/T \
-V 9_somatic_oncefiltered.vcf.gz \
-P tumor_artifact.pre_adapter_detail_metrics.txt \
-O 11_somatic_twicefiltered.vcf.gz

输出文件是一个VCF文件和相关summary文件：

11_somatic_twicefiltered.vcf.gz.summary

2.6 展示如何建立IGV以人工查看体细胞突变事件

本地window10安装IGV软件

输入6个文件，用IGV查看变异相关区域

a. PoN.vcf

b. Chr17_af-only-gnomad_grch38.vcf.gz

c. 11_somatic_twicefilted.vcf.gz

d. 2_tumor_normal_m2.bam

e. Normal.bam

f. Tumor.bam

用法从略

三、报错

3.1 报错

A USER ERROR has occurred: Bad input: Sample sample20170711-4_I6_S6_L007 is not in BAM header: [Msample20190927-1-I44_L4_D95T59, sample20170711-4_I6_S6_L007_R1]

解决办法：

需要通过AddOrReplaceReadGroups 赋值给RGSM， 然后在 Mutect2的时候，指定tumor和Normal

sample_name=SRR23455317;

#14 min
time singularity exec /data2/neoantigen/singularity/nextNEOpi_1.3.2_18734d43.sif gatk AddOrReplaceReadGroups INPUT=work/"$sample_name".sam OUTPUT=work/"$sample_name".named.sam  RGID=1 RGLB=lib1 RGPL=illumina RGPU=unit1 RGSM=1

sample_name=SRR23455318;
time singularity exec /data2/neoantigen/singularity/nextNEOpi_1.3.2_18734d43.sif gatk AddOrReplaceReadGroups INPUT=work/"$sample_name".sam OUTPUT=work/"$sample_name".named.sam  RGID=1 RGLB=lib1 RGPL=illumina RGPU=unit1 RGSM=2


#25min 
time singularity exec /data2/neoantigen/singularity/nextNEOpi_1.3.2_18734d43.sif gatk SortSam INPUT=work/"$sample_name".named.sam OUTPUT=work/"$sample_name".named.sorted.bam SORT_ORDER=coordinate

#12min
time singularity exec /data2/neoantigen/singularity/nextNEOpi_1.3.2_18734d43.sif gatk  MarkDuplicates INPUT=work/"$sample_name".named.sorted.bam OUTPUT=work/"$sample_name".named.sorted.rmDup.bam METRICS_FILE=work/"$sample_name".Input.tmp_file REMOVE_DUPLICATES=true ASSUME_SORTED=true VALIDATION_STRINGENCY=SILENT TMP_DIR=tmp/


time samtools index work/"$sample_name".named.sorted.rmDup.bam


tumorsample=SRR23455317;
normalsample=SRR23455318;
time singularity exec /data2/neoantigen/singularity/nextNEOpi_1.3.2_18734d43.sif gatk --java-options "-Xmx250g" Mutect2   -R /data4/aNEO_data/mice/Mus_musculus/UCSC/mm10/Sequence/Chromosomes/ref_mm10.fasta  -I work/"$tumorsample".named.sorted.rmDup.bam -tumor 1  -I  work/"$normalsample".named.sorted.rmDup.bam -normal 2  -O work/gatk4.vcf.gz

参考资料

• https://zhuanlan.zhihu.com/p/63899960
• https://www.jianshu.com/p/716ebbd6c97f