胰腺導管腺癌免疫相關基因預后風險模型的構建與評估

張昱 任瑞平 萬鵬 賀曉蘭

摘要：目的 構建整合胰腺導管腺癌（PDAC）分子亞型和免疫相關基因的風險評估模型。方法 以GSE71729數據集（n=145）為訓練集，整合PDAC鱗狀和非鱗狀亞型之間差異表達基因和差異免疫相關基因構建調控網絡，篩選發揮主調控鱗狀亞型的5個免疫標志基因。基于患者生存信息和免疫標志基因構建整合免疫評分（IIS）模型來預測PDAC患者的臨床預后，并在5個驗證集（n=758）中檢驗其預測效能。結果 根據IIS將PDAC患者分為高危組和低危組。在訓練集和驗證集中，高危組患者的總體生存期均短于低危組（P均<0.001）。多變量Cox回歸分析顯示IIS是PDAC的獨立預后因子（HR=2.16，95%CI=1.50～3.10，P<0.001）。結論 IIS可用于PDAC患者的危險分層，并可能成為PDAC潛在的預后標志物。

關鍵詞：胰腺導管腺癌；免疫相關基因；預后

中圖分類號： R735.9? 文獻標識碼： A? 文章編號：1000-503X（2024）03-0354-07

DOI：10.3881/j.issn.1000-503X.15736

Construction and Evaluation of a Prognostic Risk Prediction Model of Pancreatic Ductal Adenocarcinoma Based on Immune-Related Genes

ZHANG Yu，REN Ruiping，WAN Peng，HE Xiaolan

Cancer Radiotherapy and Chemotherapy Center，The Affiliated Peoples Hospital of Ningbo University，Ningbo，Zhejiang 315040，China

Corresponding author：ZHANG Yu Tel：0574-87017563，E-mail：zy95589034@163.com

ABSTRACT：Objective To construct a risk prediction model by integrating the molecular subtypes of pancreatic ductal adenocarcinoma （PDAC） and immune-related genes.Methods With GSE71729 data set （n=145） as the training set，the differentially expressed genes and differential immune-related genes between the squamous and non-squamous subtypes of PDAC were integrated to construct a regulatory network，on the basis of which five immune marker genes regulating the squamous subtype were screened out.An integrated immune score （IIS） model was constructed based on patient survival information and immune marker genes to predict the clinical prognosis of PDAC patients，and its predictive performance was tested with 5 validation sets （n=758）.Results PDAC patients were assigned into high risk and low risk groups according to the IIS.In both training and validation sets，the overall survival of patients in the high risk group was shorter than that in the low risk group （both P<0.001）.The multivariable Cox regression showed that IIS was an independent prognostic factor for PDAC （HR=2.16，95%CI=1.50-3.10，P<0.001）.Conclusion IIS can be used for risk stratification of PDAC patients and may become a potential prognostic marker for PDAC.

Key words：pancreatic ductal adenocarcinoma；immune-related gene；prognosis

Acta Acad Med Sin，2024，46（3）：354-360

胰腺導管腺癌（pancreatic ductal adenocarcinoma，PDAC）是一種常見的消化道惡性腫瘤，占胰腺癌組織學類型的90%以上，是全球第七大癌癥相關死亡原因［1］，在我國其患病人數呈逐年上升趨勢［2］。盡管近些年PDAC的治療策略有所改進，但其預后情況仍然十分嚴峻，5年存活率僅為9%［3］。由于患者個體差異及疾病的異質性，手術和化療的治療效果存在差異［4-5］。對PDAC患者進行危險分層，有助于制訂精準和個體化的治療方案，改善患者預后和臨床轉歸。高通量測序技術的不斷進步，為揭示基因表達和修飾特性提供了有力支持，從而能夠更有效地開發與PDAC預后評估緊密相關的生物標志物［6-8］。根據表達譜聚類結果，PDAC可分為鱗狀亞型、胰腺祖細胞亞型、免疫原性亞型和內/外分泌腺異常分化亞型，其中鱗狀亞型約占31%，侵襲轉移能力強，預后相對較差，中位生存期僅為13.3個月［9］。該亞型的主要分子特征包括TP53和KDM6A基因突變、TP63基因轉錄網絡上調、上皮間質轉化基因高表達以及胰腺內胚層細胞命運決定基因的甲基化［10］。識別鱗狀亞型的內在分子調控機制有助于PDAC患者的風險評估。此外，研究發現免疫微環境也與PDAC預后密切相關［11-12］。本研究擬通過整合分析鱗狀亞型PDAC中的上皮細胞間質轉化特征和免疫相關基因的相互關系，利用調控網絡推理研究主調控鱗狀亞型的免疫標志基因，構建一個風險評估模型，為PDAC患者的精準治療提供新的思路。

