生物信息學分析篩選結直腸癌靶基因及評估預后的價值

江秀玲，胡有根，湯元杰*

(1.揚州市職業大學醫學院，江蘇 揚州 225003；2.武警江蘇總隊醫院，江蘇 揚州 225003)

0 引言

結直腸癌(CRC)是第三大常見癌癥類型[1]。由于結腸癌高侵襲和高轉移特性，CRC已成為癌癥導致的死亡的主要類型之一[2]。盡管當前CRC的治療方法和手段不斷提升，但CRC患者的存活率依然很低[3]。因此，深入了解CRC分子和遺傳網絡，探索識別新的生物分子標志物極其重要。

長鏈非編碼RNA(lncRNAs)是長度大于200個核苷酸的RNA，在細胞增殖、調亡、細胞周期阻滯、細胞遷移和侵襲等多種生物學過程中發揮重要調節作用[3]。近年來，在CRC中，已發現數百種表達異常的lncRNAs，其中有發揮增強子功能的，能夠調控相鄰的編碼基因的表達，制約疾病的發生[4-8]。MicroRNA是長約22nt的短鏈RNA，能夠通過抑制目的基因的翻譯或降解，從而調節目的基因的表達。然而在實際調控過程中不僅僅是簡單的microRNA-mRNA的沉默機制，還有更為復雜的調控網絡。一些非編碼的RNA同樣存在于microRNA的結合位點，在細胞中起到miRNA海綿(miRNA sponge)的作用，進而解除miRNA對其靶基因的抑制作用，升高靶基因的表達水平，也因此構建了龐大的ceRNA網絡(ceRNETs)，這一作用機制被稱為競爭性內源RNA(ceRNA)機制[9-11]。

在本研究中，我們研究首次通過構建結直腸癌差異表達基因的ceRNA網絡，挑選出與結直腸癌密切相關的關鍵ceRNA基因，其次通過分析WGCNA網絡的關鍵基因并將關鍵基因與上述ceRNA網絡中的關鍵miRNA做交集分析，將挑選出的miRNA通過生存分析分析篩選出與生存密切相關的miRNA并通過數據庫驗證其實際意義。相關研究可以作為治療結直腸癌的新型靶向藥物以及新型結直腸癌標志物。

1 材料與方法

1.1 數據收集和預處理

通過TCGA獲取結腸癌表達數據，取表達數據較完整的樣本(449個樣本)，并篩選miRNA、mRNA以及lncRNA表達數據；通過GDC Data Portal對TCGA中結腸癌相關臨床數據進行統計，共統計439個病人，整理為標準表格。

1.2 RNA差異表達分析

使用DESeq2包分別對miRNA、mRNA、lncRNA進行差異表達分析。經過分析，獲得各RNA的差異表達列表(FDR＜0.05)。

1.3 TCGA結腸癌中差異表達mRNA的GO分析

GO主要可以分為三個主要的類群，生物學進程 (Biological Process，BP)，分 子 功 能 (Molecular Function，MF)和細胞組分(Cellular Component)。將分析得到的差異基因基于GO數據庫分別從BP、MF、CC這三個層面進行的GO注釋，得到基因參與的所有GO。采用Fisher檢驗計算每個GO的顯著性水平(P-Value)，從而篩選出差異基因、負相關基因富集的顯著性GO。

1.4 TCGA結腸癌中差異表達mRNA的KEGG Pathway分析

將分析得到的差異基因基于KEGG數據庫進行Pathway注釋，得到差異基因參與的所有Pathway Term，采用Fisher檢驗計算Pathway的顯著性水平(P-Value)，從而篩選出差異基因、負相關基因富集的顯著性Pathway Term。

1.5 ceRNA網絡構建

根據上述步驟中差異分析結果，將lncRNA、mRNA和miRNA分為上調和下調，再根據表達情況將上調的mRNA、lncRNA與下調的miRNA對應、下調的mRNA、lncRNA與上調的miRNA對應使用Miranda軟件預測mRNA、lncRNA與miRNA的結合。其對3’UTR的篩選依據主要是從序列匹配、miRNA與mRNA雙鏈的熱穩定性以及靶位點的保守性三個方面進行分析。綜合這3條原則，miRanda選取與miRNA序列互補的3’UTR所對應的排名前十位的基因，作為miRNA的候選靶基因，對于多個miRNA對應于同一靶位點的情況，miRanda使用貪心算法(Greedy Algorithm)選取其中得分最高且自由能最低的那一對。設定Score閾值為150，Energy閾值為-20，去掉重復項，并將上下調分別獲取的結果整合后進行展示。若lncRNA與mRNA與同一個miRNA結合，將此lncRNA視為相應mRNA的ceRNA。將miRNA與mRNA結合預測與miRNA與lncRNA結合預測結果按miRNA取并集。

