卵巢癌（ovarian cancer，OC）是女性生殖系统常见的恶性肿瘤之一^[1]。原发性卵巢癌的组织学类型中90%以上为上皮性卵巢癌（epithelial ovarian cancer，EOC），EOC的发展和进展与上皮组织密切相关^[2-3]。正常上皮细胞具有抗肿瘤活性，并能够通过调节细胞骨架蛋白来消除致癌转化细胞^[4]。上皮细胞是大多数人类肿瘤的来源，其恶性转化通常与细胞极性的丧失和解体密切相关，并且上皮细胞极性的破坏促进上皮-间质转化（epithelial-mesenchymal transition，EMT），这是上皮肿瘤细胞侵入周围基质的关键步骤^[5-6]。随着对上皮细胞研究的不断深入，上皮细胞相关生物标志物成为近年来的研究热点。单细胞转录组测序（single cell RNA sequencing，scRNA-seq）是在单细胞水平对转录组进行测序的1项新技术，可以研究单个细胞内的基因表达情况，能解决用组织样本测序无法解决的细胞异质性难题^[7-8]。

为了探索上皮细胞相关基因与卵巢癌的关系，笔者使用生物信息学方法，基于卵巢癌的scRNA-seq数据对卵巢癌上皮细胞进行了更精确的分析，并联合bulk RNA-seq构建了1个基于上皮细胞标记基因的风险评分模型，可能为OC患者提供潜在的治疗靶点。

1.   材料与方法

1.1   获取卵巢癌单细胞测序的转录组和临床数据

从GEO数据库（https://www.ncbi.nlm.nih.gov/geo/）获取OC单细胞测序数据（GSE118828）进行单细胞分析（包含18个样本）；获得具有完整临床信息的数据集GSE140082作为模型的验证队列（包含380个OC的bulk RNA测序样本）。从TCGA数据库（https://portal.gdc.cancer.gov/）获得379份具有生存时间及生存状态的OC样本（TCGA-OV）作为模型的试验队列，从GTEx数据库获取88份非肿瘤卵巢组织的转录组数据用于差异分析。

1.2   单细胞数据处理及细胞聚类

使用R语言（4.1.2版）中的“Seurat”和“Harmony”包对GSE118828数据集进行整理并去除批次效应。“PercentageFeatureSet”包用于确定每个细胞中线粒体基因的百分比后删除线粒体基因比例>15%的细胞，避免线粒体基因的表达影响细胞分群，同时删除状态差的细胞。以上数据标准化后使用“FindVariableFeatures”函数识别前1500个高度可变基因（highly variable features，HVGs）用于细胞分群。使用“JackStraw”函数基于HVGs进行PCA以降低维度，选择前10个PC对细胞进行聚类。应用“FindClusters”函数，根据已报道的细胞特异性标记基因对细胞簇进行注释^[9-15]。

1.3   单细胞数据的拷贝组变异（copy number variation，CNV）分析

使用“copyKAT”包进行CNV分析，识别每个样本的染色体拷贝数变异情况，鉴定良恶性细胞。选取在良恶性细胞中表达差异具有统计学意义（P < 0.05）的基因作为恶性细胞标记基因。

1.4   上皮细胞相关评分模型的构建和验证

单因素Cox分析用于识别TCGA-OV队列中与生存显著相关的恶性上皮细胞标志基因。使用“Glmnet”包整合生存时间、生存状态和基因表达量，利用Lasso-cox方法进行回归分析。设置10折交叉验证以获得最优模型。Lambda值为0.00235651033076621，确定模型公式为：RiskScore=-0.169036513490011*EPCAM+0.194701577619064*CLDN4-0.139298669105292*CXCR4+0.155318080556463*TIMP3。使用“Maxstat”包计算RiskScore的最佳截断值，基于此截断值将患者分成高、低风险2组，使用“Survival”包分析2组的预后差异。

获取患者的OS及模型风险评分，利用pROC的ROC函数进行1、3、5 a的ROC分析，并评估AUC和置信区间以获得最终的AUC结果。

1.5   TCGA队列中的功能分析

使用“Limma”包在高低风险组（H.Risk和L.Risk）间进行差异表达分析，并得到差异表达基因（differentially expressed genes，DEGs）（∣logFC∣≥1且P < 0.05）。随后，使用“clusterProfiler”包和“org.Hs.eg.db”包确定了DEGs的潜在生物学机制。基于基因本体论（gene ontology，GO）和基因组百科全书（kyoto encyclopedia of genes and genomes，KEGG）对高低风险组的DEGs进行功能富集分析，确定差异具有统计学意义的生物学功能及通路（P < 0.05）。

