基于生物信息學(xué)挖掘乳腺癌關(guān)鍵代謝基因及生存預(yù)后分析*

羅業(yè)浩，唐秀松，許棟涵，呂 挺，游香華，龐宇舟，李仁鋒

1 廣西中醫(yī)藥大學(xué)壯醫(yī)藥學(xué)院，廣西 南寧 530200； 2 澳門科技大學(xué)中醫(yī)藥學(xué)院，澳門 999078；3 湖北民族大學(xué)，湖北 恩施 445000； 4 珠海市中西醫(yī)結(jié)合醫(yī)院，廣東 珠海 519000；5 廣西中醫(yī)藥大學(xué)第一附屬醫(yī)院,廣西 南寧 530023

乳腺癌（breast cancer，BC）是女性最常見的惡性腫瘤之一，在我國(guó)該病發(fā)病率占全身惡性腫瘤總數(shù)的7%～10%，每年女性乳腺癌發(fā)病例數(shù)達(dá)到16.9萬，占全球總發(fā)病數(shù)的12.25%［1-2］。目前，乳腺癌的確切發(fā)病機(jī)制尚不分清楚，醫(yī)學(xué)界普遍認(rèn)為主要與遺傳因素、基因突變、機(jī)體免疫功能下降、神經(jīng)功能異常等有關(guān)［3-4］。關(guān)于乳腺癌的生存及預(yù)后狀態(tài)，研究顯示［5］，年齡小于35 歲（≤35 者）5 年無遠(yuǎn)處轉(zhuǎn)移生存率以及5 年總體生存率分別為83.4%，88.1%；年齡大于35 歲（≥35）則分別為88.3%，88.2%。乳腺癌的預(yù)后受自身因素（如腫瘤大小、淋巴結(jié)轉(zhuǎn)移、病理分級(jí)等）、環(huán)境因素（如不良飲食生活習(xí)慣、精神心理因素等）、遺傳因素（如家族遺傳、免疫缺陷等）、分子生物學(xué)因素如雌激素受體（estrogen receptor，ER）、孕激素受體（progesterone receptor，PR）、人類表皮生長(zhǎng)因子 受 體2（human epidermal growth factor receptor-2，HER-2）、腫瘤增 殖抗原（nucleus related antigen Ki-67、Ki-67）、腫瘤蛋白P53（Recombinant Tumor Protein p53，P53）等多種因素綜合影響［6-8］。

生物信息學(xué)（bioinformatics）是在生命科學(xué)的研究中，以計(jì)算機(jī)為工具對(duì)生物信息進(jìn)行儲(chǔ)存、檢索和分析的科學(xué)［9-10］。生物信息學(xué)作為當(dāng)前自然科學(xué)的研究重點(diǎn)之一，基因表達(dá)譜分析，代謝網(wǎng)絡(luò)分析，基因芯片設(shè)計(jì)和蛋白質(zhì)組學(xué)數(shù)據(jù)分析等，逐漸成為生物信息學(xué)中新興的重要研究領(lǐng)域［11］。因此，利用生物信息學(xué)的方法尋找乳腺癌的關(guān)鍵代謝基因并進(jìn)行生存預(yù)后分析具有重要意義，為乳腺癌的研究提供一定的依據(jù)和思路。

1 資料與方法

1.1 數(shù)據(jù)獲得

1.1.1 乳腺癌轉(zhuǎn)錄組數(shù)據(jù)的獲得 乳腺癌的轉(zhuǎn)錄組數(shù)據(jù)通過TCGA 數(shù)據(jù)庫(kù)（https：//portal.gdc.cancer.gov/）進(jìn)行檢索，在TCGA 數(shù)據(jù)庫(kù)中首先進(jìn)行Case 篩選，在Case ID 選項(xiàng)中選擇組織類型為“breast”，Program 選項(xiàng)中選擇研究類型為“TCGA”，Project 選項(xiàng)中選擇“TCGA-BRCA”，其余選項(xiàng)則按照TCGA 數(shù)據(jù)庫(kù)默認(rèn)選項(xiàng)進(jìn)行。接下來進(jìn)行File 篩選，在File 選項(xiàng)中Data Category 選擇為“transcriptome profiling”，Data Type 選擇為“Gene Expression Quantification”，Workflow Type 選擇為“HTSeq-FPKM”。按上述方法選擇好數(shù)據(jù)后進(jìn)行數(shù)據(jù)下載即可得到乳腺癌的轉(zhuǎn)錄組數(shù)據(jù)。