1.6 WGCNA網絡中關鍵基因與ceRNA網絡中關鍵基因的交集分析

使用加權的表達相關性對基因的表達數據進行共表達分析，并構建共表達網絡。接著基于加權相關性，進行層級聚類分析，并根據設定標準切分聚類結果，獲得不同的基因模塊，根據基因間表達量進行聚類分析，得到的各模塊間的相關性圖；WGNA認為基因之間的簡單的相關性不足以計算共表達，所以它把鄰接矩陣轉換為拓撲重疊矩陣，以降低噪音和假相關，獲得的新距離矩陣(Topological Overlap Matrix，TOM)。繪制TOM圖；根據上述分析共獲得8個基因共表達模塊；根據TOM值對模塊中基因的共表達關系進行分析。

獲取WGCNA與ceRNA網絡中miRNA的交集取WGCNA中差異基因最多的“turquoise(青色)”模塊進行進一步分析，將模塊中的miRNA，與ceRNA中miRNA取交集，發現20個共有miRNA。分析ceRNA及WGCNA的關鍵RNA k-core表示在一個關聯網絡中，評估基因在特定網絡中的核心程度；分別計算WGCNA共表達網絡與ceRNA網絡中各基因的k-core，并結合ceRNA預測結果，在上述步驟挑選出關鍵miRNA的基礎上，挑選WGCNA所選模板中與預后關聯顯著的lncRNA及mRNA，構建ceRNA網絡。利用Blast2Go及TopGO分析ceRNA網絡中靶點RNA的已知資料和信號通路參與情況，發現不少mRNA與結腸癌關聯性較大、或參與癌癥相關通路，呈現具體注釋結果。

1.7 生存分析

TCGA數據庫除提供組學測序數據外，還可獲取生存數據，基于TCGA中的表達數據和生存數據能夠有效地進行生存分析。對上述步驟構建的ceRNA中節點、臨床數據與總生存率(OS)進行COX回歸分析。

1.8 GEO數據庫驗證

通過GEO數據庫查找結腸癌相關轉錄組測序數據，選取結腸癌相關miRNA數據集GSE115513，利用數據集對上文分析獲得的關鍵miRNA進行差異表達分析，并繪制表達箱型圖。

2 實驗結果

2.1 結腸癌和對照樣本間差異表達的RNA

基于TCGA中獲取的結腸癌RNA表達數據，我們使用了DESeq2分別對miRNA、mRNA、lncRNA進行差異表達分析，篩選得出428個差異表達miRNA、4318個差異表達mRNA、128個差異表達lncRNA，并繪制了火山圖(圖1)。

圖1 結腸癌中差異表達RNA的火山圖

2.2 差異表達RNA的富集功能

將上述步驟中差異篩選結果進一步分析，我們將、上調mRNA和下調mRNA作為研究目標，分別從生物學進程(BP)，分子功能(MF)和細胞組分(CC)這三個角度對樣本進行GO注釋。全部差異mRNA在跨膜運輸、免疫炎癥應答、細胞因子、細胞外基質等活動中顯著富集(圖2A)；上調差異mRNA在細胞周期、有絲分裂、生長因子、細胞外基質等活動中顯著富集(圖2B)；而下調mRNA則在跨膜運輸、離子運輸、受體功能、GPC等活動中顯著富集(圖2C)。為更直觀顯示結果，我們挑選了富集最顯著的20條GO條目繪制了富集散點圖(圖3)，跨膜運輸、離子運輸和小分子代謝等基因數目占據量較大；G蛋白偶聯的嘌呤核苷酸受體通路、鋅離子細胞應答反應、脂肪酸β氧化反應Richfactor較大，而Rich factor越大，則表示富集的程度越大。

圖2 GO分析結果

圖3 GO富集性散點圖

除了對差異基因進行GO注釋之外，我們還通過KEGG Pathway分析差異RNA的富集功能。全部差異mRNA在蛋白質消化和吸收、谷氨酸能突觸通路等富集；上調差異mRNA在細胞周期通路富集；下調差異mRNA在礦物吸收、鈣離子和趨化因子信號通路富集(圖4)。

圖4 用KEGG數據庫對差異基因進行富集分析

