用生物信息學方法分析鑒定肺腺癌中的hub 基因

呂金海

（懷化學院，湖南 懷化 418008）

肺癌是常見的癌癥類型之一[1]，其分為小細胞肺癌（Small Cell LungCancer，SCLC）和非小細胞肺癌（Nonsmall Cell LungCancer，NSCLC）2 種類型。其中，NSCLC約占肺癌的85%[2]。肺腺癌（LungAdenocarcinoma，LUAD）是非小細胞肺癌的主要亞型之一。肺癌患者早期無明顯癥狀，大部分患者被發現時已經處于中晚期[3]。肺腺癌與多種風險因素有關，其成因與發病機制尚不完全清楚。

在癌癥基因組圖譜（The Cancer Genome Atlas，TCGA）、國際癌癥基因組聯盟（International Cancer Genome Consortium，ICGC）等大樣本腫瘤研究項目開展后，腫瘤樣本中生成了大量對腫瘤研究有幫助的數據[4-5]。此外，一些小規模癌癥項目的數據也被存入高通量基因表達數據庫（Gene Expression Omnibus，GEO）等數據庫[6]中，為生物信息學分析提供了一些幫助。

在各種實驗技術都存在局限性的情況下[7]，基因芯片技術和高通量測序技術的應用給予了癌癥研究一個新視角。公共數據庫中的大量數據為不同類型癌癥基因表達的綜合分析提供了便利。最近，一些研究人員將生物信息學得到的結果與臨床數據相結合，發現了用于不同類型癌癥的診斷、治療和預后的新的生物標志物[8-10]，這類新的生物標志物被稱為hub 基因。基于此，CHEN 等人[11]發現了4 個有望作為生物標記物的hub 基因。GUO 等人[12]鑒定出了9 個與腎上腺皮質癌相關的hub 基因并報道了相關通路。

本研究主要是利用綜合生物信息學探究癌癥的底層發病原理，并尋找其中可能的分子作用機制來改進診斷方法和預防措施。我們從GEO 中獲取了GSE116959[13]和GSE118370[14]2 張基因芯片，限定log2FC值大于1.00 且p值小于5.00×10-2，得到了差異表達基因（Differentially Expressed Genes，DEGs），并將其進行基因本體論（Gene Ontology，GO）和基因組的京都百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）富集分析。此后，我們組建了完整的差異表達基因的蛋白質-蛋白質互作（Protein-protein Interaction，PPI）網絡和加權基因共表達網絡（Weighted Gene Coexpression Network Analysis，WGCNA）。本文利用一系列生物信息學方法，確定了hub 基因，并使用富集分析來確定與肺腺癌相關的關鍵通路。我們進行了生存分析來探究hub 基因的表達與肺腺癌預后的關系。本研究使用了多個平臺的數據進行綜合分析，大大提高了hub 基因的可靠性。本研究結果可從分子水平進一步了解肺腺癌的發生發展，為肺腺癌的診斷、治療、監測和預后提供一些準確、有實用價值的標志物。

1 材料與方法

1.1 數據來源

GEO 數據庫是由美國國家生物技術信息中心（National Center for Biotechnology Information，NCBI）于2000 年創立的公共基因表達數據庫。我們通過GEO數據庫評估后選擇并下載了2 個基因表達譜數據集：GSE116959 和GSE118370，其中GSE116959 含有68個樣本（57 個肺腺癌樣本和11 個正常對照樣本），GSE118370 含有12 個樣本（6 個肺腺癌樣本和6 個正常對照樣本）。

1.2 方法

1.2.1 批次效應的消除和DEGs 篩選

在高通量實驗中，異構性和潛在變量被認為是偏差和變異性的主要來源。基因組實驗中最著名的潛在變異來源是批次效應。批次效應的影響可能很嚴重，甚至完全損害生物學結果。

我們首先將GEO 中得到的數據使用R 包“tidyverse”中的“left_join”命令進行合并，然后使用R包“sva”中的“ComBat”函數消除批次效應后，再選用R 包“limma”[15]比對肺腺癌組織與正常組織的數據，最后篩選出矯正后P＜5.00×10-2的數據，其中上調基因設置為“log2FC＞1.00”，下調基因設置為“log2FC＜-1.00”。

1.2.2 DEGs 的GO 分析和KEGG 分析

