卵巢癌免疫相关预后模型的构建与验证

郝苓吉 张哲 董训虎

卵巢癌是女性死亡率最高的妇科肿瘤，2020年全球新增卵巢癌病例31.4万例，死亡病例超过20.7万例，死亡率达60%［1］。肿瘤细胞减灭术后再行辅助化疗是卵巢癌的标准治疗方案。但是，卵巢癌患者确诊时大多已为晚期，且具有易浸润和转移的生物学特性，复发率极高，因此大多数患者对化疗产生耐药性，预后很差［2］，如何筛选出这部分患者对临床治疗具有重要意义。目前，包含基因表达和临床特征的大规模公开数据为研究卵巢癌预后提供了更便利和广泛的机会。

免疫系统已被证明是癌症发生和发展过程中的决定性因素［3］。近期研究表明卵巢癌是免疫原性肿瘤，肿瘤浸润淋巴细胞（tumor infiltrating lymphocytes，TILs）的存在，尤其是CD8+CD3+T细胞，与卵巢癌患者更好的预后有关［4］。此外，已有研究证实免疫系统在卵巢癌中的预后价值［5］，针对免疫检查点的靶向治疗也取得较好的效果［6］。这些证据表明，卵巢癌可以从免疫治疗中获益。因此，基于免疫相关基因的预后模型有潜力应用于卵巢癌的诊断和治疗过程。本研究通过癌症基因组图谱（TCGA）和基因表达数据库（GEO）构建卵巢癌免疫相关的预后模型，探讨其对预后分层的作用，并初步筛选模型中有潜力的预后标志物，以期为卵巢癌的精准免疫治疗提供新策略。

1 资料与方法

1.1 资料收集

本研究纳入的88例卵巢正常组织样本来自基因组织表达数据库（GTEx，https：//xenabrowser.net/data⁃pages/）；卵巢癌样本来自 TCGA（https：//xenabrowser.net/datapages/）和 GEO（https：//www.ncbi.nlm.nih.gov/geo/）数据库。纳入标准：⑴原发疾病诊断为卵巢浆液性囊腺癌；⑵包含完整的生存时间和生存结局；⑶包含完整的临床信息。经匹配，共519例卵巢癌样本纳入后续研究，其中TCGA样本366例，GSE26712样本153例，以上转录数据均经过scale标准化处理。此外，从免疫信息学数据库（https：//www.immport.org/）中收集2 483个免疫相关基因进行差异分析与后续的模型构建。

1.2 差异基因的获取和功能富集分析

提取转录数据（https：//www.immport.org/）中的免疫相关基因，通过R软件（4.0.3）的limma软件包比较GTEx和TCGA样本中具有差异表达的免疫相关基因（IRDEGs），纳入标准：log2FC＞1，FDR＜0.05。利用clusterProfiler软件包对IRDEGs进行KEGG富集分析，通过气泡图对结果进行展示，PFDR＜0.05视为差异有统计学意义。

1.3 预后模型的构建与评价

为构建预后模型，将TCGA卵巢癌样本按照7∶3的比例随机分为训练集和测试集。在训练集中，通过单因素Cox回归筛选与预后有关的IRDEGs，进而通过最小绝对收缩和选择算法（LASSO）回归分析（10折交叉验证）对模型进行精简，最后通过Cox比例风险回归分析得到最优模型，该模型遵循赤池信息准则（AIC），选取AIC值最小的模型作为最优模型，采用R软件的survminer包分析此模型的最佳截断值并对训练集进行分组（组内样本数量大于训练集总样本数的30%），采用Log⁃rank检验比较高/低风险组OS的差异并绘制生存曲线；分析1、3、5年的受试者工作特征（ROC）曲线，通过曲线下面积（AUC）和C⁃index对模型的区分度进行评价，以相同的截断值对测试集和GSE26712队列进行生存分析。然后，在训练集和测试集中分别通过单因素和多因素Cox回归分析年龄、肿瘤分级、肿瘤分期和风险值与卵巢癌预后的关系。最后，通过构建列线图［7］对1、3、5年的生存率进行预测，使用校准曲线对其进行评价（通过bootstrap法循环1 000次），校准曲线可以显示列线图的预测生存率与实际生存率之间的拟合程度，以此来评价列线图的预测准确性；最后通过5年的决策曲线分析评价该列线图在指导临床决策方面的收益情况，并与其他病理参数作比较。上述分析过程由R 4.0.3（http：//www.R⁃project.org）和SPSS 21.0软件完成，除特殊说明外，双侧P＜0.05视为差异有统计学意义。