2.3 競爭內源性RNA網絡揭示了長鏈非編碼RNA和信使RNA的競爭內源性機制

在復雜的ceRNA調控網絡中，調控機制不僅存在miRNA-mRNA沉默機制，而存在更為復雜的lncRNA-miRNA-mRNA調控網絡，一些非編碼的RNA同樣存在與microRNA的結合位點，在細胞中起到miRNA海綿(miRNA sponge)的作用，進而解除miRNA對其靶基因的抑制作用，升高靶基因的表達水平，也因此構建了龐大的ceRNA網絡。

在本研究中，我們根據miRNA篩選出miRNA的節點miRNA，然后根據節點miRNA通過miRanda軟件預測出mRNA、lncRNA與miRNA的結合，共獲得94個miRNA、1859個mRNA、117個lncRNA。通過WGCNA加權網絡相關性進一步分析，我們篩選出差異表達最多的基因模板-turquoise 青色模板，其中包含77個miRNA、45個lncRNA以及1381個mRNA(圖5A-C)。然后，將上述篩選出的miRNA取交集，獲得了20條miRNA(表1)。對下述miRNA計算其k-core值，并選出k-core值最高的關鍵miRNA hsa-mir-874、hsa-mir-92a-1、hsa-mir-940，并在此基礎上挑選WGCNA所選模板中與預后關聯顯著的lncRNA及mRNA，成功構建了ceRNA網絡(圖6)。

圖5 WGCNA分析

圖6 關鍵miRNA的ceRNA互作網絡

表1 共有miRNA在WGCNA及ceRNA中k-core

2.4 關鍵驅動基因的存活分析

本研究基于TCGA中獲取的表達數據以及生存數據分別對ceRNA節點中的miRNA，mRNA和lncRNA進行COX回歸分析，在關鍵miRNA的COX回歸分析中發現COX分析結果中與預后有顯著關系的為hsa-mir-874(圖7A)；而關鍵ceRNA 網絡中 lncRNA(LINC00858、 LINC00896、MIR503HG、PRR7-AS1、SNHG7)、mRNA(JMJD7-PLA2G4B、ANKRD23、ATP2A1、C6orf223、CPNE7、ENGASE、IZUMO4、LZUM04、LZTS3、MC1R、NPEPL1、PHF19、PLA2G4B、PPFIA4、SEZ6L2、SFXN3、SLC25A27、SLC6A20、SOX12、SRCIN1、TEX22、TMEM178B、TRAF5、VEGFA、WT1、ZGLP1、ZNF514)的COX回歸分析中均與預后有顯著關系(圖7B-C)。

圖7 Cox回歸分析基因表達與預后的關系

2.5 選擇 hsa_mir_874/LINC00896、LINC00858 進行數據庫驗證

基于上述生存預后分析，我們選擇了具有統計學意義的hsa_mir_874，在結直腸癌中屬于下調表達基因(圖8A)。GEPIA網站驗證靶基因lncRNA(MIR503HG、PRR7-AS1、LINC00896、LINC00858)在 結 腸 癌 中 的 表 達，僅 有LINC00858[15]、PRR7-AS1、LINC00896 在結直腸癌高表達 (圖8B)。GEPIA數據庫驗證靶基因mRNA(C6orf223、NPEPL1、PPFIA4、SEZ6L2、SFXN3、SOX12、TEX22、VEGFA、WT1、ZNF514)在結腸癌中表達情況，其中C6orf223、SEZ6L2、VEGFA在結腸癌中高表達(圖8C)。因此，我們推測在結直腸 癌 中 hsa_mir_874/LINC00896、LINC00858、PRR7-AS1/C6orf223、SEZ6L2、VEGFA存在ceRNA機制。

圖8 ceRNA網絡相關因子在結腸癌中的表達情況

3 討論

結直腸癌(CRC)是第三常見的癌癥，每年約136萬新病例，有70萬人死于結腸癌，是全球第四大癌癥導致死亡的癌癥類型[12]。大約10%的CRC病例是遺傳性的，高達90%是散發性的(沒有家族史或遺傳易感性)。已經確定導致CRC發展的幾個危險因素，主要包括不健康的生活習慣，如缺乏運動、吸煙、食用紅肉和加工肉類以及飲酒，甚至包括肥胖、糖尿病2型和炎癥性腸病(IBD)在內的一些疾病也與CRC[13]發病風險增加有關。目前臨床上檢測CRC的方法可以分為無創性和侵入性兩種成像技術，但最近的研究表明，它們在篩查CRC時敏感性和特異性比較低[14]。因此，越來越多研究從分子和遺傳兩個方面進行深入的研究，為結直腸癌的早期診斷尋找更高效的生物標志物。

