劉少博 黃 波

（錦州醫(yī)科大學(xué)附屬第一醫(yī)院胸外科，錦州 121000）

肺癌是全球范圍內(nèi)發(fā)病率和病死率最高的惡性腫瘤，在我國(guó)，肺癌在所有男性惡性腫瘤中發(fā)病率和病死率均位列第一，在所有女性惡性腫瘤中發(fā)病率位列第二，僅次于乳腺，病死率則位列第一［1］。肺腺癌（lung adenocarcinoma，LUAD）是目前肺癌最常見(jiàn)的病理類型，目前肺腺癌的發(fā)生率逐年增加，呈現(xiàn)出年輕化的趨勢(shì)，疾病初期癥狀少，發(fā)病迅速，病死率高且預(yù)后差，多數(shù)患者被診斷時(shí)已經(jīng)處于晚期［2-3］。當(dāng)今精準(zhǔn)醫(yī)學(xué)的發(fā)展使基因?qū)用娴闹委煾泳珳?zhǔn)，通過(guò)對(duì)肺腺癌患者進(jìn)行基因檢測(cè)，已發(fā)現(xiàn)最常見(jiàn)的肺癌驅(qū)動(dòng)基因有EGFR、ALK、ROS1 和BRAF［4］。在患小細(xì)胞肺癌的亞洲人中，EGFR 突變率可達(dá)35%～40%，基于此，近年來(lái)基因靶向治療藥物如吉非替尼、厄洛替尼和克唑替尼等廣泛用于臨床治療，免疫治療如免疫檢查點(diǎn)抑制劑PD-1/PD-L1也可通過(guò)對(duì)免疫檢查點(diǎn)的抑制來(lái)治療癌癥，對(duì)患者的生存時(shí)間及生存質(zhì)量有一定的提高［5］。不幸的是，肺腺癌的預(yù)后仍然很差，因此探索新的生物標(biāo)志物和預(yù)后基因成為精密醫(yī)學(xué)時(shí)代的研究趨勢(shì)。

目前，醫(yī)療技術(shù)和分子生物學(xué)技術(shù)都有了很大的發(fā)展，隨著基因組微陣列和高通量測(cè)序技術(shù)的進(jìn)步以及結(jié)合生物信息學(xué)分析為研究腫瘤的發(fā)生發(fā)展提供了有效方法，基因芯片和RNA 測(cè)序的廣泛應(yīng)用也極大豐富了腫瘤的相關(guān)數(shù)據(jù)。一些基于大規(guī)模、全基因組相關(guān)聯(lián)的數(shù)據(jù)庫(kù)也促進(jìn)了新生標(biāo)志物的發(fā)現(xiàn)，最常見(jiàn)的當(dāng)屬GEO和TCGA數(shù)據(jù)庫(kù)，高通量基因表達(dá)數(shù)據(jù)庫(kù)（gene expression omnibus，GEO）由美國(guó)國(guó)立生物技術(shù)信息中心（NCBI）于2000 年創(chuàng)建并維護(hù)的基因表達(dá)數(shù)據(jù)（http：www.ncbi.nlm.nih.gov/geo），收錄全世界高通量基因組數(shù)據(jù)。人類癌癥基因組圖譜（The Cancer Genome Atlas，TCGA）（http：//cancergenome.nih.gov/），包括33 種腫瘤的臨床隨訪數(shù)據(jù)和基因組學(xué)數(shù)據(jù)，因?yàn)椴煌矓?shù)據(jù)庫(kù)的內(nèi)容或多或少存在一定異質(zhì)性，綜合多個(gè)數(shù)據(jù)庫(kù)進(jìn)行生物信息學(xué)分析便可以減少樣本的異質(zhì)性和平臺(tái)差異性，將多個(gè)平臺(tái)的不同微陣列數(shù)據(jù)進(jìn)行聯(lián)合分析也可以獲得更加豐富的臨床數(shù)據(jù)。本研究通過(guò)一定篩選條件從GEO 數(shù)據(jù)庫(kù)下載3 個(gè)數(shù)據(jù)集，結(jié)合TCGA 肺腺癌數(shù)據(jù)集進(jìn)行差異基因的篩選，并對(duì)差異基因進(jìn)行加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析、富集分析、表達(dá)差異分析、生存分析等，為探討肺腺癌預(yù)后相關(guān)基因的篩選提供理論依據(jù)。

1 材料與方法

1.1 芯片數(shù)據(jù)獲取 在GEO 數(shù)據(jù)庫(kù)的檢索框里輸入關(guān)鍵詞“l(fā)ung cancer”“l(fā)ung adenocarcinoma”，條件為“homo sapiens”“expression profiling by array”。篩選標(biāo)準(zhǔn)：①標(biāo)本為L(zhǎng)UAD 組織和對(duì)應(yīng)的癌旁組織；②每個(gè)芯片數(shù)據(jù)集都包含MRNA 且數(shù)量不少于25 對(duì)，從其中選出3 組符合標(biāo)準(zhǔn)的基因表達(dá)譜數(shù)據(jù)（GSE43458、GSE27262、GSE10072）［6］。另外在TCGA數(shù)據(jù)庫(kù)中，選擇數(shù)據(jù)類別為轉(zhuǎn)錄組數(shù)據(jù)（transcriptome profiling）和原始數(shù)據(jù)（raw counts），包括535 個(gè)原發(fā)性肺腺癌樣本和59 個(gè)正常樣本（表1），然后從TCGA下載533例包括性別、年齡、生存時(shí)間、生存狀態(tài)、病理分期等與之對(duì)應(yīng)的臨床信息用于后續(xù)分析。

表1 基因芯片基本信息Tab.1 Basic information of gene chip

