一种样本级差分索引构建方法及查询方法

本发明属于生物信息学，涉及一种样本级差分索引构建方法及查询方法，尤其涉及一种通过存储样本与参考基因组差异信息实现高效存储和查询的样本级差分索引构建及查询方法。

背景技术：

1、随着高通量测序技术的快速发展，大规模群体基因组学研究已成为生物医学领域的重要方向。在基因组数据分析中，序列比对和查询是基础性操作，其效率直接决定了后续分析的质量和速度。传统的生物序列索引方法，如bwa（burrows-wheeler aligner，li, h., et al. fast and accurate short read alignment with burrows-wheelertransform. bioinformatics. 2009: 25(14), 1754-1760）和bowtie2（langmead, b. ,et al. fast gapped-read alignment with bowtie 2. nature methods, 2012: 9(4),357-359）所采用的fm-index（full-text index in minute space），为每个样本构建完整的基因组索引。然而，该方法存在以下显著缺陷，一方面是存储开销巨大。 对于每个样本，完整基因组索引的存储需求为数gb级别。当处理n个样本时，需要存储n个gb级的完整索引，存储开销呈线性增长。对于万样本级别的群体基因组学项目，总存储需求可达数十tb，不仅造成高昂的存储成本，而且查询时频繁的i/o操作成为严重的性能瓶颈。另一方面是缺乏跨样本的压缩机制。 现有方法未充分利用个体基因组与参考基因组之间的高度相似性。研究表明，人类个体基因组与参考基因组之间的一致性超过99.9%，即任意两个人类个体基因组之间仅存在约数百万个变异位点。传统方法为每个样本独立构建完整索引，未能利用这一冗余信息进行存储优化。

2、近年来，泛基因组学研究取得重要进展（ferro, e., et al. building apangenome alignment index via recursive prefix-free parsing. iscience, 2024:27(4), 109421），通过利用个体基因组与参考基因组的高度相似性，可采用差异编码的思想降低存储需求。现有差异存储格式如vcf（variant call format）虽然能够记录样本与参考基因组之间的变异信息，但其设计目标主要是变异数据的存储和交换，存在以下不足：（1）不支持从参考坐标到样本坐标的快速转换；（2）不支持样本序列的即时重建；（3）查询效率低，无法满足实时分析需求。此外，以vg（variation graph）和giraffe为代表的变异图方法虽然通过构建包含多个基因组信息的图结构实现了多基因组索引，但其构建复杂度高、内存需求大（ sirén, j., et al. haplotype-aware graph indexes. bioinformatics,2019: 35(9), 1577-1579），在大规模群体中的应用受到限制。

3、因此，本领域亟需一种既能大幅压缩存储空间，又能支持快速坐标转换和序列即时重建的样本级索引方法。

技术实现思路

1、针对现有技术中为每个样本构建完整基因组索引所导致的存储开销巨大、无法利用个体基因组与参考基因组相似性进行压缩、以及现有差异存储格式不支持快速坐标转换和序列即时重建的技术问题，本发明的目的在于提供一种样本级差分索引构建方法及查询方法。

2、为实现上述目的，本发明的技术方案为：

3、一种样本级差分索引构建方法，其步骤包括：

4、将样本基因组测序数据与参考基因组进行比对，检测样本基因组测序数据与参考基因组之间的所有变异，生成有序变异列表v；

5、基于所述有序变异列表v，计算参考基因组中每个参考位置相对于该参考位置在样本基因组测序数据中的样本坐标的累积偏移量，构建坐标映射表m；

6、提取样本基因组测序数据中满足设定条件的完整序列存入差分序列缓存c；

7、基于有序变异列表v、坐标映射表m、差分序列缓存c生成样本级差分索引。

8、优选的，样本级差分索引为d(s) = (v, m, c, p)；其中， p为样本级差分索引的元数据。

9、优选的，生成有序变异列表v的方法为：

10、采用alt-aware比对模式，将样本基因组测序数据的测序读段比对到参考基因组，生成比对结果；

11、基于所述比对结果进行变异检测，生成原始变异集合；

12、采用质量值和测序深度对所述原始变异集合进行质量过滤；

13、对过滤后的变异按参考坐标进行升序排序；

14、对相邻的复杂变异进行合并处理，生成有序变异列表v。

15、优选的，构建坐标映射表m的方法为：

16、设参考基因组任一碱基位置为整数坐标，该整数坐标对应的样本坐标为该碱基在样本基因组测序数据中的位置；定义偏移量函数o(p)为参考位置p处的累积偏移量，初始值o(0)=0；

