基于免疫細胞浸潤評分實現膀胱癌分型及預后風險評估

殷桂草，鄭生旗，張 偉，董 欣，祁樂中，李一帆

1. 揚州大學附屬醫院泌尿外科，江蘇 揚州 225000

2. 揚州大學護理學院 公共衛生學院，江蘇 揚州 225000

膀胱癌具有高度異質性，是世界上第五大常見惡性腫瘤［1］。2018年，全世界超過19.9 萬人死于膀胱癌，新發膀胱癌病例超過54.9 萬例［2］。約3/4 的膀胱癌在最初診斷時為非肌層浸潤性膀胱癌［3］，可以通過經尿道膀胱腫瘤切除聯合膀胱內灌注化療進行治療，預后一般較好；一旦病變累及膀胱肌層，或出現遠處轉移，患者的5年總存活率將從80%下降到15%以下［4-5］。因此，迫切需要開發新的預后評估策略和有效的治療方法來改善膀胱癌患者的預后。

近年來，免疫療法顯示出積極的抗腫瘤作用［6-7］。研究證實，PD-1/PD-L1 檢查點抑制劑對肌層侵襲性和轉移性膀胱癌具有較好的療效［8］。腫瘤微環境中的免疫成分在腫瘤細胞基因表達和臨床預后中發揮重要作用［9］，不同免疫細胞浸潤豐度的膀胱癌患者在基因組、轉錄組和生物學過程上存在顯著異質性［10］。因此，量化腫瘤中關鍵免疫細胞成分可以為預測腫瘤患者的免疫治療效果及預后評估提供新的方法。

隨著高通量技術的成熟和腫瘤免疫研究算法的發展，轉錄組測序可以用來描述腫瘤免疫微環境［11］。然而，目前還沒有研究系統地探究免疫浸潤評分與膀胱癌患者治療及預后的關系。本研究基于ssGSEA 算法獲得16種免疫細胞浸潤的評分，通過無監督聚類將患者分為高、低免疫浸潤兩個聚類，研究不同聚類之間免疫檢查點相關基因的表達、化療藥物敏感性和免疫治療效果的差異；通過WGCNA 發現與關鍵免疫細胞相關的核心基因，Cox 回歸篩選預后相關基因，構建預后風險評分模型，并利用多因素分析發現該風險評分是膀胱癌患者獨立的預后因素；最后，根據風險評分和臨床病理參數建立列線圖來預測患者1、3、5年的總存活率，為患者的預后評估及個性化治療提供參考。

1 資料與方法

1.1 數據收集及處理

從TCGA 數據庫（https：//portal.gdc.cancer.gov/repository）中獲取19 份正常組織標本和411份膀胱癌組織標本的RNA 測序數據及相關臨床數據。排除重復標本和患者總存活時間不足30 d 的標本。從GEO 數據庫（https：//www.ncbi.nlm.nih.gov/gds/）下載GSE13507 中標本的mRNA表達數據和患者的臨床數據。從TCIA 數據庫（https：//tcia.at/home）下載膀胱癌患者對不同免疫治療方案的免疫細胞陽性比例分數評分數據，免疫細胞陽性比例分數得分越高代表患者免疫治療的效果可能越好。

1.2 采用ssGSEA分析免疫細胞浸潤情況

利用“GSVA”包進行ssGSEA，計算正常組織和腫瘤組織16 種免疫細胞浸潤的評分，并使用R中的“limma”、“ggplot2”和“ggpubr”包對正常組織和腫瘤組織免疫浸潤評分的差異進行比較。

1.3 采用無監督聚類進行膀胱癌患者分型分析

為了進一步探討免疫浸潤評分對膀胱癌患者預后和治療的價值，使用“ConsensusClusterPlus”包基于16 種免疫細胞浸潤評分對膀胱癌患者進行聚類分型，設置maxK=9，reps=100，pItem=0.8，pFeature=1，distance=“manhattan”，clusterAlg=“pam”。同時，使用ggplot2 R 包生成無監督聚類和相應的代表性數據。之后，對不同聚類的患者進行預后分析，并進行主成分分析以驗證分型的效果。使用“Heatmap”包繪制不同分型下免疫細胞浸潤評分的熱圖。使用“limma”包、“ggplot2”包評估不同分型下免疫檢查點相關基因的表達及患者對免疫治療反應的差異。最后，使用“prophytic”包評估不同分型的膀胱癌患者對臨床藥物的敏感性。

