四种肿瘤体细胞单核苷酸突变检测方法的比较

李晓东，何小雨，陈玮，李瑞琳，赵丹，祝海栋，张裕，代闯闯，陆忠华，迟学斌,3，牛北方,4*，郎显宇*

1. 中国科学院计算机网络信息中心，北京 100190

2. 中国科学院大学，北京 100049

3. 中国科学院计算科学应用中心，北京 100190

4. 贵州大学医学院，贵州 贵阳 550025

引言

在对肿瘤长久以来的研究中，研究人员逐渐将它定义为一种“基因组疾病”[1]。肿瘤起源于某个上皮细胞，当它积累了足够的体细胞突变后，会导致其自身的增殖失调，从而扩展成一个大的细胞群，细胞群中的子代细胞也保留有与祖先细胞一致的突变。能够诱发细胞发生肿瘤变异以及促进肿瘤发展的突变可以分为两种：抑制基因 (suppressor genes) 的失活性突变和癌基因 (oncogenes) 的激活性突变。两种变异都可能导致携带它们的细胞获得显著的成长优势，从而迅速发展成肿瘤组织。

肿瘤基因组研究期望寻找到变异基因与临床表现之间的联系，从而识别出驱动肿瘤进化发展的变异，为病人提供有效的精准治疗方案。研究发现[2]，肿瘤基因组变异主要有单核苷酸变异(Single Nucleotide Variant，SNV)，插入和缺失变异(Insertion and Deletion，Indel)、结构变异 (Structure Variant，SV) 三种类型。其中肿瘤体细胞 SNV 是一种最简单的变异类型，在肿瘤的发生中扮演重要的角色[3]。SNV 可能导致对应的氨基酸改变，甚至所编码蛋白质密码子变成终止密码子，从而影响基因的功能。通过高通量测序技术，研究人员能够直接获得肿瘤和正常组织样本的原始基因序列数据。通过分析样本基因序列，并与标准参考基因组比对，寻找与肿瘤相关的 SNV，是目前肿瘤基因研究中常用的方法。

在肿瘤序列分析中，体细胞 SNV 检测是后续分析的基础。寻找到变异位点集合后，通过比对已有数据库 (如 dbSNP，COSMIC) 中的变异位点，从而寻找与疾病相关的变异基因。因此，最终分析结果的可靠性依赖于寻找到的初始变异位点集合的质量。

然而，样品降解、遗传异质性和组织污染 (杂质)等因素影响，对 SNV 的检测存在较大挑战。同时，测序过程中引入的测序误差、序列重排引入的比对误差以及碱基覆盖度不足等，也给变异位点的检测增加了困难。早期的一些检测软件如 NaiveSubtract[4]，采用统计学显著性检验方法，分别对配对的肿瘤样本和正常样本进行分析，获得样本的 SNV 候选位点集，并从肿瘤样本的候选位点集中去除在正常样本中检测到的候选位点。但是，这种简单的“减法”没有对两个样本中共有的家系多态性位点进行分析，同时也不适用于发现低等位基因频率的变异。目前，针对体细胞突变检测的算法，均是通过联合分析肿瘤和正常样本序列，从而最大程度减小因测序平台及比对算法引入的随机或系统误差。

在对变异位点检测时，研究人员希望检测到的候选位点集具有较高的准确性，同时假阳性位点数目尽可能低。由于不同方法使用的模型和评估指标不同，不同软件识别的变异位点集合存在较大的不一致性，方法准确性也存在较大差别。本文针对目前最受欢迎的四种 SNV 检测算法 Varscan2[5]、SomaticSniper[6]、Strelka[7]和 MuTect2[8]进行测试，在已知基准 SNV突变信息的数据集上，比较四种方法的表现，同时对不同方法在不同测序深度上的结果进行分析。

1 体细胞 SNV 检测分析

DNA 序列包含腺嘌呤核苷酸 (A)、鸟嘌呤核苷酸(G)、胞嘧啶核苷酸 (C)、胸腺嘧啶核苷酸 (T) 四种核苷酸。SNV 主要发生在 DNA 复制过程中，即某种核苷酸突变成为另外一种类型核苷酸，从而导致同一个位点上的基因存在两个以上的状态，称作等位基因。在二倍体基因组中，对每一个位点，相较于参考基因组，设所有可能的等位基因型组合为 Ω={AA，AB，BB} A 表示该位点基因型与参考基因组一致，B 表示该位点基因型与参考基因组不一致。AA 为纯合子参考型，AB 为杂合子突变型，BB 为纯合子突变型(在人的基因组中，理论上存在 10 种组合，即 {AA，AC，AT，AG，CC，CT，CG，TT，TG，GG}，此处简化，每个位点只包含两种可能的基因型 )。对每一个位点，肿瘤样本和正常样本中联合基因型分布如表 1所示。