1.2 數(shù)據(jù)預(yù)處理和差異表達(dá)基因的篩選 對(duì)TCGA及GEO 的數(shù)據(jù)集均使用R 軟件進(jìn)行處理，如果多個(gè)探針對(duì)應(yīng)同一個(gè)基因，則表達(dá)的平均值被認(rèn)為是該基因的表達(dá)水平，采用Benjamini-Hochberg 方法調(diào)整P值，以控制錯(cuò)誤發(fā)現(xiàn)率（FDR）。首先對(duì)GEO 數(shù)據(jù)集進(jìn)行預(yù)處理，采用Perl 語(yǔ)言對(duì)3 組原始數(shù)據(jù)集（GSE43458、GSE27263、GSE10072）進(jìn)行矩陣的合并，接下來(lái)對(duì)合并后的原始數(shù)據(jù)采用R 語(yǔ)言中Bioconductor 的R 包“Affy”中魯棒多芯片平均 算 法（RMA）（robust multichip average algorithm）進(jìn)行背景矯正、標(biāo)準(zhǔn)化和以2 為底的對(duì)數(shù)轉(zhuǎn)換，然后利用R 軟件包SVA 的combat 函數(shù)進(jìn)行批次矯正，對(duì)去除批次效應(yīng)前后的數(shù)據(jù)表達(dá)分別進(jìn)行箱線圖的繪制［7］。從TCGA 數(shù)據(jù)庫(kù)中下載的原始數(shù)據(jù)去除重復(fù)基因及其表達(dá)量之后，利用R 軟件edgeR 包的CPM 函數(shù)進(jìn)行數(shù)據(jù)的矯正及標(biāo)準(zhǔn)化處理［8］，刪除CPM（每百萬(wàn)堿基中每個(gè)轉(zhuǎn)錄本count 值）均值＜1 的樣本，并進(jìn)行以2為底數(shù)的轉(zhuǎn)換。然后對(duì)以上預(yù)處理過(guò)的兩組數(shù)據(jù)集分別使用R 軟件包Limma 篩選差異基因［9］，篩選標(biāo)準(zhǔn)為：|log2（fold-change）|＞1 以及矯正后P值（false discovery rate，F(xiàn)DR）＜0.05，對(duì)TCGA 和GEO 篩選出的差異基因分別利用R 語(yǔ)言“gplots”程序包中的“heatmap.2”函數(shù)對(duì)正常肺組織樣本和腫瘤樣本繪制聚類熱圖。利用火山圖來(lái)觀察FDR 和差異變化倍數(shù)之間的關(guān)系，并對(duì)求出的GEO和TCGA數(shù)據(jù)庫(kù)的差異基因通過(guò)在線網(wǎng)頁(yè)工具繪制韋恩圖（http：//bioinformatics.psb.ugent.be/webtools/Venn），獲取兩者共同表達(dá)的上調(diào)和下調(diào)的差異基因。

1.3 差異基因的富集分析 為了探索肺腺癌發(fā)生發(fā)展的機(jī)制，利用基因功能分析（基因本體論，gene ontology，GO）與通路分析（京都基因與基因組百科全書(shū)Kyoto encyclopedia of genes and genomes，KEGG）對(duì)基因產(chǎn)物功能進(jìn)行詳細(xì)的生物學(xué)注釋和描述。GO涵蓋了分子生物學(xué)功能（molecular function，MF）、細(xì)胞學(xué)組分（cellular components，CC）和生物學(xué)過(guò)程（biological process，BP），通過(guò)富集分析的形式全面概括了給定基因的功能信息［10］。KEGG 是整合基因組、化學(xué)和系統(tǒng)功能信息并從基因和分子網(wǎng)絡(luò)方面系統(tǒng)性分析基因功能的一個(gè)數(shù)據(jù)庫(kù)，通常用于識(shí)別功能和代謝途徑［11］。DAVID 在線分析平臺(tái)（https：//david.ncifcrf.gov/）是一個(gè)生物信息數(shù)據(jù)庫(kù)，為大規(guī)模的基因或蛋白列表提供系統(tǒng)綜合的生物功能注釋，用于從多個(gè)基因和蛋白質(zhì)集合中提取比較有意義的生物信息，使用DAVID 分別分析了上調(diào)和下調(diào)的基因在GO中的注釋并利用KEGG進(jìn)行通路分析，設(shè)定P＜0.05為顯著性基因富集。

1.4 蛋白互作網(wǎng)絡(luò)的構(gòu)建與分析 String 數(shù)據(jù)庫(kù)（http：//string-db.org/）是一個(gè)搜索已知蛋白質(zhì)之間和預(yù)測(cè)蛋白質(zhì)之間相互作用的數(shù)據(jù)庫(kù)［12］，在STRING 在線數(shù)據(jù)庫(kù)中對(duì)TCGA 和GEO 數(shù)據(jù)庫(kù)共有的DEGS 進(jìn)行了蛋白質(zhì)-蛋白質(zhì)相互作用（proteinprotein interaction，PPI）網(wǎng)絡(luò)分析，并將置信分?jǐn)?shù)＞0.9 設(shè)置為截止標(biāo)準(zhǔn)，然后將PPI 網(wǎng)絡(luò)的信息導(dǎo)入Cytoscape 3.6.0（http：//www.cytoscape.org/）中使其可視化。Cytoscape 作為生物信息分析的開(kāi)源軟件工具之一，用于可視化探索由蛋白質(zhì)、基因和其他類型相互作用組成的生物互助網(wǎng)絡(luò)，是生物信息學(xué)研究的重要工具之一［13］。使用Cytoscape 的插件Cytohubba 其中的5 種方法從DEGS 的PPI 網(wǎng)絡(luò)中篩選中樞基因，包括EPC（邊緣滲透成分）、MCC（最大團(tuán)中心性）、MNC（最大鄰域成分）、Degree（節(jié)點(diǎn)連接度）和Closeness（節(jié)點(diǎn)連接緊密度），挑選在5種計(jì)算指標(biāo)得分均出現(xiàn)的基因作為中樞基因。另外通過(guò)插件MCODE（molecular complex detection）發(fā)掘肺腺癌PPI網(wǎng)絡(luò)中不同功能的基因模塊，篩選標(biāo)準(zhǔn)設(shè)定為：Degree Cutoff=2、Node Score Cutoff=0.2、K-Core=2、Max Depth=100。篩選出其中最顯著的模塊，MCODE是通過(guò)蛋白質(zhì)復(fù)合物聚類找到緊密連接的部分，從而篩選出差異基因的基因功能模塊［14］，隨后運(yùn)用DAVID 對(duì)最顯著模塊中的基因進(jìn)行GO 和KEGG分析。

