本發(fā)明屬于系統(tǒng)生物學領域，涉及一種基于雙向隨機游走和多標簽學習的miRNA-環(huán)境因子關系預測方法。

背景技術：

眾所周知，疾病和遺傳因子和環(huán)境因子有著密切關系。一些疾病的產(chǎn)生可視為遺傳因子和環(huán)境因子之間的擾動造成的。因此，識別遺傳因子和環(huán)境因子之間關系有助于我們解析疾病復雜的機制，且進一步有助于疾病診斷，預防及其治療。

MicroRNA(miRNA)是一類新發(fā)現(xiàn)的、重要的遺傳因子。一般來說，miRNA可以通過和信使RNA的3’端非編碼區(qū)域的完全配對或者5’端非編碼區(qū)域的不完全配對，抑制信使RNA的表達，從而達到調控基因表達的目的。miRNA產(chǎn)生過程包括以下幾個過程：1)通過核糖核酸聚合酶II(RNA polymerase II)對基因的內(nèi)含子特定位點的剪切，形成初始miRNA(primary miRNA)。這種初始miRNA的長度一般為100-1000nt核苷酸。并且在5’端末尾有一個cap結構，在3’端末尾包含ploy(A)特征。2)初始miRNA經(jīng)過RNase III酶Drosha的剪切形成前體miRNA(precursor miRNA)。前體miRNA的長度一般為70nt核苷酸長度，有發(fā)夾結構，且在3’端末尾有2個突出的核苷酸。3)前體miRNA經(jīng)過RNase III酶Dicer剪切最終形成成熟體miRNA。成熟體miRNA的長度一般為19-22nt，且具有高度保守型。許多研究表明：miRNA在人類的很多重要的生命過程中起著舉足輕重的作用，包括細胞的生長、發(fā)育、增殖、分化、凋亡等。

近年來，研究發(fā)現(xiàn)：許多環(huán)境因子(EF)和miRNA的表達有著密切關系,環(huán)境因子能夠調控miRNA的表達。例如，在HepG2細胞系中，鞣花單寧(日本蛇菰屬植物提煉出來的多酚化合物)能夠調節(jié)25個miRNA表達。其中，17個miRNA為上調，8個miRNA為下調。這些環(huán)境因子和miRNA關系的擾動經(jīng)常和疾病的產(chǎn)生有著密切關系。在前列腺癌細胞中，研究人員發(fā)現(xiàn)染料木黃酮能夠促進miR-34a的表達，進而影響基因HOTAIR的表達。而該基因和癌細胞的增殖有著密切關系。因此，研究環(huán)境因子和miRNA之間的關系能夠有助于人類了解疾病產(chǎn)生機理，進而幫助人類設計有效的疾病治療方案。

隨著高通量技術不斷發(fā)展，產(chǎn)生了大量的miRNA相關數(shù)據(jù)。一些miRNA相關的數(shù)據(jù)庫也建立來存儲這些miRNA數(shù)據(jù)。miRbase數(shù)據(jù)庫是目前最權威的miRNA序列數(shù)據(jù)庫，該數(shù)據(jù)主要存儲miRNA的序列、前體結構及其miRNA對應的靶基因數(shù)據(jù)。為了保證數(shù)據(jù)的準確性，該數(shù)據(jù)庫只收集已經(jīng)被驗證過的miRNA相關的數(shù)據(jù)。miRecord數(shù)據(jù)庫為動物的miRNA靶基因數(shù)據(jù)庫，包括生物實驗驗證的靶基因和計算方法預測的靶基因。DbDEMC是一個專門存儲人類癌癥差異表達的miRNA數(shù)據(jù)，該數(shù)據(jù)收錄14種癌癥中607個差異表達的miRNA。miREnvironment數(shù)據(jù)是專門收錄miRNA和環(huán)境因子之間的關系的數(shù)據(jù)庫。該數(shù)據(jù)庫手動收集了24個物種，1242個miRNA，394個環(huán)境因子，3857條關系。

目前，在生物學中，傳統(tǒng)的預測miRNA和環(huán)境因子之間的關系的實驗方法主要包括PCR和微陣列方法。生物學家在特定組織細胞中，加入特定環(huán)境因子，進而觀測miRNA表達的差異。根據(jù)這些方法能夠準確有效的發(fā)現(xiàn)miRNA和環(huán)境因子之間關系。但是，這些方法非常費時且成本很高，尤其是在大量的生物數(shù)據(jù)面前，更顯得捉襟見肘。大規(guī)模預測miRNA和環(huán)境因子之間的關系的計算方法應運而生，利用計算方法能夠有效降低經(jīng)濟花費，提高實驗效率。利用計算方法來預測miRNA和環(huán)境因子之間的關系可以簡單分為以下幾類：

