CN111370057B - 确定样本染色体结构变异信号强度以及插入片段长度分布特征的方法及应用 - Google Patents

Abstract

本发明提出了确定样本染色体结构变异信号强度以及插入片段长度分布特征的方法及应用。具体地，涉及一种确定样本来源的方法。该方法包括：(1)将样本染色体的测序读段数据与参考基因组进行比对，确定样本染色体的测序读段数据的低质量比对率和高质量比对率；(2)基于高质量比对率所对应的读段，确定样本染色体的结构变异的比例、线粒体DNA的含量、预定长度插入片段的比例；(3)基于预定肿瘤预测模型，以及步骤(2)所获得的染色体的结构变异的比例、线粒体DNA的含量、预定长度插入片段的比例，确定样本来源的概率；(4)基于样本来源的概率，确定样本来源。

Description

确定样本染色体结构变异信号强度以及插入片段长度分布特 征的方法及应用

技术领域

本发明涉及生物信息领域，具体地，本发明涉及确定样本染色体结构变异信号强度以及统计样本插入片段长度分布的方法及应用，更具体地，本发明涉及确定样本染色体结构变异信号强度和插入片段的分布，线粒体的拷贝数的方法以及确定样本来源的方法。

背景技术

据报道，大约有10％的癌症表现出大量的染色体结构变异，染色体破碎(Chromothripsis)。染色体发生的结构变异主要有4种：

1.缺失染色体中某一片段的缺失例如，猫叫综合征是人的第5号染色体部分缺失引起的遗传病，因为患病儿童哭声轻，音调高，很像猫叫而得名。猫叫综合征患者的两眼距离较远，耳位低下，生长发育迟缓，而且存在严重的智力障碍；果蝇的缺刻翅的形成也是由于一段染色体缺失造成的。

2.重复染色体增加了某一片段果蝇的棒眼现象就是X染色体上的部分重复引起的。

3.倒位染色体某一片段的位置颠倒了180度，造成染色体内的重新排列如女性习惯性流产(第9号染色体长臂倒置)。

4.易位染色体的某一片段移接到另一条非同源染色体上或同一条染色体上的不同区域如人慢性粒白血病(第22号染色体的一部分易位到第14号染色体上造成)，夜来香也经常发生这样的变异。部分肿瘤样本是由于染色体结果变异引起的，同时也有近10％的样本存在大量的结构变量(即存在Chromothripsis)。因此通过测试样本中的染色体变异信号强度，可以区分肿瘤的来源。

同时，最近多项研究表明，基于低深度全基因组测序的cfDNA(circulating cell-free DNA，cfDNA)的片段化特征可以一定程度上反映它的细胞来源和它入血的分子机理。通常健康个体的cfDNA片段大小分布通常主要集中在167bp，这主要反映了白细胞核小体的位置分布。而癌症患者的cfDNA片段分布特征谱则与之不同，主要体现在其cfDNA整体片段大小分布比正常人的cfDNA短。通过筛选小片段cfDNA(如小于150bp)可以有效富集晚期肿瘤患者cfDNA中肿瘤来源的ctDNA(circulating tumor DNA,ctDNA)，而不用获知这些DNA片段的基因突变信息。因此，基于低深度全基因组测序的cfDNA片段化特征分析方法可以通过片段筛选来富集ctDNA，提高突变检出限，可以作为目标区域高深度测序方法的补充，被应用于科学研究和癌症筛查、用药指导、复发监控等临床场景。

发明内容

本申请是基于发明人对以下事实和问题的发现和认识作出的：

发明人发现，染色体结构变异即使在肿瘤样本中发生的概率是很低的，除此之外，其他的肿瘤样本和正常人样本并不存在差异。因此，单一地检测染色体结构变异的信号强度来预测样本的来源，其敏感性上往往不太理想。而本申请发明人创造性地将染色体结构变异信号强度与插入片段长度分布特征进行拟合，通过优化拟合流程和参数，大大提高了确定样品来源的准确性，为后续的基于样本数据的科学研究，如药物筛选或临床应用，如癌症筛查、用药指导、复发监控奠定了基础。

在本发明的第一方面，本发明提出了一种确定样本来源的方法。根据本发明的实施例，所述方法包括：(1)将样本染色体的测序读段数据与参考基因组进行比对，确定样本染色体的测序读段数据的低质量比对率和高质量比对率；(2)基于高质量比对率所对应的读段，确定样本染色体的结构变异的信号比例、线粒体DNA的含量、预定长度插入片段的比例；(3)基于预定肿瘤预测模型，以及步骤(2)所获得的染色体的结构变异的信号比例、线粒体DNA的含量、预定长度插入片段的比例，确定样本来源的概率；(4)基于样本来源的概率，确定样本来源。其中，样本染色体的结构变异的信号比例反映了样本染色体结构变异的信号强度。根据本发明实施例的方法可准确判定样本的来源，为后续的基于样本数据的科学研究，如药物筛选或临床应用，如癌症筛查、用药指导、复发监控打下了坚实的基础。

根据本发明的实施例，上述方法还可以进一步包括如下附加技术特征至少之一：

根据本发明的实施例，所述预定长度插入片段的比例包括所述预定长度插入片段的读段数在样本染色体的测序读段的总数目的比例和样本染色体预定区间的所述不同预定长度插入片段的读段数的比值。

根据本发明的实施例，所述测序读段数据的低质量比对率是通过如下方式确定的：确定任意一端测序读段的比对质量值小于30的测序读段的数目，其中，每个读段的比对质量值可由比对软件(BWA)自动输出，表示这个读段比对位置错误的概率，具体计算方法可参见Li H,Ruan J,Durbin R.(2008)Genome Research 18:1851-8；基于任意一端测序读段的比对质量值小于30的测序读段的数目与测序读段的总数目的比值，确定测序读段数据的低质量比对率。进而作为样本的比对质控。

根据本发明的实施例，所述测序读段数据的高质量比对率是通过如下方式确定的：确定测序读段两端的比对质量值均大于30的测序读段的数目；基于测序读段两端的比对质量值均大于30的测序读段的数目与测序读段的总数目的比值，确定测序读段数据的高质量比对率。

根据本发明的实施例，所述染色体结构变异的信号比例是通过如下公式确定的：