1 資料和方法

1.1 數據來源及處理

本研究納入6個公共數據集，共903例PDAC患者，其中GSE71729數據集（n=145）作為訓練集，PACA-AU（n=242）、PACA-CA（n=193）、TCGA-PAAD（n=146）、GSE85916（n=80）和GSE131050（n=97）作為驗證集。采用GEOquery包（V 1.0.7）［13］從基因表達綜合（Gene Expression Omnibus，GEO）數據庫（https：//www.ncbi.nlm.nih.gov/geo/）獲取GSE71729、GSE85916和GSE131050數據集的轉錄譜數據及相應的臨床信息。PACA-AU和PACA-CA數據集中標準化的RNA測序數據和臨床信息來自國際癌癥基因組聯盟網站（https：//icgc.org/）。采用TCGAbiolinks包（V 2.18.0）［14］從癌癥基因組圖譜（The Cancer Genome Atlas，TCGA）數據庫下載TCGA-PAAD數據集的標準化表達譜和相應臨床信息。微陣列數據中的探針ID被轉換成基因符號，若有多個探針ID對應同一種基因符號，則保留平均表達值最高的探針ID來表示該基因。GSE71729、PACA-AU以及TCGA數據集中的分子分型信息來源于Rashid等［15］的研究。

1.2 整合網絡分析

從Immport數據庫［16］獲取2498個免疫相關基因。采用R軟件limma包對PDAC樣本進行差異表達分析，以|log2FC|>0.5和校正后P<0.05為閾值篩選鱗狀和非鱗狀亞型之間差異表達的免疫相關基因。采用RTN包（V 2.10.0）［17］整合差異表達基因并構建調控網絡。采用主調控分析計算各免疫基因調控單元中的鱗狀亞型特征基因的超幾何檢驗P值，進而得到主調控鱗狀亞型的免疫標志基因。

1.3 模型構建和驗證

基于篩選的免疫標志基因和患者生存信息，采用多變量Cox比例風險回歸構建整合免疫評分（integrating immune score，IIS）模型，模型系數是基因的Cox回歸系數，基因值是FPKM。IIS風險賦分公式為：IIS=（0.0068×PLAU）+（0.2584×AQP9）+（-0.3970×COLEC12）+（0.0721×S100A2）+（0.3604×VIM），以上四分位數為閾值將患者分為高危組和低危組。采用Kaplan-Meier生存分析，在5個獨立的驗證集中對IIS模型進行性能評估。此外，納入其他臨床因素作為協變量，采用單變量和多變量Cox回歸分析檢驗IIS是否為PDAC患者的獨立預后因子。

1.4 基因集富集分析

為探究高危組中富集的信號通路，采用R軟件HTS包（V 2.3.5）［18］進行基因集富集分析，設定P<0.05為顯著性閾值，置換檢驗次數設置為1000。

1.5 免疫細胞浸潤分析

采用CIBERSORT算法表征高危組和低危組患者的腫瘤組織基因表達譜中免疫細胞的浸潤豐度。CIBERSORT是一種反卷積算法，可以估算出每個腫瘤組織樣本中22種免疫浸潤細胞的相對比例。將標準化的表達譜上傳至CIBERSORT網站（https：//cibersort.stanford.edu/），運行參數設定為1000次隨機排列，以LM22作為已知參考集，分析樣本中免疫細胞的浸潤情況。

