CN109504751B - 一种肿瘤复杂克隆结构的缺失变异识别及克隆计数方法 - Google Patents

Abstract

本发明公开了一种肿瘤复杂克隆结构的缺失变异识别及克隆计数方法，首先运用缺失识别方法（如DELLY方法）确定与缺失突变相关的所有区域范围；然后针对每一个缺失区域，采用分段扩展和最长可变拆分读段的组合策略识别缺失断点，得到飘移缺失断点候选集合；接着根据断点对的数目区分飘移缺失与固定缺失，针对固定缺失，采用贪婪策略识别缺失周围混乱的单核苷酸变异位点并逐级分离各亚克隆所属的变异位点；最后获得克隆数目并依据克隆数目信息进一步优化飘移缺失识别结果。本发明（项目编号20180550161、20180550855）不仅能够识别克隆演变中的飘移缺失突变，而且能够识别克隆分化引起的单核苷酸混乱变异情况，对克隆模式的研究提供了另一种检测途径。

Description

一种肿瘤复杂克隆结构的缺失变异识别及克隆计数方法

技术领域

本发明属于面向精准医学的数据科学技术领域，涉及一种肿瘤复杂克隆结构的缺失变异识别及克隆计数方法。

背景技术

近年来，肿瘤基因组计划(Cancer genome project)、肿瘤基因组图谱计划(TheCancer Genome Atlas，TCGA)、国际肿瘤基因组协作联盟(International cancer genomeconsortium，ICGC)及精准医疗计划等大型项目地实施，极大地推动了肿瘤基因组学的研究与发展。此外，高通量测序技术下产生的海量组学大数据为研究人类基因功能、疾病机制与精准医疗等方面奠定了研究基础，同时也为肿瘤基因组学的研究与发展提供了必要的数据支撑与保障。

肿瘤是基因突变逐步累积而形成的。基因突变是DNA序列的大分子中某些小分子或小分子片段发生变异，导致肿瘤的产生与演变，因此肿瘤中的基因突变信息是进行肿瘤研究的重要依据，也是揭示复杂疾病遗传机制的关键。其中，缺失是一种常见且功能重要的结构变异，可引起严重的基因疾病。如，较大缺失突变因其通常有基因缺失而致命；中等长度的缺失突变可导致一些人类疾病，如威廉斯氏综合征；研究发现前列腺癌、乳腺癌、卵巢癌、结肠癌和肝癌等在8号染色体上有一些较小的重叠缺失区域。肿瘤具有高度的异质性和遗传多样性，治疗中不断出现新的选择压力，从而分化出更小亚克隆产生抗体使病情持续发展。缺失突变在不同的克隆中也呈现多样性。

现有的缺失突变识别方法包括读段深度法、双末端读段匹配法、拆分读段法和从头组装法。其中，读段深度法，如CNVnator，通过较少读段深度的片段来识别缺失区域；双末端读段匹配法，如HYDRA，通过插入距离大于期望值的双末端读段信息识别缺失区域；拆分读段法，如CREST和Pindel等，通过跨越断点的读段信息识别缺失；从头组装法采用拼接所有短读段形成一条序列，再与参考序列相比对来识别缺失，如EULER-USR、SvABA和Platypus等。此外，基于以上方法的一些组合方法相继发展，如DELLY、PRISM、LUMPY、SV-Bay和SWAN等。然而，现有传统方法都没有针对性地考虑肿瘤异质性问题，所以都不适用于肿瘤复杂克隆结构中的缺失突变识别，具体原因如下。

在克隆演化过程中，对于同一缺失突变，可能存在三种情况：其一，缺失突变在不同的亚克隆中呈现较小差异性，包括长度、位置等信息的差异，此种情况称为飘移缺失，会产生许多位置相近的断点。其二，缺失突变在不同的亚克隆中自身保持较为稳定的状态，但在缺失周围会出现许多单核苷酸变异位点，且随着克隆分化变得更加复杂，使得缺失突变识别变得更加困难。此种情况称为固定缺失。其三，缺失突变及周围区域均呈现较为稳定状态，此种情况称为简单缺失，此种情况在本方法中不做考虑。

对于第一种情况，大多数传统缺失识别方法针对同一缺失区域并未考虑多种不同缺失情况存在的现象，因此仅报出一个结果。如，DELLY方法只寻找一个最大团，SWAN仅检测缺失区域的最外边界。除此之外，尽管有些方法能报出多个断点，但存在着一定局限性，如Pindel所采用的模式增长策略依赖于锚点识别，锚点识别的准确性又受测序误差与比对偏差所影响，从而产生不同程度的假阳性和假阴性识别结果。对于第二种情况，由于传统识别方法常采用拆分读段法，而被拆分读段的两部分由于受到混乱的单核苷酸变异位点影响而不能准确地映射到正确位置，从而进一步影响缺失断点的识别精度。此外，①传统识别方法中常采用一些过滤原则来减少由测序和比对误差而产生的假阳性结果。但在肿瘤复杂克隆结构中，过滤原则会将占比较低的亚克隆读段信息直接删除；②混乱的单核苷酸变异位点会引起缺失周围区域读段的比对质量下降，支持缺失断点信息的携带变异的读段将被过滤掉。