1.1.2 乳腺癌臨床數(shù)據(jù)的獲得 乳腺癌臨床數(shù)據(jù)的獲得同樣選擇TCGA 數(shù)據(jù)庫(kù)進(jìn)行篩選，按照“1.1.1”項(xiàng)中的步驟，在TCGA 數(shù)據(jù)庫(kù)中首先進(jìn)行Case 篩選，在Case ID 選項(xiàng)中選擇組織類型為“breast”，Program選項(xiàng)中選擇研究類型為“TCGA”，Project 選項(xiàng)中選擇“TCGA-BRCA”，其余選項(xiàng)則按照TCGA 數(shù)據(jù)庫(kù)默認(rèn)選項(xiàng)進(jìn)行。接下來進(jìn)行File篩選，在File 選項(xiàng)中Data Category 選擇為“clinical”，Data Format 選擇為“bcr xml”，選擇好后將數(shù)據(jù)下載即可獲得乳腺癌的臨床數(shù)據(jù)。

1.1.3 乳腺癌代謝相關(guān)基因的表達(dá)量數(shù)據(jù)獲得 首先通過GSEA 數(shù)據(jù)庫(kù)（http：//software.broadinstitute.org/gsea/index.jsp）查找相關(guān)代謝基因，下載“KEGG gene sets，gene symbols”文件，結(jié)合TCGA 數(shù)據(jù)庫(kù)所得的乳腺癌基因進(jìn)行匹配，最后獲得乳腺癌的代謝基因。

1.2 數(shù)據(jù)處理通過Perl 將獲得的乳腺癌轉(zhuǎn)數(shù)組數(shù)據(jù)和臨床數(shù)據(jù)分別運(yùn)用merge.pl、symbol.pl、getClinical.pl 等相關(guān)腳本文件運(yùn)行處理，分別獲得轉(zhuǎn)錄組數(shù)據(jù)矩陣和臨床數(shù)據(jù)矩陣。將GSEA 數(shù)據(jù)庫(kù)所得的代謝基因結(jié)合TCGA 乳腺癌的相關(guān)基因運(yùn)用相關(guān)腳本文件進(jìn)行提取乳腺癌的代謝基因表達(dá)量，運(yùn)用R×64 3.6.1 軟件進(jìn)行差異分析。獲取差異基因后，結(jié)合整理好的臨床數(shù)據(jù)運(yùn)用相關(guān)腳本文件進(jìn)行TCGA 代謝基因和生存時(shí)間合并，再次運(yùn)用R×64 3.6.1 軟件篩選預(yù)后相關(guān)的代謝基因、Lasso 回歸模型構(gòu)建、生存分析、風(fēng)險(xiǎn)曲線繪制、獨(dú)立預(yù)后分析等。

2 結(jié)果

2.1 乳腺癌轉(zhuǎn)錄組數(shù)據(jù)整理結(jié)果運(yùn)用Perl 處理后，獲得乳腺癌的轉(zhuǎn)錄組數(shù)據(jù)中，共有1222 個(gè)個(gè)樣品，其中正常樣品有113個(gè)，腫瘤樣品有1109個(gè)，進(jìn)行基因ID的轉(zhuǎn)換后得到56753個(gè)相關(guān)基因。見表1。

表1 乳腺癌轉(zhuǎn)錄組樣本相關(guān)基因

2.2 乳腺癌臨床數(shù)據(jù)整理結(jié)果在TCGA 數(shù)據(jù)庫(kù)下載好乳腺癌臨床數(shù)據(jù)后，通過Perl 運(yùn)用相關(guān)腳本文件對(duì)數(shù)據(jù)進(jìn)行整理，共得到443 個(gè)臨床數(shù)據(jù)，結(jié)果見表2。

