CN114694750B - 一种基于ngs平台的单样本肿瘤体细胞突变判别及tmb检测方法 - Google Patents

Abstract

本申请涉及生信分析技术领域，具体提供一种基于NGS平台的单样本肿瘤体细胞突变判别方法及相应TMB检测方法，本方法基于NGS平台数据得到的变异检测结果，通过公共数据库和自建国人突变数据库比对，并基于单样本计算的等位基因拷贝数构建统计分布模型的方法,进行单样本肿瘤体细胞突变的判别，以此实现TMB计算。

Description

一种基于NGS平台的单样本肿瘤体细胞突变判别及TMB检测 方法

技术领域

本申请属于生信分析领域，具体涉及一种基于NGS平台的单样本肿瘤体细胞突变判别方法及其相应的TMB检测方法。

背景技术

肿瘤样本的下一代测序（Next Generation Sequencing, NGS）被广泛用于发现生物学上重要的突变，并指导临床进行靶向治疗和用药。这些突变分为体细胞突变（Somaticmutation）和种系突变（Germline mutation），体细胞突变一般仅存在于肿瘤细胞中，不遗传给后代，而种系突变发生在受精卵中，同时存在于肿瘤细胞和正常细胞中，可以遗传给后代。细胞的癌变和肿瘤的发生发展通常与体细胞中基因突变长期积累的结果有关，但不是所有的体细胞突变都会导致细胞癌变。通常每位患者都会检测到大量的突变结果，其中只有少量突变是潜在的肿瘤驱动突变（Driver mutation）或与临床治疗相关的突变，而其余的则是乘客突变（Passenger mutation）或在功能上是良性的种系突变。这些驱动突变通过细胞信号通路和一系列复杂的生物学作用又会引发其他的基因突变，这导致肿瘤细胞中的基因突变数量往往高于正常的体细胞。肿瘤突变负荷(Tumor Mutational Burden，TMB)则是衡量肿瘤细胞中总的基因突变程度的指标，通常以每百万碱基(Mb)的基因外显子编码区包含的肿瘤体细胞单核苷酸突变（SNV）以及小片段插入缺失变异(INDEL)总数来表示。

近年来，免疫治疗在皮肤癌、膀胱癌、肺癌和肾癌以及错配修复缺陷的癌症患者中表现出令人欣喜的临床效果。免疫治疗主要控制细胞免疫反应的时间和位置，通过免疫检查点阻断增加抗肿瘤活性。细胞毒性T淋巴细胞相关抗原 4 (CTLA-4) 和程序性细胞死亡蛋白 1 (PD-1) 受体的抑制剂是目前研究广泛和临床发展最快的药物。虽然免疫治疗已在多种肿瘤治疗中获得显著的临床疗效，但这其中只有部分患者在这些治疗中获益。免疫组织化学 (IHC) 测量 PD-1/PD-L1 蛋白表达是最早被批准作为一些免疫治疗药物的伴随诊断方法，然而，PD-1/PD-L1 表达的测量在技术上具有一定的挑战性并且并不总是准确预测免疫治疗的效果。TMB作为一种新兴免疫治疗的生物标志物，在多个肿瘤类型中已被证明与患者对 CTLA-4 和 PD-1 抑制剂的反应显著相关，且TMB 与反应率的相关性比IHC检测中PD-L1 的表达量更显著。因此，TMB的准确检测对于预测免疫检查点抑制剂的疗效具有重要意义，并且可以使癌症患者有机会获得更加精准的治疗。

无论是肿瘤驱动突变或与临床治疗相关的突变检测还是TMB精确检测，准确的判别体细胞突变都至关重要。通常，区分体细胞突变与种系突变需要对肿瘤及其匹配的正常对照组织一起进行测序，在肿瘤组织中检测到但在正常对照中不存在的变异被判定为体细胞突变。然而，随着临床癌症测序变得越来越广泛，大多数病例可能无法获得匹配的正常对照样本，并且由于各种原因，组织学上正常的组织可能由于肿瘤细胞浸润以及处理过程中肿瘤细胞污染造成最终体细胞突变的判别存在偏差。因此，在没有正常对照样本的情况下准确的判别肿瘤体细胞突变，进而进行单样本TMB检测在实际临床应用中具有很大的必要性。

