聯(lián)合CGGA-TCGA數(shù)據(jù)庫挖掘膠質(zhì)瘤顯著生存相關(guān)基因

呂俊杰，張一驊，鄧 貴，劉佳欣，呂正禹，任付賓，陳麗娜

(1.哈爾濱醫(yī)科大學(xué) 生物信息科學(xué)與技術(shù)學(xué)院，哈爾濱 150081；2.哈爾濱醫(yī)科大學(xué)附屬第四醫(yī)院，哈爾濱 150081)

根據(jù)世界衛(wèi)生組織對(duì)神經(jīng)系統(tǒng)腫瘤的分類，腦膠質(zhì)瘤被分為9類[1]：星形細(xì)胞性腫瘤、少突膠質(zhì)性腫瘤、室管膜腫瘤、混合性膠質(zhì)瘤、脈絡(luò)叢腫瘤、來源未定的神經(jīng)上皮腫瘤、神經(jīng)元及神經(jīng)元與膠質(zhì)細(xì)胞混合性腫瘤、松果體瘤和胚胎性腫瘤。其中星形細(xì)胞性腫瘤中的膠質(zhì)母細(xì)胞瘤(GBM)惡性程度最高，GBM的5年生存率非常低。許多工作都揭示出一些基因的異常表達(dá)影響膠質(zhì)瘤[2-7]，也有研究從大樣本中識(shí)別到膠質(zhì)瘤的預(yù)后基因[8-12]。然而膠質(zhì)瘤作為惡性腫瘤，復(fù)雜程度高，通常預(yù)后效果不好。值得注意的是，目前雖然對(duì)膠質(zhì)瘤的診斷水平還很低，但是病人的生存時(shí)間個(gè)體差異很大[13-15]。所以嘗試從病人的生存時(shí)間差異來進(jìn)行生存分析，從而挖掘膠質(zhì)瘤生存相關(guān)基因及其可能的調(diào)控機(jī)制。此外，預(yù)后基因是否可以用于不同人群，這個(gè)工作開展得較少，因此有必要對(duì)預(yù)后基因進(jìn)行不同人群驗(yàn)證。

TCGA癌癥基因組圖譜為癌癥研究提供了豐富的資源，但是大部分樣本來源于國(guó)外，所得結(jié)果在中國(guó)推廣具有一定的差異。CGGA是由北京市神經(jīng)外科研究所、首都醫(yī)科大學(xué)附屬北京天壇醫(yī)院江濤教授團(tuán)隊(duì)建立的中國(guó)腦膠質(zhì)瘤基因組圖譜數(shù)據(jù)庫，因此，通過CGGA聯(lián)合TCGA數(shù)據(jù)庫來挖掘膠質(zhì)瘤生存相關(guān)基因，并分析其表達(dá)調(diào)控機(jī)制，對(duì)膠質(zhì)瘤尤其是中國(guó)人群膠質(zhì)瘤患者的預(yù)后分析具有重要的研究?jī)r(jià)值。

1 材料與方法

1.1 數(shù)據(jù)來源

下載CGGA數(shù)據(jù)庫mRNAseq_693 和mRNAseq_325 兩個(gè)批次的RNA表達(dá)譜數(shù)據(jù)，以及相應(yīng)的臨床數(shù)據(jù)；從哈爾濱醫(yī)科大學(xué)附屬第四醫(yī)院獲得滿足要求的13例樣本mRNAseq_13與臨床數(shù)據(jù)；TCGA中GBM和LGG共計(jì)1 162個(gè)樣本的表達(dá)譜數(shù)據(jù)集與臨床信息文件。

1.2 方法

