基于融合相似性和三部图的 circRNA 与疾病关联预测

王波，刘庭斌，张剑飞，杜晓昕，王鑫炜

(齐齐哈尔大学 计算机与控制工程学院，黑龙江 齐齐哈尔 161006)

circRNA 是具有共价闭环结构的内源性非编码RNA，最初被认为是RNA 错误剪接的副产物[1].随着高通量测序方法和生物信息学的飞速发展，多项研究已证明circRNA 与其他分子存在相互作用，参与生物体内各项生命活动的调控，例如靶基因表达、细胞增殖、免疫应答、遗传印迹、肿瘤侵袭等[2-5].与circRNA 相关的各种生物学数据存储在公共数据库中，如circBase[6]、CircR2Disease[7]、circRNADisease[8]和circ2Traits[9]等.

在“相似的circRNA 可能与相同的疾病有相似的关联”假设下，许多计算模型被用于挖掘潜在的circRNA-疾病关联，解决了传统生物实验耗时长且高成本的问题[10].这些模型可以大致分为3 类：基于信息在网络中的传播、基于机器学习和基于深度学习.Fan 等[11]提出使用异构网络的路径信息进行circRNA-疾病关联预测的KATZ 度量计算模型（KATZHCDA）.Li 等[12]提出基于网络一致性投影的计算方模型（NCPCDA），利用多源相似性和一致性投影得到预测得分矩阵.Ding等[13]结合随机游走算法和逻辑回归方法开发了名为RWLR 的计算模型.Lei 等[14]提出名为ICFCDA 的基于协作过滤推荐系统的计算模型.Deepthi 等[15]提出计算模型AE-DNN，它依赖自动编码器和深度神经网络来预测新的circRNA-疾病关联.Xiao 等[16]提出基于网络嵌入的自适应子空间学习方法（NSL2CD），同时在模型中加入综合加权图正则化项和L1 范数约束来实现投影矩阵的平滑性和稀疏性.不难看出，已有的计算模型还存在一些缺陷：1）模型使用的训练数据有限，这对模型的鲁棒性和覆盖范围有影响；2）模型主要基于单一的数据描述方法，没有将circRNA 与疾病行为信息和属性信息结合起来，全面定义circRNA 与疾病的特征，导致预测性能有限；3）研究者没有考虑编码-非编码基因-疾病关联的异质性，无法准确测量circRNA-disease关联信息.

为了改善现有计算模型不足，本研究提出基于融合相似性和三部图的circRNA 与疾病关联预测模型（prediction of circRNA and disease association based on fusion similarity and tripartite graph, FSTPGCDA）.研究工作包括1）利用数据库得到circRNA 序列信息、 circRNA-gene 关联信息、circRNA-disease 关联信息和疾病语义信息，把数据处理成circRNA-disease关联矩阵和circRNA-gene 关联矩阵.2）利用混沌博弈表示（chaotic game representation，CGR）[17]、语义相似性、Jaccard 系数[18]与拉普拉斯特征映射[19]融合相似性计算相似性.3）加权相似性得到融合相似性.4）利用circRNA-disease 关联信息和circRNA-gene 关联信息构建gene-circRNA-disease三部图[20].5）通过融合相似性方法为三部图分配初始资源，使用贪心算法进行资源分配，得出最终circRNA-disease 资源得分矩阵.6）计算预测得分并排序，进行留一交叉验证（leave-one-out crossvalidation, LOOCV）[21].

1 三部图模型

1.1 数据集

通过整合不同种类的生物关联信息，构建数据集D1、D2.在D1 中，circRNA-diseas 关联从CircFunBase 数据库[22]中下载；分别从circBase、circR2Disease 和MeSH[23]中收集circRNA 序列信息、circRNA-gene 关联信息和疾病语义信息；剔除重复后，共收集2 983 个circRNA-diseas 关联和2 318 个circRNA-gene 关联信息.在D2 中，circRNAdiseas 关联信息从circR2Cancer[24]数据库中下载；分别从circBase、miR2Disease[25]和MeSH 中收集circRNA 序列信息、circRNA-miRNA 关联信息和疾病语义信；剔除重复后，共收集到647 个circRNAdiseas 关联信息和756 个circRNA-miRNA 关联信息.数据集的关联信息及数据个数n如表1 所示.

表1 数据集关联信息Tab.1 Data set association information

1.2 相似性计算