当只有肿瘤的测序数据可用时，研究人员开发了相关的分析过滤策略来区分种系和体细胞突变。目前，最典型的方法主要是通过匹配公共人群数据库（如 dbSNP，千人基因组，gnomAD，ExAC）的方法过滤掉存在于人群数据库中的种系突变。有研究显示每个个体中都存在大量的未在目前公共人群数据库中记录的私有种系突变，而且私有种系突变率与个体的血统显著相关，同时也取决于群体在大规模测序项目中的代表性。然而，目前公共人群数据库中对于中国人群抑或东亚人群的群体数据的记录非常有限，所以这些方法并不能充分的去除中国人群私有的种系突变。因此，一方面需要扩展目前中国人群数据库的种系突变数据，另一方面需要开发新的生物信息学方法在肿瘤单样本中高灵敏度和精确度判别种系和体细胞突变。

综上所述，肿瘤单样本中体细胞突变的判别和TMB的计算受到目前公共数据库的局限性以及仍然不够完善的生物信息学过滤方法的较大影响，判别过程存在较多的假阳现象。因此本申请基于NGS平台开发了一种单样本肿瘤体细胞突变判别和TMB检测方法，该方法保证了单样本肿瘤体细胞突变判别和TMB检测的准确性

有鉴于此，提出本申请。

发明内容

为解决上述现有技术问题，本申请基于NGS平台数据得到的变异检测结果，通过公共数据库和自建中国人群突变数据库的比对策略以及基于单样本计算的等位基因拷贝数构建统计分布模型的方法进行单样本肿瘤体细胞突变的判别，并以判别结果实现TMB的计算。

因此，本申请至少包括如下几方面目的：

本申请的第一目的为提供一种基于NGS平台的单样本肿瘤体细胞突变判别方法及模型，以及模型构建方法；

本申请的第二目的为提供一种基于上述单样本肿瘤体细胞突变判别的TMB检测方法及模型，以及模型构建方法。

为实现上述目的，本申请具体提供如下技术方案：

本申请首先提供一种基于NGS平台的单样本肿瘤体细胞突变判别模型的构建方法，所述方法包括如下步骤：

1）肿瘤相关基因捕获：单样本DNA提取后通过杂交捕获肿瘤相关基因区域；

2）基因测序：NSG对肿瘤相关基因进行测序；

3）数据比对：将NGS测序数据比对到人参考基因组；

4）突变检测：基于人参考基因组比对数据分析并获取SNV与INDEL，并计算得到突变的等位基因频率（AF）；

5）拷贝数区段构建：基于正常拷贝数基线和单样本Bin区域测序深度数据划分区段，并计算区段的拷贝数；

6）国人种系突变数据库构建：基于中国人群对照样本测序数据进行种系突变检测并计算每个种系突变的等位基因频率及其突变发生频率，完成种系突变数据库构建；

7）基因突变注释：使用公共数据库和6）国人种系突变数据库对4）获得的SNV和INDEL突变进行数据库注释，分别确定先验体细胞突变和先验种系突变及其等位基因频率AF；

8）拷贝数拟合模型构建：基于5）获取的区段的拷贝数及7）获取的先验种系突变和先验体细胞突变构建拷贝数拟合模型，获得肿瘤纯度、倍性以及矫正后的区段拷贝数和次等位基因拷贝数；

9）体细胞突变判别模型构建：基于8）获得的肿瘤纯度、倍性以及区段拷贝数和次等位基因拷贝数构建体细胞突变判别模型，获得体细胞突变判别结果；

10）体细胞突变过滤：基于7）基因突变注释的结果和9）获得的体细胞突变判别结果，过滤得到可靠的体细胞突变结果。

进一步的，所述步骤3）中，所述人参考基因组为hg19基因组。

进一步的，所述步骤5）中，所述拷贝数区段构建步骤如下：

5a、构建Bin区域：将杂交捕获目标区域Bed划分成Bin；

5b、基线构建：选用正常对照样本获得每个样本的每个Bin的测序深度，建立正常拷贝数基线；

5c、区段划分：基于基线数据和肿瘤组织单样本Bin区域测序深度数据，使用循环二元分割CBS算法将每条染色体臂上拷贝数CN相近的Bin区域合并成区段Segment，并取Bin区域CN的中位值作为Segment的CN。

进一步的，所述步骤6）中，所述中国人群种系突变数据库的具体构建步骤如下：