對(duì)mRNAseq_693 和mRNAseq_325兩個(gè)批次的基因取交集，所得結(jié)果再與mRNAseq_13取并集，引用Biocoductor的sva包進(jìn)行批次校正，同時(shí)合并臨床數(shù)據(jù)。根據(jù)生存時(shí)間將樣品分為兩組，一組生存時(shí)間小于1年并且生存狀態(tài)為死亡，一組生存時(shí)間大于3年，基于R軟件，引用limma包的經(jīng)典貝葉斯方法對(duì)兩組樣品進(jìn)行差異分析，同時(shí)采用KM方法和COX方法進(jìn)行生存分析，兩種生存分析方法中均選取顯著性P值為0.001，同時(shí)計(jì)算高低表達(dá)組的5年生存率，判斷5年生存率的差異是否大于0.200；以生存時(shí)間1到3年的樣本，TCGA中膠質(zhì)瘤樣品為驗(yàn)證集進(jìn)行驗(yàn)證：對(duì)CGGA生存時(shí)間1到3年的樣本與TCGA中膠質(zhì)瘤樣品進(jìn)行同樣的KM與COX方法雙重生存分析過濾，得到生存顯著相關(guān)基因集，與差異表達(dá)基因集取交集得到共有的顯著生存相關(guān)基因集，引用survivalROC包，方法設(shè)置為KM方法，預(yù)測(cè)時(shí)間選取3年，對(duì)共有的顯著生存相關(guān)基因集繪制ROC曲線，得到AUC值大于0.700的基因集。對(duì)得到的集中基因進(jìn)行文獻(xiàn)驗(yàn)證，根據(jù)基因表達(dá)的中位值進(jìn)行高低表達(dá)分組，繪制風(fēng)險(xiǎn)曲線，采用單因素COX分析與多因素COX回歸進(jìn)行獨(dú)立預(yù)后分析；應(yīng)用GSEA軟件對(duì)目標(biāo)基因進(jìn)行富集分析，準(zhǔn)備目標(biāo)基因的表達(dá)數(shù)據(jù)文件與表型數(shù)據(jù)文件，表達(dá)數(shù)據(jù)集文件存儲(chǔ)基因表達(dá)量，表型數(shù)據(jù)文件存儲(chǔ)按照基因的中位值，將目標(biāo)基因分為高低風(fēng)險(xiǎn)兩組的標(biāo)簽l和h，得到目標(biāo)的富集結(jié)果；引用limma包采用皮爾遜相關(guān)系數(shù)法，計(jì)算基因間表達(dá)量的相關(guān)性，相關(guān)性系數(shù)cor過濾值設(shè)置為0.500，統(tǒng)計(jì)學(xué)P值的過濾標(biāo)準(zhǔn)設(shè)置為0.001，挑選正相關(guān)最顯著的20個(gè)基因(相關(guān)性系數(shù)cor為正)、負(fù)相關(guān)最顯著的20個(gè)基因(相關(guān)性系數(shù)cor為負(fù))，以樣品為橫坐標(biāo)、基因?yàn)榭v坐標(biāo)，引用R包pheatmap繪制熱圖，最后基于TCGA Spliceseq數(shù)據(jù)庫進(jìn)行選擇性剪接事件分析。

2 結(jié)果與分析

2.1 差異分析

數(shù)據(jù)分組：按生存時(shí)間小于1年并且生存狀態(tài)為死亡的樣本，與生存時(shí)間大于3年的樣本分為兩組，得到兩組分別為65個(gè)(shorter survival)與47個(gè)(longer survival)的樣本，具體信息如表1所示。

表1 樣本臨床信息

應(yīng)用R統(tǒng)計(jì)軟件limma包的經(jīng)典貝葉斯方法進(jìn)行差異分析，foldChange設(shè)置為0.500，即滿足差異倍數(shù)為1倍，Benjamini-Hochberg方法校正后padj值設(shè)置為0.050，得到10個(gè)padj值小于0.050的顯著差異基因PLAT、IGFBP2、AC062021.1、COL4A1、COL4A2、KCNJ11、TNFRSF12A、BCAT1、TAGLN2、SERPINH1。具體如表2所示。

表2 顯著差異基因

2.2 生存分析

采用KM方法和COX方法進(jìn)行雙重生存分析，發(fā)現(xiàn)這10個(gè)基因全部和生存顯著相關(guān)，具體如表3所示。

表3 生存顯著相關(guān)基因

對(duì)CGGA中生存時(shí)間1到3年的樣本與TCGA中膠質(zhì)瘤樣品進(jìn)行同樣的生存分析，得到生存顯著相關(guān)基因集，與這10個(gè)基因取交集得到共有的7個(gè)顯著生存相關(guān)基因PLAT、IGFBP2、AC062021.1、COL4A2、KCNJ11、BCAT1和SERPINH1。對(duì)得到顯著生存相關(guān)基因采用R統(tǒng)計(jì)軟件并引用survivalROC包，繪制ROC曲線，得到曲線下面積大于0.700的4個(gè)基因PLAT、IGFBP2、BCAT1和SERPINH1，具體見表4。

表4 AUC值大于0.700的生存顯著相關(guān)基因

分別對(duì)PLAT、IGFBP2、BCAT1和SERPINH1進(jìn)行單因素COX回歸分析與多因素COX獨(dú)立預(yù)后分析，結(jié)果如表5所示：在單因素與多因素COX分析中，4個(gè)基因P值均小于0.001，說明PLAT、IGFBP2、BCAT1和SERPINH1都可以作為獨(dú)立預(yù)后因子。