1.4 采用WGCNA 篩選與關鍵免疫細胞顯著相關的基因集

WGCNA 的輸入數據文件為：①TCGA 數據集中411 個膀胱癌腫瘤組織標本的轉錄組數據；②前述鑒定的關鍵免疫細胞的免疫浸潤評分數據。首先，通過中位數絕對偏差測量411 個標本中基因表達的可變性，并確定在正常組織和腫瘤組織中差異表達的基因。構建帶表型熱圖的樣本樹突圖。然后，通過軟閾值函數計算出最佳軟閾值，并通過軟閾值構造加權鄰接矩陣。鄰接矩陣進一步轉化為拓撲重疊矩陣。然后，將具有相似表達的基因聚類成一個模塊，稱為共表達模塊。計算P值和相關系數以確定共表達模塊與表型之間的關聯。最后，鑒定出與關鍵免疫細胞浸潤顯著相關的關鍵模塊，并提取模塊內的關鍵基因用于下一步研究。

1.5 功能富集分析研究關鍵基因的潛在功能

使用“ClusterProfiler”包進行富集分析，以探索關鍵基因潛在的功能，并選擇GO 和KEGG 作為背景數據庫。

1.6 構建基于關鍵基因的預測模型并驗證

首先，采用“glmnet”包對1.4 篩選出來的基因進行LASSO Cox 回歸分析，進一步確定與膀胱癌患者預后相關的基因。然后，從HPA 數據庫下載膀胱癌預后相關基因在膀胱癌組織和正常膀胱組織中的免疫組織化學結果，驗證其在組織中的蛋白質表達水平。針對無免疫組織化學結果的HOXB5和HOXB6采用qRT-PCR在臨床標本中驗證。具體方法：收集揚州大學附屬醫院8 份膀胱癌組織標本和鄰近正常組織標本，其中正常組織標本距離腫瘤組織2 cm。所有操作規程均通過揚州大學附屬醫院倫理委員會審查（2020-YKL03-G042）。所有標本均存儲在－80 ℃冰箱保存。使用TRIzol 試劑（美國Invitrogen 公司）從樣品中提取總RNA，然后用RNeasyMaxi 試劑盒（德國Qiagen 公司）純化獲得RNA。完成上述步驟后，按照miScriptⅡRT試劑盒（德國Qiagen 公司）的說明書，逆轉錄合成cDNA。使用miScriptSYBR?GreenPCR 試劑盒（德國Qiagen公司）進行實驗，反應所需儀器為LC480熒光定量PCR 儀。反應條件：95 ℃ 30 s，95 ℃10 s，55 ℃ 30 s，72 ℃ 30 s，共45 個循環。相關基因的引物序列見表1。然后，基于篩選并經驗證的膀胱癌預后相關基因構建風險評分模型：風險評分=（Coef1×mRNA1的表達量）+（Coef2×mRNA2的表達量）+…+（Coefn×mRNAn 的表達量）。其中，Coef 為相應mRNA 的LASSO Cox 回歸模型系數。通過中位風險評分將膀胱癌患者分為高、低風險組。使用R 中的“Kaplan-Meier”包評價兩組患者的預后，使用R 中的“survivalROC”包生成ROC 曲線，并使用“timeROC”包計算模型的AUC值。同時，從GEO 數據庫中下載GSE13507 數據集（包含165 例具有完整臨床隨訪數據的膀胱癌患者），對模型進行外部驗證。

表1 qRT-PCR引物設計Table 1 qRT-PCR primers

1.7 構建預測膀胱癌患者預后的列線圖

首先，結合風險評分和臨床病理特征進行單因素和多因素Cox 回歸分析，以探索風險評分模型是否為膀胱癌患者預后的獨立相關因素。然后，基于風險評分和獨立的臨床病理特征構建列線圖，并使用C 指數和校準曲線驗證列線圖的預測效能。

1.8 統計學方法

