基于圖注意力網絡的環狀RNA與疾病關聯關系預測

張瀚元,趙博偉,胡 倫*,王 磊,尤著宏

(1.中國科學院大學 中國科學院新疆理化技術研究所,新疆 烏魯木齊 830011;2.廣西科學院 大數據與智能計算研究中心,廣西 南寧 530007;3.西北工業大學 計算機學院 大數據存儲與管理工業和信息化部重點實驗室,陜西 西安 710072)

0 引 言

環狀RNA是一類收尾相連具有環狀結構的轉錄RNA,它產生于DNA轉錄過程或轉錄后的修飾[1-2],具體的產生機制還在研究中。雖然細胞內的RNA主要是以線性結構為主,但環狀RNA也大量存在,并且發現環狀RNA往往會高表達轉錄。近年來隨著高通量測序技術的發展,環狀RNA能夠通過被反向比對的雙端(two-paired)短序列識別出。數據分析表明,它們在癌癥等多種復雜疾病組織與正常組織的比對中有顯著的轉錄差異,這些有差異的環狀RNA被認為與該疾病發生和發展有關系[3]。比如,Hsa_circ_0046430在最近研究中參與miR-6785-5p/SRCIN1的ceRNA調控網絡促進結腸癌的生長[4],CircRNA DDX21則參與miR-1264/QKI的ceRNA調控網絡以弱化三陰性乳腺癌的生長[5],而利用環狀RNA基因表達數據則可以挖掘出新的胃癌標志物[6]等等。然而,通過實驗手段發現的環狀RNA與疾病的關聯關系畢竟有限,研究人員希望通過現有的研究信息和生物知識,借助機器學習和人工智能的方法,預測和挖掘環狀RNA與疾病的關聯關系[7]。

1 研究背景

建立生物基因型與表型性狀的關聯關系一直是生命科學研究的重要問題[8]。研究人員已經通過計算手段來挖掘這種關聯關系,如小RNA(microRNA)與疾病[9]、非編碼RNA(LncRNA)與疾病[10]、環狀RNA(Circular RNA,CircRNA)與疾病[11-13]的關聯關系。由于已有知識的局限,以及不同生物分子對應的疾病特征不同,目前多數有效的環狀RNA與疾病的預測方法是通過鏈路預測(Link prediction)對已知的環狀RNA與疾病關系的補全,關聯關系(Association prediction)預測可以認為是鏈路預測的一種特例[14]。主要關于鏈路預測的方法都有嘗試在環狀RNA與疾病關系預測問題上進行研究,比如KATZHCDA方法通過KATZH圖信息指標對環狀RNA與疾病的關系進行預測。KATZH指標是一種通過節點間鏈路個數來衡量節點間關系程度并用于鏈路關系的預測[15]。iCircDA-MF通過矩陣分解的方法整合環狀RNA與疾病的相關信息進行鏈路預測[16]。也有通過深度學習模型構建分類器進行相關關系的預測,如MSFCNN方法通過融合多源信息后利用兩層卷積網絡進行關系預測[17]。GIS-CDA也是一種采用了圖注意力機制的模型,但主要是利用數據融合的技術和歸納式矩陣補全[12]。以上關于圖鏈路預測的傳統方法都有應用在環狀RNA與疾病關系的預測中。AE-DNN方法通過構建編碼器(AutoEncoder)和深度神經網絡(Deep Neural Network)進行關系預測[18]。AANE-SAE[19]利用屬性網絡編碼算法(AANE)獲得淺層特征,并利用堆疊的自動編碼器(SAE)獲得深層特征,最后利用XGboost分類器進行預測。一般來說利用信息指標進行鏈路預測只局限于部分結構,無法利用到全面的圖結構信息。單純利用傳統的機器學習模型雖然也能取得較好的訓練效果,但是在驗證中相對來說具有較高的假陽性率,不利于生物實驗的驗證。矩陣分解方法的結果假陽性率低,但是偏重于已有知識的強化,發現新知識的能力較弱。

