基于生物信息學分析肝細胞癌中lncRNA-miRNA-mRNA的調控網絡

張 政，盧鴻健，李明慧，高文博，田 宇，王嘉祺，張榮花，章廣玲

(1.華北理工大學臨床醫學院/河北省醫工融合精準醫療重點實驗室，河北 唐山 063000；2.華北理工大學基礎醫學院/河北省慢性疾病基礎醫學重點實驗室，河北 唐山 063210)

肝細胞癌（hepatocellular carcinoma，HCC）是世界上常見的惡性腫瘤之一，占原發性肝癌的75%～85%。GLOBOCAN 的最新統計數據顯示[1]，2020 年全球HCC 新發病例達91 萬，死亡病例83 萬，分別居腫瘤新發病率第6 位、死亡率第3 位；我國人口占世界的18%，但新發病例占全球的45%。HCC 的發生發展過程復雜，是一種異質性較高的癌癥，大多數發現時已是中晚期，錯過了手術的最佳時間且預后較差，其5 年總生存率只有12%[2]。因此，尋找HCC 發展過程中有效的靶點和信號通路并進行相應干預，對延緩HCC 的進展或者逆轉具有重要意義。研究表明[3]，非編碼RNA（non-coding RNA，ncRNA）包括長鏈非編碼RNA（long non-coding RNA，lncRNA）和微小RNA（microRNA，miRNA）等，其可能在HCC 的發生發展中發揮了重要的作用。miRNAs 分子由18～25個bp 組成的非編碼RNA，能夠通過堿基互補配對與靶基因的3' 端非編碼區（3'untranslated region，3'UTR）特異性結合，沉默或降解靶mRNAs 從而影響其翻譯功能。lncRNAs 分子是一種長度超過200 bp 的非編碼核苷酸，含有miRNAs 的結合位點，可以作為miRNAs 的“海綿”吸附miRNAs，從而競爭性抑制miRNAs 對靶mRNAs 的作用[4]。本研究擬通過生物信息學方法探究癌癥基因組圖譜（TCGA）數據庫中HCC 組織中差異表達的基因，并構建lncRNAmiRNA-mRNA 調控網絡，為靶向診斷和治療HCC提供科學依據。

1 資料與方法

1.1 HCC 組織中差異表達基因的篩選 利用R 軟件DESeq2 包分析TCGA 數據庫中374 例HCC 組織與50 例正常肝組織中lncRNAs、miRNAs、mRNAs 的表達譜。使用R 軟件ggplot2 包對表達譜數據繪制火山圖，以|Log2FC|＞1（fold change 差異倍數大于2 倍）且P＜0.05 作為篩選標準。

1.2 miRNAs 轉錄因子的預測 利用FunRich 軟件預測差異表達miRNAs 的轉錄因子，根據其所調控的基因數和值進行排序，選取富集程度最顯著的前10 個轉錄因子進行展示。

1.3 基因功能（GO）和京都基因與基因組百科全書（KEGG）富集分析 利用R 軟件clusterProfiler 包分別分析miRNAs、mRNAs 各基因之間潛在的影響，并以調整后的P＜0.05 為標準，進行GO 和KEGG分析，分別預測其主要生物功能和參與的主要信號通路。

1.4 miRNA-mRNA 調控網絡 使用FunRich 軟件預測miRNA 的候選靶基因，結合HCC 組織中差異表達的miRNAs 和mRNAs 的FC 值，篩選具有相互作用的miRNA-mRNA 關系。

1.5 PPI 網絡構建 將篩選出的基因導入STRING 數據庫，分析蛋白之間的相互作用（protein-protein interaction，PPI），并將分析結果導入Cytoscape 中，構建可視化PPI 蛋白互作網絡圖。

1.6 lncRNA-miRNA-mRNA 調控網絡構建 利用Starbase 數據庫，預測miRNAs 的上游lncRNAs，結合HCC 中差異表達lncRNAs 的FC 值，篩選具有相互作用的lncRNA-miRNA 關系對。基于共享的miRNAs 通過Cytoscape 軟件對lncRNA-miRNAmRNA 調控網絡進行可視化。

2 結果

2.1 差異表達lncRNAs、miRNAs、mRNAs 的篩選 利用R 軟件DESeq2 包分析374 例HCC 組織和50 例正常肝 組織中lncRNAs、miRNAs、mRNAs 的表達譜，根據P＜0.05 和|Log2FC|＞1 的標準篩選篩選出HCC 組織中差異表達的lncRNAs 2850 個，其中上調2475 個，下調375 個；差異表達的miRNAs 38個，其中上調31 個，下調7 個；差異表達的mRNAs 4455 個，其中上調3199 個，下調1256 個，見圖1。

圖1 差異表達基因的火山圖

2.2 轉錄因子-miRNA 調控網絡 利用FunRich 軟件對38 個差異表達miRNAs 進行轉錄因子預測，共得到201 個轉錄因子，根據其所調控的基因數和P值，選取富集程度顯著的前10 個轉錄因子進行展示，見圖2，其基本信息見表1，前10 個轉錄因子分別為EGR1、SP1、NKX6-1、POU2F1、LHX3、SP4、MEF2A、ARID3A、MSX2、HOXA9。

