基于癌癥基因組圖譜計劃多組學數(shù)據(jù)構建膠質(zhì)母細胞瘤六基因預后模型

雷常貴，賈學淵，孫文靖

研究報告

基于癌癥基因組圖譜計劃多組學數(shù)據(jù)構建膠質(zhì)母細胞瘤六基因預后模型

雷常貴，賈學淵，孫文靖

哈爾濱醫(yī)科大學醫(yī)學遺傳學研究室，中國遺傳資源保護與疾病防控教育部重點實驗室(哈爾濱醫(yī)科大學)，哈爾濱 150081

膠質(zhì)母細胞瘤(glioblastoma, GBM)是最常見的原發(fā)性顱內(nèi)腫瘤，惡性程度極高，患者預后極差。為了識別GBM預后生物標記物，建立預后模型，本研究通過分析癌癥基因組圖譜計劃(The Cancer Genome Atlas, TCGA)數(shù)據(jù)庫中GBM的表達譜數(shù)據(jù)，篩選出不同生存期GBM患者差異基因。利用GISTIC軟件和Kaplan-Meier (KM)生存分析方法分析TCGA數(shù)據(jù)庫中的GBM拷貝數(shù)變異數(shù)據(jù)，識別影響生存的擴增基因(survival-associated amplified gene, SAG)。取短生存期組上調(diào)基因和SAG兩者的交集基因，進行單因素Cox回歸和迭代Lasso回歸篩選重要候選基因并建立預后模型；計算預后評分，根據(jù)預后評分中位數(shù)將患者分為高風險組和低風險組。用ROC曲線判斷模型的優(yōu)良，KM生存分析高低風險組預后差異，并用GEO、CGGA和Rembrandt數(shù)據(jù)庫3個外部數(shù)據(jù)集進行驗證。多因素Cox回歸分析判斷預后評分的預后獨立性。結果顯示，GBM不同生存期差異分析得到上調(diào)基因426個，下調(diào)基因65個。短生存期組上調(diào)基因與SAG交集得到47個基因。經(jīng)過篩選，最終確定六基因(、、、、、)預后模型。TCGA實驗組和3個外部驗證組模型的ROC曲線下面積均大于0.6，甚至達到0.912。KM分析顯示高低風險組的預后都存在差異(<0.05)。在多因素Cox回歸分析中，六基因預后評分是GBM患者預后的獨立影響因素(<0.05)。通過一系列分析，本研究確立了六基因(、、、、、)的GBM預后模型，模型具有很好的預測能力，可作為預測GBM患者的獨立預后標志物。

膠質(zhì)母細胞瘤；多組學數(shù)據(jù)；六基因組合；預后模型；癌癥基因組圖譜計劃

膠質(zhì)瘤(glioma)是最常見的原發(fā)性顱內(nèi)腫瘤，占所有顱內(nèi)腫瘤的80%[1]。世界衛(wèi)生組織(World Health Organization, WHO) (2016年中樞神經(jīng)系統(tǒng)腫瘤分類總結)將彌漫性膠質(zhì)瘤分為WHO II級和III級星形細胞瘤、II級和III級少突膠質(zhì)細胞瘤、IV級膠質(zhì)母細胞瘤以及兒童相關的彌漫性膠質(zhì)瘤[2]。不同級別的膠質(zhì)瘤在腫瘤病理形態(tài)(如膠原纖維含量和形態(tài)多樣性)、腫瘤發(fā)展和患者預后等方面差異很大[3]。其中膠質(zhì)母細胞瘤(glioblastoma, GBM, WHO IV)最為常見，且患者預后最差，總生存期中位數(shù)僅為12～ 14個月；但是，有3%～5%的患者可以存活超過3年，被稱為長期存活者[4]。與長期存活相關的臨床和分子因素仍然匱乏，目前在彌漫性膠質(zhì)瘤中僅有異檸檬酸脫氫酶(isocitrate dehydrogenases, IDH)突變狀態(tài)、O6-甲基鳥嘌呤-DNA 甲基轉移酶(O-6-methylguanine- DNA methyltransferase, MGMT)基因啟動子甲基化和染色體1p和19q聯(lián)合缺失被廣泛認知[4]。GBM患者不同的生存期差異也反映出關鍵腫瘤相關基因表達水平或基因組改變存在不同[5]。

拷貝數(shù)變異(copy number variations, CNV)通常被定義為涉及大規(guī)模(>1 kb)基因組DNA變化的結構變異，CNV是基因表達的顯著影響因素，可能影響多種致癌或抑癌關鍵通路的活性[6～9]。在膠質(zhì)瘤中常見的拷貝數(shù)變異包括7號染色體片段的擴增，9號和10號染色體片段的缺失，以及在小范圍內(nèi)1號和19號染色體的缺失[10]。研究顯示這些基因組結構變化主要與腫瘤相關基因的擴增，以及腫瘤抑制基因、和等的缺失或突變有關[3,11]。

各種變異對于GBM生存期的影響，目前GBM預后模型研究主要圍繞轉錄組水平開展，如構建自噬相關的4個基因(和)的GBM預后模型[12]，缺氧相關的5個基因(和)的GBM預后模型[13]，免疫狀態(tài)相關的15個基因的GBM預后模型[14]。另外通過篩選與GBM預后相關的4個miRNA (和)來建立預后模型[15]。在多組學層面，有研究通過對DNA甲基化和基因表達的綜合分析識別關鍵標志物(和)來建立預后模型[16]。由于拷貝數(shù)變異的復雜性，目前對于GBM拷貝數(shù)變異相關的預后標志物的系統(tǒng)研究尚不多見，因此探討不同生存期GBM之間的轉錄組水平差異，并結合GBM基因組層面的拷貝數(shù)變異數(shù)據(jù)，進行聯(lián)合分析，將有助于更好地識別GBM的關鍵預后分子和潛在的治療靶點。本研究對TCGA數(shù)據(jù)庫中GBM患者的基因組數(shù)據(jù)和表達譜數(shù)據(jù)聯(lián)合分析，建立了一個多基因預后模型，其在不同數(shù)據(jù)集中都能夠很好地識別預后不良的高風險患者，為GBM患者的預后風險評估提供更好的依據(jù)。

