生物信息學預測橘皮素治療結直腸癌核心基因

宋子健，李建偉，胡和智

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

1 材料與方法

1.1 數據資料收集使用DMSO溶劑和橘皮素藥物處理CRC的HCT116細胞48 h，使用TRIzol提取實驗組和溶劑對照組細胞的總RNA，使用逆轉錄試劑盒將其逆轉為cDNA，進行RNA測序。最后，將細胞樣本分為3個橘皮素實驗組和3個無橘皮素對照組。

1.2 差異基因篩選方法本研究篩選差異表達使用R語言中的程序包對橘皮素實驗組和無橘皮素對照組的RNA-seq數據進行基因差異表達分析，差異基因的篩選標準為|log2FC|>1和P<0.05。

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

1.4 lncRNA-miRNA-mRNA調控網絡構建首先，利用DIANA網站[5]獲取lncRNA-miRNA調控關系數據，然后通過miRDB網站[6]獲取miRNA-mRNA調控關系數據。利用Cytoscape3.7.2軟件[7]構建lncRNA-miRNA-mRNA調控網絡，將其中的mRNAs記為mRNA關鍵基因，以便后續研究。

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

1.6 mRNA核心基因集篩選通過STRING數據庫，收集mRNA關鍵基因集的蛋白質互作網絡(PPI)數據，PPI分數設置為0.4，并使用Cytoscape3.7.2軟件構建相應的PPI網絡[9]?；贛CODE算法和cytoHubba插件中拓撲分析方法Degree、MNC和EPC分別對PPI網絡進行分析。MCODE算法中的閾值設為：Node Density Cutoff=0.1、Node Score Cutoff=0.2、K-Core=2、Max.Depth=100;3個拓撲分析方法的閾值設置為：Degree Score≥3、MNC≥2、EPC≥4.8，最終得到mRNA核心基因集。

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

2 結果

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

圖1 RNA-seq數據的差異表達分析火山圖

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

2.3 lncRNA-miRNA-mRNA調控網絡針對5個lncRNA關鍵基因，從DIANA網站獲取651條與它們相關的lncRNA-miRNA調控關系數據，從miRDB數據庫獲取419條與它們相關的miRNA-mRNA調控關系數據?；谝陨蟽山M調控關系數據，利用Cytoscape3.7.2軟件的Merge功能構建相應的lncRNA-miRNA-mRNA調控網絡。lncRNA-miRNA-mRNA網絡由1 057個節點和909條邊組成，共有117個mRNA關鍵基因。調控網絡如圖3所示。

2.4 mRNA關鍵基因的GO和KEGG分析結果通過DAVID數據庫對117個mRNA關鍵基因進行GO富集分析，在生物過程方面，mRNA關鍵基因功能主要集中于RNA聚合酶II啟動子轉錄的調控和TOR信號通路等生物過程，表1給出了生物過程中差異性顯著的前5個條目。KEGG通路分析結果表明CRC的發生發展與PI3K-Akt信號通路密切相關。表2給出差異性顯著的前5條KEGG通路分析通路。通過GO富集分析和KEGG通路分析證實mRNA關鍵基因集與CRC密切相關。

圖2 5個lncRNA關鍵基因的生存曲線圖

表1 差異性顯著的前5個GO條目

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

2.5 mRNA核心基因集將117個mRNA關鍵基因上傳至STRING數據庫，構建對應的PPI網絡，將結果導出并保存為tsv格式。利用Cytoscape3.7.2軟件將該PPI網絡可視化，由48個節點和45條邊組成，如圖4所示。利用MCODE算法和cytoHubba插件中Degree、MNC及EPC拓撲分析方法對PPI網絡進行分析。MCODE算法分析結果如圖5所示，拓撲分析方法前10名結果如表3所示。利用Python語言的numpy程序包對MCODE算法的結果基因集和拓撲分析方法的結果基因集取交集，得到6個mRNA核心基因(FOS、GADD45A、CCND2、MYCN、BACH1和MXD1)。這些核心基因在橘皮素治療CRC中起到重要的調控作用，有成為生物標志物和藥物靶點的潛力。

表3 拓撲分析方法前10名結果

圖3 lncRNA關鍵基因對應的lncRNA-miRNA-mRNA調控網絡

圖4 mRNA關鍵基因對應的PPI網絡

2.6 mRNA核心基因與患者的預后關系利用GEPIA數據庫對6個mRNA核心基因進行生存分析，各自的生存曲線圖如圖6所示。其中，FOS、CCND2和MXD1表達水平對患者的總生存時間有著顯著影響(P<0.05)。而GADD45A、MYCN和BACH1對患者的生存率影響差異無統計學意義。

圖5 MCODE分析結果圖

3 討論

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

圖6 6個mRNA核心基因表達與患者預后的生存曲線

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