y＝x/A,

其中，y表示染色体结构变异的信号比例，

x表示异常比对测序读段的数目，

A表示样本染色体的测序读段的总数目。

需要说明的是，本申请所述的正向读段和反向读段是基于预定序列的正向测序读段和反向测序读段，当认定一个方向的测序读段为正向读段时，则另一方向的测序读段即为反向读段，因此，本申请所述的正向读段或反向读段应做广义理解，它们是指一对测序读段中所包含的两个方向的读段，而不将某一方向的读段特指为正向或反向读段。例如，当描述“所述正向读段与第一染色体的第一预定参考序列完全匹配，所述反向读段与第二预定参考序列完全匹配”时，也可以理解为“所述反向读段与第一染色体的第一预定参考序列完全匹配，所述正向读段与第二预定参考序列完全匹配”，这两种表述方式在本申请中所描述的技术方案一致。

根据本发明的具体实施例，所述异常比对测序读段是指：基于一对测序读段，所述一对测序读段包括预定序列的正向读段和反向读段，所述正向读段的一部分序列与第一染色体的第一预定参考序列完全匹配，另一部分与第二预定参考序列完全匹配，所述第二预定参考序列与所述第一预定参考序列位于同一染色体，并且二者的距离大于10000bp或所述第二预定参考序列位于不同于第一染色体的第二染色体上。所述反向读段与第三预定参考序列完全匹配，并且第三预定序列与第一预定参考序列或者第二预定参考序列在同一染色体上，且位置距离小于1000bp，所述染色体结构变异的信号比例是通过如下公式确定的：

y¹＝x¹/A，

其中，y₁表示在上述异常比对测序读段情况下的染色体结构变异的信号比例，

_x1表示上述异常比对测序读段情况下的异常比对测序读段的数目，

A表示样本染色体的测序读段的总数目。

为了便于理解，申请人将上述异常比对测序读段情况下的染色体结构变异表示为图1。

根据本发明的具体实施例，所述异常比对测序读段是指：基于一对测序读段，所述一对测序读段包括预定序列的正向读段和反向读段，所述正向读段与第一染色体的第一预定参考序列完全匹配，所述反向读段与第二预定参考序列完全匹配，所述第二预定参考序列与所述第一预定参考序列位于同一染色体并且二者的距离大于10000bp或所述第二预定参考序列位于不同于第一染色体的第二染色体上，所述染色体结构变异的比例是通过如下公式确定的：

y²＝x²/A，

其中，y₂表示上述异常比对测序读段情况下的染色体结构变异的信号比例，

_X2表示上述异常比对测序读段情况下的所述异常比对测序读段的数目，

A表示样本染色体的测序读段的总数目。

为了便于理解，申请人将上述异常比对测序读段情况下的染色体结构变异表示为图2。

根据本发明的具体实施例，所述异常比对测序读段是指：基于一对测序读段，所述一对测序读段包括预定序列的正向读段和反向读段，所述正向读段和反向读段的一部分同时比对到第一染色体的第一预定参考序列，所述正向读段和反向读段的剩余部分同时比对到第二预定参考序列完全，且比对到第二预定参考序列的起始位置一样，所述第二预定参考序列与所述第一预定参考序列位于同一染色体并且二者的距离大于10000bp或所述第二预定参考序列位于不同于第一染色体的第二染色体上，所述染色体结构变异的信号比例是通过如下公式确定的：

y³＝x³/A，

其中，y₃表示上述异常比对测序读段情况下的染色体结构变异的信号比例，

_X3表示上述异常比对测序读段情况下的异常比对测序读段的数目，

A表示样本染色体的测序读段的总数目。

为了便于理解，申请人将上述异常比对测序读段情况下的染色体结构变异表示为图3。

根据本发明的具体实施例，所述异常比对测序读段是指：基于一对测序读段，所述一对测序读段包括预定序列的正向读段和反向读段，所述正向读段或反向读段分别与第一预定参考序列的正义链完全匹配和第二预定参考序列的正义链完全匹配；或者分别与第一预定参考序列的反义链完全匹配和第二预定参考序列的反义链完全匹配；或者正向读段与第一染色体的第一预定参考序列的正义链完全匹配，反向读段与第一染色体的第二预定参考序列的反义链完全匹配，第二预定参考序列在第一染色体上的位置小于第一预定参考序列，

y⁴＝x₄/A，

其中，y₄表示上述异常比对测序读段情况下的染色体结构变异的信号比例，

_X4表示上述异常比对测序读段情况下的所述异常比对测序读段的数目，

A表示样本染色体的测序读段的总数目。

为了便于理解，申请人将上述异常比对测序读段情况下的染色体结构变异表示为图4。

根据本发明的实施例，所述线粒体DNA的含量是通过如下方式确定的：确定比对到参考线粒体基因序列的测序读段的数目；基于比对到参考线粒体基因序列的测序读段的数目与测序读段的总数目的比值，确定线粒体DNA的含量。

根据本发明的实施例，确定预定长度插入片段的读段数在样本染色体的测序读段的总数目的比例是通过如下方式进行的：将所述预定长度插入片段的长度设置为70～150bp，150～180bp，180～220bp或250～330bp。发明人发现，同一样本中，长度为70～150bp插入片段的比例(记为P70～150)与长度为180～220bp插入片段的比例(记为P180～220)之间存在着很强的负相关性，同时将癌症样本与正常样本的插入片段的比例进行t检验发现，P180～220的P-value更小，因此，统计P180～220可使结果的准确性进一步提高；同时，发明人发现，长度在250～330bp插入片段在癌症样本中也出现较多，因此，同时统计P250～330bp也可使结果的准确性进一步提高。

根据本发明的实施例，进一步包括确定峰谷间距，所述峰谷间距是指插入片段在小于150bp范围内，波峰与波谷对应长度的插入片段及其上下2bp[±2bp]的测序读段数目占样本染色体的测序读段的总数目的比例的差值。

