基于lncRNA-mRNA共表達網絡篩選胃癌分期相關的lncRNA*

哈爾濱醫科大學衛生統計學教研室(150081) 叢雨欣 張秋菊 田 偉 張 奇 李 稱 趙 敏 劉美娜

【提 要】 目的 構建lncRNA-mRNA共表達網絡，探索與胃癌分期相關的lncRNA及其調控關系，為研究胃癌的進展機制及尋找胃癌治療的潛在靶點提供依據。方法 利用TCGA數據庫，收集胃癌RNA-seq數據及臨床信息數據；采用表達數量性狀位點(eQTL)分析與加權基因共表達網絡分析(WGCNA)相結合的方法，構建lncRNA-mRNA共表達網絡模塊，結合臨床信息篩選與胃癌分期相關的模塊；采用Kruskal-Wallis秩和檢驗篩選模塊內胃癌不同分期差異表達的lncRNA。結果 獲得286例胃癌組織和30例癌旁對照組織樣本的RNA-seq數據；eQTL分析得到5118對順式作用和1 953 109對反式作用lncRNA-mRNA；2 999個lncRNA和3 884個mRNA納入WGCNA，產生25個共表達模塊，其中與胃癌分期高度相關的模塊有3個；模塊midnightblue、orange內18個樞紐lncRNA中有14個lncRNA在胃癌不同分期差異表達。結論 本研究篩選出14個與胃癌分期相關的lncRNA，這14個lncRNA可能通過調控mRNA的表達影響胃癌的進展，分析其對應的網絡調控關系，為研究lncRNA-mRNA調控機制及探索胃癌治療靶點提供了參考和方向。

lncRNA(long non-coding RNA)是一類轉錄長度超過200nt的非編碼序列，通常由RNA聚合酶II轉錄形成。由于不編碼蛋白質，lncRNA起初被認為是沒有功能的轉錄垃圾。近年來發現，lncRNA在表觀遺傳學、轉錄及轉錄后水平等方面調控基因的表達，如干擾鄰近基因的表達、作為共因子調節轉錄因子的活性、在轉錄后水平可與mRNA形成雙鏈復合物調控基因表達水平。更值得關注的是，lncRNA在細胞異常調節、誘導腫瘤發生過程中起關鍵作用[1]。多項研究表明，lncRNA 與胃癌的發生、發展、轉移和預后密切相關[2-3]。因此，探索與胃癌分期相關的lncRNA，明確其對靶基因的調控關系對胃癌機制研究及治療具有重要意義。

加權基因共表達網絡分析(weighted gene co-expression network analysis，WGCNA)作為一種高通量數據挖掘算法，已廣泛應用于篩選疾病相關基因的研究領域，如利用WGCNA構建與胃癌進展相關的lncRNA網絡[4]。本研究在此基礎上將表達數量性狀位點(expression quantitative trait loci，eQTL)分析與WGCNA相結合，收集胃癌RNA-seq數據及臨床信息數據，利用eQTL篩選出的相關程度高的lncRNA-mRNA共表達對，對應的lncRNA和mRNA納入WGCNA，構建共表達網絡模塊，獲得與胃癌分期相關的lncRNA并分析其調控關系，為胃癌進展及治療靶點研究奠定基礎。

數據與方法

1.數據來源

在TCGA(the cancer genome atlas)數據庫中，納入286例胃癌組織和30例癌旁對照組織樣本的RNA-seq reads count數據。收集匹配的臨床數據，包括年齡、性別、腫瘤TNM分期等信息。

2.分析方法

(1)預處理方法

RNA-seq數據的預處理：包括提取樣本信息、構建基因表達矩陣、將探針名轉化為基因名，最終獲得行名為樣本名、列名為基因名的矩陣；利用edge R進行差異分析，設置log2 fold change界值為1，FDR界值為0.05。

(2)表達數量性狀位點分析

表達數量性狀位點分析由數量性狀位點分析(quantitative trait loci，QTL)發展而來，QTL指定位控制數量性狀的基因在基因組中的位置[5]，eQTL是將每個基因的表達水平作為數量性狀進行分析[6]。通過全基因組mRNA表達量測序得到特定組織樣本的基因表達量，以全部SNP為自變量、每種mRNA表達量為因變量進行線性回歸，得到每一個SNP位點和mRNA表達量之間的關系。由于eQTL分析中表型數據不限于離散型數據，因此eQTL可作為lncRNA-mRNA共表達分析方法[7]。本研究將SNP數據替換為lncRNA表達量數據，分析mRNA與lncRNA表達量的相關性，獲得mRNA-lncRNA共表達對。

