CN109949861B - 肿瘤突变负荷检测方法、装置和存储介质 - Google Patents
肿瘤突变负荷检测方法、装置和存储介质 Download PDFInfo
- Publication number
- CN109949861B CN109949861B CN201910254928.5A CN201910254928A CN109949861B CN 109949861 B CN109949861 B CN 109949861B CN 201910254928 A CN201910254928 A CN 201910254928A CN 109949861 B CN109949861 B CN 109949861B
- Authority
- CN
- China
- Prior art keywords
- mutation
- indel
- mutation frequency
- site
- frequency
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Abstract
一种肿瘤突变负荷检测方法、装置和存储介质，该方法包括：获取单个测试样本的突变频率数据，包括该样本目标区域的位点突变频率，将位点突变频率与设定的阈值进行比较得到大于阈值的单核苷酸变异，并去除单核苷酸变异中的无义突变得到有效单核苷酸变异个数；获取单个测试样本的Indel突变频率数据，包括该样本目标区域的Indel突变频率，将Indel突变频率与设定的阈值进行比较得到大于阈值的Indel突变个数；根据肿瘤突变负荷的估值公式计算肿瘤突变负荷的数值，估值公式包括有效单核苷酸变异个数的权重项和Indel突变个数的权重项。本发明的方法，在不依赖公共数据库和正常配对对照样本的前提下，准确地对肺癌样本的TMB指标进行检测。
Description
技术领域
本发明涉及肿瘤检测技术领域，具体涉及一种肿瘤突变负荷检测方法、装置和存储介质。
背景技术
细胞程序性死亡蛋白1(Programmed Cell Death protein 1,PD-1)是一种通常表达于细胞表面的蛋白，通过降低免疫细胞对细胞的炎症反应而调控免疫系统，防止自身免疫的发生。PD-1的配体PD-L1可以特异性地中和PD-1，从而重新启动免疫系统对细胞的杀伤作用。这种现象又被叫做免疫检查点抑制。通过免疫检查点抑制机制(如CTLA-4和PD-L1)开发的药物在近几年的癌症治疗中被发现有着令人鼓舞的治疗效果。
细胞的癌变通常是由体细胞中的基因突变长期积累的结果，但不是所有的体细胞突变都会导致细胞癌变。目前主流观点认为，只有在驱动基因上的特异突变才赋予细胞癌变的特性，这种突变叫做驱动突变(driver mutation)。而驱动突变又会引发其他的基因突变，这导致癌细胞中的基因突变数量往往高于正常的体细胞。肿瘤突变负荷(TMB)是反映肿瘤细胞中总的基因突变程度的一个指标，通常以每百万碱基(Mb)的肿瘤基因组区域中包含的肿瘤体细胞突变总数来表示。
多个大规模临床研究发现，免疫检查点抑制剂的疗效很大程度上取决于患者癌细胞中所携带的基因突变的数量。在接受免疫检查点抑制剂治疗的患者中，TMB高与TMB低的患者对免疫疗法的疗效差异十分明显。因此，TMB的精确测量可以预测免疫检查点抑制剂的疗效，使癌症患者有机会获得更加精准的治疗。
最初的TMB采用全外显子组测序方法，对照患者的正常组织和癌组织，找出癌细胞特有的体细胞突变。可是对癌组织和正常组织同时测序成本较高，科学家开始探索只对癌症单个样本测序来测量TMB的可行性。2017年底，FDA批准的FoundationOne CDx基因检测试剂盒就采用了单样本测量TMB的方法，即只对癌症样本测序，然后通过统计方法和人群数据库信息在基因突变中确定胚系突变并过滤。这在一定程度上降低了TMB检测的成本要求。但是这类方法有几个缺陷。首先，不同癌种间TMB的差异较大，用同样的数据库过滤方法并不能保证对所有癌种的胚系基因突变进行精确的去除；其次，这种方法高度依赖公共数据库的质量和多样性，对不在数据库中的人种的突变背景无法做出精确的描述。
发明内容
本申请提供一种肿瘤突变负荷检测方法、装置和存储介质，在不依赖公共数据库和配对正常样本的前提下，准确地对肺癌样本的TMB指标进行检测。
根据第一方面，一种实施例中提供一种肿瘤突变负荷检测方法，包括如下步骤：
获取单个测试样本的突变频率数据，该突变频率数据包括该样本目标区域的位点突变频率，将上述位点突变频率与设定的位点突变频率阈值进行比较，得到位点突变频率大于上述位点突变频率阈值的单核苷酸变异，并去除上述单核苷酸变异中的无义突变得到有效单核苷酸变异个数；
获取单个测试样本的Indel突变频率数据，该Indel突变频率数据包括该样本目标区域的Indel突变频率，将上述Indel突变频率与设定的Indel突变频率阈值进行比较，得到Indel突变频率大于上述Indel突变频率阈值的Indel突变个数；
根据肿瘤突变负荷的估值公式计算肿瘤突变负荷的数值，其中上述估值公式包括上述有效单核苷酸变异个数的权重项和上述Indel突变个数的权重项。
在优选实施例中，上述肿瘤突变负荷的估值公式如下：
S/100+sgn(I)
其中，S指上述有效单核苷酸变异个数，I指Indel突变个数，sgn()为符号函数，在I大于或等于个数阈值的情况下，sgn(I)输出值为1，否则输出值为0。
在优选实施例中，上述个数阈值为2。
在优选实施例中，上述测试样本的位点突变频率是根据上述测试样本的测序深度进行校正得到的校正位点突变频率。
在优选实施例中，上述校正位点突变频率通过如下公式得到：
θadj＝θj×τ(min(1，Dj/Dlimit)，α，β)；
其中，θadj为在当前位点的校正位点突变频率，θj为在当前位点实际观测的位点突变频率，τ是以α和β为形状参数的Beta分布的累计概率分布函数，Dj为上述测试样本在当前位点的实际测序深度，Dlimit为设定的最低校正测序深度。
在优选实施例中，上述位点突变频率阈值通过如下方法确定，该方法即SNV变异训练或SNV变异统计方法：
获取一组训练样本的ACGT格式文件，该ACGT格式文件包含选定的目标区域中每个位点的位置信息、测序深度信息和突变到任意其它三种非参考碱基的突变频率；
统计所有训练样本中每种三碱基突变(mutational signature)的平均突变频率，并将其作为每种三碱基突变的先验突变频率；
从所有训练样本中提取每一位点的突变频率最大值，上述突变频率最大值满足的条件是，在同一个位点SNP的比例高于阈值p的次数在所有训练样本中至少出现设定次数n，若没有满足上述条件，该突变位点的突变频率被置换为该位点所有训练样本的突变频率平均值；将所得到的突变频率最大值或突变频率平均值乘以设定的系数值后得到的数值如果大于1，则将突变频率设为1，如果该数值小于1，则突变频率取该数值，然后将突变频率与先验突变频率进行加权，得到该突变位点的加权突变频率，作为位点后验突变频率；
求出所有训练样本在上述目标区域中各个位点的测序深度平均值；
将上述位点后验突变频率和上述测序深度平均值提供给TNER方法，在给定的显著性水平下得出上述位点突变频率阈值。
在优选实施例中，上述阈值p是0.05，上述设定次数n是10，上述系数值是5，上述给定的显著性水平是0.001。
在优选实施例中，上述先验突变频率通过如下方法确定：
获取每个训练样本中的背景突变和SNP突变，其中杂合子和纯合子SNP的突变频率分别在0.5和1处聚集并呈高斯分布，而背景突变的突变频率在0.001-0.1处聚集并呈伽马分布；
通过对上述背景突变和SNP突变的突变频率形成的混合分布进行拟合，找出混合分布的概率密度分布中背景突变与杂合子SNP突变之间概率密度分布的最低点所对应的突变频率，将该突变频率作为背景突变频率的阈值，将突变频率小于该阈值的突变作为真实背景突变；
在上述真实背景突变中对每个训练样本的相同的三碱基突变进行归类后求三碱基平均背景突变频率，然后将所有训练样本中相同的三碱基平均背景突变频率的平均值作为三碱基突变的先验突变频率。
上述三碱基突变(mutational signature)是指6种基础单碱基突变形式(A→T、A→G、A→C、C→A、C→T、C→G)与其上下文各一个碱基的组合，共有96种。
在优选实施例中，上述方法在得到位点突变频率大于上述位点突变频率阈值的单核苷酸变异之后，去除变异频率在5％以下、45％～55％之间和95％～100％之间的单核苷酸变异，再去除上述单核苷酸变异中的无义突变得到有效单核苷酸变异个数。
在优选实施例中，上述Indel突变频率阈值通过如下方法确定，该方法即Indel变异训练或Indel变异统计方法：
获取一组训练样本的Indel格式文件，该Indel格式文件包含选定的目标Indel组中每个Indel的信息，将每个Indel以染色体+位置+突变前碱基+突变类型+突变后碱基进行编码，在编码过程中只选取突变后碱基编码的第一位组成每个Indel的突变编码；
在所有训练样本中找出所有至少出现两次且突变频率都大于频率预设值的Indel编码及其对应的突变频率，每个Indel编码对应的突变频率为所有训练样本中该编码对应的突变频率最大值，将该突变频率最大值的设定倍数作为Indel检测中的上述Indel突变频率阈值，并将上述编码及其对应的上述Indel突变频率阈值以哈希表的形式保存。
在优选实施例中，上述频率预设值为1％，上述设定倍数为2倍。
在优选实施例中，上述方法在将上述Indel突变频率与设定的Indel突变频率阈值进行比较时，将没有出现在上述哈希表中的Indel排除。
在优选实施例中，上述方法在将上述Indel突变频率与设定的Indel突变频率阈值进行比较时，去除突变频率在40％～60％之间和90％～100％之间的Indel。
根据第二方面，一种实施例中提供一种肿瘤突变负荷检测装置，包括如下单元：
单核苷酸变异统计单元，用于获取单个测试样本的突变频率数据，上述突变频率数据包括该样本目标区域的位点突变频率，将上述位点突变频率与设定的位点突变频率阈值进行比较，得到位点突变频率大于上述位点突变频率阈值的单核苷酸变异，并去除上述单核苷酸变异中的无义突变得到有效单核苷酸变异个数；
Indel变异统计单元，用于获取单个测试样本的Indel突变频率数据，上述Indel突变频率数据包括该样本目标区域的Indel突变频率，将上述Indel突变频率与设定的Indel突变频率阈值进行比较，得到Indel突变频率大于上述Indel突变频率阈值的Indel突变个数；
肿瘤突变负荷计算单元，用于根据肿瘤突变负荷的估值公式计算肿瘤突变负荷的数值，其中上述估值公式包括上述有效单核苷酸变异个数的权重项和上述Indel突变个数的权重项。
根据第三方面，一种实施例中提供一种计算机可读存储介质，包括程序，该程序能够被处理器执行以实现如第一方面的方法。
本发明的肿瘤突变负荷检测方法，针对单个样本进行检测，结合了单核苷酸变异检测和插入缺失检测两个功能。本发明的方法通过学习同类正常样本(训练样本)的测序结果中单核苷酸变异和插入缺失情况，对多种癌症样本(特别是肺癌样本)中的相应突变，根据测序深度和突变频率进行检测，从而达到计算肿瘤突变负荷的目的。本发明的方法不需要参考公共数据库中的突变信息，也不需要提取成对的正常样本并进行测序。与现有的双样本肿瘤突变负荷检测技术相比，本发明的方法降低了实验操作的人力消耗和对患者的取样难度，也降低了在计算新样本肿瘤突变负荷时的计算资源消耗。与现有的单样本肿瘤突变负荷检测流程相比，本发明的方法减少了对公共数据库的依赖，并且可以在检测突变的同时根据同类样本的变异频率信息过滤掉样本中的变异假阳性。
附图说明
图1为本发明实施例中一种肿瘤突变负荷检测方法流程图。
图2为本发明实施例中Indel变异统计结果与真实全外显子组测序TMB数值间的关系。
图3为本发明实施例中一种肿瘤突变负荷检测装置结构框图。
图4为本发明实施例中在肺癌样本中TMB数值的计算结果与全外显子组测序TMB数值的线性关系，x轴为外显子组测序TMB，y轴为本方法计算的TMB数值，原点大小表示样本测序深度，原点越大，测序深度越高。
图5为本发明实施例中双样本方法在相同肺癌样本中TMB数值的计算结果与全外显子组测序TMB数值的线性关系，x轴为外显子组测序TMB，y轴为双样本方法在相同肺癌样本中计算的TMB数值，原点大小表示样本测序深度，原点越大，测序深度越高。
图6为本发明实施例中在23例不同癌症样本中TMB数值的计算结果与全外显子组测序TMB数值的线性关系，x轴为外显子组测序TMB，y轴为本方法计算的TMB数值，原点大小表示样本测序深度，原点越大，测序深度越高。
具体实施方式
下面通过具体实施方式结合附图对本发明作进一步详细说明。在以下的实施方式中，很多细节描述是为了使得本申请能被更好的理解。然而，本领域技术人员可以毫不费力的认识到，其中部分特征在不同情况下是可以省略的，或者可以由其他材料、方法所替代。
另外，说明书中所描述的特点、操作或者特征可以以任意适当的方式结合形成各种实施方式。同时，方法描述中的各步骤或者动作也可以按照本领域技术人员所能显而易见的方式进行顺序调换或调整。因此，说明书和附图中的各种顺序只是为了清楚描述某一个实施例，并不意味着是必须的顺序，除非另有说明其中某个顺序是必须遵循的。
本发明中采用的术语具体含义如下：
参考基因组：物种参考的标准基因组序列。
读长(Reads)：测序所得基因组序列片段。
BAM：一种用于存储比对信息的标准二进制文件格式。
acgt：一种记录每个位点单核苷酸变异信息的文件。
Indel：一种记录每个位点插入和缺失类型变异信息的文件。
如图1所示，本发明的一种实施例中提供一种肿瘤突变负荷检测方法，包括如下步骤：
S101：获取单个测试样本的突变频率数据，该突变频率数据包括该样本目标区域的位点突变频率，将上述位点突变频率与设定的位点突变频率阈值进行比较，得到位点突变频率大于上述位点突变频率阈值的单核苷酸变异，并去除上述单核苷酸变异中的无义突变得到有效单核苷酸变异个数；
S102：获取单个测试样本的Indel突变频率数据，该Indel突变频率数据包括该样本目标区域的Indel突变频率，将上述Indel突变频率与设定的Indel突变频率阈值进行比较，得到Indel突变频率大于上述Indel突变频率阈值的Indel突变个数；
S103：根据肿瘤突变负荷的估值公式计算肿瘤突变负荷的数值，其中上述估值公式包括上述有效单核苷酸变异个数的权重项和上述Indel突变个数的权重项。
本发明的单样本肿瘤突变负荷(TMB)检测方法，结合了单核苷酸变异(SNV)检测和插入缺失(Indel)检测两个功能。通过学习同类正常组织(训练样本)的测序结果中单核苷酸变异和插入缺失情况，对肺癌样本(测试样本)中的相应突变，根据测序深度和突变频率进行检测，从而达到计算TMB的目的。
以下分别介绍单核苷酸变异(SNV)和插入缺失(Indel)统计。
I.SNV变异统计
将癌组织样本(测试样本)特别是肺癌组织样本，与正常对照样本(训练样本)的测序读长(reads)经过比对产生的BAM格式文件用Samtools软件转化为pileup格式文件。在转化过程中，只允许测序错误和比对错误率小于0.1％的reads，对应的Phread Score和Mapping Score均为30。再将生成的pileup格式文件用sequenza-utils以默认参数转化为单核苷酸突变频率数据文件ACGT格式，或者对pileup格式文件用Varscan中的pileup2acgt工具转化为单核苷酸突变频率数据文件ACGT。ACGT格式类型的文件包含了目标区域中每个位点的位置信息、深度信息和变异/突变到任意其他三种非参考碱基的概率。
然后，(1)统计所有正常对照样本(训练样本)中每种三碱基突变(mutationalsignature)的平均突变频率，并将它作为TNER方法(参考Shibing Deng,Maruja Lira,Donghui Huang,Kai Wang,Crystal Valdez,Jennifer Kinong,Paul A.Rejto,JadwigaBienkowska,James Hardwick,Tao Xie.TNER:A Novel Bayesian Background ErrorSuppression Method for Mutation Detection in Circulating Tumor DNA,BMCBioinformatics,(2018)19:387)检验新突变的先验突变频率。(2)寻求所有正常对照样本(训练样本)中每一位点突变频率中的最大值，该最大值需要满足的条件是：在同一个位点SNP的比例高于某个阈值p(例如，0.05)的次数在所有正常对照样本(训练样本)中至少出现n(例如，10)次。如果没有满足上述条件，该突变位点的突变频率被置换为该位点在所有训练样本中的突变频率平均值。然后将所得到的突变频率，即突变频率最大值或突变频率平均值乘以设定的系数值(例如，5)后得到的数值如果大于1，则将突变频率设为1，如果该数值小于1，则突变频率取该数值，然后将突变频率与先验突变频率进行加权，得到该突变位点的加权突变频率，作为位点后验突变频率。(3)求出所有正常对照样本(训练样本)在目标区域中各个位点的测序深度平均值。(4)通过TNER方法在给定的显著性水平下(例如，0.001)结合位点后验突变频率和该位点平均测序深度，得出该位点突变频率阈值，作为检验新输入样本(测试样本)的突变频率的阈值。
在本发明一个优选的实施例中，先验突变频率通过如下方法确定：获取每个训练样本中的背景突变和SNP突变，其中杂合子和纯合子SNP的突变频率分别在0.5和1处聚集并呈高斯分布，而背景突变的突变频率在0.001-0.1处聚集并呈伽马分布；通过对背景突变和SNP突变的突变频率形成的混合分布进行拟合，找出混合分布的概率密度分布中背景突变与杂合子SNP突变之间概率密度分布的最低点所对应的突变频率，将该突变频率作为背景突变频率的阈值，将突变频率小于该阈值的突变作为真实背景突变；在真实背景突变中对每个训练样本的相同的三碱基突变进行归类后求三碱基平均背景突变频率，然后将所有训练样本中相同的三碱基平均背景突变频率的平均值作为三碱基突变的先验突变频率。
然后，根据输入样本(测试样本)的深度信息，对测试样本的位点突变频率进行校正，通过如下公式得到校正位点突变频率：
θadj＝θj×τ(min(1，Dj/Dlimit)，α，β)；
其中，θadj为在当前位点的校正位点突变频率，θj为在当前位点实际观测的位点突变频率，τ是以α和β为形状参数的Beta分布的累计概率分布函数，Dj为上述测试样本在当前位点的实际测序深度，Dlimit为设定的最低校正测序深度。
通过比较癌组织样本(测试样本)相应位置的测序突变频率(优选校正位点突变频率)和突变频率的阈值，输出测试样本中通过筛选的单核苷酸变异，即测序突变频率大于突变频率的阈值的单核苷酸变异。为提高准确度，在通过筛选的单核苷酸变异中，去除突变频率在5％以下、45％-55％之间和95％-100％之间的单核苷酸变异，然后用snpeff软件注释后得到这些单核苷酸变异中无义突变的个数。用总突变个数减去无义突变的个数，得到单核苷酸变异个数，作为单核苷酸变异统计部分的输出值参与之后的TMB检验。
II.Indel变异统计
将正常对照样本(训练样本)测序中的每个Indel以染色体+位置+突变前碱基+突变类型+突变后碱基作为编码。在编码过程中，只选取突变后碱基编码的第一位组成每个Indel的突变编码。在所有正常对照样本(训练样本)中找出所有至少出现两次并突变频率都大于频率预设值(例如，1％)的Indel编码和对应的突变频率，每个Indel编码对应的突变频率为所有正常组织样本中该编码对应突变频率的最大值。将该最大值的设定倍数(例如，2倍)作为新Indel检测中突变频率的阈值。然后，将编码和它对应的阈值以哈希表的形式储存。
在检测癌组织样本(测试样本)中新的Indel时，对于每一个癌组织样本检测出的Indel，如果该Indel出现在哈希表中并且突变频率大于阈值，该Indel就被定义为通过筛选的Indel。在通过筛选的Indel中，去除突变频率在40％-60％之间和90％-100％之间的Indel。如果癌组织样本中通过筛选的Indel数目大于或等于2个，那么定义这个样本为高Indel样本，反之则为低Indel样本。图2反映了Indel变异统计结果与真实全外显子组测序TMB数值间的关系，高Indel样本组对应的平均TMB显著高于低Indel样本组对应的平均TMB。
最终TMB的估计值结合了单核苷酸变异个数和Indel统计个数对样本TMB的判断，TMB的估计值用下式表示：
S/100+sgn(I)
其中，S指上述有效单核苷酸变异个数，I指Indel突变个数，sgn()为符号函数，在I大于或等于个数阈值(例如，2)的情况下，sgn(I)输出值为1，否则输出值为0。
本发明的方法，通过学习同类正常样本的变异信息，实现了基于测序变异频率对肺癌样本的TMB指标的准确检测。本发明的方法不需要参考公共数据库中的突变信息，也不需要提取成对的正常样本并进行测序。与现有的双样本肿瘤突变负荷检测技术相比，本发明的方法降低了实验操作的人力消耗和对患者的取样难度，也降低了在计算新样本肿瘤突变负荷时的计算资源消耗。与现有的单样本肿瘤突变负荷检测流程相比，本发明的方法减少了对公共数据库的依赖，并且可以在检测突变的同时根据同类样本的变异频率信息过滤掉样本中的变异假阳性。
本领域技术人员可以理解，上述实施方式中各种方法的全部或部分功能可以通过硬件的方式实现，也可以通过计算机程序的方式实现。当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘、光盘、硬盘等，通过计算机执行该程序以实现上述功能。例如，将程序存储在设备的存储器中，当通过处理器执行存储器中程序，即可实现上述全部或部分功能。另外，当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序也可以存储在服务器、另一计算机、磁盘、光盘、闪存盘或移动硬盘等存储介质中，通过下载或复制保存到本地设备的存储器中，或对本地设备的系统进行版本更新，当通过处理器执行存储器中的程序时，即可实现上述实施方式中全部或部分功能。
因此，对应于本发明的方法，本发明一种实施例中提供一种肿瘤突变负荷检测装置，如图3所示，包括如下单元：单核苷酸变异统计单元301，用于获取单个测试样本的突变频率数据，上述突变频率数据包括该样本目标区域的位点突变频率，将上述位点突变频率与设定的位点突变频率阈值进行比较，得到位点突变频率大于上述位点突变频率阈值的单核苷酸变异，并去除上述单核苷酸变异中的无义突变得到有效单核苷酸变异个数；Indel变异统计单元302，用于获取单个测试样本的Indel突变频率数据，上述Indel突变频率数据包括该样本目标区域的Indel突变频率，将上述Indel突变频率与设定的Indel突变频率阈值进行比较，得到Indel突变频率大于上述Indel突变频率阈值的Indel突变个数；肿瘤突变负荷计算单元303，用于根据肿瘤突变负荷的估值公式计算肿瘤突变负荷的数值，其中上述估值公式包括上述有效单核苷酸变异个数的权重项和上述Indel突变个数的权重项。
此外，本发明的一种实施例中提供一种计算机可读存储介质，包括程序，该程序能够被处理器执行以实现如第一方面的方法。
以下通过实施例详细说明本发明的技术方案，应当理解，实施例仅是示例性的，不能理解为对本发明保护范围的限制。
实施例1：
以下实施例中，训练数据：360例健康人群对照样本经过深圳裕策生物科技有限公司YuceOne Plus芯片测序生成的单核苷酸突变频率数据和110例健康人群对照样本经过深圳裕策生物科技有限公司YuceOne Plus芯片测序生成的Indel突变频率数据。测试数据：11例肺癌患者组织样本经过深圳裕策生物科技有限公司YuceOne Plus芯片测序获得的突变结果，包含单核苷酸突变频率数据和Indel突变频率数据。
将360例健康对照样本经过深圳裕策生物科技有限公司YuceOne Plus芯片测序，将测序所得读长(reads)与人参考基因组进行比对，得到BAM格式的比对结果。然后对BAM格式文件用Samtools软件转化为pileup格式文件。在转化过程中，只允许测序错误和比对错误率小于0.1％的reads，对应的Phread得分(Phread Score)和映射得分(Mapping Score)均为30。再将生成的pileup格式文件用sequenza-utils以默认参数转化为单核苷酸突变频率数据文件ACGT格式。将该ACGT文件作为SNV变异统计的输入数据，根据SNV变异统计的描述进行训练，获得每个位点突变频率阈值，用来对测试数据的目标区域进行单核苷酸变异的检测。
将110例健康对照样本经过深圳裕策生物科技有限公司YuceOne Plus芯片测序，将测序所得reads与人参考基因组进行比对，得到BAM格式的比对结果。将比对产生的BAM格式文件用Samtools软件转化为pileup格式文件。在转化过程中，只允许测序错误和比对错误率小于0.1％的reads，对应的Phread得分(Phread Score)和映射得分(Mapping Score)均为30。然后，将pileup文件用Varscan pileup2indel以默认参数转化为Indel格式。Indel格式类型的文件包含了目标区域中所有插入和缺失(Indel)的位置信息、深度信息和与具体插入缺失片段相对应的突变频率。将Indel文件根据Indel变异统计的描述进行处理，得到每个Indel突变频率阈值，用于对测试数据的目标区域进行Indel检测。
将11例肺癌组织测试样本经过深圳裕策生物科技有限公司YuceOne Plus芯片测序，将测序所得读长(reads)与人参考基因组进行比对，得到BAM格式的比对结果。(1)对BAM格式文件用Samtools软件转化为pileup格式文件，在转化过程中，只允许测序错误和比对错误率小于0.1％的reads，对应的Phread得分(Phread Score)和映射得分(MappingScore)均为30，再将生成的pileup格式文件用sequenza-utils以默认参数转化为单核苷酸突变频率数据文件ACGT格式；(2)对比对产生的BAM格式文件用Samtools软件转化为pileup格式文件，在转化过程中，只允许测序错误和比对错误率小于0.1％的reads，对应的Phread得分(Phread Score)和映射得分(Mapping Score)均为30，然后将pileup文件用Varscanpileup2indel以默认参数转化为Indel格式。对获得的单核苷酸突变频率数据和Indel突变频率数据用本发明方法计算得出TMB，然后将相同样本和与其配对的对照样本用全外显子组双样本测序方法计算TMB。将两种方法得出的TMB进行拟合，拟合结果R-square为0.76，Adjusted R-square为0.71，皮尔森相关系数为0.79(图4)。
然后，将上述11例肺癌组织测试样本和与其配对的对照样本用深圳裕策生物科技有限公司YuceOne Plus芯片双样本测序方法计算TMB(双样本测序方法计算TMB参见专利文献CN109033749A)，然后与相同样本全外显子组双样本测序方法得到的TMB进行拟合。拟合结果R-square(R2)为0.70，Adjusted R-square为0.68，皮尔森相关系数为0.84(图5)。比较本发明方法和双样本测序方法与外显子测序方法TMB的拟合结果可以发现，本发明方法可以做到在不牺牲精确度的前提下取代肺癌双样本测序TMB检测方法(表1)。
表1本发明方法和双样本方法与全外显子组测序TMB数值的相关关系
实施例2：
将23例各类癌组织测试样本经过深圳裕策生物生物科技有限公司YuceOne Plus芯片测序，对获得的单核苷酸突变频率数据和Indel突变频率数据用本发明方法计算TMB，然后将相同样本和与其配对的对照样本用全外显子组双样本测序方法计算TMB。将两种方法得出的TMB进行拟合，拟合结果R-square为0.72，Adjusted R-square为0.68，皮尔森相关系数为0.85(图6)。与实施例1中11例肺癌结果对比表明，本发明方法也可以适用于其他癌种的TMB检测。
以上应用了具体个例对本发明进行阐述，只是用于帮助理解本发明，并不用以限制本发明。对于本发明所属技术领域的技术人员，依据本发明的思想，还可以做出若干简单推演、变形或替换。
Claims (15)
1.一种肿瘤突变负荷检测方法，其特征在于，所述方法包括如下步骤：
获取单个测试样本的突变频率数据，所述突变频率数据包括该样本目标区域的位点突变频率，将所述位点突变频率与设定的位点突变频率阈值进行比较，得到位点突变频率大于所述位点突变频率阈值的单核苷酸变异，并去除所述单核苷酸变异中的无义突变得到有效单核苷酸变异个数；
获取单个测试样本的Indel突变频率数据，所述Indel突变频率数据包括该样本目标区域的Indel突变频率，将所述Indel突变频率与设定的Indel突变频率阈值进行比较，得到Indel突变频率大于所述Indel突变频率阈值的Indel突变个数；
根据肿瘤突变负荷的估值公式计算肿瘤突变负荷的数值，其中所述估值公式包括所述有效单核苷酸变异个数的权重项和所述Indel突变个数的权重项。
2.根据权利要求1所述的方法，其特征在于，所述肿瘤突变负荷的估值公式如下：
S/100+sgn(I)
其中，S指所述有效单核苷酸变异个数，I指Indel突变个数，sgn()为符号函数，在I大于或等于个数阈值的情况下，sgn(I)输出值为1，否则输出值为0。
3.根据权利要求2所述的方法，其特征在于，所述个数阈值为2。
4.根据权利要求1所述的方法，其特征在于，所述测试样本的位点突变频率是根据所述测试样本的测序深度进行校正得到的校正位点突变频率。
5.根据权利要求4所述的方法，其特征在于，所述校正位点突变频率通过如下公式得到：
θadj＝θj×T(min(l，Dj/Dlimit)，α，β)；
其中，θadj为在当前位点的校正位点突变频率，θj为在当前位点实际观测的位点突变频率，T 是以α和β为形状参数的Beta分布的累计概率分布函数，Dj为所述测试样本在当前位点的实际测序深度，Dlimit为设定的最低校正测序深度。
6.根据权利要求1所述的方法，其特征在于，所述位点突变频率阈值通过如下方法确定：
获取一组训练样本的ACGT格式文件，该ACGT格式文件包含目标区域中每个位点的位置信息、测序深度信息和突变到任意其它三种非参考碱基的突变频率；
统计所有训练样本中每种三碱基突变的平均突变频率，并将其作为每种三碱基突变的先验突变频率；
从所有训练样本中提取每一位点的突变频率最大值，所述突变频率最大值满足的条件是，在同一个位点SNP的比例高于阈值p的次数在所有训练样本中至少出现设定次数n，若没有满足所述条件，该突变位点的突变频率被置换为该位点在所有训练样本中的突变频率平均值；将所得到的所述突变频率最大值或突变频率平均值乘以设定的系数值后得到的数值如果大于1，则将突变频率设为1，如果所述数值小于1，则突变频率取所述数值，然后将所述突变频率与先验突变频率进行加权，得到该突变位点的加权突变频率，作为位点后验突变频率；
求出所有训练样本在所述目标区域中各个位点的测序深度平均值；
将所述位点后验突变频率和所述测序深度平均值提供给TNER方法，在给定的显著性水平下得出所述位点突变频率阈值。
7.根据权利要求6所述的方法，其特征在于，所述阈值p是0.05，所述设定次数n是10，所述系数值是5，所述给定的显著性水平是0.001。
8.根据权利要求6所述的方法，其特征在于，所述先验突变频率通过如下方法确定：
获取每个训练样本中的背景突变和SNP突变，其中杂合子和纯合子SNP的突变频率分别在0.5和1处聚集并呈高斯分布，而背景突变的突变频率在0.001-0.1处聚集并呈伽马分布；
通过对所述背景突变和SNP突变的突变频率形成的混合分布进行拟合，找出混合分布的概率密度分布中背景突变与杂合子SNP突变之间概率密度分布的最低点所对应的突变频率，将该突变频率作为背景突变频率的阈值，将突变频率小于该阈值的突变作为真实背景突变；
在所述真实背景突变中对每个训练样本的相同的三碱基突变进行归类后求三碱基平均背景突变频率，然后将所有训练样本中相同的三碱基平均背景突变频率的平均值作为三碱基突变的先验突变频率。
9.根据权利要求1所述的方法，其特征在于，所述方法在得到位点突变频率大于所述位点突变频率阈值的单核苷酸变异之后，去除变异频率在5％以下、45％～55％之间和95％～100％之间的单核苷酸变异，再去除所述单核苷酸变异中的无义突变得到有效单核苷酸变异个数。
10.根据权利要求1所述的方法，其特征在于，所述Indel突变频率阈值通过如下方法确定：
获取一组训练样本的Indel格式文件，该Indel格式文件包含选定的目标Indel组中每个Indel的信息，将每个Indel以染色体+位置+突变前碱基+突变类型+突变后碱基进行编码，在编码过程中只选取突变后碱基编码的第一位组成每个Indel的突变编码；
在所有训练样本中找出所有至少出现两次且突变频率都大于频率预设值的Indel编码及其对应的突变频率，每个Indel编码对应的突变频率为所有训练样本中该编码对应的突变频率最大值，将该突变频率最大值的设定倍数作为Indel检测中的所述Indel突变频率阈值，并将所述编码及其对应的所述Indel突变频率阈值以哈希表的形式保存。
11.根据权利要求10所述的方法，其特征在于，所述频率预设值为1％，所述设定倍数为2倍。
12.根据权利要求10所述的方法，其特征在于，所述方法在将所述Indel突变频率与设定的Indel突变频率阈值进行比较时，将没有出现在所述哈希表中的Indel排除。
13.根据权利要求1所述的方法，其特征在于，所述方法在将所述Indel突变频率与设定的Indel突变频率阈值进行比较时，去除突变频率在40％～60％之间和90％～100％之间的Indel。
14.一种肿瘤突变负荷检测装置，其特征在于，所述装置包括如下单元：
单核苷酸变异统计单元，用于获取单个测试样本的突变频率数据，所述突变频率数据包括该样本目标区域的位点突变频率，将所述位点突变频率与设定的位点突变频率阈值进行比较，得到位点突变频率大于所述位点突变频率阈值的单核苷酸变异，并去除所述单核苷酸变异中的无义突变得到有效单核苷酸变异个数；
Indel变异统计单元，用于获取单个测试样本的Indel突变频率数据，所述Indel突变频率数据包括该样本目标区域的Indel突变频率，将所述Indel突变频率与设定的Indel突变频率阈值进行比较，得到Indel突变频率大于所述Indel突变频率阈值的Indel突变个数；
肿瘤突变负荷计算单元，用于根据肿瘤突变负荷的估值公式计算肿瘤突变负荷的数值，其中所述估值公式包括所述有效单核苷酸变异个数的权重项和所述Indel突变个数的权重项。
15.一种计算机可读存储介质，其特征在于，包括程序，所述程序能够被处理器执行以实现如权利要求1-13中任一项所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910254928.5A CN109949861B (zh) | 2019-03-29 | 2019-03-29 | 肿瘤突变负荷检测方法、装置和存储介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910254928.5A CN109949861B (zh) | 2019-03-29 | 2019-03-29 | 肿瘤突变负荷检测方法、装置和存储介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109949861A CN109949861A (zh) | 2019-06-28 |
| CN109949861B true CN109949861B (zh) | 2020-02-21 |
Family
ID=67013275
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910254928.5A Active CN109949861B (zh) | 2019-03-29 | 2019-03-29 | 肿瘤突变负荷检测方法、装置和存储介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109949861B (zh) |
Families Citing this family (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110343748B (zh) * | 2019-08-08 | 2020-05-26 | 广州燃石医学检验所有限公司 | 基于高通量靶向测序分析肿瘤突变负荷的方法 |
| CN112840402B (zh) * | 2019-09-02 | 2024-05-07 | 北京哲源科技有限责任公司 | 获得细胞内确定性事件的方法及电子设备 |
| CN110993028B (zh) * | 2019-12-17 | 2022-03-29 | 清华大学 | 突变数据识别方法、训练方法、处理装置及存储介质 |
| CN111584002B (zh) * | 2020-05-22 | 2022-04-29 | 至本医疗科技（上海）有限公司 | 用于检测肿瘤突变负荷的方法、计算设备和计算机存储介质 |
| CN111477277A (zh) * | 2020-05-29 | 2020-07-31 | 北京优迅医学检验实验室有限公司 | 样本质量评估方法和装置 |
| CN111933219B (zh) * | 2020-09-16 | 2021-06-08 | 北京求臻医学检验实验室有限公司 | 一种分子标志物肿瘤缺失突变负荷的检测方法 |
| CN113481299B (zh) * | 2021-06-30 | 2022-05-10 | 苏州京脉生物科技有限公司 | 用于肺癌检测的靶向测序panel、试剂盒及获得靶向测序panel的方法 |
| CN114005489B (zh) * | 2021-12-28 | 2022-03-22 | 成都齐碳科技有限公司 | 基于三代测序数据检测点突变的分析方法和装置 |
| CN115424664B (zh) * | 2022-11-07 | 2023-03-10 | 北京雅康博生物科技有限公司 | 人为突变程度评估方法及装置 |
Family Cites Families (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20160002717A1 (en) * | 2014-07-02 | 2016-01-07 | Boreal Genomics, Inc. | Determining mutation burden in circulating cell-free nucleic acid and associated risk of disease |
| CN104462869B (zh) * | 2014-11-28 | 2017-12-26 | 天津诺禾致源生物信息科技有限公司 | 检测体细胞单核苷酸突变的方法和装置 |
| WO2017191076A1 (en) * | 2016-05-01 | 2017-11-09 | Genome Research Limited | Method of characterising a dna sample |
| CN106566877A (zh) * | 2016-10-31 | 2017-04-19 | 天津诺禾致源生物信息科技有限公司 | 检测基因突变的方法和装置 |
| CN108690871B (zh) * | 2018-03-29 | 2022-05-20 | 深圳裕策生物科技有限公司 | 基于二代测序的插入缺失突变检测方法、装置和存储介质 |
| CN109033749B (zh) * | 2018-06-29 | 2020-01-14 | 裕策医疗器械江苏有限公司 | 一种肿瘤突变负荷检测方法、装置和存储介质 |
| CN109411015B (zh) * | 2018-09-28 | 2020-12-22 | 深圳裕策生物科技有限公司 | 基于循环肿瘤dna的肿瘤突变负荷检测装置及存储介质 |
-
2019
- 2019-03-29 CN CN201910254928.5A patent/CN109949861B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN109949861A (zh) | 2019-06-28 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109949861B (zh) | 肿瘤突变负荷检测方法、装置和存储介质 | |
| Adie et al. | Speeding disease gene discovery by sequence based candidate prioritization | |
| US20220130488A1 (en) | Methods for detecting copy-number variations in next-generation sequencing | |
| US8725422B2 (en) | Methods for estimating genome-wide copy number variations | |
| CN108733975B (zh) | 基于二代测序的肿瘤克隆变异检测方法、装置和存储介质 | |
| KR102237923B1 (ko) | 암 검출을 위한 혈장 dna의 돌연변이 분석 | |
| US11043283B1 (en) | Systems and methods for automating RNA expression calls in a cancer prediction pipeline | |
| US20160117444A1 (en) | Methods for determining absolute genome-wide copy number variations of complex tumors | |
| Chiaromonte et al. | The share of human genomic DNA under selection estimated from human–mouse genomic alignments | |
| CN110010197B (zh) | 基于血液循环肿瘤dna的单核苷酸变异检测方法、装置和存储介质 | |
| EP2844771A1 (en) | Methods for determining absolute genome-wide copy number variations of complex tumors | |
| CN111718982A (zh) | 一种肿瘤组织单样本体细胞突变检测方法及装置 | |
| CN111916150A (zh) | 一种基因组拷贝数变异的检测方法和装置 | |
| CN111755068A (zh) | 基于测序数据识别肿瘤纯度和绝对拷贝数的方法及装置 | |
| Kato et al. | Sweepstake evolution revealed by population-genetic analysis of copy-number alterations in single genomes of breast cancer | |
| CN114530199A (zh) | 基于双重测序数据检测低频突变的方法、装置及存储介质 | |
| CN112116956A (zh) | 一种基于二代测序的肿瘤单样本tmb检测方法及装置 | |
| CN114974415A (zh) | 一种检测染色体拷贝数异常的方法和装置 | |
| CN114067909B (zh) | 一种矫正同源重组缺陷评分的方法、装置和存储介质 | |
| CN114242170B (zh) | 一种同源重组修复缺陷的评估方法、装置和存储介质 | |
| Chen et al. | DeBreak: Deciphering the exact breakpoints of structural variations using long sequencing reads | |
| CN117174178A (zh) | 一种基于二代短读长序列的单倍型距离评估方法及装置 | |
| Alsaedi | Evaluating the Application of Allele Frequency in the Saudi Population Variant Detection | |
| CN117953967A (zh) | 一种评估肿瘤纯度的方法、校正tmb的方法及设备 | |
| US20160154930A1 (en) | Methods for identification of individuals |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| TA01 | Transfer of patent application right | ||
| TA01 | Transfer of patent application right |
Effective date of registration: 20190828 Address after: 225300 Taizhou Pharmaceutical High-tech Industrial Park, Taizhou City, Jiangsu Province, Fifth Phase Standard Factory Building G129, 8-storey East and 9-storey Applicant after: Yuze Medical Devices Jiangsu Co., Ltd. Address before: 518081 Shenyan Road, Yantian District, Shenzhen City, Guangdong Province Applicant before: Shenzhen yulce Biological Technology Co., Ltd. |
|
| GR01 | Patent grant | ||
| GR01 | Patent grant |