基于长读长测序的目标区域串联重复变异分型方法和装置与流程

本发明涉及基因测序，具体而言，涉及一种基于长读长测序的目标区域串联重复变异分型方法和装置。

背景技术：

1、串联重复序列（tandom repeat，tr）是指基因组中一段短的dna序列在特定区域内连续重复排列的现象。根据重复单元的长度，串联重复一般可以分为两类：短串联重复（short tandem repeats, str），重复单元一般为2-6bp和可变数目串联重复（variablenumber tandem repeats, vntr），重复单元一般为7-100bp。这些重复序列在基因组的调控、结构、以及进化中起着重要作用，同时也因其高变异性被广泛应用于遗传学研究（如亨廷顿舞蹈症、脆性x综合征等）、法医学个体识别和种群遗传学分析等领域。

2、传统二代测序（ngs）（如illumina平台），因读长有限（通常<300bp），往往难以覆盖较长或复杂结构的tr区域，导致在tr区域的覆盖度不均、比对错误率高，进而造成拷贝数估计偏差和分型错误。随着长读长测序技术（第三代测序，tgs）的发展，特别是pacbio singlemolecule real-time (smrt)测序和oxford nanopore technologies (ont)测序平台的普及，单次测序读长可达数十甚至上百kb，能够直接跨越整个tr区域，从而在tr检测的精度和分型一致性方面展现出显著优势。

3、目前，基于tgs数据的tr检测软件已有一定发展，但各有侧重和局限：

4、straglr：主要采用无参模式，旨在检测基因组中新发或显著扩增的tr，不依赖参考tr库，因此难以针对已知的特定致病位点进行精细的定量分析和分型。

5、tricolor：同样采用无参模式，适用于全基因组范围内的新tr发现，但在针对特定已知位点的定量精度上有限，且输出结果较为宏观，缺乏单分子层面的详细解析。

6、trgt：基于已知tr区域进行检测，支持聚类分型和基于长度的分型，并能分析甲基化状态。然而，trgt主要提供每个样本的总体分型结果，无法直接输出每条测序reads的单倍型归属信息。当样本中两个单倍型在长度或序列上高度相似时，trgt的分型准确性会显著下降。

7、其他工具：如 repeathmm-long、tredparse、tandem-genotypes 等，分别在特定位点检测或全基因组扫描上有一定优势，但在长重复、复杂重复或极高测序深度下仍存在误分型、reads 分配不准确等问题。

8、尽管第三代测序技术显著改善了tr检测能力，但在针对高深度靶向pcr产物测序数据的应用场景中，现有技术仍存在明显的空白和缺陷，具体表现为：

9、1）缺乏逐条reads的精准结构解析能力：现有的主流软件大多假定覆盖深度适中，以样本整体为单位进行分型统计，无法直接对pcr富集后产生的高深度靶向测序数据进行逐条reads（read-by-read）的拷贝数（copy number, cpn）和序列结构解析。在高深度测序数据中，每一张reads都携带了独立的分子信息，忽略单分子层面的解析会导致对杂合子中稀有等位基因或复杂结构变异的漏检。此外，对于疾病样本，研究重复序列之间的打断（例如，ctg重复中出现ga的打断等，如ctgctgctggactg）也具有重要意义。有研究表明，不同疾病对应的打断特征是不同的。

10、2）单倍型分型（haplotyping）能力不足：在二倍体生物中，区分来自父本和母本的两个等位基因（即单倍型）至关重要。现有工具如trgt虽然提供分型结果，但往往无法将单倍型信息回溯到具体的测序reads上。当两个单倍型的tr长度相近或序列相似度较高时，现有算法难以准确区分，导致分型错误。而基于单分子水平的snp/indel辅助分型工具（如hiphase等）通常需结合专门的变异检测流程，缺乏针对tr区域的集成化解决方案。

11、3）高深度数据下的分型矫正机制缺失：在pcr产物的高深度测序中，由于扩增偏好性或测序误差，可能导致某些单倍型的reads比例失衡或出现未分配（unassigned）reads。现有技术缺乏有效的统计矫正算法来处理这些边缘情况（如纯合子误判、单倍型比例严重失衡等），限制了其在临床诊断中的准确性。

12、因此，亟需开发一种能够基于第三代测序数据，针对高深度测序结果（例如，pcr产物高深度测序结果），在单条reads层面精确识别串联重复拷贝数，结合多源变异信息进行单倍型分型，并具备智能矫正机制的检测方法及系统。

技术实现思路