1.6 統計學處理

采用R軟件limma包（V 3.42.2）［19］進行PDAC鱗狀和非鱗狀亞型之間的差異分析。采用R軟件plotly包對PDAC鱗狀和非鱗狀亞型之間的差異基因進行主成分分析（principal component analysis，PCA）。采用列線圖對多個影響因素進行整合并構建評分，繪制校準曲線評估列線圖的預測性能。采用R軟件（V 3.6.1）進行統計分析，P<0.05為差異有統計學意義。

2 結果

2.1 免疫標志基因的篩選

通過差異基因分析從GSE71729數據集中篩選出差異表達的95個免疫相關基因和1269個差異表達基因。PCA法篩選出5個主調控鱗狀亞型的免疫標志基因，分別為PLAU（P<0.001）、AQP9（P=0.011）、VIM（P=0.032）、S100A2（P=0.040）和COLEC12基因（P=0.041）（圖1A）。在TCGA-PAAD驗證集中，與非鱗狀亞型比較，鱗狀亞型PLAU、VIM、S100A2和COLEC12基因的表達均增高（P均<0.001）（圖1B）。

2.2 IIS模型的預測效能

高危組Bailey分型的鱗狀亞型（P<0.001）、Moffitt分型的類基底細胞亞型（P<0.001）及死亡患者（P<0.001）的比例高于低危組（圖2A）。PCA結果顯示，訓練集中IIS可以有效地區分鱗狀和非鱗狀亞型（圖2B）。此外，腫瘤分化程度越低IIS越高，且轉移灶中的IIS高于原發灶和正常組織（P均<0.001）（圖2C）。

2.3 IIS模型驗證

受試者工作特征曲線分析顯示，相較于單個免疫標志基因，IIS對患者2年生存時間的預測能力更強，其曲線下面積達到0.73（圖3A）。生存分析顯示，訓練集（HR=2.3，95% CI=1.4～3.8，P<0.001）和PACA-AU（HR=2.0，95% CI=1.4～2.9，P<0.001）、PACA-CA（HR=1.9，95% CI=1.3～2.8，P<0.001）、TCGA-PAAD（HR=1.6，95% CI=1.0～2.7，P<0.001）、GSE85916（HR=1.9，95% CI=1.1～3.4，P<0.001）、GSE131050（HR=4. 95% CI=1.7～10.0，P<0.001）驗證集中高危組患者的總體生存期均明顯低于低危組（圖3B～G）。納入性別、TNM分期和分化程度等臨床因素作為協變量，進行單變量和多變量Cox回歸分析，結果顯示IIS是PDAC的獨立預后因子（HR=2.16，95%CI=1.50～3.10，P<0.001）。

2.4 列線圖的預測效能

納入IIS與年齡、性別、TNM分期和腫瘤分級等變量構建列線圖，對患者1年和3年生存率進行預測（圖4A）；校準曲線顯示1年和3年生存率的預測校準曲線趨近標準曲線，說明該模型對PDAC患者預后的預測能力較高（圖4B）。

2.5 IIS模型多組學功能分析

基因集富集分析顯示，血管生成、上皮間充質轉化、胰島素樣生長因子受體、哺乳動物雷帕霉素靶蛋白、血管內皮生長因子、表皮生長因子受體、轉化生長因子-β、細胞周期、絲裂原活化蛋白激酶、WNT等相關通路顯著富集于高危組，而DNA修復缺陷通路則在高危組中相對下調（圖5A）。免疫細胞浸潤分析顯示，高危組中M0巨噬細胞比例顯著增加（P<0.001），低危組中單核細胞和未活化的CD4記憶性T細胞比例顯著增加（P<0.00 P=0.009）（圖5B）。

3 討論