综上所述，由于肿瘤异质性和复杂克隆结构，已有的传统缺失识别方法只适用于简单缺失识别，不能较好地适用于复杂克隆结构中的飘移缺失和固定缺失识别。然而，肿瘤精准诊疗的现实需求需要充分地考虑肿瘤异质性中的各级亚克隆引起的复杂情况，才能为临床提供有效地指导策略。

发明内容

本发明针对上述现有技术中的不足，基于对肿瘤结构异质性问题的考虑，提供一种针对肿瘤复杂克隆结构的缺失突变识别方法，并根据识别出的缺失突变信息及其周围区域的单核苷酸变异信息实现克隆数目统计。

本发明是通过以下技术方案来实现：

一种肿瘤复杂克隆结构的缺失变异识别及克隆计数方法，包括以下操作：

1)以肿瘤基因测序结果作为读段比对信息，运用DELLY方法确定克隆结构中与缺失突变相关的所有区域范围，构造缺失参考序列集合D_r；

同时为每个缺失区域筛选出与之相关的候选双末端读段，构造候选读段集合R_r；

2)针对每一个缺失区域，采用分段扩展和最长可变拆分读段的组合策略识别缺失断点，得到飘移缺失断点候选集合；

其中，通过分段扩展策略确定R_r中各缺失区域R_ri的读段正确比对位置，排除克隆演变中可能引起的可变缺失断点；最长可变拆分读段采用可容错的动态规划算法寻找最大公共子串，以此确定比对位置与缺失断点；

3)根据断点对数目区分飘移缺失与固定缺失；对于区分出的固定缺失，采用贪婪策略识别缺失周围混乱的单核苷酸变异位点并逐级分离各亚克隆所属的变异位点；

4)获得克隆数目，并依据克隆数目信息优化飘移缺失识别结果；

所述的优化是对识别出的飘移缺失区域进行过滤，若一个区域的飘移缺失个数大于识别出来的子克隆数目k，则只输出前k个飘移缺失。

所述的飘移缺失的识别为：

1)以DELLY方法检测的缺失区域为参考，将其对应的参考基因组片段作为核心，构造缺失参考序列；对于第i个缺失，假设缺失的左断点为

右断点为

分别向两侧扩展缺失区域并确定缺失参考序列d_ri，范围区间为

其中，插入间距的均值μ和标准差σ由Picard计算所得，读段长度为l；

令D_r表示缺失参考序列集合，D_r＝{d_r1,…,d_ri,…,d_rn}，其中n为缺失突变数目，d_ri表示第i个缺失参考序列；

2)将所有读段与参考基因组序列进行比对，获得BAM/SAM格式的比对结果文件，再在比对结果文件中进一步搜索获得与每个缺失相关的候选读段；识别每个缺失区域的候选双末端读段必须满足至少有一端比对结果的CIGAR标识中包括除“M”以外的字符；

令R_r表示候选读段集合，R_r＝{R_r1,…,R_ri,…,R_rn}，其中R_ri表示第i个缺失区域对应的候选读段集合；

3)使用分段扩展策略实现R_ri的重比对：

(1)局部缺失参考序列划分成固定k碱基长度的小片段，称为k-mers；相邻的k-mers间满足(k-1)个碱基重叠并将所有的k-mers存储于哈希表中；k-mers作为哈希表中的关键字，并设置k的默认值，每个k-mer在参考基因组中的比对位置作为键内容；

(2)读段划分为无覆盖等长的k-mers，作为哈希表的键；所有读段的k-mers映射于种子哈希表中且返回每个k-mer完美比对的位置列表信息；并通过种子间相邻位置关系过滤掉比对位置结果中的部分假阳性结果；

4)采用最长可变拆分读段方法识别缺失断点：

对于第i个缺失参考序列，如果在R_ri中双末端读段的两端均能映射到d_ri中，且至少一端读段的汉明距离大于2bp，则将具有最大距离的一端读段存储于集合B_d中，集合B_d用于保存满足条件的读段；

采用最大公共子串策略确定比对位置与断点，最大公共子串指读段被拆分的部分与参考序列具有最大的公共子序列；采用具有容错机制的动态规划算法识别最大公共子串；

读段在缺失参考序列中的比对位置首先由分段扩展重比对方法确定，然后分别对读段拆分的头尾两部分进行定位；对于头部分，在读段与缺失参考序列中起始位置为比对位置的等长序列之间寻找最大公共子串，直到容错率不满足为止；对于尾部分采用相同的方法，在缺失参考序列中寻找尾部的起始位置；在容错上，考虑满足头部容错数且头与尾部的总容错数范围为[0,5bp]；

统计每个断点对的读段支持数，并从候选缺失集合中删除支持数小于2的断点对；在统计时，基于克隆数目确定飘移缺失数目的上限并调整过滤规则，过滤掉由错误比对引起的缺失；还将断点对按读段支持数升序存储，当一条读段支持多个断点时，仅将此读段统计到具有最高支持数的断点对中，并相应地修改其它断点对的支持数；

识别所有缺失断点后，得到缺失集合D_s，

m为缺失数目；

为缺失片段的左断点；

为缺失片段的右断点；

