【2.6.3】使用GATK4 Mutect2 检测体细胞突变
二、常规流程
2.1 Mutect2检测体细胞突变
输入文件：
- 肿瘤样本名：HCC1143_tumor
- 配对癌旁样本名：HCC1143_normal
- PoN正常样本集VCF文件：后面给出生成PoN文件的方法
- germline资源：germline变异数据库VCF
- 目标区域文件：interval
输出文件：
- 原始的体细胞突变的VCF文件：1_somatic_m2.vcf.gz(.tbi)
- 重组装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都将考虑在内。
- 不过滤比对到不同contig上的PE reads，大概是因为BWA的PE比对策略（alt-aware mapping,unclear），如果滤掉不配对reads,会丢失一部分与备用单体型相关的SNV和Indel信息。
- -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参数
- 注意：样本交叉污染不同于tumor的normal污染和normal样本的tumor污染。目前，这个工具还不能用来计算算肿瘤的纯度。
- 如果发现了高度污染，那么第一件事情就是核查样本是否在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
参考资料