1)基于機器學習的方法

由于miRNA和環(huán)境因子的數(shù)據(jù)種類越來越多，基于機器學習模型可以利用這些數(shù)據(jù)，研究miRNA-EF關聯(lián)預測方法?；谧钚《朔ǎ珻hen等提出一種新的miRNA和環(huán)境因子之間的關系預測方法。在該方法中，采用兩個分類器分別對miRNA和環(huán)境因子進行最優(yōu)化。最終，通過線性加權方法得到每對miRNA和環(huán)境因子的打分值。

2)基于網(wǎng)絡的方法

網(wǎng)絡的方法都是基于相同的假設：相似的miRNA往往和相似的環(huán)境因子有關系。Chen等人基于隨機游走方法提出一種新的miRNA和環(huán)境因子之間的關系預測方法；該方法融合了miRNA功能相似性，miRNA拓撲相似性，EF結構相似性以及EF拓撲相似性。另外，Li等提出一種基于網(wǎng)絡推斷方法來預測miRNA-環(huán)境因子的關聯(lián)關系，該方法利用已知的miRNA和環(huán)境因子之間的關系數(shù)據(jù)，構建miRNA-環(huán)境因子二分圖。在構建的二分圖上，采用網(wǎng)絡推斷方法來預測新的miRNA和環(huán)境因子之間的關系。該方法所基于的假設為：如果兩個miRNA共享多個環(huán)境因子，那么和其中一個miRNA有關系的環(huán)境因子也很有可能和另外一個miRNA有關系。

3)基于統(tǒng)計的方法

Qiu等人提出一種miRNA-環(huán)境因子關系預測方法。該方法融合miRNA組織差異表達數(shù)據(jù)和miRNA-疾病關系數(shù)據(jù)。采用富集分析方法來識別潛在的miRNA-環(huán)境因子之間關系。

上述各類方法從不同的角度來預測潛在的miRNA-環(huán)境因子關聯(lián)。雖然上述方法取得不少成就，但是也存在一些缺陷。首先，計算miRNA和環(huán)境因子的相似性上，有的方法采用網(wǎng)絡的拓撲特征來計算miRNA和環(huán)境因子的相似性。忽略了miRNA和環(huán)境因子的生物屬性。其次，在融合不同miRNA和環(huán)境因子數(shù)據(jù)上，有的方法采用算術均值方法來融合不同的miRNA和環(huán)境因子相似性。這可能會引入其他的噪聲。另外，有的方法無法對一個新的miRNA和環(huán)境因子之間的關系進行預測。因此，有必要設計一種新的預測miRNA-環(huán)境因子關系預測方法。

技術實現(xiàn)要素：

本發(fā)明所要解決的技術問題是，針對現(xiàn)有技術的不足，提供一種基于雙向隨機游走和多標簽學習的miRNA-環(huán)境因子關系預測方法(MEI-BRWMLL)，該方法易于實施，預測準確率高。

本發(fā)明為解決技術問題所采用的技術方案如下：

一種基于雙向隨機游走和多標簽學習的miRNA-環(huán)境因子關系預測方法，包括以下步驟：

1)首先基于已知的miRNA-環(huán)境因子關系，構建miRNA-環(huán)境因子矩陣ME；然后計算miRNA相似性和環(huán)境因子相似性，構建miRNA相似性矩陣和環(huán)境因子相似性矩陣；

2)根據(jù)miRNA相似性矩陣和環(huán)境因子相似性矩陣分別構建miRNA相似性網(wǎng)絡和環(huán)境因子相似性網(wǎng)絡；再基于miRNA相似性網(wǎng)絡、環(huán)境因子相似性網(wǎng)絡和miRNA-環(huán)境因子矩陣，構建miRNA-環(huán)境因子關系網(wǎng)絡；

3)采用雙向隨機游走(Bi-random walk,BiRW)和多標簽學習方法(Multi-label learning,MLL)來預測潛在的miRNA-環(huán)境因子關系：

a)對于miRNA和環(huán)境因子都是已知的情況，采用雙向隨機游走方法，分別在miRNA相似性網(wǎng)絡和環(huán)境因子相似網(wǎng)絡上進行不同步數(shù)游走，得到miRNA-環(huán)境因子得分矩陣，分值越大表明對應的miRNA和環(huán)境因子存在關系的可能性越大；