根据本发明的实施例，确定样本染色体预定区间的所述不同预定长度插入片段的读段数的比值是通过如下方式进行的：将参考基因序列划分为多个相同长度的窗口区间，任选地，所述窗口区间的大小为100kb；确定每个窗口区间内预定长度插入片段的读段数目，任选地，预定长度插入片段的长度为100～150bp或151～220bp；确定每个窗口区间内不同预定长度插入片段的读段数目的比值。发明人发现，不同来源的样本(肿瘤VS正常)的预定长度插入片段在样本测序读段数目的总比例之间的差异小，但是稳定性好；而一定窗口区间内的不同长度预定插入片段之间的比值却差异明显，但是却受外界因素影响大(比如GC含量，序列复杂度等)，因此，综合分析预定长度插入片段在样本测序读段数目的总比例和预定区间内不同预定长度插入片段的读段数目的比值，可显著提高检测和预测结果的真实性。

根据本发明的实施例，在每个窗口区间内，进一步包括对预定长度插片片段的读段数目进行校正处理。进而进一步提高对样品来源判断的准确性。

根据本发明的实施例，所述校正处理是通过将在每个窗口区间内预定长度的插入片段的读段数目的中位值加上片段数目残差获得的。进而可消除因为每个窗口区间的GC和比对率的不同，而对预定长度插入片段的读段数目造成的影响。

根据本发明的实施例，所述片段数目的残差是通过如下方式获得的：(1)确定每个窗口区间内的GC含量和比对率，其中，每个窗口区间内的GC含量是通过如下方式确定的：统计每个窗口(bin)内的A,T,C,G碱基的数量，计算GC碱基的数量占总碱基数量的比值，即为该窗口的GC含量，每个窗口区间的比对率(mappability)是通过如下方式获得的：根据从UCSC下载的ENCODE’s mappability bigwig文件，将文件中的每个区域的mappability与划分的窗口区间比较，计算出所有区间的比对率的平均值，作为该区间的比对率；(2)将步骤(1)所获得的每个窗口区间内的GC含量和比对率进行组合和分组处理，获得每个GC含量和比对率组合对应所有窗口区间的测序读段数目的中位值；(3)基于局部加权非参数回归方法，构建GC含量和比对率组合对应窗口区间的测序读段数目中位值相对于GC含量和比对率的拟合曲线；(4)基于所述拟合曲线以及每个窗口区间内的GC含量和比对率，确定每个窗口区间内的理论插入片段数目；(5)将每个窗口区间内的预定长度的插入片段的读段数目减去步骤(4)所获得的理论插入片段数目，获得每个窗口区间内的预定长度的插入片段数目的残差。通过上述方法获得片段数目的残差可有效消除每个窗口区间内GC含量和比对率的差异。

根据本发明的实施例，进一步包括将预定区间内的预定长度插入片段的数目进行加和处理，所述加和处理处理包括分别将插入片段的长度为100～150bp的插入片段的读段数目进行加和和将插入片段的长度为151～220bp的插入片段的读段数目进行加和，

根据本发明的实施例，所述加和处理后的区间的长度可为1M，5M，10M等。根据本发明的具体实施例，所述加和处理后的区间长度为5M。进而有效提高了信号的稳定性和差异的显著性。

根据本发明的具体实施例，进一步包括将插入片段的长度为100～150bp的插入片段的读段数目加和除以插入片段的长度为151～22bp的插入片段的读段数目加和，以便获得插入片段读段数目加和的比值。

根据本发明的实施例，进一步包括对所述加和的比值进行主成分分析，根据累积方差(或者贡献率)大于某一阈值，比如80％，85％，90％，95％来决定选择的主成分个数。

根据本发明的实施例，所述预定肿瘤预测模型包括选自SVM、Lasso、GBM的至少之一。

根据本发明的实施例，所述预定肿瘤预测模型为GBM，基于ROC曲线，以及预定的敏感性和特异性，确定相应阈值。根据本发明的具体实施例，所述预定特异性为98％，对应的所述阈值为0.45。

根据本发明的实施例，在步骤(4)中，所述样本来源肿瘤的概率大于0.45，是样本来源于肿瘤患者的指示。

总之，为了便于理解，申请人将上述根据本发明实施例的确定样本来源的方法的流程总结为图5。

在本发明的第二方面，本发明提出了一种确定样本来源的系统。根据本发明的实施例，所述系统包括：比对模块，所述比对模块用于将样本染色体的测序读段数据与参考基因进行比对，确定样本染色体的测序读段数据的低质量比对率和高质量比对率；获取模块，所述获取模块用于基于高质量比对率所对应的读段，确定样本染色体的结构变异的比例、线粒体DNA的含量、预定长度插入片段的比例；确定样本来源概率模块，所述确定样本来源概率模块用于基于预定肿瘤预测模型，以及获取模块所获得的染色体的结构变异的比例、线粒体DNA的含量、预定长度插入片段的比例，确定样本来源的概率；判定模块，所述判定模块用于基于确定样本来源概率模块所确定的样本来源的概率，确定样本来源。根据本发明实施例的系统适于执行前面所述的确定样本来源的方法，进而可对样本的来源作出准确判断，为后续的基于样本数据的科学研究，如药物筛选或临床应用，如癌症筛查、用药指导、复发监控打下了坚实的基础。

为了便于理解，申请人将上述确定样本来源的系统的结构表示为图6。根据本发明的实施例，所述系统包括：比对模块100，所述比对模块100用于将样本染色体的测序读段数据与参考基因进行比对，确定样本染色体的测序读段数据的低质量比对率和高质量比对率；获取模块200，所述获取模块200用于基于比对模块100所输出的高质量比对率所对应的读段，确定样本染色体的结构变异的比例、线粒体DNA的含量、预定长度插入片段的比例；确定样本来源概率模块300，所述确定样本来源概率模块300用于基于预定肿瘤预测模型，以及获取模块200所输出的染色体的结构变异的比例、线粒体DNA的含量、预定长度插入片段的比例，确定样本来源的概率；判定模块400，所述判定模块400用于基于确定样本来源概率模块300所输出的样本来源的概率，确定样本来源。

根据本发明实施例的系统具有与前面所述的确定样本来源的方法的相同的附加技术特征，其优势和效果类似，在此不再赘述。

在本发明的第三方面，本发明提出了一种确定样本来源的电子设备。根据本发明的实施例，所述电子设备包括存储器、处理器；其中，所述处理器通过读取所述存储器中存储的可执行程序代码来运行与所述可执行程序代码对应的程序，以用于实现前面所述的确定样本来源的方法。根据本发明实施例的电子设备可用来检测待测样本是否来源于肿瘤患者，检测准确性高且成本低。

