胰腺癌中lncRNA和mRNA表達譜分析及調控網絡研究

張美晨，魏文竹，盧星宇，湯文瀛，王迪，劉躍娟

（哈爾濱醫科大學大慶校區醫學信息學院，黑龍江 大慶 163319）

0 引言

胰腺癌（pancreaticcancer，PC）因其起病隱匿、早期診斷難、手術切除困難、放化療不敏感以及預后差、死亡率高等原因，成為危害最大的惡性腫瘤之一。我國胰腺癌發病率近20年來增長了將近6倍，居惡性腫瘤死亡率的第5位，因此積極尋找對胰腺癌有更好診斷意義的臨床指標尤為重要。目前，越來越多研究表明，癌癥的發生與基因突變、缺失和表達障礙有關。本研究利用胰腺癌芯片數據，計算得到差異表達的mRNA和lncRNA，結合蛋白質互作數據，構建失調的lncRNA-mRNA網絡，通過挖掘網絡功能模塊，初步探索其作用機制，為揭示胰腺癌發生發展過程中lncRNA與其靶基因的關系和作用機制提供新線索。

1 材料與方法

1.1 胰腺癌的lncRNA和mRNA表達譜

我們從基因表達數據庫(Gene Expression Omnibus，GEO)中下載了編號為GSE61166的基因芯片數據。本研究從中選取4個胰腺癌患者和4個胰腺炎患者胰腺組織的lncRNA和mRNA表達譜數據作為研究對象。

1.2 BLAST序列比對

利用BLAST進行序列比對，在基因芯片GPL61166平臺上下載序列信息，同時在GENECODE數據庫下載人類基因序列，利用計算機軟件比對，通過BioMart數據庫，將Ensembl Gene ID轉換為Gene Symbol。然后我們做了標準化和以2為底的對數轉換，若多個探針集對應一個基因，則取其平均值作為該基因的表達量，共獲得6530個mRNA和31個lncRNA，全部數據處理流程如圖1所示。

1.3 差異基因篩選

差異基因通過基因芯片顯著性分析方法（SAM）篩選，SAM方法根據芯片數據噪音大小與表達峰度進行修改，其優點是在得到較多特征基因的同時，錯誤發現率（FDR）還保持在較低水平。利用SAM方法，閥值設定為P＜0.05，得到差異表達的820個mRNA和16個lncRNA。

1.4 相關分析

圖1 流程圖

相關分析法是對具體有依存關系的現象探討其相關方向以及程度的一種統計方法。本研究對差異表達的820個mRNA和16個lncRNA計算皮爾森相關系數，取系數大于0.8且P＜0.05的強相關對構建lncRNA-mRNA共表達網絡。

1.5 構建lncRNA-mRNA互作網絡

根據上一步得到的共表達網絡，從HPRD（Human Protein Reference Database）數據庫下載人類蛋白質和蛋白質相互作用（PPI）數據，得到人的PPI網絡，并將lncRNA-mRNA共表達網絡和PPI網絡取交集，然后對交集的mRNAs在PPI中選取直接鄰居結點，和原共表達網絡一起形成lncRNA-mRNA互作網絡。

1.6 網絡模塊挖掘

利用Cytoscape軟件中的JActiveModules插件進行模塊挖掘，獲得5個網絡功能模塊。用ClusterOne和Netmine對模塊進行聚類和進一步挖掘，尋找p-value＜0.05的顯著性lncRNA-mRNA互作對。

1.7 Gene ontology(GO)分析

通過來自R/bioConductor的clusterProfiler包完成了對模塊中編碼基因的GO terms富集度統計學分析，計算出這些基因的GO term的p-value和p-value的FDR值，定位這些基因最可能相關的GO term。

2 結果

2.1 差異表達 lncRNAs、mRNAs

對胰腺癌表達譜數據進行處理，得到差異表達的mRNAs、lncRNAs數量分別為820、16個。通過比較差異mRNAs的前5%和全部的lncRNAs在不同樣本中的表達量（見圖2、圖3），我們發現這些mRNAs與lncRNAs在胰腺癌患者中低表達。