b)對于新的miRNA，基于miRNA的相似性網(wǎng)絡和環(huán)境因子的相似性網(wǎng)絡，采用多標簽學習方法，得到新的miRNA與環(huán)境因子之間存在關系的概率值，概率值越大，兩者存在關系的可能性越大；

c)對于新的環(huán)境因子，基于miRNA的相似性網(wǎng)絡和環(huán)境因子的相似性網(wǎng)絡，采用多標簽學習方法，得到新的環(huán)境因子與miRNA之間存在關系的概率值，概率值越大，兩者存在關系的可能性越大。

所述步驟1)中，構建miRNA-環(huán)境因子矩陣ME，miRNA-環(huán)境因子矩陣每一行對應一個miRNA，每一列對應一個環(huán)境因子；若已知miRNA m_i和環(huán)境因子e_j存在關系，則ME(i,j)等于1；否則，ME(i,j)等于0；其中i＝1,2,…,m；j＝1,2,…,e；m和e分別為已知的miRNA和環(huán)境因子個數(shù)。

所述步驟1)中，構建miRNA相似性矩陣的方法為：

首先，基于miRNA序列信息和miRNA-環(huán)境因子矩陣ME，分別計算miRNA序列相似性和miRNA交互譜相似性；

然后，基于miRNA序列相似性和miRNA交互譜相似性，構建兩個的miRNA相似性矩陣；

最后，采用相似性矩陣融合方法分別將不同的miRNA相似性矩陣融合為一個miRNA相似性矩陣。

所述步驟1)中，構建環(huán)境因子相似性矩陣的方法為：

首先，基于環(huán)境因子的化學結構、化學分類系統(tǒng)信息和miRNA-環(huán)境因子矩陣ME，分別計算環(huán)境因子化學結構相似性、化學分類系統(tǒng)相似性和環(huán)境因子交互譜相似性；

然后，基于環(huán)境因子化學結構相似性、化學分類系統(tǒng)相似性和環(huán)境因子交互譜相似性構建三個環(huán)境因子相似性矩陣；

最后，采用相似性矩陣融合方法分別將不同的環(huán)境因子相似性矩陣融合為一個環(huán)境因子相似性矩陣。

所述步驟2)包括以下步驟：

首先，基于步驟1)中融合得到的miRNA相似性矩陣和環(huán)境因子相似性矩陣，分別構建miRNA相似性網(wǎng)絡和環(huán)境因子相似性網(wǎng)絡；在miRNA相似網(wǎng)絡中，節(jié)點為miRNA，邊的權值為miRNA-miRNA相似性值；在環(huán)境因子相似性網(wǎng)絡中，節(jié)點為環(huán)境因子，邊的權值為環(huán)境因子-環(huán)境因子相似性值；

然后，根據(jù)已知的miRNA-環(huán)境因子關系，構建miRNA-環(huán)境因子關聯(lián)二分圖；在miRNA-環(huán)境因子關聯(lián)二分圖中，節(jié)點分別為miRNA和環(huán)境因子，若miRNA m_i和環(huán)境因子e_j存在關系，則miRNA m_i和環(huán)境因子e_j之間存在一條邊，否則沒有邊；

最后，融合miRNA相似性網(wǎng)絡、環(huán)境因子相似性網(wǎng)絡和miRNA-環(huán)境因子關聯(lián)二分圖，構成miRNA-環(huán)境因子關系網(wǎng)絡。

所述步驟a)具體為：

a1)對融合后的miRNA相似性矩陣和融合后的環(huán)境因子相似性矩陣進行標準化，公式為：

RM＝DM^-1/2*FM*DM^-1/2

RE＝DE^-1/2*FE*DE^-1/2

其中，RM表示標準化的miRNA相似性矩陣，DM表示融合后的miRNA相似性矩陣的對角矩陣，DM對角線的值為對應的FM中一行的值的和；

RE表示標準化的環(huán)境因子相似性矩陣，DE表示融合后的環(huán)境因子相似性矩陣的對角矩陣，DE對角線的值為對應的FE中一行的值的和；

對miRNA-環(huán)境因子矩陣ME進行標準化，公式為：

A(i,j)＝ME(i,j)/sum(ME)

其中，A(i,j)表示標準化的miRNA-環(huán)境因子矩陣A中第i行j列的元素；ME(i,j)代表miRNA-環(huán)境因子矩陣ME中第i行j列的元素，sum(ME)為ME中所有元素之和；

a2)進行以下迭代：

在miRNA相似性網(wǎng)絡游走：

RD_L(t)＝a×RM×RD(t-1)+(1-a)A