6a、基于中国人群正常对照样本的NGS测序数据，优选采用Sentieon种系突变检测工具，进行种系突变的检测；

6b、基于上述6a中获得的种系突变，过滤低覆盖深度和低支持数的位点，获得自建中国人群种系突变集；

6c、基于上述6b种系突变集计算每个种系突变的等位基因频率及其在该人群中突变发生频率，完成自建种系突变数据库构建。

进一步的，所述步骤7）中，所述突变注释的具体步骤包括：

7a、使用COSMIC数据库对所述4）获得的SNV和INDEL突变进行注释，获得COSMIC突变计数，通过COSMIC突变计数来确定先验体细胞突变；

7b、使用dbSNP、gnomAD、ExAC和1000G数据库对所述4）获得的SNV和INDEL突变进行注释，获得公共人群数据库中突变的人群频率，通过公共人群数据库突变的人群频率来确定先验种系突变；

7c、使用所述6）构建的中国人群突变数据库对所述4）获得的SNV和INDEL突变进一步注释得到中国人群特有的先验种系突变。

进一步的，所述步骤8）中，所述拷贝数拟合模型的具体构建步骤如下：

8a、基于所述5）获取的Segment的CN值以及上述7）获取的先验种系突变和先验体细胞突变及其突变频率AF，（优选使用DBSCAN算法）对CN和AF作二维聚类得到CN和AF均相近的多个Segment；

8b、通过Grid-Search拟合所有可能的拷贝数CN、肿瘤纯度Purity、肿瘤倍性Ploidy和次等位基因拷贝数nB组合，按照下列公式计算每种组合的理论AF与真实AF之间距离，即差值绝对值，通过拟合距离最小的组合获得样本的Purity和Ploidy；

CN=Ploidy×Purity+2×(1-Purity)

AF=(nB×Purity+1×(1-Purity))/CN

8c、使用上述8b中的方法重新计算并校正所有Segment的拷贝数CN和次等位基因拷贝数nB。

进一步的，所述步骤9）中，所述体细胞突变判别模型的具体构建步骤如下：

9a、基于上述8）获得的Purity以及各个Segment的CN和nB结果，通过下列公式计算Segment上每个突变位点的种系突变期望AF（ExpAFgermline）和体系突变期望AF（ExpAFsomatic）；

ExpAFgermline=(Purity×M+(1-Purity))/(Purity×CN+2×(1-Purity))

ExpAFsomatic=(Purity×M)/(Purity×CN+2×(1-Purity))

其中，M表示突变位点的等位基因拷贝数，M=nB或者M=CN-nB；

9b、基于上述9a中计算的ExpAFgermline和ExpAFsomatic以及4）中获得的突变信息构建下述二项分布统计模型，分别计算位点是种系突变的概率Probability(G)和体细胞突变的概率Probability(S)；

Probability(G)=Bin(n×f,n,ExpAFgermline)

Probability(S)=Bin(n×f,n,ExpAFsomatic)

其中，n表示突变位点的测序深度，f表示突变位点的等位基因频率；

9c、基于上述9b中计算的Probability(G)和Probability(S)进行体细胞突变的判别，若（Probability(S) > α且 Probability(G) ≤α）则该位点判定为体细胞突变；若（Probability(S) ≤α且 Probability(G) ≤α且 Purity ≥ 0.2）则该位点判定为体细胞亚克隆突变；若（Probability(S) ≤α且 Probability(G) > α）则该位点判定为种系突变；其中α为判定阈值。

本申请还提供一种肿瘤突变负荷TMB评估模型的构建方法，包括上述任一所述方法的步骤，并进一步包括如下步骤：

11）基于上述任一所述步骤10）获得的体细胞突变结果及杂交捕获的基因编码区域大小，计算肿瘤突变负荷TMB值；

TMB=(MC/(Length(bp)))×1000000

其中，MC表示纳入TMB计算的体细胞突变数目，Length(bp)表示杂交捕获的基因编码区域大小。

本申请还提供一种计算机可读存储介质，其上存储有计算机程序/指令，该计算机程序/指令被处理器执行时实现上述任一所述方法的步骤。

本申请还提供一种计算机程序产品，包括计算机程序/指令，该计算机程序/指令被处理器执行时实现上述任一所述方法的步骤。

