一种基于发生甲基化C-T转换特征的全基因组甲基化位点识别方法及系统与流程

本技术属于生物信息学和表观基因组学分析，具体涉及一种基于发生甲基化c-t转换特征的全基因组甲基化位点识别方法及系统。

背景技术：

1、dna甲基化是常见的表观遗传修饰。5-甲基胞嘧啶（5mc）与基因表达调控、基因组印记、x染色体失活以及多种疾病过程有关。现有全基因组甲基化检测技术主要包括亚硫酸氢盐测序（wgbs）、tet辅助吡啶硼烷测序（taps）以及illumina 5-base dna prep等。

2、在taps和illumina 5-base技术体系中，甲基化胞嘧啶经化学处理后，在测序读段中表现为胸腺嘧啶（t）；未甲基化胞嘧啶仍保持为胞嘧啶（c）。因此，在参考基因组胞嘧啶位点上，c/t观测比例可以用于估计甲基化水平。对于负链胞嘧啶，其在参考基因组上对应为g位点，甲基化转换信号则表现为g/a关系。

3、传统亚硫酸氢盐测序的转换方向不同。该技术主要使未甲基化胞嘧啶转化为尿嘧啶，并在测序中表现为t，而甲基化胞嘧啶相对保留为c。为了适配这种转换，常规wgbs分析通常需要将读段或参考基因组进行c/t转换，构建三碱基或双链转换参考，再使用专门的甲基化比对工具分析。该路线对wgbs是合适的，但也会带来一些实际问题，例如参考构建复杂、碱基复杂度降低、重复匹配增多、唯一比对能力下降，以及与标准bam/vcf流程衔接不够直接。

4、taps和5-base数据不宜简单套用wgbs三碱基流程。原因在于，taps/5-base中未甲基化c仍为c，而发生甲基化的c表现为t。若仍按亚硫酸氢盐测序的分析思路处理，既增加流程复杂度，也容易使甲基化信号解释与实验转换方向不一致。因此，有必要建立一种面向taps/5-base转换特征的分析方法，直接利用标准参考基因组、标准bam文件和碱基堆积结果完成甲基化识别。

5、现有方法在该场景下主要存在以下不足。

6、1. 位点覆盖不完整。部分方法以变异检测结果作为甲基化信号来源，只能得到被变异检测器输出的c>t或g>a候选位点。低甲基化比例位点、弱转换信号位点，以及未达到变异候选阈值的胞嘧啶位点容易被遗漏。

7、2. 甲基化信号和变异信号容易混在一起。变异检测工具的目标是识别基因型、体细胞突变或等位基因状态，并不负责区分化学转换产生的c>t/g>a信号和真实遗传变异。高甲基化位点可能被误认为c>t或g>a变异，真实变异也可能被计入甲基化比例。

8、3. 甲基化比例不能直接由变异模型替代。甲基化水平本质上是同一胞嘧啶位点中甲基化分子和未甲基化分子的比例。常规变异检测模型输出的是基因型、等位基因频率或候选变异状态，不能稳定地给出每个胞嘧啶位点的甲基化读段数、未甲基化读段数和甲基化比例。

9、4. 正负链解释需要统一。taps/5-base数据中，正链参考c位点表现为c/t信号；负链胞嘧啶在参考基因组上表现为g/a信号。若流程只处理参考c位点，或者没有明确负链g位点的判读规则，会影响双链胞嘧啶位点的完整识别。

10、5. 上下文信息往往需要后处理。甲基化分析通常需要区分cg、chg和chh类型，并保留三碱基序列。通用碱基堆积工具只提供原始位点观测信息，缺少面向甲基化c-t转换特征的正负链判读、上下文分类、三碱基信息提取和统一结果输出。

11、6. 输出格式不统一。不同工具的结果字段差异较大，不利于后续质控、甲基化水平统计、可视化和甲基化差异分析。

12、7. 三碱基比对与常规基因组分析体系衔接不够顺畅。三碱基流程通常需要专用索引、专用比对器或转换参考。若同一样本还要做snv、indel、cnv、sv或bam质控，往往需要额外维护一套标准参考坐标下的结果。

