CN113889184B - 一种融合基因组特征的m6A甲基化局部功能谱分解方法 - Google Patents

Abstract

本发明名为一种融合基因组特征的m⁶A甲基化局部功能谱分解方法，属表观转录组学及模式识别研究领域。本发明目标是融合基因组特征对RNA甲基化局部功能谱进行高精度的分解，揭示m⁶A甲基化的详细调控机制。基于矩阵分解中的独立成分分析(ICA)方法，对m⁶A甲基化位点对应的基因组特征进行融合，以便充分考虑极强相关m⁶A位点之间的关联性。最终利用基于负熵的估计方法对m⁶A位点参与的多个调控通路进行恢复，进一步挖掘出每个调控通路对应的敏感位点，实现对m⁶A甲基化局部功能的高精度谱分解。本发明提出的算法能够为湿地实验揭示m⁶A甲基化的调控机制提供有力的参考。

Description

一种融合基因组特征的m6A甲基化局部功能谱分解方法

技术领域

本发明涉及表观转录组学及模式识别研究领域，具体涉及一种融合基因组特征的m⁶A甲基化局部功能谱分解方法。

背景技术

在表观转录组学中，N6-甲基腺嘌呤(m⁶A)作为真核生物mRNA中最普遍和最丰富的转录后RNA修饰，是指发生在碱基A第六位N原子上的甲基化。虽然最近的很多研究显示了m⁶A修饰与关键的生物学功能及病理现象有关，但是m⁶A甲基化的详细调控机制尚不明确。多项研究指出，每种m⁶A甲基化调控因子均同时调控着许多位点的甲基化水平。共享相同调控因子的位点，其甲基化水平通常与调控因子的变化一致，并在功能上表现显著相关，呈现协同甲基化模式。由于生物“个体差异性”及同一个体生物重复性样本“组内差异性”的存在，功能谱分解时更常见的现象应该是“参与同一调控功能的有可能仅是局部修饰位点，而这部分修饰位点也极有可能只在局部条件下表达”。

到目前为止，已经开发了许多方法来进行基因表达谱、DNA甲基化谱等的局部功能谱分解。这些方法的基本思想主要可以分为双聚类以及矩阵分解两类。Lazzeroni等人提出了格子模型(Plaid)，在把基因表达数据看成一个背景模型与多个双聚类叠加的基础上，使用最小化误差平方和来确定双聚类；Ihmels等人提出了迭代签名算法(ISA)，将双聚类定义为一个由基因子集和与之共调控关系最为密切的条件子集一起构成的转录模块；Chen等人提出了基于富集约束框架的加权格子模型(FBCwPlaid)，在对m⁶A甲基化局部功能谱进行分解的过程中引入了RNA表达水平，并利用基因本体论对双聚类结果进行约束；Zhang等人提出的融合同源信息的加权迭代签名算法(REW-ISA V2)，在分析过程中不仅引入了RNA表达水平，还对m⁶A位点的同源性信息进行了充分考虑。Xie等人提出的定性双聚类算法(QUBIC)，通过采用基因表达数据的定性(或半定量)测量和组合优化技术的组合，以更一般的形式解决双聚类问题。Orzechowski等人在unibic算法的基础上提出了runibic包，可以用于挖掘趋势保持的双簇。Engreitz等人基于矩阵分解中独立成分分析(ICA)的思想，将微阵列数据建模为可以共享基因的生物模块的组合，实现了对传统聚类方法的改进。Saelens等人对已有的方法进行了总结和比较，证实了基于矩阵分解思想的方法相对于双聚类方法具有更好的谱分解性能。然而，一方面，已有的基于矩阵分解思想的方法仅在基因表达谱上进行了测试分析，并没有针对RNA甲基化谱进行实际的优化或改进；另一方面，RNA甲基化谱中包含的隐含信息(例如位点的基因组特征等)并未得到充分利用，这导致目前很多方法的谱分解精度都还存在一定的改善空间。

发明内容