轉錄組學和蛋白質組學生成的數據越來越多，我們需要采用綜合策略進行分析。GO 數據庫將基因與該數據庫系統中的生物學術語緊密聯系，并把基因功能分為3 個部分：細胞組分（Cellular Component，CC）、分子功能（Molecular Function，MF）、生物過程（Biological Process，BP）[16]。KEGG 是一個包含基因表達途徑和功能的數據庫，旨在通過基因組信息和生物信息學來預測細胞中復雜的通路及生物學行為[17]。R 包“clusterProfiler”支持3 種物種，包括人類、小鼠和酵母。它提供了一種基因分類方法：groupGO，并可以自動進行基因簇富集分析[18]。我們對獲得的DEGs進行GO分析和KEGG 富集分析，并設置P=5.00×10-2為篩選閾值。

1.2.3 PPI 網絡構建與分析

我們首先將篩選出的DEGs 輸入在線數據庫檢索相互作用基因/蛋白質的搜索工具（Search Tool for the Retrieval of Interacting Genes/Proteins，STRING）數據庫[19]（https：//www.string-db.org），然后對Cytoscape[20]（https：//cytoscape.org）注釋和分析得到的PPI 網絡使用cytoHubba[21]插件進行hub 基因的篩選，最后通過MCODE[22]插件識別PPI 網絡中最重要的模塊。

1.2.4 加權基因共表達網絡的構建

WGCNA 是一種系統生物學方法，用于描述微陣列樣品中基因之間的相關模式。WGCNA 軟件包包含一組全面的功能，用于對大型高維數據集進行相關網絡分析，其使用無監督聚類來識別基因模塊，我們使用R 包“WGCNA”[23]構建加權基因共表達網絡。

軟閾值是基于近似無尺度網絡的一種準則，同時也是可以減少誤差的一種方式。它使得構建的網絡更符合冪律分布和無尺度網絡特征，并使得結果更具有生物數據特征。本文首先通過pick Soft Threshold 函數計算出軟閾值12，緊接著通過該軟閾值構建加權基因共表達網絡，然后對基因進行聚類后使用動態剪樹法識別基因模塊合并，其次將相似的基因聚類到相同顏色的模塊中，最后計算皮爾森相關系數以評估模塊與表型之間的潛在關聯性。

1.2.5 基于Oncomine 數據庫的Meta 分析

微陣列研究的一個問題是假陽性結果很普遍，而Meta 分析可以從各種各樣的微陣列數據集中識別和評估多個基因表達特征的交集，從而真實反映腫瘤的基本轉錄特征。

賽默飛公司開發了商業化的Oncomine 數據庫[24]，該數據庫是目前世界上最大的癌基因芯片數據庫和整合數據挖掘平臺。用戶可以使用Oncomine 選擇適當的研究進行比較，然后通過Meta 分析來確定多個獨立研究中顯著過表達或低表達的基因。本研究通過Oncomine 在線數據庫分析篩選出的hub 基因在肺腺癌組織與正常組織中的表達情況來驗證hub 基因的可靠性。

1.2.6 生存分析

2017 年上線的GEPIAWeb 是其基于TCGA 和基因型-組織表達（Genotype-Tissue expression，GTex）的服務器。GEPIA2（http：//GEPIA.Cancer-pku.cn）是GEPIAWeb經過更新和增強的版本，可提供更多功能[25]。GEPIA2 可根據基因表達水平進行生存分析，該功能允許用戶篩選基因在不同癌癥類型中的預后影響。為了進一步評估hub 基因在本研究中的預后重要性，我們在此平臺上進行了總體生存分析，用來檢驗每個基因的表達與總生存率之間的關系。

2 結果

2.1 DEGs 篩選

本文將GSE116959 與GSE118370 的數據進行集成并消除批次效應后，設置|log2FC|＞1（log2FC=log2病人平均表達量/正常人平均表達量），校正后p＜0.05（p大于0.05，被認為沒有統計學意義），篩選得到1433 個差異基因后對差異基因的表達模式進行了研究，圖1 顯示上調的523 個差異基因，下調的910 個差異基因。

圖1 差異基因火山圖

2.2 GO 分析與KEGG 分析