在本发明的第四方面，本发明提出了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序。根据本发明的实施例，该程序被处理器执行时实现如前面所述的确定样本来源的方法。根据本发明实施例的系统适于执行前面所述的确定样本来源的方法，进而可对样本的来源作出准确判断，为后续的基于样本数据的科学研究，如药物筛选或临床应用，如癌症筛查、用药指导、复发监控打下了坚实的基础。

附图说明

图1是根据本发明实施例的一种异常比对测序读段情况下的染色体结构变异示意图；

图2是根据本发明实施例的一种异常比对测序读段情况下的染色体结构变异示意图；

图3是根据本发明实施例的一种异常比对测序读段情况下的染色体结构变异示意图；

图4是根据本发明实施例的3种异常比对测序读段情况下的染色体结构变异示意图；

图5是根据本发明实施例的确定样本来源的方法流程图；

图6是根据本发明实施例的确定样本来源的系统的结构示意图；

图7是根据本发明实施例的染色体结构变异信号强度在癌症样本和正常样本中的分布差异，使用t检验进行比较，两组存在显著差别(P-value值<0.05)；

图8是根据本发明实施例的某一个样本的插入片段长度，图中箭头指向小于150bp的局部峰值(次高峰)；

图9是根据本发明实施例的正常样本和癌症样本在插入片段180～220之间的比例的分布差异(P-value值是通过t检验获得)；

图10是根据本实施例的主成分分析(PCA)结果；以及

图11是根据本实施的10折交叉验证，在测试集上的预测ROC曲线，其中，10折交叉验证过程中的平均特异性和敏感性见图中黑色粗线。

具体实施例方式

实施例1样本血浆分离，文库制备，上机测序；

1.血浆分离

a)准备好实验所需的仪器、试剂、耗材，高速冷冻离心机应提前预冷至4℃。

b)如果外周血样本是用EDTA抗凝管采集的，抽血之后立马放进4℃冰箱，并在2小时内进行血浆分离。如果外周血样本是用streck管等游离核酸保存管采集的，则可在常温放置，并在采血管说明书规定的时间内进行血浆分离分离。

c)记录样本信息，将采血管配平，将高速冷冻离心机换成水平转子，并设定参数：温度4℃，离心力1600g，时间10min。将采血管配平之后放置在离心机中，进行离心。

d)离心完成之后，将采血管放置在生物安全柜的离心管架上。将离心后采血管中的上清收集至新的15mL离心管中，在管壁标记样本编号以及操作时间。注意在收集上清时需要仔细操作，避免吸入白细胞。剩下的血细胞用于提取gDNA，分装至新的15mL离心管中，在管壁标记样本编号以及操作时间。

e)将高速冷冻离心机换成角转子，并设定参数：温度4℃，离心力16000g，时间10min。将装有上清的15mL离心管配平之后放置在离心机中，进行离心。

f)离心完成之后，将装有上清的15mL离心管放置在生物安全柜的离心管架上。将离心后离心管中的上清收集至新的15mL离心管中。注意在收集上清时需要仔细操作，避免吸入沉淀。这一步的目的是去除血浆当中的细胞碎片等杂质。

g)将血浆以及血细胞放置于-80℃冰箱保存，备用。

h)实验完成后，将所有物品归位，并清洁实验台面，将生物安全柜紫外灯打开，照射30min后关闭。记录详细的实验记录。

2.cfDNA提取

i)准备好实验所需的仪器、试剂、耗材。打开水浴锅，并调节温度至60℃。打开金属浴，并调节温度至56℃。确认试剂盒有效期，buffer ACB是否加有合适量的异丙醇，bufferACW1以及buffer ACW1是否加有合适量的无水乙醇。

j)记录样本编号等信息。

k)若是分离的新鲜血浆，则直接进行cfDNA提取。若血浆冻存在–80℃条件下，需将血浆样本解冻后，在16,000x g[固定角转头]的离心力以及4℃的温度条件下离心5min以去除冷冻沉淀。

l)按照表1配置所需量的ACL混合液。

表1：处理4ml样本所需的Buffel ACL以及carrier RNA(溶解于Buffer AVE)体积用量

m)转移400μl Proteinase K至装有4ml血浆的50ml离心管中。间断涡旋30s以充分混匀。

n)加入3.2ml的Buffer ACL(含有1.0μg carrier RNA)。剧烈涡旋混匀15秒。确保离心管经剧烈涡旋，以保证样本和Buffer ACL的重复混匀，从而实现高效的裂解。

o)注意：此步完成后请不要中断实验并立即进行下步的裂解孵育步骤。

p)将离心管接着60℃水浴30分钟。

q)向上述反应液中加入7.2ml的Buffer ACB。盖上管盖，间断涡旋15s以充分混匀。

r)将含有Buffer ACB的裂解液至于冰上孵育或冷藏孵育5min。

s)组装抽滤装置：把VacValve插在24孔底上，再把VacConnectors插入VacValve中，再将QIAamp Mini硅胶膜柱连接到VacConnectors上，最后把20ml扩容管插入到硅胶膜柱上。确保扩容管插入紧实以防止样本泄露。注意：将2ml收集管留下至后续空转时才使用。并在硅胶膜柱上做好样本编号的标记。VacValve可调节流速，VacConnectors可以防止污染，QIAamp Mini硅胶膜柱用于吸附DNA，扩容管用于装大体积血浆。

t)把孵育完的混合物转移至扩容管中，打开真空泵，待离心柱中的裂解液完全抽干后，关闭真空泵，打开24孔底座一侧的排气阀将压力释放到0兆帕。小心地将扩容管拆下并丢弃。

u)向QIAamp Mini硅胶膜柱中加入600μl的Buffer ACW1，关闭排气阀，并打开真空泵，进行抽滤液体。当离心柱中Buffer ACW1被抽干后，关闭真空泵，打开24孔底座一侧的排气阀将压力释放到0兆帕。

v)向QIAamp Mini硅胶膜柱中加入750μl的Buffer ACW2，关闭排气阀，并打开真空泵，进行抽滤液体。当离心柱中Buffer ACW2被抽干后，关闭真空泵，打开24孔底座一侧的排气阀将压力释放到0兆帕。

