生物信息學(xué)預(yù)測(cè)橘皮素治療結(jié)直腸癌核心基因

宋子健，李建偉，胡和智

作為第三大惡性腫瘤的結(jié)直腸癌(colorectal cancer, CRC)具有惡性程度高、病程進(jìn)展迅速、易復(fù)發(fā)和轉(zhuǎn)移等特點(diǎn)，對(duì)人類健康和生命安全構(gòu)成重大威脅[1]。目前對(duì)于CRC的分子機(jī)制和核心基因的不完全了解阻礙了對(duì)CRC的各項(xiàng)研究。橘皮素是一類天然黃酮類物質(zhì)，屬于黃酮類化合物，廣泛存在于蕓香科植物川橘果皮、酸橙果皮和柑橘莖葉中，目前已被證實(shí)具有抑制細(xì)菌和抗腫瘤等藥理作用[2]。因此，該研究利用橘皮素治療結(jié)直腸癌的RNA-seq數(shù)據(jù)預(yù)測(cè)核心基因，并利用生存分析對(duì)其進(jìn)行預(yù)后分析，為CRC的診斷及治療藥物的研制提供新的作用靶點(diǎn)。

1 材料與方法

1.1 數(shù)據(jù)資料收集使用DMSO溶劑和橘皮素藥物處理CRC的HCT116細(xì)胞48 h，使用TRIzol提取實(shí)驗(yàn)組和溶劑對(duì)照組細(xì)胞的總RNA，使用逆轉(zhuǎn)錄試劑盒將其逆轉(zhuǎn)為cDNA，進(jìn)行RNA測(cè)序。最后，將細(xì)胞樣本分為3個(gè)橘皮素實(shí)驗(yàn)組和3個(gè)無橘皮素對(duì)照組。

1.2 差異基因篩選方法本研究篩選差異表達(dá)使用R語言中的程序包對(duì)橘皮素實(shí)驗(yàn)組和無橘皮素對(duì)照組的RNA-seq數(shù)據(jù)進(jìn)行基因差異表達(dá)分析，差異基因的篩選標(biāo)準(zhǔn)為|log2FC|>1和P<0.05。

1.3 lncRNA關(guān)鍵基因篩選根據(jù)基因種類，從篩選出的差異基因集中分別提取lncRNA差異表達(dá)基因集、miRNA差異表達(dá)基因集和mRNA差異表達(dá)基因集。對(duì)lncRNA差異表達(dá)基因集采用如下方法篩選得到lncRNA關(guān)鍵基因：①提高篩選標(biāo)準(zhǔn)，以|log2FC|>2和P<0.05為新閾值，分別篩選出差異表達(dá)更為顯著的lncRNA差異表達(dá)基因集和miRNA差異表達(dá)基因集，分別記為集合A和集合B，通過歸并排序算法對(duì)兩個(gè)集合均按照P值由小到大進(jìn)行排序；②從StarBase數(shù)據(jù)庫[3]中收集與miRNA差異表達(dá)基因集B中miRNA存在調(diào)控關(guān)系的lncRNAs，記為lncRNA基因集C；③將集合A與集合C取交集，得到的lncRNA基因集記為集合D；④利用GEPIA數(shù)據(jù)庫[4]中的臨床數(shù)據(jù)對(duì)集合D中的lncRNAs進(jìn)行生存分析，得到有顯著預(yù)后價(jià)值的lncRNA關(guān)鍵基因集。顯著lncRNA能夠通過大數(shù)據(jù)網(wǎng)站(如GEPIA等)預(yù)測(cè)它們和臨床病理參數(shù)的關(guān)系，以供后續(xù)研究使用。

1.4 lncRNA-miRNA-mRNA調(diào)控網(wǎng)絡(luò)構(gòu)建首先，利用DIANA網(wǎng)站[5]獲取lncRNA-miRNA調(diào)控關(guān)系數(shù)據(jù)，然后通過miRDB網(wǎng)站[6]獲取miRNA-mRNA調(diào)控關(guān)系數(shù)據(jù)。利用Cytoscape3.7.2軟件[7]構(gòu)建lncRNA-miRNA-mRNA調(diào)控網(wǎng)絡(luò)，將其中的mRNAs記為mRNA關(guān)鍵基因，以便后續(xù)研究。