其中，a為固定參數(shù)【a的取值通過交叉驗證確定，本發(fā)明實施例中，a的取值通過交叉驗證最終確定為0.8】，RD(t-1)為游走t-1步后預測的miRNA-環(huán)境因子得分矩陣，RD_L(t)為在miRNA相似性網(wǎng)絡上游走t步后預測的miRNA-環(huán)境因子得分矩陣；RD(0)＝A；

在環(huán)境因子相似性網(wǎng)絡游走：

RD_R(t)＝a×RD(t-1)×RE+(1-a)A

其中，RD_R(t)為在環(huán)境因子相似性網(wǎng)絡上游走t步后預測的miRNA-環(huán)境因子的得分矩陣；

最終輸出游走t步后預測的miRNA-環(huán)境因子得分矩陣RD(t)：

設定隨機游走在miRNA相似性網(wǎng)絡和環(huán)境因子相似性網(wǎng)絡中的最大迭代步數(shù)分別為l和r；【l和r的取值通過交叉驗證確定，本發(fā)明實施例中，l和r的取值通過交叉驗證最終確定分別為4和2】當隨機游走在兩個網(wǎng)絡中的步數(shù)超過最大迭代參數(shù)l或r時，迭代過程終止；

a3)對于任意兩個已知的miRNA和環(huán)境因子，根據(jù)RD(t)的值判斷兩者存在關系的可能性。

所述步驟b)具體為：

b1)按如下方法計算新的miRNA m_c與已知的環(huán)境因子e_j存在關系的概率P(m_c,e_j)：

其中，和表示兩個先驗概率；和表示兩個后驗概率；

其中，e(i)表示miRNA m_c的K個最近鄰居中，本身和e_j有關系，且其對應的K個最近鄰居中，有i個miRNA和環(huán)境因子e_j有關系的miRNA的個數(shù)；

e′(i)表示miRNA m_c的K個最近鄰居中，本身和e_j有關系，且其對應的K個最近鄰居中，有i個miRNA和環(huán)境因子e_j沒有關系的miRNA的個數(shù)；

s的取值采用交叉驗證得到；

b2)根據(jù)P(m_c,e_j)的值判斷新的miRNA m_c與已知的環(huán)境因子e_j存在關系的可能性。

所述步驟c)具體為：

c1)按如下方法計算新的環(huán)境因子e_c與已知的miRNA m_i存在關系的概率P(e_c,m_i)：

其中，和表示兩個先驗概率；和表示兩個后驗概率；

其中，m(j)表示環(huán)境因子e_c的K個最近鄰居中，本身和m_i有關系，且其對應的K個最近鄰居中，有j個環(huán)境因子和miRNA m_i有關系的環(huán)境因子的個數(shù)；

m′(j)表示環(huán)境因子e_c的K個最近鄰居中，本身和m_j有關系，且其對應的K個最近鄰居中，有j個環(huán)境因子和miRNA m_i沒有關系的環(huán)境因子的個數(shù)；

r的取值采用交叉驗證得到；

c2)根據(jù)P(e_c,m_i)的值判斷新的環(huán)境因子e_c與已知的miRNA m_i存在關系的可能性。

所述s的取值為10，r的取值為5。

有益效果：

本發(fā)明通過集成雙向隨機游走和多標簽學習算法，來預測環(huán)境因子潛在的miRNA。首先，利用已知的不同生物數(shù)據(jù)信息，計算得到不同的miRNA相似性矩陣和環(huán)境因子相似性矩陣。然后，基于計算得到不同的miRNA相似性和環(huán)境因子相似性，分別創(chuàng)建不同的miRNA相似性矩陣和不同的環(huán)境因子相似性矩陣，再用相似性矩陣融合方法分別對不同的miRNA相似性矩陣和不同的環(huán)境因子相似性矩陣進行融合，減低單個數(shù)據(jù)源的得到的相似性噪聲，提高最終miRNA和環(huán)境因子相似性的可靠性。最后對不同情況下miRNA和環(huán)境因子，分別采用雙向隨機游走和多標簽學習來預測潛在的miRNA-環(huán)境因子關系，該方法能有效地挖掘環(huán)境因子潛在的miRNA。

附圖說明

圖1：本發(fā)明MEI-BRWMLL流程圖；

圖2：環(huán)境因子的度分布圖；

圖3：ClusterViz在已知數(shù)據(jù)集上識別出來的三個模塊；

圖4：在數(shù)據(jù)集上的十倍交叉驗證。

具體實施方式

以下將結合附圖和具體實施例對本發(fā)明做進一步詳細說明：

如圖1所示，本發(fā)明具體實現(xiàn)過程如下：

一.miRNA相似性計算和環(huán)境因子相似性計算