使用R 4.0.3 和GraphPad 9.0 進行統計分析。連續性變量用均數±標準差（±s）描述，采用Student’st檢驗或Wilcoxon秩和檢驗進行比較；分類變量用頻率（%）描述，采用χ2檢驗或Fisher 精確檢驗。繪制Kaplan-Meier 生存曲線和時間依賴性ROC 曲線，評估風險評分模型的預測價值。采用單因素和多因素Cox 回歸分析確定膀胱癌患者預后的獨立影響因素。P<0.05 為具有統計學意義。

2 結 果

2.1 正常組織與腫瘤組織免疫細胞浸潤評分差異及與患者預后的相關性

與正常組織標本比較，膀胱癌組織標本中B 細胞、肥大細胞、中性粒細胞、輔助性T 細胞和腫瘤浸潤淋巴細胞的浸潤評分明顯升高（均P<0.05），提示這五種免疫細胞可能在膀胱癌免疫微環境的改變中發揮重要作用，是膀胱癌的關鍵免疫細胞。見圖1。

圖1 癌癥基因組圖譜膀胱癌隊列中正常組織與腫瘤組織16種免疫細胞浸潤豐度比較Figure 1 Comparison of 16 kinds of immune cell infiltration abundance between normal tissues and tumor tissues in bladder cancer cohort from the Cancer Genome Atlas

聚類分析結果顯示，411 份膀胱癌患者標本可初步分為三類（圖2）。生存曲線分析發現，聚類2 中患者的預后顯著優于聚類1 和聚類3（均P<0.05），而聚類1和聚類3患者的預后無顯著差異（P>0.05），見圖3A。基于此，本研究將聚類1 和聚類3 合并為一類（聚類1′）重新進行生存分析。結果顯示，聚類1′和聚類2 中患者的存活率差異具有統計學意義（圖3B）。進一步行主成分分析發現，當TCGA 樣本被劃分為兩個簇時，樣本可以明顯分離（附圖1），且聚類2 患者的免疫細胞浸潤評分明顯高于聚類1′患者（附圖2）。

2.2 不同聚類患者治療情況比較

比較不同聚類之間免疫檢查點相關基因的表達發現，在聚類2 中免疫檢查點相關基因的表達均顯著升高（附圖3A）。TCIA 數據庫分析發現，聚類2 中患者對CTLA4 抑制劑、PD-1 抑制劑及兩者聯合更敏感，更能夠從免疫治療中獲益（附圖3B）。分析兩種分型患者對常用化療藥物的敏感性發現，與聚類1′比較，聚類2患者對恩貝酸、多西他賽、環巴胺和阿卡地新更為敏感（附圖3C）。

2.3 關鍵免疫細胞相關基因組篩選結果

從附圖4可以看出，當軟域值為5時，數據更符合冪律分布，平均連通性趨于穩定，數據適合進一步研究。確定了軟閾值后，根據基因間表達量的相關性構建基因聚類樹。將最小模塊數設置為10，深度靈敏度設置為 2，使用動態剪切法識別模塊。最后，將相似度小于0.8的模塊合并，共得到25 個非灰色模塊（圖4）。在正相關性模塊中，棕色模塊中除肥大細胞外，其余相關系數均大于0.7，可見其與免疫細胞的相關性系數相對較大；在負相關性模塊中，僅有黑灰色模塊，見表2。本研究后續選擇棕色模塊和黑灰色模塊，匯總兩個模塊共得到35 個基因。相關基因的表達見附圖5。

圖4 基于最優軟閾值構建基因聚類樹及共表達網絡Figure 4 Construction of gene clustering tree and coexpression network based on optimal soft threshold

圖5 風險評分模型在TCGA數據集中的驗證Figure 5 Validation of the risk scoring model in the TCGA dataset

表2 不同模塊與關鍵免疫細胞的相關性分析Table 2 Correlation analysis between different modules and key immune cells

2.4 關鍵基因富集分析結果

對35 個關鍵基因進行GO 分析結果顯示，在生物過程及分子功能分析中，35 個關鍵基因顯著參與單核細胞和淋巴細胞等免疫細胞的分化以及正向調節細胞因子的產生（附圖6）；在細胞成分分析中，相關基因主要富集于質膜的外側。KEGG 富集分析結果顯示，這些基因主要參與T細胞受體信號通路、細胞黏附分子及細胞因子-受體相互作用（附圖7）。

