CN110010197B - 基于血液循环肿瘤dna的单核苷酸变异检测方法、装置和存储介质 - Google Patents

基于血液循环肿瘤dna的单核苷酸变异检测方法、装置和存储介质 Download PDF

Info

Publication number
CN110010197B
CN110010197B CN201910255969.6A CN201910255969A CN110010197B CN 110010197 B CN110010197 B CN 110010197B CN 201910255969 A CN201910255969 A CN 201910255969A CN 110010197 B CN110010197 B CN 110010197B
Authority
CN
China
Prior art keywords
mutation
site
mutation frequency
frequency
background
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
Application number
CN201910255969.6A
Other languages
English (en)
Other versions
CN110010197A (zh
Inventor
倪帅
李淼
陈龙昀
张艳鹏
但旭
陈超
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shenzhen Yuce Biotechnology Co ltd
Original Assignee
Shenzhen Yuce Biotechnology Co ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Shenzhen Yuce Biotechnology Co ltd filed Critical Shenzhen Yuce Biotechnology Co ltd
Priority to CN201910255969.6A priority Critical patent/CN110010197B/zh
Publication of CN110010197A publication Critical patent/CN110010197A/zh
Application granted granted Critical
Publication of CN110010197B publication Critical patent/CN110010197B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/50—Mutagenesis

Landscapes

  • Bioinformatics & Cheminformatics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Genetics & Genomics (AREA)
  • Biotechnology (AREA)
  • Biophysics (AREA)
  • Chemical & Material Sciences (AREA)
  • Molecular Biology (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Analytical Chemistry (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

一种基于血液循环肿瘤DNA的单核苷酸变异检测方法、装置和存储介质，该方法包括：获取测试样本的血液循环肿瘤DNA各位点的突变数据，该突变数据包括位点突变频率；获取训练样本的每个位点背景突变频率的置信范围，该置信范围是通过对每一例训练样本中的所有三碱基突变频率和位点突变频率进行学习建模，并使用原地更新的列表对模型进行训练而得到；对测试样本的各位点的位点突变频率和模型中每个位点的背景突变频率的置信范围进行比较，输出测试样本的位点突变频率未在置信范围内的单核苷酸变异作为检测结果。本发明在大幅优化计算资源需求和检测速度的同时，提升了检测ctDNA单核苷酸突变的敏感性和准确性，满足肿瘤临床检测ctDNA单核苷酸突变可靠性需求。

Description

基于血液循环肿瘤DNA的单核苷酸变异检测方法、装置和存储 介质
技术领域
本发明涉及肿瘤检测技术领域，具体涉及一种基于血液循环肿瘤DNA的单核苷酸变异检测方法、装置和存储介质。
背景技术
循环肿瘤DNA(ctDNA)是指癌细胞死亡时被释放到患者血液中的肿瘤DNA。对ctDNA的分析有助于确定肿瘤的突变类型，同时监测肿瘤的生长。肿瘤来源的DNA可能携带与正常DNA不同的突变，从而被区分开来。然而，有时ctDNA在血液DNA中的含量极少，准确检测出突变DNA，对现有的数据分析方法提出了挑战。
近年来，DNA测序技术飞速发展。以Illumina边合成边测序技术(SBS)为代表的二代测序技术，由于价格较低、准确度较高而成为癌症基因组测序方法的首选。二代测序技术可以同时对基因组中多个区域进行测序，从而在基因层面准确地确定肿瘤的突变类型。可是，由于二代测序过程中的一些步骤如PCR扩增和荧光识别有一定的错误率，这给测序结果带来了一定的非生物来源的变异噪音。为了使癌症患者有机会获得更加精准的治疗，从背景突变噪音中分辨肿瘤来源的突变DNA显得至关重要。
ctDNA在血液DNA中的含量从0.01％到50％不等。在ctDNA含量极低时，ctDNA中携带的突变容易被测序结果的变异噪音所干扰。Aaron M Newman等人发现，二代测序的PCR扩增过程总会引起特定的碱基变异，变异集中表现为鸟嘌呤(G)到胸腺嘧啶(T)和胞嘧啶(C)到胸腺嘧啶(T)的替换。基于此发现，Aaron等人第一次提出了通过学习已知背景变异信息而降低背景突变噪音的模型iDES(integrated digital error suppression)，这是ctDNA测序应用中第一个通过学习正常样本的变异信息来获得背景突变特征的模型。经过iDES的校正，样本中不带突变噪音的位点的比例从90％提高到了98％，很大程度上提高了样本变异检测的灵敏度。
iDES有效地降低了背景变异中的噪音，而Shibing Deng等人对模型进行了更精细的优化，提出了基于学习已知背景变异信息中的连续三碱基的变异率来降低背景突变噪音的模型TNER(Tri-Nucleotide Error Reducer)。他们将单碱基突变扩展为该突变与前后各一个碱基的组合(Tri-nucleotide)，将6种变异类型扩展为96种，发现上述连续三碱基组合的突变频率在某个单碱基突变中的出现频率也并不相同。同时，Shibing Deng等人用二项分布和贝叶斯方法替代了iDES中基于高斯分布对变异次数的描述，使模型更加符合真实数据的表现。
在14组测试数据中，相比于iDES，TNER将背景中不带突变噪音的位点的比例从平均98％提高到了99％，并且将错误率从平均0.03降低到0.02。然而，TNER只适用于少量数据的训练，并没有考虑到在训练样本和测试样本增加时消耗的计算内存和时间。这导致TNER在训练样本的数量快速增加时会占用大量的内存。另外，TNER只在学习正常样本的变异信息时考虑了背景测序数据中测序深度对突变检测可信度的影响，忽略了实际检测过程中同一个样本的测序深度可能不一致的事实。这导致了在实际检测中，样本测序深度低的区域更容易出现假阳性。
单核苷酸多态性(SNP)是指在人群中所占的比例都大于1％的基因组中特定位点的单核苷酸变异，每个正常人的基因组中平均携带400-500万个SNP。TNER和iDES在统计背景碱基突变频率时没有有效地区分样本中背景突变与杂合子SNP的突变频率，导致对背景突变频率的估计偏低，影响变异检测的准确性。另外，在人群SNP位点产生的变异通常会被SNP在一部分训练样本中出现的较高变异频率所覆盖，从而iDES和TNER不能有效检出，这降低了变异检测的灵敏度。
发明内容
本发明提供一种基于血液循环肿瘤DNA的单核苷酸变异检测方法、装置和存储介质，在大幅优化计算资源需求和检测速度的同时，提升了检测ctDNA单核苷酸突变的敏感性和准确性，满足肿瘤临床检测ctDNA单核苷酸突变可靠性需求。
根据第一方面，一种实施例中提供一种基于血液循环肿瘤DNA的单核苷酸变异检测方法，包括：
获取测试样本的血液循环肿瘤DNA各位点的突变数据，上述突变数据包括位点突变频率；
获取训练样本的每个位点背景突变频率的置信范围，该置信范围是通过对每一例训练样本中的所有三碱基突变频率和位点突变频率进行学习建模，并使用原地更新的列表对模型进行训练而得到；
对上述测试样本的各位点的位点突变频率和模型中每个位点的背景突变频率的置信范围进行比较，输出测试样本的位点突变频率未在上述置信范围内的单核苷酸变异作为检测结果。
在优选实施例中，上述测试样本是肿瘤患者血液样本，上述训练样本是正常人血液样本。
在优选实施例中，上述测试样本的位点突变频率是根据上述测试样本的测序深度进行校正得到的校正位点突变频率。
在优选实施例中，上述校正位点突变频率通过如下公式得到：
θadj＝θj×T(min(l，Dj/Dlimit)，α，β)；
其中，θadj为在当前位点的校正位点突变频率，θj为在当前位点实际观测的位点突变频率，τ是以α和β为形状参数的Beta分布的累计概率分布函数，Dj为上述测试样本在当前位点的实际测序深度，Dlimit为设定的最低校正测序深度。
在优选实施例中，上述背景突变频率的置信范围由以下方法确定：
获取一组训练样本的ACGT格式文件，该ACGT格式文件包含目标区域中每个位点的位置、测序深度和突变到任意其它三种非参考碱基的突变频率；
统计所有训练样本中每种三碱基突变的平均突变频率，并将其作为每种三碱基突变的先验突变频率；
从所有训练样本中提取每一位点的突变频率平均值；在每一个位点上，将上述突变频率平均值与上述先验突变频率进行加权，得到该突变位点的加权突变频率，作为位点后验突变频率，该位点后验突变频率符合贝塔分布；
求出所有训练样本在上述目标区域中各个位点的测序深度平均值；
将上述位点后验突变频率和上述位点的测序深度平均值在给定的显著性水平下得出位点后验突变频率在贝塔分布中的置信范围，作为上述背景突变频率的置信范围。
在优选实施例中，上述给定的显著性水平是0.001。
在优选实施例中，上述三碱基突变的先验突变频率由以下方法确定：
获取每个训练样本中的背景突变和SNP突变，其中杂合子和纯合子SNP的突变频率分别在0.5和1处聚集并呈高斯分布，而背景突变的突变频率在0.001-0.1处聚集并呈伽马分布；
通过对上述背景突变和SNP突变的突变频率形成的混合分布进行拟合，找出混合分布的概率密度分布中背景突变与杂合子SNP突变之间概率密度分布的最低点所对应的突变频率，将该突变频率作为背景突变频率的阈值，将突变频率小于该阈值的突变作为真实背景突变；
在上述真实背景突变中对每个训练样本的相同的三碱基突变进行归类后求三碱基平均背景突变频率，然后将所有训练样本中相同的三碱基平均背景突变频率的平均值作为上述三碱基突变的先验突变频率。
在优选实施例中，上述方法还包括：
对上述测试样本的位点突变频率低于上述背景突变频率的置信范围下限的位点进行二次筛选。
在优选实施例中，上述二次筛选包括：
筛选上述测试样本的位点突变频率低于该位点背景突变频率的置信范围下限且高于全局突变频率阈值的单核苷酸变异。
在优选实施例中，上述全局突变频率阈值是0.005。
根据第二方面，一种实施例中提供一种基于血液循环肿瘤DNA的单核苷酸变异检测装置，包括：
测试样本数据获取模块，用于获取测试样本的血液循环肿瘤DNA各位点的突变数据，上述突变数据包括位点突变频率；
置信范围获取模块，用于获取训练样本的每个位点背景突变频率的置信范围，该置信范围是通过对每一例训练样本中的所有三碱基突变频率和位点突变频率进行学习建模，并使用原地更新的列表对模型进行训练而得到；
数据比较和输出模块，用于对上述测试样本的各位点的位点突变频率和模型中每个位点的背景突变频率的置信范围进行比较，输出测试样本的位点突变频率未在上述置信范围内的单核苷酸变异。
在优选实施例中，上述装置还包括：
二次筛选模块，用于筛选上述测试样本的位点突变频率低于该位点背景突变频率的置信范围下限且高于全局突变频率阈值的单核苷酸变异；
在优选实施例中，上述全局突变频率阈值是0.005。
根据第三方面，一种实施例中提供一种计算机可读存储介质，包括程序，该程序能够被处理器执行以实现如第一方面的方法。
本发明的方法对现有模型的训练结构和训练方式进行了优化，使用原地更新的列表来存储每一例训练样本中的所有突变频率，使同样数据下模型消耗的内存降低，使每一例新增样本消耗内存减少，在完成训练后直接计算并保存每个位点的背景突变频率置信范围，检测新的测试样本时，不需要重复计算这些背景突变频率的置信范围，使得在实际运行中，每个测试样本的检测时间大大降低。
此外，在优选实施例中，增加了根据测试样本的测序深度对测试样本的位点突变频率进行校正的步骤，降低了低测序深度区域的突变检测的假阳性；通过统计学习每个训练样本中背景突变频率与单核苷酸多态性位点突变频率的差异，重新界定选取三碱基突变频率动态阈值，增加了对背景突变频率估计的精确程度；本发明方法增加了对测试样本的位点突变频率低于背景突变频率置信范围下限的位点的二次筛选，提高了在单核苷酸多态性位点上的突变检测的灵敏度。
附图说明
图1为本发明实施例中一种基于血液循环肿瘤DNA的单核苷酸变异检测方法流程图。
图2为本发明实施例中一个真实训练数据样例中突变频率的密度分布图，其中实线代表动态界定三碱基突变频率的选取位置，虚线代表TNER方法中固定对三碱基突变频率胚系SNP突变频率的分割位置，动态地选取三碱基突变频率与胚系SNP突变频率的分割位置能够更精确地记录背景突变频率。
图3本发明实施例中优化后的RAM使用情况比较结果图，其中x轴表示训练样本的数量，y轴表示程序所占用的内存，单位为Mb，虚线为TNER方法，实线为本发明方法。
图4为本发明实施例的实际测试中两种方法(本发明方法(Optimized)和TNER方法(Original))对于每新增一个相同的训练样本时程序内存消耗的增加量比较结果图，其中x轴表示两种方法，y轴表示增加的内存，单位为Mb。
图5为本发明实施例的实际测试中样本SNV检测所需时间和节省的时间比较结果图，其中x轴表示样本，y轴表示时间，单位为秒，实线代表TNER方法在一个样本的SNV检测所需时间，虚线代表本发明方法可以节省的时间。
图6为本发明实施例中一个实际样例用两种方法检测到的变异比较结果图，其中x轴左侧为所有检测到的突变的突变频率，右侧为左侧对应突变的背景突变频率，实线连接的突变为两种方法(本发明方法和TNER方法)同时检测到的突变，虚线连接的突变为只在本发明方法中检测到的突变，可以看出本发明方法允许在背景变异频率高于样本变异频率的位点中检测单核苷酸突变，而图中背景变异频率高的位点中有75％是已知的人群SNP位点。
图7为本发明实施例中一种基于血液循环肿瘤DNA的单核苷酸变异检测装置结构框图。
图8为本发明一个实施例中比较两种方法(本发明方法和TNER方法)在189例癌症患者血液样本检测到的所有单核苷酸变异对应位点的测序深度，TNER方法(灰色)在低深度区域会富集更多突变，本发明方法拒绝了大部分低深度区域的单碱基突变，接受了更多高深度区域的单碱基突变。
具体实施方式
下面通过具体实施方式结合附图对本发明作进一步详细说明。在以下的实施方式中，很多细节描述是为了使得本发明能被更好的理解。然而，本领域技术人员可以毫不费力的认识到，其中部分特征在不同情况下是可以省略的，或者可以由其他材料、方法所替代。
另外，说明书中所描述的特点、操作或者特征可以以任意适当的方式结合形成各种实施方式。同时，方法描述中的各步骤或者动作也可以按照本领域技术人员所能显而易见的方式进行顺序调换或调整。因此，说明书和附图中的各种顺序只是为了清楚描述某一个实施例，并不意味着是必须的顺序，除非另有说明其中某个顺序是必须遵循的。
本发明中采用的术语具体含义如下：
参考基因组：物种参考的标准基因组序列。
读长(Reads)：测序所得基因组序列片段。
BAM：一种用于存储比对信息的标准二进制文件格式。
acgt：一种记录每个位点单核苷酸变异信息的文件。
Indel：一种记录每个位点插入和缺失类型变异信息的文件。
三碱基突变(mutational signature)：是指6种基础单碱基突变形式(A→T、A→G、A→C、C→A、C→T、C→G)与其上下文各一个碱基的组合，共有96种。
如图1所示，本发明的一种实施例中提供一种基于血液循环肿瘤DNA的单核苷酸变异检测方法，包括：
S101：获取测试样本的血液循环肿瘤DNA各位点的突变数据，上述突变数据包括位点突变频率；
S102：获取训练样本的每个位点背景突变频率的置信范围，该置信范围是通过对每一例训练样本中的所有三碱基突变频率和位点突变频率进行学习建模，并使用原地更新的列表对模型进行训练而得到；和
S103：对上述测试样本的各位点的位点突变频率和模型中每个位点的背景突变频率的置信范围进行比较，输出测试样本的位点突变频率未在上述置信范围内的单核苷酸变异作为检测结果。
首先，本发明的方法，对现有模型的训练结构进行了优化，分为训练阶段和测试阶段。在训练阶段，使用原地更新的列表来存储每一例训练样本(正常血液样本)中的所有突变频率，并且在完成训练后直接计算并保存每个位点的背景突变频率的置信范围，在测试样本的实际检测过程中，直接对测试样本的每个位点变异频率与相应位点的背景突变频率的置信范围进行比较。
现有方法中，假设用100个训练样本的测序数据对模型进行训练，在现有的TNER方法(Shibing Deng,Maruja Lira,Donghui Huang,Kai Wang,Crystal Valdez,JenniferKinong,Paul A.Rejto,Jadwiga Bienkowska,James Hardwick,Tao Xie.TNER:A NovelBayesian Background Error Suppression Method for Mutation Detection inCirculating Tumor DNA,BMC Bioinformatics,(2018)19:387)中，100个训练样本的选定的目标区域中的每个位点对应到其它三种非参考碱基的突变频率的矩阵会被逐一读入内存，原样汇编成一个大的数据集，然后存储在硬盘中，等待测试过程中读取并处理。
本发明中，一个结构优化的实例是，100个训练样本的数据被逐一读入内存，用一个矩阵P记录下各位点的突变频率。在下一个训练样本的数据被读入之前，释放前一个训练样本的数据所占用的内存，将下一个训练样本的数据对应位点的突变频率累加在相同的矩阵P中。最后用矩阵P各位点的累加突变频率除以训练样本的数据个数(这里是100)得到平均突变频率。
其次，本发明的方法，增加了根据测试样本的测序深度对测试样本的位点突变频率进行校正的步骤，例如，用贝塔(Beta)分布的累计分布函数作为惩罚函数对测试样本中测序深度低于某一阈值的位点的突变频率进行校正，得到校正位点突变频率。
例如，在一个具体实施例中，校正位点突变频率通过如下公式得到：
θadj＝θj×T(min(l，Dj/Dlimit)，α，β)；
其中，θadj为在当前位点的校正位点突变频率，θj为在当前位点实际观测的位点突变频率，τ是以α和β为形状参数的Beta分布的累计概率分布函数，Dj为上述测试样本在当前位点的实际测序深度，Dlimit为设定的最低校正测序深度。
需要对测试样本的位点突变频率进行校正的依据是，假设选定的目标区域测序后，某个区域A深度很低，例如只有8X，在此区域中的1条测序序列的某个位点发生1个突变，这个突变的突变频率就为1/8＝0.125。同时，某个区域B深度很高，达到8000X，在此区域中的1000条测序序列的某个位点发生突变，这个突变的突变频率同样也为1/8＝0.125。
然而，可以看出，区域B对应位点的突变频率更加可信，因为每个序列被捕获测序是相对随机的过程。在B区域的突变位点增加1条突变序列，对总体突变频率影响仅为0.0125％，而如果区域A中增加1条突变序列，对总体突变频率影响则为12.5％。说明超低深度的突变频率可信度不高。
为了降低低深度区域的类似突变频率高的噪音，需要对低深度区域的突变频率乘以一个系数Q。深度越低，对应的Q越小，比如8X的区间突变频率本来为12.5％，乘以小数0.1后，突变频率变为1.25％；20X的区间某突变的突变频率本来为5％，乘以小数0.5后，突变频率变为2.5％。降低了深度极低的区域突变频率可信度不高的问题。
再次，在本发明一个实施例中，背景突变频率的置信范围由以下方法确定：获取一组训练样本的ACGT格式文件，该ACGT格式文件包含目标区域中每个位点的位置、测序深度和突变到任意其它三种非参考碱基的突变频率；统计所有训练样本中每种三碱基突变的平均突变频率，并将其作为每种三碱基突变的先验突变频率；从所有训练样本中提取每一位点的突变频率平均值；在每一个位点上，将上述突变频率平均值与上述先验突变频率进行加权，得到该突变位点的加权突变频率，作为位点后验突变频率，该位点后验突变频率符合贝塔分布；求出所有训练样本在上述目标区域中各个位点的测序深度平均值；将上述位点后验突变频率和上述位点的测序深度平均值在给定的显著性水平(例如0.001)下得出位点后验突变频率在贝塔分布中的置信范围，作为上述背景突变频率的置信范围。
本发明的方法，通过学习训练样本的数据中三碱基背景突变频率与单核苷酸多态性位点突变频率的差异，动态界定训练样本的数据中三碱基突变频率的选取阈值。
具体而言，在本发明一个实施例中，通过以下方法确定三碱基突变的先验突变频率：获取每个训练样本中的背景突变和SNP突变，其中杂合子和纯合子SNP的突变频率分别在0.5和1处聚集并呈高斯分布，而背景突变的突变频率在0.001-0.1处聚集并呈伽马分布；通过对上述背景突变和SNP突变的突变频率形成的混合分布进行拟合，找出混合分布的概率密度分布中背景突变与杂合子SNP突变之间概率密度分布的最低点所对应的突变频率，将该突变频率作为背景突变频率的阈值，将突变频率小于该阈值的突变作为真实背景突变；在上述真实背景突变中对每个训练样本的相同的三碱基突变进行归类后求三碱基平均背景突变频率，然后将所有训练样本中相同的三碱基平均背景突变频率的平均值作为上述三碱基突变的先验突变频率。
为了更清楚的说明这一点，图2示出了一个真实训练数据样例中突变频率的密度分布图。可以看出，绝大多数突变的突变频率集中在-3附近(对应突变频率为0.001左右)。右边两个峰代表突变频率在0.5和1时的SNP的密度分布。两个峰之外的其它所有突变理论上都是背景突变。
要获取到所有的背景突变，才能够正确地估计平均的三碱基背景突变频率，从而分辨背景突变和肿瘤中的真实突变。因此，获取所有的背景突变非常重要。虚线代表现有的TNER方法中固定对三碱基突变频率胚系SNP突变频率的分割位置，固定取值0.1。这样会漏判一些突变频率高于0.1的背景突变。而实线代表动态界定三碱基突变频率的选取位置，动态地选取三碱基突变频率与胚系SNP突变频率的分割位置，更精确地识别所有的背景突变所对应的三碱基突变频率。
最后，本发明的方法，对测试样本的位点突变频率低于背景突变频率的置信范围下限的位点进行二次筛选。具体而言，在一个实施例中，筛选测试样本的位点突变频率低于该位点背景突变频率的置信范围下限且高于全局突变频率阈值(例如0.005)的单核苷酸变异。这对检测人群SNP在单个患者中发生的突变有较好的检测效果。
为了更清楚的说明这一点，结合现有技术解释如下：现有的TNER方法对于突变的筛选法则比较简单，突变频率大于阈值K，则判断为肿瘤细胞中的真实突变。在计算突变背景阈值K时，联合了三碱基突变频率T和该位点的突变频率S。但是在训练阶段，位点的突变频率S会被人群携带SNP的比例所影响。例如，在100个样本中计算位点突变频率S时，如果其中45个样本在相同的位点带有纯合子SNP(A→T)，那么根据TNER算法的统计方法，该位点在100个样本中的位点(A→T)突变频率就大约为45％，结合三碱基突变频率T后，阈值K在99.99％置信度下的置信区间的置信上限确定为49％。这样的阈值无法检测出肿瘤样品中突变频率低于49％的真实突变。
本发明的方法，经过优化以后，对位点背景阈值K在99.99％置信度下的对应的最小置信区间也做了界定(例如43％)。如果在某例测试样本ctDNA中该突变(A→T)频率为5％，虽然超过了芯片中所有位点的平均变异频率(假设为0.67％)，但是小于当前位点阈值K的置信区间的最大值。会继续比较阈值K所对应的置信区间的最小值，如果阈值K最小值也无法解释观测到的变异频率，那么输出该变异作为SNP位点发生的体细胞变异。
本发明的方法对现有模型的训练结构和训练方式进行了优化，使用原地更新的列表来存储每一例训练样本中的所有突变频率，使同样数据下模型消耗的内存(RAM)降低(图3)，使每一例新增样本消耗内存减少2/3以上(图4)，在完成训练后直接计算并保存每个位点的背景突变频率置信范围，检测新的测试样本时，不需要重复计算这些背景变异频率的置信范围，使得在实际运行中，对每个测试样本的检测时间节省约80％(图5)。
此外，在优选实施例中，增加了根据测试样本的测序深度对测试样本的位点突变频率进行校正的步骤，降低了低测序深度区域的突变检测的假阳性；通过统计学习每个训练样本中背景突变频率与单核苷酸多态性位点突变频率的差异，重新界定选取动态阈值，增加了对背景突变频率估计的精确程度(图2)。由于在单核苷酸多态性位点上的背景突变频率较高，在大多数情况下会覆盖真实的低频突变，本发明的方法增加了对测试样本的位点突变频率低于背景变异频率置信下限的位点的二次筛选，提高了在单核苷酸多态性位点上的突变检测的灵敏度(图6)。
本领域技术人员可以理解，上述实施方式中各种方法的全部或部分功能可以通过硬件的方式实现，也可以通过计算机程序的方式实现。当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘、光盘、硬盘等，通过计算机执行该程序以实现上述功能。例如，将程序存储在设备的存储器中，当通过处理器执行存储器中程序，即可实现上述全部或部分功能。另外，当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序也可以存储在服务器、另一计算机、磁盘、光盘、闪存盘或移动硬盘等存储介质中，通过下载或复制保存到本地设备的存储器中，或对本地设备的系统进行版本更新，当通过处理器执行存储器中的程序时，即可实现上述实施方式中全部或部分功能。
因此，对应于本发明的方法，本发明一种实施例中提供一种基于血液循环肿瘤DNA的单核苷酸变异检测装置，如图7所示，包括：测试样本数据获取模块701，用于获取测试样本的血液循环肿瘤DNA各位点的突变数据，上述突变数据包括位点突变频率；置信范围获取模块702，用于获取训练样本的每个位点背景突变频率的置信范围，该置信范围是通过对每一例训练样本中的所有三碱基突变频率和位点突变频率进行学习建模，并使用原地更新的列表对模型进行训练而得到；数据比较和输出模块703，用于对上述测试样本的各位点的位点突变频率和模型中每个位点的背景突变频率的置信范围进行比较，输出测试样本的位点突变频率未在上述置信范围内的单核苷酸变异。
此外，本发明的一种实施例中提供一种计算机可读存储介质，包括程序，该程序能够被处理器执行以实现如本发明的基于血液循环肿瘤DNA的单核苷酸变异检测方法。
以下通过实施例详细说明本发明的技术方案，应当理解，实施例仅是示例性的，不能理解为对本发明保护范围的限制。
实施例1
本实施例中，训练数据：10例健康人群血液样本经过深圳裕策生物科技有限公司的靶向(Target)芯片测序生成的ACGT文件；测试数据：189例癌症患者血液样本经过深圳裕策生物科技有限公司的靶向(Target)芯片测序生成的ACGT文件。
将10例健康人群血液样本，经过深圳裕策生物科技有限公司的靶向(Target)芯片进行测序后所得基因组序列片段与人参考基因组进行比对，得到BAM格式的比对结果。然后对BAM格式文件用Samtools软件转化为pileup格式文件。在转化过程中，只允许测序错误和比对错误率小于0.1％的reads，对应的Phread得分(Phread Score)和映射得分(MappingScore)均为30。再将生成的pileup格式文件用sequenza-utils以默认参数转化为单核苷酸突变频率数据文件ACGT格式。将该ACGT文件作为置信范围获取模块的输入数据，根据置信范围的确定方法，获取训练样本的每个位点背景突变频率的置信范围。
189例癌症患者血液样本，经过深圳裕策生物科技有限公司的靶向(Target)芯片进行测序后所得基因组序列片段与人参考基因组进行比对，得到BAM格式的比对结果。然后对BAM格式文件用Samtools软件转化为pileup格式文件。在转化过程中，只允许测序错误和比对错误率小于0.1％的reads，对应的Phread得分(Phread Score)和映射得分(MappingScore)均为30。再将生成的pileup格式文件用sequenza-utils以默认参数转化为单核苷酸突变频率数据文件ACGT格式。将该ACGT文件用本发明方法进行单核苷酸变异检测。
最终得到的189例癌症患者血液样本中编号为18A01602XJ03的样本中的突变检测结果(图6)和所有189例癌症患者血液样本中检测到的单核苷酸变异突变频率与TNER方法检测到的单核苷酸变异突变频率的分布进行比较(图8)。
在编号为18A01602XJ03样本中，本发明方法通过对该样本的位点突变频率低于背景突变频率的置信范围下限的位点进行二次筛选，成功地独立检测出8个背景变异频率高于样本变异频率的单核苷酸变异，而TNER方法没有检测出。其中6个为已知的SNP，并且均在于癌症高度相关的基因上(表1)。BRCA2是关键的DNA错配修复基因，BRCA2基因突变可能会增加乳腺癌、卵巢癌的患病风险。EGFR基因编码表皮细胞生长因子受体，通常会在癌细胞的表面大量出现，促进癌细胞生长。KDR基因编码血管内皮生长因子受体-2，在部分肿瘤中高表达，促进血管内皮细胞分裂、增殖，诱导肿瘤血管增生；PTCH也作为抑癌基因被研究。可以发现，本方法能够在一些背景变异频率高的关键SNP位点上检测单核苷酸变异，提高了单核苷酸变异检测的灵敏度。
表1
变异位点 dbSNP编码 变异所在的基因
Chr13:32890572 rs1799943 BRCA2基因
Chr13:32911888 rs1801406 BRCA2基因
Chr13:32929232 rs1799955 BRCA2基因
chr4:55961159 rs2219471 KDR基因
chr7:55214348 rs2017454 EGFR基因
chr9:98229389 rs2066829 PTCH1基因
比较所有189例癌症患者血液样本检测到的单核苷酸变异对应位点的测序深度的密度分布可以发现，TNER方法在低深度区域富集了更多突变，而低深度区域的测序结果可信度偏低，增加了单核苷酸变异检测假阳性的风险。本发明方法拒绝了大部分低深度区域的单碱基突变，接受了更多高深度区域的单碱基突变(图8)。在降低了假阳性风险的同时，提高了检测的灵敏度。
综合以上实施例，本发明方法克服了现有技术难以在训练数据中背景突变频率偏高的区域检测单核苷酸突变的问题，提高了单核苷酸变异检测的灵敏度。对训练数据中背景突变频率可信度偏低的区域的突变频率进行了更精确地估计，提高了单核苷酸突变检测的特异性。
以上应用了具体个例对本发明进行阐述，只是用于帮助理解本发明，并不用以限制本发明。对于本发明所属技术领域的技术人员，依据本发明的思想，还可以做出若干简单推演、变形或替换。

Claims (12)

1.一种基于血液循环肿瘤DNA的单核苷酸变异检测方法，其特征在于，所述方法包括：
获取测试样本的血液循环肿瘤DNA各位点的突变数据，所述突变数据包括位点突变频率；其中，所述测试样本的位点突变频率是根据所述测试样本的测序深度进行校正得到的校正位点突变频率，所述校正位点突变频率通过如下公式得到：
Figure 612157DEST_PATH_IMAGE001
，
其中，θ adj 为在当前位点的校正位点突变频率，θ j 为在当前位点实际观测的位点突变频率，τ是以α和β为形状参数的Beta分布的累计概率分布函数，Dj为所述测试样本在当前位点的实际测序深度，Dlimit为设定的最低校正测序深度；
获取训练样本的每个位点背景突变频率的置信范围，该置信范围是通过对每一例训练样本中的所有三碱基突变频率和位点突变频率进行学习建模，并使用原地更新的列表对模型进行训练而得到；
对所述测试样本的各位点的位点突变频率和模型中每个位点的背景突变频率的置信范围进行比较，输出测试样本的位点突变频率未在所述置信范围内的单核苷酸变异作为检测结果。
2.根据权利要求1所述的单核苷酸变异检测方法，其特征在于，所述测试样本是肿瘤患者血液样本，所述训练样本是正常人血液样本。
3.根据权利要求1所述的单核苷酸变异检测方法，其特征在于，所述背景突变频率的置信范围由以下方法确定：
获取一组训练样本的ACGT格式文件，该ACGT格式文件包含目标区域中每个位点的位置、测序深度和突变到任意其它三种非参考碱基的突变频率；
统计所有训练样本中每种三碱基突变的平均突变频率，并将其作为每种三碱基突变的先验突变频率；
从所有训练样本中提取每一位点的突变频率平均值；在每一个位点上，将所述突变频率平均值与所述先验突变频率进行加权，得到该突变位点的加权突变频率，作为位点后验突变频率，该位点后验突变频率符合贝塔分布；
求出所有训练样本在所述目标区域中各个位点的测序深度平均值；
将所述位点后验突变频率和所述位点的测序深度平均值在给定的显著性水平下得出位点后验突变频率在贝塔分布中的置信范围，作为所述背景突变频率的置信范围。
4.根据权利要求3所述的单核苷酸变异检测方法，其特征在于，所述给定的显著性水平是0.001。
5.根据权利要求3所述的单核苷酸变异检测方法，其特征在于，所述三碱基突变的先验突变频率由以下方法确定：
获取每个训练样本中的背景突变和SNP突变，其中杂合子和纯合子SNP的突变频率分别在0.5和1处聚集并呈高斯分布，而背景突变的突变频率在0.001-0.1处聚集并呈伽马分布；
通过对所述背景突变和SNP突变的突变频率形成的混合分布进行拟合，找出混合分布的概率密度分布中背景突变与杂合子SNP突变之间概率密度分布的最低点所对应的突变频率，将该突变频率作为背景突变频率的阈值，将突变频率小于该阈值的突变作为真实背景突变；
在所述真实背景突变中对每个训练样本的相同的三碱基突变进行归类后求三碱基平均背景突变频率，然后将所有训练样本中相同的三碱基平均背景突变频率的平均值作为所述三碱基突变的先验突变频率。
6.根据权利要求1所述的单核苷酸变异检测方法，其特征在于，所述方法还包括：
对所述测试样本的位点突变频率低于所述背景突变频率的置信范围下限的位点进行二次筛选。
7.根据权利要求6所述的单核苷酸变异检测方法，其特征在于，所述二次筛选包括：
筛选所述测试样本的位点突变频率低于该位点背景突变频率的置信范围下限且高于全局突变频率阈值的单核苷酸变异。
8.根据权利要求7所述的单核苷酸变异检测方法，其特征在于，所述全局突变频率阈值是0.005。
9.一种基于血液循环肿瘤DNA的单核苷酸变异检测装置，其特征在于，所述装置包括：
测试样本数据获取模块，用于获取测试样本的血液循环肿瘤DNA各位点的突变数据，所述突变数据包括位点突变频率；其中，所述测试样本的位点突变频率是根据所述测试样本的测序深度进行校正得到的校正位点突变频率，所述校正位点突变频率通过如下公式得到：
Figure 945050DEST_PATH_IMAGE001
，
其中，θ adj 为在当前位点的校正位点突变频率，θ j 为在当前位点实际观测的位点突变频率，τ是以α和β为形状参数的Beta分布的累计概率分布函数，Dj为所述测试样本在当前位点的实际测序深度，Dlimit为设定的最低校正测序深度；
置信范围获取模块，用于获取训练样本的每个位点背景突变频率的置信范围，该置信范围是通过对每一例训练样本中的所有三碱基突变频率和位点突变频率进行学习建模，并使用原地更新的列表对模型进行训练而得到；
数据比较和输出模块，用于对所述测试样本的各位点的位点突变频率和模型中每个位点的背景突变频率的置信范围进行比较，输出测试样本的位点突变频率未在所述置信范围内的单核苷酸变异。
10.根据权利要求9所述的单核苷酸变异检测装置，其特征在于，所述装置还包括：
二次筛选模块，用于筛选所述测试样本的位点突变频率低于该位点背景突变频率的置信范围下限且高于全局突变频率阈值的单核苷酸变异。
11.根据权利要求10所述的单核苷酸变异检测装置，其特征在于，所述全局突变频率阈值是0.005。
12.一种计算机可读存储介质，其特征在于，包括程序，所述程序能够被处理器执行以实现如权利要求1-8中任一项所述的单核苷酸变异检测方法。
CN201910255969.6A 2019-03-29 2019-03-29 基于血液循环肿瘤dna的单核苷酸变异检测方法、装置和存储介质 Active CN110010197B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910255969.6A CN110010197B (zh) 2019-03-29 2019-03-29 基于血液循环肿瘤dna的单核苷酸变异检测方法、装置和存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910255969.6A CN110010197B (zh) 2019-03-29 2019-03-29 基于血液循环肿瘤dna的单核苷酸变异检测方法、装置和存储介质

Publications (2)

Publication Number Publication Date
CN110010197A CN110010197A (zh) 2019-07-12
CN110010197B true CN110010197B (zh) 2021-07-20

Family

ID=67169321

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910255969.6A Active CN110010197B (zh) 2019-03-29 2019-03-29 基于血液循环肿瘤dna的单核苷酸变异检测方法、装置和存储介质

Country Status (1)

Country Link
CN (1) CN110010197B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110867207B (zh) * 2019-11-26 2021-07-30 北京橡鑫生物科技有限公司 验证ngs变异检测方法的评估方法及评估装置
CN113186255A (zh) * 2021-05-12 2021-07-30 深圳思勤医疗科技有限公司 基于单分子测序检测核苷酸变异方法与装置
WO2023284260A1 (zh) * 2021-07-12 2023-01-19 广州燃石医学检验所有限公司 基于血液测序的肿瘤内异质性的评估方法及其用于预测免疫疗法的应答
CN114242158B (zh) * 2022-02-21 2022-05-13 臻和(北京)生物科技有限公司 ctDNA单核苷酸变异位点检测方法、装置、存储介质及设备
CN115410649B (zh) * 2022-04-01 2023-03-28 北京吉因加医学检验实验室有限公司 一种同时检测甲基化和突变信息的方法及装置
CN115440299A (zh) * 2022-08-25 2022-12-06 中国科学院心理研究所 确定背景微生物的方法、设备、介质和程序产品
CN115424664B (zh) * 2022-11-07 2023-03-10 北京雅康博生物科技有限公司 人为突变程度评估方法及装置
CN116356001B (zh) * 2023-02-07 2023-12-15 江苏先声医学诊断有限公司 一种基于血液循环肿瘤dna的双重背景噪声突变去除方法

Family Cites Families (19)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2013166517A1 (en) * 2012-05-04 2013-11-07 Complete Genomics, Inc. Methods for determining absolute genome-wide copy number variations of complex tumors
KR101591475B1 (ko) * 2014-07-16 2016-02-03 고려대학교 산학협력단 레일리 산란을 이용한 순환 종양 DNA(ctDNA)의 종양 특이적 돌연변이와 후성유전학적 변이의 동시 검출방법
WO2016090584A1 (zh) * 2014-12-10 2016-06-16 深圳华大基因研究院 确定肿瘤核酸浓度的方法和装置
CN105734122A (zh) * 2014-12-31 2016-07-06 深圳市作物分子设计育种研究院 Simm法快速定位突变性状相关基因
CN105063208B (zh) * 2015-08-10 2018-03-06 北京吉因加科技有限公司 一种血浆中游离的目标dna低频突变富集测序方法
US10364468B2 (en) * 2016-01-13 2019-07-30 Seven Bridges Genomics Inc. Systems and methods for analyzing circulating tumor DNA
CN106021994B (zh) * 2016-05-13 2019-03-26 万康源(天津)基因科技有限公司 一种肿瘤突变位点筛选及互斥基因挖掘的方法
CN106022001B (zh) * 2016-05-13 2018-09-18 万康源(天津)基因科技有限公司 一种肿瘤突变位点筛选及互斥基因挖掘的系统
CN106650312B (zh) * 2016-12-29 2022-05-17 浙江安诺优达生物科技有限公司 一种用于循环肿瘤dna拷贝数变异检测的装置
CN106778073B (zh) * 2017-01-19 2019-09-06 北京吉因加科技有限公司 一种评估肿瘤负荷变化的方法和系统
CN108517360A (zh) * 2017-02-27 2018-09-11 北京医院 一种循环肿瘤游离dna突变检测质控品及其制备方法
CN107423578B (zh) * 2017-03-02 2020-09-22 北京诺禾致源科技股份有限公司 检测体细胞突变的装置
WO2018204657A1 (en) * 2017-05-04 2018-11-08 The Johns Hopkins University Detection of cancer
CN108154010B (zh) * 2017-12-26 2018-10-19 东莞博奥木华基因科技有限公司 一种ctDNA低频突变测序数据分析方法和装置
CN108733975B (zh) * 2018-03-29 2021-09-07 深圳裕策生物科技有限公司 基于二代测序的肿瘤克隆变异检测方法、装置和存储介质
CN108875302B (zh) * 2018-06-22 2022-02-22 广州漫瑞生物信息技术有限公司 一种检测细胞游离肿瘤基因拷贝数变异的系统和方法
CN109033749B (zh) * 2018-06-29 2020-01-14 裕策医疗器械江苏有限公司 一种肿瘤突变负荷检测方法、装置和存储介质
CN109022553B (zh) * 2018-06-29 2019-10-25 裕策医疗器械江苏有限公司 用于肿瘤突变负荷检测的基因芯片及其制备方法和装置
CN109411015B (zh) * 2018-09-28 2020-12-22 深圳裕策生物科技有限公司 基于循环肿瘤dna的肿瘤突变负荷检测装置及存储介质

Also Published As

Publication number Publication date
CN110010197A (zh) 2019-07-12

Similar Documents

Publication Publication Date Title
CN110010197B (zh) 基于血液循环肿瘤dna的单核苷酸变异检测方法、装置和存储介质
US8725422B2 (en) Methods for estimating genome-wide copy number variations
TWI814753B (zh) 用於標靶定序之模型
US20220130488A1 (en) Methods for detecting copy-number variations in next-generation sequencing
CN109949861B (zh) 肿瘤突变负荷检测方法、装置和存储介质
WO2019204360A1 (en) Systems and methods for determining tumor fraction in cell-free nucleic acid
KR20220034923A (ko) 서열-특정 오류(sse)를 유발시키는 서열 패턴을 식별하기 위한 심층 학습-기반 프레임워크
WO2019046804A1 (en) IDENTIFICATION OF FALSE POSITIVE VARIANTS USING A MODEL OF IMPORTANCE
CN115083529A (zh) 一种检测样本污染率的方法及装置
CN114530199A (zh) 基于双重测序数据检测低频突变的方法、装置及存储介质
Duan et al. Common copy number variation detection from multiple sequenced samples
CN115631790A (zh) 单细胞转录组测序数据的体细胞突变提取方法及装置
US20200105374A1 (en) Mixture model for targeted sequencing
O’Fallon et al. Algorithmic improvements for discovery of germline copy number variants in next-generation sequencing data
Huang et al. Robust analysis of allele-specific copy number variations from scRNA-seq data with XClone
CN114974416B (zh) 一种检测相邻多核苷酸变异的方法及装置
CN116434830B (zh) 基于ctDNA多位点甲基化的肿瘤病灶位置识别方法
CN115762641B (zh) 一种指纹图谱构建方法及系统
CN116153394A (zh) 检测snv的装置和方法
O’Fallon et al. Algorithmic improvements for discovery of germline copy number variants in next-generation sequencing data
CN114708905A (zh) 基于ngs的染色体非整倍体检测方法、装置、介质和设备
CN115713107A (zh) 用于变体识别的神经网络
Zeng et al. NONNEGATIVE LEAST SQUARE–A NEW LOOK INTO SAGE DATA
CN117941002A (zh) 染色体和亚染色体拷贝数变异检测
CN115966259A (zh) 一种基于逻辑回归建模的样本同源性检测校验方法及系统

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
GR01 Patent grant
GR01 Patent grant