1、本发明的主要目的在于提供一种基于长读长测序的目标区域串联重复变异分型方法，以解决现有技术中靶向高深度测序数据分型不准、缺乏单分子结构可视化及无法精准归属单倍型的问题。

2、为了实现上述目的，根据本发明的第一个方面，提供了一种基于长读长测序的目标区域串联重复变异分型方法，该目标区域串联重复变异分型方法包括：

3、s1，获取目标区域的长读长测序数据，将上述长读长测序数据比对至参考基因组，生成比对文件，筛选出位于上述目标区域内的测序reads，得到目标测序reads；

4、s2，对每条上述目标测序reads进行序列解析，识别上述目标区域内串联重复单元的序列，得到上述串联重复单元的位置和拷贝数；

5、s3，基于上述比对文件，检测上述目标区域内的多源变异信息，得到多源变异数据；

6、s4，利用单倍型构建算法，结合上述拷贝数和上述多源变异数据，对每条上述目标测序reads进行单倍型归属判定，并根据归属将每条上述目标测序reads归类为第一单倍型、第二单倍型或未分配类别，得到单倍型分型初步结果；

7、s5，基于上述第一单倍型、第二单倍型和上述未分配类别的目标测序reads比例分布和/或拷贝数分布特征，对上述未分配类别中的目标测序reads进行重新分配，得到单倍型分型最终结果；

8、其中，上述多源变异信息包括如下任意多种：单核苷酸变异、小片段插入缺失和结构变异；上述串联重复单元为一种或多种。

9、进一步地，上述s2包括如下步骤：

10、s21，构建上述串联重复单元的正则表达式或序列模式匹配规则；

11、s22，利用上述正则表达式或序列模式匹配规则对上述目标测序reads进行扫描，定位所有匹配到的上述串联重复单元的位置；

12、s23，统计每条上述目标测序reads中匹配到的串联重复单元的数量，得到上述拷贝数。

13、进一步地，上述s1中，上述长读长测序数据为pcr产物的测序数据时，在筛选出位于上述目标区域内的上述测序reads之后，得到上述目标测序reads之前，还需要基于引物序列匹配度对上述测序reads进行有效性校验，保留测通的测序reads，从上述测通的测序reads中获得上述目标测序reads。

14、进一步地，从上述测通的测序reads中获得上述目标测序reads的步骤包括：去除上述测通的测序reads两端的引物序列和侧翼序列，得到截取后reads，上述截取后reads即为上述目标测序reads。

15、进一步地，上述s3包括如下步骤：

16、s31，使用单核苷酸变异检测工具对上述比对文件进行分析，提取上述目标区域内的单核苷酸变异位点信息和小片段插入缺失位点信息；

17、s32，使用结构变异检测工具对上述比对文件进行分析，提取上述目标区域内的结构变异位点信息；

18、s33，整合上述单核苷酸变异位点信息、小片段插入缺失位点信息和上述结构变异位点信息，得到上述多源变异数据。

19、进一步地，上述单核苷酸变异检测工具为deepvariant或longshot，上述结构变异检测工具为sawfish、cutesv或clair3。

20、进一步地，上述s4包括如下步骤：

21、s41，将上述拷贝数和上述多源变异数据进行合并，生成综合变异vcf文件；

22、s42，将上述综合变异vcf文件输入至上述单倍型构建算法中，对每条上述目标测序reads进行单倍型相位推断，输出每条目标测序reads所属的单倍型类别。

23、进一步地，上述s5包括如下步骤：

24、s51，统计上述第一单倍型、上述第二单倍型和上述未分配类别中上述目标测序reads的数量，计算各自占总有效测序reads的比例；

25、s52，判断上述第一单倍型的比例和上述第二单倍型的比例是否均大于或等于第一预设阈值；若是，则判定分型结果稳定，不进行重新分配，将当前分型结果作为上述单倍型分型最终结果；若否，执行s53；

26、s53，判断上述第一单倍型或上述第二单倍型的比例大于或等于第二预设阈值，且上述未分配类别的比例大于或等于第三预设阈值；若是，则将上述未分配类别重新校正为上述第一单倍型或第二单倍型中比例较小的单倍型类别，并将上述比例较小的单倍型类别重新矫正为未分配类别，得到新的分型结果，将上述新的分型结果记为上述单倍型分型最终结果；若否，执行s54；

27、s54，获取所有的上述目标测序reads中串联重复单元的拷贝数的分布直方图，识别出前两个最高峰值对应的拷贝数值区间，以这两个拷贝数值区间的中间点作为拷贝数阈值；

28、s55，根据上述拷贝数阈值，将所有的上述目标测序reads区分为两个单倍型；计算上述两个单倍型的比例，判断上述两个单倍型的比例是否为（4~10）:1，若是，则上述两个单倍型为上述单倍型分型最终结果，若否，执行s56；