1.miRNA-環(huán)境因子矩陣構建

基于已知的miRNA-環(huán)境因子關系，本發(fā)明首先構建miRNA-環(huán)境因子矩陣ME。其中，miRNA-環(huán)境因子矩陣每一行對應一個miRNA，每一列對應一個環(huán)境因子。若miRNA m_i和環(huán)境因子e_j存在關聯(lián)關系，ME(i,j)等于1；否則，ME(i,j)等于0。

2.miRNA相似性計算

本發(fā)明miRNA相似性計算包括兩個部分：miRNA序列相似性計算和miRNA交互譜相似性計算。本發(fā)明采用Emboss-Needle工具來計算兩條成熟體miRNA序列的相似性。Emboss-Needle參數(shù)則是按照指定的參數(shù)(Matrix＝EDNAfull,Gap open＝10,Gap extend＝0.5)。

miRNA交互譜相似性計算是基于假設：相似miRNA一般都會有相似的交互模式，本發(fā)明定義兩個miRNA m_i和m_j之間的交互譜相似性為：

miR_Gip(m_i,m_j)＝exp(-γ_m||IP(m_i)-IP(m_j)||²)

其中，m_i和m_j分別代表第i個miRNA和第j個miRNA；i，j＝1,2,…,m；m表示已知的miRNA的總個數(shù)；IP(m_i)代表miRNA-環(huán)境因子矩陣中的第i行。

3.環(huán)境因子相似性計算

環(huán)境因子可以分為兩種：化學分子環(huán)境因子和非化學分子環(huán)境因子。對于化學分子環(huán)境因子，本發(fā)明采用化學結構相似性、化學分類系統(tǒng)信息相似性和交互譜相似性來計算其相似性。對于非化學分子環(huán)境因子，只采用交互譜相似性來計算相似性?；瘜W結構相似性計算是利用SIMCOMP工具來實現(xiàn)的，該工具是通過比對兩個化學分子的共同子結構大小來比較兩個化學分子的化學結構相似性，即共同的子結構越大，兩個化學分子的化學結構相似性越高?；瘜W分類系統(tǒng)信息(ATC：anatomical therapeutic chemical)是基于其在各個組織中化學，藥理及其作用等屬性，其提供為非結構信息。本發(fā)明采用Lin相似性計算方法來計算兩個化學分子之間的化學分類系統(tǒng)信息相似性，該方法基于節(jié)點的層次結構，化學分子的注釋信息重疊越多，相似性越高。

對于所有的環(huán)境因子(包括化學分子環(huán)境因子和非化學分子環(huán)境因子)，本發(fā)明采用以下公式計算兩個環(huán)境因子e_i和e_j之間的交互譜相似性：

EF_Gip(e_i,e_j)＝exp(-γ_e||IP(e_i)-IP(e_j)||²)

其中，e_i和e_j分別代表第i個環(huán)境因子和第j個環(huán)境因子；i，j＝1,2,…,e；e表示已知的環(huán)境因子的總個數(shù)，IP(e_i)代表miRNA-環(huán)境因子矩陣中的第i列。

二.相似性矩陣融合

從上面部分可知，本發(fā)明已經(jīng)計算得到2種miRNA的相似性數(shù)據(jù)和3種環(huán)境因子相似性數(shù)據(jù)。以miRNA為例子，首先，本發(fā)明對miRNA的兩種相似性矩陣按照下面方式分別進行標準化：

其中，W表示miRNA相似性矩陣，W(i,j)為其第i行第j列的元素；NM^v₀(i,j)表示經(jīng)過標準化后的miRNA相似性矩陣NM^v₀第i行第j列的元素。

為了度量相似性矩陣的局部一致性，本發(fā)明基于前面計算得到的miRNA相似性，對每個miRNA選擇相似性最大的前K個miRNA作為其最近鄰居集合，從而構建局部相似性矩陣：

其中，N_i表示第i個miRNA的最近鄰居集合。

對于miRNA每種數(shù)據(jù)類型，相似性矩陣融合過程是通過迭代更新相似性矩陣：

其中，SM^v表示數(shù)據(jù)類型為v的miRNA的局部相似性網(wǎng)絡；本發(fā)明中v＝1,2分別表示miRNA序列局部相似性矩陣和miRNA交互譜局部相似性矩陣；n代表是數(shù)據(jù)類型的總的個數(shù)，在這里，n為2；NM^v_u表示第u(u＝1,2,3…)次迭代更新的相似性矩陣；本發(fā)明中k＝1,2分別指代miRNA序列相似性矩陣和miRNA交互譜相似性矩陣；迭代計算NM^v_u，直到兩次迭代之間的差值NM^v_u-NM^v_u-1＜ε，結束迭代。本發(fā)明設定ε＝10e^-6。

