基于生物信息學構建前列腺癌單細胞層面lncRNA-miRNA-mRNA調控網絡*

莫泳鋒，謝遠亮，陳蔓瑩，覃洪宇，葉 雨**

(1.廣西醫科大學第二附屬醫院，廣西南寧 530007；2.廣西醫科大學附屬腫瘤醫院，廣西南寧 530021)

前列腺癌(Prostate Cancer，PCa)是常見的泌尿系統惡性腫瘤之一。據統計，2018年全世界有將近130萬PCa新發病例和35.9萬PCa死亡病例。PCa的發病率為13.5%，死亡率為6.7%，分別居全球男性腫瘤疾病第二位和第五位[1]。隨著社會的發展，人口逐漸向老齡化過渡，飲食健康問題及生活方式的改變使PCa的發病形勢仍十分嚴峻[2,3]。因此，開展PCa的機制研究，不斷深入探索腫瘤的發病機制，從分子機制上取得進一步突破,對臨床診治工作有必要且關鍵的作用。

生物信息學分析是基于高通量測序結果研究腫瘤發病機制和鑒定預后生物標志的重要方法。lncRNA是內源性非編碼長RNA(長度200-100 000個核苷酸)，轉錄本主要位于細胞核中，其穩定性和保守性較差，但特異性強[4]。lncRNA參與包括PCa在內的多種類型癌癥的細胞增殖、遷移、侵襲和凋亡等生物學過程[5-8]。miRNA是內源性非編碼微小RNA(長度約22個核苷酸)，在細胞和體液(包括血漿或血清)中相對穩定[9,10]。在大多數情況下，miRNA與靶信使RNA的3′未翻譯區相互作用，在轉錄后誘導mRNA降解和抑制翻譯[11]。有研究提出了競爭性內源性RNA(Competing Endogenous RNA，ceRNA)，即編碼RNA和非編碼RNA可以通過競爭相同的miRNA反應元件(miRNA Response Elements，MREs)來相互調節彼此的表達水平[12-14]。ceRNA的提出意味著在lncRNA-mRNA之間存在著通信調節網絡，通過競爭相同的miRNA序列來調節彼此的表達，形成一個復雜多樣的lncRNA-miRNA-mRNA調控網絡從而調節生物學過程和腫瘤發生，如轉錄調控、轉錄后調控、蛋白修飾等[15]。近年來雖然對ceRNA的研究越來越深入并取得了一定的成果，但是lncRNA-miRNA-mRNA網絡對PCa的具體調控機制還不明確。

單細胞RNA測序(Single cell RNA sequencing，scRNA-seq)是一種可以在單細胞水平上探索復雜組織細胞內轉錄組信息的新技術[16]。與批量RNA測序(Bulk RNA-seq)相比，scRNA-seq可以反映單個細胞的轉錄組表達水平，為細胞分子機制研究提供更精細且可靠的研究方法。lncRNA-miRNA-mRNA網絡的調控機制在bulk RNA-seq的支持下已得到越來越多學者的認可，并證實其在腫瘤的發生發展中起著重要作用。目前，關于scRNA-seq與bulk RNA-seq技術結合構建的lncRNA-miRNA-mRNA網絡研究仍較少，將二者結合構建單細胞水平的lncRNA-miRNA-mRNA調控網絡可為腫瘤致病機制研究提供新方向，同時有望獲得潛在的預后生物標志物，為PCa的診斷和治療提供新思路。

1 資料與方法

1.1 數據獲取