2.   结果

2.1   高质量细胞的识别及主成分分析

通过计算，发现GSE118828数据集中测序深度与基因数目呈正相关，所有细胞中均未测到线粒体基因，见图1A～1B。消除在<10个细胞中表达的基因和表达<200个基因的细胞，共得到3066个高质量细胞。获取了在样本中波动最大的1500个高变基因（HVGs），见图1C。使用HVGs进行主成分分析，可视化结果显示前15个PC均在虚线上方且靠近y轴，越靠近y轴表明实际基因与理论基因的差值越小（P < 0.05），见图1D。

图  1  卵巢癌的单细胞RNA测序分析

A：线粒体基因与测序深度的关系（线粒体基因含量为0）和测序深度与基因数目呈正相关关系；B：线粒体基因含量；C：前1500个高变基因；D：前15个主成分。

Figure  1.  Single cell RNA sequencing analysis of ovarian cancer

下载: 全尺寸图片 幻灯片

2.2   鉴定上皮细胞类群

将3 066个细胞样本根据前10个PC进行T-SNE聚类，共得到13个亚群，见图2A。差异分析表明各亚群之间基因表达具有显著差异P，见图2C。对细胞亚群进行注释共得到8个主要细胞群（上皮细胞、T细胞、B细胞、单核细胞、成纤维细胞、组织干细胞、平滑肌细胞和内皮细胞），上皮细胞在其中占了相对较大的比例，见图2B。筛选出56个在上皮细胞中表达显著的基因用于后续分析（P < 0.05）。

图  2  单细胞数据的降维和聚类

A：细胞的T-SNE图，显示细胞群；B：细胞注释的T-SNE图；C：细胞群标记物的相对表达热图（仅显示前10名）。

Figure  2.  Dimensionality reduction and clustering of single cell data

下载: 全尺寸图片 幻灯片

2.3   鉴定良恶性细胞并筛选恶性基因集

CNV将细胞样本区分为1 530个恶性细胞和1 536个非恶性细胞，两类细胞间基因表达具有明显差异P，见图3A。筛选出恶性和非恶性细胞间的657个差异基因作为恶性细胞基因集（P < 0.05）。恶性细胞基因集与上皮细胞基因集取交集共得到45个基因。对以上基因在TCGA-OV和GTEx中的表达进行差异分析，得到36个在肿瘤组织和正常卵巢组织中差异显著的恶性上皮细胞标记基因（MECRGs），见图3B。

图  3  单细胞CNV分析及预后基因筛选

A：拷贝组变异分析热图；B：MECRGs在良恶性患者中的表达热图；C：基因表达和OS之间单变量Cox回归分析结果的森林图。

Figure  3.  Single cell CNV analysis and prognosis gene screening

下载: 全尺寸图片 幻灯片

2.4   筛选预后相关基因

对MECRGs进行单因素Cox回归分析，得到4个与预后显著相关的基因（P < 0.05），分别为EPCAM、CLDN4、CXCR4、TIMP3，其中EPCAM、CXCR4为保护性因素，CLDN4、TIMP3为预后危险因素，见图3C。

2.5   预后模型的构建及验证

TCGA-OV作为试验组，Lasso-cox分析显示以上4个基因均参与模型构建，见图4A～4B。其中EPCAM、CXCR4的表达量随着风险评分的增加呈现下调趋势，CLDN4、TIMP3的表达量随着风险评分的增加呈现上调趋势，见图4C。K-M生存分析显示H.Risk组患者的OS明显低于L.Risk组患者（P < 0.05），见图4D。ROC曲线显示AUC值在1 a、3 a、5 a分别达到0.55、0.60、0.62，见图4E。

图  4  TCGA-OV中的多基因风险评分构建与验证

A～B：Lasso-cox回归分析；C：TCGA队列中OS状态、OS和风险评分的分布、模型基因的表达热图；D：TCGA-OV队列高、低风险组患者OS的K-M曲线（P < 0.001）；E：TCGA-OV队列ROC曲线；F：GEO队列中OS状态、OS和风险评分的分布、模型基因的表达热图；G：GEO队列高低风险组患者OS的K-M曲线（P < 0.001）；H：GEO队列ROC曲线。

Figure  4.  Construction and verification of polygene risk score in 4TCGA-OV

下载: 全尺寸图片 幻灯片