1.5 TCGA 數(shù)據(jù)集DEGS 的WGCNA 分析及關(guān)鍵基因的確定 加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（WGCNA）是從全基因組表達(dá)中理解基因功能和基因關(guān)聯(lián)的一種重要方法，可用于檢測(cè)高度相關(guān)基因的共表達(dá)模塊（module-membership，MM）以及與臨床特征相關(guān)的模塊（gene-significance，GS），為預(yù)測(cè)共表達(dá)基因的功能和發(fā)現(xiàn)在人類疾病中起關(guān)鍵作用的基因提供了很好的見(jiàn)解［15-17］。此外，轉(zhuǎn)錄組學(xué)中另一個(gè)強(qiáng)大的分析是差異基因表達(dá)分析，它為研究基因組調(diào)控的分子機(jī)制和發(fā)現(xiàn)實(shí)驗(yàn)組與對(duì)照組之間表達(dá)水平的定量變化提供了方法，這種基因表達(dá)的差異可以發(fā)現(xiàn)特定疾病的潛在生物標(biāo)志物［18］。因此，采用兩種方法，將WGCNA 和差異基因表達(dá)分析的結(jié)果結(jié)合起來(lái)，可以高度提高相關(guān)基因的識(shí)別能力。使用R 軟件WGCNA 包［19］對(duì)TCGA 數(shù)據(jù)集的差異基因構(gòu)建共表達(dá)網(wǎng)絡(luò)，首先計(jì)算差異基因各個(gè)基因之間的Pearson 系數(shù)使其轉(zhuǎn)化為相似矩陣，通過(guò)WGCNA包的pick soft threshold 函數(shù)自動(dòng)進(jìn)行網(wǎng)絡(luò)拓?fù)浞治鲞x擇軟閾值β，β 可以強(qiáng)調(diào)基因之間強(qiáng)弱相關(guān)性。確定β 后相似矩陣轉(zhuǎn)化為鄰接矩陣，再將鄰接矩陣轉(zhuǎn)換為拓?fù)渲丿B矩陣（TOM），設(shè)置模塊最小基因數(shù)為50，剪切高度為0.25，通過(guò)層次聚類使表達(dá)相近的基因置于同一基因模塊，并將閾值設(shè)置為20 000以消除異常值，利用動(dòng)態(tài)混合切割方法，將表達(dá)模式類似的基因分到不同的模塊中。得到這些數(shù)據(jù)后，計(jì)算基因模塊和表型（癌組織和正常樣本）的Pearson 相關(guān)系數(shù)，選擇與腫瘤發(fā)生密切相關(guān)的基因模塊，用GO和KEGG分析挖掘目標(biāo)模塊所參與的生物學(xué)功能，然后利用基因和模塊的相關(guān)性和基因與臨床性狀的相關(guān)性進(jìn)行顯著模塊核心基因的挖掘。如果模塊中一個(gè)基因同時(shí)具有較大的MM和GS，則該基因被認(rèn)為是模塊中的核心基因，將MM＞0.7 和GS＞0.35 定義為候選的核心基因，然后利用Cytohubba 篩選的中樞基因與模塊篩選的核心基因取交集，并將交集中的基因定義為最終的關(guān)鍵基因。

1.6 關(guān)鍵基因的生存分析及差異分析的表達(dá)Kaplan-Meier plotter 是基于EGA、TCGA 和GEO 數(shù)據(jù)庫(kù)評(píng)估大量基因?qū)ι嬗绊懙某Ｓ镁W(wǎng)站工具，利用Kaplan-Meier plotter 驗(yàn)證9 個(gè)關(guān)鍵基因與肺癌患者預(yù)后總生存率的關(guān)系。GEPIA（http：//gepia.cancer-pku.cn/）是一個(gè)在線的基因表達(dá)譜動(dòng)態(tài)數(shù)據(jù)分析數(shù)據(jù)庫(kù)，可用于分析癌癥和正常組織之間的表達(dá)差異以及總生存率，進(jìn)一步驗(yàn)證關(guān)鍵基因的mRNA表達(dá)水平［20］。HPA（https：//www.proteinatlas.org/）提供了大量人類蛋白質(zhì)的表達(dá)譜，呈現(xiàn)為大多數(shù)人類組織的免疫組織化學(xué)（IHC）等實(shí)驗(yàn)數(shù)據(jù)的蛋白質(zhì)表達(dá)譜數(shù)據(jù)庫(kù)［21］。用免疫組化法（IHC）從人蛋白圖譜數(shù)據(jù)庫(kù)（HPA）中檢測(cè)肺腺癌與正常組織之間生存相關(guān)基因的蛋白表達(dá)。