對篩選出來的差異基因做GO 分析和KEGG 分析。GO 分析的結果顯示，BP 主要圍繞細胞外基質組織、細胞外結構組織和有絲分裂的細胞核分裂；CC 主要表現在含膠原的細胞外基質、分泌顆粒內腔和膠原三聚體等結構；MF 主要體現在細胞外基質結構組成、糖胺聚糖結合和肽聚合。KEGG 分析顯示，差異基因主要富集在補體和凝血級聯、細胞外基質（Extracellular matrix，ECM）受體相互作用和黏著作用的通路中，具體見表1。

表1 富集分析結果表

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

使用STRING 分析DEGs 間的相互作用。結果顯示共有716 個DEGs 網絡節點和3759 條邊被用于構建PPI 網絡。圖2 為使用cytoHubba 中的Degree算法初步篩選出20 個hub 基因：CDK1、CDC20、CCNA2、CCNB1、KIF11、CCNB2、KIF20A、KIF2C、GNG11、BUB1、FPR2、TOP2A、CENPE、PIK3R1、CDCA8、BUB1B、UBE2C、MAD2L1、AGT、ASPM。MCODE 可檢測可能代表分子復合物的大型PPI 網絡中的密集連接區域。該方法通過局部鄰域密度和從局部密集的種子蛋白向外遍歷的頂點加權，根據給定的參數隔離密集的區域。圖3為使用MCODE 鑒定出了3 個最重要的基因簇。

圖2 hub 基因網絡圖

圖3 MCODE 模塊圖

2.4 加權基因共表達網絡的構建與分析

本文以12 為軟閾值構建網絡模塊，通過對模塊的進一步分析，計算模塊間的差異，并構建基因模塊樹狀圖（見圖4）。本文使用動態剪切樹法識別基因模塊，并將相似度較高的模塊合并，最終得到18 個模塊（見圖5）。18 個不同帶狀矩形框代表18 個基因集模塊。圖4 顯示樹狀圖上的每一個峰都對應1 個基因，相似的基因被聚類到相同帶狀矩形框的模塊中。表2顯示數據集分為肺腺癌和正常2 種表型。相關系數大于0.6 的模塊與性狀顯著相關。

表2 模塊-性狀相關性表

圖4 基因模塊樹狀圖

圖5 基因矩形模塊圖

綜合WGCNA 的基因模塊樹狀圖、基因矩形模塊圖、模塊-性狀相關性表和PPI 網絡分析的結果，我們初步篩選出17 個hub 基因：CDK1、CDC20、CCNA2、CCNB1、KIF11、CCNB2、KIF20A、KIF2C、BUB1、TOP2A、CENPE、PIK3R1、CDCA8、BUB1B、UBE2C、MAD2L1、ASPM，并進行了進一步的分析和驗證。

2.5 Oncomine 數據庫差異表達分析

使用Oncomine 數據庫，設置Analysis Type：Cancer vs.NormalAnalysis；Cancer Type：LungCancer；DataType：mRNA 選取肺腺癌患者和正常人的數據進行meta 分析，詳細數據見表3。除沒有研究發現PIK3R1 在肺腺癌患者與正常人中的表達量有差異外，所有hub 基因在肺腺癌與正常人中的表達量都具有統計學差異（p＜0.05）。

表3 Oncomine 數據庫meta 分析結果

2.6 癌癥與正常人hub 基因相關的生存分析

總生存時間是指從隨機化開始至因任何原因引起死亡的時間，總生存期的延長可以體現確切的臨床獲益，是抗腫瘤藥物最可靠的療效評價指標。表4 生存分析結果將病人依據表達水平的中位數分為高表達組和低表達組，對得到的hub 基因進行生存分析以探究基因表達量與總生存時間之間的相關性。所有hub 基因的高表達都與肺腺癌患者的總體生存時間縮短相關。這提示了這些hub 基因在不同程度上與肺腺癌的病理和生理有關，可以作為監測肺腺癌預后的潛在生物標志物以判斷肺腺癌的嚴重程度、預測患者的生存時間，或是作為治療靶點。

3 結論