對于每一種數(shù)據(jù)類型的miRNA相似性矩陣，都可以得到一個融合后的相似性矩陣。本發(fā)明通過算術平均值計算所有數(shù)據(jù)類型融合后的miRNA相似性矩陣FM，其定義如下：

對于環(huán)境因子的三種相似性矩陣，利用相同的相似性矩陣融合方法也可以得到一個融合后的環(huán)境因子相似性矩陣FE。

三.基于雙向隨機游走和多標簽學習方法預測miRNA-環(huán)境因子相互關系

本發(fā)明采用兩種方法來預測miRNA-環(huán)境因子之間關聯(lián)關系：雙向隨機游走來預測已知的miRNA和環(huán)境因子之間潛在的關聯(lián)關系和多標簽學習方法來預測新的miRNA(環(huán)境因子)和環(huán)境因子(miRNA)之間的關聯(lián)關系。本發(fā)明采用這兩種方法的理由包括以下幾點：首先，以前的研究表明雙向隨機游走在預測已知的樣本之間相互作用關系中取得很好的效果。但是其對于新的樣本之間關聯(lián)關系效果不是特別突出。然而，多標簽學習方法在新的樣本之間關系預測上優(yōu)勢很大。因此，本發(fā)明通過融合兩種方法優(yōu)勢來提高miRNA-環(huán)境因子之間關系。

(1)基于雙向隨機游走方法預測已知的miRNA-環(huán)境因子之間潛在相互關系

本發(fā)明預測miRNA-環(huán)境因子之間相互關聯(lián)關系基于的假設為：相似的miRNA一般都會和相似的環(huán)境因子之間有關系。其具體的實施步驟如下：

首先，本發(fā)明采用拉普拉斯標準化方法分別對上部分計算得到的融合后的miRNA相似性矩陣和融合后的環(huán)境因子相似性矩陣進行標準化。以miRNA為例，其標準化的定義如下：

RM＝DM^-1/2*FM*DM^-1/2

其中，DM表示融合后的miRNA相似性矩陣的對角矩陣。DM對角線的值為對應的FM中一行的值的和。通過類似的過程也可以得到標準化的環(huán)境因子相似性矩陣。

對于miRNA-環(huán)境因子矩陣ME，其標準化矩陣A定義為：

A(i,j)＝ME(i,j)/sum(ME)

其中，ME(i,j)代表miRNA-環(huán)境因子矩陣ME中第i行j列的元素，sum(ME)為ME中所有元素之和。

至此，本發(fā)明得到標準化的miRNA相似性矩陣(RM)，環(huán)境因子相似性矩陣(RE)和miRNA-環(huán)境因子關聯(lián)矩陣(A)。然后，本次采用雙向隨機游走方法來預測已知的miRNA和環(huán)境因子之間潛在的關聯(lián)關系?；趍iRNA相似性矩陣和環(huán)境因子相似性矩陣，構建miRNA相似性網(wǎng)絡和環(huán)境因子相似性網(wǎng)絡。其中，兩個網(wǎng)絡中，節(jié)點分別為miRNA和環(huán)境因子，邊的權重分別為miRNA-miRNA相似性值和環(huán)境因子-環(huán)境因子相似性值?？紤]到miRNA相似性矩陣和環(huán)境因子相似性矩陣可能有不同的結構和拓撲特征，隨機游走在兩個網(wǎng)絡中的最優(yōu)的迭代步數(shù)可能不一致。因此，本發(fā)明設定兩個參數(shù)l和r分別來表示隨機游走在miRNA相似性網(wǎng)絡和EF相似性網(wǎng)絡中最大迭代步數(shù)；l和r的取值通過交叉驗證確定，本發(fā)明實施例中，l和r的取值通過交叉驗證最終確定分別為4和2。其迭代過程可表示為：

在miRNA相似性網(wǎng)絡游走：

RD_L(t)＝a×RM×RD(t-1)+(1-a)A

其中，a為固定參數(shù)【a的取值通過交叉驗證確定，本發(fā)明實施例中，a的取值通過交叉驗證最終確定為0.8】，RD(t-1)為游走t-1步后預測的miRNA-環(huán)境因子得分矩陣，RD_L(t)為在miRNA相似性網(wǎng)絡上游走t步后預測的miRNA-環(huán)境因子得分矩陣；RD(0)＝A；

在環(huán)境因子相似性網(wǎng)絡游走：

RD_R(t)＝a×RD(t-1)×RE+(1-a)A