1.7 預(yù)后模型的構(gòu)建和驗(yàn)證 Cox 回歸模型是一種以生存時(shí)間和生存結(jié)局為變量，可同時(shí)分析多種因素對(duì)生存期影響的半?yún)?shù)回歸模型，將從TCGA網(wǎng)站下載的患者臨床數(shù)據(jù)，去除總生存率缺少的數(shù)據(jù)后將表達(dá)和生存數(shù)據(jù)合并，然后將數(shù)據(jù)集隨機(jī)平均分為訓(xùn)練集和驗(yàn)證集，使用訓(xùn)練集建立模型并在驗(yàn)證集進(jìn)行驗(yàn)證，將篩選的GEO 樣本和TCGA 樣本均存在差異的479 個(gè)基因，利用訓(xùn)練集中的生存數(shù)據(jù)使用R 軟件“survival”生存分析軟件包進(jìn)行單變量Cox比例風(fēng)險(xiǎn)回歸分析得到與預(yù)后顯著相關(guān)的基因（P＜0.01）［22］，然后通過(guò)glmnet 程序包以生存狀態(tài)為應(yīng)變量，篩選出的基因表達(dá)值作為反應(yīng)變量進(jìn)行1 000 次Lasso 回歸分析對(duì)基因個(gè)數(shù)進(jìn)行降維處理，從而降低模型的誤差獲得廣義的線性模型［23］，而后進(jìn)行多因素Cox 比例風(fēng)險(xiǎn)回歸分析，獲得風(fēng)險(xiǎn)基因并構(gòu)建風(fēng)險(xiǎn)預(yù)后模型［24］。該模型使用疾病風(fēng)險(xiǎn)評(píng)分作為預(yù)后狀態(tài)的預(yù)測(cè)因子，疾病風(fēng)險(xiǎn)評(píng)分由多變量Cox 比例風(fēng)險(xiǎn)回歸分析的參數(shù)β 和樣本中每個(gè)基因的表達(dá)量確定［25］。利用預(yù)后模型分別對(duì)驗(yàn)證集和訓(xùn)練集進(jìn)行風(fēng)險(xiǎn)評(píng)分的計(jì)算，依據(jù)風(fēng)險(xiǎn)指數(shù)的中位數(shù)分別將驗(yàn)證集和訓(xùn)練集分為高、低風(fēng)險(xiǎn)組，結(jié)合生存信息繪制生存曲線得出高、低風(fēng)險(xiǎn)表達(dá)生存狀況，評(píng)價(jià)模型預(yù)測(cè)效果是否顯著（P＜0.05），在這個(gè)過(guò)程中使用的統(tǒng)計(jì)方法是對(duì)數(shù)秩檢驗(yàn)。使用R軟件“survival ROC”包計(jì)算時(shí)間依賴的受試者工作曲線（ROC 曲線）評(píng)估回歸模型在1 年、3 年、5 年生存期的預(yù)測(cè)能力［26］，AUC＞0.5時(shí)而且越接近1，預(yù)后越好。利用生存時(shí)間和基因風(fēng)險(xiǎn)模型分別繪制散點(diǎn)圖和高低風(fēng)險(xiǎn)熱圖，并通過(guò)驗(yàn)證集驗(yàn)證回歸模型在預(yù)測(cè)肺腺癌患者生存預(yù)后的價(jià)值和穩(wěn)定性，以此來(lái)證明得到的風(fēng)險(xiǎn)評(píng)分是合理的。此外，為了使模型更有效地應(yīng)用于臨床過(guò)程，將臨床信息（性別、年齡、分期）納入預(yù)后模型，剔除臨床資料缺失的樣本，共獲得480份樣本，利用這些樣本風(fēng)險(xiǎn)評(píng)分和臨床信息進(jìn)行列線圖的繪制。

2 結(jié)果

2.1 篩選差異表達(dá)基因 經(jīng)過(guò)對(duì)3 組GEO 基因芯片進(jìn)行合并及數(shù)據(jù)標(biāo)準(zhǔn)化之后共有104個(gè)正常肺樣本和163 個(gè)肺腺癌樣本，進(jìn)行批次矯正用以消除GSE43458、GSE27262 和GSE10072 的批次效應(yīng)（圖1A），然后在合并后的GEO 微陣列數(shù)據(jù)集中得到337 個(gè)顯著下調(diào)基因和154 個(gè)顯著上調(diào)基因（圖1B、C），從包含59 個(gè)正常樣本和535 個(gè)肺腺癌樣本的TCGA 數(shù)據(jù)集中得到2 101 個(gè)下調(diào)基因1 481 個(gè)上調(diào)基因（圖1E、F）。將兩個(gè)數(shù)據(jù)集取交集得到148 個(gè)上調(diào)基因和331個(gè)下調(diào)基因（圖1D）。

圖1 差異表達(dá)基因熱圖及火山圖Fig.1 Heatmap and volcano map of DEGs

2.2 差異基因的GO 及KEGG 分析 將篩選出的479 個(gè)差異基因通過(guò)DAVID 進(jìn)行功能和途徑的富集，利用GO 分析，將所有差異基因同時(shí)富集到BP、CC、MF 這3 種生物學(xué)關(guān)系中，結(jié)果表明：148 個(gè)上調(diào)的差異基因主要參與核分裂、有絲分裂姐妹染色單體分離、核仁染色體分離以及細(xì)胞外基質(zhì)組織等生物過(guò)程，其產(chǎn)物主要參與有絲分裂的紡錘體、膠原三聚體復(fù)合體、中間體、染色體上的著絲粒等細(xì)胞組分，發(fā)揮絲氨酸內(nèi)肽酶活性、血小板衍生生長(zhǎng)因子結(jié)合、蛋白酶結(jié)合、金屬內(nèi)肽酶活性、糖胺聚糖結(jié)合、絲氨酸水解酶及肽酶活性等生物學(xué)分子功能（圖2A）。涉及的信號(hào)通路主要包括：細(xì)胞周期、蛋白質(zhì)的消化吸收、ECM-受體相互作用、P53 信號(hào)通路、卵母細(xì)胞的減數(shù)分裂、孕酮介導(dǎo)的卵母細(xì)胞成熟、IL-17 和松弛素信號(hào)通路等（圖2B）。331 個(gè)下調(diào)的DEGs 涉及的生物學(xué)過(guò)程主要包括：血管系統(tǒng)發(fā)育生成的調(diào)節(jié)、阿米巴樣細(xì)胞遷移、細(xì)胞-基質(zhì)黏附、組織和上皮細(xì)胞的遷移、負(fù)調(diào)控生長(zhǎng)以及對(duì)糖皮質(zhì)激素的反應(yīng)；涉及的細(xì)胞學(xué)組分主要包括：含膠原蛋白的細(xì)胞外基質(zhì)、細(xì)胞-細(xì)胞連接、膜筏、質(zhì)膜的外側(cè)、黏著斑、細(xì)胞-底物連接、血小板α 顆粒等；參與的分子生物學(xué)功能主要包括：酰胺結(jié)合、肽結(jié)合、糖胺聚糖結(jié)合、細(xì)胞因子結(jié)合、生長(zhǎng)因子結(jié)合、跨膜受體蛋白激酶活性、淀粉樣蛋白-β 結(jié)合、轉(zhuǎn)化生長(zhǎng)因子-β 結(jié)合、跨膜受體蛋白絲氨酸/蘇氨酸激酶活性（圖2C）。KEGG 信號(hào)通路主要包括細(xì)胞因子-細(xì)胞因子受體相互作用、細(xì)胞黏附分子、血管平滑肌收縮、補(bǔ)體和凝血級(jí)聯(lián)、cAMP信號(hào)通路等（圖2D）。

圖2 差異基因的GO和KEGG富集分析Fig.2 Enrichment analysis of differentially expressed genes by GO and KEGG