為了能夠提高預測的能力,就需要引入更多生物知識及其關系網絡來提取特征信息,比如構建生物知識的異構網絡等[20]。隨著近年來圖表示學習(graph represent learning)算法的發展,圖表示學習在人類社會網絡鏈路預測的相關問題上取得了較好的結果[21]。一些圖表示學習方法被用于環狀RNA與疾病關聯關系的預測,如Lei通過隨機游走算法實現特征的提取,并利用K鄰接聚類的方法實現環狀RNA與疾病關聯關系的預測[22];本課題組發表的iGRLCDA通過因子圖卷積網絡(factor Graph Convolution Network)在異構圖上提取特征[23],利用隨機森林分類器實現環狀RNA與疾病關聯關系的預測,取得了較好的結果。理論上,圖卷積網絡也可以直接做鏈路預測[14],但是不容易訓練成功?？紤]到環狀RNA與疾病的關系中大部分關系未知,所以iGRLCDA利用因子圖卷積網絡在主要的圖結構上對節點分類(node classification)。依據節點分類模型提取出所有節點的特征,最后依據分類器實現鏈路關系預測。在iGRLCDA的設計過程中,發現對傳統機器學習方法進行調優的過程比較費時且需要一定技巧,希望設計一種具有自適應且綜合性能良好的模型來實現環狀RNA與疾病關聯關系的預測。深度學習模型無疑具有較好的自適應性,但目前對于環狀RNA與疾病關系預測深度學習AE-DNN模型[18]部分性能并不出眾,反映非平衡數據性能的MCC指標為0.58,低于iGRLCDA[23]的0.714 6。此外,在驗證集上AE-DNN模型的AUC為0.85,也低于iGRLCDA[23]的0.928 7。在實現自動編碼器(AutoEncoder,AE)與深度全連接神經網絡(Deep Neural Network,DNN)的基礎上,嵌入圖注意力機制(Graph Attention Network,GAT)[24],實現了GAT-AE-DNN結構的端到端的深度學習模型GATECDA,在環狀RNA與疾病預測的CircR2Disease數據集中[25],其綜合性能AUC得分為0.961 8,MCC關系為0.757 6。GATECDA采用端到端的GAT-AE-DNN深度學習模型,具有自適應性、易于泛化和拓展等特點,訓練過程也更容易。

2 基于圖表示學習方法的預測

基于圖表示學習方法進行特征提取并預測關聯關系的基礎在于從圖中學習相應的知識并將圖結構信息融合入圖中節點的特征。相較于傳統上只利用節點內部的信息,圖表示學習可以利用節點有聯系的不同節點的特征來強化自身以反映與相關節點的聯系。以環狀RNA參與的ceRNA調控網絡為例,如果只考慮其自身的序列信息,那么可能在表示中無法反映出環狀RNA通過吸附miRNA來調節LncRNA的關系。但利用圖表示學習方法提取特征后,所提取的特征來源于環狀RNA自身,但也能把現有的調控關系反映出來。

目前,主要的圖表示學習方法有矩陣分解的方法、隨機游走的方法、圖神經網絡的方法等。其中圖注意力網絡(Graph attention networks,GATs)也是圖神經網絡中一種主要的方法[21,24],在多個同質數據集上的鏈路預測中取得了較好的性能。研究中首先建立異構的環狀RNA與疾病關系的網絡。所謂異構是因為環狀RNA或疾病在各自向量空間內存在關系圖,如圖1所示,需要在不同向量空間表述的節點關系中挖掘關聯關系。比如關系圖G=(u,v),其中的u與v分別表示不同類型的節點,它們各自在自身的向量空間存在不同的維度u_feature和v_feature。已經知道部分u與v之間存在聯系,因此構成了異構關系圖。圖表示學習方法實質就是在考慮異構關系圖G的結構上把u_feature和v_feature映射到同一個空間成為node_feature,該node_feature可以區分整體關系圖G中不同節點的類別。

隨后,u與v之間的已知關系(u,v)->R為預測的正樣本集,隨機產生的關系(u,v)->R*為預測的負樣本,正負樣本具有相同的大小N(N=739)并一同作為大小為2N的訓練集。在訓練集上采用五折交叉驗證。此外,為了驗證不同模型的性能,從訓練集中拿出n(n=50)個關系作為驗證集。最后,將提取的節點特征聯系起來利用分類器進行預測。圖1展示了GATECDA的整體流程,從異構生物知識中獲得環狀RNA與疾病的特征,并用深度模型預測關聯關系。