表1 差異表達前10 名轉錄因子的基本信息

圖2 差異表達miRNAs 轉錄因子的預測

2.3 miRNAs 的GO 和KEGG 富集分析 對38 個差異表達的miRNAs 進行GO 和KEGG 功能富集分析，其生物過程主要涉及血管內皮細胞遷移的調節、血管內皮細胞增殖調控與新生血管生成、內皮細胞增殖的調節、血管內皮細胞增殖與新生血管生成；分子功能主要涉及mRNA 結合參與轉錄后基因沉默、RNA 聚合酶Ⅱ復合物結合、RNA 聚合酶核心酶結合、基礎轉錄機制結合；KEGG 分析主要富集于miRNAs 在癌癥中的作用，見圖3。

圖3 差異表達miRNAs 的GO 和KEGG 功能富集分析

2.4 miRNA-mRNA 調控網絡 通過FunRich 軟件預測差異表達的38 個miRNAs 的下游靶基因，共得到1530 個靶基因，然后與差異表達的4455 個mRNAs取交集獲得286 個候選靶基因，見圖4。通過Starbase 和MiRDB 數據庫比對分析，獲得14 302 對miRNA-mRNA，與差異表達的miRNA-mRNA 匹配后，最后篩選出交叉表達負調控的12 個miRNAs 和108 個mRNAs 作為后續研究對象。

圖4 差異表達miRNAs 的下游靶基因與差異表達mRNAs 的韋恩圖

2.5 108 個mRNAs 的蛋白PPI 調控網絡及富集分析利用STRING 數據庫和Cytoscape 軟件構建篩選出的108 個基因的PPI 網絡，見圖5。圖中節點代表mRNAs，有相互作用關系的mRNAs 之間有連線，共有44 個節點和134 條邊；節點越大、顏色越深則介數中心度越大，邊的粗細代表連接評分，邊越粗，mRNAs 間的互作關系越大。

圖5 108 個mRNAs 的蛋白互作網絡圖

2.6 108 個mRNAs 的GO 與KEGG 富集分析 利用GO 功能分析篩選出的108 個mRNAs 的生物功能，生物過程主要涉及肌肉適應、MAP 激酶活性的調節、胚胎器官發育、蛋白絲氨酸/蘇氨酸激酶活性的正調控；分子功能主要涉及MAP 激酶酪氨酸/絲氨酸/蘇氨酸磷酸酶活性、MAP 激酶磷酸酶活性、激酶調節活性、DNA 結合轉錄激活物活性和RNA 聚合酶Ⅱ特異性。KEGG 富集分析108 個mRNAs 參與的信號通路，主要涉及調節干細胞多能性的信號通路、催乳素信號通路、MAPK 信號通路、EGFR 酪氨酸激酶抑制劑耐藥性，見圖6。

圖6 108 個mRNAs 的GO 和KEGG 富集分析圖

2.7 構建lncRNA-miRNA-mRNA 調控網絡 利用Starbase 數據庫，分析篩選出的12 個miRNAs 的上游lncRNAs 共221 個，結合差異表達lncRNAs 的FC 值，篩選出5 對lncRNA-miRNA，包括2 個miRNAs（miR-23A 和miR-146B）與5 個lncRNAs（GAS5、MCM3AP-AS1、LINC00173、KCNQ1OT1、SNHG7）；結合篩選出的2 個miRNAs 的與組織中差異表達mRNAs 的FC 值，共匹配42 對miRNAmRNA（包括40 個mRNAs）。通過Cytoscape 軟件對lncRNA-miRNA-mRNA 調控網絡進行可視化處理，見圖7。

圖7 HCC 組織中lncRNA-miRNA-mRNA 的調控網絡

3 討論

HCC 是由肝細胞惡性病變引起的癌癥，患者發現時大多數已是中晚期，具有發病隱匿、進展快、易耐藥、預后差以及易復發等特點，且現有的治療效果并不理想[5]。因此尋找開發有效的HCC 診斷和預后的生物標記物具有重要意義。越來越多的研究表明，lncRNAs 和miRNAs 等非編碼RNAs 與HCC 的發病密切相關。LncRNA n339260 通過降低miR-30e-5p的表達使TP53INP1 的表達增加，從而促進肝細胞癌的進展[6]。LncRNA ZEB1-AS1 通過調節miR-1224-5p/MAP4K4 軸，調節持續沙眼衣原體感染中線粒體介導的HeLa 細胞凋亡[7]。抑制lncRNA GAS6-AS2 通過調節miR-136-5p/OXSR1 軸能夠在體外和體內水平減輕敗血癥相關的急性腎損傷[8]。然而在HCC 中尚沒有報道完整的lncRNA-miRNA-mRNA 調控網絡，因此深入研究lncRNAs 和miRNAs 在HCC 發生發展中的作用及調控機制具有重要意義。