表2 乳腺癌臨床數(shù)據(jù)

2.3 乳腺癌臨床數(shù)據(jù)整理結(jié)果將GSEA 數(shù)據(jù)庫(kù)所得的代謝基因和TCGA 數(shù)據(jù)庫(kù)所得的乳腺癌基因通過相關(guān)腳本文件處理后，得到乳腺癌代謝基因的表達(dá)量，共得到944個(gè)數(shù)據(jù)，見表3。運(yùn)用R×64 3.6.1軟件進(jìn)行差異分析，結(jié)果見表4，圖1。

圖1 代謝基因差異分析火山圖

表3 乳腺癌代謝基因表達(dá)量

表4 代謝基因差異分析結(jié)果

2.4 TCGA 代謝基因和生存時(shí)間合并結(jié)果首先對(duì)乳腺癌的臨床數(shù)據(jù)進(jìn)行整理［12］，去掉生存時(shí)間（futime）為空的或生存時(shí)間（futime）小于30 天的樣品，同樣去掉生存狀態(tài)（fustat）為空的樣品，結(jié)合“2.3”項(xiàng)所得的差異代謝基因運(yùn)用相關(guān)腳本文件進(jìn)行生存時(shí)間合并，結(jié)果得到119 個(gè)相關(guān)基因。見附表1。

2.5 Lasso 模型構(gòu)建結(jié)果對(duì)所得的預(yù)后相關(guān)代謝基因運(yùn)用R×64 3.6.1 軟件中g(shù)lmnet 包及survival包構(gòu)建Lasso模型［13-14］，其中risk score值的計(jì)算公式為：risk score=（Coefficient-Gene1×expression of Gene1）+（Coefficient-Gene2×expression of Gene2）+…+（Coefficient genen×expression Genen）。見附表2。

2.6 預(yù)后相關(guān)代謝基因篩選結(jié)果對(duì)“2.4”項(xiàng)所得的差異代謝基因運(yùn)用單因素Cox 回歸、Pvalue＜0.05進(jìn)行篩選，HR值大于1表明該基因?yàn)楦唢L(fēng)險(xiǎn)基因，P＜0.05 表示該基因與預(yù)后相關(guān)，運(yùn)用R×64 3.6.1軟件中survival包進(jìn)行預(yù)后相關(guān)代謝基因篩選［15-16］。見表5，圖2。

圖2 預(yù)后相關(guān)代謝基因森林圖

表5 預(yù)后相關(guān)代謝基因篩選結(jié)果

2.7 生存分析結(jié)果根據(jù)“2.5”項(xiàng)所得的風(fēng)險(xiǎn)值（risk score），根據(jù)風(fēng)險(xiǎn)值的中位值大小判斷病人的高低風(fēng)險(xiǎn)后進(jìn)行生存分析，以P＜0.05 表示該病人在高低風(fēng)險(xiǎn)組間，其生存具有差異［17］，運(yùn)用R×64 3.6.1 軟件中survival 包及survminer 包繪制生存曲線圖。如圖3。

圖3 生存曲線圖

2.8 風(fēng)險(xiǎn)分析結(jié)果根據(jù)“2.5”項(xiàng)所得生存風(fēng)險(xiǎn)曲線運(yùn)用R×64 3.6.1軟件中pheatmap包繪制風(fēng)險(xiǎn)分析［18］，結(jié)果顯示隨著風(fēng)險(xiǎn)值的增大，死亡的病人逐漸增多。見圖4—6。

圖4 風(fēng)險(xiǎn)熱圖

圖5 風(fēng)險(xiǎn)曲線圖

圖6 生存狀態(tài)圖

2.9 獨(dú)立預(yù)后分析結(jié)果根據(jù)“2.7”項(xiàng)所得風(fēng)險(xiǎn)分析結(jié)果，運(yùn)用R×64 3.6.1 軟件中survival 包分別進(jìn)行單因素及多因素獨(dú)立預(yù)后分析，以風(fēng)險(xiǎn)值（risk score）P＜0.001 表示該風(fēng)險(xiǎn)值可以作為獨(dú)立預(yù)后因素。見圖7—8。