其中，RD_R(t)為在環(huán)境因子相似性網(wǎng)絡上游走t步后預測的miRNA-環(huán)境因子的得分矩陣；

最終的輸出游走t步后預測的miRNA-環(huán)境因子得分矩陣RD(t)：

當隨機游走在兩個網(wǎng)絡中的步數(shù)超過最大迭代參數(shù)l或r時，迭代過程終止。

對于任意兩個已知的miRNA和環(huán)境因子，RD(t)的值代表兩者存在關系的可能性，表示已知的miRNA和環(huán)境因子之間潛在的關聯(lián)關系。

(2)基于多標簽學習方法預測新的miRNA-環(huán)境因子之間潛在相互關系

本發(fā)明將miRNA-環(huán)境因子之間關系預測問題轉化為概率事件。假設需要預測miRNA m_i和環(huán)境因子e_j之間關系。若在已知的miRNA中，與某個miRNA最相似的K個miRNA中有10個miRNA和環(huán)境因子e_j有關系，則認定這個miRNA有70％概率和環(huán)境因子e_j存在關聯(lián)關系。并且，在與miRNA m_i最相似的K個miRNA中有10個miRNA和環(huán)境因子e_j存在關聯(lián)關系，那么miRNA m_i也有70％概率和環(huán)境因子e_j有關系。

假定有m個已知的miRNA和e個環(huán)境因子，則miRNA-環(huán)境因子矩陣ME的維度為m行和e列。miRNA m_i和環(huán)境因子e_j之間存在關系的概率定義為P_ij，代表miRNA m_i和環(huán)境因子e_j之間存在關系的可信程度。當已知miRNA m_i和環(huán)境因子e_j之間存在關系時，P_ij＝ME(i,j)＝1；否則，當已知miRNA m_i和環(huán)境因子e_j之間不存在關系時，P_ij＝ME(i,j)＝0。給定一個新的miRNA m_c，其與環(huán)境因子e_j存在關系的概率記為P(m_c,e_j)。基于miRNA相似性，選取與miRNA m_c相似性最大的K個miRNA，即選擇miRNA m_c的K個最近鄰居，N_c表示miRNA m_c的K個最近鄰居。根據(jù)訓練樣本的統(tǒng)計信息，應用條件概率和貝葉斯決策理論來預測測試樣本。本發(fā)明按如下方法計算P(m_c,e_j)：

其先驗概率可利用已知的m個miRNA來估計：

表示和環(huán)境因子e_j有關的miRNA的先驗概率，m表示已知的miRNA的個數(shù)；

后驗概率可以按照以下方法計算：

其中，s的取值通過交叉驗證得到，本實施例最后s取值為10；e(i)表示miRNA m_c的K個最近鄰居中，本身和e_j有關系，且其對應的K個最近鄰居中，有i個miRNA和環(huán)境因子e_j有關系的miRNA的個數(shù)；

e′(i)表示miRNA m_c的K個最近鄰居中，本身和e_j有關系，且其對應的K個最近鄰居中，有i個miRNA和環(huán)境因子e_j沒有關系的miRNA的個數(shù)。

四、實驗驗證

1.網(wǎng)絡分析

miRNA-環(huán)境因子網(wǎng)絡包含224個miRNA，124個環(huán)境因子和729條miRNA-環(huán)境因子關系。圖2表示miRNA-環(huán)境因子關系中環(huán)境因子的度的分布。行坐標代表度，縱坐標代表miRNA的個數(shù)。從圖2中可以發(fā)現(xiàn)大多數(shù)的環(huán)境因子的度都為1。也就是只有1個miRNA和這些環(huán)境因子有關系。其中，吉西他濱的度最高，有56個miRNA與其有關系。

本發(fā)明利用Cytoscape的ClusterViz插件對miRNA-環(huán)境因子關系網(wǎng)絡進行聚類分析。圖3表示利用ClusterViz聚類得到的3個模塊。三角形表示環(huán)境因子，矩形表示miRNA。從這些模塊可以發(fā)現(xiàn)：環(huán)境因子可以調控功能相似的miRNA。例如，模塊3中4個環(huán)境因子(DDT，E2，BPA和電離輻射)和let-7家族都有關系。

2.評價指標

為了驗證本發(fā)明的有效性，本發(fā)明采用十倍交叉驗證來測試MEI-BRWMLL的預測性能。將已知的miRNA-環(huán)境因子已知的關系隨機分為十份，隨機選取一份作為測試集，剩下額九份作為訓練集。在每次交叉驗證實驗中，將測試集里的關系數(shù)據(jù)集中關系刪除，利用剩下的九份訓練集中的已知信息來預測測試集中關系。若刪除后，某個環(huán)境因子沒有一個miRNA與其關聯(lián)，則采用多標簽學習方法來進行預測；否則，采用雙向隨機游走方法進行預測。從而，可以獲得刪除的miRNA-環(huán)境因子得分/存在關系的概率。