PDAC是胰腺癌中最常見的一種類型，其病情隱匿且進展迅速，預后極差。全身化療是針對進展期和復發PDAC患者最主要的治療手段，但由于嚴重的不良反應和患者低耐受性，患者臨床收益有限［20］。隨著分子生物學的發展以及腫瘤免疫微環境的深入研究，以程序性細胞死亡蛋白1/ 程序性死亡配體1抗體和細胞毒性T淋巴細胞抗原4抑制劑為代表的免疫療法在肺癌、黑色素瘤、泌尿系統腫瘤中表現出較好的效果，顯著改善了患者預后［21］。然而對于PDAC免疫療法收效甚微，目前的證據表明處于高度免疫抑制狀態的腫瘤微環境是阻礙免疫療法有效性的主要因素。不同的基質和免疫浸潤細胞類型與PDAC患者的預后相關［22］。因此，探究腫瘤微環境中的免疫狀態對于PDAC整體評估和治療策略的制訂具有重要意義。近年來，癌癥分子分型研究揭示了PDAC的分子特征，基于表達譜的特征聚類分析，PDAC被分為鱗狀亞型、胰腺祖細胞亞型、免疫原性亞型和內/外分泌腺異常分化亞型，其中鱗狀亞型惡性程度最高［9］。本研究篩選了鱗狀亞型和非鱗狀亞型之間的差異表達基因和差異免疫基因，構建了調控網絡，并篩選出主調控鱗狀亞型的5個免疫標志基因，分別是PLAU、AQP9、VIM、S100A2和COLEC12基因。這些基因廣泛參與免疫信號傳導和腫瘤進展。據報道，PLAU在胃癌中呈高表達，其相關基因參與轉化生長因子-β、耐藥形成等信號通路，免疫細胞浸潤狀態削弱，提示患者預后不良［23］。AQP9是存在于多種免疫細胞中的一類水甘油通道蛋白，高表達于乳腺癌、結腸癌和肺癌中，與較差預后顯著相關，且參與調節CD4+T細胞、CD8+T細胞、中性粒細胞等免疫細胞的浸潤［24］。VIM是一種間充質狀態標志物，在多種轉移性腫瘤中呈高表達［25］。S100A2屬于S100蛋白家族，是一類高度保守的EF-手型鈣結合蛋白［26］，其表達失調與多種癌癥的發生發展過程密切相關［27］。Li等［28］研究發現高表達S100A2通過調節糖酵解重編程促進結直腸癌細胞增殖。COLEC12是一種清道夫受體，作為固有免疫系統的模式識別分子，與宿主防御相關的多種功能有關［29］，參與白細胞募集和癌癥轉移［30］。本研究成功構建IIS模型，實現對PDAC患者的危險分層，不同于以往的單基因或多基因標志物，該模型結合免疫系統和轉錄組水平的分子亞型特征，有效減少了腫瘤異質性的干擾。經IIS危險分層后，生存分析結果發現高危組患者的總體生存期明顯低于低危組，且在5個驗證集中的表現一致。納入其他臨床因素作為協變量，單變量和多變量Cox回歸分析證實IIS是PDAC預后的獨立預測因子。

本研究構建了基于公共數據庫的PDAC預后模型，盡管納入了多中心樣本并進行了深入分析和驗證，但仍存在一些局限性。首先，轉錄組測序實驗操作流程復雜，獲取的數據需要依賴專門的計算資源和軟件進行處理和解析；其次，不同測序平臺數據的標準不一，限制了模型的臨床適用性；最后，由于本研究為回顧性分析，樣本量有限，統計效力不足。未來工作需結合更多前瞻性數據和分子生物學實驗以進一步優化和完善該模型。

綜上，本研究結果表明，基于PDAC亞型分類和免疫相關基因構建的IIS模型，在PDAC患者危險分層方面具有一定的應用價值，有助于個體化治療的實施。

利益沖突 所有作者聲明無利益沖突