1.2.1 融合相似性计算 本研究的数据集存在稀疏问题，使得计算的过程时间长，算法时间复杂度高.与传统相似性度量方法相比，Jaccard 相似性能够改善余弦相似性[26]只考虑单一变量而忽略其他信息量的弊端，适合在稀疏度过高的数据中使用.拉普拉斯特征映射是基于图的降维算法，在降维后仍能保持原有的数据结构.本研究将拉普拉斯特征映射和Jaccard 结合进行相似性计算.基本思路如下.

1）以计算样本关联矩阵相似性为例，a、b关联矩阵记为M，元素Mi,j=1 为第i个a样本和第j个b样本存在关联.将M拆分成2 个矩阵A和B，操作方式如下：将M每 行中值为1 的列下标分别作为对应的b样本的属性，构成行向量，该行向量的元素数量nB表示b样本的数量.例如，如果第i行中值 为1 的 列 下 标 为 {j1,j2,j3} ，那 么Ai=[0,0,0,···,0,1,0,···,0,1,0,···,0,0,0]中1 的位置分别对应第j1,j2,j3个b样本的属性，其余位置为0.将所有nA个向量按行堆叠，得到nA×nB的矩阵，其中nA为a样本的数量.矩阵B同理可得.其中Ai,j=1 为 第i个a样 本 具有第j个 属 性，Bi,j=1 为 第i个b样本具有第j个属性.

2）使用拉普拉斯特征映射将A和B映射到特征空间中，得到kA个a样本的特征向量和kB个b样本的特征向量.使用Jaccard 系数计算这些特征向量之间的相似性，即计算所有非零元素的索引集合的Jaccard 系数的平均值，即

式 中：J(A,B) 为 矩 阵A和B的 Jaccard 系 数，φA(ai) 为将a样本i映射到特征空间中得到的特征向量，φB(aj) 为 将b样本j映射到特征空间中得到的特征向量.

1.2.2 circRNA 相似性计算 现有序列比对算法只能量化位置信息或非线性信息，能够将这2 类信息结合的算法鲜少.为此基于CGR 的方法利用Pearson 相关系数[27]来量化位置与非线性信息之间的相似性和差异性.1）将CGR 空间划分为Ng网 格(Ng=2s×2s;在本研究中s=3 )，网格表示为

2）分别对每个网格中横坐标x和纵坐标y进行累加，若点在网格内，量化位置信息为

3）计算每个网格Zi的评分量化非线性信息，其中Ni为第i个网络的评分，

4）每个网格被描述为3 个属性, 并融合属性构造描述第i个circRNAc(i).通过Pearson 相关系数确定序列相似度(c(i),c(j)).

式中： C ov 为协方差，D为方差.circRNA 序列相似度矩阵大小为n×n.

如果RNA 影响同一种人类疾病，它们的功能往往是相似的[28].从circRNA-gene 关联矩阵中利用融合相似性计算得出circRNA-gene 之间circRNA 的 相 似 性.从circRNA-disease 关 联 矩阵中利用上述的融合相似性计算得出circRNAdisease 之间circRNA 的相似性.

通过从不同角度分析circRNA 的特征，可以得到3 个相似矩阵，包括和.为了解决由于数据和可用性不足导致的circRNA 之间差异较小的问题，引入序列相似度和circRNAgene 从不同的角度来描述circRNA，使模型更具信息量.相似性信息的完成是通过组合3 个矩阵来完成的.融合相似度定义如下：

1.2.3 disease 相似性计算 疾病语义相似.根据MeSH 数据库的语义信息将疾病表示为有向无环图（directed acyclic graph，DAG）.DAG 中的节点代表疾病，边代表疾病之间的关系.如果疾病在病理上相似，则DAG 的更多部分将被共享[29].计算疾病贡献值的模型为

式中：n(DAGs(q)) 为疾病q的D AGs的数量定义，并将所有疾病的数量定义为nd.疾病d(i) 与疾病d(j)的语义相似度评分描述为

式中：Nd(i)为疾病d(i) 的DAG 出现的所有疾病.病理相似的疾病往往与功能相似的circRNAs

1.2.4 gene 相似性计算 从circRNA-gene 关联矩阵中利用融合相似性计算得出circRNA-gene 之间gene 的相似性.

1.3 gene-circRNA-disease 三部图

1.3.1 三部图介绍 由circRNA-disease 和circRNA-gene构建三部图，加入circRNA 相似性和disease 相似性来对节点进行资源分配.