1.1 变异检测算法

本文分析的四种变异检测算法，主要使用经过质控、比对、去重、重校验后的 BAM 类型文件作为输入数据。BAM 文件是 SAM 文件的二进制格式，文件中包含每个 read 片段的原始序列，碱基质量，比对到参考序列的坐标信息，比对质量等信息。算法通过联合分析参考基因组、肿瘤和正常组织样本基因组，计算各位点的相关特征信息，采用不同模型，确定每个位点在不同样本中基因型，从而确定在肿瘤样本中该位点是否为体细胞单核苷酸变异位点 (图 1)。

1.1.1 Mutect2

Mutect2 算法基于隐马尔可夫模型的贝叶斯方法，设计了两个贝叶斯分类器。对变异位点的检测分为四步：(1) 去除低质量序列数据；(2) 使用第一个贝叶斯分类器对肿瘤样本进行变异检测；(3) 过滤由于相关测序过程中产生的假阳性误差，这类误差很难被误差模型捕获；(4) 使用第二个贝叶斯分类器，对识别的变异位点进行区分，判断为体细胞变异 (SNV) 或家系变异 (Germline variant)。

1.1.2 SomaticSniper

SomaticSniper 算法使用贝叶斯概率模型，基于参考序列，计算每种基因型的先验概率分布。通过每个位点在肿瘤和正常样本中的观测数据，计算位点联合基因型的后验概率，采用对数运算转换为每个位点的体细胞变异得分。转换后的分值在 0-255 之间，表示位点的变异程度，用户可指定相应的体细胞变异阈值，以过滤变异分数较低的候选位点。

图1 假设的肿瘤-正常样本比对示意图，等位基因数(Allelic Counts)， 表示该位点与参考基因一致的 read 数，表示覆盖该位点的总 read 数Fig. 1 The comparison diagram of hypothesized tumor-normal sample

1.1.3 Strelka

Strelka 首先搜索正常和肿瘤样本 BAM 文件中发生插入缺失的序列，对其重排。算法使用重排后的序列数据，基于贝叶斯概率模型，计算出候选位点在正常样本中最可能的基因型和正常基因型状态下体细胞变异概率。模型中使用到的序列信息包括碱基质量和比对质量，链偏差，位点为变异位点的先验概率以及正常样本中杂合型的预期率。由于测序的正常样本中一般包含有家系变异，肿瘤样本包含真正体细胞变异的肿瘤组织和正常组织，因此 Strelka 使用等位基因频率而不是二倍体基因型计算体细胞变异概率。对寻找到的原始 SNV 和 Indel 变异位点，Strelka 依据配置参数，以剔除假阳性位点。

1.1.4 Varscan2

Varscan2 采用启发式和统计学算法，分别生成包含 SNV 变异和 Indel 变异的 VCF 文件。相较于其他算法，Varscan2 使用经 samtools 程序基于 BAM 文件生成的 mpileup 类型文件作为输入[9]。在每一个位点，程序同时读取肿瘤和正常样本中的数据后，对位点的深度进行标准化，然后进行启发式配对比较。当一个位点存在某种碱基出现频率超过阈值，该位点基因型为纯合型，否则为杂合型。对在肿瘤和正常组织样本中基因型不匹配的位点，统计两个样本中参考支持型 (reference-supporting) 和变异支持型 (variantsupporting) 的序列数，进行单尾 Fisher 精确检验。如果检验的 P 值小于设定的显著性阈值，算法依据正常样本的基因型来确定该位点的变异类型，当正常样本为纯合型，分类为体细胞变异 (SNV)，当正常样本为杂合性，分类为杂合性缺失 (LOH)。如果大于显著性阈值，则为家系变异 (Germline variant)。在生成的VCF 文件中，标注每个位点的最大概率变异类型。

1.2 数据特征说明

上述四种变异检测算法均是通过统计待分析位置的数据信息，选择可用于分析的有效特征。不同的变异检测算法使用的数据特征各有不同，根据数据特征的特点，可分为以下六种类别，如表 2 所示。

(1) 序列测序深度：正常样本和肿瘤样本中覆盖某一位点的短 read 数；