w)向QIAamp Mini硅胶膜柱中加入750μl的无水乙醇溶液，关闭排气阀，并打开真空泵，进行抽滤液体。当离心柱中无水乙醇被抽干后，关闭真空泵，打开24孔底座一侧的排气阀将压力释放到0兆帕。关闭真空泵电源。

x)盖上QIAamp Mini硅胶膜柱并从真空支管上取下后放置到干净的2ml收集管中，将VacConnector丢弃。收集管在全速条件(20,000x g；14,000rpm)下离心3min。

y)将QIAamp Mini硅胶膜柱放置到新的2ml收集管中，开盖并置于56℃条件下的金属浴上干燥10min至硅胶膜彻底干燥。

z)将QIAamp Mini硅胶膜柱取出后放置到干净的1.5ml洗脱管(试剂盒自带)中，并将使用过的2ml的收集管丢弃。

aa)向QIAamp Mini硅胶膜柱中硅胶膜的中央小心加入55μl的Nuclease-freewater。盖上管盖后在室温孵育3min。

bb)将洗脱管置于小型离心机中全速(20,000x g；14,000rpm)离心1min来洗脱cfDNA。

cc)质量标准与评估

Qubit HS定量：取1μLcfDNA使用

dsDNA HS Assay Kit定量，记录浓度。

Agilent 2100检测：测定cfDNA片段分布。

dd)实验完成后，将所有物品归位，并清洁实验台面，将生物安全柜紫外灯打开，照射30min后关闭。记录详细的实验记录。

3.cfDNA文库构建

ee)建库前准备

i.从4℃冰箱取出纯化DNA所用的磁珠(AMPureXP beads,Beckman)，室温平衡30min再使用。

ii.从-20℃冰箱内取出End Repair&A-Tailing Buffer和End Repair&A-TailingBufferenzyme mix试剂，置于冰盒上解冻，待用。

iii.将要建库的cfDNA样本名称、取样日期、DNA浓度记录在实验记录本上，并编写好编号，方便之后操作。

iv.取相应数量的200μL PCR管，写好编号(管盖和管壁都标注编号)。

v.按cfDNA建库起始量10ng≤X≤100ng标准计算每个cfDNA样本所需要的DNA溶液体积，记录在实验记录本上，并取相应的体积置于对应的200μL PCR管内。

vi.向每个200μL PCR管内加入适量的Nuclease-Free water，使终体积达到50μL。

vii.注：在建库过程中配制所有反应体系应遵循如下规则：若样本少于四个，不需配制混合体系，每个样本独立加入反应体系中的每种成分溶液；若超过四个样本，则将反应体系中每个成分溶液按所需用量的105％配制混合体系，然后逐一加入各个样本中。

ff)末端修复&加A

i.按照表2所示，配制末端修复&加A反应体系。

表2：

ii.向每个200μl PCR管内加入10μL上述末端修复反应体系，混匀后低速离心，设定PCR仪，程序如下表3。

表3：

iii.将反应体系从PCR仪中取出，放置在小黄板上，并进行接头连接反应。

gg)接头连接反应体系

i.按照表4所示，配制接头连接反应体系。

表4：

成分	1个反应体系	8个反应体系(过量5％)
			PCR-级水(PCR-grade water)	5μL	42μL
连接缓冲液(Ligation Buffer)	30μL	252μL
			DNA连接酶(DNA Ligase)	10μL	84μL
总体积(Total volume)	45μL	378μL

ii.向每个反应管中加入45μL上述反应体系，温和混合均匀，低速离心。

iii.根据input DNA量加入适量的adapter，具体DNA:adapter如下表5，每个反应管各加入5μL adapter。另外根据测序要求，每个样本加入不同的adapter，使得同一个lane中不会出现两个样本使用同一个adapter的情况，记录好每个样本使用的adapter信息。

表5：

iv.混合均匀，并放入PCR仪中，设定温度20℃，反应15min。

hh)DNA纯化

i.配制80％乙醇(例如配置50mL 80％乙醇：40mL无水乙醇+10mL Nuclease-freeWater)，80％乙醇应现用现配。

ii.准备相应数量的1.5mL样本管，并做好相应的标记。

iii.将事先在室温平衡好的磁珠充分震荡混匀，并向每个管中分装88μL。

iv.将上述加了adapter的DNA与磁珠混匀。室温静置10min。

v.将1.5mL样本管置于磁力架上，进行磁珠吸附，直至溶液澄清。

vi.小心移除上清液，再加入200μL 80％乙醇，将样本管水平旋转360度，静置30s后弃上清液。(此过程，离心管一直保持在磁力架上。)

vii.重复步骤上述步骤一次。

viii.应将所有残留的酒精溶液移除。打开管盖，常温下干燥磁珠，挥发乙醇，以免过多乙醇影响后续反应体系中酶的效果。注意：不可过分干燥磁珠，否则会导致DNA不容易从磁珠上洗脱下来，造成产量损失。当磁珠表面不再有光泽时即为干燥完成。

ix.每个样本管内加入21μL Nuclease-Free water，重悬浮磁珠，充分混匀后室温静置5min。

x.准备一批新的200μL PCR管，管盖管壁标注对应的样本编号。

xi.将样本管置于磁力架，进行磁珠吸附，直至溶液澄清后，将上清液转移至对应编号的PCR管中，作为PCR实验的模板。

ii)文库扩增

i.按照表6所示，配制文库扩增反应体系。

表6：

ii.每个0.2mL样本管内加入30μL Pre-PCR扩增反应体系，温和混匀并低速离心，放入PCR仪中反应。

iii.将PCR仪设定如下程序，PCR cycle应根据input DNA量适当调整，见表7。

表7：

iv.循环数选择参考表格8。

表8：

Input DNA量(ng)	PCR cycle
		X>50ng	4
25ng<X≤50ng	5
		10ng<X≤25ng	6
X≤10ng	7

v.Pre-PCR反应结束后，开始进行文库纯化。

jj)文库纯化

i.准备相应数量的1.5mL样本管，并做好相应的标记。

ii.将事先在室温平衡好的磁珠充分震荡混匀，并向每个管中分装50μL。

iii.将上述加了adapter的DNA与磁珠混匀。室温静置10min。