2.3 蛋白互助網(wǎng)絡(luò)的構(gòu)建及中樞基因鑒定 基于String 數(shù)據(jù)庫(kù)利用Cytoscape 軟件對(duì)差異表達(dá)基因進(jìn)行PPI網(wǎng)絡(luò)的構(gòu)建（圖3A），包括478個(gè)節(jié)點(diǎn)和816個(gè)邊緣，首先使用5 種方法分析前30 位基因，取共有的基因?yàn)長(zhǎng)UAD 的中樞基因，得到的19 個(gè)中樞基因分別為：ASPM、AURKA、CENPF、CEP55、DLGAP5、KIF4A、MELK、NCAPG、NDC80、NEK2、NUSAP1、PBK、PRC1、PTTG1、RRM2、TOP2A、TTK、KIF20A 和TPX2（表2）。利用Cytoscape 的插件MCODE 獲得最顯著的模塊（圖3B），可見(jiàn)中樞基因都位于最顯著模塊而且都為上調(diào)基因。GO 富集分析表明，在生物過(guò)程中，該模塊的基因主要在細(xì)胞分裂和有絲分裂核分裂以及染色體分離中富集；細(xì)胞組分分析表明，基因在紡錘體、染色體、中間體中明顯富集；分子功能分析表明，這些基因主要參與ATP 和部分蛋白質(zhì)的結(jié)合（圖3C）。KEGG 分析表明這些基因主要參與細(xì)胞周期和卵母細(xì)胞減數(shù)分裂（圖3D）。

表2 多種CytoHubba方法中樞基因的排序Tab.2 Sequencing of central genes by various cytohubba methods

圖3 蛋白互助網(wǎng)絡(luò)的可視化及最顯著模塊的分析Fig.3 Visualization of PPI network and analysis of most significant modules

2.4 關(guān)鍵基因的篩選 利用TCGA 數(shù)據(jù)集中提取的3 582個(gè)差異基因表達(dá)譜，選取軟閾值β=3建立基因調(diào)控網(wǎng)絡(luò)（圖4D），結(jié)果顯示綠松石色模塊與正常樣本表型相關(guān)系數(shù)最大為0.82，藍(lán)色模塊與肺腺癌樣本表型相關(guān)系數(shù)最大為0.54（圖4A），另外根據(jù)各模塊間的Pearson 相關(guān)系數(shù)也發(fā)現(xiàn)藍(lán)色和綠松石色一致性最大，因此選擇藍(lán)色模塊為目的模塊，模塊中MM＞0.7 和GS＞0.35 的基因定義為核心基因，綠松石色和藍(lán)色基因分布如圖4B、C。另外，經(jīng)過(guò)cytoscape篩選的19個(gè)中樞基因均位于藍(lán)色模塊，GO（圖4E）和KEGG（圖4F）分析結(jié)果表明，藍(lán)色模塊與有絲分裂、染色體分離、細(xì)胞周期、DNA 的轉(zhuǎn)錄復(fù)制、p53 信號(hào)通路以及卵母細(xì)胞的減數(shù)分裂等關(guān)系更為密切，可能與癌細(xì)胞過(guò)度增殖有關(guān)，其模塊內(nèi)的基因可能對(duì)藥物開(kāi)發(fā)有重要的作用。核心基因和PPI網(wǎng)絡(luò)中識(shí)別的中樞基因共有的基因作為最終的關(guān)鍵基因，分別為ASPM、CEP55、DLGAP5、KIF4A、MELK、NEK2、RRM2、TOP2A、TPX2。

圖4 WGCNA分析與最顯著模塊基因富集分析Fig.4 WGCNA analysis and most significant module gene enrichment analysis

2.5 關(guān)鍵基因的預(yù)后分析及表達(dá)差異 在PPI 網(wǎng)絡(luò)和WGCNA 共同篩選獲得了9 個(gè)關(guān)鍵基因，這些基因可能在肺腺癌的發(fā)生發(fā)展進(jìn)程中起關(guān)鍵作用，利用Kaplan-Meier 曲線分析得出這9 個(gè)關(guān)鍵基因?qū)颊叩目偵鏁r(shí)間有著顯著影響（P＜0.01，圖5），為了進(jìn)一步驗(yàn)證，利用人類蛋白圖譜數(shù)據(jù)庫(kù)獲得癌癥和正常組織中9種基因蛋白水平的免疫組織化學(xué)染色圖像，結(jié)果表明除ASPM 無(wú)相關(guān)數(shù)據(jù)之外，其余基因在LUAD 中均有顯著上調(diào)（圖6A），另外通過(guò)GEPIA 數(shù)據(jù)庫(kù)分析上述基因在基因水平上肺腺癌與癌旁樣本之間均存在顯著差異且均在肺腺癌組織中呈現(xiàn)高表達(dá)狀態(tài)（圖6B），進(jìn)一步說(shuō)明這些基因在肺腺癌的發(fā)生發(fā)展中有一定作用，提示這些基因可能成為預(yù)后的分子標(biāo)志物和治療靶點(diǎn)。

圖5 9個(gè)hub基因的總生存率（OS）分析Fig.5 Overall survival（OS）analysis of 9 hub genes

圖6 驗(yàn)證核心基因表達(dá)水平Fig.6 Validate expression level of critical genes