与现有技术相比，本申请至少具有如下优势：

(1) 本申请开发了一种基于NGS平台的单样本肿瘤体细胞突变判别和TMB检测方法，在无配对样本情况下也能够准确地检测肿瘤样本的体细胞突变和TMB。该方法一方面节省了配对样本的实验、测序和分析时间及成本，另一方面可以在无配对样本检测情况下准确的为癌症患者临床决策提供信息并扩大治疗选择。

(2) 本发明基于中国人群NGS数据构建了种系突变数据库，有效补充了目前公共人群数据库中中国人群的种系突变信息，有助于准确的标记中国人群私有的种系突变，保证了单样本体细胞突变判别和TMB检测的准确性。

(3) 本发明通过拷贝数拟合模型构建了基于肿瘤纯度和等位基因拷贝数的单样本体细胞突变判别统计分布模型。该模型方法能够准确的标记公共数据库如dbSNP、gnomAD等未记录的罕见私有种系突变并进行体细胞主克隆突变和亚克隆突变的分类，去除假阳性体细胞突变。

附图说明

为了更清楚地说明本申请具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本申请的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1 单样本体细胞突变判别和TMB检测流程图；

图2 自建中国人群种系突变数据库流程图；

图3 单样本TMB检测与配对样本TMB检测结果的一致性；

图4 单样本TMB检测与WES TMB检测结果的一致性。

具体实施方式

下面将结合附图对本申请的技术方案进行清楚、完整地描述，显然，所描述的实施例是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。

以下术语或定义仅仅是为了帮助理解本申请而提供。这些定义不应被理解为具有小于本领域技术人员所理解的范围。

除非在下文中另有定义，本申请具体实施方式中所用的所有技术术语和科学术语的含义意图与本领域技术人员通常所理解的相同。虽然相信以下术语对于本领域技术人员很好理解，但仍然阐述以下定义以更好地解释本申请。

如本申请中所使用，术语“包括”、“包含”、“具有”、“含有”或“涉及”为包含性的(inclusive)或开放式的，且不排除其它未列举的元素或方法步骤。术语“由…组成”被认为是术语“包含”的优选实施方案。如果在下文中某一组被定义为包含至少一定数目的实施方案，这也应被理解为揭示了一个优选地仅由这些实施方案组成的组。

在提及单数形式名词时使用的不定冠词或定冠词例如“一个”或“一种”，“所述”，包括该名词的复数形式。

本申请中的术语“大约”、“大体”表示本领域技术人员能够理解的仍可保证论及特征的技术效果的准确度区间。该术语通常表示偏离指示数值的±10％，优选±5％。

此外，说明书和权利要求书中的术语第一、第二、第三、(a)、(b)、(c)以及诸如此类，是用于区分相似的元素，不是描述顺序或时间次序必须的。应理解，如此应用的术语在适当的环境下可互换，并且本申请描述的实施方案能以不同于本申请描述或举例说明的其它顺序实施。

下面为具体的实施例。

实验例、本申请方法体系建立

如图1所示，本申请的单样本体细胞突变判别和TMB计算包括以下步骤：

1、从肿瘤组织单样本提取DNA并通过杂交捕获（比如使用人肿瘤多基因检测Panel）肿瘤相关基因区域捕获肿瘤相关基因；

2、通过NSG平台对捕获的肿瘤相关基因进行测序，获得高通量测序原始数据并进行质控，去除低质量测序序列；

3、将质控后的测序数据比对到人参考基因组hg19上，并去除重复序列；

4、突变检测：基于人参考基因组比对数据分析并获取单核苷酸位点变异（SNV）与小片段插入缺失变异（INDEL），并计算得到突变的等位基因频率（AF）；

5、拷贝数区段构建：

5a、构建Bin区域：将多基因检测Panel目标区域（Bed）划分成Bin；

5b、基线构建：选用正常对照样本获得每个样本的每个Bin的测序深度，建立正常拷贝数基线；

5c、区段划分：基于基线数据和肿瘤组织单样本Bin区域测序深度数据，使用循环二元分割（Circular Binary Segmentation，CBS）算法将每条染色体臂上拷贝数（CN）相近的Bin区域合并成区段（Segment），并取Bin区域CN的中位值作为Segment的CN；

6、自建中国人群种系突变数据库：按照图 1 流程进行中国人群种系突变数据库构建：

