朱東月,朱 平

(江南大學理學學院,中國江蘇 無錫 214122)

乳腺癌是全球女性發(fā)病率最高的癌癥之一,相關數(shù)據顯示大約每18 秒就有一名女性被確診為乳腺癌[1]。乳腺癌是女性死亡的主要原因之一,占全球癌癥死亡人數(shù)的14%[2]。目前,對于乳腺癌的診斷仍缺乏高效敏感的分子標志物[3]。因此,探索相關生物標志物對乳腺癌的診斷尤為重要。

癌癥的大多數(shù)生物標志物主要是基于高通量基因組學數(shù)據分析獲得的。RNA 測序技術的發(fā)展有助于揭示基因組的復雜性。已有資料顯示,超過90%的基因組可以被轉錄,但只有1%～2%的轉錄產物可以編碼蛋白質,而大多數(shù)轉錄產物是長鏈非編碼RNA(long noncoding RNA,lncRNA)和微 RNA (microRNA,miRNA)[4]。lncRNA 是一類重要的非編碼RNA,長度約為200 個核苷酸,由RNA 聚合酶Ⅱ轉錄,經選擇性剪切加工而成熟。lncRNA 在復雜的生命過程中扮演著重要的調控角色,如參與調節(jié)免疫反應、惡性細胞的增殖和轉移[5]。為了揭示 lncRNA 的功能,Zhang 等[6]通過構建lncRNA-mRNA 網絡,發(fā)現(xiàn)了7 個與心肌梗死發(fā)生和復發(fā)相關的lncRNA。miRNA 是一類重要的非編碼小RNA,其長度約為22 個核苷酸,可以通過與堿基配對來指導復合物識別mRNA 的3＇端非翻譯區(qū),從而對靶基因起到調控的作用[7]。競爭性內源RNA(competing endogenous RNA,ce-RNA)假說主要是指不同類型的RNA 之間通過競爭一個或多個miRNA 的結合位點來調控基因表達[8]。ceRNA 在很多癌癥的預后和治療中都發(fā)揮著重要的作用。在胃癌中,lncRNA HOTAIR 通過結合 miR-331-3P 抑制 HER2 的表達,且 HOTAIR/HER2 的表達水平與晚期胃癌具有顯著相關性[9];在腎透明細胞癌中,PTENP1 作為競爭性內源RNA 抑制癌癥的發(fā)展[10];在乳腺癌中,linc-ROR通過調節(jié)miRNA 促進癌細胞的發(fā)展和轉移[11]。

本文利用差異 mRNA、lncRNA 和 miRNA 構建ceRNA 網絡,并結合網絡的拓撲性質和RNA的表達數(shù)據計算該網絡節(jié)點的權重,然后基于隨機森林的改進算法篩選出一組RNA,該組RNA在識別乳腺癌方面具有較高精度和效率,結合其生物功能分析發(fā)現(xiàn),該組RNA 可作為乳腺癌的生物標志物。

1 材料與方法

1.1 RNA差異表達分析

