綜合生物信息學分析揭示乳腺癌細胞衰老的預后及免疫特征

王科江 李雪森 （西南醫科大學基礎醫學院腫瘤醫學研究所，瀘州 646000）

乳腺癌是全球范圍內高度流行的癌癥。過去二十年，早期診斷工具和治療方法進步使乳腺癌病死率與1990年相比降低了1/3［1-2］。但乳腺癌仍是女性患者癌癥死亡的主要原因，威脅全球婦女的健康。越來越多的證據表明乳腺癌具有高度異質性，局部腫瘤微環境對腫瘤發展至關重要，包括乳腺癌起始、進展、轉移及耐藥［3］。

細胞衰老是衰老的關鍵過程之一，是衰老與癌癥間的紐帶［4］。衰老與腫瘤的聯系復雜，目前知之甚少。細胞衰老在腫瘤發生和發展過程中具有“雙刃劍”作用。一方面，衰老細胞進入永久性細胞周期停滯的背景下，衰老確保組織動態平衡并防止腫瘤發生［5-6］。另一方面，衰老細胞通過促進多種類型免疫抑制細胞積累及各種炎癥相關信號分子和細胞因子激活改變腫瘤微環境，廣泛影響腫瘤微環境和腫瘤生長［7-8］。值得注意的是，衰老和細胞衰老可改變免疫細胞健康狀況，并最終影響癌癥治療效果，尤其是免疫治療。但細胞衰老與腫瘤微環境的關系尚不清楚，細胞衰老相關基因在評價乳腺癌免疫浸潤和臨床預后方面的價值有待進一步研究。

為系統評估乳腺癌患者細胞衰老與預后的相關性，本研究建立了基于細胞衰老相關基因的新風險模型，并探討了其作為預測預后生物標志物的潛在重要性，基于細胞衰老相關特征深入分析了風險亞群、免疫檢查點和免疫細胞滲透間的關系，為乳腺癌免疫治療策略及細胞衰老調控機制提供了新的見解。

1 資料與方法

1.1 數據來源與下載 從TCGA網站（https：//portal.gdc.cancer.gov/）下載乳腺癌RNA-seq轉錄組數據和相關臨床信息，對mRNA count數據進行標準化，將標準化的count數據轉化為TPM，進一步對TPM數據進行log2對數轉化用于后續分析。將表達數據與臨床數據相匹配，去除重復樣本，共113例正常樣本和1 050例腫瘤樣本納入分析（表1），從GEO數據庫（https：//www.ncbi.nlm.nih.gov/geo/）下載并整理具有預后信息的327例乳腺癌患者隊列（GSE20685）作為外部驗證集。細胞衰老相關基因來自CellAge數據庫（https：//genomics.senescence.info/cells/），該數據庫收錄了278個經實驗驗證的細胞衰老相關基因。

表1 TCGA-BRCA隊列臨床信息基線資料Tab.1 Baseline clinical information of TCGA-BRCA cohort

1.2 差異顯著的細胞衰老相關基因鑒定及富集分析 R軟件“DEseq2”包鑒別差異表達的細胞衰老相關基因［9］，cut-off值設為|log2FC|＞1，Padj＜0.05，利用R軟件“clusterProfiler”包對差異表達的細胞衰老相關基因進行GO和KEGG通路富集分析［10］。

1.3 細胞衰老相關基因預后風險模型構建及外部數據集驗證 TCGA中1 050個乳腺癌患者樣本作為訓練集，采用單因素Cox回歸分析評估TCGA中差異表達的細胞衰老相關基因與患者總生存期（OS）的關系，篩選P＜0.05的所有基因進行進一步分析。采用LASSO分析選擇OS相關的最佳預后基因。通過多因素Cox回歸分析得到最優模型，計算風險評分，風險評分=基因表達量1×Coef1+基因表達量2×Coef2+……+基因表達量n×Coefn（Coef：基因在多因素Cox回歸分析中的回歸系數，n：預后相關細胞衰老相關基因總數）。根據風險評分中位值將患者分為高、低風險組。將GSE20685的327個乳腺癌患者隊列作為外部驗證集，將基因表達數據帶入TCGA訓練集風險評分公式計算各樣本風險評分進行驗證。利用“survival”“survminer”“pheatmap”包繪制TCGA訓練集和GSE20685外部驗證集風險評分曲線圖、風險評分散點圖、生存曲線圖、高低風險評分組細胞衰老相關基因表達熱圖。