1.5 基因功能注釋和通路富集分析為了更深層次了解橘皮素在治療CRC中的調(diào)控功能，對(duì)mRNA關(guān)鍵基因集進(jìn)行基因本體論(gene ontology,GO)分析和京都基因組百科全書(Kyoto encyclopedia of genes and genomes,KEGG)通路分析，從而找到差異基因的分子功能、參與的主要生物過程以及它們所屬的代謝通路等。DAVID網(wǎng)站[8]是一個(gè)生物信息數(shù)據(jù)庫，它整合了生物學(xué)數(shù)據(jù)和分析工具，主要用于GO富集分析和KEGG通路分析。將mRNA關(guān)鍵基因集導(dǎo)入到DAVID數(shù)據(jù)庫中進(jìn)行分析，相關(guān)閾值設(shè)置為P<0.05、Kappa Score=0.5、Min Level=5和Max Level=8。

1.6 mRNA核心基因集篩選通過STRING數(shù)據(jù)庫，收集mRNA關(guān)鍵基因集的蛋白質(zhì)互作網(wǎng)絡(luò)(PPI)數(shù)據(jù)，PPI分?jǐn)?shù)設(shè)置為0.4，并使用Cytoscape3.7.2軟件構(gòu)建相應(yīng)的PPI網(wǎng)絡(luò)[9]?；贛CODE算法和cytoHubba插件中拓?fù)浞治龇椒―egree、MNC和EPC分別對(duì)PPI網(wǎng)絡(luò)進(jìn)行分析。MCODE算法中的閾值設(shè)為：Node Density Cutoff=0.1、Node Score Cutoff=0.2、K-Core=2、Max.Depth=100;3個(gè)拓?fù)浞治龇椒ǖ拈撝翟O(shè)置為：Degree Score≥3、MNC≥2、EPC≥4.8，最終得到mRNA核心基因集。

1.7 統(tǒng)計(jì)學(xué)處理通過GEPIA數(shù)據(jù)庫對(duì)mRNA核心基因集進(jìn)行在線生存分析，驗(yàn)證其有效性。數(shù)據(jù)集限定為COAD和READ數(shù)據(jù)集，時(shí)間軸單位設(shè)置為月；基因表達(dá)差異采用t檢驗(yàn)，在CRC中表達(dá)量與預(yù)后的關(guān)系采用Log-rank檢驗(yàn)，以P<0.05表示差異有統(tǒng)計(jì)學(xué)意義。核心mRNA能夠通過大數(shù)據(jù)網(wǎng)站(如GEPIA等)預(yù)測(cè)它們和臨床病理參數(shù)的關(guān)系，以供后續(xù)研究使用。

2 結(jié)果

2.1 差異基因篩選RNA測(cè)序后共有21 460條基因數(shù)據(jù)，根據(jù)閾值P<0.05，|log2FC|>1對(duì)基因數(shù)據(jù)進(jìn)行差異表達(dá)分析，共得到2 614個(gè)差異表達(dá)基因，上調(diào)差異表達(dá)基因1 711個(gè)，下調(diào)差異表達(dá)基因903個(gè)。其中，含有197個(gè)lncRNA差異表達(dá)基因，128個(gè)miRNA差異表達(dá)基因，1 938個(gè)mRNA差異表達(dá)基因。差異表達(dá)基因的火山圖如圖1所示，圖中紅色為上調(diào)基因，綠色為下調(diào)基因。

圖1 RNA-seq數(shù)據(jù)的差異表達(dá)分析火山圖