表5 獨(dú)立預(yù)后分析

進(jìn)一步經(jīng)文獻(xiàn)驗(yàn)證，PLAT、IGFBP2和BCAT1均與膠質(zhì)瘤生存顯著相關(guān)[15-17]，經(jīng)文獻(xiàn)證實(shí)均可以做獨(dú)立預(yù)后因子，繼續(xù)挖掘SERPINH1與膠質(zhì)瘤的生存相關(guān)性，根據(jù)SERPINH1的表達(dá)中位值進(jìn)行高低表達(dá)分組，繪制風(fēng)險(xiǎn)曲線，如圖1(a)CGGA中SERPINH1風(fēng)險(xiǎn)曲線。從圖1(b)TCGA中SERPINH1風(fēng)險(xiǎn)曲線可以看出，高低表達(dá)組的5年生存率顯著不同，低風(fēng)險(xiǎn)組5年生存率約80%，高風(fēng)險(xiǎn)組約30%，說明SERPINH1與膠質(zhì)瘤生存顯著相關(guān)。

圖1 CGGA與TCGA中SERPINH1風(fēng)險(xiǎn)曲線

2.3 ROC曲線評(píng)估預(yù)測(cè)能力

分別基于CGGA和TCGA數(shù)據(jù)庫繪制SERPINH1基因的3年和5年的ROC曲線，評(píng)估SERPINH1的預(yù)測(cè)能力，AUC值介于0.500和1.000之間，如圖2(a)所示：在CGGA中3年和5年的AUC值均大于0.700。由圖2(b)可知：在TCGA中5年生存率預(yù)測(cè)準(zhǔn)確性為0.691，3年生存率預(yù)測(cè)準(zhǔn)確性可達(dá)0.807，說明SERPINH1可以作為獨(dú)立預(yù)后因子。

圖2 ROC曲線評(píng)估

2.4 GSEA富集分析

應(yīng)用GSEA_4.0.0進(jìn)行基因集富集分析，對(duì)h組與l組的GO功能類按照P值大小進(jìn)行排序，選取 h 組最顯著的前 5 個(gè) GO 功能類與 l 組最顯著的前 5 個(gè) GO 功能類進(jìn)行繪圖。圖3展示的是SERPINH1主要顯著富集于哪些功能類。發(fā)現(xiàn)這些基因主要富集于代謝、細(xì)胞外基質(zhì)、多細(xì)胞有機(jī)大分子代謝過程、谷氨酸能突觸傳遞、后腦細(xì)胞分化、突觸前膜、去小腦皮層形成等過程，說明這些基因可能參與調(diào)控腦細(xì)胞形成與分化相關(guān)通路。已有研究表明這些通路的異常改變導(dǎo)致了膠質(zhì)瘤的發(fā)生[18]。

圖3 GSEA 富集分析結(jié)果展示

2.5 共表達(dá)分析

橫坐標(biāo)按SERPINH1的表達(dá)量從小到大進(jìn)行排序，即從左到右SERPINH1的表達(dá)量越來越高，引用R包pheatmap繪制熱圖(圖4)。從圖4可以看出，前20個(gè)基因隨著目標(biāo)基因的表達(dá)增高，表達(dá)量也越來越高，后20個(gè)基因的表達(dá)量則隨著目標(biāo)基因的表達(dá)量增高而減少。從圖4還可以看出，SERPINH1與IGFBP2顯著正相關(guān)，而IGFBP2是已知的膠質(zhì)瘤致病基因[11-13]。

圖4 相關(guān)性熱圖

基于通路富集發(fā)現(xiàn)這些基因主要富集于代謝、谷氨酸能突觸傳遞、后腦細(xì)胞分化、突觸前膜、去小腦皮層形成等過程，說明這些基因參與調(diào)控腦細(xì)胞發(fā)育形成與分化過程。這些通路在膠質(zhì)瘤中發(fā)生了異常改變，說明本文的研究結(jié)果與已知研究結(jié)果相一致[18]。可能是這些基因的異常表達(dá)導(dǎo)致通路的異常，從而導(dǎo)致疾病的發(fā)生。

2.6 選擇性剪接事件分析