如果m>1，为飘移缺失，否则为固定缺失；若飘移缺失数目超过亚克隆数目n，则以支持数最高的n个断点对作为最终识别结果。

所述的逐级分离策略识别每个亚克隆中的单核苷酸变异位点SNVs，包括以下操作：

1)按照断点对数目区分固定缺失与飘移缺失，再根据固定缺失周围的SNVs数目区分固定缺失与简单缺失；设定SNVs数目的阈值为σ_s；

2)构建每个固定缺失的SNV参考序列，它包括

和

两个子序列，len值的设定必须满足一定的长度，由此确保获得足够的读段信息；

候选读段集合R_r中的读段依次与对应的SNV参考序列重比对；如果双末端读段的两端均与SNV参考序列相匹配，将具有最大汉明距离的一端存储于集合B_s中；否则，当只有一端读段与SNV参考序列相匹配时，将不匹配的另一端存储于集合O_s中；

3)使用B_s集合中的读段来校正D_s中的断点；采用最长可变拆分读段方法确定读段在每个候选SNV参考序列中的位置，最可信的断点对由具有最大读段支持数的候选SNV参考序列报出；

4)在校正断点之后使用集合O_s和B_s中的读段识别SNV；

对于集合O_s首先通过双末端读段的插入距离确定比对位置，插入距离服从均值为μ，标准差为σ的正态分布，插入距离应在[μ-3σ,μ+3σ]范围内；假定读段比对位置为P_s，如果未匹配读段是正向，它的位置范围为[P_s-(μ-l)-3σ,P_s-(μ-l)+3σ]；

否则，未匹配读段是反向，它的位置范围是[P_s+(μ-l)-3σ,P_s+(μ-l)+3σ]；具有最小汉明距离的位置被选作为读段的最终比对位置，之后逐个碱基比对识别SNVs；

5)将固定缺失周围的SNVs分离到不同的亚克隆中，其中SNV的支持数必须大于设置的阈值，将包含混乱SNVs的区域命名为S。

所述分离SNVs的方法如下：

1)按比对位置升序存储S中的所有SNVs；

2)将每个读段中的SNVs与S中的SNVs相比较，识别出最后一代亚克隆的读段，如果一条读段包含所有SNVs，则其必定来自于最后一代亚克隆；

3)移除最后一代克隆中的读段并修改S；如果读段移除后，支持数小于阈值的SNV均被视为最后一代亚克隆；

4)使用2)估计移除的读段数目，如果此值小于阈值，将其直接删除；

5)迭代执行2)到4)，直至S为空或执行次数大于假设的最大克隆数为止，实现克隆计数，获得克隆数目。

与现有技术相比，本发明具有以下有益的技术效果：

本发明方法主要针对复杂克隆结构中缺失突变及其周围的混乱SNVs进行识别，并通过分离SNVs实现克隆计数；显著提高了缺失变异识别能力，有益于理解克隆发生与演变机制。

本发明方法首先使用DELLY方法确定缺失区间，提高了算法效率，进一步采用分段扩展重比对策略确定读段比对位置，使用哈希表实现读段的快速存储与查找，并引入相邻过滤思想快速筛选出错误比对位置，减少了验证过程的执行步骤，极大地提高了算法运行速度并提高了种子比对精度。

进一步的，本发明方法能够充分解决肿瘤异质性中的复杂情况，采用最长可变拆分读段方法识别飘移缺失断点，采用可容错的动态规划算法，通过寻找最大公共子串分别对拆分读段的头尾两部分确定比对位置，并采用过滤规则及读段支持数降低了假阳性识别结果。

进一步的，本发明方法解决了克隆演化过程中变异位点愈加复杂的情况，现有方法没有考虑多级克隆中的SNVs识别问题；本方法以不同克隆中SNVs分布具有不同特性的思想为出发点，自下而上逐级分离克隆中的SNVs，实现了克隆数目的精确统计。

总之，本方法适用于肿瘤异质性中的复杂克隆情况，不仅能够识别克隆演变中的飘移缺失突变，而且能够识别克隆分化引起的单核苷酸混乱变异情况，对克隆模式的研究提供了另一种检测途径，对临床诊疗与用药提供指导策略。

附图说明

图1为本发明方法流程示意图；

图2-1a～图2-4c为本发明方法(SV-Del)识别飘移缺失的仿真实验结果与现有方法识别结果的对比图；其中，图2-1a～图2-1c为识别不同缺失长度的对比结果，图2-2a～图2-2c为识别不同飘移区间的对比结果，图2-3a～图2-3c为识别不同覆盖度的对比结果，图2-4a～图2-4c为识别不同亚克隆数目和占比的对比结果；

图3-1a～图3-3为本发明方法(SV-Del)识别单核苷酸变异位点的仿真实验与现有方法的对比结果图；其中，图3-1a～图3-1b为不同亚克隆数目的对比结果图，图3-2为不同覆盖度的对比结果图，图3-3为不同变异频率的对比结果图；

图4a～图4c为本发明方法分离SNVs及克隆计数的实验结果图，其中，图4a为不同亚克隆数目的实验结果图，图4b为不同覆盖度的实验结果图，图4c为不同变异频率的实验结果图。

具体实施方式

下面结合附图对本发明做进一步详细描述，所述是对本发明的解释而不是限定。