2.2 lncRNA關(guān)鍵基因集以閾值|log2FC|>2且P<0.05為標(biāo)準(zhǔn)對(duì)197個(gè)lncRNA差異表達(dá)基因再次進(jìn)行篩選，得到116個(gè)lncRNA基因，命名為集合A；以閾值|log2FC|>2且P<0.05為標(biāo)準(zhǔn)對(duì)miRNA差異表達(dá)基因集進(jìn)行篩選，得到92個(gè)miRNA基因，命名為集合B；從StarBase數(shù)據(jù)庫中共收集到65個(gè)與集合B存在調(diào)控關(guān)系lncRNA基因，命名為集合C；對(duì)集合A和集合C取交集，得到32個(gè)lncRNA基因，命名為集合D；利用GEPIA數(shù)據(jù)庫進(jìn)行生存分析，共得到5個(gè)具有顯著預(yù)后價(jià)值的lncRNA基因(MALAT1、NEAT1、LINC00342、LINC01133、LINC00662)，記為lncRNA關(guān)鍵基因集，生存曲線如圖2所示。由圖2可知，5個(gè)lncRNA關(guān)鍵基因的LogrankP均<0.05，MALAT1(LogrankP=0.022)、NEAT1(LogrankP=0.013)、LINC00342(LogrankP=0.035)、LINC01133(LogrankP=0.036)和LINC00662(LogrankP=0.045),這證明5個(gè)lncRNA關(guān)鍵基因均對(duì)患者的生存周期有顯著影響。

2.3 lncRNA-miRNA-mRNA調(diào)控網(wǎng)絡(luò)針對(duì)5個(gè)lncRNA關(guān)鍵基因，從DIANA網(wǎng)站獲取651條與它們相關(guān)的lncRNA-miRNA調(diào)控關(guān)系數(shù)據(jù)，從miRDB數(shù)據(jù)庫獲取419條與它們相關(guān)的miRNA-mRNA調(diào)控關(guān)系數(shù)據(jù)?；谝陨蟽山M調(diào)控關(guān)系數(shù)據(jù)，利用Cytoscape3.7.2軟件的Merge功能構(gòu)建相應(yīng)的lncRNA-miRNA-mRNA調(diào)控網(wǎng)絡(luò)。lncRNA-miRNA-mRNA網(wǎng)絡(luò)由1 057個(gè)節(jié)點(diǎn)和909條邊組成，共有117個(gè)mRNA關(guān)鍵基因。調(diào)控網(wǎng)絡(luò)如圖3所示。

2.4 mRNA關(guān)鍵基因的GO和KEGG分析結(jié)果通過DAVID數(shù)據(jù)庫對(duì)117個(gè)mRNA關(guān)鍵基因進(jìn)行GO富集分析，在生物過程方面，mRNA關(guān)鍵基因功能主要集中于RNA聚合酶II啟動(dòng)子轉(zhuǎn)錄的調(diào)控和TOR信號(hào)通路等生物過程，表1給出了生物過程中差異性顯著的前5個(gè)條目。KEGG通路分析結(jié)果表明CRC的發(fā)生發(fā)展與PI3K-Akt信號(hào)通路密切相關(guān)。表2給出差異性顯著的前5條KEGG通路分析通路。通過GO富集分析和KEGG通路分析證實(shí)mRNA關(guān)鍵基因集與CRC密切相關(guān)。

圖2 5個(gè)lncRNA關(guān)鍵基因的生存曲線圖

表1 差異性顯著的前5個(gè)GO條目

表2 差異顯著的前5個(gè)KEGG通路分析條目

2.5 mRNA核心基因集將117個(gè)mRNA關(guān)鍵基因上傳至STRING數(shù)據(jù)庫，構(gòu)建對(duì)應(yīng)的PPI網(wǎng)絡(luò)，將結(jié)果導(dǎo)出并保存為tsv格式。利用Cytoscape3.7.2軟件將該P(yáng)PI網(wǎng)絡(luò)可視化，由48個(gè)節(jié)點(diǎn)和45條邊組成，如圖4所示。利用MCODE算法和cytoHubba插件中Degree、MNC及EPC拓?fù)浞治龇椒▽?duì)PPI網(wǎng)絡(luò)進(jìn)行分析。MCODE算法分析結(jié)果如圖5所示，拓?fù)浞治龇椒ㄇ?0名結(jié)果如表3所示。利用Python語言的numpy程序包對(duì)MCODE算法的結(jié)果基因集和拓?fù)浞治龇椒ǖ慕Y(jié)果基因集取交集，得到6個(gè)mRNA核心基因(FOS、GADD45A、CCND2、MYCN、BACH1和MXD1)。這些核心基因在橘皮素治療CRC中起到重要的調(diào)控作用，有成為生物標(biāo)志物和藥物靶點(diǎn)的潛力。