17、对于按参考坐标升序排列的变异记录，依次处理每条变异记录，根据变异类型计算偏移量变化量δ = len(alt) - len(ref)；其中，len(alt)为变异记录中替代等位基因的序列长度； len(ref)为变异记录中参考等位基因的序列长度；

18、计算参考基因组中参考位置

p的累积偏移量；其中，表示变异记录的参考坐标，表示变异记录中参考等位基因的序列长度，表示变异记录中替代等位基因的序列长度；

19、采用分层索引结构构建坐标映射表m，用于存储锚点的累积偏移量及粗粒度区间内的变异增量；其中，以固定间隔l在参考基因组上设置锚点生成多个粗粒度区间，每个位置上的锚点记录该位置的参考坐标及其对应的累积偏移量，根据各锚点构建粗粒度锚点索引；在每个粗粒度区间内构建细粒度索引，存储该粗粒度区间内所有变异记录的完整信息。

20、优选的，所述设定条件包括：（a）包含结构变异的区域；（b）与参考基因组序列一致性低于设定阈值的区域；（c）样本特异性插入序列；所述差分序列缓存c中的每个条目包含区域起始坐标、区域终止坐标、序列类型标识、序列数据。

21、一种基于所构建样本级差分索引的查询方法，其特征在于，

22、接收查询请求，获取目标参考坐标x；

23、在坐标映射表m的粗粒度锚点索引中定位目标参考坐标x所在的粗粒度区间[l,r]，获取该粗粒度区间[l, r]起始锚点处的基准偏移量o_base；

24、在该粗粒度区间[l, r]对应的细粒度索引中，在区间[l, x)内通过二分查找定位目标参考坐标x之前的最后一个变异记录v_last；

25、从基准偏移量o_base出发，累加从该粗粒度区间[l, r]起始锚点l处的锚点到v_last之间所有变异导致的偏移量变化，得到精确累积偏移量o_accurate；

26、计算样本坐标y = x + o_accurate；

27、若y为空则返回无对应位置，否则返回样本坐标y。

28、一种基于所构建样本级差分索引的样本序列重建方法，其步骤包括：

29、根据待重建样本序列获取参考基因组中目标区域[region_start, region_end]的参考序列s_ref；

30、从变异列表v中获取位于该目标区域[region_start, region_end]内的所有变异记录，按参考坐标升序排列；

31、初始化当前参考位置cur = region_start，样本序列构建列表seq_list = []；

32、依次处理该目标区域[region_start, region_end]内每条变异记录v，根据变异类型分别应用对应的序列编辑操作；

33、处理完所有变异后，若cur ≤ region_end，则将参考序列s_ref中剩余部分添加到seq_list；

34、将seq_list中的所有片段拼接，输出完整的样本序列。

35、一种样本级差分索引构建系统，其特征在于，包括

36、有序变异列表生成模块，用于将样本基因组测序数据与参考基因组进行比对，检测样本基因组测序数据与参考基因组之间的所有变异，生成有序变异列表v；

37、坐标映射表构建模块，用于基于所述有序变异列表v，计算参考基因组中每个参考位置相对于该参考位置在样本基因组测序数据中的样本坐标的累积偏移量，构建坐标映射表m；

38、差分序列缓存模块，用于提取样本基因组测序数据中满足设定条件的完整序列存入差分序列缓存c；

39、索引生成模块，用于基于有序变异列表v、坐标映射表m、差分序列缓存c生成样本级差分索引。

40、一种计算设备，其特征在于，包括：处理器、存储有计算机程序的存储器，所述计算机程序被处理器运行时，执行上述的方法。

41、一种计算机可读存储介质，其特征在于，存储指令，当所述指令在计算机上运行时，使得计算机执行上述的方法。

42、本发明的优点如下：

43、1. 存储效率显著提升。 本发明通过仅存储样本与参考基因组的差异信息，将存储复杂度从o(n × gb)降低到o(gb + n × mb)。对于万样本级别的群体基因组学项目，存储需求可从数十tb降低至数百gb，压缩比达90%以上。

44、2. 坐标转换速度快。 本发明采用分层索引结构，坐标转换的时间复杂度为o(logn)，平均转换时间小于1毫秒，支持实时查询需求。

45、3. 序列重建灵活高效。 本发明支持从参考序列出发，通过应用变异列表即时重建任意区域样本序列，无需存储完整的样本基因组序列。

46、4. 增量更新友好。 新增样本时只需构建该样本的差分索引，无需重建全局索引，具有良好的可扩展性。

47、5. 兼容现有生物信息学生态系统。 本发明可与bwa、gatk、minimap2等主流工具链无缝集成，降低用户的学习和迁移成本。