1.4 卵巢癌预后标志物的筛选

为探讨风险模型中单基因的预后作用，通过TCGA数据库对模型中的每个基因进行生存分析，初步筛选具有统计学意义的关键基因，以P＜0.05视为差异有统计学意义。

2 结果

2.1 卵巢癌中IRDEGs的筛选和功能分析

比较GTEx和TCGA样本中免疫相关基因的转录数据后共得到495个IRDEGs，见图1A。对IRDEGs进行KEGG富集分析，共有146个条目得到富集，排名前10位的条目分别是细胞因子⁃细胞因子受体相互作用（PFDR＜0.001）、EB病毒感染、甲型流感、趋化因子信号通路、病毒蛋白与细胞因子和细胞因子受体的相互作用、Th17细胞分化、自然杀伤细胞介导的细胞毒性、脂质与动脉粥样硬化、卡波西肉瘤相关疱疹病毒感染、抗原处理和呈递（均PFDR＜0.001），见图1B。

图1 IRDEGs的筛选与KEGG富集分析Fig.1 Screening of the IRDEGs and KEGG enrichment analysis

2.2 免疫相关预后模型的构建

将TCGA队列随机分为训练集（n=258）和测试集（n=108），样本信息见表1。单因素Cox回归分析显示，训练集中有26个IRDEGs与预后有关（P＜0.05，图2A）；为避免过度拟合，26个IRDEGs进一步进行LASSO回归分析得到20个关键的IRDEGs（图2B～C）；最后经多因素Cox回归分析确定了包含11 个 IRDEGs（C5AR1、CX3CR1、CXCL11、CXCL13、IGF1、IL27RA、NFKBIB、PENK、PI3、PSMC1和PSME3）的预后模型，该模型的AIC为1 391.44，C⁃index为0.69（图3）。计算公式如下：Risk score=C5AR1×0.226+CX3CR1×0.149+CXCL11×（-0.310）+CXCL13×（-0.192）+IGF1×0.222+IL27RA×0.277+NFKBIB×0.124+PENK×0.212+PI3×0.273+PSMC1×（-0.153）+PSME3×0.188。

图2 卵巢癌免疫相关预后模型的构建Fig.2 Construction of the immune⁃related prognostic signature in ovary cancer

图3 11⁃基因预后模型的森林图Fig.3 Forest plot of the 11⁃gene signature

表1 TCGA样本的基线资料Tab.1 Baseline characteristics of TCGA samples

根据上述公式计算样本的风险值并以最佳截断值（risk score=0.23）为界分别对训练集、测试集和GSE26712队列进行生存分析，同时绘制ROC曲线，结果如图4所示。在训练集中，低风险组的中位生存时间为3.17年（95%CI：0.17～10.59年），高风险组的中位生存时间为2.02年（95%CI：0.06～6.18年），差异有统计学意义（P＜0.001），见图4A；ROC曲线显示其1、3、5年的AUC分别是0.67、0.71、0.75，见图4D。以相同的截断值对测试集和GSE26712队列进行生存分析发现，在测试集中低风险组的中位生存时间（2.99年，95%CI：0.15～11.03年）优于高风险组（2.25年，95%CI：0.04～5.83年），差异有统计学意义（P=0.002），见图4B；1、3、5年的AUC分别为0.68、0.62、0.60，见图4E。在GSE26712队列中低风险组的中位生存时间（3.5年，95%CI：0.32～12.29年）优于高风险组的中位生存时间（2.26年，95%CI：0.12～9.16年），差异有统计学意义（P=0.011），见图 4C；1、3、5年的 AUC 分别是 0.63、0.66、0.63，见图4F。