因?yàn)檫x擇性剪接是腫瘤的關(guān)鍵致病因素之一，大約95%以上的基因會(huì)發(fā)生可變剪接[19-20]，所以，基因SERPINH1的高表達(dá)與生存相關(guān)，可能是因?yàn)樗目勺兗艚邮录c生存顯著相關(guān)，并且某個(gè)剪接事件是腫瘤發(fā)生的高風(fēng)險(xiǎn)因素，所以基于TCGA Spliceseq數(shù)據(jù)庫對(duì)SERPINH1基因的選擇性剪接進(jìn)行分析，從https://bioinformatics.mdanderson.org/TCGASpliceSeq/下載LGG和GBM的選擇性剪接數(shù)據(jù)，該網(wǎng)站提取了TCGA各種癌癥的選擇性剪接數(shù)據(jù)，選取樣品PSI值大于75%(Percentage of samples with PSI value of 75)、最小PSI范圍(樣本之間的增量)[Minumun PSI range(delta across samples)]和最小PSI標(biāo)準(zhǔn)偏差(Minumum PSI standard deviation)都設(shè)置為0，提取出SERPINH1的數(shù)據(jù)，利用R包impute對(duì)空值進(jìn)行補(bǔ)缺，PSI均值小于0.05的過濾掉，選擇性剪接在所有樣品中波動(dòng)太小的(sd值小于0.01)也刪掉。發(fā)現(xiàn)SERPINH1發(fā)生的選擇性剪接事件類型為AP、AT和ES，分別對(duì)SERPINH1選擇性剪接事件進(jìn)行單因素COX分析，結(jié)果如表6所示：和生存顯著相關(guān)的剪接事件為發(fā)生于1、4.1和2.2號(hào)外顯子的AP和ES。其中發(fā)生于4.1號(hào)外顯子的AP(SERPINH1|17866|AP)是高風(fēng)險(xiǎn)選擇性剪接事件(HR>1)。SERPINH1的異常剪接可能是導(dǎo)致膠質(zhì)瘤發(fā)生或者不良預(yù)后的重要因素。

表6 生存顯著相關(guān)的選擇性剪接事件

3 討論與結(jié)論

由于膠質(zhì)瘤患者的生存時(shí)間個(gè)體差異較大，所以嘗試從病人的生存時(shí)間差異來進(jìn)行生存分析。采用KM與COX方法進(jìn)行雙重生存分析，得到同時(shí)滿足兩種方法顯著性P值均小于0.05的10個(gè)基因：PLAT、IGFBP2、AC062021.1、COL4A1、COL4A2、KCNJ11、TNFRSF12A、BCAT1、TAGLN2和SERPINH1。KM方法把基因看成離散變量，COX方法則視基因?yàn)檫B續(xù)變量，從結(jié)果可以看出，這10個(gè)基因基于KM方法的顯著性P值均為0，小于0.050；基于COX方法的顯著性P值也都遠(yuǎn)遠(yuǎn)小于0.050，可見這些基因都與生存顯著相關(guān)；經(jīng)TCGA驗(yàn)證得到兩個(gè)數(shù)據(jù)庫共有的7個(gè)顯著生存相關(guān)基因PLAT、IGFBP2、AC062021.1、COL4A2、KCNJ11、BCAT1和SERPINH1。

通過引用survivalROC包繪制ROC曲線評(píng)估7個(gè)基因的預(yù)后預(yù)測(cè)性能，發(fā)現(xiàn)PLAT、IGFBP2、BCAT1和SERPINH1的AUC值均大于0.7，其他3個(gè)基因的AUC值也都大于0.6，說明獨(dú)立預(yù)后預(yù)測(cè)性能較好。通過單因素COX分析與多因素COX回歸分析，發(fā)現(xiàn)PLAT、IGFBP2、BCAT1和SERPINH1的顯著性P值都小于0.010，說明這4個(gè)基因均可以作為獨(dú)立預(yù)后因子。

BCAT1、PLAT和IGFBP2均有文獻(xiàn)證實(shí)可以作為膠質(zhì)瘤獨(dú)立預(yù)后因子，所以繼續(xù)對(duì)SERPINH1進(jìn)行獨(dú)立預(yù)后分析，根據(jù)SERPINH1表達(dá)的高低分組，基于CGGA與TCGA數(shù)據(jù)庫繪制基因風(fēng)險(xiǎn)曲線可以看出：高低表達(dá)組的5年生存率顯著不同，低風(fēng)險(xiǎn)組5年生存率約80%，高風(fēng)險(xiǎn)組約30%；從生存曲線可以看到，SERPINH1在CGGA數(shù)據(jù)庫中的分類性能更優(yōu)于在TCGA中的分類性能。這可能預(yù)示著SERPINH1更適合作為中國(guó)膠質(zhì)瘤人群的獨(dú)立預(yù)后因子。