2.6 預(yù)后模型的構(gòu)建 將表達(dá)和生存數(shù)據(jù)合并后的494 個(gè)TCGA 數(shù)據(jù)集樣本分為訓(xùn)練集和驗(yàn)證集，為保證能預(yù)測(cè)出有效的預(yù)后模型，首先使用訓(xùn)練集的生存數(shù)據(jù)對(duì)479 個(gè)差異基因進(jìn)行單因素Cox 比例風(fēng)險(xiǎn)回歸分析，共鑒定出34個(gè)對(duì)預(yù)后有顯著影響的基因（P＜0.01），然后通過(guò)Lasso 回歸分析，可以得到19 個(gè)基因進(jìn)行后續(xù)分析，進(jìn)一步使用多變量Cox 比例風(fēng)險(xiǎn)回歸分析，共獲得12 個(gè)風(fēng)險(xiǎn)基因（圖7A），分別 為CA4、ENO1、FBLN5、FZD4、INAVA、NEK2、RRAS、SEMA5A、TIMP1、TMPRSS11E、EFNB2、AKAP12，進(jìn)行風(fēng)險(xiǎn)預(yù)后模型的構(gòu)建，即Risk score=（0.001×ENO1）-（0.208×CA4）+（0.006×FBLN5）+（0.041×FZD4）+（0.055×INAVA）+（0.075×NEK2）+（0.006×RRAS）+（0.083×SEMA5A）+（0.001×TIMP1）+（0.013×TMPRSS11E）+（0.018×EFNB2）+（0.006×AKAP12），通過(guò)風(fēng)險(xiǎn)得分算出高低風(fēng)險(xiǎn)組，分別在訓(xùn)練集和驗(yàn)證集進(jìn)行生存分析，得出低風(fēng)險(xiǎn)組的患者生存狀況明顯優(yōu)于高風(fēng)險(xiǎn)組（圖7B、C）。使用ROC 曲線對(duì)模型的預(yù)測(cè)性能進(jìn)行評(píng)估，結(jié)果可見(jiàn)：訓(xùn)練集中使用ROC 曲線對(duì)風(fēng)險(xiǎn)模型的預(yù)測(cè)AUC 分別為0.785、0.748、0.771（圖7D～F），驗(yàn)證集中得出AUC 分別為0.736、0.706、0.621（圖7G～I(xiàn)），另外可從生存時(shí)間和風(fēng)險(xiǎn)評(píng)分繪制的散點(diǎn)圖中看出，隨著風(fēng)險(xiǎn)得分的增加，死亡的患者也增加，存活時(shí)間相對(duì)減少，由此可見(jiàn)模型有相對(duì)較好的預(yù)測(cè)能力（圖8）。

圖7 基因風(fēng)險(xiǎn)模型的構(gòu)建Fig.7 Construction of gene risk model

圖8 風(fēng)險(xiǎn)模型得分與生存時(shí)間、臨床信息的關(guān)系Fig.8 Relationship between risk model score and survival time and clinical information

3 討論

腫瘤的發(fā)生發(fā)展涉及多個(gè)環(huán)節(jié)、因素和階段，而細(xì)胞周期的改變是驅(qū)使細(xì)胞向惡性轉(zhuǎn)化的關(guān)鍵一步，只有突破細(xì)胞周期的調(diào)控才可以抑制腫瘤的發(fā)生發(fā)展。隨著高通量測(cè)序技術(shù)和基因微陣列的高速發(fā)展，可以檢測(cè)到一些基因的改變與疾病的關(guān)系，為疾病的診斷及預(yù)后提供一定的理論幫助，由于不同平臺(tái)或者數(shù)據(jù)集中小樣本會(huì)存在局限性，本文通過(guò)多個(gè)數(shù)據(jù)集進(jìn)行整合，分別通過(guò)PPI 網(wǎng)絡(luò)和WGCNA 共表達(dá)分析進(jìn)行關(guān)鍵基因的挖掘，PPI 網(wǎng)絡(luò)是基于互助的蛋白質(zhì)網(wǎng)絡(luò)，WGCNA 是基于基因之間的相關(guān)性構(gòu)造的網(wǎng)絡(luò)，兩者相結(jié)合為新的預(yù)后基因的篩選提供了巨大潛能，首先對(duì)3 組GEO 數(shù)據(jù)集和TCGA 數(shù)據(jù)集進(jìn)行標(biāo)準(zhǔn)化處理，之后將3 組GEO數(shù)據(jù)集進(jìn)行合并和批次矯正。通過(guò)生物信息學(xué)分析，共得到479 個(gè)差異基因（上調(diào)148 個(gè)、下調(diào)331 個(gè)），GO 分析表明主要與細(xì)胞分裂增殖、周期調(diào)控、減數(shù)分裂和有絲分裂核分裂以及染色體分離等生物過(guò)程相關(guān)，主要參與組成紡錘體、染色體、中間體等細(xì)胞組分并參與ATP 和部分蛋白質(zhì)的結(jié)合；KEGG 分析表明這些基因主要參與細(xì)胞周期和卵母細(xì)胞減數(shù)分裂。最終確定了9 個(gè)與LUAD 患者預(yù)后明顯相關(guān)的關(guān)鍵基因，分別為ASPM、CEP55、DLGAP5、KIF4A、MELK、NEK2、RRM2、TOP2A、TPX2。

細(xì)胞增殖是癌癥的特征，而惡性表型特征不受控制的基礎(chǔ)就是細(xì)胞周期的去調(diào)控，癌癥遺傳學(xué)已經(jīng)表明，生長(zhǎng)信號(hào)網(wǎng)絡(luò)中的過(guò)度激活突變，加上腫瘤抑制蛋白功能的喪失，推動(dòng)了癌基因的增殖，細(xì)胞周期引擎位于復(fù)雜的致癌信號(hào)網(wǎng)絡(luò)的匯合點(diǎn)下游，是腫瘤診斷和治療的重要靶點(diǎn)，它的失控是所有癌癥細(xì)胞異常增殖的核心［27］。9 個(gè)關(guān)鍵基因多通過(guò)紡錘體和中心體形成來(lái)參與影響細(xì)胞周期的進(jìn)程，在人類多種惡性腫瘤中發(fā)現(xiàn)了異常表達(dá)水平，有可能成為抗癌治療的靶點(diǎn)。人類異常紡錘體樣小頭畸形相關(guān)蛋白ASPM 產(chǎn)物多位于紡錘體和中心體，主要使細(xì)胞有絲分裂時(shí)紡錘體向兩極運(yùn)動(dòng)，并且維持細(xì)胞質(zhì)的均等分裂［28］，在多種癌癥中高表達(dá)。相關(guān)研究顯示ASPM 在膠質(zhì)母細(xì)胞瘤、前列腺癌中的表達(dá)水平與腫瘤的病理分級(jí)及臨床分期密切相關(guān)［29］。最新研究顯示ASPM 在肺腺癌中高表達(dá)，并與生存率、臨床分期及預(yù)后相關(guān)［30］。中心體相關(guān)蛋白CEP55 主要功能為錨定微管聚合相關(guān)蛋白和參與紡錘體形成，并與中心體相偶聯(lián)，磷酸化后發(fā)揮調(diào)控細(xì)胞周期的作用，達(dá)到對(duì)細(xì)胞增殖的調(diào)控［31-32］，研究表明CEP55 的高表達(dá)可以促進(jìn)癌癥的增殖、遷移和侵襲，例如乳腺癌，前列腺癌，腎癌等［33-37］。JING等［38］發(fā)現(xiàn)CEP55在非小細(xì)胞肺癌組織中的表達(dá)顯著增加，并且其過(guò)度表達(dá)與患者的不良預(yù)后相關(guān)。DLGAP5 是一種有絲分裂紡錘體蛋白，促進(jìn)微管蛋白聚合物的形成，在紡錘體組配中起重要作用，可作為信號(hào)分子具有重要的生物學(xué)功能［39-40］。BRANCHI 等［41］研究顯示，DLGAP5 的下調(diào)導(dǎo)致結(jié)直腸癌的侵襲和遷移潛能顯著降低。染色體相關(guān)驅(qū)動(dòng)蛋白KIF4A 是一種基于微管的運(yùn)動(dòng)蛋白，是染色體濃縮和分離機(jī)制的重要組成部分，在有絲分裂的多個(gè)步驟中發(fā)揮作用，并對(duì)調(diào)節(jié)后期紡錘體、胞質(zhì)分裂、中間帶形成和胞質(zhì)分離期間染色體的完整性發(fā)揮重要作用，腫瘤中高表達(dá)可增強(qiáng)肝細(xì)胞癌、口腔癌和乳腺癌的增殖和侵襲［42-46］。相關(guān)研究顯示KIF4A 可作為肺癌的預(yù)后生物標(biāo)志物和治療靶點(diǎn)［47］。