29、s56，将上述目标测序reads随机分配为比例为1:1的两个单倍型，则上述两个单倍型为上述单倍型分型最终结果；

30、其中，上述第一预设阈值、第二预设阈值和第三预设阈值均为30%~35%。

31、进一步地，在得到上述单倍型分型最终结果后，上述目标区域串联重复变异分型方法进一步包括：

32、s6，根据上述单倍型分型最终结果，针对每个单倍型，基于其包含的上述目标测序reads的串联重复单元序列生成共识序列，计算每个单倍型的拷贝数统计指标和序列纯度，并生成可视化结果。

33、进一步地，上述s6包括如下步骤：

34、s61，提取每个单倍型的所有上述目标测序reads的上述串联重复单元的序列，生成上述共识序列，将每条上述目标测序reads的上述串联重复单元的序列与上述共识序列进行比对，统计完全匹配的比例，得到上述单倍型的上述序列纯度；

35、s62，计算每个单倍型的所有上述目标测序reads的上述拷贝数的众数、中位数和平均数；

36、s63，筛选出拷贝数位于上述众数附近的预设数量条上述目标测序reads，展示上述目标测序reads中串联重复单元的序列结构和变异信息。

37、为了实现上述目的，根据本发明的第二个方面，提供了一种基于长读长测序的目标区域串联重复变异分型装置，该装置包括：

38、比对筛选模块，被设置为获取目标区域的长读长测序数据，将上述长读长测序数据比对至参考基因组，生成比对文件，筛选出位于上述目标区域内的测序reads，得到目标测序reads；

39、解析模块，被设置为对每条上述目标测序reads进行序列解析，识别上述目标区域内串联重复单元的序列，得到上述串联重复单元的位置和拷贝数；

40、检测模块，被设置为基于上述比对文件，检测上述目标区域内的多源变异信息，得到多源变异数据；

41、单倍型分型第一模块，被设置为利用单倍型构建算法，结合上述拷贝数和上述多源变异数据，对每条上述目标测序reads进行单倍型归属判定，并根据归属将每条上述目标测序reads归类为第一单倍型、第二单倍型或未分配类别，得到单倍型分型初步结果；

42、单倍型分型第二模块，被设置为基于上述第一单倍型、第二单倍型和上述未分配类别的目标测序reads比例分布和/或拷贝数分布特征，对上述未分配类别中的目标测序reads进行重新分配，得到单倍型分型最终结果；

43、其中，上述多源变异信息包括如下任意多种：单核苷酸变异、小片段插入缺失和结构变异；上述串联重复单元为一种或多种。

44、进一步地，上述解析模块包括：

45、构建单元，被设置为构建上述串联重复单元的正则表达式或序列模式匹配规则；

46、定位单元，被设置为利用上述正则表达式或序列模式匹配规则对上述目标测序reads进行扫描，定位所有匹配到的上述串联重复单元的位置；

47、拷贝数统计单元，被设置为统计每条上述目标测序reads中匹配到的串联重复单元的数量，得到上述拷贝数。

48、进一步地，上述长读长测序数据为pcr产物的测序数据时，上述比对筛选模块还包括：校验单元，被设置为在筛选出位于上述目标区域内的上述测序reads之后，得到上述目标测序reads之前，还需要基于引物序列匹配度对上述测序reads进行有效性校验，保留测通的测序reads，从上述测通的测序reads中获得上述目标测序reads。

49、进一步地，上述校验单元还包括：截取模块，被设置为去除上述测通的测序reads两端的引物序列和侧翼序列，得到截取后reads，上述截取后reads即为上述目标测序reads。

50、进一步地，上述检测模块包括：第一分析单元，被设置为使用单核苷酸变异检测工具对上述比对文件进行分析，提取上述目标区域内的单核苷酸变异位点信息和小片段插入缺失位点信息；第二分析单元，被设置为使用结构变异检测工具对上述比对文件进行分析，提取上述目标区域内的结构变异位点信息；整合单元，被设置为整合上述单核苷酸变异位点信息、小片段插入缺失位点信息和上述结构变异位点信息，得到上述多源变异数据。

51、进一步地，上述单核苷酸变异检测工具为deepvariant或longshot，上述结构变异检测工具为sawfish、cutesv或clair3。

52、进一步地，上述单倍型分型第一模块包括：合并单元，被设置为将上述拷贝数和上述多源变异数据进行合并，生成综合变异vcf文件；单倍型相位推断单元，将上述综合变异vcf文件输入至上述单倍型构建算法中，对每条上述目标测序reads进行单倍型相位推断，输出每条目标测序reads所属的单倍型类别。

