基于生物信息學挖掘乳腺癌關鍵代謝基因及生存預后分析*

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

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

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

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

1 資料與方法

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

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

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

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

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

2 結果

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

表1 乳腺癌轉錄組樣本相關基因

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

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

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

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

表3 乳腺癌代謝基因表達量

表4 代謝基因差異分析結果

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

2.5 Lasso 模型構建結果對所得的預后相關代謝基因運用R×64 3.6.1 軟件中glmnet 包及survival包構建Lasso模型［13-14］，其中risk score值的計算公式為：risk score=（Coefficient-Gene1×expression of Gene1）+（Coefficient-Gene2×expression of Gene2）+…+（Coefficient genen×expression Genen）。見附表2。

2.6 預后相關代謝基因篩選結果對“2.4”項所得的差異代謝基因運用單因素Cox 回歸、Pvalue＜0.05進行篩選，HR值大于1表明該基因為高風險基因，P＜0.05 表示該基因與預后相關，運用R×64 3.6.1軟件中survival包進行預后相關代謝基因篩選［15-16］。見表5，圖2。

圖2 預后相關代謝基因森林圖

表5 預后相關代謝基因篩選結果

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

圖3 生存曲線圖

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

圖4 風險熱圖

圖5 風險曲線圖

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

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

圖7 多因素獨立預后分析森林圖

3 討論

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

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

圖8 單因素獨立預后分析森林圖

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