通過SERPINH1基因在CGGA與TCGA中的3年和5年生存率的ROC曲線可以看出，在CGGA中AUC值都大于0.700，在TCGA中5年生存率曲線AUC值為0.691接近于0.700，而3年生存率曲線下面積為0.807，說明基因SERPINH1可以作為獨(dú)立預(yù)后因子，預(yù)測(cè)性能比較好。發(fā)現(xiàn)SERPINH1在CGGA中3年和5年的生存率預(yù)測(cè)準(zhǔn)確性近似，在TCGA中3年和5年的生存率預(yù)測(cè)準(zhǔn)確性差別比較明顯，這種差異可能是因?yàn)槿朔N的差異造成的，尚需進(jìn)一步驗(yàn)證。

通過基因富集分析與共表達(dá)分析，發(fā)現(xiàn)PLAT、IGFBP2、BCAT1和SERPINH1主要富集于代謝與腦發(fā)育形成等過程。已有研究表明這些通路在膠質(zhì)瘤中發(fā)生了異常改變，可能是這些基因的異常表達(dá)導(dǎo)致通路異常，從而導(dǎo)致疾病的發(fā)生。并且我們發(fā)現(xiàn)IGFBP2與SERPINH1顯著正相關(guān)，IGFBP2是膠質(zhì)瘤已知的致病基因，這說明SERPINH1可能通過與IGFBP2互作共同影響膠質(zhì)瘤的預(yù)后，后續(xù)工作將進(jìn)一步基于OMIM數(shù)據(jù)庫挖掘已知膠質(zhì)瘤致病基因與PLAT、IGFBP2、BCAT1和SERPINH1的關(guān)系，并基于STRING數(shù)據(jù)庫，繪制基因間的蛋白質(zhì)互作網(wǎng)絡(luò)，并對(duì)網(wǎng)絡(luò)進(jìn)行分析研究。

為進(jìn)一步挖掘?qū)е耂ERPINH1和膠質(zhì)瘤顯著生存相關(guān)的分子調(diào)控機(jī)制，基于TCGASpliceseq數(shù)據(jù)庫分析了基因SERPINH1的選擇性剪接事件，發(fā)現(xiàn)SERPINH1選擇性剪接事件主要是ES和AP，經(jīng)過單因素COX回歸分析發(fā)現(xiàn)發(fā)生2.2號(hào)外顯子的ES剪接事件，1號(hào)和4.1號(hào)外顯子的AP事件與生存顯著相關(guān)，其中SERPINH1發(fā)生1號(hào)外顯子的AP剪接事件是高風(fēng)險(xiǎn)選擇性剪接事件，說明基因SERPINH1的異常剪接或者調(diào)控SERPINH1剪接調(diào)控因子的異常會(huì)導(dǎo)致膠質(zhì)瘤的發(fā)生或不良預(yù)后。后續(xù)可以構(gòu)建SERPINH1剪接調(diào)控網(wǎng)絡(luò)，研究生存相關(guān)基因異常剪接以及剪接調(diào)控因子的異常對(duì)膠質(zhì)瘤的致病性，進(jìn)一步深入研究膠質(zhì)瘤的選擇性剪接調(diào)控機(jī)制。同時(shí)，也發(fā)現(xiàn)了SERPINH1在CGGA和TCGA中預(yù)測(cè)能力的差異，后續(xù)將進(jìn)一步研究基于不同人種的膠質(zhì)瘤獨(dú)立預(yù)后因子的分子機(jī)制。

結(jié)合差異分析、生存分析、共表達(dá)分析、ROC曲線評(píng)估預(yù)測(cè)能力，并通過TCGA驗(yàn)證，最終得到PLAT、IGFBP2、BCAT1和SERPINH1等4個(gè)可以作為獨(dú)立預(yù)后標(biāo)志因子的基因。分析顯示這4個(gè)基因不僅可以作為獨(dú)立的預(yù)后因子，并且這4個(gè)基因在不同人群中都具有較好的獨(dú)立預(yù)后預(yù)測(cè)能力，這說明本文找到的預(yù)后基因具有潛在的應(yīng)用價(jià)值。對(duì)SERPINH1選擇性剪接事件進(jìn)行分析，結(jié)果顯示選擇性剪接事件可能是SERPINH1影響膠質(zhì)瘤患者預(yù)后的重要因素，這為進(jìn)一步理解膠質(zhì)瘤預(yù)后基因的特征提供了依據(jù)。