1.4 預后模型性能評估 在TCGA隊列中對風險評分和相關臨床參數進行單因素和多因素Cox回歸分析，確定風險評分和相關臨床參數是否可作為乳腺癌患者OS的獨立預測因子。納入臨床病理參數和風險評分利用R軟件“regplot”包生成預后列線圖預測TCGA中乳腺癌患者一年、三年和五年OS，并采用校準曲線驗證預測的生存概率與實際觀察結果的吻合度。采用ROC曲線評價風險評分及相關臨床參數預測率。

1.5 GSEA基因集富集分析 GSEA可通過將預定義的基因集與特定表型進行比對找出預定義基因集相關表型。根據細胞衰老關鍵預后基因風險評分中位值進行分組，并與C2.CP.KEGG.v7.2基因集、Hallmark Gene Sets（癌癥特征基因集合）基因集進行比對，C2.CP.KEGG.v7.2基因集和Hallmark收集自分子簽名數據庫（MSigDB），利用R軟件“clusterProfiler”包進行分析［10］，探索細胞衰老關鍵預后基因在乳腺癌中可能的作用機制。顯著富集的基因集符合以下標準：P＜0.05，FDRq＜0.25。

1.6 免疫微環境及免疫細胞浸潤分析 利用“CIBERSORT”包反卷積算法評估所有TCGA乳腺癌病例中免疫細胞浸潤比例，并繪制免疫細胞浸潤比例豐度圖和22種免疫細胞相關性熱圖［11］。根據細胞衰老關鍵預后基因風險評分中位值進行分組，評估22種免疫細胞在高、低風險組的浸潤比例和免疫檢查點基因表達。利用“ESTIMATAE”包進行免疫微環境分析，計算免疫評分，評估風險評分與免疫微環境的關系［12］。

1.7 體細胞突變分析 體細胞突變數據來源于TCGA（https：//portal.gdc.cancer.gov/），數據類型為“Masked_Somatic_Mutation”，利用“maftools”包分析TCGA體細胞突變數據［13］，計算基因突變頻率。

1.8 統計學分析 所有統計分析和數據繪圖均采用R軟件（4.1.2）進行。

2 結果

2.1 TCGA數據庫中差異表達的細胞衰老相關基因鑒定及差異表達基因功能富集分析 278個細胞衰老相關基因中有68個基因在腫瘤組織和正常樣本中的表達存在顯著差異，其中36個基因表達上調，32個基因表達下調（圖1A、B）。通過功能富集分析探討上述68個細胞衰老差異表達基因在TCGA-BRCA隊列中的生物學功能和顯著富集通路，GO及KEGG富集分析結果顯示，差異基因顯著富集于細胞周期、細胞衰老相關信號通路及生物學進程（圖1C、D）。

圖1 細胞衰老相關基因差異分析及GO/KEGG富集分析Fig.1 Difference analysis of senescence related genes and GO/KEGG enrichment analysis

2.2 細胞衰老相關基因預后模型構建及驗證 在TCGA隊列中對上述68個細胞衰老差異基因進行單因素Cox回歸分析，結果顯示8個基因與乳腺癌患者OS相關（P＜0.05，圖2A）。為避免過度擬合，對8個有預后價值的基因進行LASSO回歸（圖2B、C），對LASSO回歸確定的7個預后相關基因進行多因素Cox回歸，確定了WT1、IFNG、TP63、IGFBP6、CPEB1 5個細胞衰老關鍵預后基因的最佳預后模型，風險評分=WT1×0.162 094+IFGN×（-0.320 910）+TP63×（-0.120 800）+IGFBP6×（-0.141 510）+CPEB1×0.385 455，預后模型森林圖如圖2D所示。根據中位風險評分將患者分為高、低風險組，風險評分散點圖、風險評分曲線圖顯示隨著風險評分升高，相同時間點患者死亡率越高（圖3A、B），熱圖展示了5個基因在高、低風險組的表達（圖3C），生存曲線顯示高風險組預后較差（圖3D）。為驗證預后模型，對來自GEO數據庫的GSE20685乳腺癌隊列按照與TCGA隊列相同的風險評分計算公式計算各樣本風險評分，根據TCGA隊列中位評分將GSE20685驗證隊列分為高、低風險組，結果顯示構建的模型具有普適性，在GSE20685隊列中得到了驗證（圖3E～H）。