本发明作为以下2018年度辽宁省自然科学基金计划项目的科技成果：

(1)项目编号：20180550161，项目名称“肿瘤克隆模式检测及突变易感性关联分析方法研究”；

(2)项目编号：20180550855，项目名称“MicroRNA-223-3调节帕金森神经炎症的相关机制研究”。

本发明基于高通量的二代测序数据，针对肿瘤中的复杂克隆结构提出了一种缺失突变的启发式识别方法。为了提高算法效率，首先运用缺失识别方法(如DELLY方法)确定与缺失突变相关的所有区域范围；然后针对每一个缺失区域，采用分段扩展和最长可变拆分读段的组合策略识别缺失断点，得到飘移缺失断点候选集合；接着根据断点对的数目区分飘移缺失与固定缺失；针对固定缺失，采用贪婪策略识别缺失周围混乱的单核苷酸变异位点并逐级分离各亚克隆所属的变异位点；最后获得克隆数目并依据克隆数目信息进一步优化飘移缺失识别结果。

本发明公开的一种针对肿瘤复杂克隆结构的缺失变异识别及克隆计数方法，包括以下操作：

1)以肿瘤基因测序结果作为读段比对信息，运用DELLY方法确定克隆结构中与缺失突变相关的所有区域范围，构造缺失参考序列集合D_r；

同时为每个缺失区域筛选出与之相关的候选双末端读段，构造候选读段集合R_r；

2)针对每一个缺失区域，采用分段扩展和最长可变拆分读段的组合策略识别缺失断点，得到飘移缺失断点候选集合；

其中，通过分段扩展策略确定R_r中各R_ri(第i个缺失区域)的读段正确比对位置，排除克隆演变中可能引起的可变缺失断点；最长可变拆分读段采用可容错的动态规划算法寻找最大公共子串，以此确定比对位置与缺失断点；

3)根据断点对数目区分飘移缺失与固定缺失；对于区分出的固定缺失，采用贪婪策略识别缺失周围混乱的单核苷酸变异位点并逐级分离各亚克隆所属的变异位点；

4)获得克隆数目，并依据克隆数目信息优化飘移缺失识别结果：对识别出的飘移缺失区域进行过滤，若一个区域的飘移缺失个数大于识别出来的子克隆数目k，则只输出前k个飘移缺失。

进一步的，包括以下操作：

S1、输入数据集包括两部分：

S11、候选缺失区域集合

