发明内容
本发明所要解决的技术问题是如何准确鉴定或检测ctDNA测序数据中的突变检测结果为真实突变或测序噪音。
为了解决上述技术问题,本发明首先提供了检测ctDNA测序数据中噪音水平的装置,所述装置可包括如下模块:
A1)ctDNA测序数据处理模块:用于获得健康受试者的ctDNA测序数据,经过质控与人类参考基因组进行比对得到比对文件;
A2)突变频率检测模块:用于基于所述比对文件,获得所述ctDNA测序数据中每个位点的测序深度、每个位点的点突变reads支持数、每个位点的插入或缺失突变reads支持数,基于所述每个位点的测序深度和所述每个位点的点突变reads支持数计算所述每个位点的点突变频率,基于所述每个位点的测序深度和所述每个位点的插入或缺失突变reads支持数计算所述每个位点的插入或缺失突变频率;
A3)点突变噪音水平模型构建模块:用于基于每个位点的点突变特征分组信息和所述每个位点的点突变频率获得所述每个位点的点突变噪音水平;所述每个位点的点突变特征分组信息包括每个位点的三核苷酸序列和每个位点的比对次数;
A4)插入或缺失突变噪音水平模型构建模块:用于基于每个位点的插入或缺失突变特征分组信息和所述每个位点的插入或缺失突变频率获得所述每个位点的插入或缺失突变噪音水平;所述每个位点的插入或缺失突变特征分组信息包括每个位点的重复序列元件信息和每个位点的染色体片段重复信息;
A5)噪音水平输出模块:用于输出所述ctDNA测序数据的噪音水平;所述噪音水平包括所述每个位点的点突变(SNV)噪音水平和/或所述每个位点的插入或缺失突变(InDel)噪音水平。
上文所述位点可为一个核苷酸,具体可为测序产品Panel范围内包含的位点。所述每个位点的三核苷酸序列包括所述每个位点在参考基因组上的核苷酸及其上下游各一个核苷酸组成的三个核苷酸序列。所述每个位点的三核苷酸序列可由人类参考基因组文件获得。所述每个位点的比对次数可为唯一比对或多次比对。所述每个位点的比对次数可由人类参考基因组的Mappability文件获得。所述每个位点的重复序列元件信息表示所述每个位点是否属于重复序列区域以及属于具体哪个重复区域类型。所述每个位点的重复序列元件信息可由人类基因组的重复序列元件数据文件获得。所述每个位点的染色体片段重复信息表示所述每个位点是否属于大片段重复区域。所述每个位点的染色体片段重复信息可由人类基因组的染色体片段重复数据文件获得。
上述装置中,A3)所述点突变噪音水平模型构建模块可包括如下模块:
A3-1)点突变特征整合分组模块:用于基于所述每个位点的三核苷酸序列和所述每个位点的比对次数,将具有相同三核苷酸序列且具有相同比对次数的位点分为一组得到点突变特征相同分组;
A3-2)点突变先验背景噪音水平获得模块:用于基于所述点突变特征相同分组的每个位点的点突变频率和A3-1)所述点突变特征相同分组使用式1计算获得所述点突变特征相同分组的平均点突变频率:
式1中:λ tm 表示所述点突变特征相同分组内所有位点的平均点突变频率,μ ptm 代表A3-1)所述点突变特征相同分组内每个位点的点突变频率,n代表A3-1)所述点突变特征相同分组内位点的数量;
A3-3)点突变噪音水平获得模块:用于基于所述每个位点所属点突变特征相同分组的平均点突变频率和所述每个位点的点突变频率使用式2和式3获得所述每个位点的点突变噪音水平:
式2和式3中:p表示测序产品Panel范围内某个位点的信息,λ tm 表示位点p所属点突变特征相同分组的所有位点的平均点突变频率;μ p 表示位点p的点突变频率;ω 1 代表权重;ε p 表示位点p的点突变噪音水平。
上述装置中,所述Mappability文件可来源于UCSC网站(网址:https://genome.ucsc.edu/cgi-bin/hgTables/hgsid=1510554375_HANSVqcS3xxfWO2Ui7wYz8FFIlzA,参考基因组版本:Hg19;100bp reads长度比对版本)。
上述装置中,A4)所述插入或缺失突变噪音水平模型构建模块可包括如下模块:
A4-1)插入或缺失突变特征整合分组模块:用于基于所述每个位点的重复序列元件信息和所述每个位点的染色体片段重复信息,将具有相同的重复序列元件信息且具有相同染色体片段重复信息的位点分为一组得到插入或缺失突变特征相同分组;
A4-2)插入或缺失突变先验背景噪音水平获得模块:用于基于所述每个位点的插入或缺失突变频率和A4-1)所述插入或缺失突变特征相同分组使用式4计算获得所述插入或缺失突变特征相同分组的平均插入或缺失突变频率:
式4中:λ rsl 表示所述插入或缺失突变特征相同分组内所有位点的平均插入或缺失突变频率,μ prsl 代表A4-1)所述插入或缺失突变特征相同分组内每个位点的插入或缺失突变频率,m代表A4-1)所述插入或缺失突变特征相同分组内位点的数量;
A4-3)插入或缺失突变噪音水平获得模块:用于基于所述每个位点所属插入或缺失突变特征相同分组的平均插入或缺失突变频率和所述每个位点的插入或缺失突变频率使用式5和式6获得所述每个位点的插入或缺失突变噪音水平:
式5和式6中:p表示Panel范围内某个位点的信息;ω 2 代表权重,λ rsl 表示位点p所属插入或缺失突变特征相同分组内所有位点的平均插入或缺失突变频率,μ pl 表示位点p的插入或缺失突变频率,ε pl 表示位点p的插入或缺失突变噪音水平。
上述装置中,所述重复序列元件信息包含DNA重复元件、长插入片段重复元件、短插入片段重复元件、低复杂度区域重复元件、简单重复区域重复元件、长片段重复区域元件、非重复区域、其它重复元件共计8种类型。所述染色体片段重复信息包含位于非大片段重复区域和位于大片段重复区域两种类型。所述人类基因组的重复序列元件数据文件(Repeats文件)可来源于UCSC网站(网址:https://genome.ucsc.edu/cgi-bin/hgTables,参考基因组版本:Hg19,选择RepeatMasker)。所述人类基因组的染色体片段重复数据文件(Segmental Duplication文件)可来源于UCSC网站(网址:https://genome.ucsc.edu/cgi-bin/hgTables,参考基因组版本:Hg 19,选择Segmental Dups)。
上文所述测序数据可为测序产品对应Panel的测序数据。
为了解决上述技术问题,本发明还提供了鉴定ctDNA测序数据中的突变检测结果为真实突变或测序噪音的装置,所述装置可包括如下模块:
C1)ctDNA 测序数据质控和比对模块:用于获得待测样本ctDNA测序数据,经过质控与人类参考基因组进行比对得到比对文件;
C2)ctDNA 测序数据突变分析模块:用于基于所述比对文件获得所述待测样本ctDNA测序数据中突变位点的测序深度、突变reads支持数和突变频率,得到待测样本ctDNA测序数据中的突变检测结果;
C3)统计学检验模块:用于基于所述突变检测结果和噪音水平,使用统计学假设检验的方法确定所述突变检测结果为真实突变或测序噪音。所述噪音水平可利用上文所述的装置获取。
上述装置中,所述突变可为SNV和/或InDel。
上述装置中,所述噪音水平可为SNV噪音水平和/或InDel噪音水平。
上述装置中,所述统计学假设检验可使用式7和式8计算假设检验的统计学显著水平:
式7和式8中:p表示测序产品Panel范围内某个位点的信息,δ p 是位点p的噪音水平(点突变噪音水平,插入或缺失突变噪音水平);d p 是位点p的测序深度;λ p 是位点p期望观测到的噪音突变reads支持数;k p 是位点p在所述待测样本中的突变(点突变,插入或缺失突变)reads支持数;j代表突变reads支持数取值,j∈(0,k p ) ;P(k p , λ p )表示k p 大于等于λ p 的统计学P值大小。
为了解决上述技术问题,本发明还提供了检测ctDNA测序数据中噪音水平的方法,所述方法可包括如下步骤:
B1)ctDNA测序数据处理:获得健康受试者的ctDNA测序数据,经过质控与人类参考基因组进行比对得到比对文件;
B2)突变频率检测:基于所述比对文件,获得所述ctDNA测序数据中每个位点的测序深度、每个位点的点突变reads支持数、每个位点的插入或缺失突变reads支持数,基于所述每个位点的测序深度和所述每个位点的点突变reads支持数计算所述每个位点的点突变频率,基于所述每个位点的测序深度和所述每个位点的插入或缺失突变reads支持数计算所述每个位点的插入或缺失突变频率;
B3)点突变噪音水平模型构建:基于每个位点的点突变特征分组信息和所述每个位点的点突变频率获得所述每个位点的点突变噪音水平;所述每个位点的点突变特征分组信息包括每个位点的三核苷酸序列和每个位点的比对次数;
B4)插入或缺失突变噪音水平模型构建:基于每个位点的插入或缺失突变特征分组信息和所述每个位点的插入或缺失突变频率获得所述每个位点的插入或缺失突变噪音水平;所述每个位点的插入或缺失突变特征分组信息包括每个位点的重复序列元件信息和每个位点的染色体片段重复信息;
B5)噪音水平输出:输出所述ctDNA测序数据的噪音水平;所述噪音水平包括所述每个位点的点突变(SNV)噪音水平和/或所述每个位点的插入或缺失突变(Indel)噪音水平。
上文所述位点可为一个核苷酸,具体可为测序产品Panel范围内包含的位点。所述每个位点的三核苷酸序列包括所述每个位点在参考基因组上的核苷酸及其上下游各一个核苷酸组成的三个核苷酸序列。所述每个位点的三核苷酸序列可由人类参考基因组文件获得。所述每个位点的比对次数可为唯一比对或多次比对。所述每个位点的比对次数可由人类参考基因组的Mappability文件获得。所述每个位点的重复序列元件信息表示所述每个位点是否属于重复序列区域以及属于具体哪个重复区域类型。所述每个位点的重复序列元件信息可由人类基因组的重复序列元件数据文件获得。所述每个位点的染色体片段重复信息表示所述每个位点是否属于大片段重复区域。所述每个位点的染色体片段重复信息可由人类基因组的染色体片段重复数据文件获得。
上述方法中,B3)所述点突变噪音水平模型构建可包括如下步骤:
B3-1)点突变特征整合分组:基于所述每个位点的三核苷酸序列和所述每个位点的比对次数,将具有相同三核苷酸序列且具有相同比对次数的位点分为一组得到点突变特征相同分组;
B3-2)点突变先验背景噪音水平获得:基于所述点突变特征相同分组的每个位点的点突变频率和B3-1)所述点突变特征相同分组使用式1计算获得所述点突变特征相同分组的平均点突变频率:
式1中:λ tm 表示所述点突变特征相同分组内所有位点的平均点突变频率,μ ptm 代表B3-1)所述点突变特征相同分组内每个位点的点突变频率,n代表B3-1)所述点突变特征相同分组内位点的数量;
B3-3)点突变噪音水平获得:基于所述每个位点所属点突变特征相同分组的平均点突变频率和所述每个位点的点突变频率使用式2和式3获得所述每个位点的点突变噪音水平:
式2和式3中:p表示Panel范围内某个位点的信息;λ tm 表示位点p所属点突变特征相同分组的所有位点的平均点突变频率;μ p 表示位点p的点突变频率;ω 1 代表权重;ε p 表示位点p的点突变噪音水平。
上述方法中,所述Mappability文件可来源于UCSC网站(网址:https://genome.ucsc.edu/cgi-bin/hgTables/hgsid=1510554375_HANSVqcS3xxfWO2Ui7wYz8FFIlzA,参考基因组版本:Hg19;100bp reads长度比对版本)。
上述方法中,B4)所述插入或缺失突变噪音水平模型构建可包括如下步骤:
B4-1)插入或缺失突变特征整合分组:基于所述每个位点的重复序列元件信息和所述每个位点的染色体片段重复信息,将具有相同重复序列元件信息且具有相同染色体片段重复信息的位点分为一组得到插入或缺失突变特征相同分组;
B4-2)插入或缺失先验背景噪音水平获得:基于所述每个位点的插入或缺失突变频率和B4-1)所述插入或缺失突变特征相同分组使用式4计算获得所述插入或缺失突变特征相同分组的平均插入或缺失突变频率:
式4中:λ rsl 表示所述插入或缺失突变特征相同分组内所有位点的平均插入或缺失突变频率,μ prsl 代表B4-1)所述插入或缺失突变特征相同分组内每个位点的插入或缺失突变频率,m代表B4-1)所述插入或缺失突变特征相同分组内位点的数量;
B4-3)插入或缺失突变噪音水平获得:基于所述每个位点所属插入或缺失突变特征相同分组的平均插入或缺失突变频率和所述每个位点的插入或缺失突变频率使用式5和式6获得所述每个位点的插入或缺失突变噪音水平:
式5和式6中:p表示Panel范围内某个位点的信息;ω 2 代表权重,λ rsl 表示位点p所属插入或缺失突变特征相同分组内所有位点的平均插入或缺失突变频率,μ pl 表示位点p的插入或缺失突变频率,ε pl 表示位点p的插入或缺失突变噪音水平。
上述方法中,所述重复序列元件信息包含DNA重复元件、长插入片段重复元件、短插入片段重复元件、低复杂度区域重复元件、简单重复区域重复元件、长片段重复区域元件、非重复区域、其它重复元件共计8种类型。所述染色体片段重复信息包含位于非大片段重复区域和位于大片段重复区域两种类型。所述人类基因组的重复序列元件数据文件(Repeats文件)可来源于UCSC网站(网址:https://genome.ucsc.edu/cgi-bin/hgTables,参考基因组版本:Hg19,选择RepeatMasker)。所述人类基因组的染色体片段重复数据文件(Segmental Duplication文件),可来源于UCSC网站(网址:https://genome.ucsc.edu/cgi-bin/hgTables,参考基因组版本:Hg 19,选择Segmental Dups)。
上文所述测序数据可为测序产品对应Panel的测序数据。
为了解决上述技术问题,本发明还提供了鉴定ctDNA测序数据中的突变检测结果为真实突变或测序噪音的方法,所述方法可包括如下步骤:
D1)ctDNA 测序数据质控和比对:获得待测样本ctDNA测序数据,经过质控与人类参考基因组进行比对得到比对文件;
D2)ctDNA 测序数据突变分析:基于所述比对文件获得所述待测样本ctDNA测序数据中突变位点的测序深度、突变reads支持数和突变频率,得到待测样本ctDNA测序数据中的突变检测结果;
D3)统计学检验分析:用于基于所述突变检测结果和噪音水平,使用统计学假设检验的方法确定所述突变检测结果为真实突变或测序噪音;所述噪音水平可利用上文所述的方法获取。
上述方法中,所述突变可为SNV和/或InDel。
上述方法中,所述噪音水平可为SNV噪音水平和/或InDel噪音水平。
上述方法中,所述统计学假设检验可使用式7和式8计算假设检验的统计学显著水平:
式7和式8中:p 表示测序产品Panel范围内某个位点的信息,δ p 是位点p的噪音水平(点突变噪音水平,插入或缺失突变噪音水平);d p 是位点p的测序深度;λ p 是位点p期望观测到的噪音突变reads支持数;k p 是位点p在所述待测样本中的突变(点突变,插入或缺失突变)reads支持数;j代表突变reads支持数取值,j∈(0,k p ) ;P(k p , λ p )表示k p 大于等于λ p 的统计学P值大小。
为了解决上述技术问题,本发明还提供了存储有计算机程序的计算机可读存储介质,所述计算机可读存储介质可为检测ctDNA测序数据中噪音水平的计算机可读存储介质或鉴定ctDNA测序数据中的真实突变或测序噪音的计算机可读存储介质。所述检测ctDNA测序数据中噪音水平的计算机可读存储介质或鉴定ctDNA测序数据中的真实突变或测序噪音的计算机可读存储介质的计算机程序可使计算机执行如上文所述的方法的步骤或运行上文所述装置的模块。
上文所述的装置和/或上文所述方法的下述任一种应用也属于本发明的保护范围:
E1)在制备筛查或辅助筛查癌症患者的产品中的应用;
E2)在开发或制备癌症个体化治疗的产品中的应用。
针对上述问题,本发明的目的是提供一种基于位点的噪音水平来检测突变的方法。
本发明的技术对于单肿瘤样本ctDNA数据检测SNV和InDel,无需肿瘤患者配对正常样本的数据,这使得减少试验及测序成本投入。
为实现上述目的,本发明采取以下技术方案:本发明由于采取以上技术方案,其具有以下优点:
1. 本发明对于ctDNA样本测序数据检测点突变(SNV)和/或插入缺失突变(InDel),无需提供正常对照样本,减少成本投入。
2. 本发明利用健康个体数据集,构建位点的背景噪音水平,可以通过估计噪音水平与真实观测肿瘤样本的突变频率统计学检验差异确定真实突变信号。相较于现有技术所有位点使用统一的阳性判断值,本发明使用位点的阳性判断值,对于区分低频的噪音/突变信号更灵敏。
3. 本发明针对测序数据内不同位点模拟出特异的噪音水平,对于评估位点的最低检测限(Limit of Detection, LoD)具有明显优势。
4. 本发明操作简单,资源消耗小,速度快。因其利用健康个体数据集构建噪音水平模型,分析真实肿瘤样本时,若使用Panel测序产品、试验条件等不变情况下,无需重新构建模型,直接基于噪音水平与真实观测突变频率进行统计学检验即可。
具体实施方式
下面结合具体实施方式对本发明进行进一步的详细描述,给出的实施例仅为了阐明本发明,而不是为了限制本发明的范围。以下提供的实施例可作为本技术领域普通技术人员进行进一步改进的指南,并不以任何方式构成对本发明的限制。
下述实施例中的实验方法,如无特殊说明,均为常规方法,按照本领域内的文献所描述的技术或条件或者按照产品说明书进行。下述实施例中所用的材料、试剂等,如无特殊说明,均可从商业途径得到。
实施例1、基于健康个体的ctDNA测序数据构建位点的噪音水平模型。
本发明基于北京泛生子基因科技有限公司收集的40例健康个体样本(样本年度血液检测报告无异常),使用Onco Sonar产品(CT_DNA_170实体瘤基因突变检测试剂盒(可逆末端终止测序法),北京泛生子基因科技有限公司,产品货号RSN113)依据标准试验流程,经DNA提取、文库构建、探针捕获、上机测序等步骤得到40个无已知癌症的健康人个体样本的ctDNA原始测序数据。本实施例包括对原始数据预处理,获得突变频率、确定突变分组特征、构建噪音水平模型等步骤。
1.原始数据预处理。
原始测序数据为FASTQ格式,对40例健康个体样本的ctDNA原始测序数据进行质控(包括数据量、去重前/后深度、Q30和捕获效率等)、比对等标准数据分析流程。其中质控、比对软件分别为Trimmomatic(v0.33,http://www.usadellab.org/cms/index.php/page=trimmomatic)和BWA(v0.7.10-r789,https://bio-bwa.sourceforge.net/)。
质控软件Trimmomatic需要设置参数,提供测序使用接头序列,滑动窗口剪切窗口大小为4bp,滑动窗口内碱基平均质量值阈值不低于15;原始数据中测序reads的起始端切除质量值低于3的碱基;测序reads的末端切除质量值低于3的碱基,其余参数无明确说明即使用软件默认参数。ctDNA样本质控合格标准为:Q30≥80%,去除重复序列(简称“去重”)前测序深度≥15000X,去重率≥2.5。BWA软件将质控合格的测序数据比对到人类参考基因组,参考基因组版本为Hg19,获得健康个体样本有效比对文件BAM文件(包含有效比对数据),用于后续分析。
2. 测序产品Panel范围内检测位点的突变频率获得。
基于BAM文件统计Panel范围内(Onco Sonar产品,Panel大小约0.3M)每个位点的测序深度及检出4种可能核苷酸A/C/T/G 点突变reads支持数,计算得到每个位点4种核苷酸的点突变频率(未检测到点突变的频率记为0)。具体为将BAM文件使用Samtools(v1.7,http://www.htslib.org/doc/samtools.html)软件中的mpileup命令生成pileup格式文件,经程序(北京泛生子基因科技有限公司,ctDNA-SNV-InDel-pipeline),对Panel范围内每个位点的测序深度(即每个位点的测序总reads数)及突变reads支持数进行统计,利用每个位点发生点突变reads支持数除以该位点的测序深度即可计算得到该位点发生点突变的突变频率。同理,基于BAM文件统计Panel范围内每个位点的测序深度(位点测序总reads数)及检出插入或缺失突变reads支持数,利用每个位点插入或缺失突变reads支持数除以该位点的测序深度计算得到该位点插入或缺失突变的突变频率。
3. 筛选特征分组并构建位点的点突变的噪音水平模型。
3.1筛选点突变特征分组,确定每个组别点突变先验背景噪音水平。
经美国加州大学圣克鲁兹分校开发和维护的数据库(以下简称UCSC数据库)下载人类参考基因组的Mappability文件(网址:https://genome.ucsc.edu/cgi-bin/hgTables/hgsid=1510554375_HANSVqcS3xxfWO2Ui7wYz8FFIlzA)(基因组版本:Hg19;100bpreads长度比对版本),该文件给出某个基因组序列位置在全基因组范围内是唯一比对或多次比对。基于参考基因组文件(版本:Hg19),提取位点的三核苷酸序列(tri-nucleotidecontext)(即某个位点在参考基因组上的核苷酸序列及其上游一个核苷酸和下游一个核苷酸组成的三个核苷酸序列)。三核苷酸序列及比对次数(即唯一比对还是多次比对)作为后续构建模型的特征。
将每个位点的三核苷酸序列和比对次数组合进行特征分组,具体为将具有相同三核苷酸序列且具有相同比对次数的位点分为一组,组内包括的位点即为点突变特征相同的位点。根据不同三核苷酸序列及比对次数最多可分为128个组(每个位点上核苷酸有A/T/C/G 4种可能形式,三个核苷酸排列最多有43=64种;比对次数分为2种情况:唯一比对和多次比对;三核苷酸序列特征与比对次数特征可最多分为64×2=128个分组,每一个组被称为点突变特征相同分组。基于上文内容将Panel测序范围内的位点进行分组(每个位点对应的三核苷酸序列及比对次数是唯一的,所以每个位点只会属于唯一分组),根据步骤2得到的位点的点突变频率,将点突变特征相同分组内所有位点(特征相同位点)的点突变频率取平均值来计算该分组的平均点突变频率(式1),该平均点突变频率作为该分组的先验背景噪音水平:
式1中:λ tm 表示点突变特征相同分组内所有位点的平均点突变频率,μ ptm 代表点突变特征相同分组内每个位点的点突变频率,n代表点突变特征相同分组内位点的数量。
3.2构建位点的点突变噪音水平模型。
结合每个位点所在点突变特征相同分组及其对应的点突变先验背景噪音水平(步骤3.1获得的平均点突变频率)及步骤2检测到的每个位点的点突变频率构建模型,计算当前位点发生点突变的噪音水平:
式2和式3中:p表示测序产品Panel范围内某个位点的信息,如:chr7: 55249071;λ tm 表示位点p所属点突变特征相同分组的所有位点的平均点突变频率;μ p 表示位点p的点突变频率;ω 1 代表权重;ε p 表示位点p的点突变噪音水平。
4. 筛选插入或缺失突变特征分组并构建位点的插入或缺失突变的噪音水平模型。
4.1 筛选插入或缺失突变的分组特征,确定每个组别插入或缺失突变先验背景噪音水平。
经UCSC数据库下载人类基因组的重复序列元件数据文件(Repeats文件,网址:https://genome.ucsc.edu/cgi-bin/hgTables,基因组版本Hg19,选择RepeatMasker)和人类基因组的染色体片段重复数据文件(Segmental Duplication文件,网址:https://genome.ucsc.edu/cgi-bin/hgTables,基因组版本Hg 19,选择Segmental Dups)。Repeats文件提供某个基因组位点是否属于重复序列区域及属于具体重复区域的类型,分为DNA、LINE、SINE、Low_complexity、Simple_repeats、LTR、NonRepeats、Other等共计8种不同类型(DNA类型代表DNA重复元件;LINE类型表示长插入片段重复元件;SINE类型表示短插入片段重复元件;Low complexity类型表示低复杂度区域重复元件;比如G-rich,AT-rich区域;Simple repeats类型代表简单重复区域重复元件,如(TTAGGG)n,(GATATA)n;LTR类型表示长片段重复区域元件;NonRepeats类型表示位于非重复区域;Other类型表示其它重复元件,具体类型详细信息可参考UCSC数据库)。Segmental Duplication文件提供某个基因组位点是否属于大片段重复区域(染色体重复区域可分为位于非大片段重复区域和位于大片段重复区域)。位点的重复序列元件信息(Repeats信息)及其染色体片段重复信息(Segmental Duplication信息)作为后续构建插入或缺失突变噪音水平模型的特征。
将具有相同Repeats信息(即属于上述8种类型的同一种,例如都是DNA重复元件)和Segmental Duplication信息的位点分为一组,称为插入或缺失突变特征相同分组,组内包括的位点即为插入或缺失突变特征相同的位点。Repeats信息可分为8组;SegmentalDuplication信息可分为2组,Repeats信息与Segmental Duplication信息组合可最多分为8×2=16个分组,每一个组被称为插入或缺失突变特征相同分组。将Panel范围内的位点进行分组(每个位点的Repeats信息及Segmental Duplication信息是唯一的,所以每个位点只会属于唯一分组)。基于步骤2得到的位点插入或缺失突变的突变频率,将插入或缺失突变特征相同分组内所有位点的插入或缺失突变频率取平均值来计算其平均突变频率(式4),该信息作为该分组的插入或缺失突变的先验背景噪音水平:
式4中:λ rsl 表示插入或缺失突变特征相同分组内所有位点的平均插入或缺失突变频率,μ prsl 代表插入或缺失突变特征相同分组内每个位点的插入或缺失突变频率,m代表插入或缺失突变特征相同分组内位点的数量。
4.2 构建位点的插入或缺失突变的噪音水平模型。
结合每个位点分组及对应的插入或缺失突变的先验背景噪音水平(步骤4.1获得的特征相同分组的平均突变频率)及步骤2统计得到的每个位点的插入或缺失突变频率使用式5和式6构建模型,计算当前位点发生插入或缺失突变的噪音水平。
式5和式6中:p表示测序产品Panel范围内某个位点的信息,ω 2 代表权重,λ rsl 表示位点p所属插入或缺失突变特征相同分组内所有位点的平均插入或缺失突变频率,μ pl 表示位点p的插入或缺失突变频率;ε pl 表示位点p的插入或缺失突变噪音水平。
位点的点突变噪音水平模型和位点的插入或缺失突变的噪音水平模型共同构成ctDNA NGS测序数据噪音水平模型。
实施例2、鉴定ctDNA 测序数据中的突变检测结果为真实突变或测序噪音的流程。
对于待测样本,基于与实施例1相同的标准试验流程,经DNA提取、文库构建、探针捕获和上机测序等流程得到ctDNA测序数据;再经过数据预处理、比对等步骤,对于质量控制合格样本,统计位点的测序深度及突变reads支持数和突变频率,得到待测样本的ctDNA测序数据的突变检测结果。根据实施例1中构建的ctDNA测序数据噪音水平模型,对待测样本内发生突变的位点,使用统计学检验方法来鉴定样本ctDNA 测序数据中的突变(点突变、插入或缺失突变)检测结果是真实突变还是测序噪音。
将样本ctDNA NGS测序数据使用实施例1中的方法进行质控及比对到参考基因组(参考基因组版本:Hg19)的步骤,得到BAM文件;基于BAM文件统计突变位点测序深度、突变支持数(支持突变的reads数),计算得到位点的突变频率;基于位点的信息(即发生突变位点的染色体号及基因组位置)及突变的类型(突变是点突变还是插入或缺失突变类型)匹配实施例1中构建模型获得的位点的噪音水平(点突变类型匹配点突变噪音水平,插入或缺失突变匹配插入或缺失突变的噪音水平)。根据统计学假设检验思想,H0:该位点突变为测序噪音;H1:该位点突变为真实突变;根据突变位点的测序深度、突变支持数以及该位点的突变噪音水平计算位点统计学显著水平(P-value)(式7和式8);若P-value<alpha(阳性阈值),则应拒绝H0,接受H1,即检测到的突变为真实突变;反之,接受H0,即检测到的突变为测序噪音。
式7和式8中:p 表示测序产品Panel范围内某个位点的信息,δ p 是位点p的噪音水平(点突变噪音水平,插入或缺失突变噪音水平);d p 是位点p的测序深度,λ p 是位点p的噪音突变reads支持数;k p 是位点p在所述待测样本中突变(点突变,插入或缺失突变)reads支持数;j代表突变reads支持数取值,j∈(0,k p );P(k p , λ p )表示k p 大于等于λ p 的统计学P值大小。
本实施例中突变阳性阈值为0.0001。
实施例3、本发明基于肿瘤样本ctDNA测序数据对突变检测结果的验证。
基于肿瘤样本ctDNA 测序数据将本发明实施例1和2中建立的突变检测结果判定方法与现有技术进行比较,评估指标是已知阳性变异(即点突变、插入或缺失突变)检出数量。选用北京泛生子基因科技有限公司收集10例既往临床检测肿瘤样本(包括5例肺癌和5例肠癌样本),按照Onco Sonar产品试验流程,经DNA提取、等比例混样,文库构建、探针捕获和上机测序等流程得到NGS测序数据。
采用现有技术(即使用samtools软件mpileup命令统计ctDNA测序数据内发生变异位点总深度、突变支持数及突变频率。真实突变判断标准:突变频率大于等于0.1%,突变支持数大于等于4,参考网址http://www.htslib.org/doc/1.7/samtools.html)分别对10例样本测序数据单独进行分析(即不包含DNA等比例混样过程)。检测结果(见图2)显示:阳性突变位点总数648个,其中低频突变(突变频率<0.5%)位点 486个,低频突变中,基因突变频率AF(Allele Frequency)≤0.1%位点192个,0.1<AF≤0.2%位点96个,0.2<AF≤0.3%位点96个,0.3<AF≤0.5%位点102个。
分别采用现有技术和实施例2建立的方法对混样后的测序数据进行分析,并对突变结果进行比较。通过比较已知阳性突变位点检出突变比例发现,本发明的突变检测方法在低频突变(低于0.5%)检出上具有优势,性能可提升约26%。
计算方法:如图2所示,针对某个频率分组,利用本发明方法检出突变比例减去现有方法检出突变比例的差值除以现有方法检出突变比例即当前分组性能提升比例;具体计算步骤,AF≤0.1%分组性能提升比例:(0.197-0.0989)/0.0989=99%,0.1<AF≤0.2% 分组性能提升比例:(0.5-0.3437)/0.3437=45%,0.2<AF≤0.3%分组性能提升比例:(0.8333-0.7395)/0.7395=13%,0.3<AF≤0.5% 分组性能提升比例:(0.8725-0.8828)/0.8823=-1%,0.5<AF≤2.5%分组性能提升比例:(1-1)/1=0%,AF>2.5% 分组性能提升比例:(1-1)/1=0%。基于上述所述每个频率分组性能提升的比例,再计算6个频率分组性能提升的平均值即为总体性能提升比例:(99%+45%+13%-1%+0%+0%)/6=26%。
综上,本发明基于ctDNA测序数据中位点的测序深度、突变reads支持数及该位点噪音水平计算位点统计学显著水平P-value,不再依赖于统一划定的阳性突变的判断标准,尤其适用于低频突变的检出,不会因为未达到阈值从而导致漏检;本发明无需提供正常对照样本,减少成本投入及缩短分析周期。因此,本发明可应用于肿瘤样本检测低频的突变,从而为肿瘤患者提供个性化用药指导。
以上对本发明进行了详述。对于本领域技术人员来说,在不脱离本发明的宗旨和范围,以及无需进行不必要的实验情况下,可在等同参数、浓度和条件下,在较宽范围内实施本发明。虽然本发明给出了特殊的实施例,应该理解为,可以对本发明作进一步的改进。总之,按本发明的原理,本申请欲包括任何变更、用途或对本发明的改进,包括脱离了本申请中已公开范围,而用本领域已知的常规技术进行的改变。