1 材料與方法

1.1 數(shù)據(jù)收集

運用TCGAbiolinks包[17]下載癌癥基因組圖譜計劃(The Cancer Genome Atlas, TCGA)數(shù)據(jù)庫(https:// portal.gdc.cancer.gov/)中GBM患者的RNA-seq表達數(shù)據(jù)(GBM. RNAseq HtSeq Count, level 3)，DNA拷貝數(shù)變異數(shù)據(jù)(GBM. SNP6.0 array)，同時下載與患者相關臨床數(shù)據(jù)。提取的臨床數(shù)據(jù)信息包括總生存期(overall survival time, OS.time)、年齡、性別、基因突變狀態(tài)等。我們用樣本匹配RNA-seq數(shù)據(jù)，CNV數(shù)據(jù)和對應的臨床數(shù)據(jù)提取了139個GBM患者數(shù)據(jù)用于分析[18]。

我們一共采用了3組數(shù)據(jù)作為驗證模型。從美國國家生物技術信息中心(National Center for Bio-technology Information，NCBI)的高通量芯片表達譜數(shù)據(jù)庫(Gene Expression Omnibus database, GEO, https://www.ncbi.nlm.nih.gov/geo/)中獲取膠質(zhì)瘤的GSE16011數(shù)據(jù)集，數(shù)據(jù)集中包含159個GBM樣本，從中獲取150個包含總生存期和生存狀態(tài)信息的GBM樣本作為驗證組進行后續(xù)分析和模型驗證。從中國腦膠質(zhì)瘤圖譜數(shù)據(jù)庫(Chinese Glioma Genome Atlas, CGGA, http://www.cgga.org.cn/)下載了Part C的mRNAseq_325表達譜數(shù)據(jù)和臨床數(shù)據(jù)，我們將GBM的88個樣本數(shù)據(jù)作為驗證組進行后續(xù)分析和模型驗證。此外，我們還從美國腦腫瘤分子數(shù)據(jù)庫(The Repository of Molecular Brain Neoplasia Data, Rembrandt, http://gliovis.bioinfo.cnio.es/)提取146個GBM樣本的表達譜數(shù)據(jù)和臨床數(shù)據(jù)作為驗證組進行后續(xù)分析和模型驗證。總生存期定義為從患者診斷日期開始至死亡日期截止。TCGA中GBM數(shù)據(jù)和GEO中GSE16011數(shù)據(jù)下載日期截止于2020年9月25日，CGGA數(shù)據(jù)和Rembrandt數(shù)據(jù)下載日期截止于2021年1月18日。

1.2 分析流程

本研究分析流程見圖1所示。

1.3 不同生存期GBM患者表達譜差異分析

對于TCGA中139個GBM表達譜數(shù)據(jù)，將每個患者的總生存期小于1年(OS.time<1 year)作為短生存期患者，大于3年(OS.time>3 year)作為長生存期患者進行分組。運用Deseq2包[19]進行差異分析，篩選<0.05，并且絕對值Log2(fold change)>1作為短生存期組與長生存期組相比差異的顯著基因。最終，通過火山圖展示短生存期患者和長生存期患者的差異基因。

圖1 分析流程圖

1.4 GBM拷貝數(shù)變異數(shù)據(jù)分析

腫瘤中重要靶點的基因組識別(genomic iden-tification of significant targets in cancer, GISTIC)是一種用于識別更有可能引發(fā)腫瘤發(fā)病機制的變異區(qū)域的算法[20]。利用GISTIC 2.0軟件對139個GBM患者CNV數(shù)據(jù)進行分析，選擇0.1作為拷貝數(shù)增加和缺失的閾值，-value=0.25作為顯著性閾值，確定了GBM患者中所有基因的五種離散拷貝數(shù)狀態(tài)(?2, ?1, 0, 1, 2)。以基因的(1, 2)拷貝數(shù)作為基因擴增組，(?2, ?1, 0)拷貝數(shù)作為基因非擴增組。

1.5 生存分析

Kaplan-Meier(KM)生存分析[21]用于分析基因與患者預后的關系。我們將GISTIC 2.0分析得到不同擴增狀態(tài)的基因進行KM生存分析，篩選<0.05作為影響生存的擴增基因(survival-associated amplifica-tion gene, SAG)。

1.6 單因素Cox回歸和迭代Lasso回歸篩選候選基因

將差異分析中短生存期組上調(diào)基因和SAG取交集基因；通過survival包對交集基因進行單因素Cox回歸分析，篩選<0.10的基因。為了進一步篩選重要的預后基因，避免單因素Cox回歸分析的過度擬合問題，通過glmnet包[22]對單因素Cox回歸分析結果中得到的基因進行500次迭代Lasso回歸分析。選擇被Lasso重復抽取大于400次的變量作為最終的關鍵預后基因[23]。

1.7 多因素Cox回歸分析與預后模型構建

迭代Lasso回歸分析篩選出的關鍵預后基因，基因迭代頻次的大小反應了基因的重要程度。根據(jù)迭代頻次大小排序，依次通過survival包納入單個基因進行多因素Cox比例回歸模型的構建，篩選接受者操作特性曲線(receiver operating characteristic curve, ROC曲線)下面積最大的基因集作為最佳預后模型。該步驟通過timeROC包[24]來完成。并用predict函數(shù)預測每個患者風險得分，將風險得分的中位數(shù)進行劃分高低風險組進行KM生存分析。并用GEO、CGGA和Rembrandt數(shù)據(jù)庫3個外部數(shù)據(jù)集來驗證模型的優(yōu)良。

1.8 富集分析

利用cluster profiler包[25]對短生存期組上調(diào)基因進行GO功能富集和KEGG通路富集分析，以確定在3個類別如生物學過程(biological process, BP)、細胞成分(cellular component, CC)和分子功能(mole-cular function, MF)中過度表達的功能富集以及過度表達的KEGG通路。對于該分析，F(xiàn)DR<0.05被認為具有統(tǒng)計學意義?；蚋患治?gene set enrich-ment analysis, GSEA)以不同生存期的差異分析后所得到的差異基因的Log2(fold change)值排序為輸入文件，使用分子特征數(shù)據(jù)集(molecular signatures database, MSigDB)中GO基因集(C5)，KEGG基因集(C2)和特征基因集(HALLMARK)注釋基因富集情況。GSEA通過fgsea包[26]完成，標準化富集得分(norma-lized enrichment score, NES)和校正值用來評價基因集富集效果，其中校正<0.05為有統(tǒng)計學意義。

1.9 加權基因共表達網(wǎng)絡構建

基于TCGA中GBM的RNA-seq表達譜數(shù)據(jù)，篩選出方差最大的前5000個基因，利用加權基因共表達網(wǎng)絡(weighted gene co-expression network ana-lysis, WGCNA)包[27]對這5000個基因的表達矩陣進行分析。WGCNA簡要步驟如下：①基因模塊識別：選擇合適的軟閾值β構建網(wǎng)絡，最佳軟閾值需達到無尺度擬合指數(shù)0.9以上，且需要網(wǎng)絡的平均連通性較高。以此軟閾值構建拓撲重疊矩陣(topological overlap matrix, TOM)。WGCNA將具有高拓撲重疊指數(shù)的一組基因定義為同一模塊，對基因進行聚類，并按照混合動態(tài)剪切樹的方法確定基因模塊。②模塊篩選：將模塊進行主成分分析，即各模塊第一主成分與臨床表型進行皮爾森相關分析，得到模塊與表型的相關系數(shù)即模塊顯著性(module significance, MS)及值，<0.05被認為有統(tǒng)計學意義。

我們運用WGCNA包構建了高低風險組的加權基因共表達網(wǎng)絡。WGCNA具體參數(shù)如下：當使用0.9作為相關系數(shù)閾值時，我們選擇的軟閾值能力是4，并且選擇10作為模塊中的最小基因數(shù)。為了合并可能的相似模塊，我們將0.2定義為切割高度的閾值。

1.10 統(tǒng)計分析

利用R 3.5.3軟件進行數(shù)據(jù)下載、統(tǒng)計分析和相應的數(shù)據(jù)可視化。包括TCGAbiolinks包下載TCGA中GBM的RNA-seq表達譜數(shù)據(jù)、CNV數(shù)據(jù)和臨床數(shù)據(jù)。利用Deseq2包對不同生存期進行差異分析。通過survival包進行KM生存分析、單變量、多變量Cox比例風險分析和模型構建，并利用survminer包進行相應的可視化。通過glmnet包進行迭代Lasso回歸分析。利用ggplot2包和pheatmap包繪制風險因子關聯(lián)圖。

2 結果與分析

2.1 GBM不同生存期患者表達譜差異分析

用Deseq2包對TCGA中139例GBM的短生存期組(OS.time<1 year) 79例和長生存期組(OS.time> 3 year) 7例患者的表達譜數(shù)據(jù)進行差異分析，根據(jù)篩選標準<0.05且絕對值Log2(fold change)>1，我們一共篩選得到491個mRNA基因(附表1)，其中上調(diào)基因426個，下調(diào)基因65個(圖2)。為了探索短生存期組上調(diào)基因參與的生物學過程，進行了GO和KEGG富集分析。根據(jù)GO分析結果顯示上調(diào)基因主要分布于細胞外基質(zhì)、含膠原蛋白的細胞外基質(zhì)、內(nèi)質(zhì)網(wǎng)內(nèi)腔等組織；參與了細胞外組織、白細胞遷移、膠原組織等生物過程；主要與受體結合、生長因子結合等生物過程有關(圖3A)。與上調(diào)基因相關的KEGG通路，包括細胞因子–受體相互作用，PI3K-Akt信號通路，腫瘤中的蛋白多糖，IL-17信號通路，腫瘤壞死因子信號通路等(圖3B)。

為了彌補由于表達量閾值的影響，可能會遺漏某些差異的富集通路。我們對不同生存期差異分析的差異基因進行GSEA分析，按照校正后值從大到小排列，在GO基因集中，短生存期組顯著富集了與細胞外基質(zhì)、含膠原的細胞外基質(zhì)和細胞外結構等生物過程(附圖1A)。在KEGG基因集中，短生存期組顯著富集了趨化因子信號通路、細胞因子受體相互作用等，與單獨選取上調(diào)基因富集分析結果部分重疊(附圖1B)。在HALLMARK基因集中，短生存期組顯著富集了多種腫瘤相關的信號通路。如上皮間質(zhì)型轉化(epithelial-mesenchymal transition, EMT)、TNFA，KRSA信號通路、缺氧等(附圖1C)。通過對不同生存期組間差異表達基因的GSEA分析，進一步確認了上述在短生存期組上調(diào)基因富集分析相關的信號通路和生物過程。

圖2 TCGA中GBM患者不同生存期基因差異表達分析

差異基因火山圖分析，紅點代表上調(diào)的差異基因，藍點代表下調(diào)的差異基因，灰點代表沒有差異的基因。Log2(fold change)為差異分析中基因的差異倍數(shù)，OS.time<1 year比OS.time>3 year。

圖3 TCGA中GBM短生存期組上調(diào)基因GO和KEGG富集分析

A：GO功能富集分析結果。BP：生物學過程；CC：細胞組分；MF：分子功能。B：KEGG通路富集分析結果。

2.2 GISTIC和KM分析CNV數(shù)據(jù)

對139例GBM的CNV數(shù)據(jù)進行GISTIC分析，得到所有基因離散化的拷貝數(shù)狀態(tài)。將所有基因分成擴增組和非擴增組進行KM生存分析，根據(jù)<0.05得到2405個SAG。將短生存期組上調(diào)基因和SAG交集得到的47個候選基因，用于后續(xù)分析(附表2)。

2.3 構建GBM預后風險模型

為了進一步篩選預后基因，我們將得到的47個候選基因，進行單因素Cox比例風險分析，篩選<0.10的基因。篩選得到18個基因見表1所示。對上述基因集進行迭代Lasso回歸，進行500次迭代處理，篩選頻次大于400的基因。獲得12個基因依據(jù)頻次大小分別納入多因素Cox回歸模型，計算每次納入基因的ROC曲線下的面積。選取面積最大的基因集合作為TCGA實驗組的最佳模型。最終，我們確定了六基因(、、、、、)組合，六基因組合在多因素Cox回歸模型中ROC曲線下的面積最大為0.912 (圖4，A和B)。用predict函數(shù)預測TCGA實驗組每個患者風險得分，根據(jù)風險得分的中位值(1.116)，我們將139例中69例納入高風險組，70例納入低風險組。從風險因子關聯(lián)圖中可以看到，高風險組的死亡人數(shù)明顯更多，存活人數(shù)更少(圖4C)。6個基因表達量在患者中的分布表明，高風險組患者的6個基因的表達量比低風險組的更高。并且六基因模型的風險評分分組的KM分析顯示，高低風險組存在明顯差異。高風險組預后較差，<0.001具有統(tǒng)計學意義(圖4D)。

表1 單因素Cox回歸篩選P<0.10的基因

HR為風險比(hazard ratio)用于描述相對危險度的指標。

Log2(fold change)為差異分析中基因的差異倍數(shù)，OS.time<1 year比OS.time>3 year。

圖4 TCGA中GBM實驗組六基因預后模型建立

A：根據(jù)迭代Lasso回歸的迭代次數(shù)納入回歸模型，選取曲線下面積的基因集合作為最佳模型。AUC (area under curve)被定義為ROC曲線下與坐標軸圍成的面積。B：六基因預后模型的ROC曲線下面積為0.912。C：六基因預后模型風險因子關聯(lián)圖。上圖：預后模型中高風險(紅色)和低風險(藍色) GBM患者的風險評分分布；中圖：散點圖顯示了模型中GBM患者的生存狀況。紅點表示死亡的患者，藍點表示存活的患者；下圖：模型中6個基因的表達量熱圖，橫坐標為基因表達量，縱坐標為GBM樣本。D：高風險組和低風險組的KM生存分析。

2.4 外部數(shù)據(jù)驗證六基因預后模型

為了對六基因預后模型進行驗證，我們首先將六基因模型應用于GEO腦膠質(zhì)瘤的GSE16011表達譜驗證組，其ROC曲線下面積為0.825 (圖5A)。說明六基因預后模型效能好，對于GBM的預后具有很好的預測能力。用predict函數(shù)預測GSE16011驗證組每個患者風險得分，根據(jù)風險得分的中位值(1.0347)，我們將150例中74例納入高風險組，76例納入低風險組。GSE16011驗證組的風險因子關聯(lián)圖同TCGA實驗組中具有一致性，高風險組的死亡人數(shù)多，存活人數(shù)少(圖5B)。6個基因表達量隨著風險評分增高而增高。KM生存分析同樣顯示了高風險組預后較差(圖5C)。同樣我們將六基因模型應用于CGGA驗證組，結果顯示，CGGA驗證組的ROC曲線下面積為0.636 (附圖2A)，說明六基因預后模型在CGGA驗證組有中等效能，對于GBM的預后具有一定的預測能力。CGGA驗證組的風險因子關聯(lián)圖同TCGA實驗組中具有一致性，高風險組的死亡人數(shù)多，存活人數(shù)少(附圖2B)。用predict函數(shù)預測CGGA驗證組每個患者風險得分，根據(jù)風險得分的中位值(0.998)我們將88例中44例納入高風險組，44例納入低風險組。CGGA驗證組(附圖2C) KM生存分析顯示了高風險組預后較差的結果。進而我們將六基因模型應用于Rembrandt驗證組，結果顯示，其ROC曲線下面積為0.845 (附圖3A)，說明六基因預后模型效能好，對于GBM的預后具有很好的預測能力。Rembrandt驗證組的風險因子關聯(lián)圖同TCGA實驗組中具有一致性，高風險組的死亡人數(shù)多，存活人數(shù)少(附圖3B)。用predict函數(shù)預測Rembrandt驗證組每個患者風險得分，根據(jù)風險得分的中位值(0.986)，我們將146例中73例納入高風險組，73例納入低風險組。Rembrandt驗證組(附圖3C) KM生存分析同樣顯示了高風險組預后較差的結果。值均小于0.05，具有統(tǒng)計學意義。從3個驗證組結果看，六基因預后模型對實驗組和外部數(shù)據(jù)驗證組都有較好的預測作用，該模型可以對GBM患者的生存狀態(tài)進行預測。

圖5 GEO中GBM GSE16011驗證組的六基因預后模型

A：六基因預后模型在GSE16011驗證組的ROC曲線下面積為0.825。B：GSE16011驗證組的六基因預后模型風險因子關聯(lián)圖。上圖：預后模型中高風險(紅色)和低風險(藍色)GBM患者的風險評分分布；中圖：散點圖顯示了模型中GBM患者的生存狀況。紅點表示死亡的患者，藍點表示存活的患者；下圖：模型中6個基因的表達量熱圖，橫坐標為基因表達量，縱坐標為GBM樣本。C：GSE16011驗證組的高風險組和低風險組的KM生存分析。

2.5 評估六基因模型預后獨立性

為了評估六基因模型的預后獨立性，我們結合患者的臨床指標包括：六基因預后評分分組、年齡、性別和基因突變狀態(tài)等進行多因素Cox回歸比例風險模型。結果見表2所示。在TCGA GBM實驗組中，多因素Cox回歸分析顯示，六基因評分風險分組(HR=2.39, 95%CI=1.582～3.617,=0.0004)和年齡(HR=1.02, 95%CI=1.005～1.04,=0.01056)與生存顯著相關，并且均為危險因素，HR均大于1。在GEO GSE16011驗證組中，多因素Cox回歸分析顯示，六基因評分風險分組(HR=1.46, 95%CI=1.012～ 2.125,=0.04285)和年齡(HR=1.03, 95%CI=1.02～ 1.052,=0.0001)與生存顯著相關，同樣均為危險因素，HR大于1。在CGGA驗證組中，多因素Cox回歸分析顯示，六基因評分風險分組(HR=1.92, 95%CI=1.189～3.102,=0.00766)與生存顯著相關，為危險因素，HR大于1。在Rembrandt驗證組中，多因素Cox回歸分析顯示，六基因評分風險分組(HR=1.76, 95%CI=1.18～2.61,=0.005)與生存顯著相關，為危險因素，HR大于1。綜上所述，六基因預后模型可以作為一個獨立于其他臨床因素的預后指標，用于GBM的風險預后。

表2 多因素Cox回歸比例風險模型

2.6 加權共表達網(wǎng)絡分析

通過WGCNA包將RNA-seq表達譜數(shù)據(jù)中方差最大的前5000基因用于加權共表達網(wǎng)絡分析。當使用0.9作為相關系數(shù)閾值時，軟閾值被選為4(圖6A)。通過WGCNA分析，構建了13個共表達模塊(圖6B)，與性狀相關分析顯示黑色模塊與高風險組最相關，模塊顯著性(MS)為0.33，<0.05 (圖6C)。我們對黑色模塊包含的206個基因進行GO功能和KEGG富集分析。GO功能富集分析結果顯示，模塊內(nèi)的基因主要分布在細胞外基質(zhì)，含膠原蛋白的細胞外基質(zhì)，內(nèi)質(zhì)網(wǎng)內(nèi)腔等組織，參與了細胞外組織、白細胞遷移、與氧反應等生物過程；主要有與受體相關，生長因子相關，胞外基質(zhì)結合等分子功能(圖7A)。與模塊內(nèi)的基因相關的KEGG通路，包括細胞因子–細胞因子受體相互作用、腫瘤壞死因子信號通路、焦點粘著、類風濕性關節(jié)炎、IL?17信號通路等(圖7B)。與前面的上調(diào)基因富集分析在GO功能和KEGG通路上具有高度一致性。

3 討論

GBM作為最常見、惡性程度最高的膠質(zhì)瘤亞型，雖然總體預后差，但是具有廣泛分子異質(zhì)性和患者的生存期差異大的特點。目前研究表明可以作為預測GBM預后的相關分子生物標志物，包括基因突變、染色體1p/19q聯(lián)合缺失和基因啟動子甲基化等。具有上述相關分子改變的GBM患者對化療敏感并且具有預后較好的特點。然而具有擴增、重排、、突變和融合或點突變的分子改變患者提示預后不佳的效果。實際上，目前已有研究通過篩選一系列與GBM相關的功能核心基因來建立預后模型，如自噬、缺氧、免疫等方面相關分子。上述特征是基于單一轉錄組學數(shù)據(jù)的生物學過程發(fā)展起來的。一個層面數(shù)據(jù)的改變很難解釋疾病的整個發(fā)生過程。因此，本研究對TCGA數(shù)據(jù)庫中GBM患者的RNA-seq表達譜數(shù)據(jù)和CNV數(shù)據(jù)聯(lián)合分析，建立六基因(、、、、、)預后模型，發(fā)現(xiàn)TCGA GBM實驗組中該風險模型是很好的預后預測因子，能夠很好地識別預后不良高風險的GBM患者；并且在GEO中GSE16011驗證組、CGGA驗證組和Rembrandt驗證組顯示該模型具有良好的穩(wěn)定性和可重復性。相較與以往研究，多組學聯(lián)合分析更有利于全面評估患者的預后風險，同時也為后續(xù)的分析研究提供了思路?；诨蚪M和轉錄組構建的預測模型能夠很好地識別預后不良高風險的GBM患者，并且能夠在不同數(shù)據(jù)集中得到很好的驗證，結果具有一致性，故能為GBM患者的預后風險評估提供更好的依據(jù)。

圖6 WGCNA加權共表達網(wǎng)絡識別高風險組相關模塊

A：基于TCGA中GBM RNA-seq表達數(shù)據(jù)的軟閾值的確定。B：基因聚類樹狀圖。樹形圖下方的色塊表示由動態(tài)樹切割方法識別的模塊，模塊內(nèi)的基因高度相關。C：模塊與高低風險組相關性的熱圖。色塊中上面的數(shù)字代表相關性，下面數(shù)字代表值，紅色呈正相關，綠色呈負相關。

圖7 高風險組相關模塊內(nèi)的基因GO和KEGG富集分析

A：高風險組相關模塊內(nèi)的基因GO功能富集分析結果。BP：生物學過程；CC：細胞組分；MF：分子功能。B：高風險組相關模塊內(nèi)的基因通路富集分析結果。

六基因模型區(qū)分的高風險組相關基因的富集分析結果，與上調(diào)基因的富集分析結果核心成分相重疊，主要集中在細胞外基質(zhì)、白細胞遷移、細胞外組織、受體結合和生長因子相關等功能；細胞因子受體相互作用，腫瘤壞死因子信號通路，IL?17信號通路等。并且在GSEA分析HALLMARK基因集中，短生存期組顯著富集了腫瘤EMT信號通路，缺氧等其他經(jīng)典腫瘤相關信號通路。細胞外基質(zhì)相關生物過程通常與腫瘤發(fā)生EMT甚至侵襲遷移的惡性表型密切相關[28]。因此，六基因模型可以獲得與GBM預后相關的核心細胞生物學過程，細胞組分，分子功能及相關通路。進一步說明了六基因組合模型可以較好地用于預測GBM患者的預后。

在膠質(zhì)瘤以往的研究報道中，6個基因中有4個基因與膠質(zhì)瘤密切相關?；蚓幋a一種含有同源結構域的轉錄因子，在胚胎神經(jīng)發(fā)育過程中是必不可少的。LEF1-AS1通過上調(diào)EN2參與GBM的惡性發(fā)展[29]?；蚓幋a的蛋白質(zhì)是一種血小板衍生生長因子，屬于CXC趨化因子家族，這種生長因子是一種有效的中性粒細胞趨化劑和激活劑。它被證明可以刺激各種細胞過程，包括DNA合成、有絲分裂、糖酵解等，在發(fā)育和治療耐藥中發(fā)揮作用。PPBP趨化因子可能是未來GBM研究的一個有希望的靶點[30]，該基因編碼的蛋白也被預測作為肺癌早期診斷的標志物[31]。基因編碼富含亮氨酸重復序列蛋白61。在人類蛋白質(zhì)圖譜(Human Protein Atlas, HPA, https://www.proteinatlas. org/)數(shù)據(jù)庫中，免疫組化結果顯示LRRC61在高級別膠質(zhì)瘤組織和低級別膠質(zhì)瘤組織中表達高于正常大腦皮層組織(附圖4，A～C)。基因編碼SEL1L家族成員3，在HPA數(shù)據(jù)庫中，免疫組化結果顯示SEL1L3也在高級別膠質(zhì)瘤組織和低級別膠質(zhì)瘤組織中表達高于正常大腦皮層組織(附圖4，D～F)。另外兩種基因中，基因是羧肽酶A/B亞家族的一員，與第7號染色體上的另外3個家族成員位于一個簇中。羧肽酶是一類含鋅的外肽酶，催化羧基末端氨基酸的釋放，并被合成為酶原，通過蛋白水解被激活。該基因可能參與組蛋白高乙酰化途徑[32]，有文獻報道，通過CPA4激活STAT3和ERK1/2信號通路促進細胞生長并預測結直腸癌預后不良[33]?；蚓幋aDNA損傷誘導轉錄樣蛋白4，有研究表明DDIT4L可能是通過心臟mTOR信號傳導病理應激到自噬的重要傳感器，并且DDIT4L可以作為心血管疾病的治療靶點，在自噬和mTOR信號通路起主要作用[34]。

綜上所述，根據(jù)我們的研究，提出六基因(、、、、、)預后模型可為GBM患者的預后風險評估提供更好的依據(jù)。其包含已知與GBM相關的基因和未知基因，值得我們在今后的臨床樣本中驗證，并進行更多的臨床學與功能學研究。

附錄：

附加材料詳見文章電子版www.chinagene.cn。

致謝：

感謝哈爾濱醫(yī)科大學傅松濱教授在指導文章寫作和審閱文章過程中給予的幫助。

[1] Schwartzbaum JA, Fisher JL, Aldape KD, Wrensch M. Epidemiology and molecular pathology of glioma., 2006, 2(9): 494–503.

[2] Su CL, Li L, Chen XW, Zhang J, Shen NX, Wang ZX, Yang SQ, Li J, Zhu WZ, Wang CY. Summary of classification of central nervous system tumors in WHO 2016., 2016, 31(7): 570–579.蘇昌亮, 李麗, 陳小偉, 張巨, 申楠茜, 王振熊, 楊時騏, 李娟, 朱文珍, 王承緣. 2016年WHO中樞神經(jīng)系統(tǒng)腫瘤分類總結. 放射學實踐, 2016, 31(7): 570–579.

[3] Parsons DW, Jones S, Zhang XS, Lin JCH, Leary RJ, Angenendt P, Mankoo P, Carter H, Siu IM, Gallia GL, Olivi A, McLendon R, Rasheed BA, Keir S, Nikolskaya T, Nikolsky Y, Busam DA, Tekleab H, Diaz LA, Hartigan J, Smith DR, Strausberg RL, Marie SKN, Shinjo SMO, Yan H, Riggins GJ, Bigner DD, Karchin R, Papadopoulos N, Parmigiani G, Vogelstein B, Velculescu VE, Kinzler KW. An integrated genomic analysis of human glioblastoma multiforme., 2008, 321(5897): 1807–1812.

[4] Krex D, Klink B, Hartmann C, von Deimling A, Pietsch T, Simon M, Sabel M, Steinbach JP, Heese O, Reifenberger G, Weller M, Schackert G, German Glioma Network. Long-term survival with glioblastoma multiforme., 2007, 130: 2596–2606.

[5] Gerber NK, Goenka A, Turcan S, Reyngold M, Makarov V, Kannan K, Beal K, Omuro A, Yamada Y, Gutin P, Brennan CW, Huse JT, Chan TA. Transcriptional diversity of long-term glioblastoma survivors., 2014, 16(9): 1186–1195.

[6] Asiedu MK, Thomas CF, Dong J, Schulte SC, Khadka P, Sun ZF, Kosari F, Jen J, Molina J, Vasmatzis G, Kuang R, Aubry MC, Yang P, Wigle DA. Pathways impacted by genomic alterations in pulmonary carcinoid tumors., 2018, 24(7): 1691–1704.

[7] Chang J, Tan WL, Ling ZQ, Xi RB, Shao MM, Chen MJ, Luo YY, Zhao YJ, Liu Y, Huang XC, Xia YC, Hu JL, Parker JS, Marron D, Cui QH, Peng LN, Chu JH, Li HM, Du ZL, Han YL, Tan W, Liu ZH, Zhan QM, Li Y, Mao WM, Wu C, Lin DX. Genomic analysis of oesophageal squamous-cell carcinoma identifies alcohol drinking- related mutation signature and genomic alterations., 2017, 8: 15290.

[8] Liang L, Fang JY, Xu J. Gastric cancer and gene copy number variation: emerging cancer drivers for targeted therapy., 2016, 35(12): 1475–1482.

[9] Wang Z, Wang LX. Gene amplification in lung cancer., 2010, 33(3): 155–159.王箴, 王良旭. 基因擴增在肺癌中的研究進展. 國際遺傳學雜志, 2010, 33(3): 155–159.

[10] Ohgaki H, Kleihues P. Genetic alterations and signaling pathways in the evolution of gliomas., 2009, 100(12): 2235–2241.

[11] Jin WX, Yang TM, Dong XY, Hua ZC, Xu XX. P53 mutation in human brain gliomas: Comparison of loss of heterozygosity with mutation frequency., 2000, 22(6): 357–360.金衛(wèi)新, 楊天明, 董雪吟, 華子春, 徐賢秀. 腦膠質(zhì)瘤中P53基因突變及其與染色體17p雜合性丟失的比較. 遺傳, 2000, 22(6): 357–360.

[12] Wang YL, Zhao WJ, Xiao Z, Guan GF, Liu X, Zhuang MH. A risk signature with four autophagy-related genes for predicting survival of glioblastoma multiforme., 2020, 24(7): 3807–3821.

[13] Lin WZ, Wu SH, Chen XC, Ye YL, Weng YL, Pan YH, Chen ZJ, Chen L, Qiu XX, Qiu SF. Characterization of hypoxia signature to evaluate the tumor immune microenvironment and predict prognosis in glioma groups., 2020, 10: 796.

[14] Wang JJ, Wang H, Zhu BL, Wang X, Qian YH, Xie L, Wang WJ, Zhu J, Chen XY, Wang JM, Ding ZL. Development of a prognostic model of glioma based on immune-related genes., 2021, 21(2): 116.

[15] Hermansen SK, S?rensen MD, Hansen A, Knudsen S, Alvarado AG, Lathia JD, Kristensen BW. A 4-miRNA signature to predict survival in glioblastomas., 2017, 12(11): e0188090.

[16] Jia DY, Lin W, Tang HL, Cheng YF, Xu KW, He YS, Geng WJ, Dai QX. Integrative analysis of DNA methylation and gene expression to identify key epigenetic genes in glioblastoma., 2019, 11(15): 5579– 5592.

[17] Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, Sabedot TS, Malta TM, Pagnotta SM, Castiglioni I, Ceccarelli M, Bontempi G, Noushmehr H. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data., 2016, 44(8): e71.

[18] Li X, Li MW, Zhang YN, Xu HM. Common cancer genetic analysis methods and application study based on TCGA database., 2019, 41(3): 234–242.李鑫, 李夢瑋, 張依楠, 徐寒梅. 常用腫瘤基因分析方法及基于TCGA數(shù)據(jù)庫的分析應用. 遺傳, 2019, 41(3): 234–242.

[19] Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2., 2014, 15(12): 550.

[20] Beroukhim R, Getz G, Nghiemphu L, Barretina J, Hsueh T, Linhart D, Vivanco I, Lee JC, Huang JH, Alexander S, Du JY, Kau T, Thomas RK, Shah K, Soto H, Perner S, Prensner J, Debiasi RM, Demichelis F, Hatton C, Rubin MA, Garraway LA, Nelson SF, Liau L, Mischel PS, Cloughesy TF, Meyerson M, Golub TA, Lander ES, Mellinghoff I, Sellers WR. Assessing the significance of chromosomal aberrations in cancer: methodology and application to glioma., 2007, 104(50): 20007–20012.

[21] Stalpers LJA, Kaplan EL. Edward L. Kaplan and the Kaplan-Meier survival curve., 2018, 33(2): 109–135.

[22] Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent., 2010, 33(1): 1–22.

[23] Sveen A, ?gesen TH, Nesbakken A, Meling GI, Rognum TO, Liest?l K, Skotheim RI, Lothe RA. ColoGuidePro: a prognostic 7-gene expression signature for stage III colorectal cancer patients., 2012, 18(21): 6001–6010.

[24] Blanche P, Dartigues JF, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks., 2013, 32(30): 5381–5397.

[25] Yu GC, Wang LG, Han YY, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters., 2012, 16(5): 284–287.

[26] Korotkevich G, Sukhov V, Budin N, Shpak B, Artyomov M, Sergushichev A. Fast gene set enrichment analysis., 2019: 060012.

[27] Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis., 2008, 9: 559.

[28] Zeng Y, Xiao HL, Guo QN. Research progress on the mechanism of epithelial mesenchymal transition in glioblastoma., 2019, 26(8): 533–538.曾英, 肖華亮, 郭喬楠. 膠質(zhì)母細胞瘤上皮間葉轉化機制研究進展. 診斷病理學雜志, 2019, 26(8): 533–538.

[29] Zeng S, Zhou C, Yang DH, Xu LS, Yang HJ, Xu MH, Wang H. LEF1-AS1 is implicated in the malignant development of glioblastoma via sponging miR-543 to upregulate EN2., 2020, 1736: 146781.

[30] Achyut BR, Shankar A, Iskander ASM, Ara R, Angara K, Zeng P, Knight RA, Scicli AG, Arbab AS. Bone marrow derived myeloid cells orchestrate antiangiogenic resistance in glioblastoma through coordinated molecular networks., 2015, 369(2): 416–426.

[31] Du Q, Li EC, Liu YE, Xie WL, Huang C, Song JQ, Zhang W, Zheng YJ, Wang HL, Wang Q. CTAPIII/CXCL7: a novel biomarker for early diagnosis of lung cancer., 2018, 7(2): 325–335.

[32] Ross PL, Cheng I, Liu X, Cicek MS, Carroll PR, Casey G, Witte JS. Carboxypeptidase 4 gene variants and early-onset intermediate-to-high risk prostate cancer., 2009, 9: 69.

[33] Pan HD, Pan JX, Ji L, Song SB, Lv H, Yang ZR, Guo YB. Carboxypeptidase A4 promotes cell growth via activating STAT3 and ERK signaling pathways and predicts a poor prognosis in colorectal cancer., 2019, 138: 125–134.

[34] Simonson B, Subramanya V, Chan MC, Zhang AF, Franchino H, Ottaviano F, Mishra MK, Knight AC, Hunt D, Ghiran I, Khurana TS, Kontaridis MI, Rosenzweig A, Das S. DDiT4L promotes autophagy and inhibits pathological cardiac hypertrophy in response to stress., 2017, 10(468): eaaf5967.

附錄

附圖1 TCGA中GBM不同生存期差異基因GSEA分析結果

Supplementary Fig. 1 GSEA analysis genes of GBM patients with different survival time in TCGA

A：GSEA分析GO基因集中，短生存期組顯著富集的結果。B：GSEA分析KEGG基因集中，短生存期組和長生存期顯著富集的結果。C：GSEA分析HALLMARK基因集中，短生存期組和長生存期顯著富集的結果。NES為normalized enrichment score，標準化富集得分。

附圖2 CGGA驗證組的六基因預后模型

Supplementary Fig. 2 Six-gene prognostic model of CGGA validation group

A：六基因預后模型在CGGA驗證組的ROC曲線下面積為0.636。B：CGGA驗證組的六基因預后模型風險因子關聯(lián)圖。上圖為預后模型中高風險(紅色)和低風險(藍色) GBM患者的風險評分分布；中圖：散點圖顯示了模型中GBM患者的生存狀況。紅點表示死亡的患者，藍點表示存活的患者；下圖為模型中6個基因的表達量熱圖，橫坐標為基因表達量，縱坐標為GBM樣本。C：CGGA驗證組的高風險組和低風險組的KM生存分析。

附圖3 Rembrandt驗證組的六基因預后模型

Supplementary Fig. 3 Six-gene prognostic model of Rembrandt validation group

A：六基因預后模型在Rembrandt驗證組的ROC曲線下面積為0.845。B：Rembrandt驗證組的六基因預后模型風險因子關聯(lián)圖。上圖為預后模型中高風險(紅色)和低風險(藍色) GBM患者的風險評分分布；中圖：散點圖顯示了模型中GBM患者的生存狀況。紅點表示死亡的患者，藍點表示存活的患者；下圖為模型中6個基因的表達量熱圖，橫坐標為基因表達量，縱坐標為GBM樣本。C：Rembrandt驗證組的高風險組和低風險組的KM生存分析。

附圖4 預后模型的基因在膠質(zhì)瘤組織和正常組大腦皮層組織中的蛋白表達情況

Supplementary Fig. 4 Protein expression of prognostic model genes in glioma tissue

A：LRRC61在正常大腦皮層組織中的免疫組化結果；B：LRRC61在低級別膠質(zhì)瘤組織中的免疫組化結果；C：LRRC61在高級別膠質(zhì)瘤組織中的免疫組化結果；D：SEL1L3在在正常大腦皮層組織中的免疫組化結果；E：SEL1L3在低級別膠質(zhì)瘤組織組織中的免疫組化結果；F：SEL1L3在高級別膠質(zhì)瘤組織中的免疫組化結果，圖片數(shù)據(jù)來源于人類蛋白質(zhì)圖譜數(shù)據(jù)庫(HPA)。

附表1 TCGA中GBM患者不同生存期的差異基因

續(xù)附表1

續(xù)附表1

續(xù)附表1

續(xù)附表1

續(xù)附表1

續(xù)附表1

Log2(fold change)為差異分析中基因的差異倍數(shù)，OS.time<1 year比OS.time>3 year。

Diff type為差異基因類型，UP為上調(diào)基因，DOWN為下調(diào)基因。

附表2 短生存期組上調(diào)基因和SAG交集結果

Supplementary Table 2 The intersection of up-regulated genes in short survival group and SAG

KM_為KM生存分析的值。為不同生存期組患者差異分析的值。

Establish six-gene prognostic model for glioblastoma based on multi-omics data of TCGA database

Changgui Lei, Xueyuan Jia, Wenjing Sun

Glioblastoma (GBM) is the most common primary intracranial tumor with extremely high malignancy and poor prognosis.In order to identify the GBM prognostic biomarkers and establish a prognostic model, we analyzed the expression profile data of GBM in The Cancer Genome Atlas (TCGA) database as the experimental group. First, we identified the differentially expressed genes of different survival periods among the GBM patients. The GISTIC software and Kaplan Meier (KM) survival curve were used to analyze the copy number variation of GBM to identify the survival-associated amplified gene (SAG). We selected the intersection genes of up-regulated ones in short survival group and SAG, performed univariate Cox regression and iterative Lasso regression with them to identify the important candidate genes and establish a prognostic model. Based on the model, the prognostic score was calculated. The patients were divided into high-risk and low-risk groups according to the median prognostic score. Meanwhile ROC curve was used to evaluate the validity of the model, applying the KM survival analysis of the high-risk and low-risk groups. Multivariate Cox regression analysis was used to determine the independence of the prognostic score. All the data were verified with three external datasets: GEO GSE16011, CGGA, and Rembrandt. The results showed that differential expression analysis of different survival periods of GBM identified 426 up-regulated genes and 65 down-regulated genes in the TCGA GBM dataset. The intersection of up-regulated genes in short survival group and SAG yielded 47 genes. After the screening, the six-gene combination (,,,,,) prognostic model was finally determined. The area under ROC curve of the model in TCGA experimental group and three external validation group were all greater than 0.6, even reaching 0.912. KM analysis showed that the prognosis of the high-risk and low-risk groups was significant different (<0.05). In the multivariate Cox regression analysis, the six-gene prognostic score was an independent factor influencing the prognosis of GBM patients (<0.05).In summary, this study established a prognostic model of six-gene (,,,,,) for GBM. This six-gene model has good predictive ability and could be used as an independent prognostic marker for GBM patients.

glioblastoma; multi-omics data; six-gene combination; prognostic model; The Cancer Genome Atlas

2020-12-11;

2021-03-20

教育部“創(chuàng)新團隊發(fā)展計劃”項目(編號IRT1230)資助[Supported by Program for Changjiang Scholars and Innovative Research Team in University (No. IRT1230)]

雷常貴，在讀碩士研究生，專業(yè)方向：醫(yī)學遺傳學。E-mail: leichanggui@hrbmu.edu.cn

賈學淵，博士，副教授，研究方向：醫(yī)學遺傳學。E-mail: jiaxueyuan@hrbmu.edu.cn

雷常貴和賈學淵并列第一作者。

孫文靖，博士，教授，研究方向：醫(yī)學遺傳學。E-mail: sunwj@ems.hrbmu.edu.cn

10.16288/j.yczz.20-428

2021/4/7 16:51:07

URI: https://kns.cnki.net/kcms/detail/11.1913.R.20210407.1059.002.html

(責任編委: 何順民)