TCGA 數(shù)據庫[12](https：//cancergenome.nih.gov/)提供了各種基因組的高通量數(shù)據,包括mRNA、lncRNA 和miRNA 的表達數(shù)據。首先,從TCGA數(shù)據庫中下載與乳腺癌有關的mRNA、lncRNA 和miRNA 的表達矩陣。其中mRNA、lncRNA 表達數(shù)據包括113 個正常樣本和1 109 個乳腺癌樣本,miRNA 表達數(shù)據包含76 個正常樣本和752 個乳腺癌樣本。然后,利用R 軟件中的edgeR 包[13]對RNA 數(shù)據進行標準化和差異分析,以獲得差異mRNA、lncRNA 和miRNA。在篩選差異表達mRNA 時,將篩選參數(shù)設置為：|log2FC|＞2,P＜0.001; 在篩選差異lncRNA 和差異miRNA 時,將篩選參數(shù)設置為：|log2FC|＞1.5,P＜0.05。

1.2 ceRNA網絡的建立

本文研究的ceRNA 網絡包含lncRNA-miRNA和miRNA-mRNA 兩種靶向關系。miRcode 數(shù)據庫[14](http：//www.mircode.org/)覆蓋了完整的GENECODE 注釋的轉錄組,可用于預測lncRNA 和miRNA 之間的靶向關系。通過該數(shù)據庫,獲得了本研究所需的lncRNA-miRNA 關系對。此外,使用 miRDB(http：//mirdb.org)、miRTarBase (http：//mirtarbase.mbc.nctu.edu.tw/)和 TargetScan (http：//www.targetscan.org/vert_72/)在線分析工具[15～17]預測mi-RNA 的靶基因mRNA。為了提高生物信息分析的可靠性,本文利用3 個數(shù)據庫中重疊的miRNA-mRNA 關系對進行后續(xù)的生物學分析?；讷@得的 lncRNA-miRNA 和 miRNA-mRNA 關系對,利用Cytoscape 軟件[18]將ceRNA 網絡可視化。

1.3 ceRNA網絡節(jié)點權重的計算

邊聚集系數(shù)(edge clustering coefficient,ECC)是由Watts 等[19]于1998年提出的,用于衡量網絡中一個節(jié)點與其鄰接節(jié)點之間的親疏程度,是描述網絡拓撲性質的重要特征之一。根據Watts 等[19]對于邊聚集系數(shù)的定義,ceRNA 網絡的邊聚集系數(shù)E(x,y)可以用下式表示：

式中N(x)表示節(jié)點x 鄰接節(jié)點的集合,N(y)表示節(jié)點y 鄰接節(jié)點的集合。由上式可以看出,E(x,y)是一個局部指標,用于衡量網絡中x 和y 兩個節(jié)點之間的親疏程度。E(x,y)的取值范圍為[0,1],其值越大表示兩個節(jié)點之間的關系更為緊密,屬于同一模塊的可能性越大。

皮爾遜相關系數(shù)(Pearson correlation coefficient,PCC)可用于篩選ceRNA 網絡中的mRNA-lncRNA關系對[6]。該系數(shù)可用如下公式表示：

其中n 表示乳腺癌樣本數(shù),Xi表示在i 樣本中x的表達量,Yi表示在i 樣本中y 的表達量,分別表示在樣本中x 和y 的平均表達量。使用值描述x 和y 相關性的強弱。的取值范圍為[0,1],該值越大,mRNA 和lncRNA 的相關性越強。

在ceRNA 網絡中,ECC 值和PCC 值較高的節(jié)點對所在模塊的調控能力較強,該節(jié)點更傾向成為生物標志物[20]。因此,當衡量一個節(jié)點在網絡中的調控程度時,需要綜合考慮ceRNA 網絡的拓撲性質和基因的表達相關性。利用下式定義每個節(jié)點的權重值：

1.4 基于隨機森林的逐步特征選擇方法

隨機森林(random forest,RF)已經廣泛用于癌癥生物標志物的分類和預測[21～22]。該算法可以避免過擬合且不受數(shù)據維度的影響,具有較好的分類精度和預測效果,但不可以降低變量子集的維度。本文將RF 改進為基于隨機森林的逐步特征選擇(stepwise feature selection based on random forest,SFS-RF)算法,用于篩選高效率和低維度的乳腺癌生物標志物,具體步驟為：1) 運用具有5-折交叉驗證的隨機森林算法計算單個變量的分類精度,并按照精度由高到低進行排序; 2) 將每個變量逐次放入到變量集中,如果增加一個變量可以使得分類精度提高,則保留該變量; 反之,則終止計算。算法運行結束時,即可篩選出分類精度最高且數(shù)目最小的變量集,以及相應的分類精度。

將本文采用的SFS-RF 方法與LASSO (least absolute shrinkage and selection operator)[23]、主成分分析(principal component analysis,PCA)[24]進行比較,并利用準確率(accuracy,ACC)、靈敏度(sensitivity,SN)、特異性(specificity,SP)和馬修斯系數(shù)(Matthews correlation coefficient,MCC) 4 個評價指標衡量變量集的分類能力。ACC、SN、SP 和MCC的計算公式如下所示：

式中TP 表示真陽性,即預測為正樣本實際也是正樣本的數(shù)量; TN 表示真陰性,即預測為負樣本實際也是負樣本的數(shù)量; FP 表示假陽性,即預測為正樣本實際為負樣本的數(shù)量; FN 表示假陰性,即預測為負樣本實際為正樣本的數(shù)量。