采用DELLY方法(https://www-korbel.embl.de/software.html)识别与缺失突变相关的所有断点集合，每一对断点被视为一个缺失区域。

S12、候选读段信息集合

采用BWA工具(https://sourceforge.net/projects/bio-bwa/files/)获得读段比对信息，然后为每个缺失区域筛选出与之相关的候选双末端读段，这些读段必须满足至少存在一端比对结果的CIGAR标识中包括除‘M’以外的其它字符。

S2、识别飘移缺失，具体识别方法如下：

S21、克隆演变可能引起可变的缺失断点，由此采用分段扩展重比对策略确定读段比对的正确位置。

S22、采用最长可变拆分读段方法，可容错的动态规划算法寻找最大公共子串，以此确定缺失断点。

S23、过滤假阳性断点。

S3、识别固定缺失周围的单核苷酸变异位点(SNVs)。

S31、固定缺失周围的SNVs识别。

S32、将SNVs有效分离到所属亚克隆中，实现克隆数目的统计。

S33、以克隆数目为依据进一步优化飘移缺失的识别结果。

下面分别对各个操作步骤进行详细的说明。

飘移缺失的识别

本发明方法首先采用分段扩展重比对策略确定读段比对位置。首先以DELLY方法检测的缺失区域为参考，将其对应的参考基因组片段为核心，构造缺失参考序列。对于第i个缺失，具体方法如下：

假设缺失的左断点为

右断点为

分别向两侧扩展缺失区域并确定缺失参考序列d_ri，范围区间为

其中，插入间距的均值μ和标准差σ由Picard(http://sourceforge.net/projects/picard/)计算所得，读段长度为l。

令D_r表示缺失参考序列集合，D_r＝{d_r1,…,d_ri,…,d_rn}，其中n为缺失突变数目，d_ri表示第i个缺失参考序列。

然后，搜索BAM/SAM文件获得与每个缺失相关的候选读段。识别每个缺失区域的候选双末端读段必须满足至少有一端比对结果的CIGAR标识中包括除‘M’以外的其它字符。

令R_r表示候选读段集合，R_r＝{R_r1,…,R_ri,…,R_rn}，其中R_ri表示第i个缺失区域对应的候选读段集合。

使用分段扩展策略实现R_ri的重比对，具体方法如下：

(1)局部缺失参考序列划分成固定k碱基长度的小片段，称为k-mers。相邻的k-mers间满足(k-1)个碱基重叠并将所有的k-mers存储于哈希表中。k-mers作为哈希表中的关键字，在此设k的默认值为10bp，每个k-mer在参考基因组中的比对位置作为键内容。

(2)读段划分为无覆盖等长的k-mers，作为哈希表的键。所有读段的k-mers映射于种子哈希表中且返回每个k-mer完美比对的位置列表信息。

如果仅通过k-mer的映射位置信息来确定读段的比对位置，则会产生假阳性在内的多个比对结果。故需要验证比对结果，从而极大的降低了算法效率。因此，在本发明方法中，引入相邻过滤思想将错误比对位置快速有效地筛出。相邻过滤的基本思想是：一个k-mer的位置列表中种子的可能位置必须满足在允许的编辑距离下读段中其它相邻k-mer的比对位置在参考序列中也是相邻的。此方法通过种子间相邻位置关系能够有效地过滤掉部分假阳性结果，减少了验证过程的执行步骤，提高了种子正确比对的精度。

具体地，在本发明方法中，将第一个k-mer的位置作为读段的候选位置。对于一个读段，假定第i个k-mer有多个比对位置，此读段也将有多个比对位置。在此种情况下，采用所有k-mers投票确定读段的候选位置，每个候选位置的k-mers支持数为投票数，且返回最小汉明距离允许范围内的位置作为最终位置(汉明距离表示两个序列中对应碱基位不同的数目)。存储读段候选位置信息包括至少σ_s个完美比对的k-mers，σ_s为支持数阈值，它的取值与读段长度相关，在此σ_s默认值为5。本发明方法中不考虑较小的缺失。

本发明的方法然后采用最长可变拆分读段方法识别缺失断点，具体方法如下：

对于第i个缺失参考序列，如果在R_ri中双末端读段的两端均能映射到d_ri中，且至少一端读段的汉明距离大于2bp，则将具有最大距离的一端读段存储于集合B_d中(集合B_d用于保存满足条件的读段)。

当一个读段被一个断点切分时，可在缺失参考序列中确定被拆分两部分的位置，继而获得断点对。但单核苷酸变异位点(SNVs)和测序误差影响被拆分读段两部分位置的准确检测，所以本发明方法采用最大公共子串策略确定比对位置与断点。最大公共子串指读段被拆分的部分与参考序列具有最大的公共子序列。

在此，采用具有容错机制的动态规划算法识别最大公共子串。动态规划的表达式如下：

其中，A[i]与B[j]分别表示序列A和B中的第i个和第j个碱基，k表示容错碱基数。DP[k][i][j]表示容错为k且以A[i]和B[j]为结束的最大公共子串长度。

读段在缺失参考序列中的比对位置首先由分段扩展重比对方法确定，然后分别对读段拆分的头尾两部分进行定位，对于头部分，在读段与缺失参考序列中起始位置为比对位置的等长序列之间寻找最大公共子串，直到容错率不满足为止。对于尾部分采用相同的方法，在缺失参考序列中寻找尾部的起始位置。在容错上，首先考虑满足头部容错数且头与尾部的总容错数范围为[0,5bp]。

本方法统计每个断点对的读段支持数，并从候选缺失集合中删除支持数小于2的断点对。当头/尾部片段较短时，可能会产生多个比对位置，将会获得多个缺失。但仅有一个缺失是正确的，其余均为错误比对引起的缺失应被过滤掉。为了降低此种情况的假阳性错误率，应基于克隆数目确定飘移缺失数目的上限并相应地调整过滤规则。

此外，肿瘤中有多个共存克隆，亚克隆占比使得覆盖度不同，高覆盖度的亚克隆中可能会存在一些比低覆盖度的亚克隆支持数更高的假阳性断点。此时，低覆盖度的亚克隆中识别的断点可能被过滤掉。为了解决此问题，首先将断点对按读段支持数升序存储，当一条读段支持多个断点时，仅将此读段统计到具有最高支持数的断点对中，并相应地修改其它断点对的支持数。按照相同方法处理所有的缺失，得到缺失集合D_s，D_s＝{d_sj}_j＝1:m，

m为缺失数目。如果m>1，为飘移缺失，否则为固定缺失。若飘移缺失数目超过亚克隆数目n，则输出支持数最高的n个断点对作为最终识别结果。

克隆计数

当肿瘤中存在缺失突变时，在其周围区域可能出现许多SNVs，且随着克隆演变而愈加混乱，使得缺失突变的识别变得更加困难。目前已有的SNV识别方法，如常用工具GATK，当覆盖度较低时无法精确地识别SNVs。此外，现有方法无法有效区分亚克隆中的SNVs。但是，每个亚克隆中的SNVs呈现不同的特征，它的识别更有利于实现克隆计数和缺失检测。

基于上述原因，本方法试图通过固定缺失周围的SNVs识别，并将SNVs有效分离到所属亚克隆中，实现克隆数目的统计，最后以克隆数目为依据进一步优化飘移缺失的识别结果。

本方法采用逐级分离策略识别每个亚克隆中的SNVs，具体方法如下：

(1)按照上一步检测的断点对数目区分固定缺失与飘移缺失，再根据固定缺失周围的SNVs数目区分固定缺失与简单缺失。设定SNVs数目的阈值为σ_s，在此不考虑小于阈值的简单缺失情况。

(2)构建每个固定缺失的SNV参考序列，它包括

和

两个子序列，len的值由用户设定，但其必须满足一定的长度确保获得足够的读段信息。

然后，候选读段集合R_r中的读段依次与对应的SNV参考序列重比对。如果双末端读段的两端均与SNV参考序列相匹配，将具有最大汉明距离的一端存储于集合B_s中；否则，当只有一端读段与SNV参考序列相匹配时，将不匹配的另一端存储于集合O_s中。

(3)由于DELLY报出的断点信息可能存在偏差，所以需要进一步验证并获得正确的断点。在此，使用B_s中的读段校正D_s中的断点。采用最长可变拆分读段方法确定读段在每个候选SNV参考序列中的位置，最可信的断点对由具有最大读段支持数的候选SNV参考序列报出，以此校正DELLY报出的错误断点并获得最优的SNV参考序列。

(4)在校正断点之后使用O_s和B_s中的读段识别SNV。由于集合O_s中的读段位置未知，首先通过双末端读段的插入距离确定比对位置。插入距离服从均值为μ，标准差为σ的正态分布，绝大多数的插入距离应在[μ-3σ,μ+3σ]范围内。假定读段比对位置为P_s，如果未匹配读段是正向，它的位置范围为[P_s-(μ-l)-3σ,P_s-(μ-l)+3σ]；否则，未匹配读段是反向，它的位置范围是[P_s+(μ-l)-3σ,P_s+(μ-l)+3σ]。具有最小汉明距离的位置被选作为读段的最终比对位置，之后逐个碱基比对识别SNVs。

(5)将固定缺失周围的SNVs分离到不同的亚克隆中。为了避免测序误差的干扰，SNV的支持数必须大于用户设置的阈值。将包含混乱SNVs的区域命名为S。

分离SNVs的方法如下：

步1：按比对位置升序存储S中的所有SNVs；

步2：将每个读段中的SNVs与S中的SNVs相比较，识别出最后一代亚克隆的读段，包含所有SNVs的读段必定来自于最后一代亚克隆。

步3：移除最后一代克隆中的读段并修改S。如果读段移除后，支持数小于阈值的SNV被视为最后一代亚克隆。

步4：使用步2估计移除的读段数目，如果此值小于阈值，它可能是由测序误差所引起，将其直接删除。

步5：迭代执行步2到步4，直至S为空或执行次数大于10(假设最大克隆数为10)为止。实现克隆计数，获得克隆数目。

由于存在单点变异及测序误差，识别出的飘移缺失可能会存在假阳性错误。当利用缺失周围混乱的单核苷酸变异位点判断出克隆数目后，可进一步验证和校正飘移缺失结果。即根据子克隆数目对识别出的飘移缺失区域进行过滤，若一个区域的飘移缺失个数大于识别出来的子克隆数目k，则只输出前k个飘移缺失，实现结果优化。

为了测试本发明方法SV-Del的实用性与有效性，下面通过仿真实验分别与四个较常用方法SWAN、Pindel、VarScan2和GATK进行比较分析。仿真实验中，随机选取人类参考基因组(hg19)中的几条染色体，每条染色体中均嵌入包含飘移缺失、固定缺失和简单缺失在内的三种缺失。为了避免缺失间的相互干扰，实验中限制任意两个缺失间距至少大于读段长度5个碱基的长度。使用仿真工具TNSim(http://github.com/lnmxgy/TNSim)分别生成正常基因组和肿瘤基因组的读段，并在读段中设置1％测序误差，读段长度为100bp，双末端读段的插入距离服从正常分布N(500bp,49bp)。设置默认参数：覆盖度为200×，缺失飘移区间为[0,5]bp，缺失长度为[10,1100]bp。每组实验中根据需要改变相应的参数，实验采用的评价指标包括精度、召回率和F1-measure。

1、识别飘移缺失的仿真实验

1)改变缺失长度。分别设置缺失长度变化区间为[10,50)bp、[50,100)bp和[100,1100]bp，亚克隆数目为3且所占比例为5:3:2。实验结果如图2-1a～c所示。

2)改变缺失的飘移区间。从图2-1a～c可见，GATK和VarScan2两种方法仅适用于识别较短的缺失。在以下实验中，仅与SWAN和Pindel两种方法做比较分析。设置缺失的飘移区间分别为[0,5]bp、[5,10]bp和[10,20]bp；亚克隆比例为2:3:5。实验结果如图2-2a～c所示。

3)改变覆盖度。覆盖度会影响不同比例亚克隆中的缺失识别精度，设置4个亚克隆，占比为5:3:1:1，分别在覆盖度为100×、150×、200×和250×下进行测试，实验结果如图2-3a～c所示。

4)改变克隆结构。设置克隆数和所占比例分别为3(6:3:1)、4(3:3:2:2)和5(4:2:2:1:1)，结果如图2-4a～c所示。