隨著基因組學技術水平的不斷提升，lncRNAs被認為是基因組的一種重要功能表達。LncRNAs是一種長度大于200個核苷酸的RNA，可通過與DNA、RNA和蛋白質分子相互作用來調節疾病的生理和病理過程[15,16]。越來越多的lncRNA被報道在CRC進展中發揮重要作用。例如有研究發現在CRC組織和細胞系中lncRNA CPS1-IT1表達顯著下調，并且CPS1-IT1表達低的患者預后較差。進一步體外實驗發現CPS1-IT1可顯著抑制細胞增殖、遷移和侵襲能力，加速細胞凋亡，從而抑制上皮間充質轉化(epi- mesenchymal transition,EMT)，在CRC中發揮抑瘤作用[17]。

一些lncRNA作為競爭內源性RNA(competing endogenous RNAs, ceRNAs)在基因表達調控中發揮重要作用。ceRNA假說描述了一個復雜的由miRNA介導的轉錄后調控網絡，即通過共享一個或多個miRNA響應元件(MREs)，蛋白編碼RNA和非編碼RNA相互競爭與miRNA結合，進而相互調節表達[9]。lncRNA下調導致更多的miRNA分子能夠自由結合到含有相同MREs的mRNA上，從而導致其蛋白表達水平降低，相反的，lncRNA過表達導致miRNA與mRNA結合的分子減少，其蛋白表達水平則升高。盡管ceRNA假說得到了很多研究者的支持，但與結直腸癌相關的ceRNA機制探索較少。若能發現有效的結直腸癌ceRNA網絡作用機制，則能更深入從分子層面分析結直腸癌發展機制，發現更多新型生物標志物。

在本研究中，我們基于TCGA數據庫中可獲得的結直腸癌表達數據以及臨床數據，分析了結直腸癌中表達差異的mRNA、miRNA和lncRNA。使用Miranda軟件以及WGCNA加權基因共表達網絡獲取關鍵miRNA，構建關鍵miRNA的ceRNA互作網絡；同時，我們根據預后分析，選擇與生存最密切相關的miRNA hsa_mir_874進行后續驗證，已有文獻報道，hsa-mir-874可抑制結腸癌進展，增加細胞凋亡，增加化療敏感性[18]。我們研究發現hsa_mir_874的靶基因lncRNA為LINC00858、LINC00896、MIR503HG、PRR7-AS1均與生存預后密切相關。這其中，LINC00858已有研究報道能夠促進骨肉瘤增殖和非小細胞肺腺癌發展[19,20]，并有研究證明其在結直腸癌組織中表達量高。而hsa_mir_874的靶基因mRNA C6orf223、IZUMO4、NPEPL1、PPFIA4、SEZ6L2、SFXN3、SOX12、TEX22、VEGFA、WT1、ZNF514 均與結腸癌預后相關，只有C6orf223、SEZ6L2、VEGFA在結腸癌中高表達。其中越來越多地研究表明VEGFA也在結直腸癌中過表達，并且干擾VEGFA的低表達能夠有效抑制結直腸癌的發展[21,22]。通過生存預后分析，我們推測lncRNA LINC00858、LINC00896、PRR7-AS1與基因C6orf223、SEZ6L2、VEGFA共享一個或多個hsa_mir_874位點，隨后我們在GEPIA數據庫驗證了這一猜想。

本研究首次通過WGCNA關鍵基因分析與ceRNA網絡篩選出結直腸癌中的靶向lncRNA與生物標志物，發現lncRNA LINC00858、LINC00896、PRR7-AS1通過 ceRNA 機制影響結直腸癌的預后，但其具體的作用機制還有待深入探索，這也提示我們，LINC00858、 LINC00896、PRR7-AS1有望作為結直腸癌診斷和治療的新型生物標志物或結直腸癌新型靶向藥物。

綜上所述，我們提出了一種更加多角度構建ceRNAs網絡的方法，并確定 LINC00896、LINC00858、PRR7-AS1可以作為結直腸癌治療生物標志物的關鍵基因，INC00858、LINC00896通過ceRNA機制影響結直腸癌的預后。往后研究應更加重視ceRNA網絡的構建篩選關鍵基因或RNA競爭內源性相互作用的驗證，以探尋更為有效的生物標志物。