6a、基于3722例正常对照样本的NGS测序数据，采用Sentieon种系突变检测工具进行种系突变的检测；

6b、基于上述6a中获得的种系突变，过滤掉低覆盖深度和低支持数的位点，获得自建中国人群种系突变集；

6c、基于上述6b种系突变集计算每个种系突变的等位基因频率及其在该人群中突变发生频率，完成自建种系突变数据库构建。

7）突变注释：

7a、上述4）获得的SNV和INDEL突变使用COSMIC数据库进行注释，获得COSMIC突变计数，通过COSMIC突变计数来确定先验体细胞突变；

7b、上述4）获得的SNV和INDEL突变使用dbSNP、gnomAD、ExAC和1000G数据库进行注释，获得公共人群数据库中突变的人群频率，通过公共人群数据库突变的人群频率来确定先验种系突变；

7c、上述4）获得的SNV和INDEL突变使用上述6）构建的中国人群突变数据库进一步注释过滤得到中国人群特有的先验种系突变；

8）拷贝数拟合模型构建：

8a、基于上述5）获取的Segment的CN值以及上述7）获取的先验种系突变和先验体细胞突变及其突变等位基因频率AF，使用DBSCAN算法对CN和AF作二维聚类得到CN和AF均相近的多个Segment；

8b、通过Grid-Search拟合所有可能的拷贝数（CN）、肿瘤纯度（

）、肿瘤倍性（Ploidy）和次等位基因拷贝数（nB）组合，按照下列公式计算每种组合的理论AF与真实AF之间距离（差值绝对值），通过拟合距离最小的组合获得样本的Purity和Ploidy。

CN=Ploidy×Purity+2×(1-Purity)

AF=(nB×Purity+1×(1-Purity))/CN

8c、使用上述8b中的方法重新计算并校正所有Segment的拷贝数（CN）和次等位基因拷贝数（nB）。

9）体细胞突变判别模型：

9a、基于上述8）获得的Purity以及各个Segment的CN和nB结果，通过下列公式计算Segment上每个突变位点的种系突变期望AF（ExpAFgermline）和体系突变期望AF（ExpAFsomatic），其中M表示突变位点的等位基因拷贝数，M=nB或者M=CN-nB；

ExpAFgermline=(Purity×M+(1-Purity))/(Purity×CN+2×(1-Purity) )

ExpAFsomatic=(Purity×M)/(Purity×CN+2×(1-Purity) )

9b、基于上述9a中计算的ExpAFgermline和ExpAFsomatic以及4）中获得的突变信息构建下述二项分布统计模型，分别计算位点是种系突变的概率Probability(G)和体细胞突变的概率Probability(S)，其中n表示突变位点的测序深度，f表示突变位点的等位基因频率。

Probability(G)=Bin(n×f,n,ExpAFgermline)

Probability(S)=Bin(n×f,n,ExpAFsomatic)

9c、基于上述9b中计算的Probability(G)和Probability(S)进行体细胞突变的判别，若（Probability(S) > α且 Probability(G) ≤α）则该位点判定为体细胞突变；若（Probability(S) ≤α且 Probability(G) ≤α且 Purity ≥ 0.2）则该位点判定为体细胞亚克隆突变；若（Probability(S) ≤α且 Probability(G) > α）则该位点判定为种系突变；其中α为判定阈值。

10）体细胞突变过滤：基于上述7）的突变注释结果和9）体细胞突变判别模型，过滤得到可靠的体细胞突变结果。

11）TMB计算：基于上述10）获得的体细胞突变结果以及人肿瘤多基因检测Panel覆盖编码区域的大小，肿瘤突变负荷TMB值通过如下方式计算得到；

TMB=MC/(Length(bp))×1000000

其中，MC表示纳入TMB计算的体细胞突变数目，Length(bp)表示人肿瘤多基因检测Panel覆盖编码区域的大小。

实施例1

本实施例中阳性标准品具体体细胞突变信息如下表1所示，其中体细胞突变的频率主要分布在1%-2%。本实施例共20例HD789为gDNA标准品，模拟福尔马林固定和石蜡包埋（FFPE）组织样本，来检测本方法的检测性能。所有样本均使用人肿瘤多基因检测Panel靶向捕获的方式获得DNA测序数据，且靶向捕获范围覆盖标准品中的阳性位点，测序深度1500X。