从上述实验结果分析可见，本发明的方法在不同参数下均表现出较好的性能，且对于克隆结构中的飘移缺失识别效果较好。

2、识别SNVs的实验

1)识别SNVs。分别在覆盖度为100×和200×下进行测试，克隆数目及占比分别为3(6:3:1)、4(3:3:2:2)和5(4:2:2:1:1)，实验结果如图3-1a～b所示。改变覆盖度为100×、150×、200×和250×，亚克隆数为4且占比为5:3:1:1，实验结果如图3-2所示。调整固定缺失周围的SNVs变异频率为1％、2％和3％，实验结果如图3-3所示。从结果可见，VarScan2和GATK的过滤规则在低覆盖度下降低了F1-Measure值，本发明的方法SV-Del识别单核苷酸变异位点的F1-Measure值不受克隆数目影响，而只与覆盖度相关。当覆盖度达到一定水平时，SV-Del均可有较好的SNV识别效果。此外，SNVs较大的变异率不同程度的影响方法识别性能，但是SV-Del比其它方法具有更好的性能。

2)分离SNVs，克隆计数。实验结果如图4a～c所示。从结果可见，克隆数目影响SNVs分离性能，当克隆覆盖度较低时，较少的读段支持数使得SNVs分离困难。当覆盖度达到一定水平后，分离性能保持稳定。变异频率越高分离精度越低。