圖2 細胞衰老相關基因預后風險模型構建Fig.2 Construction of a prognostic risk model for cell senescence related genes

圖3 TCGA訓練集和GSE20685外部驗證集高、低風險評分組患者基因表達及生存狀態評估Fig.3 Evaluation of gene expression and survival status of patients in TCGA training set and GSE20685 external validation set with high and low risk score

2.3 預后模型性能評估及列線圖建立 為進一步驗證預后模型的準確性，研究風險評分是否可作為TCGA-BRCA隊列的OS獨立預后因素，對臨床特征和風險評分進行單因素Cox回歸分析和多因素Cox回歸分析，結果顯示僅風險評分、年齡、淋巴結轉移分期（N_stage）在單因素Cox回歸分析和多因素回歸分析中均為獨立預后因素（P＜0.001，圖4A、B）。納入風險評分、年齡、淋巴結轉移分期3個獨立預后因子構建預后列線圖（圖4C），校準曲線提示該列線圖預測生存率與隊列實際生存率擬合度較好（圖4D）。繪制多重ROC曲線比較風險評分和臨床參數預測患者預后準確率，結果顯示，風險評分具有較好的預測能力（圖4E）。

圖4 預后風險評分模型的獨立預后價值Fig.4 Independent prognostic value of prognostic risk score model

2.4 GSEA富集分析 為進一步探討風險模型分類亞組間基因功能和通路差異，根據風險評分中位值將患者分為高、低風險組進行GSEA富集分析，選擇顯著富集的前10個基因集進行展示，結果顯示腫瘤免疫相關特征基因集和信號通路在高風險組中顯著下調，如KEGG數據集中的抗原加工與呈遞、趨化因子信號通路、細胞因子及細胞因子受體相互作用（圖5A），Hallmark基因集（癌癥特征基因集合）中的IL2_STAT5信號通路、IL6_JAK_STAT3信號通路，提示高風險患者可能處于免疫抑制狀態（圖5B）。

圖5 風險得分相關基因顯著富集于腫瘤免疫相關通路Fig.5 Genes related to risk score were significantly enriched in tumor immune-related pathways

2.5 高、低風險患者免疫微環境分析 基于GSEA分析結果推測高、低風險患者腫瘤免疫微環境可能存在差異。應用CIBERSORT反卷積算法進一步確認風險評分與免疫成分的關系，構建乳腺癌患者22種免疫細胞豐度圖（圖6A），分析22種免疫細胞的相關性（圖6B）。進一步分析高風險患者和低風險患者22種免疫細胞浸潤比例，結果顯示CD8 T細胞、M0巨噬細胞、M1巨噬細胞、M2巨噬細胞比例存在顯著差異（P＜0.001，圖6C）。產生活性氧的M1巨噬細胞是腫瘤生長抑制因子，而產生IL-10和TGF-β的M2巨噬細胞可促進腫瘤生長［14］。許多腫瘤中，如結直腸癌、乳腺癌、卵巢癌或膀胱癌，CD8 T細胞浸潤是很好的預后標志［15］。值得注意的是，本研究中高風險患者M2巨噬細胞顯著多于低風險患者，CD8 T細胞、M1巨噬細胞顯著少于低風險患者（P＜0.001，圖6C），表明基于風險評分對乳腺癌患者進行風險分層得到的結果與既往研究相符，風險評分可作為提示腫瘤微環境免疫細胞浸潤和預測患者預后的指標。為評估高、低風險患者對免疫治療的影響，探討風險評分與公認的免疫檢查點基因的相關性，結果顯示低風險患者PDCD1、CD274、CTLA4、CD86、LAG3、HAVCR2、TIGIT表達高于高風險患者（圖6D）；提示低風險患者免疫檢查點基因表達高于高風險患者，低風險患者可能有更好的免疫治療反應。進一步利用R軟件包ESTIMATE根據基因表達計算各腫瘤中各患者免疫評分，比較高、低風險組免疫評分差異及風險評分與免疫評分的相關性，結果顯示高風險組免疫評分顯著低于低風險患者，且風險評分與免疫評分呈負相關（圖6E、F）。提示高風險患者腫瘤微環境傾向于免疫抑制狀態，隨著風險評分升高免疫評分降低，同時高風險患者伴隨促腫瘤生長的腫瘤相關巨噬細胞高度浸潤及免疫檢查點基因表達顯著降低，可能影響高風險患者免疫治療應答，尤其是免疫檢查點抑制劑。