表1 HD789标准品体细胞突变信息

20例标准品样本测序数据，经过测序数据质控后均使用BWA（v0.7.17）MEM算法进行人参考基因组比对分析，去除重复序列后得到的BAM格式文件作为输入文件，进行本检测方法的单样本体细胞突变检测和判别分析。

20例阳性标准品检出结果如表2所示，所有180个阳性体细胞突变位点，基于本检测方法均能正确的判别和检出。

表2 阳性标准品检测性能评估结果

注：TP=True Positive，表示真阳；FN=False Negative，表示假阴。

对比例1 与具有配对样本的体细胞突变检测比较。

本对比例中采用来源于泛癌种检测的180例临床回顾样本，每例临床样本均进行肿瘤组织和匹配的对照样本的人肿瘤多基因检测Panel靶向捕获测序，测序深度1500X。180例临床样本测序数据，经过测序数据质控后均使用BWA（v0.7.17）MEM算法进行人参考基因组比对分析，去除重复序列后得到的BAM格式文件作为输入文件。

180例临床样本均分别进行具有配对样本的体细胞突变检测和判别分析以及基于本检测方法的单样本体细胞突变检测和判别分析。其中，以配对样本体细胞突变检测结果为标准，主要定义如下：如果突变位点同时在匹配的正常样本中检出并具有显著差异的等位基因频率，则判定为种系突变，若只有肿瘤样本中检出则判定为为体细胞突变。为了进行本检测方法与标准的配对样本检测方法的性能比较，主要进行下述3个方面测试：

1）本检测方法的单样本体细胞突变检测和判别与标准的配对样本检测方法在已知肿瘤体细胞热点突变检测和判别中的性能比较；

2）本检测方法的单样本体细胞突变检测和判别与标准的配对样本检测方法在所有体细胞突变检测和判别中的性能比较；

3）本检测方法的单样本TMB检测与标准的配对样本TMB检测方法的一致性评估。

结果如表3所示，180例临床样本中，配对样本检测出的128个体细胞热点突变均能使用本检测方法正确的检测和判别，体细胞热点检出准确率100%。同时，针对180例临床样本中所有的体细胞突变的评估显示，以配对样本检出结果为标准，本检测方法的灵敏度（Sensitivity=TSC/(TSC+FGC)）达到90.64%，阳性预测值（PPV=TSC/(TSC+FSC))达到90.04%。

表3 184例泛癌临床样本单样本体细胞突变检测与配对样本模式下体细胞突变检测一致性

注：TSC=True Somatic Count，真阳体细胞突变数目；FSC=False Somatic Count，假阳体细胞突变数目；FGC=False Germline Count，假阳种系突变数目。

图3结果显示，其中Paired-TMB为配对样本模式下TMB检测结果，TumorOnly-TMB为本检测方法单样本TMB检测结果，180例泛癌临床样本Paired-TMB与TumorOnly-TMB的相关性R²达到0.9946。

该结果表明，对比例1基于本检测方法的体细胞突变判别分析与配对样本的体细胞突变检测结果一致性超过90%，尤其是对于体细胞热点突变可以达到100%的一致性；同时，对比例1基于本检测方法的TMB检测准确性非常高，单样本TMB检测结果与配对样本TMB检测结果相关性超过99%。因此，基于本检测方法能够准确地检测肿瘤样本的体细胞突变和TMB，在单样本的情况下同样可以为癌症患者临床决策提供信息并扩大治疗选择。

对比例2 与WES-TMB金标准的比较

本对比例中选取24例泛癌肿瘤样本，提取肿瘤样本的DNA并纯化，采用“人类多基因突变联合检测试剂盒（可逆末端终止测序法）（SimcereDx Onco500）”进行DNA片段化、末端修复、接头连接和文库扩增；使用杂交探针对文库进行杂交捕获，捕获产物进行扩增、定量和标准化后应用NextSeq550Dx测序仪进行PE100双端测序，获得肿瘤样本目标区域的测序数据，测序深度1500X。24例肿瘤样本测序数据，经过数据质控后使用BWA（v0.7.17）MEM算法进行人参考基因组比对分析，比对结果通过分子标签（UMI）进行碱基矫正和去重复后进行SNV/INDEL的突变检测；基于SNV/INDEL突变结果采用本检测方法进行单样本体细胞突变判别分析并计算TMB。同时，24例肿瘤样本进行具有配对样本的全外显子组测序（WES），WES数据通过数据质控，BWA（v0.7.15）MEM算法进行人参考基因组比对和VarDict（v1.5.7）进行配对样本模式的SNV/INDEL突变检测，保留检测范围内高置信度的体细胞突变进行TMB计算。其中，以配对样本全外显子组测序检测的TMB作为金标准，评估基于本检测方法的单样本TMB检测结果与临床公认的WES TMB检测结果的一致性。

