基于子空间学习的m7G与疾病的关联预测方法

基于子空间学习的m7g与疾病的关联预测方法
技术领域
1.本发明涉及生物信息学研究领域，具体是用基于子空间学习的m7g与疾病的关联预测方法，为表观转录组层面探究肿瘤发生机制和致病通路提供参考。本发明有望为以m7g修饰位点为靶点的肿瘤治疗策略提供新思路。


背景技术：

2.7-甲基鸟苷(n7-methylguanosine,m7g)是一种在鸟苷(guanosine,g)的n7位置出现的一种rna修饰，是信使rna(message rna,mrna)的5’端帽结构(5’capping)的一部分，广泛存在于多个物种的各类rna 中。近年来，学界研究表明m7g在生物学过程及正常生理功能中发挥重要的调节作用。其中，研究发现 m7g通过rna的转录、剪切、加工、翻译和衰变等过程参与多种疾病甚至是恶性肿瘤的发生和转移。shaheen 等人发现m7g转移酶复合物mettl1-wdr4上的突变会导致trna上m7g位点数量减少，从而诱发原发性侏儒症和脑发育畸形。林等人发现敲除掉wdr4基因会破坏胚胎干细胞的分化，并导致神经系统紊乱。此外，邓等人发现mettl1基因的沉默会引发多能干细胞上m7g位点数量减少，使中胚层的分化和血管发育受到了严重抑制。另外，最近在肺癌、结肠癌和肝癌细胞中也报道了m7g位点的异常。因此，如果能够进一步确定与某一疾病相关的m7g修饰位点，则将为以m7g位点为靶点的肿瘤治疗策略提供新思路。虽然湿实验能够提供较为准确的m7g-疾病关联关系，但是其过程较为复杂且耗费大量的人力物力财力。所以，我们以现有的记录m7g-疾病关联的生物信息数据库为依托，用计算方法预测未知的m7g-疾病关联，为生物学家下一步湿实验的开展提供关键的m7g位点。
3.虽然目前解决m7g-疾病关联预测的算法很少，但是从数学角度上看，实现m7g-疾病关联预测的本质是填充m7g位点与疾病的关联矩阵，即矩阵填充问题。近年来，子空间学习作为一种捕获矩阵潜在结构的方法为解决矩阵填充问题提供了新思路，在生物信息学领域引起了广泛的关注。jolliffe等人提出的主成分分析(principal component analysis，pca)是应用最广泛的子空间学习算法，并在高斯噪声下取得了显著的效果。pca认为给定的数据来自一个低维线性子空间和一个噪声空间，通过求解核范数最小化和 frobenius范数最小化问题来学习去噪的低维线性子空间，以此恢复矩阵中的缺失元素。但在实际情况下，给定的观测数据并不一定源于同一个低维完备子空间，而是源于多个。为了从噪声中捕获多个低维完备子空间的结构特性，刘等人提出了低秩表示(low-rank representation,lrr)模型。lrr认为在一组样本中，每个样本都可以表示为子空间中基的线性组合。进一步地，lrr以子空间自身为字典，旨在找到所有样本的低秩表示。在样本划分的过程中，lrr将原始的高维空间分割成了多个低维子空间。lrr模型是一种核范数最小化模型(nuclear norm minimization,nnm)，它采用核范数来对秩函数进行凸松弛。作为一个凸问题， nnm能够得到一个全局最优解，但是仍旧偏离原本秩最小化问题的解。这归因于两个原因：一是核范数并不是秩函数的最佳近似，二是用奇异值阈值算法(singular value thresholding,svt)时，无论奇异值大与小，svt对他们的惩罚和收敛程度相同。但是大的奇异值包含更多
信息，小的奇异值则含有很多噪声。所以，理论上应该对小的奇异值加以更大的惩罚，使其以较快的速度收敛；对大的奇异值加以更小的惩罚，使其以较慢的速度收敛。近年来，schatten-p范数逐渐应用在矩阵填充领域。相比于核范数，schatten-p范数是一个对秩函数更紧的非凸近似。nie等人最先提出基于schatten-p范数的lrr模型，即schatten-p范数最小化模型。其在矩阵填充场景下，获得了比传统lrr模型更高的预测准确率和更好的鲁棒性。接着，张等人提出了求解schatten-p范数最小化问题的广义矩阵软阈值算子(generalized matrix soft thresholding, gmst)。gmst能够让小的奇异值收敛速度更快，大的奇异值收敛速度变慢。实验表明，相比于svt， gmst能够更好地刻画子空间的结构。尽管schatten-p范数最小化问题能够更有效地挖掘子空间的全局结构特性，但是却忽略了子空间的局部结构特性。另外，稀疏子空间聚类算法(sparse subspace clustering,ssc)可以挖掘子空间的局部稀疏表达特性。形式上，ssc将每个数据点表示成子空间中其他数据点的线性组合，以此划分高维空间，但ssc无法捕获关联矩阵的整体性质。有效的m7g-疾病关联预测算法有待提出，且作为目前取得较好效果的子空间学习算法，如何在矩阵存在缺失元素的前提下，明晰其关联模式，仍是亟待改善的问题。