圖6 高、低風險患者免疫微環境分析Fig.6 Analysis of immune microenvironment in patients with high and low risk

2.6 高、低風險患者基因組突變分析 結果展示了高、低風險亞組前30個高頻突變基因，高風險組460個樣本中420個樣本發生了突變（91.3%），低風險組466個樣本中422個樣本發生了突變（90.56%），其中高風險患者抑癌基因TP53突變頻率（38.00%）高于低風險（31.00%，圖7A、B紅框部分）。為進一步說明TP53基因在兩組隊列中的突變差異，使用“Maftool”包的“mafCompare”函數比較兩組間TP53差異突變情況，該函數對兩個隊列生成的2×2關聯表進行fisher檢驗，找出差異突變基因，高風險組中TP53基因173個樣本發生了突變，低風險組中143個樣本發生了突變，OR＜1（OR=0.735）且P＜0.05（P=0.027），表明高風險患者抑癌基因TP53突變頻率更高。

圖7 高低風險亞組基因突變瀑布圖Fig.7 Waterfall map high and low risk subgroups

3 討論

乳腺癌是世界最常見的惡性腫瘤，據最近癌癥數據統計，乳腺癌占所有新診斷惡性腫瘤的30%以上，盡管乳腺癌在診斷、手術和藥物開發方面取得了一定進展，但治療反應不充分及復發和轉移，乳腺癌治療仍面臨嚴峻挑戰［16］。此外，多數治療方法均針對癌細胞開發，忽略了腫瘤微環境周圍的非癌癥成分，導致反應結果較差，腫瘤微環境包括各種非惡性細胞和非細胞成分，包括免疫系統成分、血細胞、內皮細胞、脂肪細胞和間質成分，可能以不同方式改變腫瘤行為［17］。衰老是一個復雜的生物過程，具有細胞自主性和旁分泌效應，對腫瘤微環境有重大影響［7］。但乳腺癌中衰老細胞如何與腫瘤免疫細胞相互作用及其在評估腫瘤免疫浸潤性和臨床預后方面的價值尚未見報道。

本研究分析了乳腺癌中細胞衰老相關基因的表達模式、預后價值及對腫瘤免疫微環境的影響，使用LASSO-Cox構建了新的基于TCGA數據集中5個衰老特征表達的生存預測模型，該模型在GEO數據庫的公開數據集GSE20685得到了很好的驗證。值得注意的是，本研究發現單因素及多因素Cox回歸提示風險評分是乳腺癌患者的獨立預后因素，根據風險評分中位值分組的高、低風險患者免疫狀態具有高度異質性，高風險患者伴隨促腫瘤生長的腫瘤相關巨噬細胞高度浸潤以及免疫檢查點基因表達顯著降低，具有良好預后提示的CD8 T細胞浸潤豐度降低，且高風險患者基于ESITIMATE算法得出的免疫評分顯著低于低風險患者。提示細胞衰老相關基因預測模型指導乳腺癌患者預后具有積極作用，可與特異性免疫檢查點因子或腫瘤微環境聯合作為預測免疫檢查點抑制劑反應的生物標志物。

