《Frontiers in Immunology》：Construction and validation of a novel diagnostic model for esophageal squamous cell carcinoma: an integrated analysis of multi-omics data
引言
食管癌是全球范围内严重威胁人类健康的恶性肿瘤，发病率排名第七，死亡率排名第六。2020年全球癌症统计数据显示，中国的新发食管癌病例和死亡人数分别占全球的53.70%和55.35%。其中，食管鳞状细胞癌（ESCC）占中国病例的95%以上，与西方以食管腺癌（EAC）为主的情况形成鲜明对比。遗憾的是，ESCC患者在早期通常没有明显症状，多数在出现吞咽困难或转移症状后才被确诊，大部分患者确诊时已处于中晚期，总体五年生存率仅为6%-15%。若能早期发现并及时治疗，患者生存率将显著提高。
目前，针对ESCC可靠的早期肿瘤标志物和诊断模式仍然明显不足。传统的肿瘤标志物，如鳞状细胞癌抗原、癌胚抗原和细胞角蛋白19片段，其特异性或敏感性不足以作为ESCC的独立诊断生物标志物。MicroRNAs在早期ESCC检测中显示出潜力，但现有的诊断模型仍有局限。因此，开发新的诊断方法以实现ESCC的早期检测显得尤为重要。
大多数肿瘤诊断模型主要使用来自癌症基因组图谱（TCGA）数据库的数据构建，而该数据库的样本绝大多数来自欧洲血统人群。然而，ESCC的高危人群主要集中在亚洲，特别是中国和日本。仅依赖TCGA数据开发ESCC诊断模型存在明显局限性。本研究使用的GSE23400、GSE17351和GSE20347数据集主要包含来自中国和日本人群的样本。利用这些数据集构建ESCC诊断模型并鉴定抑癌基因，可为ESCC的预防和治疗策略提供有价值的见解。
在本研究中，我们基于GEO数据库筛选出五个与ESCC显著相关的基因，构建了诊断模型和列线图，并通过动物模型中肿瘤与正常组织之间基因差异表达模式的实验证据，以及人ESCC组织免疫荧光分析对保护因子进行了补充验证。此外，通过整合单细胞RNA分析工具，我们深入探究了基因可能参与的生物学过程，为ESCC患者的临床治疗提供了潜在的见解。
材料与方法
研究设计
本研究的流程如图1所示。为获取ESCC样本，我们检索并整合了三个GEO数据集（GSE23400、GSE17351和GSE20347），并移除了批次效应。我们进行了差异表达基因（DEGs）分析，以识别ESCC与癌旁正常组织之间的差异表达基因。采用加权基因共表达网络分析（WGCNA）来描绘与ESCC关联最强的基因模块，并通过WGCNA和DEGs的交集来精炼ESCC相关基因集。为识别关键通路和生物学过程，我们对合并数据集、DEGs和WGCNA结果进行了基因集富集分析（GSEA）。使用LASSO回归识别潜在的生物标志物，随后通过逻辑回归构建ESCC诊断模型，并在TCGA ESCC数据集中进行外部验证。为表征肿瘤与癌旁正常组织的免疫微环境，进行了免疫浸润分析。特别检查了保护基因SORBS2对免疫微环境的影响，并通过GSEA进一步探索其功能。为验证候选生物标志物在正常和肿瘤组织之间的差异表达，我们建立了小鼠ESCC模型并进行了蛋白质免疫印迹（WB）实验。进行了单细胞RNA分析以探索肿瘤组织的细胞类型及保护基因的潜在功能。最后，对人ESCC标本进行免疫荧光染色，以确认SORBS2蛋白的表达和定位。
数据获取与预处理
批量转录组数据集
本研究中的所有数据集均从GEO数据库下载。GSE23400数据集包含53个癌旁正常食管组织和53个ESCC组织。GSE17351数据集包含5个正常组织和5个ESCC组织。GSE20347数据集中，共有17个癌旁正常组织和17个ESCC组织。
批次效应消除与数据合并
应用稳健多阵列平均（RMA）技术来校正背景噪声、标准化信号并计算表达值。利用R软件环境中的inSilicoMerging包整合上述三个数据集。随后实施策略以去除批次效应，从而得到正常和食管组织的一致表达矩阵。
ESCC相关基因的鉴定
差异表达基因
我们评估了GSE23400、GSE17351和GSE20347中不同样本的基因表达水平分布。应用主成分分析（PCA）和样本层次聚类树来探索组间差异和组内样本重复性。异常样本从后续分析中移除。使用limma包筛选DEGs，截断标准为P< 0.05 且 |log2FC| > 0.585，并进一步通过pheatmap和ggplot2包在热图和火山图中可视化。
差异表达基因的基因本体论和京都基因与基因组百科全书富集分析
使用clusterProfiler R包进行基因本体论（GO）分析。对于功能分析，参数设定显著性水平为P< 0.05，从而识别分子功能（MFs）、生物过程（BPs）和细胞组分（CCs）。此外，使用clusterProfiler包设定P< 0.05进行功能分析。该方法能够通过KEGG功能富集分析揭示与可能作为治疗靶点的基因相关的富集信号通路。
加权基因共表达分析
使用WGCNA R包创建涵盖ESCC和正常组织样本中所有基因的共表达网络。为进一步分析，选择了方差高达50%的基因。通过确定基因间的皮尔逊相关系数生成共表达矩阵。随后，使用公式 amn= |cmn|β制定邻接矩阵，其中amn表示基因m和n之间的邻接性，cmn是皮尔逊相关系数，β代表软阈值功率。然后，将该加权邻接矩阵转换为拓扑重叠测度矩阵，以评估基因的网络连接性。通过应用平均连锁层次聚类生成矩阵的聚类树状图。设定最小基因模块大小为30以识别合适的模块，并设定0.25的阈值以合并相似模块。
针对ESCC保护基因的基因集富集分析
在本研究中，使用分子特征数据库（MSigDB）集合，通过clusterProfiler包执行GSEA。目的是识别ESCC或正常组织之间通路的显著差异。ESCC或正常组织之间的区别作为分析的表型标签，置换次数配置为1000。其余参数保持默认设置。调整后P< 0.05 且 FDR < 0.25的通路被认为是显著富集的。
针对ESCC保护基因的免疫细胞浸润分析
CIBERSORT是一种计算工具，能够基于基因表达数据推断基因表达谱并计算异质细胞群体内各种细胞类型的相对丰度。我们使用CIBERSORT评估了22种不同免疫细胞类型的相对比例。此外，我们应用两样本t检验来探索正常组与ESCC组之间，以及低表达SORBS2组与高表达SORBS2组之间免疫细胞浸润的差异。
新型诊断模型的构建与性能评估
样本量计算
为达到80%的检验效能，双侧α为0.05，假设比值比（OR）为2.1，我们使用R中的powerMediation包计算所需样本量。假设ESCC组与正常对照组比例为1:1，计算所需总样本量为245人。
新型诊断模型的构建与性能评估
我们利用合并GEO队列中126个ESCC组织样本及其癌旁正常组织的基因表达数据开发了一种新的诊断模型，两组共包含252个样本。最初，我们使用glmnet包中的LASSO技术选择与ESCC密切相关的基因以防止模型过拟合，并进行50折交叉验证。随后，基于LASSO识别的生物标志物，使用多元逻辑回归构建诊断模型，并使用rms包以列线图形式图形化表示。为量化过拟合，我们在训练集（n = 173）内进行了1000次Bootstrap重采样，并计算了乐观校正R2。此外，我们从TCGA数据库中提取了所有ESCC样本，共91个样本，作为外部验证集。从临床元数据中提取TCGA-ESCC样本的种族信息，自我报告的种族分为亚洲人、白人、黑人或非裔美国人、或未知。我们通过受试者工作特征（ROC）曲线分析评估模型在训练集和外部验证集上的预测性能，评估指标包括曲线下面积（AUC）值、精确度、特异性、敏感性、阳性预测值（PPV）和阴性预测值（NPV）。我们还使用校准曲线分析来评估预测概率与ESCC实际发生频率之间的一致性，并采用决策曲线分析（DCA）来评估列线图的临床效用。
小鼠ESCC组织中标志基因差异表达的验证
ESCC细胞培养
小鼠ESCC细胞系AKR购自Cellverse有限公司。细胞在补充有10%胎牛血清和1%青霉素-链霉素的DMEM中培养，直至形成汇合单层。当细胞生长到适合传代的密度时，用胰蛋白酶消化并以1:3的比例传代，并在37°C和5% CO2条件下培养。离心后，将4 × 10^6个AKR细胞重悬于PBS和基质胶（4:1比例）的混合物中，其中PBS购自中国碧云天生物技术，基质胶购自美国康宁。
ESCC小鼠模型的构建
本研究中进行的所有动物实验均获得了中国医学科学院肿瘤医院动物实验伦理委员会的批准（批准号：NCC2024A559）。将八只5周龄雄性C57BL/6小鼠（北京华阜康生物科技股份有限公司）适应7天，然后随机分为两组，每组四只：ESCC组和正常组。在ESCC组中，将1×10^6个AKR细胞皮下注射到每只小鼠的右腋窝区域，以建立小鼠ESCC模型。从第三天开始测量肿瘤大小。使用公式计算肿瘤体积：V = 1/2 × (长 × 宽^2)，并在尸检时记录肿瘤重量。接种后第10天，在CO2室中对小鼠实施安乐死。立即从ESCC小鼠收集肿瘤组织，并从对照组收集正常食管组织。
5种生物标志物表达的蛋白质免疫印迹验证
将速冻的ESCC肿瘤或癌旁正常上皮（20 mg，每组n = 4）在冰上切碎，并在含有蛋白酶和磷酸酶抑制剂以及3毫米钢珠的150 μL冰冷RIPA裂解液中裂解。组织在60 Hz下匀浆5个循环，并在4°C下以12,000 g离心20分钟。用BCA试剂盒对上清液进行定量，并用裂解液调整至2 μg/μL。等量蛋白质（20 μg）与5× Laemmli缓冲液混合，在100°C下变性15分钟，上样到10% SDS-PAGE凝胶（含5%浓缩胶）上，在120 V下电泳，直至溴酚蓝到达凝胶底部。通过湿转法将蛋白质转移到0.22 μm PVDF膜上，用5%脱脂牛奶的TBST溶液在室温下封闭1小时，并在4°C下与一抗孵育过夜：抗BID（1:1000）、抗CBX3（1:1000）、抗ECT2（1:2000）、抗KIF14（1:1000）、抗SORBS2（1:500）。TBST洗涤三次，每次5分钟后，将膜与HRP标记的山羊抗兔IgG（1:5000）在室温下孵育1小时。使用ECL显色，并使用六一仪器化学发光成像仪成像。在ImageJ中量化条带密度，归一化至β-肌动蛋白，并使用GraphPad Prism绘制为平均值±标准误。
单细胞RNA测序数据获取与分析
从GEO检索公共scRNA-seq数据集（GSE196756），并使用Seurat进行处理。原始读数使用Cell Ranger比对到GRCh38。符合以下条件的细胞未纳入本研究：（1）检测到的基因计数在1,000–35,000范围外的细胞被丢弃；（2）检测到的RNA计数<200或>7,500的细胞被移除；（3）线粒体转录本超过其总基因计数10%的任何细胞被过滤掉。使用Harmony校正批次效应。通过PCA和UMAP执行降维。聚类使用Louvain算法。通过典型标记物注释细胞类型。差异表达分析识别了簇之间的DEGs。通过CellChat预测细胞-细胞通讯。
免疫荧光测定
本研究中使用所有人组织的获取均获得了中国医学科学院肿瘤医院伦理委员会的批准（批准号：23/510-4253）。将福尔马林固定石蜡包埋（FFPE）切片（5 μm）在二甲苯中脱蜡，通过梯度乙醇（100%→95%→80%→70%；各3分钟）再水化，并在10 mM柠檬酸盐缓冲液（pH 6.0，95°C，15分钟）中微波加热。冷却至室温后，切片用5%牛血清白蛋白（BSA/PBS）在25°C下封闭30分钟。与一抗——转胶蛋白/SM22多克隆抗体（1:200）和SORBS2多克隆抗体（1:500）——在4°C潮湿避光条件下共孵育过夜。PBS洗涤三次后，切片与Alexa Fluor 488标记的山羊抗小鼠IgG（1:500）和Alexa Fluor 594标记的山羊抗兔IgG（1:500）在室温避光条件下孵育1.5小时。PBS充分洗涤后，用含DAPI的抗淬灭封片剂封片，在4°C黑暗中固化过夜，并使用配备Plan-Apochromat 63×/1.40油镜的Zeiss LSM 900共聚焦显微镜成像。顺序扫描防止光谱串扰。使用ImageJ对每个样品≥5个随机选择的视野进行荧光强度定量分析，并进行背景扣除。
统计分析
所有生物信息学分析均在R v4.3.1中进行。关键软件包包括：用于差异表达分析的limma；用于共表达网络构建的WGCNA；用于数据可视化的ggplot2；用于批次效应校正的inSilicoMerging。实验数据以平均值±标准差表示。对于两组间数据比较，应用t检验，而对于多组间比较，则使用单因素方差分析。使用GraphPad Prism 10软件进行统计分析和绘图。P值 < 0.05表示具有统计学显著性差异。
结果
差异表达基因分析揭示113个显著上调和173个显著下调基因
如图2A、B所示，去除批次效应前，各数据集间的样本分布存在显著差异，箱线图突显了批次效应的存在。去除批次效应后，各数据集间的数据分布变得更加一致。通过倍数变化（FC）和P值过滤（|log2FC| > 0.585 且 P< 0.05），在鳞癌组织中总共发现113个基因显著上调，173个基因显著下调（图2C、D）。
差异表达基因和合并数据集的基因集富集分析揭示了关键通路和生物过程
我们对DEGs的基因进行了GO和KEGG富集分析。在生物过程方面，主要富集于细胞分化、细胞粘附和有丝分裂细胞周期G2/M转换的正调控（图3A）。在分子功能方面，主要富集于细胞骨架蛋白结合、泛素样蛋白连接酶结合和生长因子结合（图3B）。在细胞组分方面，主要富集于细胞外基质组分、细胞外空间和细胞外基质组分（图3C）。在KEGG通路方面，主要富集于血小板活化、半胱氨酸和蛋氨酸代谢以及TGF-β信号通路（图3D）。
然后，对合并数据集进行GSEA。在生物过程方面，p53类介导的信号转导、DNA依赖性DNA复制和mRNA转运在ESCC中均显著上调（图3E）。在分子功能方面，snoRNA结合、多巴胺受体结合和苏氨酸型肽酶活性在ESCC中均显著上调（图3F）。在细胞组分方面，宿主细胞组分、浓缩染色体和核孔在ESCC中均显著上调（图3G）。在KEGG通路方面，DNA复制、细胞周期和剪接体在ESCC中均显著上调（图3H）。
加权基因共表达网络分析揭示与ESCC最相关的蓝色模块，包含13个基因
使用WGCNA包建立基因共表达网络。生成了描绘无标度拓扑模型拟合系数的散点图（图4A），红色五角星标记了超过无标度拓扑模型拟合指数R2 > 0.85的初始点，对应的软阈值β为4，R2值为0.88。使用此软阈值β将相关矩阵转换为邻接矩阵后，组装了拓扑重叠矩阵和基因层次聚类树状图。使用共表达聚类树状图描绘了各模块之间的距离（图4B），从而识别出5个基因模块，如图4C所示，突出了不同模型之间相关性的基因聚类树状图。如图4D所示，蓝色模块与食管鳞状细胞癌的相关性最高（相关指数：0.82，P< 0.001）。该基因模块包含13个基因，包括ANP32E、ATAD2、BID、CBX3、CCNB1、DTL、ECT2、GMPS、KIF14、MCM10、NDC1、NETO2和SORBS2。
对蓝色模块中同时也是差异表达基因的基因进行基因集富集分析
我们将WGCNA识别的蓝色模块中的基因与DEGs中的相关基因取交集，得到了13个基因（图5A）。基于已识别模块与DEGs重叠的这13个基因，我们进行了PPI网络分析以评估它们的相互作用（图5B）。
随后，我们对WGCNA识别的蓝色模块中的基因进行了GSEA。在生物过程方面，主要富集于细胞周期、有丝分裂细胞周期G1/S转换和有丝分裂细胞周期（图5C）。在分子功能方面，主要富集于酶结合、细胞骨架衔接蛋白活性和死亡受体结合（图5D）。在细胞组分方面，主要富集于细胞骨架、核膜和非膜结合细胞器（图5E）。在KEGG通路方面，主要富集于p53信号通路、细胞周期和鞘脂信号通路（图5F）。
ESCC诊断模型的构建与外部验证
选择蓝色模块中的13个基因作为LASSO回归分析的潜在预测因子，LASSO系数路径图如图6A所示。共有5个非零系数基因，包括BID（系数0.37）、CBX3（系数0.03）、ECT2（系数0.42）、KIF14（系数0.38）和SORBS2（系数-0.12）。LASSO回归模型中的系数在图6B中可视化。使用这五个基因开发了ESCC诊断模型，并通过多元逻辑回归进行分析。在多元逻辑回归分析中，SORBS2成为ESCC的统计学显著预测因子（表1）。我们的预测模型以便于使用的列线图形式可视化（图6C）。
我们使用TCGA作为外部验证来评估诊断模型的性能。训练集和测试集的AUC值分别为0.99和0.98（图6D、E），诊断模型的性能见表2，表明模型具有良好的准确性。校准曲线表明模型校准可接受，实际频率与预测概率之间具有良好的一致性（图6F、G）。训练集中的R2为0.724；Bootstrap内部验证产生的R2为0.6663，而外部验证集（TCGA）的R2为0.702。最终，我们生成了DCA曲线，以证明我们的模型在临床环境中的实际效用（图6H、I）。按血统进行的分层分析表明，在亚洲亚组（n=25；AUC = 0.95，R2 = 0.85，Brier = 0.04）中区分度一致，在白人亚组（n=57；AUC = 0.98，R2 = 0.84，Brier = 0.04）中性能稳定。由于黑人或非裔美国人病例数较少（n = 2），且未知/混合组的组成可能具有异质性（n = 7），为了避免不可靠的估计和潜在的偏倚，这些亚组的性能指标未在表2中呈现。
小鼠ESCC模型中标志基因表达的实验验证
为了确认五个基因在ESCC和正常组织之间的表达差异，我们建立了小鼠ESCC模型（图7A）。我们发现，与正常小鼠相比，ESCC模型小鼠的体重没有显著增加（图7B）。ESCC模型的肿瘤生长趋势如图7C所示。正常食管组织和ESCC组织的比较如图7D所示。苏木精-伊红（HE）染色显示，与正常小鼠相比，ESCC小鼠模型的组织结构遭到破坏，其特征是细胞的大小和形状极不规则，细胞排列紊乱（图7E）。随后，我们进行了蛋白质免疫印迹以评估蛋白水平，结果显示，与癌旁正常食管组织相比，ESCC组织中BID、CBX3、ECT2、KIF14蛋白表达显著增加，而SORBS2蛋白表达降低（图7F、G）。
免疫细胞浸润分析揭示抗肿瘤免疫受损与SORBS2在ESCC免疫微环境中的低表达相关
我们使用CIBERSORT算法分析了ESCC的免疫微环境（IME），该算法可分析22种不同的免疫细胞类型。结果显示，在这22种细胞类型中，正常组和ESCC组之间共有11种微环境细胞表现出显著差异（图8A）。值得注意的是，四种淋巴细胞群——记忆B细胞、CD8+T细胞、滤泡辅助T细胞和静息肥大细胞——在IME内的浸润相对于正常对应物显著减弱。这种免疫表型景观表明抗肿瘤免疫受到严重损害。此外，我们分析了SORBS2表达与ESCC的IME中免疫细胞浸润之间的相关性，揭示记忆B细胞、CD8+T细胞、滤泡辅助T细胞和静息肥大细胞在低表达SORBS2组相对于高表达SORBS2组的浸润减弱（图8B）。因此，IME中SORBS2的低表达与ESCC抗肿瘤免疫功能受损相关。
基因集富集分析揭示了SORBS2在免疫反应和细胞信号通路中的功能
为了探索SORBS2的功能意义，我们进行了GSEA。GO富集分析显示，SORBS2的生物过程主要参与由循环免疫球蛋白介导的体液免疫反应、调节吞噬作用的免疫反应细胞表面受体信号通路以及激活免疫反应的细胞表面受体信号通路（图9A）。在细胞组分方面，SORBS2富集于核糖核蛋白复合物、T细胞受体复合物和受体复合物（图9B），而在分子功能方面，SORBS2与核糖体的结构组成、免疫受体活性和信号受体活性相关（图9C）。利用KEGG富集分析，SORBS2与细胞粘附分子、剪接体和神经活性配体-受体相互作用通路显示出显著关联（图9D）。
单细胞分析揭示SORBS2主要在肌成纤维细胞和血管内皮细胞簇中表达
对来自GSE196756数据库的六个样本进行了单细胞测序分析，结果显示nCount_RNA和nFeature_RNA之间存在正相关（图10A）。方差图在所有细胞中识别出250个基因，其中2000个高变基因用红色标记，并标注了前10个基因（图10B）。图10C、D显示细胞被分为10个簇。图10E表明，十个已识别的标记基因用于显示在不同的簇中，而热图则显示了每种细胞类型的前5个标记基因（图10F）。基于关键标记基因的表达，这些簇被分类为B细胞、癌细胞、CD8 T细胞、DC细胞、成纤维细胞、巨噬细胞、肌成纤维细胞和血管内皮细胞（图10G、H）。图10I说明了五个基因在细胞类型间的分布。图10J-L表明SORBS2