(2) 碱基质量：在测序时，每测一个碱基会给出一个相应的质量值，用于衡量测序准确度，通常使用的 Phred 碱基质量值公式为，Q-score=-10×log10P，P为预估的碱基错误检查率；

(3) 链偏差：采用双端测序时，由于二代测序的读长限制，某些位置的碱基可能只在正向链或反向链上被测到；

(4) 序列比对质量：将样本序列比对到参考基因组时，由比对软件生成，用于衡量比对准确度；

(5) 位点位置值：位点位于所在序列的位置信息。位于所在序列中间时值为 1，位于所在序列两端值为0。越接近序列两端，由测序和比对过程产生误差的概率越大；

(6) 先验概率与检验水平：在对肿瘤和样本序列的比较中，检测方法采用统计假设检验，在用户指定的显著水平下，判断检测位点的基因型是否匹配或位点是否为体细胞变异位点。

1.3 候选 SNV 位点过滤

控制假阳性率是体细胞 SNV 变异检测方法面临的主要难点。理想状态下，变异识别算法应具有较高的敏感性，在最大程度发现真实变异的同时，保证检测到的变异为假阳性变异的概率最低。上述四种检测算法使用了多种信息特征来计算每个候选位点的体细胞变异概率。然而，由于 DNA 序列数据中的复杂因素影响，每个算法使用的特征因素并不完全相同，初始得到的候选变异位点集合中存在大量的假阳性变异位点。因此需要对候选位点进行过滤，以减少假阳性。

Varscan2 和 SomaticSniper 需要结合后续的流程，进行过滤。对初始得到的候选位点，提取位点的坐标信息，采用 bamreadcount[10]计算每个单核苷酸位点的详细信息，剔除不满足指定条件阈值的候选位点。Strelka 和 Mutect2 算法内置过滤流程。用户在算法运行前，需确定相关阈值参数，算法分别输出过滤前的变异位点集和过滤后的高可信度变异位点集。

2 实验数据

对于突变检测类方法性能的衡量，困难在于，没有一个已知完全变异信息的标准数据集。因此，无法对变异识别方法准确率和召回率进行准确评估。为解决这个问题，研究人员尝试使用两种途径构造可以作为衡量标准的数据集 (表 3)。

表2 四种变异检测算法使用的特征明细Table 2 The detailed characteristics of the four detection methods

途径 1：使用多测序平台对数据进行多次测序。对同一个样本，使用不同建库方法和不同测序平台进行测序。通过整合多个测序结果，得到高置信度的原始序列文件。美国国家标准和技术研究院 (NIST)，对人类基因组计划中的 NA12787 (NIST RM 8398) 样本，采用多种测序平台和建库方法，对样本进行测序，获得不同种类的数据[11]。通过不同数据间相互验证、补充，从而获得该样本的各类型变异信息。该数据通常可作为家系变异 (Germline variant)，插入缺失变异 (Indel) 和结构变异 (SV) 检测及全基因组组装的标准数据。

途径 2：使用序列数据模型软件，人工合成相应的变异数据。目前采用的对变异数据的模拟有两种方法：(1) 基于参考基因组，模拟读长和变异位点配置信息，合成突变基因组；(2) 基于已有的经过比对的正常序列数据，在其中一些混入预先设定的变异，设定的变异包括随机变异和已发现的相关疾病变异，以模拟肿瘤组织数据。通过上述方法，研究人员能够控制样本的准确变异信息，以对不同算法的性能进行有效评估。由于第一种方法无法模拟真实数据中变异位点的非随机分布，非独立误差和拷贝数变异等特性，目前分析上常采用第二种方法。

本文使用的测试数据来自由国际癌症基因组协会 (ICGC) 和癌症基因组图谱 (TCGA) 癌症基因组计划联合举办的变异识别挑战项目[12]。正常组织的基因序列数据是通过 Illumina HiSeq 2000 测序仪产生的双端 (paired-end) 测序数据，平均读长约为 101bp，平均覆盖度约为 60X。使用 GRCh37/hg19 基因组作为参考基因组，采用 BWA[13]软件进行序列比对，通过bamSurgeons[14]软件模拟产生对应的肿瘤基因序列数据。数据包含确定的 SNV 变异信息，可对每个检测算法的结果进行有效评估。

表3 实验数据描述Table 3 Experimental data description

表4 不同变异识别方法识别的相同候选位点比例Table 4 The proportion of the same candidate sites identi fi ed by different methods

3 实验结果及分析