以上结果表明本方法适用于肿瘤异质性中的复杂克隆情况，不仅能够识别克隆演变中的飘移缺失突变，而且能够识别克隆分化引起的单核苷酸混乱变异情况，对克隆模式的研究提供了另一种检测途径。

以上给出的实施例是实现本发明较优的例子，本发明不限于上述实施例。本领域的技术人员根据本发明技术方案的技术特征所做出的任何非本质的添加、替换，均属于本发明的保护范围。

Claims (7)

1.一种非医疗诊断的肿瘤复杂克隆结构的缺失变异识别及克隆计数方法，其特征在于，包括以下操作：

1)以肿瘤基因测序结果作为读段比对信息，运用DELLY方法确定克隆结构中与缺失突变相关的所有区域范围，构造缺失参考序列集合D_r；

同时为每个缺失区域筛选出与之相关的候选双末端读段，构造候选读段集合R_r；

2)针对每一个缺失区域，采用分段扩展和最长可变拆分读段的组合策略识别缺失断点，得到飘移缺失断点候选集合；

其中，通过分段扩展策略确定R_r中各缺失区域R_ri的读段正确比对位置，排除克隆演变中可能引起的可变缺失断点；最长可变拆分读段采用可容错的动态规划算法寻找最大公共子串，以此确定比对位置与缺失断点；

3)根据断点对数目区分飘移缺失与固定缺失；对于区分出的固定缺失，采用贪婪策略识别缺失周围混乱的单核苷酸变异位点并逐级分离各亚克隆所属的变异位点；

4)获得克隆数目，并依据克隆数目信息优化飘移缺失识别结果；

所述的优化是对识别出的飘移缺失区域进行过滤，若一个区域的飘移缺失个数大于识别出来的子克隆数目k，则只输出前k个飘移缺失。

2.如权利要求1所述的非医疗诊断的肿瘤复杂克隆结构的缺失变异识别及克隆计数方法，其特征在于，所述的与缺失突变相关的所有区域范围的确定为：

1)采用DELLY方法识别与缺失突变相关的所有断点集合，每一对断点被视为一个缺失区域，形成候选缺失区域集合；

2)采用BWA比对工具获得读段比对信息，然后为每个缺失区域筛选出与之相关的候选双末端读段，形成候选读段信息集合；集合中的读段必须满足至少存在一端比对结果的CIGAR标识中包括除“M”以外的字符。

3.如权利要求1或2所述的非医疗诊断的肿瘤复杂克隆结构的缺失变异识别及克隆计数方法，其特征在于，飘移缺失的识别为：

1)以DELLY方法检测的缺失区域为参考，将其对应的参考基因组片段作为核心，构造缺失参考序列；对于第i个缺失，假设缺失的左断点为

右断点为

分别向两侧扩展缺失区域并确定缺失参考序列d_ri，范围区间为

其中，插入间距的均值μ和标准差σ由Picard计算所得，读段长度为l；

令D_r表示缺失参考序列集合，D_r＝{d_r1，…，d_ri，…，d_rn}，其中n为缺失突变数目，d_ri表示第i个缺失参考序列；

2)将所有读段与参考基因组序列进行比对，获得BAM/SAM格式的比对结果文件，再在比对结果文件中进一步搜索获得与每个缺失相关的候选读段；识别每个缺失区域的候选双末端读段必须满足至少有一端比对结果的CIGAR标识中包括除“M”以外的字符；

令R_r表示候选读段集合，R_r＝{R_r1，…，R_ri，…，R_rn}，其中R_ri表示第i个缺失区域对应的候选读段集合；

3)使用分段扩展策略实现R_ri的重比对：

(1)局部缺失参考序列划分成固定k碱基长度的小片段，称为k-mers；相邻的k-mers间满足(k-1)个碱基重叠并将所有的k-mers存储于哈希表中；k-mers作为哈希表中的关键字，并设置k的默认值，每个k-mer在参考基因组中的比对位置作为键内容；

(2)读段划分为无覆盖等长的k-mers，作为哈希表的键；所有读段的k-mers映射于种子哈希表中且返回每个k-mer完美比对的位置列表信息；并通过种子间相邻位置关系过滤掉比对位置结果中的部分假阳性结果；

4)采用最长可变拆分读段方法识别缺失断点：

对于第i个缺失参考序列，如果在R_ri中双末端读段的两端均能映射到d_ri中，且至少一端读段的汉明距离大于2bp，则将具有最大距离的一端读段存储于集合B_d中，集合B_d用于保存满足条件的读段；