模型基因在GSE140082验证集中的上下调趋势与在试验集中一致，见图4F。K-M生存分析显示H.Risk组患者的OS明显低于L.Risk组（P < 0.05），见图4G。ROC曲线显示AUC值在1a、3a分别为0.51、0.58，见图4H，表明随着观察时间的延长，模型准确性逐渐提高。

2.6   TCGA队列中上皮细胞相关标记的功能富集分析

根据风险评分的最佳截断值-0.450193847032549将TCGA-OV为高、低风险2组，差异分析显示H.Risk组中311个基因上调，174个基因下调。DEGs的GO、KEGG富集分析表明，2组之间的差异与细胞粘附、细胞外基质、神经活性配体-受体相互作用、钙信号通路、转化生长因子-β相关，见图5A～5B。

图  5  路径分析的主要结果

A～B：在高低风险组中差异基因的GO和KEGG富集分析。

Figure  5.  Main results of path analysis

下载: 全尺寸图片 幻灯片

3.   讨论

3.1   上皮相关风险模型具有一定预后价值

由于诊断和治疗的延误，OC患者不可避免地出现不良预后。因此，有效的预后生物标志物仍然是OC患者迫切需要的。OC的发生和发展与上皮组织之间存在密切关联，因此，需要进一步研究上皮细胞相关基因和OC患者预后之间的潜在相关性。

这项研究首次建立了1个基于上皮细胞相关基因的OC风险评分模型，笔者用严格的标准和可信的算法创建了该模型。为证明风险模型的实用价值，本研究将TCGA-OV和GSE140082分别作为试验集和验证集，并根据风险评分将患者分为高、低风险组，以验证所构建的风险模型的预后价值。本究结果显示在试验集和验证集中低风险组患者的预后均优于高风险组。

3.2   模型基因可能成为卵巢癌的潜在靶点

本研究构建的模型包括4个基因（EPCAM、CLDN4、CXCR4和TIMP3）。EPCAM表达能通过阻止细胞间黏附、促进免疫逃逸或激活致癌基因来促进癌症侵袭。既往研究表明，EPCAM在原发食管癌中高表达，在早期食道癌中EPCAM表达减少会诱导EMT，从而促进癌症进展；而在肝癌中，EPCAM上调有促进血管生成的作用，EPCAM阳性的肝细胞癌患者预后更差、复发风险更高^[16-19]。CLDN4在各种类型肿瘤中广泛出现高表达或表达缺失。当CLDN4相关蛋白发生变化时，将影响细胞间的通透性，可引起疾病的发生^[20]。既往研究表明CLDN4在生殖系统肿瘤如卵巢癌、宫颈癌中呈现异常的高表达或低表达，CLDN4基因表达增加可以促进卵巢癌细胞的EMT进程，与卵巢癌的不良预后相关^[21-23]，该基因在卵巢癌中可能有着重要的应用价值。CXCR4是1种已被证实的治疗靶点，当激活受体的信号途径调控失常时，CXCR4可引起癌细胞生长及扩散。而且CXCR4的高表达还与许多癌症亚型的不良预后和化疗抵抗相关，部分是通过增强癌症和基质之间的相互作用^[24]。TIMP3是1种与细胞外基质紧密结合的蛋白质，TIMP3对MMP的调控可以抑制肿瘤的生长、肿瘤细胞的侵袭和迁移。然而，TIMP3的调控机制目前尚不清楚^[25]。综上所述，EPCAM、CLDN4、CXCR4和TIMP3基因与肿瘤的发生及进展密切相关，但在卵巢癌中的研究相对较少。因此，对以上4个基因继续深入研究有可能为卵巢癌的靶向治疗提供新的依据。

高、低风险组间差异基因的功能富集分析也表明，差异基因主要与细胞粘附、细胞外基质、神经活性配体-受体相互作用、钙信号通路、转化生长因子-β 相关。所有这些信号都被报道促进了肿瘤的发展。在前列腺癌中，细胞粘附性的增加会诱导肿瘤细胞发生上皮间质转化，促进细胞侵袭和转移^[26]。在膀胱癌中，转化生长因子β1在肿瘤早期抑制肿瘤增值，在晚期则通过促进上皮间质转化、诱导免疫抑制与肿瘤微环境形成，从而促进肿瘤进展^[27]。因此，参与模型的基因也可能通过调节以上通路来影响卵巢癌的发生与进展。

3.3   本研究存在的局限性

本研究仍然存在一定的局限性。首先，一些临床信息的不可获得性阻碍了笔者充分探索模型和临床特征之间的关系，例如治疗的细节、肿瘤的病理学细节。其次，虽然这项研究使用了外部队列来验证模型的可靠性，但是仍需要更多的数据来进行验证。