基于融合相似性和三部圖的 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 關聯，從資源分配方面進行致病機制的具體分析.