图4 训练集、测试集和GSE26712队列的生存曲线和ROC曲线Fig.4 Survival curves and ROC curves in the training cohort,testing cohort and GSE26712 cohort

2.3 预后模型和临床参数的关系

为了进一步评估该模型的独立预后效能，本研究通过单因素和多因素Cox回归分析了年龄、肿瘤分期、肿瘤分级和风险值与卵巢癌预后的关系（详见表2），结果显示，在训练集（HR=2.58，95%CI：2.15～3.25，P＜0.001）和测试集（HR=1.68，95%CI：1.17～2.41，P=0.005）中，risk score均能够视为独立的危险因素。

表2 单因素和多因素Cox回归评价预后模型和临床参数的关系Tab.2 Univariable and multivariable Cox regression to evaluate the relationship between prognostic model and clinical parameters

2.4 列线图的构建与评价

根据上述多因素Cox回归的结果，本研究在训练集中进行列线图的构建（图5A），以预测卵巢癌患者1、3、5年生存率，列线图显示风险值与卵巢癌患者OS呈负相关。校准曲线提示该列线图的预测生存率与队列的实际生存率的拟合度较好（图5B）。此外，5年决策曲线发现，与年龄、肿瘤分期、肿瘤分级相比，列线图在指导临床决策方面具有更明显的净收益（图5C）。

图5 训练集中列线图的建立和评价Fig.5 Construction and evaluating of the nomogram in the training cohort

2.5 卵巢癌中关键基因的初步筛选

通过TCGA数据库对预后模型中的11个基因进行生存分析，以初步筛选卵巢癌潜在的预后标志物。结果显示，CXCL11（P＜0.001）和CXCL13（P＜0.001）表达越高，卵巢癌患者预后越好；而C5AR1（P＜0.001）、CX3CR1（P＜0.001）、IL27RA（P=0.008）、PENK（P=0.020）、PI3（P＜0.001）和PSME3（P=0.003）表达越高，卵巢癌者的预后越差。见图6。

图6 预后模型中关键基因的生存分析Fig.6 Survival analysis of the hub genes in the prognostic model

3 讨论

近年来，随着免疫治疗在肿瘤领域取得的巨大进展，其作为癌症的治疗手段或标准化疗的替代方案已被广泛接受［8］，基于免疫特性的预后模型也成为预测癌症患者风险的研究重点之一［9］。本研究通过多个公开的卵巢癌队列，构建并验证了免疫相关基因的预后模型，结果发现该模型可以有效对卵巢癌患者进行危险分层，预测效能优于传统的临床病理因素，且其中的关键基因可能是卵巢癌潜在的预后标志物。首先，本研究从数据库中收集免疫相关基因，分析比较得出卵巢癌中的IRDEGs，富集分析发现这些IRDEGs大多参与了病毒感染、细胞因子和细胞因子受体的相互作用以及趋化因子信号通路。已有临床研究表明，辅助抗病毒治疗的抗增殖和细胞毒性作用可保护癌症患者免受复发感染，从而使患者受益［10］。在卵巢癌中，INGERSLEV等［11］报道，与正常组织相比，上皮性卵巢癌组织中的EB病毒感染率更高，有可能是致癌因素。此外，细胞因子和细胞因子受体的相互作用可以广泛参与炎症、血管生成和趋化过程，抑制肿瘤发生和发展［12］。趋化因子是一类由细胞分泌的小细胞因子或信号蛋白，具有诱导附近反应细胞定向趋化的能力。在本研究构建的IRDEGs中，CXCL11和CXCL13隶属于趋化因子家族并且与卵巢癌患者较好的预后有关［13］。而另外一个IRDEG CX3CR1是趋化因子CX3CL1的单一受体，在抗病毒和抗肿瘤免疫中发挥重要作用［14］。以上结果说明本研究构建的IRDEGs可能通过参与趋化因子信号通路从而在肿瘤的进展过程中发挥至关重要的作用。

