雷常貴，賈學(xué)淵，孫文靖

研究報告

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

雷常貴，賈學(xué)淵，孫文靖

哈爾濱醫(yī)科大學(xué)醫(yī)學(xué)遺傳學(xué)研究室，中國遺傳資源保護(hù)與疾病防控教育部重點(diǎn)實(shí)驗(yàn)室(哈爾濱醫(yī)科大學(xué))，哈爾濱 150081

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

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

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

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

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

1 材料與方法

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

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

我們一共采用了3組數(shù)據(jù)作為驗(yàn)證模型。從美國國家生物技術(shù)信息中心(National Center for Bio-technology Information，NCBI)的高通量芯片表達(dá)譜數(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樣本作為驗(yàn)證組進(jìn)行后續(xù)分析和模型驗(yàn)證。從中國腦膠質(zhì)瘤圖譜數(shù)據(jù)庫(Chinese Glioma Genome Atlas, CGGA, http://www.cgga.org.cn/)下載了Part C的mRNAseq_325表達(dá)譜數(shù)據(jù)和臨床數(shù)據(jù)，我們將GBM的88個樣本數(shù)據(jù)作為驗(yàn)證組進(jìn)行后續(xù)分析和模型驗(yàn)證。此外，我們還從美國腦腫瘤分子數(shù)據(jù)庫(The Repository of Molecular Brain Neoplasia Data, Rembrandt, http://gliovis.bioinfo.cnio.es/)提取146個GBM樣本的表達(dá)譜數(shù)據(jù)和臨床數(shù)據(jù)作為驗(yàn)證組進(jìn)行后續(xù)分析和模型驗(yàn)證?？偵嫫诙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患者表達(dá)譜差異分析

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

圖1 分析流程圖

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

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

1.5 生存分析

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

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

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

1.7 多因素Cox回歸分析與預(yù)后模型構(gòu)建

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

1.8 富集分析

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

1.9 加權(quán)基因共表達(dá)網(wǎng)絡(luò)構(gòu)建

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

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

1.10 統(tǒng)計(jì)分析

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

2 結(jié)果與分析

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

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

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

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

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

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

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

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

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

2.3 構(gòu)建GBM預(yù)后風(fēng)險模型

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

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

HR為風(fēng)險比(hazard ratio)用于描述相對危險度的指標(biāo)。

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

圖4 TCGA中GBM實(shí)驗(yàn)組六基因預(yù)后模型建立

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

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

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

圖5 GEO中GBM GSE16011驗(yàn)證組的六基因預(yù)后模型

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

2.5 評估六基因模型預(yù)后獨(dú)立性

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

表2 多因素Cox回歸比例風(fēng)險模型

2.6 加權(quán)共表達(dá)網(wǎng)絡(luò)分析

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

3 討論

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

圖6 WGCNA加權(quán)共表達(dá)網(wǎng)絡(luò)識別高風(fēng)險組相關(guān)模塊

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

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

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

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

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

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

附錄：

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

致謝：

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

[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)腫瘤分類總結(jié). 放射學(xué)實(shí)踐, 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.王箴, 王良旭. 基因擴(kuò)增在肺癌中的研究進(jìn)展. 國際遺傳學(xué)雜志, 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ù)庫的分析應(yīng)用. 遺傳, 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ì)母細(xì)胞瘤上皮間葉轉(zhuǎn)化機(jī)制研究進(jìn)展. 診斷病理學(xué)雜志, 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分析結(jié)果

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

A：GSEA分析GO基因集中，短生存期組顯著富集的結(jié)果。B：GSEA分析KEGG基因集中，短生存期組和長生存期顯著富集的結(jié)果。C：GSEA分析HALLMARK基因集中，短生存期組和長生存期顯著富集的結(jié)果。NES為normalized enrichment score，標(biāo)準(zhǔn)化富集得分。

附圖2 CGGA驗(yàn)證組的六基因預(yù)后模型

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

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

附圖3 Rembrandt驗(yàn)證組的六基因預(yù)后模型

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

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

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

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

A：LRRC61在正常大腦皮層組織中的免疫組化結(jié)果；B：LRRC61在低級別膠質(zhì)瘤組織中的免疫組化結(jié)果；C：LRRC61在高級別膠質(zhì)瘤組織中的免疫組化結(jié)果；D：SEL1L3在在正常大腦皮層組織中的免疫組化結(jié)果；E：SEL1L3在低級別膠質(zhì)瘤組織組織中的免疫組化結(jié)果；F：SEL1L3在高級別膠質(zhì)瘤組織中的免疫組化結(jié)果，圖片數(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交集結(jié)果

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)新團(tuán)隊(duì)發(fā)展計(jì)劃”項(xiàng)目(編號IRT1230)資助[Supported by Program for Changjiang Scholars and Innovative Research Team in University (No. IRT1230)]

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

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

雷常貴和賈學(xué)淵并列第一作者。

孫文靖，博士，教授，研究方向：醫(yī)學(xué)遺傳學(xué)。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

(責(zé)任編委: 何順民)