结果如图4所示，其中WES-TMB为全外显子组（WES）的TMB检测结果，TumorOnly-TMB为本检测方法应用“人类多基因突变联合检测试剂盒（可逆末端终止测序法）（SimcereDxOnco500）”的单样本TMB检测结果，24例泛癌临床样本WES-TMB与TumorOnly-TMB的相关性R²达到0.9717。

以全外显子组（WES）的TMB检测结果作为临床公认的金标准，该结果表明，对比例2中将本检测方法应用于SimcereDx Onco500的单样本TMB检测结果与WES TMB的一致性达到97.17%。由此可知，基于本检测方法的单样本TMB检测与金标准WES TMB的检测结果一致性非常高。同时，本检测方法可以应用在“人类多基因突变联合检测试剂盒（可逆末端终止测序法）（SimcereDx Onco500）”进行单样本TMB检测，在无对照样本的情况下同样可以为癌症患者临床决策提供信息并扩大治疗选择。

前述对本申请的具体示例性实施方案的描述是为了说明和例证的目的。这些描述并非想将本申请限定为所公开的精确形式，并且很显然，根据上述教导，可以进行很多改变和变化。对示例性实施例进行选择和描述的目的在于解释本申请的特定原理及其实际应用，从而使得本领域的技术人员能够实现并利用本申请的各种不同的示例性实施方案以及各种不同的选择和改变。本申请的范围意在由权利要求书及其等同形式所限定。

Claims (6)

1.一种基于NGS平台的单样本肿瘤体细胞突变判别模型的构建方法，其特征在于，所述方法包括如下步骤：

1）肿瘤相关基因捕获：单样本DNA提取后通过杂交捕获肿瘤相关基因区域；

2）基因测序：NSG对肿瘤相关基因进行测序；

3）数据比对：将NGS测序数据比对到人参考基因组；

4）突变检测：基于人参考基因组比对数据分析并获取SNV与INDEL，并计算得到突变的等位基因频率AF；

5）拷贝数区段构建：基于正常拷贝数基线和单样本Bin区域测序深度数据划分区段，并计算区段的拷贝数；

6）国人种系突变数据库构建：基于中国人群对照样本测序数据进行种系突变检测并计算每个种系突变的等位基因频率及其突变发生频率，完成种系突变数据库构建；

7）基因突变注释：使用公共数据库和6）国人种系突变数据库对4）获得的SNV和INDEL突变进行数据库注释，分别确定先验体细胞突变和先验种系突变及其突变等位基因频率AF；

8）拷贝数拟合模型构建：基于5）获取的区段的拷贝数及7）获取的先验种系突变和先验体细胞突变构建拷贝数拟合模型，获得肿瘤纯度、倍性以及矫正后的区段拷贝数和次等位基因拷贝数；

9）体细胞突变判别模型构建：基于8）获得的肿瘤纯度、倍性以及区段拷贝数和次等位基因拷贝数构建体细胞突变判别模型，获得体细胞突变判别结果；

10）体细胞突变过滤：基于7）基因突变注释的结果和9）获得的体细胞突变判别结果，过滤得到可靠的体细胞突变结果；

所述步骤7）中，所述突变注释的具体步骤包括：

7a、使用COSMIC数据库对所述4）获得的SNV和INDEL突变进行注释，获得COSMIC突变计数，通过COSMIC突变计数来确定先验体细胞突变；

7b、使用dbSNP、gnomAD、ExAC和1000G数据库对所述4）获得的SNV和INDEL突变进行注释，获得公共人群数据库中突变的人群频率，通过公共人群数据库突变的人群频率来确定先验种系突变；

7c、使用所述6）构建的中国人群突变数据库对所述4）获得的SNV和INDEL突变进一步注释得到中国人群特有的先验种系突变；