13、8. 亚硫酸氢盐测序本身可能带来dna损伤、覆盖不均一和序列复杂度下降等问题。对于taps/5-base数据，继续沿用亚硫酸氢盐流程并不合适。

技术实现思路

1、目的：鉴于以上技术问题中的至少一项，本技术提供一种基于发生甲基化c-t转换特征的全基因组甲基化位点识别方法及系统，从bam碱基堆积出发进行分析，对全基因组c/g位点连续扫描，按taps/5-base转换规则解释c>t和g>a信号，并同时输出甲基化信号计数、未甲基化信号计数、上下文分类和三碱基序列。

2、本技术采用的技术方案为：

3、第一方面，本技术提供了一种基于发生甲基化c-t转换特征的全基因组甲基化位点识别系统，包括：

4、序列比对模块，用于：将甲基化测序读段比对至原始参考基因组，得到染色体编号和基因组位置，并生成排序去重去低质量后的bam文件；

5、碱基堆积信息获取模块，用于：对bam文件进行全基因组碱基堆积，得到各位点的碱基堆积字符串和参考碱基；

6、碱基堆积字符串解析模块，用于：根据参考碱基解析碱基堆积字符串，统计各位点a、c、g、t各类型碱基的读段支持数，得到碱基计数结果；

7、甲基化信号提取模块，用于：根据参考碱基类型，按taps或5-base甲基化c-t转换规则从碱基计数结果中提取甲基化信号和未甲基化信号，得到链方向、甲基化信号计数和未甲基化信号计数；

8、变异混杂位点识别与处置模块，用于：基于链方向，根据参考碱基类型和外部小变异vcf文件或已知变异数据库识别c>t或g>a变异混杂位点，并对变异混杂位点进行标记或过滤；

9、甲基化上下文分类及三碱基序列提取模块，用于：根据链方向和参考基因组判断甲基化上下文类型并提取三碱基序列。

10、在一些实施例中，所述的系统还包括结果输出模块，用于：输出甲基化位点识别结果；其中，所述甲基化位点识别结果包括：染色体编号、基因组位置、链方向、甲基化信号计数、未甲基化信号计数、甲基化上下文类型和三碱基序列。

11、在一些实施例中，所述序列比对模块使用通用序列比对工具进行比对，直接使用未经化学修饰转换的原始参考基因组文件，不使用三字母或通配符基因组。

12、在一些实施例中，所述碱基堆积字符串解析模块按照以下碱基堆积字符串编码规则解析：

13、字符“.“和“,“计为与参考碱基相同的碱基类型；

14、大写a/c/g/t和小写a/c/g/t计为对应错配碱基类型；

15、字符“^“及其后一个映射质量字符跳过不计；

16、字符“$“跳过不计；

17、字符“+“或“-“及其后的插入或缺失描述整体跳过不计；

18、字符“*“跳过不计。

19、在一些实施例中，所述甲基化信号提取模块的提取规则为：

20、当参考碱基为c时，读段中t碱基计数结果为甲基化信号计数，c碱基计数结果为未甲基化信号计数，该位点标记为正链；

21、当参考碱基为g时，读段中a碱基计数结果为甲基化信号计数，g碱基计数结果为未甲基化信号计数，该位点标记为负链。

22、在一些实施例中，所述变异混杂位点识别与处置模块具体用于：

23、接收外部小变异vcf文件或已知变异数据库，作为注释层；

24、对参考碱基为c且变异等位基因为t的位点，标记为正链c>t变异混杂位点；对参考碱基为g且变异等位基因为a的位点，标记为负链g>a变异混杂位点；

25、对于变异混杂位点，根据预设模式保留标记或过滤去除。

26、进一步地，对于变异混杂位点，根据预设模式保留标记或过滤去除，包括：

27、若预设模式为标记模式，保留变异混杂位点的甲基化信号计数和甲基化上下文类型，同时输出变异混杂标记、变异等位基因和过滤状态；

28、若预设模式为过滤模式，在生成高置信甲基化位点识别结果时去除变异混杂位点，减少真实变异对甲基化比例的影响。