53、进一步地，上述单倍型分型第二模块包括：第一统计单元，被设置为统计上述第一单倍型、上述第二单倍型和上述未分配类别中上述目标测序reads的数量，计算各自占总有效测序reads的比例；第一分配单元，被设置为判断第一单倍型的比例和第二单倍型的比例是否均大于或等于第一预设阈值；若是，则判定分型结果稳定，不进行重新分配，将当前分型结果作为单倍型分型最终结果；若否，执行第二分配单元；上述第二分配单元，被设置为判断上述第一单倍型或上述第二单倍型的比例大于或等于第二预设阈值，且上述未分配类别的比例大于或等于第三预设阈值；若是，则将未分配类别重新校正为上述第一单倍型或第二单倍型中比例较小的单倍型类别，并将比例较小的单倍型类别重新矫正为未分配类别，得到新的分型结果后，将上述新的分型结果记为单倍型分型最终结果；若否，依次执行获取单元和第三分配单元；上述获取单元，被设置为获取所有的上述目标测序reads中串联重复单元的拷贝数的分布直方图，识别出前两个最高峰值对应的拷贝数值区间，以这两个拷贝数值区间的中间点作为拷贝数阈值；上述第三分配单元，被设置为根据上述拷贝数阈值，将所有的上述目标测序reads区分为两个单倍型；计算上述两个单倍型的比例，判断上述两个单倍型的比例是否为（4~10）:1，若是，则上述两个单倍型为单倍型分型最终结果，若否，执行第四分配单元；上述第四分配单元，被设置为将上述目标测序reads随机分配为比例为1:1的两个单倍型，则上述两个单倍型为单倍型最终结果；其中，上述第一预设阈值、第二预设阈值和第三预设阈值均为30%~35%。

54、进一步地，上述装置还包括：输出模块，被设置为根据上述单倍型分型最终结果，针对每个单倍型，基于其包含的上述目标测序reads的串联重复单元序列生成共识序列，计算每个单倍型的拷贝数统计指标和序列纯度，并生成可视化结果。

55、进一步地，上述输出模块包括：

56、比对单元，被设置为提取每个单倍型的所有上述目标测序reads的上述串联重复单元的序列，生成上述共识序列，将每条上述目标测序reads的上述串联重复单元的序列与上述共识序列进行比对，统计完全匹配的比例，得到上述单倍型的上述序列纯度；

57、计算单元，被设置为计算每个单倍型的所有上述目标测序reads的上述拷贝数的众数、中位数和平均数；

58、展示单元，被设置为筛选出拷贝数位于上述众数附近的预设数量条上述目标测序reads，展示上述目标测序reads中串联重复单元的序列结构和变异信息。

59、为了实现上述目的，根据本发明的第三个方面，提供了一种计算机可读存储介质，上述存储介质包括存储的程序，其中，在上述程序运行时控制上述存储介质所在设备执行上述的基于长读长测序的目标区域串联重复变异分型方法。

60、为了实现上述目的，根据本发明的第四个方面，提供了一种处理器，上述处理器用于运行程序，其中，上述程序运行时执行上述的基于长读长测序的目标区域串联重复变异分型方法。

61、应用本发明的技术方案，通过获取目标区域的长读长测序数据并比对，筛选得到目标测序reads；对每条目标测序reads进行序列解析，识别串联重复单元的位置和拷贝数；基于比对文件检测多源变异信息（包括单核苷酸变异、小片段插入缺失和结构变异中的至少两种）；利用单倍型构建算法，结合拷贝数和多源变异数据，对每条目标测序reads进行单倍型归属判定，将其归类为第一单倍型、第二单倍型或未分配类别；最后基于两类单倍型及未分配类别的reads比例分布和/或拷贝数分布特征，对未分配类别中的目标测序reads进行重新分配，得到单倍型分型最终结果。

62、上述方法通过逐条reads的cpn检测与单倍型矫正，克服了现有软件在高深度测序中无法针对已知motif进行单分子层面定量及分型准确性的瓶颈，显著提升了复杂串联重复区域（如长重复或高变异区域）的分型精度；此外，该方法能够有效区分杂合子与纯合子样本，并通过计算motif纯度及基于单倍型的序列结构可视化，为遗传疾病研究、法医学个体识别及群体遗传分析提供了高分辨率、高一致性的分子分型依据，解决了传统二代测序读长限制导致的检测偏差及现有三代分析工具在靶向高深度数据中reads分配不准确的问题。