圖6 風險評分模型在GSE13507數據集中的驗證Figure 6 Validation of the risk scoring model in GSE13507 dataset

圖7 基于風險評分構建膀胱癌患者預后相關列線圖Figure 7 Prognostic nomogram of bladder cancer patients based on risk score

2.5 膀胱癌預后風險評分模型的構建與驗證

基于之前篩選出來的35個基因，進一步通過LASSO Cox 回歸篩選出4 個與患者預后相關的基因GPR171、HOXB3、HOXB5和HOXB6。通過HPA數據庫比較相關基因的蛋白表達水平發現，GPR171 在腫瘤組織中高表達，HOXB3 在腫瘤組織中低表達（附圖8）；通過臨床組織標本進行驗證發現，HOXB5和HOXB6mRNA 在腫瘤組織中均高表達（分別是正常組織的3.75 和2.87 倍，均P<0.01）?；贕PR171、HOXB3、HOXB5和HOXB6構建膀胱癌患者的預后風險評分模型：風險評分=（－0.3654×GPR171 表達水平）+（－0.2885×HOXB3表達水平）+（－0.2885×HOXB5 表達水平）+（0.3087×HOXB6表達水平）。采用上述公式計算樣本的風險評分范圍為－2.41~－0.03，評分中位數為－0.83，根據中位風險評分將膀胱癌患者分為高風險組和低風險組，四個基因在高低風險組的表達見附圖9。

圖8 列線圖模型用于預測患者1、3、5年總存活率的校準圖Figure 8 Calibration plots of nomogram models for predicting 1， 3， or 5 year overall survival rates

Kaplan-Meier 生存曲線提示，高風險組患者預后較差，風險評分與預后呈負相關（圖5A）。ROC 曲線分析結果顯示，風險評分預測患者1、3、5年總生存時間的曲線下面積分別為0.735、0.765、0.799，均超過0.7（圖5B）。在外部數據集GSE13507 中驗證評分模型的預測價值，結果與前述一致，低風險組預后好（圖6A）；隨時間變化的ROC 曲線顯示，1年的AUC 為0.667，3年的AUC 為0.613，5年時AUC 為0.629（圖6B）。提示基于關鍵免疫細胞相關基因構建的預后風險評分模型對于膀胱癌患者預后預測的準確性較好。

2.6 膀胱癌預后評估列線圖的構建和驗證

單因素Cox 回歸分析結果顯示，年齡、T 分期、病理分期和風險評分與患者的預后顯著相關；多因素分析顯示，風險評分、T 分期和病理分期是膀胱癌患者預后的獨立危險因素。見表3。結合風險評分和臨床參數構建列線圖，C 指數為0.72（圖7）；校準曲線顯示，模型的預測結果與實際結果基本重合（圖8），提示列線圖對膀胱癌患者1、3、5年總存活率具有較好的預測效能。

表3 膀胱癌患者預后相關危險因素分析結果Table 3 Analysis of prognostic risk factors in patients with bladder cancer

3 討 論

膀胱癌是世界范圍內發病率較高且最致命的惡性腫瘤之一，腫瘤的高復發風險使得膀胱癌患者的預后較差，給社會和患者造成了沉重的經濟負擔［2］。目前，手術結合膀胱內灌注治療是非肌層浸潤性膀胱癌的主要治療手段［12］。隨著對膀胱癌的免疫病理機制和腫瘤免疫微環境的進一步了解，基于腫瘤免疫微環境調控的免疫治療已成為膀胱癌患者新的治療選擇［13］。腫瘤免疫微環境由許多具有抗腫瘤或促腫瘤活性的不同免疫細胞亞群組成。既往研究發現，巨噬細胞作為先天免疫效應細胞，在免疫治療中起著至關重要的作用，與免疫檢查點阻斷應答有關［14-15］。Jiang等［16］通過巨噬細胞相關基因構建并驗證了預測膀胱癌患者預后的模型。除了巨噬細胞之外，CD8+T細胞、NK細胞和B細胞等多種免疫細胞參與機體免疫。這些細胞通過不同途徑促進或抑制機體免疫微環境的形成，各個細胞間存在協同或抑制效應，從而在腫瘤免疫中發揮作用。因此，系統研究免疫細胞在膀胱癌的作用具有重要意義。