本研究构建了卵巢癌免疫相关的预后模型，该模型具有中等程度的预测效能，所构建的列线图能够较好地预测卵巢癌患者的OS，且优于传统的临床指标。值得注意的是，无论是在训练集还是在测试集，年龄、肿瘤分期与分级等传统指标均不能够独立预测卵巢癌患者的预后，因此本研究最终选择预后模型进行列线图构建，有理由推测在更大样本量的研究背景下，该模型可以结合临床参数构建更加准确的列线图，同时也提示本研究结论仍需在更独立的队列中进行验证的必要性。近年来，陆续有学者关注卵巢癌的免疫相关预后模型研究，如CAO等［15］利用TCGA数据库成功构建21⁃免疫基因的预后模型，其预测5年总体生存率的AUC值为0.749。也有学者［16⁃17］利用单样本基因集富集分析（ssGSEA）的方法评价卵巢癌样本和免疫细胞的关系，并进行打分，根据样本得分情况将队列分为高/低免疫评分组，从而进行后续的建模，预测5年生存率的AUC分别为0.724和0.704。本研究通过Immport网站提取转录组中免疫相关的基因并进行后续的模型构建，预测卵巢癌5年OS的AUC值为0.75，同时该模型进一步在测试集和GSE26712队列中得到验证，其具有中等程度（AUC＞0.7）的预测效能，且较21⁃免疫基因模型更简化，但不足的是该模型未能在更大样本量的临床队列中进一步验证与评价。

本研究还对模型中的关键基因进行生存分析，结果发现，其中8个基因与卵巢癌患者预后有关，包括C5AR1、CX3CR1、CXCL11、CXCL13、IL27RA、PENK、PI3和PSME3。CXCL11位于人类染色体4q21.2上，属于ELRCXC趋化因子家族，其受体为CXCR3和CXCR7，癌细胞可以通过自分泌产生CXCL11或通过调节微环境中的肿瘤基质细胞释放CXCL11［18］。CXCL13可以在高级别浆液性卵巢癌中形成免疫活性肿瘤微环境并增强PD⁃1检查点阻断的疗效［19］。C5AR1和体内T细胞的归巢以及免疫治疗有关［20］，而CX3CR1的表达受缺氧诱导因子⁃1α调控并能促进卵巢癌细胞的上皮⁃间充质转化［21］。PSME3被报道能促进胰腺癌细胞分泌TGFB1从而诱导胰腺星状细胞增殖［22］，但是与卵巢癌的关系尚未见报道。PI3是一种乳清酸性蛋白，在卵巢浆液性囊腺癌中，PI3可以通过MAPK信号通路促进细胞有丝分裂，且与更差的预后有关［23］。IL27RA则被认为是促炎症细胞因子，在野生型幼稚CD4+T细胞中，能通过STAT1的磷酸化诱导T⁃bet表达，说明IL27RA对Th1细胞的免疫应答至关重要［24］。本研究结果也提示上述基因有潜力成为卵巢癌的预后标志物，但其确切的作用有待进一步功能和机制研究探索。

尽管本研究的结果具有潜在的实质性临床意义，但仍有以下问题需要注意及探讨：⑴本研究纳入了多个队列进行研究，试图对模型进行更充分验证，但是由于平台的不同及人群异质性，基因表达水平易受采样偏差影响；⑵本研究纳入的病例数量有限，研究结果仍需更大规模的临床前瞻性研究验证；⑶本研究构建的预后模型由多个基因组成，因此有必要开展进一步的生物学实验探索其在卵巢癌中的功能和具体机制。

总之，本研究基于IRDEGs构建卵巢癌预后模型并筛选预后标志物，有助于对患者进行分层并预测生存，这些发现为卵巢癌的免疫和靶向治疗提供了新的见解。