eQTL分析的優點是可以區分順式作用關系和反式作用關系，有利于具體分析lncRNA對mRNA的調控作用。順式調控[8]是指lncRNA對染色體臨近位置(小于1Mb)的mRNA的表達調控，位于編碼蛋白上下游的lncRNA可能與啟動子或共表達基因的其他順式作用元件有交集，從而在轉錄或者轉錄后水平對基因的表達進行調控。反式調控是指對不同染色體或者染色體遠端位置的mRNA進行表達調控，調控關系不受空間距離的限定，因此trans-eQTL數量龐大，且可能存在假陽性相關。基于此，本研究選取eQTL分析獲得的全部順式作用共表達對和前5%具有顯著性的反式作用共表達進行WGCNA分析，構建lncRNA與mRNA共表達網絡。

(3)加權基因共表達網絡分析

加權基因共表達網絡分析是一種在高通量基因表達數據中，利用系統生物學思想，計算基因表達相關性，構建基因共表達模塊，進而發現具有生物學意義模塊的高通量數據挖掘算法[9]。該算法首先計算各基因間Pearson相關系數，構建加權鄰接矩陣：

aij=power(Sij，β)≡|Sij|β

其中aij代表基因i與基因j的鄰接系數；Sij代表基因i與基因j的Pearson相關系數。β為軟閾值，軟閾值的確定應滿足使共表達網絡服從無尺度網絡分布，即出現連接度為k的節點的對數lgk與該節點出現的概率的對數lg[P(k)]呈負相關，且R2應>0.8。

基因間的鄰接矩陣轉換為拓撲矩陣(topological matrix，TOM)，將某個基因與所有直接、間接相關的基因間的關系聯系起來，基于拓撲重疊性計算基因與基因間的相異度，根據相異度利用動態剪切樹法劃分基因的功能模塊，完成共表達網絡的基礎構建。繼而計算模塊特征值(即該模塊的第一主成分)與臨床表型信息的相關系數，篩選出與胃癌分期相關的模塊作為關鍵模塊，分析關鍵模塊內的lncRNA及其共表達的mRNA。根據網絡模塊節點中度的分析，確定模塊中的樞紐lncRNA，分析其在胃癌分期是否差異表達。由于目前已知功能的lncRNA微乎其微，這種分析策略有效縮小候選lncRNA的范圍，有助于尋找與胃癌分期相關的lncRNA。

(4)統計分析