iv.将1.5mL样本管置于磁力架上，进行磁珠吸附，直至溶液澄清。

v.小心移除上清液，再加入200μL 80％乙醇，将样本管水平旋转360度，静置30s后弃上清液。(此过程，离心管一直保持在磁力架上。)

vi.重复步骤上述步骤一次。

vii.应将所有残留的酒精溶液移除。打开管盖，常温下干燥磁珠，挥发乙醇，以免过多乙醇影响后续反应体系中酶的效果。注意：不可过分干燥磁珠，否则会导致DNA不容易从磁珠上洗脱下来，造成产量损失。当磁珠表面不再有光泽时即为干燥完成。

viii.每个样本管内加入35μL Nuclease-Free water，重悬浮磁珠，充分混匀后室温静置5min。

ix.准备一批新的离心管，管盖上标注所属项目，取样日期，样本名称；管壁上标注接头信息，建库日期，浓度。

x.将1.5mL样本管置于磁力架上，进行磁珠吸附，直至溶液澄清后，将上清液转移至对应的新的写有样本信息的1.5mL离心管。

xi.取1ul样本测浓度，1ul样本使用Agilent 2100测定文库片段大小，并记录相应信息。

xii.样本放入相对应项目的冻存盒内，置于-20℃保存。

xiii.实验完成后，将所有物品归位，并清洁实验台面，将超净工作台紫外灯打开，照射30min后关闭。记录详细的实验信息。

4.文库pooling

kk)准备好实验所需的仪器、试剂、耗材。

ll)按照测定的浓度以及所需要测定的数据量，计算pooling体积。

mm)取一个新的1.5ml离心管，做好标记。按照计算的pooling体积进行pooling。

nn)混合均匀之后，测定浓度，并记录信息。

oo)实验完成后，将所有物品归位，并清洁实验台面。

5.上机测序

将上述pooling好的文库用Tris-HCl以及NaOH进行稀释变性，然后进行上机测序。

实施例2

一、按照实施例1的方法，完成对样本的建库测序，获得下机数据，过滤掉低质量等reads后，使用比对软件(bwa)将这些测序reads比对到人的参考基因组上(hg19)。

二、对比对后的bam文件进行统计：

1.过滤只有一端比对的reads和重复reads(用samtools或者picard标记的重复reads)

2.统计低质量比对率：低质量比对即为任意一端reads的比对质量值小于30的比对reads，将这些reads全部加起来，除以总的reads数目，即为低质量比对率；

3.对比对质量值高(大于30)的异常比对，统计支持结构变异的reads信号强度：包括1)如下图4所示，3种比对方向异常的读段(reads)的数目统计；正常情况下reads的两端，一段比对到正义链(F:forward)，一段比对到反义链(reverse)，且比对到正义链的比对位置，小于比对到反义链的比对位置，俗称FR比对。因此，当两条reads同时比对到正义链(FF)，同时比对到反义链(RR)，以及一端比对到反义链的位置小于另一端比对到正义链的位置(RF)，这三种比对都是可能是因为结构变异引起的，命名为“FF_RR_RF”。2)除了上面的所述，如图2所示：reads的两端比对到不同染色体上，或者比对到同一条染色体，但是两端的距离>10000碱基，正常的距离是小于1000bp，命名为“Diff_Chr_Pos_PE”。3)如图1所示：一端的reads，从中间分开成两部分比对,分别比对到基因组上的两个不同位置(不同染色体，或者两段距离>10000bp)，这种比对称为“剪接比对”。另外一条reads的正常比对到人的染色体上，且与剪接比对的任意一部分相距不超过1000bp,命名为“PE_Sclip”；4)如图3所示：reads的两端比对起始位置接近(<150),且两条reads都是剪接比对，同时两条reads的剪接断点是在同一个位置，命名为“TwoSpliceMap”。按照上面所述的统计每种异常reads的数目，除以样本总的比对reads数目得到的ratio，为染色体结构异常信号强度值。基于实施例1，发明人分析进400例癌症和正常样本，使用t检验，发现染色体结构变异信号强度在癌症样本的强度，显著高于正常样本(P-value值<0.05)，结果如图7所示。本实施例中样本的染色体异常信号值，见表9。

表9：

4.对于比对质量高(>30)的正常比对reads，统计插入片段长度(正常比对到染色体上reads两端的距离)分布。很多研究报道，来源癌症肿瘤细胞的游离DNA片段的插入片段长度，小于来源正常细胞的游离DNA片段的插入片段长度。如图8所示，某一实施例样本的插入片段分布图。根据这个分布不，发明人统计插入片段在70～150,180～220,250～300之间的比例，记为P70_150,P180_220,P250_300。同时，通过比较发现P70_150与P180_220之间存在很强的负相关性。因此，选择其中的一个作为癌症和正常的预测特征P180_220(如图9：癌症和正常样本之间的t检验，P180_220的p-value更小，)。同时，如图8所示，在小于150bp的部分，存在小的波峰和波谷(图中箭头所示)，不同样本都存在同样的线下，因此，发明人统计每一个次高峰(峰值对应的插入片段长度:81 92 102 112 122 134)与对应的波谷之前的差值(波谷对应的插入片段长度：84 96 106 116 126 137。将6个差值加起来，命名为“峰谷间距”；

5.基于参考文献(Jiang,P.,et al.,.Proc Natl Acad Sci U S A,2015.112(11):p.E1317-25.)所述，癌症样本的线粒体的DNA拷贝数，显著性高于正常样本含有的线粒体DNA拷贝数。发明人用正常比对到人的线粒体上的reads数目，除以总的比对reads数目，作为衡量线粒体DNA的含量；近400例样本的统计结果如下表10。

表10：

6.同时发明人将整个基因组均匀的分成100kb大小的区域(bins)，统计每个区间插入片段长度在100到150bp的reads数目，记为“短片段数目”，同时统计每个区间的插入片段在151到220bp的reads数目，记为“长片段数目”。考虑到每个区域的GC含量和比对率(Mappability)不一样，因此发明人使用局部加权非参数回归参数(loess)分别对短片段数目和长片段数目进行校正。具体过程如下：

a)bins的过滤包括：1)mappability>0.6；2)N的比例<0.5；3)不在从UCSC上下载的region文件wgEncodeDacMapabilityConsensusExcludable.bed和wgEncodeDukeMapabili tyRegionsExcludable.bed；4)过滤掉X，Y染色体；