圖2 lncRNA和mRNA表達

2.2 網絡的全局特征

構建的lncRNA—mRNA共表達網絡共包含了826個結點和7176條邊，涉及到16個lncRNA和810個mRNA，利用Cytoscape可視化網絡（如圖4）對共表達網絡和PPI網絡整合并得到一個包含2373個結點和8626條邊，涉及到16個lncRNA和2357個mRNA的lncRNA—mRNA互作網絡。通過拓撲屬性分析(如圖5)，我們發現MIR7-3HG具有較大的度585，平均鄰居結點12，較小的聚類系數0.447，大多數的路徑長度在2-5之間。

圖3 lncRNA和mRNA直方圖

圖4 lncRNA-mRNA共表達網絡

圖5 網絡度分布

2.3 網絡功能模塊

對lncRNA—mRNA互作網絡進行模塊挖掘，得到5個模塊。我們發現LINC00670與CCL2關系對在四個模塊中都出現，LINC00670是一個存在轉錄因子CTCF綁定位點的長鏈非編碼RNA，CTCF蛋白在印記調控區域和分化甲基化區域1和MAR3結合可抑制胰島素樣生長因子2基因[14]。另外，我們對與LINC00670有互作關系的靶基因做GO功能富集分析發現，其靶基因以最顯著的形式富集在了胰液分泌通路上。

在模塊一中，LINC00670的介數達到0.021,遠高于平均介數值0.008。CCL2趨化因子及其受體可加強腫瘤細胞的增殖和侵襲能力，并聚集免疫細胞等機制以促進腫瘤的發生發展，CCL2-CCR2信號通道在腫瘤細胞的轉移及增殖過程中發揮著重要的作用，尤其是在腫瘤細胞轉移的早晚期[15]。LINC00670同時還與IL6和NFASC相連。其中，IL6（白細胞介素-6），是炎癥介質網絡的關鍵成分。有研究表明IL6與癌細胞增殖以及許多胰腺癌病人的惡病質過程有關。在模塊二中受體活性修飾蛋白1（RAMP1）可以與降鈣素受體（CTR）直接作用產生胰淀粉樣酶(AMY)受體，在胰腺腫瘤引起的胰腺導管阻塞、胰腺膿腫、胰腺損傷中AMY高表達。有研究表明ADM能夠增強胰腺癌細胞的侵襲性。與LINC00671相聯的IAPP在II型糖尿病中會引起β細胞凋亡。用一種具有中和抗體的細胞特異性促進劑明顯減少了病變的活性巨噬細胞積累，從而減少纖維化和病灶的生長。在模塊四中CKM,DCX,IL6,IAPP都在模塊一、二、三中出現，CKM,DCX,IL6同時連接LINC00670。有研究表明在胰腺癌中從UCM中更快地釋放出DOX，在凋亡的外圍細胞中可能會停止內吞作用和胞吐作用，有研究表明組織和血液的MMP9增加反映了胰腺癌與糖尿病相關。

2.4 模塊的功能分析

對模塊中編碼基因進行GO功能富集分析結果見圖6。通過來自R/bioConductor的clusterProfiler數據包，我們完成了模塊中編碼基因的GO term富集度統計學分析，閥值選取為Bonferroni校正P值 ＜0.01。

圖6第一排從左至右分別是模塊一、三、五，第二排從左至右分別是模塊二和四。得到的五個模塊中都有與炎癥顯著性功能相關的注釋；模塊1中基因功能注釋到炎癥、增殖、凋亡；幾個模塊里還有關于侵襲轉移的功能注釋。其中多個功能有重復出現。如:taxis、regulation of developmental process、multi-organism cellular process。 通 過GO注釋的模塊功能與免疫反應有關。

3 結論