所有分析通過R 3.6.1和Cytoscape 3.7.2實現。①edgeR包篩選差異表達基因；②篩選出的基因通過MatrixEQTL包進行eQTL分析，獲得lncRNA-mRNA共表達對；③利用WGCNA包實現共表達網絡的構建并篩選與胃癌分期相關的關鍵模塊；④將模塊內基因共表達網絡的權重信息導入Cytoscape 3.7.2軟件，篩選模塊樞紐基因并繪制基因共表達網絡圖；⑤對樞紐lncRNA進行Kruskal-Wallis秩和檢驗，分析各樞紐lncRNA在不同胃癌分期是否差異表達；⑥模塊內基因映射至在線網站DAVID(http：//david-d.ncifcrf.gov/)中，進行GO和KEGG富集分析。

結 果

1.預處理結果

差異分析獲得4 767個差異表達的mRNA，其中2 466個mRNA上調，2 301個mRNA下調；3 542個差異表達的lncRNA，其中2 767個lncRNA上調，775個lncRNA下調。對差異基因進行火山圖的可視化，見圖1。

圖1 差異表達基因火山圖

2.構建胃癌lncRNA-mRNA共表達網絡

以年齡、性別為協變量對差異基因進行eQTL分析，獲得5 118對順式作用共表達對和1 953 109對反式作用共表達對，取全部的順式作用共表達對和前5%具有顯著性的反式作用共表達對，去除重復基因，最終獲得2 999個lncRNA和3 884個mRNA進入下一步的WGCNA分析。

選取軟閾值β=4，通過動態剪切樹法進行模塊初步識別并合并相似模塊，設置每個基因網絡模塊最少的基因數目為50，模塊合并閾值為0.3，最終得到25個基因網絡模塊，灰色模塊是無法聚類到其他任何模塊的基因集合，見圖2。

圖2 加權基因共表達網絡模塊聚類圖

根據各模塊的特征向量分別計算每個模塊與胃癌分期的相關性，繪制模塊樣本性狀相關性熱圖，見圖3。結果顯示，midnightblue、orange、yellow、purple、blue、royalblue、green共七個模塊與胃癌分期的相關性有統計學意義(P<0.05)。

圖3 臨床信息與模塊相關性熱圖

3.篩選與胃癌分期相關lncRNA

選擇與胃癌分期相關性較強的前三個模塊，midnightblue、orange、yellow模塊，各模塊按度的大小排序前50的基因為樞紐基因，見圖4。三個模塊中分別包含10個、8個、3個樞紐lncRNA和40個、42個、47個樞紐mRNA。秩和檢驗結果獲得14個與胃癌分期相關lncRNA：midnightblue模塊有AP002954.4、AC002331.1、LINC01272、RP11-44K6.4、RP11-638I2.9、AC069363.1、LINC01094，orange模塊有RP11-443C10.1、RP11-13P5.2、AC093850.2、RP11-576I22.2、AC007750.5、RP11-95H3.1、RP11-867G23.1；其中AP002954.4在T3、T4期表達水平相較于T1期上調，AC002331.1、RP11-44K6.4在T4期表達水平相較于T1、T2期上調，AC069363.1在T4期表達水平相較于T1期上調，RP11-638I2.9在T4期表達水平相較于T2期上調，LINC01272在T2期表達水平相較T1期上調，LINC01094、RP11-443C10.1、RP11-13P5.2、AC093850.2、RP11-576I22.2、AC007750.5、RP11-95H3.1、RP11-867G23.1在T2、T3、T4表達水平相較于T1期上調，差異均具有統計學意義，見表1。與14個lncRNA共表達的mRNA共有59個，其中有56個mRNA秩和檢驗有統計學意義。

表1 樞紐lncRNA秩和檢驗結果

圖4 midnightblue(a)、orange(b)、yellow(c)模塊的共表達網絡圖

4.富集分析

lncRNA-mRNA共表達模塊中的59個mRNA進行GO和KEGG分析，結果見表2。GO分析結果顯示，模塊基因主要涉及到免疫應答(immune response)、防御反應(defense response)、趨化作用(chemotaxis)和炎性反應(inflammatory response)等生物過程；KEGG分析結果顯示，模塊基因主要富集于細胞因子-受體相互作用通路(cytokine-cytokine receptor interaction)、Toll樣受體信號通路(toll-like receptor signaling pathway)、趨化因子信號通路(chemokine signaling pathway)、黏著斑通路(focal adhesion)等，這些通路與機體炎性反應有關。

表2 富集分析結果

討 論

本研究首次將WGCNA與eQTL相結合應用在癌癥數據中，篩選與胃癌分期相關的lncRNA。由于目前大部分lncRNA功能未知，通過分析lncRNA與mRNA表達相關性，利用已知功能的mRNA推導lncRNA的功能成為一種重要的研究策略。相較于既有的單純通過WGCNA構建與胃癌進展相關的lncRNA網絡的研究[4]，本研究結合eQTL的優勢在于，一方面篩選出相關程度較高的lncRNA與mRNA進行共表達網絡分析，提高lncRNA功能推斷的準確性；另一方面將lncRNA與mRNA的共表達區分為順式共表達對與反式共表達對，有利于后續進一步研究生物學上lncRNA與mRNA的調控關系。本研究獲得14個與胃癌分期相關的lncRNA，總體上隨著胃癌分期的進展，lncRNA表達水平呈上升趨勢。與14個lncRNA共表達的59個mRNA中56個mRNA在不同分期表達水平差異有統計學意義，這說明eQTL與WGCNA相結合的方法能夠有效篩選出表達相關性高的lncRNA-mRNA，篩選出的lncRNA可能通過調控mRNA的表達在胃癌的發展進程中具有關鍵作用，可作為潛在的生物標志物。研究發現，midnightblue模塊的樞紐lncRNA LINC01272在肺鱗癌與癌旁組織中顯著差異表達，且可以作為診斷肺鱗癌早期與晚期的標志物[10]；orange模塊的樞紐lncRNA AC093850.2在肺鱗癌、肺腺癌、乳腺癌及胰腺癌中表達水平均高于癌旁組織[11]。目前仍未有研究發現這些樞紐lncRNA與胃癌分期有關，本研究結果為深入研究胃癌分期相關lncRNA提供了參考和依據。

lncRNA通過調控mRNA的轉錄、穩定性等多種方式影響靶基因的表達，lncRNA-mRNA共表達網絡是研究lncRNA功能和調控機制的重要方式。GO分析和KEGG分析結果顯示lncRNA-mRNA共表達網絡中的mRNA主要參與免疫應答、防御反應、炎性反應等生物過程，可以推斷這些樞紐lncRNA可能通過調控mRNA影響機體的炎性反應、免疫應答等，在一定程度影響胃癌的進展過程。已有多項研究發現，非編碼RNA可以通過影響腫瘤和免疫細胞中免疫調節分子的表達影響抗腫瘤免疫應答、參與炎性反應過程[12]。midnightblue模塊中的樞紐lncRNA LINC01272已被證實與炎性反應有關，Wang Sen等的研究發現LINC01272在炎癥性腸病患者中的表達水平顯著高于健康對照，是炎癥性腸病的潛在診斷標志物[13]。

綜上，利用eQTL和WGCNA相結合的方法構建lncRNA-mRNA共表達網絡，探索與胃癌分期相關的lncRNA及其調控關系。篩選出的14個與胃癌分期相關的lncRNA可為研究胃癌的進展機制提供依據，其對應的網絡調控關系為后續更深入研究lncRNA-mRNA調控機制及探索胃癌治療靶點提供了參考和方向。