图1 基于融合相似性和三部图的circRNA 与疾病关联预测模型流程图Fig.1 Flow chart of circRNA and disease association prediction model based on fusion similarity and tripartite graph

1.3.2 三部图资源推荐 三分网络资源分配的贪心算法[30]流程如下.1）对每个节点初始化资源.根据相似性矩阵SC、SD、SG，为gene-circRNAdisease 中每个节点分配资源，分别为节点C、D、G，分配初始资源RC、RD、RG.例如，对于特定的ci，位于疾病dj上的初始资源定义为

2）对于每个节点，计算与其相邻的节点之间的收益值：

式中： G aini,j为节点i和节点j之间的收益值，wi,j为 节 点i和 节 点j之 间 的 相 似 性，Ri、Rj分 别 为 节 点i和节点j目前拥有的资源数量.3）选择最优节点，即从未被选择的节点中, 选择与已选择节点之间的收益值最大的节点加入已选择节点集合中.4）更新资源，根据新加入的节点，更新其相邻节点的资源数量：

构建三部图时存在的孤立节点的处理过程：在贪心算法的过程中，比较孤立节点与已有资源的节点，并将资源分配给孤立节点的邻居节点，以提高整个网络的连通性.计算每个节点的收益值：

式中：v为 当前节点，Nv为v的 邻居节点集合，sim(v,u)为 节 点v和u之 间 的 相 似 度，ru为 节 点u已 经 分 配 的资源值.G ainv的计算结果代表如果将资源分配给节点v整个网络的收益会增加多少.

2 实验结果与分析

2.1 评估指标

FSTPGCDA 推断circRNA-disease 之间潜在关联的性能使用L O O C V 进行评估.在每个LOOCV 实验的步骤中，将每个已知的circRNAdisease 关联从训练样本中依次移除，作为测试样本，其他已知关联作为模型学习的训练样本.定义TP 为被预测为正类的样本，TN 为被预测为负类的负样本，FP 为被预测为正类的负样本，FN 为被预测为负类的正样本.根据ROC 曲线下的相应面积计算AUC（area under curve）[31].根据PR 曲线计算曲线下面积AUPR（area under precisionrecall curve）[32].

式中：n为正负样本总数，p为精密度.

2.2 模型预测能力评估

自身效果对比采用AUC、AUPR、TPR、精密度、F1 评分和MCC 评估指标，分别用LOOCV 和5、10 折对比.对比结果如表2 所示.可以看出，各评估指标在交叉验证中差异不超0.1%，该模式具有较好的鲁棒性.

表2 所提模型在不同测试方法下的评估指标对比Tab.2 Comparison of indicators for proposed model at different test methods

2.3 本研究模型与其他模型的比较

进行FSTPGCDA 与KATZHCDA、iCDA-CMG[33]、DMFCDA[34]、KGANCDA[35]和DRGCNCDA[36]的模型性能对比实验.不同模型的ROC 曲线如图2所示.KATZHCDA、iCDA-CMG、DMFCDA、KGANCDA 和DRGCNCDA 的AUC 分别为84.69%、86.25%、88.61%、87.14%和93.99%，FSTPGCDA 的AUC 为97.01%，优于其他5 种方法.不同方法的PR 曲线如图3 所示.KATZHCDA、iCDA-CMG、DMFCDA、KGANCDA 和DRGCNCDA 的AUPR分别为69.44%、79.43%、55.33%、1.20%和2.07%，FSTPGCDA 的AUPR 为86.26%，优于其他5 种方法.

图2 不同模型的ROC 对比Fig.2 Comparison of ROC for different models

图3 不同模型的PR 对比Fig.3 Comparison of PR for different models

2.4 不同相似性在本研究模型中的对比实验

为了验证不同相似性在FSTPGCDA 的效果，将融合相似性与Jaccard 相似性、余弦相似性、高斯核相似性和Pearson 相似性对比，不同相似性在dataset1 的相应ROC 曲线如图4 所示.Jaccard相似性、余弦相似性、高斯核相似性和Pearson 相似性的AUC 分别为96.16%、92.42%、91.47%和80.25%，融合相似性的AUC 为97.01%，优于其他相似性计算.

图4 所提模型不同相似性的ROC 对比Fig.4 ROC comparison of different similarity for proposed model

2.5 案例研究