隨著免疫檢查點抑制劑在臨床中的運用，肌層浸潤性膀胱癌患者能夠獲得更好的預后。PD-1（PDCD1）與腫瘤組織產生的PD-L1（CD274）結合，導致有限的宿主免疫應答。PD-L1 抑制劑通過增加浸潤的CD8+T 細胞水平，顯示出有效的抗腫瘤免疫應答［17］。T 淋巴細胞在癌癥的“免疫監視”中發揮著重要作用，不僅通過表達多態性抗原受體以識別特異性抗原，還擁有效應功能以及記憶特征［18-19］。因此，本研究以16 種免疫細胞為基礎，通過ssGSEA 算法量化了TGCA 膀胱癌隊列中免疫細胞的豐度，系統地探討了免疫細胞評分與膀胱癌患者治療及預后的關聯。通過無監督聚類將膀胱癌基于免疫細胞評分分為兩個亞型，其中聚類2 患者的預后明顯優于聚類1′患者。有報道稱，PD-L1 的高表達水平與膀胱癌腫瘤組織的惡性程度及不良預后顯著相關，且此類患者術后復發率較高［15］。本文資料顯示，PDCD1、CD28和CTLA4等重要免疫檢查點相關基因在聚類2 顯著上調，而免疫檢查點相關基因的表達與腫瘤免疫微環境顯著相關，低表達的免疫檢查點相關基因可能促進腫瘤免疫逃逸［11］。阿特珠單抗是目前唯一批準用于膀胱癌治療的PD-L1 抑制劑，但只有少數患者從免疫治療中受益［20］。本研究開發的免疫評分聚類分型方法根據膀胱癌患者對PD-1 或CTLA4 抑制劑的敏感性進行分層發現，聚類2 對PD-1 抑制劑和CTLA4 抑制劑聯合治療反應更強，表明聚類2 是免疫有利亞群，聚類2 聚類中的患者更可能從阿特珠單抗治療中獲益。此外，本研究還發現聚類2 對于部分化療藥物的半抑制濃度小于聚類1′，提示聚類2 的膀胱癌患者更容易從這些化療藥物中獲得臨床療效，避免更高劑量的化療藥物所致各種不良反應。