本研究模型中共納入5個細胞衰老相關基因：WT1、IFNG、TP63、IGFBP6、CPEB1。研究發現，原癌基因KRAS誘導的腫瘤模型中，WT1的RNA干擾顯著降低了腫瘤負擔，誘導了細胞衰老，并抑制了腫瘤細胞增殖［18］。WT1還通過調節肌動蛋白活性影響細胞骨架重排，參與細胞運動，提示WT1在癌細胞侵襲和遷移中可能發揮作用［19］。IFNG屬于Ⅱ型干擾素家族成員，由免疫細胞如T細胞和NK細胞產生，通過激活效應免疫細胞和增強抗原呈遞，在抗菌、抗病毒和抗腫瘤反應中發揮關鍵作用［20］。臨床研究發現結直腸癌患者外周血單個核細胞中TNF和IFNG均受到抑制，復發性結直腸癌患者外周血單個核細胞IFNG表達降低［21］。TP63作為TP53、TP63和TP73家族基因成員之一，在轉錄調控、誘導細胞凋亡、細胞周期停滯、衰老、代謝重編程和干細胞維持方面發揮重要作用［22］。此外，TP63的兩個亞型TAp63和DNp63在干細胞維持、腫瘤轉移中發揮不同調節作用［23］。研究發現TAp63在13個腫瘤發展特征中始終表現出腫瘤抑制模式，而DNp63在腫瘤發生發展中作為腫瘤抑制基因或癌基因具有雙重作用［24］。IGFBP6是一種IGF結合蛋白，可延長IGF半衰期，并可抑制或刺激IGF在細胞培養中的促生長作用，IGFBP6可能在結腸癌發展過程中作為一種抑癌基因，低表達IGFBP6可作為結腸癌的獨立預后生物標志物［25］。另一項研究發現對TMZ（替莫唑胺）敏感的膠質瘤細胞分泌的IGFBP6驅動旁分泌機制在體內外抑制表達IGF1-R的TMZ耐藥細胞增殖，耐TMZ膠質瘤細胞中IGFBP6過表達延長了生存期［26］。CEBP1是編碼胞質多聚腺苷酸化元件結合蛋白家族成員，研究發現CPEB1是腫瘤細胞（Hela）進入M期和細胞增殖所必需的，單獨去除其家族成員CPEB4對細胞周期的影響很小，但CPEB4與CPEB1共同耗盡時可觀察到對細胞增殖的協同抑制［27］。相反的是，研究發現CPEB1表達下調促進乳腺癌細胞增殖、轉移和細胞周期進展，且CPEB1通過介導SIRT1/SOX2通路調控乳腺癌發生［28］，但CPEB1也能在腫瘤細胞中發揮促增殖作用，前提是CPEB1首先在原代細胞中誘導細胞衰老［29］，提示CPEB1在腫瘤中的調控作用可能與衰老微環境相關。進一步證實細胞衰老相關基因在腫瘤中的調控是雙向的，發揮抑癌或促癌作用的關鍵可能取決于腫瘤微環境差異。

另外，本研究雖提供了一定理論基礎和研究參考，但也存在一定局限性。利用公共數據庫的回顧性數據對現有模型進行開發和驗證，尚需更多前瞻性研究數據驗證其臨床應用，上述發現還需進一步研究證實。其次，這些細胞衰老關鍵預后基因在預后模型中調節乳腺癌進程的潛在機制尚不清楚，其生物學功能需進一步探索。最后，由于細胞衰老在癌癥中的作用很大程度上被忽視，因此更廣泛地了解癌癥、衰老和免疫環境間的聯系非常重要。本研究確定并驗證了基于5個細胞衰老相關基因的細胞衰老相關信號，作為免疫細胞浸潤指標對乳腺癌具有獨立預后價值。因此，識別影響腫瘤免疫應答的細胞衰老相關基因并進一步研究其調控機制可能有助于危險分層，并為提高乳腺癌對免疫治療的反應提供靶點。