MELK 是一種細(xì)胞周期依賴性的絲/蘇氨酸蛋白激酶，在有絲分裂期間參與細(xì)胞周期、胞質(zhì)分裂、mRNA 剪接和細(xì)胞凋亡，是治療多種癌癥的理想治療靶點(diǎn)，在癌細(xì)胞存活中起著不可或缺的作用［48-49］。其高表達(dá)與人類星形細(xì)胞瘤和前列腺癌的惡性程度相關(guān)并且與乳腺癌患者的不良預(yù)后相關(guān)［50-51］。目前研究表明MELK是小細(xì)胞肺癌一個(gè)有前途的治療靶點(diǎn)，其抑制劑OTS167 可作為一類新的抗SCLC 藥物進(jìn)行臨床評(píng)估［52］。NEK2 是位于中心體的絲氨酸/蘇氨酸激酶，通過(guò)參與有絲分裂中心體的復(fù)制和紡錘體的裝配對(duì)細(xì)胞的分裂增殖進(jìn)行調(diào)節(jié)［53-54］。表達(dá)失調(diào)會(huì)造成染色體不穩(wěn)定（CIN）和非整倍體，這也是許多腫瘤的標(biāo)志性變化［55-56］。據(jù)報(bào)道，NEK2表達(dá)增加與腫瘤進(jìn)展有關(guān)，在多種腫瘤中顯著表達(dá)并對(duì)預(yù)后產(chǎn)生不良影響，如胰腺導(dǎo)管腺癌、前列腺癌，結(jié)腸癌［57-59］。ZHONG 等［60］研究表明NEK2 可能是非小細(xì)胞肺癌患者預(yù)后不良的更有效的腫瘤增殖標(biāo)志物。RRM2 是DNA 合成和修復(fù)的限速酶，是細(xì)胞凋亡的重要調(diào)控基因，已被報(bào)道是膠質(zhì)瘤中具有功能意義的潛在預(yù)后生物標(biāo)志物［61］，在非小細(xì)胞肺癌和細(xì)胞系中異常上調(diào)預(yù)示著預(yù)后不良，有研究顯示敲除RRM2通過(guò)內(nèi)在途徑導(dǎo)致頭頸鱗狀細(xì)胞癌和非小細(xì)胞肺癌細(xì)胞系的凋亡［62］。拓?fù)洚悩?gòu)酶IiαTOP2A 是在轉(zhuǎn)錄過(guò)程中控制和改變DNA 拓?fù)錉顟B(tài)的酶，參與了多種惡性腫瘤細(xì)胞的有絲分裂過(guò)程［63］。miR-144-3p 通過(guò)靶向TOP2A 抑制膠質(zhì)瘤細(xì)胞的生長(zhǎng)并促進(jìn)其凋亡［64］。在乳腺癌中與erbb2 同時(shí)缺失或擴(kuò)增，很可能是預(yù)測(cè)蒽環(huán)類藥物受益患者亞群的有用標(biāo)志物［65］。TOP2A 的高表達(dá)與非小細(xì)胞肺癌中癌細(xì)胞的增殖和侵襲以及干擾密切相關(guān)［66］。已被廣泛用作NSCLC 的獨(dú)立預(yù)后因子，其高表達(dá)與NSCLC 患者的不良預(yù)后相關(guān)［67］。靶向非洲爪蟾驅(qū)動(dòng)蛋白樣蛋白2TPX2 是一種微管相關(guān)蛋白，參與紡錘體的組裝并維持其結(jié)構(gòu)穩(wěn)定，調(diào)節(jié)有絲分裂的關(guān)鍵點(diǎn)，在多種人類癌癥中過(guò)度表達(dá)，并促進(jìn)癌癥發(fā)展。有報(bào)道顯示在前列腺癌中敲除TPX2 能誘導(dǎo)細(xì)胞周期靜止和凋亡并且降低細(xì)胞的侵襲能力和抑制細(xì)胞的增殖。TPX2 沉默通過(guò)調(diào)節(jié)PI3K/AKT 信號(hào)抑制肺腺癌和肝細(xì)胞癌增殖［68-69］。其高表達(dá)與非小細(xì)胞肺癌的不良預(yù)后有關(guān)，可能為預(yù)后相關(guān)基因［70］。