b)根据每个bin的GC值：统计每个窗口(bin)内A，T,C,G碱基的数量；以及G和C的数量。GC所占的比值，为该窗口的GC含量。

c)Mappability计算：根据从UCSC下载的ENCODE’s mappability bigwig文件，将文件中的每个region的mappability与bin比较，计算出每个bin里面所有region的mappability的平均值，作为该bin的mappability值。

d)每个区间的数目，相对于bins的长度校正(除以该bin非N的比例)；

e)将每个bin的GC和mappability组合，并按照它们的组合进行分组，同时计算每个GC和mappability组合对应所有bins的reads数目中位数。

f)使用广义交叉验证的方法(loess)，构建GC和mappability相对于长片段或者短片段的数目的拟合曲线。最后针对每个bin，根据其对应的GC含量和mappability，以及上面拟合的曲线，计算出该区域对应的理论片段数目，用该区间统计到的片段数目减去理论片段数目得到片段数目的残差。

g)使用该样本的长片段或者短片段数目的中位值加上残差值，作为该区域最后的校正值；并将相邻片段加起来，最终计算出每5M一个区域的长片段数目校正值和短片数目校正值；

h)对每个5M区间的片段数目进行过滤，要求过滤掉区间里面片段数目显著性小于3倍标准差的区间，最终得到537个5M区间。

i)对于过滤后的每个区间，用短片段数目除以长片段数目得到每个区间的片段比值。如下表11所述部分样本，部分区间的片段比值。

表11：

7.将上面所有样本在每5M区间的片段比值进行主成分分析(PCA)，对数据进行降维处理(见图10)，要求累积方差>80％，因此发明人选择前10个主成分。实施例样本(部分)对应的前10个主成分的值，如表12。

表12：

实施例3：

基于实施例2得到：(1)支持染色体结构变异的信号强度；(2)样本的线粒体含量比率；(3)整个样本插入片段在180～220和250～300之间的比例，以及小于150bp的波峰与波谷之间的“峰谷间距”。(3)每5M区间的短片段与长片段数目的比值，经主成分分析降维后，取前10个主成分的值。

将样本的这些统计值作为特征向量输入，使用机器学习的方法(比如：SVM，Lasso，GBM)，并基于上面近400例癌症和正常样本，采用10折交叉验证，测试肿瘤预测的效果。将样本平均分成10分，依次利用其中的9份数据作为训练集，建立肿瘤预测模型。剩余的一份作为训练集，用来衡量模型预测效果。并计算出对于每个测试集的AUC值(定义为ROC曲线下与坐标轴围成的面积)，详见图11.其中黑色的粗线表示10折交叉验证的均值。其中平均AUC值为0.864(LASSO的10折交叉验证AUC平均值是0.871，SVM是0.867)，其中最大的AUC值为0.932；

基于上面选择的模型，构建预测模型，对所有样本(包括部分第三方独立验证样本)进行肿瘤预测，确定所有样本来源肿瘤的概率。最终基于ROC曲线，取98％特异性下对应的p-value值作为cut-off值：0.43。

同理，发明人也可以基于对正常人和手术后的肿瘤病人进行实施例1和实施例2的统计结果，基于机器学习构建肿瘤复发预测模型。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。

Claims (25)

1.一种确定样本来源的方法，其特征在于，包括：

(1)将样本染色体的测序读段数据与参考基因组进行比对，确定样本染色体的测序读段数据的低质量比对率和高质量比对率；

(2)基于高质量比对率所对应的读段，确定样本染色体的结构变异的信号比例、线粒体DNA的含量、预定长度插入片段的比例；

(3)基于预定肿瘤预测模型，以及步骤(2)所获得的染色体的结构变异的信号比例、线粒体DNA的含量、预定长度插入片段的比例，确定样本来源的概率；

(4)基于样本来源的概率，确定样本来源，

其中，所述染色体的结构变异的信号比例是通过如下公式确定的：

y＝x/A,

其中，y表示染色体的结构变异的信号比例，

x表示异常比对测序读段的数目，

A表示样本染色体的测序读段的总数目；

所述预定长度插入片段的比例包括所述预定长度插入片段的读段数在样本测序读段的总数目的比例和样本染色体预定区间的不同预定长度插入片段的读段数的比值。

2.根据权利要求1所述的方法，其特征在于，所述测序读段数据的低质量比对率是通过如下方式确定的：

确定任意一端测序读段的比对质量值小于30的测序读段的数目；

基于任意一端测序读段的比对质量值小于30的测序读段的数目与测序读段的总数目的比值，确定测序读段数据的低质量比对率。

3.根据权利要求1所述的方法，其特征在于，所述测序读段数据的高质量比对率是通过如下方式确定的：

确定测序读段两端的比对质量值均大于30的测序读段的数目；

基于测序读段两端的比对质量值均大于30的测序读段的数目与测序读段的总数目的比值，确定测序读段数据的高质量比对率。

4.根据权利要求1所述的方法，其特征在于，所述异常比对测序读段是指：基于一对测序读段，所述一对测序读段包括预定序列的正向读段和反向读段，所述正向读段的一部分序列与第一染色体的第一预定参考序列完全匹配，另一部分与第二预定参考序列完全匹配，所述第二预定参考序列与所述第一预定参考序列位于同一染色体，并且二者的距离大于10000bp或所述第二预定参考序列位于不同于第一染色体的第二染色体上，所述反向读段与第三预定参考序列完全匹配，并且第三预定参考序列与第一预定参考序列或者第二预定参考序列在同一染色体上，且位置距离小于1000bp。

5.根据权利要求1所述的方法，其特征在于，所述异常比对测序读段是指：基于一对测序读段，所述一对测序读段包括预定序列的正向读段和反向读段，所述正向读段与第一染色体的第一预定参考序列完全匹配，所述反向读段与第二预定参考序列完全匹配，所述第二预定参考序列与所述第一预定参考序列位于同一染色体并且二者的距离大于10000bp或所述第二预定参考序列位于不同于第一染色体的第二染色体上。