圖7 多因素獨(dú)立預(yù)后分析森林圖

3 討論

本次研究共得到了3 個(gè)與乳腺癌相關(guān)的代謝基因，分別為POLR2K、NMNAT2、SUCLA2。POLR2K 指的是RNA聚合酶Ⅱ亞基K，該基因是一種蛋白質(zhì)編碼基因，其相關(guān)途徑包括RNA聚合酶I啟動(dòng)子逃逸和RNA 聚合酶Ⅲ轉(zhuǎn)錄起始，POLR2K 編碼RNA 聚合酶Ⅱ的最小亞基，它是三種RNA 聚合酶的常見亞基［19］。有研究推測(cè)POLR2K 的上調(diào)可能促進(jìn)聚合酶Ⅲ的組裝，從而促進(jìn)細(xì)胞增殖和癌癥發(fā)展，并且可作為乳腺癌的危險(xiǎn)因素［20-22］。NMNAT2指的是煙酰胺核苷酸腺苷酸轉(zhuǎn)移酶2，該基因是一種蛋白質(zhì)編碼基因，其主要功能是催化β-煙酰胺單核苷酸和三磷酸腺苷（adenosine-triphosphate，ATP）合成煙酰胺腺嘌呤二核苷酸［23-24］。研究表明，NMNAT2 的表達(dá)與腫瘤的浸潤(rùn)深度和腫瘤分期（tumor node metastasis，TNM）有關(guān)［25］。SUCLA2是一種蛋白質(zhì)編碼基因，這是三羧酸循環(huán)的重要組成部分。SCS-A 水解ATP，將琥珀酸酯轉(zhuǎn)化為琥珀酰-CoA［26］。目前研究發(fā)現(xiàn)，該基因的突變是線粒體DNA 耗竭綜合征的主要原因，腦和骨骼肌是主要受累器官［27］。遺憾的是，目前尚未見相關(guān)文獻(xiàn)報(bào)道其與癌癥或乳腺癌相關(guān)，這為我們后期的研究提供一種思路。

在生存分析（圖3）中，共有545個(gè)樣品被納入研究，結(jié)果發(fā)現(xiàn)隨著時(shí)間（year）的延長(zhǎng)，生存的患者數(shù)量逐漸減少，高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組最長(zhǎng)生存時(shí)間均為24年。在生存風(fēng)險(xiǎn)熱圖（圖5）中，我們發(fā)現(xiàn)POLR2K 基因是最為明顯的高表達(dá)基因，提示該基因可能與乳腺癌的生存預(yù)后密切相關(guān)。在獨(dú)立預(yù)后因素分析中，其中單因素獨(dú)立預(yù)后分析（圖8）發(fā)現(xiàn)年齡（age）、狀態(tài)（stage）、TNM 分期這3 個(gè)因素可作為單因素的獨(dú)立預(yù)后，而在多因素獨(dú)立預(yù)后分析中，年齡（age）與多因素預(yù)后密切相關(guān)。

圖8 單因素獨(dú)立預(yù)后分析森林圖

盡管我們已經(jīng)使用生物信息學(xué)技術(shù)對(duì)大量樣本鑒定出了可能影響乳腺癌預(yù)后的潛在代謝基因，但本次研究仍然存在某些局限性。首先，樣本有些部分缺乏一些臨床隨訪信息。其次，僅使用生物信息學(xué)分析獲得的結(jié)果是不夠的，需要通過實(shí)驗(yàn)驗(yàn)證來確認(rèn)。因此，需要對(duì)更大的樣本進(jìn)行進(jìn)一步的遺傳和實(shí)驗(yàn)研究并進(jìn)行實(shí)驗(yàn)驗(yàn)證。總之，本次研究基于生物信息學(xué)挖掘了乳腺癌的相關(guān)代謝基因，并對(duì)關(guān)鍵代謝基因進(jìn)行了生存預(yù)后分析，為乳腺癌的治療研究提供了一種參照方法和思路。