3 實驗結果與分析

3.1 實驗環境及參數設置

研究工作在一臺雙路Intel至強2365V2處理器的工作站上實現,內存為96 GB。在實現過程中,實際使用內存不超過16 GB,主要在屬性節點的特征提前上花費較多。GATECDA模型采用python 3.7語言實現,模型利用tensorflow 2.7張量流計算框架和keras深度學習框架封裝構建,GAT層的實現采用了dgl圖神經網絡工具包。

3.2 數據集

考慮通過環狀RNA的序列信息相似性,疾病關系的語義信息相似性和由已知的環狀RNA與疾病關系信息相似性來建立異構網絡。其中,環狀RNA序列信息源自circBase[26]數據庫中基于hg19基因組的推測的環狀RNA選擇性剪切序列。疾病關系的語義信息采用引用字典Mesh的關系獲得[27]。環狀RNA與疾病關系信息由CircR2Disease數據庫[25]中經過實驗驗證的關系獲得。部分因數據庫環狀RNA的id對應不上的序列也可以由CircR2Disease數據庫[25]提供的基因組位置或對應的基因Symbol獲得。一共獲得739個環狀RNA與疾病關系作為正樣本集,涉及到661個環狀RNA和100種疾病。在這個關系中,還存在65 261個未標注的環狀RNA與疾病的隨機關系,隨機從里面取得739個作為負樣本集。最后從1 478個正負樣本關系中取出50個關系作為驗證集,剩余的1 428個關系作為訓練集。

3.3 環狀RNA與疾病的特征提取

根據獲得的數據信息,可以構建三組節點間相似關系信息,包括環狀RNA與疾病、環狀RNA與環狀RNA、疾病與疾病。

(1)環狀RNA與疾病關聯:所有從CircR2Disease[25]的739個環狀RNA與疾病關系,涉及到661個環狀RNA和100種疾病,可以構成661×100的關系矩陣RD,其中有關系為1,否則為0。從該關系矩陣就可以通過Gaussian Interaction Profile (GIP)方法獲得單個環狀RNA或疾病的特征向量。GIP方法也是藥物與疾病關系等預測中常使用的方法[28],可以通過函數SE(p(i),p(j))從關系矩陣中兩個表示為0-1向量V(p)獲得節點i與j的相似性,如公式(1)。

(1)

(2)

其中,V(p(i))-V(p(j))表示兩個0-1向量間的差異,通過L2范式獲得差異的距離,乘以歸一化因子θ后獲得e指數的冪。最后,通過冪指數函數SE可以獲得0-1關系矩陣RD中任意兩個節點間的相似性,進而原來稀疏的0-1關系矩陣就轉化為稠密關系。其中環狀RNA或疾病可以獲得761個維度的特征。

(2)環狀RNA與環狀RNA相似性:可以獲得環狀RNA的序列信息,并通過序列相似性獲得環狀RNA與環狀RNA的661×661的相似矩陣CC。由此,可生成單個環狀RNA的特征向量。這里的環狀RNA的相似性由skip-gram結構的word2vec生成[29]。因為RNA序列結構的復雜性,RNA序列的作用區域可能局限于內部的短序列片段中,直接獲取兩條RNA序列的相似性不能反映它們相互作用的關系[30]。word2vec模型在自然預言處理中廣泛使用,它通過一個單詞在上下文中的出現關系來挖掘其特征表示。在生物序列的挖掘中,定義6-mer,如“ACCATC”為一個單詞w。

(3)

word2vec在該任務中是尋找參數Θ使得所有屬于語料T中每個句子S內單詞W的聯合概率乘積最大。在訓練中語料T包括13 000條環狀RNA序列。

(4)