通過本項目的初步研究，找到與胰腺癌發生密切相關的的五個模塊，而這些模塊又分別與炎癥、免疫、細胞增殖和凋亡等有關。表明mRNA和lncRNA通過多種途徑參與胰腺癌的發生及發展過程，在多個模塊中發現lncRNA00670與胰腺癌密切相關。隨著mRNA和lncRNA研究技術的不斷發展和成熟，利用生物信息學分析將會發掘更多的功能性mRNA和lncRNA，為診斷和治療胰腺癌等惡性病疾病提供新方向。

圖6 GO功能富集注釋

[1] Liu R C, Traverso L W. Diagnostic laparoscopy improves staging of pancreatic cancer deemed locally unresectable by computed tomography[J]. Surgical Endoscopy And Other Interventional Techniques, 2005, 19(5): 638-642.

[2] Lin Q J, Yang F, Jin C, et al. Current status and progress of pancreatic cancer in China[J]. World Journal of Gastroenterology:WJG, 2015, 21(26): 7988.

[3] Li N, Yang M, Shi K, et al. Long non-coding RNA HOXA11-AS in human cancer: A meta -analysis[J]. Clinica Chimica Acta,2017,474:165-170.

[4] Bister V, Skoog T, Virolainen S, et al. Increased expression of matrix metalloproteinases-21 and-26 and TIMP-4 in pancreatic adenocarcinoma[J]. Modern Pathology, 2007, 20(11): 1128 -1140.

[5] Kayed H, Jiang X, Keleg S, et al. Regulation and functional role of the Runt-related transcription factor-2 in pancreatic cancer[J].British journal of cancer, 2007, 97(8): 1106-1115.

[6] Namkoong H, Shin S M, Kim H K, et al. The bone morphogenetic protein antagonist gremlin 1 is overexpressed in human cancers and interacts with YWHAH protein[J]. BMC cancer, 2006, 6(1): 74.

[7] Zhu Y, Xu H, Chen H, et al. Proteomic analysis of solid pseudopapillary tumor of the pancreas reveals dysfunction of the endoplasmic reticulum protein processing pathway[J]. Molecular &Cellular Proteomics, 2014, 13(10): 2593-2603.

[8] Wang Y, Yu S, Huang D, et al. Cellular Prion Protein Mediates Pancreatic Cancer Cell Survival and Invasion through Association with and Enhanced Signaling of Notch1[J]. The American journal of pathology, 2016, 186(11): 2945-2956.

[9] Qin C F, Zhao F L. Long non-coding RNA TUG1 can promote proliferation and migration of pancreatic cancer via EMT pathway[J].European Review for Medical and Pharmacological Sciences, 2017,21: 2377-2384.

[10] Tahira A C, Dazzani B, Reis E M, et al. Long noncoding intronic RNAs are differentially expressed in primary and metastatic pancreatic cancer[J]. Molecular cancer, 2011, 10(1): 141.

[11] Ma C, Nong K, Zhu H, et al. H19 promotes pancreatic cancer metastasis by derepressing let-7’s suppression on its target HMGA2-mediated EMT[J]. Tumor Biology, 2014, 35(9): 9163-9169.

[12] Saung M T, Zheng L. Current standards of chemotherapy for pancreatic cancer[J]. Clinical therapeutics, 2017.

[13] Yi X, Zai H, Long X, et al. Krüppel-like factor 8 induces epithelialto-mesenchymal transition and promotes invasion of pancreatic cancer cells through transcriptional activation of four and a half LIM-only protein 2[J]. Oncology letters, 2017, 14(4): 4883-4889.

[14] Huang X, Ta N, Zhang Y, et al. Microarray Analysis of the Expression Profile of Long Non-Coding RNAs Indicates lncRNA RP11-263F15. 1 as a Biomarker for Diagnosis and Prognostic Prediction of Pancreatic Ductal Adenocarcinoma[J]. Journal of Cancer, 2017, 8(14): 2740.

[15] Li H, Wang X, Fang Y, et al. Integrated expression profiles analysis reveals novel predictive biomarker in pancreatic ductal adenocarcinoma[J]. Oncotarget, 2017, 8(32): 52571.