6.根据权利要求1所述的方法，其特征在于，所述异常比对测序读段是指：基于一对测序读段，所述一对测序读段包括预定序列的正向读段和反向读段，所述正向读段和反向读段的一部分同时比对到第一染色体的第一预定参考序列，所述正向读段和反向读段的剩余部分同时比对到第二预定参考序列，且比对到第二预定参考序列的起始位置一样，所述第二预定参考序列与所述第一预定参考序列位于同一染色体并且二者的距离大于10000bp或所述第二预定参考序列位于不同于第一染色体的第二染色体上。

7.根据权利要求1所述的方法，其特征在于，所述异常比对测序读段是指：基于一对测序读段，所述一对测序读段包括预定序列的正向读段和反向读段，所述正向读段或反向读段分别与第一预定参考序列的正义链完全匹配和第二预定参考序列的正义链完全匹配；或者分别与第一预定参考序列的反义链完全匹配和第二预定参考序列的反义链完全匹配；或者正向读段与第一染色体的第一预定参考序列的正义链完全匹配，反向读段与第一染色体的第二预定参考序列的反义链完全匹配，第二预定参考序列在第一染色体上的位置小于第一预定参考序列。

8.根据权利要求1所述的方法，其特征在于，所述线粒体DNA的含量是通过如下方式确定的：

确定比对到参考线粒体基因序列的测序读段的数目；

基于比对到参考线粒体基因序列的测序读段的数目与测序读段的总数目的比值，确定线粒体DNA的含量。

9.根据权利要求1所述的方法，其特征在于，确定预定长度插入片段的读段数在样本染色体的测序读段的总数目的比例是通过如下方式进行的：

将所述预定长度插入片段的长度设置为70～150bp，150～180bp，180～220bp或250～330bp。

10.根据权利要求1所述的方法，其特征在于，进一步包括确定峰谷间距，所述峰谷间距是指插入片段在小于150bp范围内，波峰与波谷对应长度的插入片段及其上下2bp[±2bp]的测序读段数目占样本染色体的测序读段的总数目的比例的差值。

11.根据权利要求1所述的方法，其特征在于，确定样本染色体预定区间的所述不同预定长度插入片段的读段数的比值是通过如下方式进行的：

将参考基因序列划分为多个相同长度的窗口区间；

确定每个窗口区间内不同预定长度插入片段的读段数目；

确定每个窗口区间内不同预定长度插入片段的读段数目的比值。

12.根据权利要求11所述的方法，其特征在于，所述窗口区间的大小为100kb。

13.根据权利要求11所述的方法，其特征在于，所述预定长度插入片段的长度为100～150bp或151～220bp。

14.根据权利要求11所述的方法，其特征在于，在每个窗口区间内，进一步包括对预定长度插片片段的读段数目进行校正处理。

15.根据权利要求14所述的方法，其特征在于，所述校正处理是通过将在每个窗口区间内预定长度的插入片段的读段数目的中位值加上片段数目残差获得的。

16.根据权利要求15所述的方法，其特征在于，所述片段数目的残差是通过如下方式获得的：

(1)确定每个窗口区间内的GC含量和比对率；

(2)将步骤(1)所获得的每个窗口区间内的GC含量和比对率进行组合和分组处理，获得每个GC含量和比对率组合对应窗口区间的测序读段数目的中位值；

(3)基于局部加权非参数回归方法，构建GC含量和比对率组合对应窗口区间的测序读段数目中位值相对于GC含量和比对率的拟合曲线；

(4)基于所述拟合曲线以及每个窗口区间内的GC含量和比对率，确定每个窗口区间内的理论插入片段数目；

(5)将每个窗口区间内的预定长度的插入片段的读段数目减去步骤(4)所获得的理论插入片段数目，获得每个窗口区间内的预定长度的插入片段数目的残差。

17.根据权利要求11所述的方法，其特征在于，进一步包括将预定区间内的预定长度插入片段的数目进行加和处理，所述加和处理包括分别将插入片段的长度为100～150bp的插入片段的读段数目进行加和和将插入片段的长度为151～220bp的插入片段的读段数目进行加和。

18.根据权利要求17所述的方法，其特征在于，所述加和处理后的区间的长度为5M。

19.根据权利要求17所述的方法，其特征在于，进一步包括将插入片段的长度为100～150bp的插入片段的读段数目加和除以插入片段的长度为151～220bp的插入片段的读段数目加和，以便获得插入片段读段数目加和的比值。

20.根据权利要求19所述的方法，其特征在于，进一步包括对所述加和的比值进行主成分分析。

21.根据权利要求1所述的方法，其特征在于，所述预定肿瘤预测模型包括选自SVM、Lasso、GBM的至少之一。

22.根据权利要求1所述的方法，其特征在于，所述预定肿瘤预测模型为GBM，基于ROC曲线、预定的敏感性或特异性，确定相应阈值。

23.根据权利要求22所述的方法，其特征在于，所述预定的特异性为98％，所述阈值为0.45。

24.根据权利要求22所述的方法，其特征在于，在步骤(4)中，所述样本来源肿瘤的概率大于0.45，是样本来源于肿瘤患者的指示。

25.一种确定样本来源的系统，其特征在于，包括：

比对模块，所述比对模块用于将样本染色体的测序读段数据与参考基因进行比对，确定样本染色体的测序读段数据的低质量比对率和高质量比对率；

获取模块，所述获取模块用于基于高质量比对率所对应的读段，确定样本染色体的结构变异的比例、线粒体DNA的含量、预定长度插入片段的比例；

确定样本来源概率模块，所述确定样本来源概率模块用于基于预定肿瘤预测模型，以及获取模块所获得的染色体的结构变异的比例、线粒体DNA的含量、预定长度插入片段的比例，确定样本来源的概率；

判定模块，所述判定模块用于基于确定样本来源概率模块所确定的样本来源的概率，确定样本来源；

其中，所述染色体的结构变异的比例是通过如下公式确定的：

y＝x/A,

其中，y表示染色体的结构变异的比例，

x表示异常比对测序读段的数目，

A表示样本染色体的测序读段的总数目；

所述预定长度插入片段的比例包括所述预定长度插入片段的读段数在样本测序读段的总数目的比例和样本染色体预定区间的不同预定长度插入片段的读段数的比值。

Images