李偉華 張廣鳳 馬驪驪 何品 王雯 李頂夫

乳腺癌是全球女性最常見(jiàn)的惡性腫瘤之一，其發(fā)病率逐年上升，2020年已超過(guò)肺癌成為全球第一大癌癥[1]。乳腺癌是一組高度異質(zhì)性的疾病，每個(gè)病例的臨床表現(xiàn)、形態(tài)學(xué)、分子特性、治療效果及預(yù)后差異很大[2]。臨床上，常用腫瘤分期、組織學(xué)分級(jí)和分子亞型作為評(píng)估乳腺癌患者的預(yù)后因素。然而，這些臨床病理特征并不能準(zhǔn)確地預(yù)測(cè)患者預(yù)后的信息。這可能導(dǎo)致對(duì)患者預(yù)后的判斷不準(zhǔn)確，一些低風(fēng)險(xiǎn)患者可能接受不必要或過(guò)度的治療，而其他高風(fēng)險(xiǎn)患者可能因治療不當(dāng)而面臨復(fù)發(fā)或轉(zhuǎn)移[3]。因此，迫切需要尋找新的分子標(biāo)記物來(lái)預(yù)測(cè)乳腺癌患者的預(yù)后，從而有利于患者的精確治療。

腫瘤通過(guò)代謝為腫瘤細(xì)胞提供重要的能量來(lái)源，如以“有氧糖酵解”或“Warburg效應(yīng)”為代表的糖代謝，可以增加大分子和中間產(chǎn)物的合成以維持腫瘤的增殖[4]。在腫瘤細(xì)胞內(nèi)，產(chǎn)生了許多代謝的變化以滿足腫瘤的能量和合成要求，這也成為腫瘤的一個(gè)重要特征[5]。本研究構(gòu)建了一個(gè)結(jié)合多個(gè)乳腺癌代謝相關(guān)基因的模型來(lái)預(yù)測(cè)乳腺癌患者的預(yù)后。通過(guò)對(duì)癌癥基因組圖譜計(jì)劃(The Cancer Genome Atlas Program，TCGA)數(shù)據(jù)庫(kù)中乳腺癌代謝相關(guān)基因進(jìn)行差異分析，篩選出代謝相關(guān)的差異基因，通過(guò)lasso回歸去掉過(guò)度擬合的基因用作預(yù)后模型的構(gòu)建。對(duì)高低風(fēng)險(xiǎn)兩組進(jìn)行生存比較，對(duì)風(fēng)險(xiǎn)得分進(jìn)行了單因素和多因素Cox回歸分析，計(jì)算出每個(gè)基因的風(fēng)險(xiǎn)系數(shù)和在樣本中的表達(dá)情況。根據(jù)基因表達(dá)量分析不同風(fēng)險(xiǎn)分組的樣本的生存狀態(tài)。

1 材料與方法

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