(3)疾病與疾病相似性: 建立疾病與疾病100×100的相似關系,就可以獲得單個疾病100個維度的特征信息。疾病與疾病的相似關系源自MeSH數據庫。作為醫學引用詞典,MeSH數據庫通過分析大量醫學論文的引用關系提供了醫學主題詞關系。利用醫學主題詞關系,基于王等人[31]發表的方法,可以獲得關于疾病間的相似關系。醫學主題詞關系構建了有向無環圖(DAG)?？梢杂浤骋患膊參與的DAG(d)=(d,N(d),E(d)),其中N(d)表示與某一疾病相關的所有節點,包括疾病或者癥狀;E(d)表示與之涉及的所有邊。如果在DAG(d)中還存在另一疾病s,那么可以通過如下公式計算疾病d與疾病s的關系:

(5)

在公式(5)中如果疾病d與疾病s關聯,那么它們的關系為1,否則找出疾病d到疾病s所有共同關聯的子節點數量,作為它們之間的關系。在復雜疾病中,疾病d的影響力為所有與之有關疾病的關系的累加和,定義如下:

(6)

有了以上(6)的信息,可以定義兩個疾病間的互信息SS1:

在公式(7)中,兩兩疾病間的相似關系可以理解為與它們相關所有節點的關系除以兩個疾病的整體影響。但是有些疾病可能影響的節點少,但它卻很重要,于是設計了另一個指標DCd(s):

(8)

其中,num(contain(DAG(d),s))表示DAG(d)圖中包含疾病s的數量,num(diseases)表示所有的疾病。這樣關聯數量少的疾病DC的分就越高。于是,有了第二個衡量疾病關系的互信息SS2:

(9)

最后,將SS1與SS2共同考慮得到SS=0.5*SS1+0.5*SS2,作為最后疾病之間的語義相似關系。

3.4 GATECDA模型的實現

在GATECDA的實現如圖2所示。首先,構建了環狀RNA與疾病的初始特征,計算環狀RNA與疾病之間關聯關系的相似性,疾病的語義相似性和環狀RNA的序列相似性。其次,GATECDA加入了圖注意力網絡(Graph attention networks,GATs)提取環狀RNA與疾病異質關系圖中的特征表示。最后,將得到的環狀RNA與疾病的特征表示送入AE-DNN深度學習模型進行關系預測,其中包含了自動編碼器(AutoEncoder,AE)和深度神經網絡(Deep Neural Network,DNN)。筆者認為GAT起到了特征提取與融合的作用,AE起到了特征降維的作用,DNN起到了分類器的作用。單層圖注意力網絡GAT也是由數個神經元組成的單元,一般不超過三層,比圖卷積網絡更容易達到訓練效果[24]。相比圖卷積網絡是一種淺層的神經網絡結構,因為本身屬于神經網絡,所以可以嵌入到深度學習模型中。

圖2 GATECDA模型深度學習模型的結構

模型首先接受生物知識圖G及其節點特征。圖G可以認為是一個M*N的二部圖(bipartie graph)。M可以認為是所有的環狀RNA,而N為疾病,同時M和N各自的特征也被作為參數。圖注意力網絡在接受數據后完成了以下工作:

Wupdatenode=[sigmoid(X*[F(j),F(i)])]

(10)

(11)

(12)

F*(i)=LeakyReLU(α*F(i))

(13)

其中,j表示i節點的所有鄰接節點。Wupdatenode構成了輸入層的神經網絡,X*[F(j),F(i)]為該層輸入的數據,其中X為自定義特征矩陣,[F(j),F(i)]表示i和j的聯合特征向量。在學習一遍所有節點后,希望單個節點更新后在整體中起到最大作用,這里用α體現特征的更新,F*是更新后的特征。此外,作為一種隨機過程,每更新一輪被認為是1個頭(head)的注意力,更新k次為多個頭(k-heads)的注意力,在GATECDA中k為8。最后,所有1至k次的特征更新都被均方和作為最后的特征,如公式(14):

(14)