为了评估FSTPGCDA 的实用价值，进行膀胱癌案例研究，膀胱癌在circFunBase 和circR2Disease数据集中有足够的数据，能够避免模型缺陷导致的偏差.根据相应的预测得分，通过按降序排列选择前15 个分数，并通过PubMed 进行验证，验证结果为PMID 号.膀胱癌是发生在膀胱黏膜上的恶性肿瘤，是泌尿系统最常见的恶性肿瘤，占中国泌尿生殖系肿瘤发病率的第一位[37].膀胱癌筛查模型的研究，对膀胱癌早期发现和高危人群预警具有重要意义.在癌前病变阶段进行筛检，早诊早治，可降低膀胱癌发病率和病死率.进一步研究膀胱癌与circRNA 之间关联有助于提高膀胱癌的诊断和治疗水平.选择预测分数前15 名的CircRNA 进行验证，有14 个得到验证.如表3 所示，hsa_circ_0001946（排名第1）对应的CDR1 基因，与Purkinje 细胞质抗原34 和62 kd 反应的抗Yo(I 型)自身抗体在一例膀胱移行细胞癌并发副肿瘤性小脑变性和抗Yo 抗体反应的患者的血清和脑脊液中被发现.肿瘤切除后抗体滴度下降[38].hsa_circ_0028173（排名第8）在膀胱癌细胞中，这些 DEmRNA 在甘油酯代谢、p53 信号通路和卵母细胞减数分裂中显着富集.circRNA 相互作用对可能在BC 中发挥重要作用[39].hsa_circ_0000144（排名第9）下调环状RNA hsa_circ_0000144 通过刺激miR-217 和抑制RUNX2 表达抑制膀胱癌进展[40].以CDR1 基因为例进行进一步分析，验证该基因是否与膀胱癌相关.如图5 所示，在研究中，将所有膀胱癌患者样本分为高表达组和低表达组，通过生存分析看到CDR1 基因高表达组膀胱癌患者的生存天数相对较短.图中，TS为生存时间，PS为生存概率.如图6 所示，进一步的结果表明，这些基因在癌症样本中的表达明显低于正常样本.图中，R为每百万份转录数.基于以上结果，最终得出这些基因的表达与膀胱癌患者的生存时间和临床病理特征显著负相关.此外，BLCA 富集分析也显示，CDR1 基因低表达组对人类来说主要在蛋白质消化吸收、EMC 受体相互作用、心肌病、癌症中枢碳代谢、黑色素瘤等疾病过程中富集，如图7 所示.图中，ER 为富集率.

图5 CDR1 基因在胃癌患者的生存分析图Fig.5 Survival analysis of CDR1 gene in patients with gastric cancer

图6 CDR1 基因在正常和肿瘤样本中的分化表达Fig.6 Differentiation and expression of CDR1 gene in normal and tumor sample

图7 免疫缺陷中富集的基因集Fig.7 Gene set enriched in immune deficiency

表3 前15 个与膀胱癌有关联的circRNATab.3 First 15 circRNAs associated with bladder cancer

3 结 语

在生物医学研究中，预测circRNA 与疾病关联有利于理解疾病的发病机制，进一步提高疾病诊断、治疗、预后和预防的质量.本研究提出新的计算模型FSTPGCDA，通过整合实验验证的circRNA 序列信息、circRNA-disease 关联信息、circRNA-gene 关联信息和疾病语义信息来识别潜在的circRNA-disease 关联.基于gene-circRNA-disease三部图的资源分配模型，能够更好地描述编码非编码基因疾病关联的异质性，丰富资源分配过程中的生物信息.先用各种相似性和融合相似性计算得到各自相似性，通过加权融合相似性得到融合相似性矩阵，解决了数据高度稀疏性，以便更好地预测得分.利用融合相似性矩阵为三部图分配初始资源，利用贪心算法为三部图进行资源分配，产生推荐该算法有效地减少了资源分配过程中的不可知偏差.在LOOCV、5 折和10 折对比实验中，不同评估指标的评估结果表明，相比其他参与对比的模型，FSTPGCDA 具有较好的预测能力和鲁棒性.案例研究的分析进一步证明，FSTGPCDA 有助于在实践中识别潜在的circRNA疾病关联.下一步计划整合gene-disease 关联或采用的其他生物信息，增加更多的经过实验验证的circRNA-disease 关联，从资源分配方面进行致病机制的具体分析.