采用最大公共子串策略确定比对位置与断点，最大公共子串指读段被拆分的部分与参考序列具有最大的公共子序列；采用具有容错机制的动态规划算法识别最大公共子串；

读段在缺失参考序列中的比对位置首先由分段扩展重比对方法确定，然后分别对读段拆分的头尾两部分进行定位；对于头部分，在读段与缺失参考序列中起始位置为比对位置的等长序列之间寻找最大公共子串，直到容错率不满足为止；对于尾部分采用相同的方法，在缺失参考序列中寻找尾部的起始位置；在容错上，考虑满足头部容错数且头与尾部的总容错数范围为[0，5bp]；

统计每个断点对的读段支持数，并从候选缺失集合中删除支持数小于2的断点对；在统计时，基于克隆数目确定飘移缺失数目的上限并调整过滤规则，过滤掉由错误比对引起的缺失；还将断点对按读段支持数升序存储，当一条读段支持多个断点时，仅将此读段统计到具有最高支持数的断点对中，并相应地修改其它断点对的支持数；

识别所有缺失断点后，得到缺失集合D_s，D_s＝{d_sj}_j＝1:m，

m为缺失数目；

为缺失片段的左断点；

为缺失片段的右断点；

如果m>1，为飘移缺失，否则为固定缺失；若飘移缺失数目超过亚克隆数目n，则以支持数最高的n个断点对作为最终识别结果。

4.如权利要求3所述的非医疗诊断的肿瘤复杂克隆结构的缺失变异识别及克隆计数方法，其特征在于，所述的过滤掉部分假阳性结果的操作为：

将第一个k-mer的位置作为读段的候选位置；对于一个读段，假定第i个k-mer有多个比对位置，此读段也将有多个比对位置；采用所有k-mers投票确定读段的候选位置，每个候选位置的k-mers支持数为投票数，且返回最小汉明距离允许范围内的位置作为最终位置；存储读段候选位置信息包括至少σ_s个完美比对的k-mers，σ_s为与读段长度相关的支持数阈值。

5.如权利要求3所述的非医疗诊断的肿瘤复杂克隆结构的缺失变异识别及克隆计数方法，其特征在于，所述的动态规划的表达式为：

其中，A[i]与B[j]分别表示序列A和B中的第i个和第j个碱基，k表示容错碱基数；DP[k][i][j]表示容错为k且以A[i]和B[j]为结束的最大公共子串长度。

6.如权利要求1或3所述的非医疗诊断的肿瘤复杂克隆结构的缺失变异识别及克隆计数方法，其特征在于，所述的逐级分离策略识别每个亚克隆中的单核苷酸变异位点SNVs，包括以下操作：

1)按照断点对数目区分固定缺失与飘移缺失，再根据固定缺失周围的SNVs数目区分固定缺失与简单缺失；设定SNVs数目的阈值为σ_s；

2)构建每个固定缺失的SNV参考序列，它包括

和

两个子序列，len值的设定必须满足一定的长度，由此确保获得足够的读段信息；

候选读段集合R_r中的读段依次与对应的SNV参考序列重比对；如果双末端读段的两端均与SNV参考序列相匹配，将具有最大汉明距离的一端存储于集合B_s中；否则，当只有一端读段与SNV参考序列相匹配时，将不匹配的另一端存储于集合O_s中；

3)使用B_s集合中的读段来校正D_s中的断点；采用最长可变拆分读段方法确定读段在每个候选SNV参考序列中的位置，最可信的断点对由具有最大读段支持数的候选SNV参考序列报出；

4)在校正断点之后使用集合O_s和B_s中的读段识别SNV；

对于集合O_s首先通过双末端读段的插入距离确定比对位置，插入距离服从均值为μ，标准差为σ的正态分布，插入距离应在[μ-3σ，μ+3σ]范围内；假定读段比对位置为P_s，如果未匹配读段是正向，它的位置范围为[P_s-(μ-l)-3σ，P_s-(μ-l)+3σ]；

否则，未匹配读段是反向，它的位置范围是[P_s+(μ-l)-3σ，P_s+(μ-l)+3σ]；具有最小汉明距离的位置被选作为读段的最终比对位置，之后逐个碱基比对识别SNVs；

5)将固定缺失周围的SNVs分离到不同的亚克隆中，其中SNV的支持数必须大于设置的阈值，将包含混乱SNVs的区域命名为S。

7.如权利要求6所述的非医疗诊断的肿瘤复杂克隆结构的缺失变异识别及克隆计数方法，其特征在于，分离SNVs的方法如下：

1)按比对位置升序存储S中的所有SNVs；

2)将每个读段中的SNVs与S中的SNVs相比较，识别出最后一代亚克隆的读段，如果一条读段包含所有SNVs，则其必定来自于最后一代亚克隆；

3)移除最后一代克隆中的读段并修改S；如果读段移除后，支持数小于阈值的SNV均被视为最后一代亚克隆；

4)使用2)估计移除的读段数目，如果此值小于阈值，将其直接删除；

5)迭代执行2)到4)，直至S为空或执行次数大于假设的最大克隆数为止，实现克隆计数，获得克隆数目。

Images