表3 拓?fù)浞治龇椒ㄇ?0名結(jié)果

圖3 lncRNA關(guān)鍵基因?qū)?yīng)的lncRNA-miRNA-mRNA調(diào)控網(wǎng)絡(luò)

圖4 mRNA關(guān)鍵基因?qū)?yīng)的PPI網(wǎng)絡(luò)

2.6 mRNA核心基因與患者的預(yù)后關(guān)系利用GEPIA數(shù)據(jù)庫對(duì)6個(gè)mRNA核心基因進(jìn)行生存分析，各自的生存曲線圖如圖6所示。其中，F(xiàn)OS、CCND2和MXD1表達(dá)水平對(duì)患者的總生存時(shí)間有著顯著影響(P<0.05)。而GADD45A、MYCN和BACH1對(duì)患者的生存率影響差異無統(tǒng)計(jì)學(xué)意義。

圖5 MCODE分析結(jié)果圖

3 討論

CRC的發(fā)生與外界環(huán)境、行為方式、遺傳等多種因素密切相關(guān)，盡管目前對(duì)CRC的研究已取得了較大進(jìn)步，但是CRC的預(yù)后仍然效果不佳。隨著測(cè)序技術(shù)的飛速發(fā)展和生物信息技術(shù)的不斷突破，CRC的分子機(jī)制研究成為了當(dāng)前的一個(gè)熱點(diǎn)，尋找CRC診斷及預(yù)后的生物標(biāo)志物和藥物靶點(diǎn)為CRC的診療提供了新的思路。橘皮素已被證實(shí)具有抑制細(xì)菌和抗腫瘤等藥理作用，但對(duì)橘皮素治療CRC的分子機(jī)制卻不甚了解。因此，本文以橘皮素治療CRC的RNA-seq數(shù)據(jù)為研究對(duì)象，通過生物信息學(xué)分析方法，篩選出117個(gè)mRNA關(guān)鍵基因。GO富集分析和KEGG通路分析結(jié)果表明，關(guān)鍵基因均富集在與CRC相關(guān)的功能和通路上。其中，KEGG分析結(jié)果表明PI3K-Akt信號(hào)通路與CRC密切相關(guān)。研究[10]表明，SPOCK1在CRC細(xì)胞系中過表達(dá)，沉默SPOCK1可逆轉(zhuǎn)CRC細(xì)胞中的EMT過程，顯著減弱了遷移/侵襲，抑制體外增殖和體內(nèi)腫瘤的生長。敲除SPOCK1明顯降低了HCT116細(xì)胞中p-PI3K和p-Akt的蛋白表達(dá)水平。此外，mRNA關(guān)鍵基因還顯著富集在乙型肝炎、肺結(jié)核、胰腺癌、甲型流感、前列腺癌等多種疾病的相關(guān)信號(hào)通路，提示橘皮素對(duì)多種疾病具有治療作用，為今后的相關(guān)研究提供了新的思路。

圖6 6個(gè)mRNA核心基因表達(dá)與患者預(yù)后的生存曲線

通過生存分析，從mRNA關(guān)鍵基因中篩選出3個(gè)與CRC預(yù)后密切相關(guān)的核心基因(FOS、CCND2和MXD1)。其中，F(xiàn)OS和CCND2已被文獻(xiàn)[11-12]證實(shí)與CRC有密切關(guān)系。有研究[13]表明，MXD1參與了乳腺癌癌細(xì)胞的增殖和轉(zhuǎn)移過程，在乳腺癌組織中MXD1表達(dá)顯著下調(diào)，并影響了乳腺癌患者的預(yù)后。因此，課題組推斷MXD1也極有可能與CRC的發(fā)生和發(fā)展相關(guān)。綜上所述，3個(gè)核心基因有成為生物標(biāo)志物和藥物靶點(diǎn)的可能，對(duì)CRC的發(fā)病機(jī)制及治療提供了新的思路，也為CRC藥物靶點(diǎn)研究提供重要參考。