針對測試集中每對被刪除關系的miRNA-環(huán)境因子【即實際存在關系的miRNA-環(huán)境因子】，對預測得到的miRNA-環(huán)境因子得分/存在關系的概率設定閾值，如果得分/概率大于這個閾值，則被認為是一個true positive(TP，表示miRNA-環(huán)境因子實際存在關系，預測也存在關系)；否則為false negative(FN，表示miRNA-環(huán)境因子實際存在關系，而預測不存在關系)。另外，對于已知沒有關系的miRNA-環(huán)境因子【即實際不存在關系的miRNA-環(huán)境因子】，若得分/概率大于這個閾值，被認為是一個false positive(FP，表示miRNA-環(huán)境因子實際不存在關系，而預測存在關系)；否則為true negative(TN，表示miRNA-環(huán)境因子實際不存在關系，而預測也不存在關系)。從而進一步計算出True-positive rate(TPR)和False-positive rate(FPR)。TPR和FPR的計算公式如下：

通過改變閾值，可以計算得到不同的TPR和FPR，進而畫出ROC曲線，并且計算該曲線下方的面積可以得到AUC值，AUC值被用來表示全局預測性能。本發(fā)明將MEI-BRWMLL和另外兩個miRNA-環(huán)境因子關系預測方法(miREFScan和miREFRWR)比較。圖4表示三種方法在miRNA-環(huán)境因子關系預測性能。從圖4可以發(fā)現(xiàn)，MEI-BRWMLL的AUC值為0.8208。另外兩種方法miREFScan和miREFRWR的AUC值分別為0.7963和0.7905。實驗結果表明：本發(fā)明要優(yōu)于現(xiàn)在的兩種方法。

3.案例分析

3,3-二吲哚基甲烷(DIM)來源于吲哚-3-甲醇消化產(chǎn)物，是一種存在于十字花科植物如西蘭花，抱子甘藍，卷心菜和甘藍菜的化合物。越來越多的研究表明：3,3-二吲哚基甲烷和許多癌癥有著密切關系。例如，3,3-二吲哚基甲烷能夠抑制結腸癌細胞中基因HDAC1，HDAC2和HDAC3的表達。表1為MEI-BRWMLL預測的和3,3-二吲哚基甲烷有關系的前15個預測的miRNA。在這預測的15個miRNA中，9個miRNA在最近的文獻被驗證和3,3-二吲哚基甲烷有關系。排名第1的為hsa-mir-146a，根據(jù)文獻表明3,3-二吲哚基甲烷能夠誘導hsa-mir-164a表達，且hsa-mir-164a能夠抑制基因MTA-2，IRAK-1和NFKB的表達，進而減少胰腺癌細胞轉移。排名第2的為hsa-mir-16。有文獻表明在腦CD4+T細胞中，3,3-二吲哚基甲烷促進hsa-mir-16的表達，從而抑制靶基因EAE的表達。最近的研究表明3,3-二吲哚基甲烷和排名6，8，12的hsa-mir-181a，hsa-mir-125b，hsa-mir-34a有關系。3,3-二吲哚基甲烷抑制這些miRNA表達，從而可以誘導肝臟中caspase-2基因表達。排名第9的為hsa-mir-200b。研究表明3,3-二吲哚基甲烷可以上調乳腺癌的SKBR3細胞系中的hsa-mir-200b。排名第11為hsa-mir-221，研究表明在胰腺癌總，3,3-二吲哚基甲烷可以改變hsa-mir-221的表達，從而改變PTEN，p27，p57和PUMA基因表達。最近研究發(fā)現(xiàn)，在前列腺癌中，3,3-二吲哚基甲烷可以上調hsa-let-7e的表達，從而抑制基因EZH2的表達。根據(jù)文獻表明在乳腺癌中，3,3-二吲哚基甲烷和赫賽汀共同作用來促進hsa-mir-200c的表達，從而抑制其靶基因的表達。另外，本發(fā)明也發(fā)現(xiàn)了一些新的miRNA，如hsa-mir-24，hsa-mir-155，hsa-mir-233，hsa-mir-181b，hsa-mir-126和hsa-mir-222。雖然這些miRNA分子機制仍然未知，這可以留給生物學家通過生物實驗方法來驗證期生物功能。

表1預測的miRNA