為了更好地將這些相關免疫浸潤評分應用于臨床，本研究通過WGCNA 篩選出與關鍵免疫細胞顯著相關的基因，GO 富集分析顯示這些基因主要參與免疫細胞的分化與細胞因子產生的正向調節。KEGG 富集分析表明這些基因主要參與T 細胞受體信號通路。通過LASSO Cox 回歸進一步篩選出四個與膀胱癌患者預后相關的基因。其中，GPR171 是一種G 蛋白偶聯受體，由神經內分泌途徑蛋白衍生的BigLEN 是GPR171 的內源性配體。有研究顯示，GPR171 可調節攝食和焦慮行為［21］，并在多種腫瘤中被證實是T 淋巴細胞的特征基因［22-23］。研究發現，GPR171在腫瘤浸潤淋巴細胞中上調，且其在黑色素瘤患者中的表達可以預測患者接受免疫治療的反應［24］。GPR171 在T淋巴細胞中的表達動力學與其他免疫檢查點（如PD1、TIGIT、CD112R 和TIM36）非常相似，在T 淋巴細胞中GPR171 可上調并通過信號轉導抑制T淋巴細胞介導的免疫反應。此外，GPR171 信號的破壞可促進T 淋巴細胞介導的抗腫瘤免疫［25］。HOXB3、HOXB5 和HOXB6 都屬于HOXB 家族。HOXB基因家族是一個長度約為180 bp的同源片段，編碼60 個氨基酸的同源結構域，與轉錄調控相關。HOXB 基因家族成員在腫瘤細胞增殖、遷移和凋亡以及腫瘤的形成和發展中均起著至關重要的作用［26］。已有研究表明，多個HOXB 基因家族成員在膀胱癌、膽管癌、子宮內膜癌和乳腺癌種差異表達［27］，如ISL LIM 同源盒蛋白1 和LIM 同源盒蛋白5 在膀胱腫瘤中表達上調，并與腫瘤的發展及轉移密切相關［28］。HOXB3是一種蛋白編碼基因，與急性髓系白血病的一個生物學亞群顯著相關。乳腺癌和肺腺癌患者HOXB3 的表達與患者的預后相關。乳腺癌患者的HOXB3表達水平較低，往往預后較差［29］，而在肺腺癌中高表達HOXB3 的患者預后不良。研究發現，HOXB3 與惡性腫瘤的增殖、侵襲和轉移顯著相關［30］。在膀胱癌中，HOXB3與吉西他濱的耐藥性相關，并可作為膀胱癌患者預后的生物標志物［29］。HOXB5 能夠與特定的DNA 側結合，從而影響器官發育和腫瘤發生中的關鍵信號通路［31］。據報道，HOXB5在腫瘤中表達失調，如HOXB5在視網膜母細胞瘤細胞表達上調，并通過增加基質金屬蛋白酶的產生促進腫瘤細胞的遷移和侵襲［32］。此外，HOXB5還具有通過激活Wnt/β-連環蛋白途徑促進非小細胞肺癌惡性進展的能力［31］。在乳腺癌中，Zhang 等［33］報道HOXB5 的過表達增強了乳腺癌細胞的侵襲能力。研究證實，HOXB5在膀胱癌組織和細胞系中過度表達，并促進膀胱癌細胞的增殖和遷移［34］。HOXB6 在結直腸腫瘤中差異表達并可作為患者預后的生物標志物［35］。研究表明，敲低HOXB6的表達能夠抑制肝細胞癌細胞的遷移、侵襲和增殖［36］。此外，Huang等［37］報道了Hsa_circ_0007031 通過海綿miR-196a-5p 調節HOXB6 來促進骨肉瘤細胞的增殖和遷移。目前關于HOXB6 與膀胱癌的研究還甚少，需要進一步研究其機制。

本研究基于上述關鍵基因構建了風險評分模型，并通過外部驗證證實了該風險評分模型的預測效能。結合臨床特征構建的列線圖也展示出對膀胱癌患者預后較好的預測價值。與其他基因集構建的膀胱癌預后風險評估模型比較，本研究構建的模型在預測膀胱癌患者5年存活率上有著更高的效能［38-39］。此外，相比Liang等［40］構建的15個基因的模型，本研究構建的模型納入的基因數較少，檢驗難度小，便于臨床應用。

綜上所述，本文系統研究了膀胱癌患者免疫細胞浸潤狀態，并基于關鍵免疫細胞相關基因構建了預后風險評分模型，對膀胱癌患者的個體化治療具有重要意義。但是，由于樣本量有限，本研究中構建的模型還需通過前瞻性、多中心、更大樣本的研究進一步驗證，并通過開展相關功能體內外實驗探索關鍵免疫細胞和核心基因在膀胱癌中發揮作用的相關機制。

本文附圖見電子版。

志謝研究得到國家自然科學基金（82002675）、江蘇省科技計劃青年基金（BK2020938）、揚州市重點研發計劃社會發展項目（YZ2020110）、揚州市軟科學研究計劃（YZ2022267）、江蘇省博士后研究資助計劃（2020Z268）、揚州大學高層次人才研究啟動基金（2019LYF）支持

AcknowledgementsThis work was supported by the National Natural Science Foundation of China (82002675),Science and Technology Plan Youth Fund Project of Jiangsu Province (BK2020938), Key R&D Plan-Social Development Project of Yangzhou City (YZ2020110), Soft Science Research Plan of Yangzhou City (YZ2022267), Postdoctoral Research Funding Program of Jiangsu Province (2020Z268),and Yangzhou University High-level Talent Research Startup Fund (2019LYF)

利益沖突所有作者均聲明不存在利益沖突

Conflict of InterestsThe authors declare that there is no conflict of interests

?The author(s) 2024. This is an open access article under the CC BY-NC-ND 4.0 License (https://creativecommons.org/licenses/by-nc-nd/4.0/)