本发明目标是揭示m⁶A甲基化的详细调控机制，建立了一种融合基因组特征的m⁶A甲基化局部功能谱分解方法。我们基于矩阵分解中的独立成分分析(ICA)方法，对m⁶A甲基化位点对应的基因组特征进行融合，以便充分考虑不同m⁶A位点之间的关联性。最终利用基于负熵的估计方法对m⁶A位点参与的多个调控通路进行恢复，进一步挖掘出每个调控通路对应的敏感位点，实现对m⁶A甲基化局部功能的高精度谱分解。融合基因组特征的m⁶A甲基化局部功能谱分解方法的具体实施步骤是：

步骤(1)：基于国家生物信息中心(National Center for BiotechnologyInformation,NCBI)的序列阅读档案(Sequence Read Archive，SRA)数据库获取原始测序样本，每组样本均包括对应的输入控制数据及免疫沉淀数据，提取上述两组数据中m⁶A位点的读段数。

步骤(2)：构建各组样本的m⁶A位点甲基化水平矩阵P_N×M＝{p_ij}和表达水平矩阵W_N×M＝{w_ij}，1≤i≤N，1≤j≤M，对W做最大最小归一化使其取值范围为[0,1]。其中，N为样本数目，M为m⁶A位点数目，p_ij和w_ij分别对应于第j个m⁶A位点在第i个样本中的甲基化水平和表达水平。

步骤(3)：基于m⁶A位点的基因组特征，利用斯皮尔曼秩相关系数构建相关系数绝对值大于0.8的极强相关m⁶A位点之间的相关性矩阵，表示为G_M×M＝{g_ij}，1≤i,j≤M。

步骤(4)：构建融合基因组特征的m⁶A甲基化局部功能谱模型，如式(1)所示。

其中A_N×N＝{a_ij}表示待求解的混淆矩阵；S_N×M＝{s_ij}表示调控通路中m⁶A位点的作用矩阵，其中s_ij对应于第j个m⁶A位点在第i个调控通路中的作用程度；α用于调节相关性矩阵的融合权重，其取值范围为0≤α≤1。

步骤(5)：引入辅助矩阵X'＝(P*W)·[αI+(1-α)G]^-1，利用主成分分析对X'白化处理，如式(2)所示。

其中V_N×N＝{v_ij}表示白化矩阵，D_N×N＝{d_ij}是由辅助矩阵X'中各行间的方差矩阵的特征值构成的对角矩阵，E_N×N＝{e_ij}是相应的特征向量的正交矩阵。

步骤(6)：经过白化处理后，方法转化为求解Z＝VAS中的矩阵A及矩阵S。令Y＝WZ，其中W＝A^-1V^-1，采用基于负熵的估计方法，对Y中的独立成分逐一求解。假设当前求解的是y_i，对应于Y中的第i行(即第i个独立成分)，1≤i≤N，本方法中采用的负熵近似表达式如(3)所示。

J(y_i)＝[E{G(y_i)}-E{G(v)}]² (3)

其中，v是与y_i具有相同方差的高斯变量，非二次函数G如式(4)所示。

其中，a为常数且1≤a≤2。求解y_i的具体步骤如下：

第一步，假设当前求解的是对应第i个独立成分中的第k个解混元素w_ik，1≤i,k≤N。求式(3)的梯度，可得负熵近似表达的梯度算法如式(5)所示。

其中，r＝E{G(y_i)}-E{G(v)}，W_i·是矩阵W中的第i个行向量(解混向量)，Z_k·是矩阵Z中的第k个行向量，函数g是非二次函数G的导数，g(y)＝tanh(ay)。w_ik的更新如式(6)所示。

其中，为w_ik更新后的值，函数g'是函数g的导数(非二次函数G的二次导数)，g'(y)＝a(1-tanh²(ay))。

第二步，同理可更新W_i·中的其他(N-1)个解混元素，重复第一步及第二步，直至W_i·收敛，W_i·的收敛条件如式(7)所示。

其中，ε是人为预设的一个很小的正数，是当前轮更新后的第i个解混向量，W_i·是上一轮的第i个解混向量。

步骤(7)：对解混向量W_(i+1)·正交化处理，以提取第i+1个独立成分y_i+1，如式(8)所示。