TCGA(https：//www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga)下載了113例正常對(duì)照和1109例腫瘤樣本的mRNA表達(dá)數(shù)據(jù)，并下載與之對(duì)應(yīng)的患者臨床信息用作預(yù)后分析。從京都基因與基因組百科全書(Kyoto Encyclopedia of Genes and Genomes，KEGG)(https：//www.genome.jp/kegg/)下載41條代謝相關(guān)的通路，將其中的944個(gè)基因提取出來(lái)作為代謝相關(guān)基因。

1.2 預(yù)后相關(guān)的基因篩選

通過(guò)limma包對(duì)KEGG中提取的代謝基因進(jìn)行差異分析，F(xiàn)DR<0.05，│logFC│>1作為篩選標(biāo)準(zhǔn)，并通過(guò)Cox回歸篩選出P<0.05的預(yù)后基因。隨機(jī)抽取70%的基因作為訓(xùn)練集，其余樣本作為驗(yàn)證集，進(jìn)一步通過(guò)Lasso回歸防止過(guò)度擬合，13個(gè)代謝相關(guān)預(yù)后基因用作預(yù)后模型構(gòu)建。

1.3 構(gòu)建預(yù)后模型

利用Cox回歸得到的代謝預(yù)后相關(guān)基因權(quán)重進(jìn)行預(yù)后模型構(gòu)建，通過(guò)風(fēng)險(xiǎn)得分中位值將樣本分為高低風(fēng)險(xiǎn)兩組，并利用survival包及survminer包繪制生存曲線。利用pheatmap包進(jìn)行風(fēng)險(xiǎn)得分與生存狀態(tài)之間的熱圖繪制。

風(fēng)險(xiǎn)得分=代謝相關(guān)風(fēng)險(xiǎn)基因表達(dá)值(FPKM)×風(fēng)險(xiǎn)系數(shù)

1.4 基因集富集分析(Gene set enrichment analysis，GSEA)

GSEA可以通過(guò)將預(yù)定義的基因集與特定表型進(jìn)行比對(duì)從而找出與預(yù)定義基因集相關(guān)的表型。根據(jù)代謝相關(guān)預(yù)后基因的風(fēng)險(xiǎn)得分進(jìn)行分組，并與KEGG的基因集進(jìn)行比對(duì)，對(duì)代謝相關(guān)預(yù)后基因在乳腺癌中可能的作用機(jī)制進(jìn)行探索。

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

所有數(shù)據(jù)的統(tǒng)計(jì)分析均使用R軟件(版本4.0.2)，P值多重檢驗(yàn)校正采用Benjamini & Hochberg(BH)方法，P<0.05認(rèn)為具有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 代謝相關(guān)預(yù)后基因篩選

通過(guò)從TCGA數(shù)據(jù)庫(kù)下載得到113例正常對(duì)照和1 109例乳腺癌患者樣本的mRNA表達(dá)譜數(shù)據(jù)。從中提取KEGG通路中代謝相關(guān)的基因用作進(jìn)一步差異分析。最終篩選出168個(gè)代謝相關(guān)的差異基因(圖1)。并通過(guò)單因素Cox回歸確定與預(yù)后相關(guān)的14個(gè)基因(圖2)。

圖1 代謝相關(guān)的差異表達(dá)基因Figure 1 Metabolism-related differentially expressed genesNote：A.Differential gene expression；B.Volcano map showing the differential gene distribution.

圖2 森林圖 Cox回歸單因素分析篩選出7個(gè)高風(fēng)險(xiǎn)基因和7個(gè)保護(hù)性基因Figure 2 Forest map Cox regression single factor analysis screened out 7 high-risk genes and 7 protective genes

2.2 Lasso回歸篩選關(guān)鍵基因及預(yù)后模型的構(gòu)建

進(jìn)一步通過(guò)lasso回歸去掉過(guò)度擬合的基因，最終確定了13個(gè)基因(NMNAT2、NT5E、QPRT、UGP2、MTHFD2、TSTA3、TYMP、ALDH2、ALDH1A1、IDO1、IL4I1、INPP1和ENPP6)用作預(yù)后模型的構(gòu)建(圖3)。并計(jì)算出每個(gè)基因的風(fēng)險(xiǎn)系數(shù)(表1)。展示每個(gè)基因在樣本中的表達(dá)情況(圖4A)。根據(jù)基因表達(dá)量將樣本分為高低風(fēng)險(xiǎn)兩組，展示不同風(fēng)險(xiǎn)分組的樣本的生存狀態(tài)。隨著風(fēng)險(xiǎn)得分的增加，相同時(shí)間節(jié)點(diǎn)，患者死亡率增高(圖4B，C)。

圖3 比例風(fēng)險(xiǎn)模型構(gòu)建Figure 3 The construction of proportional risk modelNote：A.Adjust the cross-validation of parameter selection in the proportional hazard model；B.LASSO coefficient spectrum of 14 genes in breast cancer.

表1 風(fēng)險(xiǎn)基因的風(fēng)險(xiǎn)系數(shù)

圖4 高低風(fēng)險(xiǎn)分組基因表達(dá)及生存狀態(tài)Figure 4 The high-and low-risk group gene expression and survival status of breast cancer patientsNote：A.The expression of risk genes in different groups；B.The distribution of risk scores in different groups；C.The distribution of survival status of patients in different groups.

2.3 預(yù)后模型的檢驗(yàn)

通過(guò)對(duì)高低風(fēng)險(xiǎn)兩組進(jìn)行生存比較，發(fā)現(xiàn)高風(fēng)險(xiǎn)組患者的生存率更低(圖5)。并對(duì)風(fēng)險(xiǎn)得分進(jìn)行了單因素和多因素的Cox回歸分析，結(jié)果表明，利用代謝相關(guān)基因構(gòu)建的風(fēng)險(xiǎn)得分能夠作為患者的一個(gè)獨(dú)立危險(xiǎn)因素去預(yù)測(cè)患者的預(yù)后(圖6)，并利用ROC曲線下面積去檢驗(yàn)了風(fēng)險(xiǎn)得分去預(yù)測(cè)患者預(yù)后是較為準(zhǔn)確的(圖7)。

圖5 高風(fēng)險(xiǎn)組乳腺癌患者生存時(shí)間更短Figure 5 A shorter survival time in breast cancer patients from the high-risk group

圖6 單因素和多因素Cox回歸確定風(fēng)險(xiǎn)得分為獨(dú)立預(yù)后因素Figure 6 Univariate and multivariate Cox regression determined the risk score as an independent prognostic factor

圖7 風(fēng)險(xiǎn)評(píng)分對(duì)預(yù)后具有較高的預(yù)測(cè)率Figure 7 The risk score had a high predictive rate for prognosis

此外，根據(jù)風(fēng)險(xiǎn)得分進(jìn)行GSEA分析(圖8)，發(fā)現(xiàn)其顯著富集到代謝相關(guān)的通路(KEGG_ARGININE_AND_PROLINE_METABOLISM、KEGG_GLYCOLYSIS_GLUCONEOGENESIS、KEGG_PYRIMIDINE_METABOLISM、KEGG_STARCH_AND_SUCROSE_METABOLISM)。

圖8 風(fēng)險(xiǎn)得分相關(guān)基因顯著富集在代謝相關(guān)通路Figure 8 Risk score-related genes were significantly enriched in metabolic pathways

3 討論

腫瘤細(xì)胞為適應(yīng)其對(duì)能量和物質(zhì)的需求而發(fā)生了許多改變，而代謝變化成為其一個(gè)重要特征，對(duì)代謝相關(guān)基因的研究來(lái)判斷腫瘤的特征和預(yù)后被越來(lái)越多的人所重視[6]。因此，本研究把代謝相關(guān)基因作為乳腺癌的預(yù)后的判斷。本研究從TCGA數(shù)據(jù)庫(kù)下載了乳腺癌患者和對(duì)照樣本的mRNA表達(dá)譜數(shù)據(jù)。從中提取KEGG通路中代謝相關(guān)的基因用作進(jìn)一步差異分析。最終篩選出168個(gè)代謝相關(guān)的差異基因。并通過(guò)單因素Cox回歸確定與預(yù)后相關(guān)的13個(gè)基因。

用于構(gòu)建預(yù)后模型的13個(gè)基因包括NMNAT2、NT5E、QPRT、UGP2、MTHFD2、TSTA3、TYMP、ALDH2、ALDH1A1、IDO1、IL4I1、INPP1和ENPP6。NMNAT2對(duì)合成煙酰胺腺嘌呤二核苷酸(NAD+)、調(diào)節(jié)細(xì)胞代謝具有重要的作用[7]。NMNAT2前體缺失的神經(jīng)元減弱了小鼠他莫昔芬對(duì)葡萄糖代謝和胰島素耐受的影響，最后證實(shí)NMNAT2參與脂質(zhì)和糖代謝[8]。NT5E是一種糖基磷脂酰肌醇錨定的細(xì)胞表面蛋白，在細(xì)胞外單磷酸腺苷(AMP)分解代謝為腺苷的過(guò)程中起著關(guān)鍵作用，有研究表明NT5E基因甲基化與乳腺癌的發(fā)生有關(guān)，并與乳腺癌的不良預(yù)后因素相關(guān)[9]。NT5E表達(dá)變化和相關(guān)代謝組學(xué)變化也與腫瘤對(duì)包括鉑類藥物在內(nèi)的幾種化療藥物的敏感性相關(guān)[10]。QPRT催化色氨酸轉(zhuǎn)移生成煙酸單核苷酸，這是關(guān)鍵輔酶煙酰胺腺嘌呤二核苷酸(NAD)從頭生物合成的前體，對(duì)細(xì)胞呼吸鏈中的氧化磷酸化至關(guān)重要[11]。乳腺癌中DSCAM-AS1基因敲除降低了QPRT表達(dá)，從而抑制乳腺癌的生長(zhǎng)、遷移和侵襲[12]。UGP2是一種由508個(gè)氨基酸殘基組成的酶，在糖原生物合成中起著至關(guān)重要的作用[13]。MTHFD2是一種線粒體單碳代謝酶，已證明和包括乳腺癌在內(nèi)的多種惡性腫瘤發(fā)展密切相關(guān)[14]。TSTA3巖藻糖合成酶在糖代謝的過(guò)程中發(fā)揮了重要的作用[15]。TSTA3在乳腺癌中是一種致癌基因，miR-125a-5p/miR-125b抑制乳腺癌細(xì)胞中TSTA3的表達(dá)，TSTA3通過(guò)調(diào)節(jié)CXCR4的表達(dá)來(lái)控制細(xì)胞增殖和侵襲[16]。TYMP是胸苷酸合成酶，通過(guò)促進(jìn)胸苷酸合成而調(diào)節(jié)代謝過(guò)程。在乳腺癌的研究中發(fā)現(xiàn)TYMP高表達(dá)和化療藥卡培他濱反應(yīng)顯著相關(guān)[17]。ALDH2是一種乙醛脫氫酶，ALDH2多態(tài)性對(duì)乙醛代謝產(chǎn)生影響，具有ALDH2的Lys/Lys基因型的亞裔女性中，乳腺癌患病風(fēng)險(xiǎn)增加[18]。ALDH1A1可促進(jìn)維生素A和視黃醇合成全反式維甲酸，ALDH1A1敲除的小鼠對(duì)飲食誘導(dǎo)的肥胖具有抵抗力[19]。ALDH1A1降低乳腺癌細(xì)胞內(nèi)pH值，促進(jìn)TAK1磷酸化，激活NFκB信號(hào)促進(jìn)了腫瘤的生長(zhǎng)[20]。IDO1是一種含血紅素的酶，催化色氨酸分解為犬尿氨酸。IDO1已證明是女性癌癥中重要一個(gè)免疫相關(guān)基因，參與了抗腫瘤免疫過(guò)程，并與突變負(fù)荷相關(guān)[21]。IL4I1催化苯丙氨酸氧化脫氨基生成苯丙酮酸，同時(shí)生成H2O2和NH3。IL4I1通過(guò)重塑乳腺癌細(xì)胞中的糖代謝來(lái)加強(qiáng)免疫治療的潛力[22]。INPP1是一種使游離多磷酸肌醇去磷酸化的酶，參與磷脂酰肌醇信號(hào)通路的6種關(guān)鍵分子合成[23]。ENPP6對(duì)含有膽堿的化合物(如GPC22和O-磷酰膽堿N-?；掖及?3)具有磷酸二酯酶活性，以產(chǎn)生磷酸膽堿[24]。

進(jìn)一步計(jì)算每個(gè)基因的風(fēng)險(xiǎn)系數(shù)，展示每個(gè)基因在樣本中的表達(dá)情況。根據(jù)基因表達(dá)量將樣本分為高低風(fēng)險(xiǎn)兩組，展示不同風(fēng)險(xiǎn)分組樣本的生存狀態(tài)。隨著風(fēng)險(xiǎn)得分的增加，相同時(shí)間節(jié)點(diǎn)，患者死亡率增高。最后對(duì)高低風(fēng)險(xiǎn)兩組進(jìn)行生存比較，發(fā)現(xiàn)高風(fēng)險(xiǎn)組患者的生存率更低。對(duì)風(fēng)險(xiǎn)得分進(jìn)行了單因素和多因素的Cox回歸分析，結(jié)果表明，利用代謝相關(guān)基因構(gòu)建的風(fēng)險(xiǎn)的得分能夠作為患者的一個(gè)獨(dú)立危險(xiǎn)因素去預(yù)測(cè)患者的預(yù)后，利用ROC曲線下面積檢驗(yàn)了風(fēng)險(xiǎn)得分預(yù)測(cè)患者預(yù)后的結(jié)果是較為準(zhǔn)確的。根據(jù)風(fēng)險(xiǎn)得分進(jìn)行GSEA分析，發(fā)現(xiàn)風(fēng)險(xiǎn)得分相關(guān)基因顯著富集到代謝相關(guān)的通路。

綜上所述，通過(guò)篩選得到的代謝相關(guān)基因顯著影響著乳腺癌患者的代謝情況，進(jìn)而對(duì)患者的預(yù)后產(chǎn)生影響。通過(guò)單因素和多因素Cox回歸分析，得到影響乳腺癌患者預(yù)后的代謝相關(guān)基因，富集分析表明這些基因通過(guò)多種通路影響乳腺癌患者的預(yù)后情況，這為研究乳腺癌患者代謝與生存的關(guān)系提供了證據(jù)。研究構(gòu)建的預(yù)后模型可以幫助臨床醫(yī)師對(duì)患者的預(yù)后提供有益的參考。