注意力的思想與word2vec一致,就是每個節點都朝著在整體背景中最顯著去改變。而多頭的概念與主成分分析(PCA)的概念相似。所以認為多頭注意力網絡起到了特征提取與融合的作用。隨后的AE-DNN模型由自動編碼器(AutoEncoder,AE)和深度神經網絡DNN(Deep Neural Network)組成,是深度學習中的經典模型,在很多機器翻譯任務中都有較為出色的表現。AE層接受稀疏的數據,在不斷收窄的多層網絡中實現信息的融合、壓縮與標準化,之后又以多層變寬的網絡壓縮后的數據還原回輸入數據。AE具有降維的作用,在GATECDA中,如圖2(2)把兩層GAT網絡得到的1 522維的特征壓縮為128維的特征。經過AE處理過的數據又被送入深度神經網絡6層神經網絡構建的DNN進行關聯關系的分類預測,如圖2(3)。在所有的AE-DNN層中,都使用了Batch normalization和dropout機制。Batch normalization是一種歸一化方法,可以減小異常數據的干擾。dropout機制是在每一層反饋梯度時,只更新一定比例的神經元,該模型訓練時dropout的值為0.3。Batch normalization和dropout機制都是為了防止模型過擬合,提高模型泛化能力。

3.5 評估指標

在取得對預測結果評估矩陣的真陽性率(True Positive,TP)、真陰性率(True Negative,TN)、假陽性率(False Positive,FP)、假陰性率(False Negative,FN)后,采用了準確率(Acc.)、敏感度(Sen.)、精準率(Pre.)、F1打分(F1)和Matthews關系(MCC)來較全面地評估模型的性能,這些也是機器學習領域的主流評價方法,如下:

(15)

(16)

(17)

(18)

MCC=

(19)

在五折交叉驗證的測試下,衡量受試者工作特征曲線(ROC)下面積(AUC)也是機器學習領域里衡量模型性能的主要指標。通過模型在逐一增長的測試集上預測結果真陽性率(TPR)與假陽性率(FPR)的平面坐標位置,就可以做出ROC曲線。

3.6 模型能力評估

為了評估GATECDA模型的能力,在CircR2Disease數據集上進行五折交叉驗證,即將訓練集劃分為5等份,進行五次訓練。每次以其中四份進行訓練,一份進行測試(285個樣本)。圖3展示了GATECDA模型的訓練過程的ROC曲線及AUC值。GATECDA模型的五折交叉驗證平均AUC值為0.961 8,每次的AUC值分別為0.947 6,0.952 0,0.963 7和0.979 5。其綜合性能在表1中體現,平均準確率為87.53%,敏感度為93.62%,精準度為83.80%,F1打分為88.35%,MCC關系為0.757 6, 精準度-召回曲線下面積AUPRC為0.903 2,ROC曲線下面積AUC為0.961 8。

表1 GATECDA在CircR2Disease數據集上五折交叉驗證

圖3 GATECDA模型在CircR2Disease 數據集生成的ROC曲線

3.7 不同預測模型比較

比較了已經發表的環狀RNA與疾病關聯關系預測的幾種方法在CircR2Disease[25]數據集上五折交叉驗證中的AUC值, 見表2。 它們包括基于圖表示學習方法GATECDA、iGRLCD[23]和GIS-CDA[12],深度學習模型AE-DNN[18]與AANE-SAE[19],以上模型在文中研究背景中均有介紹。通過比較可以看出,GATE-CDA在五折交叉驗證中平均的AUC為0.961 8,高于iGRLCDA[23]的0.928 7和AE-DNN[18]的0.930 3。對于衡量不平衡數據集上性能的MCC值,GATECDA的0.757 6,也高于AE-DNN的0.583 6和iGRLCDA模型的0.714 6。其中GIS-CDA與GATECDA模型都采用了圖注意力機制,不過GIS-CDA是先用編碼器融合不同維度的特征后再使用圖注意力機制, GATECDA模型首先使用圖注意力機制而不是進行編碼的信息融合,因而比GIS-CDA模型的AUC略高。筆者認為在設計異構網絡模型時,越能完整和直接地利用圖結構信息,越有利于模型的預測。GATECDA不足在于實現的圖注意力機制(CAT)是一種淺學習[14,24],對于以后更大規模數據集或知識圖譜上能力提升空間不如圖卷積網絡(GCN)模型[21]

表2 不同預測模型的比較

3.8 不同分類器比較