其中是由已经估计出的i个解混向量W_1·,W_2·,...,W_i·构成的矩阵，它的列依次对应这i个向量。

步骤(8)：重复步骤(6)及(7)完成Y中所有独立成分的提取，利用柯尔莫可洛夫-斯米洛夫检验(Kolmogorov–Smirnov test,KS test)对各独立成分进行正态性检验。

步骤(9)：Y＝WZ，估计Y中各独立成分对应的概率密度，提取95％置信区间以外的m⁶A甲基化位点集合作为敏感位点，因而每个独立成分可提取出两个m⁶A甲基化局部功能模块。

综上，一种融合基因组特征的m⁶A甲基化局部功能谱分解方法相比现有算法，具有如下优点：一方面，本发明基于独立成分分析及中心极限定理的基本思想，将谱分解问题巧妙地转化为求解最大非高斯性的问题，为m⁶A甲基化数据提供了更合理的模型；另一方面，本发明引入m⁶A甲基化位点的基因组特征，在优化过程中充分考虑了极强相关m⁶A位点之间的相互影响，使获得的独立成分更加贴近真实的调控通路，即算法具有更高的谱分解精度。

附图说明

图1是融合基因组特征的m⁶A甲基化局部功能谱分解方法的算法模型流程图

具体实施方式

为了进一步解释本次发明的具体内容和优点，以下是具体实施方式和附图的详细说明。

为了验证本算法在真实数据集上的性能，本实验从NCBI数据库中获得10个公开的人类m⁶A MeRIP-Seq数据集中的32个样本，使用WHISTLE对m⁶A位点进行估计、DESeq2对各m⁶A位点的甲基化水平进行测定。经过上述处理后，共获得了包含32种样本下69446个m⁶A位点的数据，所有这些样本都包含IP和INPUT两组数据。由于部分原始位点的估计精度较低，需要对甲基化位点进行适当的过滤。过滤步骤如下：

1、如果表达水平低于8，或者同一位点(IP+INPUT)样本的读段数低于50，则甲基化水平将被屏蔽为NA。在所有32个条件中，如果出现太多缺失值(NA计数>10)，则应删除位点；

2、由于不同条件之间的差异越大表明功能越明显，因此将保留不同条件下甲基化水平值的中值绝对偏差大于0.4的位点。

经过过滤后，仅保留了32种样本下24086个m⁶A位点的数据。根据过滤后的数据，构建位点的甲基化水平矩阵P以及表达水平矩阵W。进一步通过R中的GenomicFeature包提取出24086个m⁶A位点对应的基因组特征，计算两两位点之间的斯皮尔曼秩相关系数，仅保留极强相关位点之间的相关性，构建矩阵G。在求解过程中，本方法所依赖的预设参数仅包含相关性矩阵的融合程度α、最大迭代次数epoch以及收敛条件ε，初始化参数可以通过网格寻优的方式确定，其中：α＝0.4,epoch＝8000,ε＝5e^-15。获得敏感m⁶A位点后提取其对应的基因，基于包含相等数目大小的随机基因组计算其平均相对提升率(Average relativepromotion rate,ARPR)，以此来评价本次实施例算法的性能。不同相关性矩阵的融合程度α对本发明提出的算法(Independent Component Analysis of Fusion Genome Features,FGFICA)性能的影响如表1所示。此外，本发明提出的算法FGFICA和Lazzeroni等人提出的算法(Plaid)、等人提出的算法(ISA)、Zhang等人提出的算法(REW-ISA V2)、Xie等人提出的算法(QUBIC2)，Orzechowski等人提出的算法(runibic)的性能比较如表2所示。

表1：不同α值下FGFICA的实验结果

表2：基于真实数据集各方法的实验结果

最后说明：上述实施方式是为了更好的说明本发明的思路，绝不是对本发明的限制，凡是根据本发明本质内容所做的等效的替换、修改或补充，均应包含在本发明的保护范围之内。

Claims (2)

1.一种融合基因组特征的m⁶A甲基化局部功能谱分解方法，其包括以下步骤：