同時,用受試者操作特征曲線(receiver operator characteristic curve,ROC 曲線)評價變量集的分類精度[25]。曲線下的面積(area under curve,AUC)用于數(shù)值化地反映算法的分類能力[25],其取值范圍在[0,1],當AUC 值越接近1 時,模型的預測準確率越高。

1.5 富集分析

利用BiNGO 插件[26]對生物標志物調控的RNA進行GO (gene ontology)功能富集分析。GO 功能富集主要是從分子功能(molecular function,MF)、生物過程(biological process,BP)和細胞組分(cellular component,CC) 3 個層面探究RNA 的生物功能。利用 KOBAS 數(shù)據庫[27](http：//kobas.cbi.pku.edu.cn/)對生物標志物調控的RNA 進行KEGG (kyoto encyclopedia of genes and genomes)通路富集分析。在后續(xù)的分析中,P＜0.05 的生物功能或者通路被認為具有統(tǒng)計學意義。

2 結果

2.1 差異表達的RNA

利用edgeR 包在正常樣本和癌癥樣本中篩選出的差異表達mRNA、lncRNA 和miRNA 分別有847 個、438 個、152 個。其中,上調 mRNA 有 532個、下調 mRNA 有 315 個;上調 lncRNA 有 341個、下調lncRNA 有97 個; 上調miRNA 和下調 miRNA 依次有 128 個、24 個。

2.2 乳腺癌ceRNA網絡

從miRcode 在線數(shù)據庫中獲得了99 個差異lncRNA-miRNA 相互作用關系對,其中包含27個差異 lncRNA 和 21 個差異 miRNA。利用miRDB、miRTarBase 和 TargetScan 數(shù)據庫對差異miRNA 的靶基因mRNA 進行預測,獲得差異miRNA-mRNA 相互作用關系對39 個,其中包括20 個差異miRNA 和 27 個差異mRNA。最終,通過lncRNA-miRNA-mRNA 關系對建立ceRNA 網絡,該網絡由 24 個 lncRNA、20 個 miRNA 和 27個 mRNA 組成(圖 1)。在 ceRNA 網絡中,節(jié)點表示差異mRNA、lncRNA 或者miRNA,邊表示它們之間的聯(lián)系。

2.3 SFS-RF方法篩選的生物標志物

權重不僅考慮了不同RNA 在ceRNA 網絡中的鄰接節(jié)點聚集情況,同時也考慮了RNA 之間表達的相關性。通過公式(3)計算ceRNA 網絡節(jié)點的權重,發(fā)現(xiàn)的差異RNA 有27 個。用SFS-RF 方法對這27 個RNA 進行篩選,獲得乳腺癌生物標志物,以用于高精度識別乳腺癌樣本。最終從27個差異RNA 中篩選出一組生物標志物CHL1-AS2、LINC00466 和 LINC00337。

以這3 個RNA 作為識別乳腺癌樣本的變量時,評價指標 ACC、SN、SP、MCC 和 AUC 的值依次為 0.98、0.85、0.98、0.84 和 0.98 (ROC 曲線見圖2A),說明該組RNA 對識別乳腺癌樣本具有很高的靈敏度和特異性。lncRNA、miRNA 和mRNA 均由基因轉錄而來,因此進一步利用GEO 數(shù)據庫[28]驗證該組RNA 對乳腺癌樣本的識別情況。從GEO 數(shù)據庫中下載基因芯片GSE75333 (包含3個正常樣本和3 個乳腺癌樣本),并用ROC 曲線檢驗這3 個差異lncRNA 識別乳腺癌的靈敏度和特異性。結果顯示：CHL1-AS2、LINC00466 和 LINC00337 具有較高的識別精度(AUC=0.90)。此外,為了探究該組生物標志物在正常樣本和乳腺癌樣本中的表達情況,利用箱線圖將它們在兩組樣本中的表達水平進行可視化。結果如圖2B～D 所示,CHL1-AS2 在正常樣本中顯著高表達(P=6.8E-07);而LINC00466 和LINC00337 在正常樣本中低表達,在癌癥樣本中顯著高表達(P 依次為4.3E-59、7.4E-63)。這說明該組RNA 的異常表達可能與乳腺癌的發(fā)病機制相關。

2.4 方法比較

腫瘤基因表達具有多變量、高維度和高冗余的特點,為了獲得有效且可靠的信息,LASSO 回歸和PCA 已廣泛運用于生物標志物的篩選[29～30]。根據TCGA 數(shù)據庫中下載的乳腺癌RNA 數(shù)據,利用SFS-RF、LASSO 和PCA 方法篩選生物標志物,用 ACC、SN、SP、MCC 和 ROC 曲線衡量它們篩選生物標志物的分類精度。3 種方法篩選出的生物標志物數(shù)量及分類精度見表1。表1 直觀地顯示：在分類精度相近的情況下,SFS-RF、LASSO 和PCA 方法篩選出來的生物標志物數(shù)量依次為3、19 和9,顯然SFS-RF 方法篩選出的生物標志物最少。在乳腺癌的研究中,通過本文提出的方法,只需檢測3 個基因的轉錄水平就能達到與其他方法(如LASSO 回歸和PCA)相近的識別精度。因此,SFS-RF 方法在篩選生物標志物時具有更高的效率。

圖1 乳腺癌ceRNA 網絡正方形表示lncRNA; 圓表示miRNA; 三角形表示mRNA; 綠色節(jié)點表示下調RNA; 紅色節(jié)點表示上調RNA。Fig.1 Breast cancer ceRNA networkSquares represent lncRNA; Circles represent miRNA; Triangles represent mRNA; Green nodes indicate down-regulated RNAs;Red nodes indicate up-regulated RNAs.

2.5 生物功能分析

在乳腺癌ceRNA 網絡中,LINC00466、CHL1-AS2 和LINC00337 通過競爭性結合15 個miRNA調控32 個差異RNA 在乳腺癌中的表達(圖3)。為了探究這些差異RNA 對乳腺癌發(fā)病機制的影響,對這些RNA 進行GO 和KEGG 富集分析。GO 分析結果顯示,這些RNA 顯著富集在染色體的形成、細胞增殖的正調控、蛋白激酶活性、姐妹染色單體的分離、調節(jié)磷酸鹽代謝過程等387 個功能類別,其中在生物過程(BP)和分子功能(MF)中排名前10 的條目如圖4A 所示。KEGG 通路分析結果表明,它們主要富集在TNF 信號通路、PI3KAkt 信號通路、MAPK 信號通路、細胞衰老和乳腺癌等15 條通路,排名前10 的富集通路見圖4B。從上述結果可以看出,該組生物標志物通過調控ceRNA 網絡中RNA 的表達,參與乳腺腫瘤細胞的增殖、遺傳物質的合成和信號傳導等生物過程。

圖2 生物標志物的ROC 曲線和差異表達(A) ROC 曲線; (B～D) 生物標志物在正常樣本和乳腺癌樣本中的差異表達。NP：正常樣本; TP：癌癥樣本。Fig.2 ROC curve and differential expression of biomarkers(A)ROC curve;(B～D)Differential expression of biomarkers in normal and cancer samples.NP：Normal samples;TP：Cancer samples.

表1 SFS-RF、LASSO 和 PCA 方法的比較Table 1 Comparison of SFS-RF,LASSO and PCA methods

值得注意的是,圖 3 顯示 hsa-mir-145 與hsa-mir-183 通過一個差異lncRNA C5orf17 發(fā)生間接相互作用。由此,本文獲得了一個有意義的miRNA—lncRNA—miRNA 調控軸,即 hsa-mir-183—C5orf17—hsa-mir-145。該調控軸是生物標志物 LINC00466、CHL1-AS2 和 LINC00337 相互調控的中間渠道。

3 討論

乳腺癌是女性中發(fā)病率最高的惡性腫瘤之一,通過構建其加權ceRNA 網絡,可以挖掘出與乳腺癌發(fā)生發(fā)展較為緊密的RNA,并且深入探討這些RNA 參與的生物過程及其具有的生物功能,可以為乳腺癌的發(fā)病機制探究提供一定的理論支撐。本研究通過對乳腺癌的ceRNA 網絡加權提取了27 個差異RNA,然后采用SFS-RF 方法篩選出一組可以作為乳腺癌生物標志物的RNA,即LINC00466、CHL1-AS2 和 LINC00337。GEO 數(shù)據集驗證結果顯示,這3 個差異lncRNA 在識別乳腺癌樣本方面具有很高的可靠性和準確性(AUC=0.90)。

在這組RNA 中,LINC00466 已被證實可以作為乳腺癌的競爭性內源RNA,并參與乳腺腫瘤的發(fā)生、發(fā)展、增殖和轉移等過程[31]。CHL1-AS2 和LINC00337 是本文發(fā)現(xiàn)的兩個新的乳腺癌競爭性內源RNA。目前,暫未發(fā)現(xiàn)這兩種RNA 在乳腺癌識別中的應用價值,但已有研究表明它們對子宮內膜異位癥[32]、肺腺癌[33]和胃癌[34]的發(fā)生有一定的影響。Zhang 等[32]使用定量聚合酶鏈反應(qPCR)在異位子宮內膜中發(fā)現(xiàn)CHL1-AS2 的表達水平明顯高于正常子宮內膜,故認為這可能與子宮內膜異位癥的發(fā)生有關。LINC00337 位于1 號染色體,由5 個外顯子組成,在肺腺癌ceRNA 網絡中,LINC-00337 通過競爭性結合hsa-mir-373 和hsa-mir-519 調控PBK 和KIF23 的表達,進而調節(jié)免疫系統(tǒng),最終影響患者的預后[33]。Hu 等[34]發(fā)現(xiàn)胃癌樣本中LINC00337 的表達水平明顯高于正常樣本,并且其可利用BZH2 作為媒介抑制p21 表達,以促進胃癌細胞的增殖。因此,CHL1-AS2 和LINC00337 可能是乳腺癌潛在的生物標志物,值得進一步研究。

圖3 生物標志物的ceRNA 子網絡Fig.3 ceRNA subnetwork of three biomarkers

在GO 富集分析結果中,我們發(fā)現(xiàn)了幾類與腫瘤細胞增殖、遺傳物質合成及蛋白質代謝密切相關的功能注釋,如濃縮染色體、細胞增殖的正調控、蛋白激酶活性和正調控高分子代謝等功能條目,說明 LINC00466、CHL1-AS2 和 LINC00337 通過競爭性地結合miRNA 調控多個RNA 參與細胞增殖、分化和遺傳物質的合成等生物過程。在KEGG富集分析中,多條通路已被證實與乳腺癌的發(fā)生相關。腫瘤壞死因子(tumor necrosis factor,TNF)可以直接殺傷腫瘤細胞,相關研究表明當TRADD 低表達時,TNF 信號通路殺傷腫瘤細胞的作用減弱[35]。MALAT1 低表達可激活PI3K-Akt 信號通路,進而調節(jié)乳腺腫瘤細胞的轉移[36]。MAPK 信號通路對乳腺癌的發(fā)生發(fā)展有著重要的影響,研究顯示：EGFR 過表達激活MAPK 信號轉導通路,活化原癌基因c-fos,刺激細胞惡性增殖[37]。KEGG通路富集結果與GO 富集結果相互驗證,說明在生物功能和生物通路中,研究這些生物標志物是有意義的。同時,本文獲得了一個有意義的調控軸 hsa-mir-183—C5orf17—hsa-mir-145,其可作為 LINC00466、CHL1-AS2 和 LINC00337 相互調控的中間渠道,且hsa-mir-183[38]和hsa-mir-145[39]已被證實與乳腺癌的發(fā)生、發(fā)展和預后等方面有著密切的聯(lián)系。這些表明LINC00466、CHL1-AS2和LINC00337 可以作為生物標志物,用于高精度識別乳腺癌樣本。

圖4 GO 和KEGG 富集分析Fig.4 GO and KEGG enrichment analyses

總的來講,本文構建了乳腺癌加權ceRNA 網絡,并將其用于研究RNA 分子與乳腺癌發(fā)病機制之間的聯(lián)系。其次,利用SFS-RF 方法,篩選出了一組可用作乳腺癌生物標志物的RNA——LINC00466、CHL1-AS2 和 LINC00337。其中,CHL1-AS2 和LINC00337 作為首次發(fā)現(xiàn)的標志物,為探究乳腺癌ceRNA 調控機制提供了新的思路。