在对上述四种方法进行测试过程中，Mutect2 和strelka 使用算法默认的参数。SomaticSniper 设定碱基质量参数，对初始识别的候选位点，采用程序开发者推荐的过滤方法进行过滤，获得高置信度变异位点集。Varscan2 使用基于算法默认参数，并依据数据集样本纯度，对相应参数进行调整，获得最佳的高置信度候选位点集。详细算法运行参数参见补充材料 1。

3.1 不同方法检测的变异位点数量存在显著差异

在三组正常-肿瘤全基因组样本对上的检测结果显示，四种检测算法识别的体细胞突变数量存在显著差异。其中，strelka 检测到的候选位点数目最少，Mutect2 和 SomaticSniper 检测到的候选位点数目较多。算法检测变异位点数量与样本真实变异位点数量呈正相关 (表 4)。

为确定不同算法检测结果的一致性，我们对每个样本的四个候选集交叉比较，计算同时被两种算法检测的位点比例。如图 2 所示，四种算法中，Mutect2、SomaticSniper 与其他算法的一致性较低，Strelka 的一致性最高。Strelka 检测到的位点，多数也被其他方法检测到。但同时，其他方法检测到的位点也被 Strelka 检测到的整体比例最低。由此说明，四种算法中 Strelka 的检测标准最为严苛，检测的候选位点数量也最少。

3.2 四种检测算法的性能评估

通过与各样本的真实变异位点比较，我们对各个算法的准确率、召回率进行衡量， 比较检测算法的综合性能。四种检测方法中，Varscan2 和 Strelka 的综合性能优于另外两种算法，SomaticSniper 的性能最低。Mutect2 损失一定的准确度，但能够检测到最多的真实变异。受样本纯度的影响，识别算法在 IS2 数据集的整体准确率表现相对数据集 IS1 和 IS3 更低，此时 Strelka 的准确率显著高于其他几种算法 (表 5)。

图2 四种变异检测方法检测到的候选位对比Fig. 2 The comparison of candidate positions number detected by four different methods

表5 四种检测算法在各个样本集上的检测性能评估Table 5 Four detection algorithms are evaluated on each sample set

图3 四种检测算法的交集位点数目和准确率Fig. 3 The intersection number and accuracy of the four detection algorithms

图4 四种检测算法的并集位点数目和召回率Fig. 4 The recall rate and convergence number of the four algorithm

图5 三组样本中真实变异位点的测序深度分布，横轴表示位点的测序深度Fig. 5 Sequencing depth distribution of real variation sites in three groups of samples

图6 在三组样本上四种检测方法的结果对比，虚线表示当前样本低测序深度部分真实变异数，蓝色折线表示检测算法检测到的低测序深度个数，黄色柱图表示各算法检测结果与真实变异一致的位点数Fig. 6 The detection effect of four detection methods at low-sequencing depth site

对四种检测方法得到的候选位点进行统计发现，被所有方法均检测到的位点数分别为 3563，4870，5768，交集中真实变异位点的准确率显著高于单个方法 (图 3)。四种方法的并集数分别为 8289，11511，11921，并集的召回率显著高于任意单一方法 (图 4)。

表6 四种检测算法对低测序深度变异位点的检测性能比较Table 6 The performance of different algorithms for low sequencing depth

3.3 测序深度对检测算法准确率的影响

在各个算法模型中，序列测序深度是所有检测算法均衡量的一个因素。提取三组样本中真实位点分别在对应肿瘤和正常样本中的测序深度数，在显著水平0.05 条件下，采用 Fisher 检验，真实变异位点在肿瘤和正常组织中的测序深度分布无显著差异。

在三组基准变异数据集中，我们选取在肿瘤和正常样本中测序深度均小于平均测序深度 50% 的真实变异位点，作为低测序深度位点，其个数分别为 119，120，400，占总变异位点的比例分别为3.33%，2.77%，5.06%。分析四种检测算法在低测序深度部分的检测结果，与真实变异位点比对发现，SomaticSniper 检测到的总位点数目最多，但也引入了最多的假阳性位点，Strelka 检测到的数目最少。Mutect2 检测到的真实低测序深度变异位点数目最多，同时引入的假阳性位点也较少 (图 5-6)。

对低测序深度部分检测算法的性能评估，各算法的准确率和召回率远小于其平均水平。分析结果表明，当测序深度较低时，变异频率较低的等位基因很难被识别到，从而对位点基因型的确定存在较大的干扰，影响算法的准确性。四种算法中，Mutect2 对低深度部分位点的检测具有最好的准确率和召回率。当样本数据的测序深度较低时，使用 Mutect2 会获得较好的检测结果 (表 6)。