步骤(1)：基于国家生物信息中心的序列阅读档案数据库获取原始测序样本，每组样本均包括对应的输入控制数据及免疫沉淀数据，提取上述两组数据中m⁶A位点的读段数；

步骤(2)：构建各组样本的m⁶A位点甲基化水平矩阵P_N×M＝{p_ij}和表达水平矩阵W_N×M＝{w_ij}，1≤i≤N，1≤j≤M，对W做最大最小归一化使其取值范围为[0,1]，其中，N为样本数目，M为m⁶A位点数目，p_ij和w_ij分别对应于第j个m⁶A位点在第i个样本中的甲基化水平和表达水平；

步骤(3)：基于m⁶A位点的基因组特征，利用斯皮尔曼秩相关系数构建相关系数绝对值大于0.8的极强相关m⁶A位点之间的相关性矩阵，表示为G_M×M＝{g_bd}，1≤b,d≤M；

步骤(4)：构建融合基因组特征的m⁶A甲基化局部功能谱模型，如式(1)所示；

其中A_N×N＝{a_ef}表示待求解的混淆矩阵，a_ef对应于矩阵A中第e行第f列下的元素取值；S_N×M＝{s_op}表示调控通路中m⁶A位点的作用矩阵，其中s_op对应于矩阵S中第o行第p列下的元素取值，式(1)中矩阵S’由S和G加权求和得到；α用于调节相关性矩阵的融合权重，其取值范围为0≤α≤1；

步骤(5)：引入辅助矩阵X'＝(P*W)·[αI+(1-α)G]^-1，其中I表示单位矩阵，利用主成分分析对X'白化处理，如式(2)所示：

其中V表示白化矩阵，Z为经白化处理后的待分析矩阵，D是由辅助矩阵X'中各行间的方差矩阵的特征值构成的对角矩阵，E是相应的特征向量的正交矩阵；

步骤(6)：经过白化处理后，方法转化为求解Z＝VAS中的矩阵A及矩阵S，令Y＝WZ，其中W＝A^-1V^-1，采用基于负熵的估计方法，对Y中的独立成分逐一求解；假设当前求解的是y_h，对应于Y中的第h行，本方法中采用的负熵近似表达式如(3)所示；

J(y_h)＝[E{G(y_h)}-E{G(v)}]² (3)

其中，E{}表示计算数学期望值，v是与y_h具有相同方差的高斯变量，非二次函数G如式(4)所示：

其中，log cosh为双曲余弦函数的对数，a为常数且1≤a≤2；求解y_h的具体步骤如下：

第一步，假设当前求解的是对应第l个独立成分中的第k个解混元素w_lk，1≤l,k≤N；求式(3)的梯度，可得负熵近似表达的梯度算法如式(5)所示：

其中，r＝E{G(y_h)}-E{G(v)}，W_l·是矩阵W中的第l个解混向量，Z_k·是矩阵Z中的第k个行向量，w_lk的更新如式(6)所示；

其中，为w_lk更新后的值；

第二步，同理可更新W_l·中的其他(N-1)个解混元素，重复第一步及第二步，直至W_l·收敛，W_l·的收敛条件如式(7)所示；

其中，ε是人为预设的一个很小的正数，是当前轮更新后的第l个解混向量；

步骤(7)：对解混向量W_(l+1)·正交化处理，以提取矩阵Y中第l+1行，即y_l+1，如式(8)所示：

其中是由已经估计出的l个解混向量W_1·,W_2·,...,W_l·构成的矩阵，它的列依次对应这l个向量；

步骤(8)：重复步骤(6)完成独立成分矩阵Y中所有独立成分的提取，利用柯尔莫可洛夫-斯米洛夫检验对各独立成分进行正态性检验；

步骤(9)：Y＝WZ，估计Y中各独立成分对应的概率密度，提取95％置信区间以外的m⁶A甲基化位点集合作为敏感位点，因而每个独立成分可提取出两个m⁶A甲基化局部功能模块。

2.根据权利要求1所述的融合基因组特征的m⁶A甲基化局部功能谱分解方法得到的独立成分矩阵Y。