所述步骤8）中，所述拷贝数拟合模型的具体构建步骤如下：

8a、基于所述5）获取的Segment的CN值以及所述7）获取的先验种系突变和先验体细胞突变及其突变等位基因频率AF，对CN和AF作二维聚类得到CN和AF均相近的多个Segment；

8b、通过Grid-Search拟合所有可能的拷贝数CN、肿瘤纯度Purity、肿瘤倍性Ploidy和次等位基因拷贝数nB组合，按照下列公式计算每种组合的理论AF与真实AF之间距离，即差值绝对值，通过拟合距离最小的组合获得样本的Purity和Ploidy；

CN=Ploidy×Purity+2×(1-Purity)

AF=(nB×Purity+1×(1-Purity))/CN

8c、使用上述8b中的方法重新计算并校正所有Segment的拷贝数CN和次等位基因拷贝数nB；

所述步骤9）中，所述体细胞突变判别模型的具体构建步骤如下：

9a、基于上述8）获得的Purity以及各个Segment的CN和nB结果，通过下列公式计算Segment上每个突变位点的种系突变期望的AF：ExpAFgermline，和体系突变期望的AF：ExpAFsomatic；

ExpAFgermline=(Purity×M+(1-Purity))/(Purity×CN+2×(1-Purity) )

ExpAFsomatic=(Purity×M)/(Purity×CN+2×(1-Purity) )

其中，M表示突变位点的等位基因拷贝数，M=nB或者M=CN-nB；

9b、基于上述9a中计算的ExpAFgermline和ExpAFsomatic以及4）中获得的突变信息构建下述二项分布统计模型，分别计算位点是种系突变的概率Probability(G)和体细胞突变的概率Probability(S)；

Probability(G)=Bin(n×f,n,ExpAFgermline)

Probability(S)=Bin(n×f,n,ExpAFsomatic)

其中，n表示突变位点的测序深度，f表示突变位点的等位基因频率；

9c、基于上述9b中计算的Probability(G)和Probability(S)进行体细胞突变的判别，若Probability(S) > α，且 Probability(G) ≤ α则该位点判定为体细胞突变；若Probability(S) ≤ α，Probability(G) ≤ α，且 Purity ≥ 0.2则该位点判定为体细胞亚克隆突变；若Probability(S) ≤ α，且 Probability(G) > α则该位点判定为种系突变；其中α为判定阈值。

2.根据权利要求1所述的构建方法，其特征在于，所述步骤3）中，所述人参考基因组为hg19基因组。

3.根据权利要求1-2任一所述的构建方法，其特征在于，所述步骤5）中，所述拷贝数区段构建具体步骤如下：

5a、构建Bin区域：将杂交捕获的目标区域Bed划分成Bin；

5b、基线构建：选用正常对照样本获得每个样本的每个Bin的测序深度，建立正常拷贝数基线；

5c、区段划分：基于基线数据和肿瘤组织单样本Bin区域测序深度数据，使用循环二元分割CBS算法将每条染色体臂上拷贝数CN相近的Bin区域合并成区段Segment，并取Bin区域CN的中位值作为Segment的CN。

4.根据权利要求1-2任一所述的构建方法，其特征在于，所述步骤6）中，所述中国人群种系突变数据库的具体构建步骤如下：

6a、基于中国人群正常对照样本的NGS测序数据，进行种系突变的检测；

6b、基于上述6a中获得的种系突变，过滤低覆盖深度和低支持数的位点，获得自建中国人群种系突变集；

6c、基于上述6b种系突变集计算每个种系突变的等位基因频率及其在该人群中突变发生频率，完成自建种系突变数据库构建。

5.一种肿瘤突变负荷TMB评估模型的构建方法，其特征在于，包括权利要求1-4任一所述方法，并进一步包括如下步骤：

11）基于权利要求1-4任一所述步骤10）获得的体细胞突变结果及杂交捕获的基因编码区域大小，计算肿瘤突变负荷TMB值；

TMB=(MC/(Length(bp)))×1000000

其中，MC表示纳入TMB计算的体细胞突变数目，Length(bp)表示杂交捕获的基因编码区域大小。

6.一种计算机可读存储介质，其上存储有计算机程序/指令，其特征在于，该计算机程序/指令被处理器执行时实现权利要求 1-4任一所述方法的步骤。

Images