1）共篩選出肺腺癌相關16 個hub 基因。從高通量GEO 數據庫中獲得基因表達譜序列號，并使用R中的“limma”包對DEGs 進行鑒定，篩選發現1433 個DEGs（上調523 個，下調910 個）。為確定差異基因所涉及的生物學過程及信號通路，我們進行了GO 富集分析和基因組的KEGG 富集分析，發現DEGs 主要與細胞外基質組織、細胞外結構組織和有絲分裂的細胞核分裂有關，主要富集在補體和凝血級聯、ECM受體相互作用和黏著作用的通路上。采用PPI 分析得到了由716 個DEGs 網絡節點和3759 條邊構成的PPI 網絡，并在其中初步篩選出20 個hub 基因。通過WGCNA分析從DEGs 中得到了18 個模塊，對hub 基因進行鑒定，篩選出17 個hub 基因。使用Oncomine 數據庫對hub 基因在肺腺癌患者與正常人的表達情況進行Meta 分析，共鑒定出16 個hub 基因與肺腺癌相關。16個hub 基因分別為BUB1B、CDK1、CDC20、CCNA2、CCNB1、KIF11、CCNB2、KIF20A、KIF2C、BUB1、TOP2A、CENPE、CDCA8、UBE2C、MAD2L1、ASPM。

2）16 個hub 基因通常在肺腺癌患者中過表達。分析表明肺腺癌表型16 個hub 基因表達高于正常表型，據KEGG 分析，這些hub 基因主要與細胞周期、細胞衰老和p53 信號通路有關。GO 分析也證明它們顯著參與了細胞核分裂，細胞器分裂和染色體分離過程。根據KEGG 分析和GO 分析的結果，這些基因與染色體的變化有關，可能在腫瘤的發生和發展中起到重要作用。

3）hub 基因高表達的組總體生存時間低于低表達的組。癌癥與正常人hub 基因相關的生存分析發現，在肺腺癌患者中，這些hub 基因高表達的組總體生存時間下降。因此，這16 個hub 基因可能與肺腺癌的病理和生理密切相關，也可能是潛在的生物標志物。

4 展望

2 份關于LUAD 的數據集在篩選后一共發現1433 個差異表達基因，通過GO 和KEGG 等注釋工作發現其DEGs 功能板塊聯系緊密，共涉及53 條生物學通路。這53 條通路和16 個hub 基因為LUAD 發生和發展提供了更詳細的分子機制，有望作為潛在的生物標志物和治療靶點。雖然肺腺癌的治療在近年有所改善，但它仍然是最具攻擊性和致命性的腫瘤類型之一[26]。肺腺癌分子機制的識別對肺腺癌的診斷、治療和預后至關重要。微陣列基因表達譜被廣泛用于篩選參與腫瘤發生的DEGs，為臨床提供有價值的信息[27]。

癌癥基因普查是一項正在進行的工作，目的是對那些包含與癌癥有因果關系的突變的基因進行分類，并解釋這些基因的功能異常如何導致癌癥[28]。在16個hub 基因中，BUB1B 編碼參與紡錘體檢查點功能的激酶，并在細胞分裂過程中控制染色體分離[29]；CDK1編碼的蛋白質對于真核細胞G1/S 和G2/M 期的過渡至關重要[30]；CDC20 表達的抑制可降低細胞增殖并誘導G2/M 細胞周期阻滯[31]；CCNA2 通過整合素αVβ3信號通路促進肺癌細胞的侵襲和遷移[32]；CCNB1 編碼一種有絲分裂調控蛋白[33]；KIF11 在紡錘體功能中起關鍵作用[34]；CCNB2 編碼細胞周期蛋白，表達下調能抑制肺腺癌的發生[35]；KIF20A 敲除的可改變肺癌細胞表型并調控JNK 通路[36]；KIF2C 編碼微管解聚和染色體分離的調節因子[37]；BUB1 編碼絲氨酸/蘇氨酸蛋白激酶與激活紡錘體啟動子中發揮作用[38]；TOP2A是公認的致癌基因；CENPE 可以促進肺腺癌細胞的增殖[39]；CDCA8 的表達能刺激的乳腺癌細胞生長和存活[40]；UBE2C 編碼E2 泛素結合蛋白參與了細胞癌變[41]；MAD2L1 是有絲分裂紡錘體裝配檢查點的1 個組件[42]；ASPM 通過調節CDK4 參與肺鱗癌的發生發展[43]。這16 個hub 基因與細胞周期、細胞核分裂、染色體分離等細胞分裂活動相關，但它們能否作為肺腺癌預測、診斷、個體化預防、治療和預后的生物標志物有待進一步實驗驗證。