作者貢獻聲明 張昱：研究選題設計、數據分析、起草和修改文章；任瑞平：數據收集和整理、數據分析、起草論文；萬鵬：數據分析、修訂論文學術內容；賀曉蘭：研究選題設計、文稿審閱和格式修訂

參 考 文 獻

［1］Bray F，Ferlay J，Soerjomataram I，et al.Global cancer statistics 2018：GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries[J].CA Cancer J Clin，2018，68（6）：394-424.DOI：10.3322/caac.21492.

［2］馮程程，許傳志，何杰宇，等.1990-2019年中國20～84歲人群胰腺癌發病趨勢分析[J].中華腫瘤防治雜志，2022，29（18）：1323-1329.DOI：10.16073/j.cnki.cjcpt.2022.18.04.

［3］DeSantis CE，Ma J，Gaudet MM，et al.Breast cancer statistics，2019[J].CA Cancer J Clin，2019，69（6）：438-451.DOI：10.3322/caac.21583.

［4］Springfeld C，Ferrone CR，Katz MHG，et al.Neoadjuvant therapy for pancreatic cancer[J].Nat Rev Clin Oncol，2023，20（5）：318-337.DOI：10.1038/s41571-023-00746-1.

［5］Lo W，Zureikat A.Neoadjuvant therapy in pancreatic cancer：a review and update on recent trials[J].Curr Opin Gastroenterol，2022，38（5）：521-531.DOI：10.1097/MOG.0000000000000874.

［6］Chen H，Yu Y，Zhou L，et al.Cuproptosis-related LncRNAs signature as biomarker of prognosis and immune infiltration in pancreatic cancer[J].Front Genet，2023，14：1049454.DOI：10.3389/fgene.2023.1049454.

［7］Zhang S，Yang J，Wu H，et al.Establishment of a 7-gene prognostic signature based on oxidative stress genes for predicting chemotherapy resistance in pancreatic cancer[J].Front Pharmacol，2023，14：1091378.DOI：10.3389/fphar.2023.1091378.

［8］Chen K，Wang Q，Liu X，et al.Single cell RNA-seq identifies immune-related prognostic model and key signature-SPP1 in pancreatic ductal adenocarcinoma[J].Genes，2022，13（10）：1760.DOI：10.3390/genes13101760.

［9］Collisson EA，Bailey P，Chang DK，et al.Molecular subtypes of pancreatic cancer[J].Nat Rev Microbiol，2019，16（4）：207-220.DOI：10.1038/s41575-019-0109-y.

［10］Giordano M，Cacciato Insilla A，Campani D.Molecular subtypes of pancreatic ductal adenocarcinoma［M］//Esposito I，Karamitopoulou-Diamantis E.Pathology of the pancreas.Cham：Springer International Publishing，2022：136-140.DOI：10.1007/978-3-319-28845-1_5533-1.

［11］Falcomatà C，Brthel S，Schneider G，et al.Context-specific determinants of the immunosuppressive tumor microenvironment in pancreatic cancer[J].Cancer Discov，2023，13（2）：278-297.DOI：10.1158/2159-8290.CD-22-0876.

［12］Sherman MH，Beatty GL.Tumor microenvironment in pancreatic cancer pathogenesis and therapeutic resistance[J].Annu Rev Pathol，2023，18：123-148.DOI：10.1146/annurev-pathmechdis-031621-024600.

［13］Davis S，Meltzer PS.GEOquery：a bridge between the Gene Expression Omnibus （GEO） and BioConductor[J].Bioinformatics，2007，23（14）：1846-1847.DOI：10.1093/bioinformatics/btm254.

［14］Colaprico A，Silva TC，Olsen C，et al.TCGAbiolinks：an R/Bioconductor package for integrative analysis of TCGA data[J].Nucleic Acids Res，2016，44（8）：e71.DOI：10.1093/nar/gkv1507.