4 结束语

通过对肿瘤-正常样本基因序列分析，发现准确的体细胞点突变位点对肿瘤的研究和癌症的个性化治疗至关重要。在本文中，我们验证了四种变异检测方法在三组模拟的正常-肿瘤全基因组序列上的突变识别效果，从不同维度对检测方法的性能进行了评估。尽管四种方法均检测到大量相同的 SNV 位点，尤其是真实的变异位点，但通过对三组基准数据的总体分析表明，VarScan2 在检测高质量 SNV方面表现最为出色，在准确率和召回率上较其他方法更为均衡；Strelka 对变异位点的识别标准最为严苛，故能保证较高的准确性，但会遗漏较多真实的变异位点；Mutect2 能够检测到最多的真实变异，但同时也会引入较多的假阳性位点，这对后续变异位点的验证工作造成较大的干扰；SomaticSniper 与Mutect2 特点相似，但准确性上表现更低。当样本纯度较低时，Strelka 具有更好的表现。此外，样本数据的测序深度会显著影响各算法的检测效果。针对低测序深度样本，我们更推荐使用 Mutect2 进行变异位点检测。研究人员可根据研究目的，结合各算法的不同特点，采用多种方法的组合，能够获得更多或更准确的候选位点集。

[1]Strausberg R L, Simpson A J G, Old L J, et al. Oncogenomics and the development of new cancer therapies[J].Nature,2004, 429(6990): 469-474

[2]Strausberg R L, Simpson A J G.Whole-genome cancer analysis as an approach to deeper understanding of tumour biology[J]. British journal of cancer, 2010,102(2):243-248

[3]Carlson B.SNPs-A shortcut to personalized medicine[J].Genetic Engineering & Biotechnology News, 2008,28(12): 12-12.

[4]DePristo M A, Banks E, Poplin R, et al.A framework for variation discovery and genotyping using next-generation DNA sequencing data[J].Nature genetics, 2011, 43(5):491-498.

[5]Koboldt D C, Zhang Q, Larson D E, et al.VarScan 2:somatic mutation and copy number alteration discovery in cancer by exome sequencing[J].Genome research, 2012,22(3): 568-576.

[6]Saunders C T, Wong W S W, Swamy S, et al.Strelka:accurate somatic small-variant calling from sequenced tumor–normal sample pairs[J].Bioinformatics, 2012,28(14): 1811-1817.

[7]Cibulskis K, Lawrence M S, Carter S L, et al. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples[J].Nature biotechnology,2013, 31(3): 213-219.

[8]Larson D E, Harris C C, Chen K, et al.SomaticSniper:identification of somatic point mutations in whole genome sequencing data[J].Bioinformatics, 2011, 28(3):311-317.

[9]Li H.A statistical framework for SNP calling,mutation discovery,association mapping and population genetical parameter estimation from sequencing data[J].Bioinformatics,2011, 27(21): 2987-2993.

[10]Li H.Toward better understanding of artifacts in variant calling from high-coverage samples[J].Bioinformatics,2014, 30(20): 2843-2851.

[11]Zook J M, Catoe D, McDaniel J, et al. Extensive sequencing of seven human genomes to characterize benchmark reference materials[J].Scientific Data, 2016,3(25).

[12]Ewing A D, Houlahan K E, Hu Y, et al.Combining tumor genome simulation with crowdsourcing to benchmark somatic single-nucleotide-variant detection[J].Nature methods, 2015, 12(7): 623-630.

[13]Li H, Durbin R.Fast and accurate long-read alignment with Burrows–Wheeler transform[J].Bioinformatics, 2010,26(5): 589-595.

[14]Xu H, DiCarlo J, Satya R V, et al.Comparison of somatic mutation calling methods in amplicon and whole exome sequence data[J].BMC genomics, 2014,15(1): 244-253.

[15]Abeel T, Helleputte T, Van de Peer Y, et al.Robust biomarker identification for cancer diagnosis with ensemble feature selection methods[J].Bioinformatics,2009, 26(3): 392-398.

[16]Alioto T S, Buchhalter I, Derdak S, et al.A comprehensive assessment of somatic mutation detection in cancer using whole-genome sequencing[J].Nature communications,2015, 6.

[17]Guo Y, Li J, Li C I, et al. The effect of strand bias in Illumina short-read sequencing data[J].BMC genomics,2012, 13(1): 666-676.