已有研究表明，轉錄因子調控miRNAs 的表達，且這些轉錄因子在HCC 的發病過程中具有重要的調控作用。本研究預測了調控HCC 中差異表達miRNAs 的上游轉錄因子共201 個，根據值篩選出差異顯著的前10 個轉錄因子分別為EGR1、SP1、NKX6-1、POU2F1、LHX3、SP4、MEF2A、ARID3A、MSX2、HOXA9，其中EGR1 差異最明顯。EGR1 是一種轉錄因子，可能參與腫瘤細胞的增殖、侵襲和轉移以及腫瘤血管生成，與癌癥的發生和發展密切相關[9]。有報道表明[10]，EGR1 可能以轉錄因子方式調控hsa-miR-146b 的表達。EGR1 在miR-195 的上游啟動子區有結合位點，且調控miR-195 的表達，EGR1 可抑制miR-195/AKT3 軸并抑制胃癌細胞凋亡[11]。EGR1 還能夠通過負調控miR-491-5p 促進細胞增殖、生長和遷移，并抑制細胞凋亡，從而抑制胃癌的進程[12]。因此，研究轉錄因子對miRNAs 的調控作用可在系統層次上對疾病的發生提供理論依據。

近年來，伴隨著高通量測序技術的發展產生了大量的基因數據，這為進一步研究HCC 的發生發展及預后監測提供了有力的基礎。本研究首先分析TCGA 數據庫中差異表達基因的表達譜，得到的差異表達的lncRNAs 2850 個、miRNAs 38 個、mRNAs 4455 個。對38 個差異表達的miRNAs 進行GO 和KEGG 功能富集分析，顯示主要涉及血管內皮細胞遷移和增殖、新生血管生成、內皮細胞增殖的調節等生物過程，mRNA 結合參與轉錄后基因沉默等分子功能，且KEGG 分析顯示這些miRNAs 主要在癌癥中發揮作用。通過FunRich 軟件預測得到的miRNAs 靶基因與差異表達的mRNAs 交叉匹配，對篩選出的108 個mRNAs 進行GO 分析發現，主要在胚胎器官發育、MAP 激酶活性的調節等生物學過程富集，DNA 結合轉錄激活物活性和RNA 聚合酶Ⅱ特異性、激酶調節活性等分子功能富集。KEGG 分析發現，靶基因參與的信號通路主要涉及調節干細胞多能性的信號通路、MAPK 信號通路等。因此可以得出在HCC 中差異表達的mRNAs 其生物功能主要與細胞增殖、形態分化和細胞遷移有關。同樣地，在PPI 網絡的Top 10 基因中KLF4、NRG1、FOSB、MAP2K1 也參與了上述生物學功能以及相關信號通路。研究發現[13]，miR-146B 通過靶向KLF4 調節肝星狀細胞活化。NRG1 是鱗狀細胞癌分化的關鍵調節因子，抑制腫瘤細胞角化和終末鱗狀分化并促進腫瘤細胞增殖[14]。TGF-β1通過誘導FOSB 促進前列腺癌細胞遷移和侵襲[15]。CircRNA ZFR 通過上調MAP2K1 促進HCC 的增殖[16]。

本研究尋找與miRNAs 相匹配的上游lncRNAs和下游mRNAs 構建HCC 中lncRNA-miRNA-mRNA 的調控網絡，其中包含2 個關鍵的miRNAs（miR-23A、miR-146B）以及與其匹配的5 個lncRNAs（GAS5、MCM3AP-AS1、LINC00173、KCNQ1OT1、SNHG7）和40 個mRNAs（ROBO1、LIN28B 等）。研究表明[17]，miR-23A 是HCC 中下調最顯著的miRNA，可能是改善肝癌患者索拉非尼反應性的潛在靶點。miR-146B 通過靶定PTP1B 抑制胃癌細胞的增殖并促進其凋亡[18]。miR-146B 通過NF-κB（核因子κB）途徑靶向TLR4（Toll 樣受體4），抑制膽囊癌腫瘤的發生和轉移[19]。以上研究表明，miR-23A、miR-146B具有抑制腫瘤的作用。此外，lncRNA GAS5 靶向miR-23A 抑制CCl4 誘導的肝纖維化[20]。LncRNA SNHG7 通過調節miR-146B-5p/ROBO1 軸促進胰腺癌進展[21]。LncRNA SNHG17 通過抑制miR-23A-3p 來調節CXCL12 介導的血管生成，從而促進結直腸腺癌細胞的增殖和遷移[22]。LncRNA KCNQ1OT1通過miR-146A-5p/ACER3 軸抑制HCC 的放射敏感性并促進腫瘤發生[23]。

綜上所述，lncRNA-miRNA-mRNA 調控網絡的構建，有望闡明HCC 的發生發展機制，為HCC 的臨床診斷和預后監測提供有價值的標記物。但是lncRNA-miRNA-mRNA 調控網絡的功能以及各差異基因的功能，還有待通過體內外實驗和臨床實驗進一步證實。