本研究建立了一個(gè)用于預(yù)測(cè)患者生存率的預(yù)后模型，該模型包含12 個(gè)關(guān)鍵基因，分別為CA4、ENO1、FBLN5、FZD4、INAVA、NEK2、RRAS、SEMA5A、TIMP1、TMPRSS11E、EFNB2、AKAP12，碳酸酐酶ⅳ（CA4）是人類12 種活性同工酶的一種，其低表達(dá)可以促進(jìn)癌細(xì)胞的增殖，據(jù)報(bào)道CA4 是一種新的結(jié)直腸癌腫瘤抑制因子，可以作為結(jié)直腸癌復(fù)發(fā)的獨(dú)立生物標(biāo)志物［71-73］。在模型中系數(shù)最大，說(shuō)明CA4 是LUAD 中一個(gè)非常重要的預(yù)后因素，對(duì)判斷患者預(yù)后具有重要的參考價(jià)值。烯醇酶1（ENO1）作為一種糖酵解酶，在葡萄糖代謝中起著關(guān)鍵作用，并導(dǎo)致許多癌癥的腫瘤進(jìn)展，新的研究證明通過(guò)PI3K/AKT 途徑促進(jìn)非小細(xì)胞肺癌的糖酵解、增殖、遷移和侵襲［74-75］。FBLN5 是Fibulin 蛋白家族成員之一，其表達(dá)水平與肺癌等多種腫瘤的發(fā)生相關(guān)，并能夠影響腫瘤的增殖侵襲及預(yù)后進(jìn)展，因此有可能成為腫瘤診斷新的分子標(biāo)志物。研究發(fā)現(xiàn)FBLN5 能夠通過(guò)特殊機(jī)制調(diào)控腫瘤微環(huán)境從而調(diào)控腫瘤的發(fā)生［76］。FZD4是卷曲基因家族的成員，據(jù)報(bào)道，腫瘤抑制劑miR-493通過(guò)抑制FZD4的表達(dá)來(lái)抑制癌細(xì)胞的生長(zhǎng)和遷移能力［77］，而且已證實(shí)FZD4的敲除導(dǎo)致膀胱癌細(xì)胞遷移和侵襲顯著減少［78］。先天免疫激活因子INAVA，是一種已知為克羅恩病風(fēng)險(xiǎn)基因的蛋白質(zhì)編碼基因［79］，通過(guò)對(duì)肺腺癌患者染色體基因的整體分析，發(fā)現(xiàn)INAVA 在肺腺癌的發(fā)展和進(jìn)展中發(fā)揮重要作用［28］。已被證實(shí)，INAVA 通過(guò)上調(diào)基質(zhì)金屬蛋白酶9的表達(dá)促進(jìn)甲狀腺乳頭狀癌和肝癌侵襲性［80］。RRAS 基因的研究較少，功能、機(jī)制尚未被充分了解。SEMA5A是存在于無(wú)脊椎動(dòng)物和脊椎動(dòng)物中的跨膜蛋白，在多種癌癥中高表達(dá)并且與預(yù)后有關(guān)［81-82］。已被證明在試管內(nèi)能促進(jìn)胃癌細(xì)胞系的遷移和侵襲［83］。癌組織中SEMA5A 在轉(zhuǎn)錄和翻譯水平的下調(diào)與非吸煙女性非小細(xì)胞肺癌患者的低存活率有關(guān)［84］。TIMP1是基質(zhì)金屬蛋白酶的抑制酶，其功能與基質(zhì)金屬蛋白酶（MMPs）相反，有研究表明TIMP1 與大多數(shù)實(shí)體癌的侵襲和轉(zhuǎn)移潛能密切相關(guān)，并且在肺癌中表達(dá)異常并可作為其侵襲轉(zhuǎn)移的潛在的分子標(biāo)志［85］。TMEM185A 可通過(guò)下調(diào)EGFR/AKT 信號(hào)通路，使細(xì)胞在凋亡刺激下對(duì)凋亡敏感，從而抑制食管鱗狀細(xì)胞癌的發(fā)展［86］，但最近發(fā)現(xiàn)在膀胱癌患者中又顯著上調(diào)，與膀胱癌患者的整體存活率顯著相關(guān)［87］。EFNB2 是一種膜錨定配體，屬于受體酪氨酸激酶，能促進(jìn)膠質(zhì)瘤和黑色素瘤的細(xì)胞遷移、侵襲和血管生成，也是卵巢癌和食管鱗狀細(xì)胞癌的不良預(yù)后指標(biāo)［88-91］。EFNB2 基因的敲除抑制結(jié)腸直腸癌細(xì)胞的生長(zhǎng)，逆轉(zhuǎn)了惡性表型并削弱了耐藥性［92］。α-激酶錨蛋白12（AKP12）是一種細(xì)胞支架蛋白，其表達(dá)在多種惡性腫瘤中被抑制，是一種潛在的抑癌基因，不僅能夠抑制腫瘤發(fā)生，而且能夠抑制腫瘤轉(zhuǎn)移。研究證明AKAP12 基因與腫瘤患者的預(yù)后呈正相關(guān)［93-94］；在肺腺癌腫瘤組織中表達(dá)顯著低于正常肺組織，在癌組織中有淋巴結(jié)轉(zhuǎn)移的低于不伴有淋巴結(jié)轉(zhuǎn)移的，而且肺癌臨床分級(jí)越高，AKAP12 的表達(dá)量越低［95］。

綜上所述，本研究利用GEO 數(shù)據(jù)庫(kù)的微陣列數(shù)據(jù)與來(lái)自TCGA 的RNA 測(cè)序數(shù)據(jù)進(jìn)行整合，以確定中樞基因和更重要的關(guān)鍵基因。最后確定了9個(gè)與LUAD 發(fā)病機(jī)制和進(jìn)展相關(guān)的關(guān)鍵基因。這些基因在肺癌的細(xì)胞周期及其異常行為中都起著重要作用，表明這些基因在LUAD 治療以及預(yù)后中具有巨大的潛力。此外，我們進(jìn)行了生存分析，并建立了一個(gè)Cox 比例風(fēng)險(xiǎn)模型來(lái)識(shí)別預(yù)后的生物標(biāo)志物。構(gòu)建了一個(gè)由12 個(gè)基因組成的預(yù)測(cè)總生存率的基因標(biāo)志。這些結(jié)果將為進(jìn)一步研究LUAD 的發(fā)病機(jī)制和藥物治療提供參考。然而，我們的研究所有數(shù)據(jù)為公共數(shù)據(jù)庫(kù)的數(shù)據(jù)分析并使用其他數(shù)據(jù)庫(kù)和臨床數(shù)據(jù)進(jìn)行了驗(yàn)證，但缺乏實(shí)驗(yàn)驗(yàn)證仍然是本研究的局限性，需要進(jìn)一步的實(shí)驗(yàn)研究來(lái)證實(shí)從生物信息學(xué)分析得到的預(yù)測(cè)。