從GDC(https://portal.gdc.cancer.gov/)中獲取TCGA數據庫PCa的RNA-Seq測序文件，共獲得499例腫瘤組織和52例正常前列腺組織。通過GEO數據庫(https://www.ncbi.nlm.nih.gov/geo/)獲取GSE157703中的2例前列腺癌單細胞測序數據。

1.2 構建表達矩陣

通過R_4.1.1軟件(https://www.r-project.org/)，使用rtracklayer_1.50.0和SummarizedExperiment_1.20.0命令集構建人類基因信息比對數據庫。用data.table_1.14.0、tidyr_1.1.3和dplyr命令集讀取PCa的RNA-Seq測序文件進行整理獲取全轉錄組表達矩陣，然后與基因信息庫進行基因匹配，獲取轉錄組的基因屬性，通過基因屬性分別獲取lncRNA和mRNA的表達矩陣。

1.3 差異基因的篩選

用R語言的DESeq2_1.30.1和edgeR_3.32.1命令集對構建的lncRNA和mRNA表達矩陣進行整理分析，分別得到lncRNA和mRNA差異表達基因(Differentially Expressed Genes，DEGs)，使用plot命令集繪制靶基因的mRNA火山圖，用pheatmap_1.0.12命令集繪制Top 100-DEGs的熱圖。DEGs的篩選條件為|log2FC|>1且adj.P<0.05，log2FC表示取差異表達倍數以2為底的對數，adj.P表示矯正的P值。

1.4 lncRNA-miRNA-mRNA網絡的構建

通過Mircode數據庫(http://www.mircode.org/)，運用perl軟件對lncRNA-miRNA關系對進行預測。將獲取的miRNA通過TargetScan數據庫(https://www.targetscan.org/vert_72/)、miRTarBase數據庫(https://mirtarbase.cuhk.edu.cn/)和miRDB數據庫(http://mirdb.org/)分別預測miRNA的靶基因，取交集得到最終的miRNA靶基因，然后與整合的DEGs取交集，最終得到多條lncRNA-miRNA-mRNA軸，據此構建一個lncRNA-miRNA-mRNA調控網絡，并通過Cytoscape_3.8.1進行可視化。

1.5 功能富集分析

對構建的lncRNA-miRNA-mRNA網絡中的靶基因進行生物學功能及通路富集分析。首先對lncRNA-miRNA-mRNA網絡中的mRNA運用R語言的org.Hs.eg.db_3.12.0包進行ID轉換；然后通過clusterProfiler_3.18.1、enrichplot_1.10.2和ggplot2_3.3.3包進行基因本體論(Gene Ontology，GO)以及京都基因與基因組百科全書(Kyoto Encyclopedia of Genes and Genomes，KEGG)分析，以PvalueCutoff=0.05、QvalueCutoff=0.05為篩選條件，結果以氣泡圖和圈圖顯示；最后使用STRING (https://string-db.org/)數據庫對網絡中的mRNA進行蛋白互作分析，設定Medium confidence>0.4，結果使用Cytoscape_3.8.1將數據可視化。

1.6 單細胞測序數據的整合以及構建單細胞水平調控網絡

對獲取的單細胞測序文件，首先通過R語言的Seurat_4.0.1、ggplot2_3.3.3、clustree_0.4.3和cowplot_1.1.1包等命令集，合并整理構建單細胞表達矩陣、矯正批次效應以及質量控制；然后進行降維，根據Top 2 000差異基因進行聚類和分群；最后根據mark基因對各細胞亞群進行定義，通過拷貝數變異(Copy Number Variation,CNV)定義腫瘤細胞亞群，并進一步驗證lncRNA-miRNA-mRNA網絡中的靶基因在各細胞亞群上的表達情況，同時構建單細胞水平的lncRNA-miRNA-mRNA調控網絡。

1.7 預后分析與免疫組織化學驗證

利用R軟件的cgdsr、DT、survival等數據包，獲取494例PCa的臨床數據及轉錄組數據，根據基因的中位表達情況，將其分為目的基因高表達組和低表達組，并分析無病生存期(Disease-Free Survival，DFS)，根據疾病有無復發或死亡將其分為復發組和死亡組或無復發組和死亡組，并分析基因的表達情況，通過ggpubr包進行可視化。然后將基因輸入在線數據庫人類蛋白質圖譜(https://www.proteinatlas.org/)，參數選擇前列腺癌病理圖譜，獲取目的基因免疫組織化學的鏡下圖像并導出，通過圖像處理軟件Adobe Illustrator可視化。

2 結果與分析

2.1 差異表達基因的篩選與鑒定

通過獲取TCGA的PCa測序數據，經過篩選得到3 013個DEmRNA，其中上調1 282個，下調1 731個[圖1(a)]。由圖1(b)可知，共得到1 101個DElncRNA，其中上調443個，下調658個。取前100個差異基因作熱圖分析，結果顯示差異基因在腫瘤和正常組織中存在差異表達(圖2)。

圖2 差異基因熱圖分析結果Fig.2 Heat map analysis results of differential gene

2.2 lncRNA-miRNA-mRNA網絡的構建與分析

將得到的DElncRNA通過Mircode數據庫進行miRNA預測，得到1 935個lncRNA-miRNA關系對，然后將miRNA通過TargetScan、miRTarBase和miRDB數據庫預測mRNA，與mRNA DEGs取交集得到97個miRNA的靶基因，最終納入構建調控網絡的基因共175個，包括51個lncRNA、27個miRNA和97個mRNA。結果通過調控網絡圖進行可視化(圖3)。

Blue diamond represents lncRNA,red triangle represents miRNA,green oval represents mRNA,and the linear connection represents that they have regulatory relationships

2.3 GO生物功能富集與KEGG信號通路富集分析

根據GO分析和KEGG分析，以P<0.05，FDR<0.05為篩選條件，并按照富集基因個數從大到小排列，將排名在前15的富集結果和KEGG分析結果用氣泡圖和圈圖表示(圖4)。GO分析結果顯示，lncRNA-miRNA-mRNA調控網絡的靶基因在腺狀細胞遷移、化學突觸傳遞的調節、跨突觸信號的調節、細胞形態分化與發生的調節、上皮細胞遷移和組織遷移等生物學過程富集[圖4：(a)(c)(e)]。KEGG分析結果顯示，miRNA在癌癥、軸突導向、胞間的黏附、肌動蛋白細胞骨架調節、Rap1、MAPK、PI3K-Akt和Ras等信號通路富集[圖4:(b)(d)(f)]。

圖4 GO和KEGG分析結果Fig.4 Analysis results of GO and KEGG

2.4 蛋白質互作網絡的構建與分析

使用STRING數據庫分析lncRNA-miRNA-mRNA調控網絡靶基因的蛋白互作關系，通過Cytoscape軟件繪制，得到蛋白-蛋白互作(PPI)網絡圖，圓圈越大則互作對越多。由圖5可知，共得到58個蛋白質的互作網絡，其中Top 10基因分別為BDNF、FGF2、EZH2、MET、KIT、PIK3R1、PTGS2、RRM2、DLGAP5和CEP55。

Circle represents proteins,the larger the circle,the higher the interaction score;and the linear connection represents the interaction between proteins

2.5 單細胞測序數據的整合及構建單細胞水平調控網絡

使用R軟件將獲取的GSE157703中的2例前列腺癌單細胞測序數據進行整合分析，得到6 803個細胞、24 639個基因的單細胞矩陣，最終定義了B細胞、T細胞、成纖維細胞、肥大細胞、內皮細胞、上皮細胞、平滑肌細胞、髓系細胞和腫瘤細胞共9個細胞亞群(圖6)。

Sample 1 and Sample 2 are sequenced samples of cancer tissue from two different PCa patients

為驗證lncRNA-miRNA-mRNA調控網絡靶基因在各個亞群的表達情況，取平均表達10%以上的亞群作為單細胞水平調控網絡構建的標準[圖7：(a)(b)]，據此可以得到B細胞表達靶基因PEG10，調控網絡中有25個lncRNA、1個miRNA和1個mRNA，共25條lncRNA-miRNA-mRNA軸[圖7(c)]；得到T細胞表達靶基因DUSP2，調控網絡中有9個lncRNA、2個miRNA和1個mRNA，共18條lncRNA-miRNA-mRNA軸[圖7(d)]；得到成纖維細胞表達靶基因PDPN、HAS2、MAP1B和INMT，調控網絡中有32個lncRNA、5個miRNA和4個mRNA，共70條lncRNA-miRNA-mRNA軸[圖7(e)]；得到肥大細胞表達靶基因KIT和PTGS2，調控網絡中有6個lncRNA、1個miRNA和2個mRNA，共12條lncRNA-miRNA-mRNA軸[圖7(f)]；得到內皮細胞表達靶基因GJA1、SNCG、TGFBR3、DUSP5、INMT和MET，調控網絡中有46個lncRNA、6個miRNA和6個mRNA，共158條lncRNA-miRNA-mRNA軸[圖7(g)]；得到平滑肌細胞表達靶基因RBPMS2、SNCG、MAP1B和ARC，調控網絡中有39個lncRNA、5個miRNA和4個mRNA，共77條lncRNA-miRNA-mRNA軸[圖7(h)]；得到上皮細胞表達靶基因EGLN3，調控網絡中有11個lncRNA、3個miRNA和1個mRNA，共23條lncRNA-miRNA-mRNA軸[圖7(i)]；得到髓系細胞表達靶基因SLC7A11和DUSP5，調控網絡中有31個lncRNA、3個miRNA和2個mRNA，共43條lncRNA-miRNA-mRNA軸[圖7(j)]；得到腫瘤細胞表達靶基因TRIM36、TRIM29、TP53INP1、ITGA2、PDLIM5、AMOT、PRRG4和MET等，調控網絡中有45個lncRNA、10個miRNA和8個mRNA，共188條lncRNA-miRNA-mRNA軸[圖7(k)]。

圖7 lncRNA-miRNA-mRNA調控網絡在單細胞亞群中的表達情況Fig.7 Expression of lncRNA-miRNA-mRNA regulatory network in single cell subpopulation

2.6 PCa腫瘤細胞相關靶基因預后分析及驗證

通過R軟件分別獲取ITGA2、MET、PDLIM5、PRRG4以及TP53INP1基因的DFS及其分組表達情況(圖8)。結果顯示，ITGA2、PDLIM5H和PRRG4的DFS存在差異且均具有統計學意義(P<0.05)，上述基因的低表達組有更高的腫瘤復發率或死亡率。為了進一步驗證調控網絡的靶基因在PCa單細胞水平上的腫瘤細胞群的表達情況，通過在線數據庫分別驗證ITGA2、PDLIM5和PRRG4基因在腫瘤細胞中的蛋白水平表達情況，結果顯示上述靶基因在PCa患者的腫瘤組織病理中均有不同程度的表達。

3 討論

隨著測序技術的不斷發展，研究重點逐步從編碼RNA向非編碼RNA轉變。隨著近年來對腫瘤lncRNA-miRNA-mRNA軸的不斷深入研究，越來越多的研究表明lncRNA-miRNA-mRNA與PCa的發生和發展密切相關。體外實驗發現，lncRNA KCNQ1OT1/miR-211-5p/CHI3L1調控PCa細胞的增殖、遷移、侵襲[17]。Luo等[18]的研究發現，lncRNA TTN-AS1/miR-193a-5p在PCa中起到抑制細胞凋亡的作用。Jiang等[19]的研究發現，lncRNA HOXD-AS1/miR-361-5p/FOXM1參與PCa腫瘤的轉移。也有研究提出lncRNA HOTAIR通過競爭性吸附miR-590-5p，參與調節STAT3，并參與PCa細胞對多西他賽的耐藥[20]。上述研究結果表明，lncRNA-miRNA-mRNA調控網絡參與了調節PCa的增殖、遷移、侵襲和凋亡等多個生物學過程以及化療藥物的抵抗。以往構建的lncRNA-miRNA-mRNA調控網絡是基于bulk測序來進行的，具有宏觀的指導意義，但更精確的微觀調控關系仍待進一步研究。單細胞測序技術的提出，使測序結果能精確到單個細胞水平，為腫瘤分子機制的研究提供了新方法。本研究通過對TCGA-PRAD的bulk測序和GEO中GSE157703的scRNA數據進行整合分析，以此構建單細胞水平的lncRNA-miRNA-mRNA調控網絡，這些調控網絡可以為PCa的分子機制研究提供理論依據。

本研究通過對TCGA測序數據進行處理，得到3 013個DEmRNA和1 101個DElncRNA，篩選出51個lncRNA、27個miRNA和97個mRNA來構建lncRNA-miRNA-mRNA調控網絡。通過GO分析發現，調控網絡的靶基因在腺狀細胞的遷移、細胞形態分化與發生調節、上皮細胞遷移和組織遷移等生物學過程富集。KEGG分析發現，調控網絡的靶基因在胞間黏附、肌動蛋白細胞骨架調節等信號通路富集。因此可以得出生物功能主要與細胞形態分化、胞間連接和細胞遷移有關，信號通路參與了細胞黏附、細胞連接的形成。同樣地，在PPI的Top 10基因中BDNF、FGF2、MET和KIT也參與了上述生物功能以及相關信號通路。遠處轉移是PCa患者死亡的重要原因之一，研究PCa轉移的機制對于改善其預后至關重要[21]。研究發現lncRNA HOXD-AS1通過miR-361-5p/FOXM1參與PCa患者的轉移[19]。本研究結果提示lncRNA-miRNA-mRNA網絡參與了PCa的細胞連接、遷移以及黏附等相關生物學功能和通路的調控，可見其在轉移相關分子機制的調控中有重要作用。

根據標記基因，在單細胞測序的聚類定義后，可得到B細胞、T細胞、成纖維細胞、肥大細胞、內皮細胞、上皮細胞、平滑肌細胞、髓系細胞和腫瘤細胞共9個細胞亞群的lncRNA-miRNA-mRNA調控關系。下面將其分為免疫細胞(B細胞、T細胞和髓系細胞)、上皮細胞、間質細胞(成纖維細胞、平滑肌細胞、肥大細胞和內皮細胞)以及腫瘤細胞4大類進行討論。

免疫細胞的lncRNA-miRNA-mRNA調控網絡主要通過調控PEG10、SLC7A11、DUSP2和DUSP5基因參與調節腫瘤的免疫微環境。有研究發現，lncRNA OIP5-AS1通過miR-128-3p/SLC7A11信號傳導抑制長期接觸鎘的前列腺癌鐵凋亡[22]。本研究發現，SLC7A11在髓系細胞群表達較高，提示lncRNA-miRNA-SLC7A11信號軸可通過調控PCa腫瘤微環境中的髓系細胞參與調控PCa細胞凋亡。

上皮細胞的lncRNA-miRNA-mRNA調控網絡主要通過調控EGLN3基因參與腫瘤相關微環境的調節。miR-1205可以通過靶向EGLN3的3′-UTR位點，下調其轉錄活性參與調控PCa細胞增殖[23]。本研究發現，EGLN3主要在上皮細胞中表達，提示lncRNA-miRNA-EGLN3可能參與PCA腫瘤微環境中的上皮細胞群的相關增殖。

間質細胞的lncRNA-miRNA-mRNA調控網絡主要通過調控GJA1、SNCG、TGFBR3、DUSP5、INMT、MET、KIT、PTGS2、PDPN、HAS2、MAP1B、INMT、RBPMS2、SNCG、MAP1B和ARC基因參與腫瘤相關微環境的調節。研究發現，上調的lncRNA DLX6-AS1通過miR-497-5p/SNCG軸在PCa中起到促進細胞增殖、抑制凋亡從而調控PCa進展的作用[24]。相關研究表明，在細胞實驗中MiR-124-3p通過直接吸附PTGS2抑制AKT/NF-κB通路，從而抑制PCa細胞的細胞活力、增殖、遷移、侵襲并促進凋亡[25]。這些結果提示本研究構建的lncRNA-miRNA-mRNA調控網絡，主要通過調控PCa的間質細胞群從而間接促進腫瘤細胞的遷移和侵襲等過程。

腫瘤細胞的lncRNA-miRNA-mRNA調控網絡主要參與TRIM36、TRIM29、TP53INP1、ITGA2、PDLIM5、AMOT、PRRG4和MET基因的調控。研究發現，lncRNA625可通過靶向miR-432調控TRIM29，從而調節Wnt/β-連環蛋白通路，抑制PCa細胞增殖，促進凋亡[26]。另外，lncRNA AGAP2-AS1/miR-195-5p/PDLIM5在PCa細胞增殖、遷移和侵襲以及腫瘤生長過程發揮調控作用[27]。通過預后分析得出腫瘤細胞lncRNA-miRNA-mRNA調控網絡中的靶基因ITGA2、PDLIM5H和PRRG4可以預測PCa的DFS，這些基因在蛋白水平也得到表達驗證，提示其可作為預后的觀測指標之一。

4 結論

本研究成功構建了在PCa單細胞水平下包括腫瘤細胞在內的9個細胞亞群的數百條精確的lncRNA-miRNA-mRNA調控網絡，并發現該調控網絡不僅參與PCa的發生發展調控，同時在臨床預后預測方面也有重要的意義。本研究構建的多條lncRNA-miRNA-mRNA調控軸雖然為PCa的基礎研究提供了新方向，但是目前在各個細胞亞群的具體生物學功能尚未得知，還需通過體外實驗進一步證實。