29、在一些实施例中，所述甲基化上下文分类及三碱基序列提取模块具体用于：

30、对于正链参考c位点，读取参考基因组中该位点下游第1位碱基和第2位碱基，三碱基序列为"c"加参考基因组中该位点下游第1位碱基和第2位碱基；

31、对于负链参考g位点，读取参考基因组中该位点上游第1位碱基和第2位碱基并取其反向互补碱基，三碱基序列为"c"加参考基因组中该位点下游第1位碱基和第2位碱基的反向互补碱基；

32、若该位点下游第1位碱基为g，分类为cg类型；

33、若该位点下游第1位碱基非g且下游第2位碱基为g，分类为chg类型；

34、若该位点下游第1位碱基和第2位碱基均非g，分类为chh类型。

35、第二方面，本技术提供了一种基于发生甲基化c-t转换特征的全基因组甲基化位点识别方法，基于第一方面所述的系统，包括：

36、将甲基化测序读段比对至原始参考基因组，得到染色体编号和基因组位置，并生成排序去重去低质量后的bam文件；

37、对bam文件进行全基因组碱基堆积，得到各位点的碱基堆积字符串和参考碱基；

38、根据参考碱基解析碱基堆积字符串，统计各位点a、c、g、t各类型碱基的读段支持数，得到碱基计数结果；

39、根据参考碱基类型，按taps或5-base甲基化c-t转换规则从碱基计数结果中提取甲基化信号和未甲基化信号，得到链方向、甲基化信号计数和未甲基化信号计数；

40、基于链方向，根据参考碱基类型和外部小变异vcf文件或已知变异数据库识别c>t或g>a变异混杂位点，并对变异混杂位点进行标记或过滤；

41、根据链方向和参考基因组判断甲基化上下文类型并提取三碱基序列。

42、第三方面，本技术提供了一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现第一方面所述方法。

43、有益效果：本技术提供的基于发生甲基化c-t转换特征的全基因组甲基化位点识别方法及系统，不依赖变异检测结果发现甲基化位点，适用于taps和illumina 5-base dnaprep等甲基化测序数据分析。具有以下优点：

44、（1）位点覆盖更完整。本技术不以变异检测结果作为甲基化位点来源，而是基于参考基因组c/g位点和bam碱基堆积逐位点扫描。不考虑位点深度和质量要求，即可输出甲基化信号计数，低比例甲基化位点不会因为未进入vcf而丢失。

45、（2）正负链判读更清楚。本技术明确规定参考c位点以t为甲基化信号、c为未甲基化信号；参考g位点以a为甲基化信号、g为未甲基化信号。该规则可同时覆盖双链胞嘧啶位点。

46、（3）更符合甲基化比例的性质。本技术直接输出甲基化读段数和未甲基化读段数，并据此计算甲基化比例，不依赖基因型或突变置信度等变异检测输出。

47、（4）变异混杂可控。本技术利用外部vcf或已知变异数据库对c>t/g>a混杂位点进行标记或过滤，减少真实遗传变异被误当作甲基化转换信号的风险。

48、（5）上下文信息一次生成。本技术在输出甲基化信号计数的同时给出cg、chg、chh类型和三碱基序列，便于后续cpg筛选、非cpg分析、差异甲基化分析和功能注释。

49、（6）输出格式更利于后续分析。本技术把染色体、位置、链方向、计数、比例、上下文、三碱基序列和混杂状态放在同一结果表中，便于后续甲基化水平和差异甲基化分析。

50、（7）不需要三碱基转换参考。本技术减少了转换参考构建、参数适配和结果回译步骤，部署成本较低。

51、（8）与常规基因组分析更容易衔接。本技术使用标准四碱基参考和标准bam文件，可与常规变异检测、覆盖度分析和注释流程共享同一坐标体系。

52、（9）更适合taps/5-base数据。本技术按甲基化c表现为t、未甲基化c保持为c的实验特征进行解释，避免把亚硫酸氢盐测序的转换模型直接套用于taps/5-base数据。

53、（10）适合大规模全基因组处理。本技术可以逐行读取pileup结果，不需要一次性加载全基因组所有位点，内存压力较小。