［15］Rashid NU，Peng XL，Jin C，et al.Purity independent subtyping of tumors （PurIST），a clinically robust，single-sample classifier for tumor subtyping in pancreatic cancer[J].Clin Cancer Res，2020，26（1）：82-92.DOI：10.1158/1078-0432.CCR-19-1467.

［16］Bhattacharya S，Andorf S，Gomes L，et al.ImmPort：disseminating data to the public for the future of immunology[J].Immunol Res，2014，58：234-239.DOI：10.1007/s12026-014-8516-1.

［17］Fletcher MNC，Castro MAA，Wang X，et al.Master regulators of FGFR2 signalling and breast cancer risk[J].Nat Commun，2013，4（1）：2464.DOI：10.1038/ncomms3464.

［18］Wang X，Terfve C，Rose JC，et al.HTSanalyzeR：an R/Bioconductor package for integrated network analysis of high-throughput screens[J].Bioinformatics，201 27（6）：879-880.DOI：10.1093/bioinformatics/btr028.

［19］Ritchie ME，Phipson B，Wu DI，et al.limma powers differential expression analyses for RNA-sequencing and microarray studies[J].Nucleic Acids Res，2015，43（7）：e47-e47.DOI：10.1093/nar/gkv007.

［20］Von Hoff DD，Ervin T，Arena FP，et al.Increased survival in pancreatic cancer with nab-paclitaxel plus gemcitabine[J].N Engl J Med，2013，369（18）：1691-1703.DOI：10.1056/NEJMoa1304369.

［21］Aroldi F，Zaniboni A.Immunotherapy for pancreatic cancer：present and future[J].Immunotherapy，2017，9（7）：607-616.DOI：10.2217/imt-2016-0142.

［22］Mahajan UM，Langhoff E，Goni E，et al.Immune cell and stromal signature associated with progression-free survival of patients with resected pancreatic ductal adenocarcinoma[J].Gastroenterology，2018，155（5）：1625-1639.e2.DOI：10.1053/j.gastro.2018.08.009.

［23］Ai C，Zhang J，Lian S，et al.FOXM1 functions collaboratively with PLAU to promote gastric cancer progression[J].J Cancer，2020，11（4）：788.DOI：10.7150/jca.37323.

［24］Liu X，Xu Q，Li Z，et al.Integrated analysis identifies AQP9 correlates with immune infiltration and acts as a prognosticator in multiple cancers[J].Sci Rep，2020，10（1）：20795.DOI：10.1038/s41598-020-77657-z.

［25］Baumgart E，Cohen MS，Neto BS，et al.Identification and prognostic significance of an epithelial-mesenchymal transition expression profile in human bladder tumors[J].Clin Cancer Res，2007，13（6）：1685-1694.DOI：10.1158/1078-0432.CCR-06-2330.

［26］Chen H，Xu C，Liu Z.S100 protein family in human cancer[J].Am J Cancer Res，2014，4（2）：89-115.

［27］Wolf S，Haase-Kohn C，Pietzsch J.S100A2 in cancerogenesis：a friend or a foe[J].Amino acids，201 41：849-861.DOI：10.1007/s00726-010-0623-2.

［28］Li C，Chen Q，Zhou Y，et al.S100A2 promotes glycolysis and proliferation via GLUT1 regulation in colorectal cancer[J].FASEB J，2020，34（10）：13333-13344.DOI：10.1096/fj.202000555R.

［29］Ma YJ，Hein E，Munthe-Fog L，et al.Soluble collectin-12 （CL-12） is a pattern recognition molecule initiating complement activation via the alternative pathway[J].J Immunol，2015，195（7）：3365-3373.DOI：10.4049/jimmunol.1500493.

［30］Elola MT，Capurro MI，Barrio MM，et al.Lewis x antigen mediates adhesion of human breast carcinoma cells to activated endothelium.Possible involvement of the endothelial scavenger receptor C-type lectin[J].Breast Cancer Res Treat，2007，101：161-174.DOI：10.1007/s10549-006-9286-9.

（收稿日期：2023-06-20）