技术实现要素：

4.本发明目标是建立基于子空间学习的m7g-疾病关联预测方法，简称spblrsr。以“相似的位点会引发相似的疾病”这个假设为依托，本发明利用已知的m7g-疾病关联信息、m7g位点相似性信息和疾病相似性信息，从结构缺失、关联模式不清晰的高维异构子空间出发，习得体现其统计及结构特性的完备低维子空间，进而恢复不完备的高维关联矩阵，提取潜在的疾病相关m7g位点。其具体实施步骤是：
5.步骤(1)：从m7gdisai数据库中获得m7g-疾病关联矩阵h
sd
(741*177)，m7g位点化学相似性矩阵 h
ss_chemical
(741*741)，m7g位点碱基累计频率(cumulative nucleotide frequency,cnf)相似性矩阵h
ss_cnf
(741*741)和疾病相似性矩阵h
dd
(177*177)。其中：
[0006][0007]
步骤(2)：构建m7g位点综合相似性矩阵h
ss
：
[0008]hss
＝αh
ss_chemical
(1-α)h
ss_cnf
,0≤α≤1
ꢀꢀ
(1)
[0009]
其中，α是m7g位点化学相似性和m7g位点碱基累计频率相似性的组合系数。
[0010]
步骤(3)：构建m7g-疾病异构矩阵。m7g-疾病关联矩阵的稀疏度很高，难以保证算法的正常求解。同时，高稀疏度使得本身就存在缺失元素的矩阵的关联模式更不清晰。为了明晰关联模式，并保证算法能够正常求解，构建m7g-疾病异构矩阵，如(2)式所示。
[0011][0012]hsdt
是h
sd
的转置。m为m7g位点的个数，n为疾病的个数，则m7g-疾病异构矩阵h的维度是(m n)* (m n)。由(2)可知，h具有良好的性质。首先，h是一个对称且半正定的矩阵，其奇异值是正实数，并且等于特征值。且其左奇异向量和右奇异向量相等，等于特征向量。半正定矩阵奇异值和奇异向量的特殊性大大的减少了计算的复杂度，加快了算法运行的速
度。其次，m7g-疾病异构矩阵h的稀疏度比关联矩阵h
sd
小很多，为算法的正常求解提供保证。最后，h中的缺失值只出现在h
sd
和其转置中，因此可将 h
sd
的填充问题转化为h的填充问题。
[0013]
1.步骤(4)：建立spblrsr模型，意在从不完备的高维子空间中习得能够保存其统计和结构特性的完备低维子空间，从而预测矩阵中的缺失关联，为生物学家下一步湿实验的开展提供关键的m7g位点。这里用x来表示完备低秩异构矩阵，用c来表示稀疏自表达矩阵，将schatten-p范数最小化模型和ssc模型结合起来，建立目标函数(3)：
[0014][0015]
其中，x,c遵从上述定义，分别是预测出的完备低秩异构矩阵和稀疏的自表达矩阵。λ是平衡系数。 ||x-xc||
f2
考虑了矩阵恢复过程中可能出现的噪声，diag(c)＝0的约束是为了避免在c的求解过程中出现平凡解c＝i。p
ω
是投影算子，其定义如下：
[0016][0017]
步骤(5)：引入连个辅助矩阵m和a，得到模型(3)的等价形式(4)。进而，采用交替方向乘子法(alternative direction method of multiplier,admm)求解，具体过程如下：
[0018][0019]
上式的增广拉格朗日函数如(5)所示：
[0020][0021]
其中，μ是惩罚因子，y1和y2拉格朗日乘子。在第k 1次循环，spblrsr需要交替求解子问题x
k 1
,m
k 1
, a
k 1
,c
k 1
,y
1k 1
和y
2k 1
.
[0022]
x的更新：固定mk,ak,ck,y
1k
,y
2k
,最小化增广拉格朗日函数：
[0023][0024]
记在f(x)中对x求偏导，极小值x
*
的解析解如(6)所示。
[0025]
x
*
＝(y
2k
μmk)(i-((ak)
t
ak) ak(ak)
t
μi)-1
ꢀꢀ
(7)
[0026]
进而，加上数值范围和投影矩阵的限制，x
k 1
的解析解如(8)所示：
[0027]
x
k 1
＝l
[0,1]
(x
*
),p
ω
(x
k 1
)＝p
ω
(h)
ꢀꢀ
(8)
[0028]
其中l
[0,1]
定义为：
[0029]
[0030]
m的更新：固定x
k 1
,ak,ck,y
1k
,y
2k
,最小化增广拉格朗日函数：
[0031][0032]
上述模型是schatten-p范数最小化问题，其解析解由gmst算子给出：
[0033][0034]
其中，是gmst算子。u,v分别是矩阵的左、右奇异矩阵，ξ＝(ξ1,ξ2,...,ξr)是的奇异值。是广义向量软阈值算子(generalized vector soft thresholding,gvst)，它的基本构成是广义软阈值算子(generalized soft thresholding,gst)：
[0035][0036]
对于每一个gst算子有：
[0037][0038]
其中，且通过求解得到。
[0039]
a的更新：固定x
k 1
,m
k 1
,ck,y
1k
,y
2k
,最小化增广拉格朗日函数，对a求偏导，a的解析解如(11)所示：
[0040]ak 1
＝((x
k 1
)
t
x
k 1
μi)-1
((x
k 1
)
t
x
k 1-y
1k
μck)
ꢀꢀ
(11)
[0041]
c的更新：固定x
k 1
,m
k 1
,a
k 1
,y
1k
,y
2k
,最小化增广拉格朗日函数,转化成l1范数最小化的问题，其解析解可由软阈值算子表示：
[0042][0043]ck 1
＝c
k 1-diag(c
k 1
)
[0044]
其中，是软阈值算子。
[0045]
y1，y2的更新:固定x
k 1
,m
k 1
,a
k 1
,c
k 1
,由梯度下降法得：
[0046][0047]
其中，δ是学习率。
[0048]
步骤(6)：重复上述过程直至算法收敛。完备低秩m7g-疾病异构矩阵x中的完备低秩关联矩阵x
sd
为最终所求。
附图说明
[0049]
图1是基于子空间学习的m7g与疾病的关联预测模型流程图。图2是spblrsr，mf，svt及rpca方法在“留一疾病交叉验证”机制下的核密度估计图。针对177个疾病，各算法在该机制下会产生177个auc值。该图展示了各方法预测出的177个auc的分布情况,说明了各算法为疾病提供相关m7g位点的能力。
具体实施方式
[0050]
为了进一步解释本发明的具体内容和优点，以下是具体实施方式和附图的详细说明。
[0051]
为了说明本算法的有效性，我们在两种交叉验证框架下比较spblrsr和其他三种在矩阵填充领域广泛应用的算法，分别是矩阵分解算法(matrix factorization,mf),奇异值阈值算法(singular value thresholding,svt)和鲁棒性主成分分析算法(robust principal component analysis,rpca)。首先，我们在十折交叉验证(10-fold cross validation,10-fold cv)实验中比较他们的预测准确度。m7gdisai记录了768对 m7g位点和疾病的关联，这些关联涉及741个m7g位点和177个疾病。我们把130389个未知关联样本当作候选样本(candidate set)。另外，把768个已知关联以近乎相同的规模随机分成独立的十份，每次都使用九份作为训练集(training set)，剩下一份作为测试集(test set)。算法在训练集上训练参数，在测试集上和候选集上应用，得到预测分数。接着，我们采用auc(area under roc curve)和aupr(area under pr curve)来衡量四个算法的预测能力，结果如表1所示。最后，用wilcoxon signed-rank test来检验算法预测结果差异的显著性。表2展示了mf，svt和rpca与spblrsr算法预测结果的差异显著性。
[0052]
表1：10-fold cv下各算法的预测效果。
[0053][0054]
表2：10-fold cv下mf,svt和rpca与spblrsr的差异显著性。
[0055][0056]
实验表明，spblrsr的预测能力显著地优于mf,svt和rpca。
[0057]
其次，对于一种疾病而言，更普遍的情况是不知道任何与它相关的m7g位点。这里，
我们使用“留一疾病交叉验证”法(leave-one-disease-out-cross-validation，lodocv)来验证spblrsr为新疾病预测相关的m7g位点的能力。同10-fold cv一样，我们把130389个未知关联当作候选样本(candidate set)。不同的是，对于疾病d，我们将关联矩阵中所有已知与d相关的m7g位点都作为测试样本，而将m7g位点与其他疾病之间的已知关联视作为训练样本。最后，用训练集训练算法，再将训练好的算法作用在测试集和候选集上。我们对所有的177种疾病，都做上述的实验，进而获得177个auc值。图2展示了各算法lodocv实验的结果。从图2中可以看到，mf、svt和rpca产生的auc得分主要分布在0.3 和0.6之间，而spblrsr得到的auc得分集中在0.6和0.85之间。说明spblrsr在lodocv实验中优于mf,svt和rpca，也说明spblrsr在某疾病没有相关m7g位点存在的前提下，能够更准确地预测出与之相关的m7g位点。
[0058]
综上，spblrsr比mf、svt和rpca在预测m7g-疾病关联上更有效，表现出了更高的预测准确性和为疾病提供相关m7g位点的能力。
[0059]
最后说明：上述实施方式是为了更好的说明本发明的思路，绝不是对本发明的限制，凡是根据本发明本质内容所做的等效的替换、修改或补充，均应包含在本发明的保护范围之内。