比較GATECDA和不同分類器模型在驗證集上的預測能力。其中KNN、RF、XGboost和SVM為iLearnPlus工具[32]封裝好的分類器。GATECDA是該文提出的端到端圖注意力網絡、自動編碼器與深度神經網絡結合的深度學習模型(GAT_AE_DNN),其中AE是自動編碼器加輸出層的分類器,DNN是深度神經網絡分類器。SVM是支持向量機(Support Vector Machine),KNN是K鄰接分類器(K-nearest Neighbor),RF是隨機森林分類器(Random Forest),XGboost是極限學習分類器(Extreme Gradient boost)。以上所有模型都在1 428個正負關系構成的訓練集上加以訓練,并在獨立劃分出的50個樣本的驗證集上做性能比較。從圖4中可以看出,在驗證集樣本上GATECDA的AUC最高為0.972 6, XGboost的AUC值為0.895 0,KNN為0.733 3,RF為0.640 8, SVM為0.667 2。

圖4 不同分類器模型在驗證集上的ROC曲線

3.9 特征消融實驗

為了分析圖結構的已有知識信息與節點屬性信息對模型能力的貢獻,設計了特征消融實驗,見表3。研究中,GATECDA模型既使用已有知識構建圖G,也采用節點屬性特征,得到的預測結果AUC為0.961 8,AUPR為0.903 2。GATECDA-F是GATECDA模型只包含圖結構信息,得到的預測結果AUC為0.582 7,AUPR為0.785 7。GATECDA-G是GATECDA模型只包含節點屬性特征,得到的預測結果AUC為0.491 5,AUPR為0.732 8。最后為該結果符合預期,圖注意力網絡在考慮圖結構和節點屬性特征時可以強化特征信息。

表3 特征消融實驗

4 案例研究

通過GATECDA從661個環狀RNA和100種疾病的65 261個未標注潛在組合中預測3 743個關聯關系,約占未標注總數的5.7%。表4列出預測結果排名前30的關聯關系,并且通過文獻檢索查到相關CircRNA或其所在基因在以前的生物實驗中有發現與相關疾病存在聯系。在預測的結果得到的前30個環狀RNA與疾病的關聯關系中,其中有25個關聯能夠在最近醫學文獻中被發現存在關聯。預測結果可以幫助研究人員縮小篩查范圍,盡快找到與疾病相關的關鍵標志物。實驗中獲得的差異信息很多,一般的方法是做富集分析或是在基因共表達網絡尋找關鍵基因。如果結合已有知識對環狀RNA與疾病的關聯關系預測可以為尋找關鍵基因和疾病標志物提供一種新的角度。

表4 預測排名前30個環狀RNA與疾病的關系及文獻檢索

續表4

5 挑戰與發展

筆者認為,目前采用圖表示學習提取特征進行環狀RNA與疾病關聯關系預測的方法比其他方法能獲得較好的綜合性能。針對目前取得的進展,一方面需要利用更豐富的生物網絡知識,即利用復雜異構網絡實現對任意環狀RNA與疾病的預測,同時保持驗證中較低的假陽性率。從這一點上看,GATECDA的基礎在于已有知識的補全,因而更適合于降低假陽性率的新知識的挖掘。另一方面,研究環狀RNA與疾病關系的預測最初也是想實現環狀RNA、調控分子、生物過程、生物性狀到疾病完整鏈路的預測,但相關的知識和計算方法以前達不到一定的積累。隨著圖神經網絡、圖表示學習和生物信息等方法在相關方面的進展,關聯關系預測方法與生物知識的不斷積累,圖表示學習的方法能夠在與大規模知識圖譜不斷結合與發展。利用GATECDA多頭注意力機制和易于訓練的特點,在多目標的二部圖(bipartite graph)結構中應當會比較適用,挖掘出中間的調控過程,實現鏈路預測。

6 結束語

環狀RNA與疾病關聯關系預測模型在利用圖表示學習機制后性能有所提升,圖神經網絡與深度學習結合的模型更易于訓練與泛化。筆者認為利用人工智能技術挖掘已有生命科學知識進行相關的預測,其結果可以有助于解釋在高通量實驗中發現的大量異常信息,為研究人員推薦出與研究背景相關的關鍵信息,這將加快和提高相關領域的研究進展。