《Frontiers in Cell and Developmental Biology》：TRIM26 deficiency drives gastric cancer lymph node metastasis via TGF-β signaling activation and modulates gemcitabine response
1 引言
胃癌（GC）是全球范围内主要的消化道恶性肿瘤，其发病率和死亡率分别位居第五和第四位，对全球健康构成重大负担。该疾病具有明显的地域分布特征，东亚地区（包括中国、日本和韩国）的发病率和死亡率显著高于西方国家。尽管内镜筛查、手术技术和靶向治疗的进步改善了早期患者的预后，但大多数病例在诊断时已处于晚期，治疗选择有限。因此，晚期胃癌的5年生存率仍低于30%，显著低于其他胃肠道恶性肿瘤。这一严峻的预后凸显了阐明胃癌转移机制，特别是淋巴结转移（LNM）机制的迫切性，这对于改善患者预后至关重要。
LNM是胃癌扩散的主要和初始途径，严重影响TNM分期、手术计划和辅助治疗。临床证据一致表明，转移性淋巴结的负荷与总生存期（OS）相关，LNM阳性患者的复发率显著更高，长期结局更差。预后差异显著：广泛淋巴结阳性病例的5年生存率降至20%以下，而淋巴结阴性患者则超过60%。因此，理解LNM的分子基础对于改进风险分层和指导精准治疗至关重要。从机制上讲，胃癌淋巴扩散是一个由肿瘤细胞与肿瘤微环境（TME）之间动态串扰驱动的多步骤过程。关键事件包括上皮-间质转化（EMT），它增强了细胞的运动性和侵袭性；通过蛋白酶分泌进行的细胞外基质重塑；以及为扩散开辟途径的血管生成/淋巴管生成。免疫逃逸进一步支持转移，因为肿瘤相关巨噬细胞、调节性T细胞和其他免疫抑制成分分泌细胞因子以抑制抗肿瘤 immunity。此外，代谢重编程——如谷氨酰胺和脂质代谢的失调——为支持转移进展提供能量和生物合成前体。因此，胃癌LNM并非源于孤立的遗传改变，而是多种信号通路和细胞相互作用的汇聚。虽然TGF-β、PI3K/AKT和Wnt/β-catenin等通路是公认的贡献者，但调控LNM的核心调控网络仍未完全阐明，临床上可应用的分子驱动因素仍然缺乏。
三联基序（TRIM）家族E3泛素连接酶在泛素化中起着核心作用，调节多种过程，包括蛋白质降解、DNA损伤修复、免疫信号传导和细胞周期进程。TRIM蛋白包含特征性的RING、B-box和卷曲螺旋结构域，在癌症发病机制中表现出双重功能。例如，TRIM28和TRIM59通常通过增强增殖和转移而作为癌蛋白发挥作用，而TRIM16和TRIM62则通过抑制生长和促进凋亡而作为肿瘤抑制因子发挥作用。其中，TRIM26是一种具有组织特异性作用的上下文依赖性调节因子。它通过促进p53降解和阻碍凋亡在肝细胞癌中显示出致癌活性，而在抗病毒 immunity 中，它抑制STING-IRF3通路。蛋白质组学和转录组学分析进一步将TRIM26的异常表达与各种癌症中的肿瘤侵袭和免疫逃逸联系起来。尽管有这些见解，TRIM26在胃癌中的功能——特别是其对LNM的影响——在很大程度上仍未探索，这突出了其作为阐明胃癌新转移机制的关键靶点的潜力。
单细胞RNA测序（scRNA-seq）通过实现对细胞异质性、功能状态和发育轨迹的高分辨率分析，改变了肿瘤生物学——这些能力超出了批量转录组分析的范围。在胃癌中，scRNA-seq揭示了肿瘤进化和微环境动态，例如免疫抑制细胞组成的变化、癌症相关成纤维细胞的激活以及淋巴结中转移前生态位的重塑。当与WGCNA、机器学习（LASSO、随机森林、SVM-RFE）、伪时间推断和细胞间通讯分析等分析方法结合时，scRNA-seq为系统识别驱动转移进展的关键基因和信号网络提供了一个强大的框架。
本研究旨在阐明TRIM26在驱动LNM和调节药物反应中的双重作用。通过对单细胞和批量转录组数据的整合分析，我们将描绘以TRIM26为中心的调控网络，表征其对转移和肿瘤微环境重塑的功能影响，并 pinpoint 与TRIM26相关的化疗药物。这些见解将批判性地评估TRIM26作为临床生物标志物和治疗靶点的潜力。
2 材料与方法
2.1 数据获取与预处理
scRNA-seq数据集来源于GEO数据库（登录号：GSE163558），包含三个原发性GC组织、一个癌旁正常组织和六个转移性GC样本——包括两个LNM。从中选择三个原发性GC和两个LNM组织进行下游分析以确保可比性。
使用Seurat R包进行数据预处理。初始质量控制涉及过滤表达少于三个基因的细胞和在少于50个细胞中检测到的基因。这些步骤有效减少了技术噪音和背景信号，为后续分析步骤产生了高质量的表达矩阵。
批量RNA-seq训练数据来自TCGA-STAD项目，仅限于具有可用转录组谱（基因表达定量）的胃源性腺瘤和腺癌。年龄、性别和种族等临床变量未用作排除标准。基因水平表达值通过对多探针测量值取平均进行合并，并移除无检测表达的样本。数据使用Toil流程统一处理，转换为FPKM然后TPM格式，最后进行log2转换以提高可比性。整理相应的临床数据，并排除年龄小于18岁或生存时间少于30天的患者以增强临床相关性。
2.2 scRNA-seq数据的处理与可视化
使用Seurat（v4.0.0）处理scRNA-seq数据。应用两步过滤策略去除低质量细胞和潜在的双重态。排除在少于3个细胞中检测到的基因和表达少于250个基因的细胞。基于以下阈值进行进一步的质量控制：（1）保留检测到500–6000个基因的细胞；（2）包括UMI >1,000的细胞，排除UMI计数前3%的细胞；排除线粒体基因含量>35%的细胞，并移除线粒体表达最高的前2%细胞；（3）移除核糖体RNA比例最高和最低各1%的细胞；（4）保留总RNA计数>1,000的细胞，排除总RNA计数前3%的细胞。为进一步通过去除潜在双重态（单个液滴中捕获多个细胞）确保数据质量，我们使用DoubletFinder R包（v2.0.3）进行双重态检测，估计双重态形成率为10%。预测的双重态从后续分析中移除。过滤后，使用NormalizeData函数对表达数据进行标准化，并通过FindVariableFeatures选择高变基因。进行主成分分析（PCA）进行降维，并使用Harmony算法以“orig.ident”作为协变量校正批次效应。使用诊断图评估收敛性。基于前30个主成分应用UMAP进行可视化。使用FindNeighbors和FindClusters进行细胞聚类，并通过Idents分配簇身份。
2.3 单细胞RNA测序分析与细胞类型注释
使用标准分析流程处理单细胞RNA测序数据。基于检测到的基因数量和线粒体基因表达过滤低质量细胞。在标准化和识别高变基因后，通过主成分分析（PCA）进行降维，然后使用基于图的聚类算法进行无监督聚类。
根据先前胃癌单细胞研究中报道的成熟规范标记基因的表达对细胞亚群进行注释。具体而言，上皮细胞通过EPCAM和KRT19的高表达来识别；恶性上皮细胞通过升高的EPCAM、KRT8和KRT18表达进一步表征。T细胞通过CD3D和CD3E表达定义，CD4+和CD8+T细胞亚群分别通过CD4和CD8A识别。B细胞基于MS4A1和CD79A表达进行注释，而浆细胞通过MZB1和XBP1的高水平识别。髓系细胞通过LST1和TYROBP表达表征，巨噬细胞通过CD68和C1QC定义，树突状细胞通过FCER1A和CLEC10A定义。成纤维细胞使用COL1A1和COL1A2识别，内皮细胞通过PECAM1和VWF识别。
CopyKAT分析显示EPCAM阳性上皮簇表现出广泛的非整倍性，证实了它们的恶性性质。因此，这些恶性上皮细胞被定义为胃癌细胞并用于下游分析。
2.4 使用hdWGCNA鉴定GC细胞中与LNM相关的基因
为了识别GC中与LNM相关的关键基因模块，我们进行了高维加权基因共表达网络分析（hdWGCNA）。使用SetupForWGCNA准备高质量表达矩阵，并通过MetacellsByGroups（k = 25, max_shared = 10）生成元细胞以提高数据稳健性。对元细胞矩阵进行标准化、标准化，并使用PCA和Harmony进行整合，随后进行UMAP可视化。在通过TestSoftPowers确定最佳软阈值功率后构建共表达网络。使用ModuleEigengenes计算模块特征基因以总结表达模式。使用FindMarkers（P < 0.05）进行原发和LNM样本之间的差异表达分析。候选LNM相关基因定义为显著差异表达基因（DEG）与相关WGCNA模块中识别的枢纽基因的交集。
2.5 基于机器学习的LNM相关基因筛选
为了识别与LNM相关的稳健基因特征，我们采用了多步骤特征选择方法。首先，使用单变量Cox回归筛选具有预后意义的转移相关基因。随后，应用三种不同的机器学习算法来细化基因集：（1）通过glmnet包进行LASSO回归，执行带正则化的特征选择；（2）SVM-RFE（e1071包）基于支持向量机权重进行递归特征消除；（3）随机森林（randomForest包）通过平均减少Gini重要性对基因排序。最终候选基因定义为所有三种方法共同选择的特征的交集，然后通过生存分析和表达模式评估进行验证。
2.6 表达模式、诊断和预后价值的评估
系统评估了候选基因在TCGA队列中正常胃粘膜、原发性GC和LNM组织中的表达谱。为了评估候选基因的预后相关性，进行了单变量和多变量Cox比例风险回归分析。将相关的临床病理变量，包括年龄、性别、肿瘤分期和组织学分级，纳入多变量模型以调整潜在的混杂因素，从而确定候选基因是否作为独立的预后指标。
2.7 下游信号通路分析
使用两种互补的计算方法研究下游信号通路：将TCGA-STAD样本按关键基因的中位表达分为高表达组和低表达组。使用limma（∣log2FC∣ >1, FDR <0.05）识别DEG。同时，对前25%变异最大的基因进行WGCNA以构建共表达网络，并优化软阈值功率以实现无标度拓扑。选择与LNM显著相关的基因模块（P < 0.05, ∣相关系数∣ >0.3），并将其与DEG的重叠用于通过clusterProfiler（P < 0.05）进行KEGG通路富集分析。将Mfuzz软聚类应用于TCGA样本中前5,000个高变异基因（c = 50个簇）。使用χ2检验或逻辑回归测试表达动态与LNM状态的关联。对来自LNM相关簇的核心基因进行KEGG富集分析以识别表型相关通路。
2.8 伪时间轨迹分析
使用Monocle2和Monocle3中的伪时间分析重建GC LNM期间的细胞状态转换。在Monocle2中，创建CellDataSet，随后估计大小因子和离散度。保留平均表达≥0.1的基因，并选择前200个高变基因进行DDRTree降维。使用orderCells沿轨迹对细胞排序，BEAM分析识别分支依赖性差异表达基因（P < 0.05），特别关注转录因子和信号成分。为了验证，使用基于UMAP的降维应用Monocle3，使用learn_graph进行轨迹推断，并使用graph_test评估空间自相关性（Moran's I）。这种双框架方法确保了驱动转移进展的转录程序的稳健识别。
2.9 细胞间通讯分析
使用CellChat表征肿瘤微环境内的细胞间信号网络。利用其整理的配体-受体数据库并结合蛋白质-蛋白质相互作用网络，我们计算了通讯概率并量化了通路水平的信息流。根据其信号模式将细胞群体分类为不同的功能角色（发送者、接收者、介导者和影响者）。原发和LNM组织之间的比较分析识别了转移相关的细胞间通讯改变。
2.10 转录因子活性推断
通过整合DoRothEA数据库与Viper算法计算推断转录因子（TF）活性。使用高置信度调控相互作用（置信度水平A–C）构建调控子，计算每个细胞的活性分数，并通过PCA和UMAP可视化。跨细胞群体的差异TF活性分析通过层次聚类（pheatmap）进行和可视化。从显著改变的TF构建蛋白质-蛋白质相互作用网络，并通过Cytoscape中的MCODE识别枢纽调节因子。
2.11 代谢重编程分析
使用scMetabolism在单细胞分辨率下量化代谢通路活性，该工具应用ssGSEA和AUCell算法对关键过程进行评分，包括糖酵解、脂肪酸代谢、氨基酸代谢和TCA循环。在原发和LNM样本之间以及候选基因的高表达组和低表达组之间进行比较分析。在这些比较中一致改变的途径被识别为与转移进展相关的核心代谢特征。
2.12 分子对接
进行分子对接以评估TRIM26与化疗药物吉西他滨和紫杉醇之间的潜在结合。从小分子PubChem获得三维结构，并使用ChemOffice优化。从RCSB PDB数据库检索TRIM26的晶体结构，并使用PyMOL 2.6去除非蛋白质组分进行准备。使用AutoDock 1.5.6，添加氢原子并将可旋转键分配给配体。使用AutoDock Vina进行对接模拟，结合位点特定坐标定义搜索空间。选择具有最有利结合能的构象，其中低于-5.0 kcal/mol和-7.0 kcal/mol的值分别表示实质性和高亲和力结合。在Discovery Studio 2019和PyMOL 2.6中可视化蛋白质-配体相互作用。
2.13 分子动力学模拟
使用GROMACS 2022进行分子动力学（MD）模拟以评估TRIM26与吉西他滨和紫杉醇复合物的动态稳定性。TRIM26结构使用amber14sb力场建模，而配体参数通过AutoFF服务器从GAFF2导出。每个蛋白质-配体复合物在TIP3P水盒（1.0 nm边界）中溶剂化，并用离子中和。使用Particle Mesh Ewald方法计算静电相互作用，截止值为1.0 nm。通过LINCS算法约束键，并使用Verlet leap-frog积分器应用1 fs时间步长。
系统在能量最小化（3,000步最陡下降，2000步共轭梯度）后，在NPT条件（310 K）下进行100 ns生产模拟。分析轨迹的均方根偏差（RMSD）、均方根涨落（RMSF）、回转半径（Rg）、溶剂可及表面积（SASA）和氢键计数，以评估构象稳定性和结合相互作用。
2.14 细胞培养
人GC细胞系AGS和HGC-27购自北京协和医学院细胞资源中心，在含有10%胎牛血清（Gibco, 美国, #10099-141）和1%青霉素-链霉素（Gibco, 美国, #15140122）的DMEM或RPMI-1640培养基（HyClone, 美国, #SH30809.01, SH30022.01）中培养。细胞在37°C、5% CO2的湿润环境中维持，每3天以1:3的比例传代。
2.15 CCK-8实验
使用Cell Counting Kit-8（Dojindo, 日本, #CK04）测量细胞活力。将细胞以2–5 × 103个细胞/孔的密度接种在96孔板中（每组3–5个重复），并使其贴壁12–24小时。用浓度为0.2、0.4、0.6和0.8 μmol/L的吉西他滨（MCE, 美国, # HY-17026）处理后，向每孔加入10 μL CCK-8试剂并孵育4小时。使用酶标仪（BioTek Synergy H1）记录450 nm处的吸光度。最后，使用GraphPad Prism软件通过四参数逻辑回归模型计算IC50值。模型方程为y = A2 + (A1 - A2)/(1 + (x/IC50) ^ p)，其中y表示细胞活力，x是抑制剂浓度，A1和A2分别对应最大和最小活力值，p代表斜率因子。
2.16 慢病毒转染
TRIM26过表达的慢病毒购自吉玛基因有限公司（中国上海）。当AGS和HGC27细胞的汇合度达到约20%时，以1:30的感染复数（MOI）将慢病毒转导到细胞中，并加入polybrene以提高感染效率。转导后24小时更换培养基，转导后48小时加入含嘌呤霉素的培养基进行筛选。在嘌呤霉素选择后存活的细胞被认为是稳定转导的克隆，并通过Western blot分析验证转导效率。
2.17 集落形成实验
将细胞以500–1,000个细胞/孔的密度接种在6孔板中（每组三个重复孔），并培养10–14天。将形成的集落用4%多聚甲醛（Beyotime, 中国 #P0099）固定，用结晶紫（Beyotime, #C0121）染色，并在成像后定量。
2.18 Transwell实验
通过Transwell小室实验（Corning, 美国, # 353097）评估细胞迁移能力。简言之，将总共3 × 104个悬浮于400 μL无血清培养基中的细胞接种到小室的上室。对于侵袭实验，在接种细胞前，将Matrigel以1:15的比例稀释并预涂于上室。同时，下室填充含有20% FBS的完全培养基作为趋化剂。孵育24小时后，小心擦去膜上表面未迁移的细胞。将粘附在膜下表面的迁移细胞用4%多聚甲醛固定，结晶紫染色，最后在倒置显微镜（Nikon, 日本东京）下观察和拍照。
2.19 Western blot
使用含有磷酸酶和蛋白酶抑制剂（Beyotime, #P0013, P1045）的RIPA裂解缓冲液裂解细胞。离心后收集上清液，并通过BCA试剂盒（Beyotime, #P0011）测定蛋白质浓度。然后，将蛋白质样品与上样缓冲液混合，并在95°C下变性5分钟后再上样。通过SDS-PAGE电泳（NCM Biotech, 中国, #P2012）实现蛋白质（30 μg）分离，然后转移至膜（GE Healthcare Life Science, 德国, #10600023）。在封闭非特异性结合位点后，将膜在4°C下与一抗孵育过夜。第二天，洗涤后，将膜与相应的二抗孵育，最后使用化学发光法成像（NCM Biotech, #10,100）。Western blot分析使用的抗体如下：TRIM26（Proteintech, 美国, #27013-1-AP），GAPDH（HUABIO, 中国, #ET1601-4），TGF-β1（HUABIO, #HA721143），P-Smad2/3（abcam, #ab254407），Smad2/3（abcam, #ab207447）。
2.20 统计分析
使用GraphPad Prism 9.0进行所有数据的统计分析。对于服从正态分布的数据，应用Student's t检验；对于非正态分布数据，改用Wilcoxon秩和检验。*P < 0.05, P < 0.01, 和*P < 0.001。
3 结果
3.1 GC免疫微环境的单细胞转录组图谱
经过严格的质量控制，从原发性GC病变和匹配的淋巴结转移中保留了22,450个高质量细胞，表达24,577个基因。数据标准化、主成分分析（前30个PC）和基于UMAP的非线性降维识别出18个转录不同的细胞簇。数据集包括11,442个来自原发肿瘤的细胞和11,008个来自转移淋巴结的细胞，揭示了位点之间的组成差异。跨簇的转录异质性明显，簇0、2、4和7显示高全局基因表达，而簇6、8和14显示相对较低的活动性。
基于CopyKAT的拷贝数改变（CNA）推断识别了与EPCAM阳性簇大部分重叠的非整倍体群体，证实了它们的恶性上皮起源。使用规范标记基因结合多方法注释（SingleR、CellMarker数据库和手动整理），我们将18个亚簇分类为八个主要谱系，包括上皮细胞、免疫细胞、成纤维细胞和内皮细胞。比较分析揭示了这些细胞类型在原发和转移位点之间相对丰度的变化，强调了淋巴扩散过程中肿瘤微环境的动态重塑。
3.2 使用hdWGCNA识别与LNM相关的基因模块
为了识别与LNM相关的上皮基因模块，我们使用软阈值功率4对单细胞表达矩阵进行了加权基因共表达网络分析（hdWGCNA）。这揭示了14个不同的共表达模块。其中，M3和M5模块与上皮细胞身份表现出强相关性，并共同包含127个在上皮亚簇中显著富集的枢纽基因。尽管存在一些模块间共表达，但M3和M5保持了相对独立的表达模式，并且它们与上皮的强关联通过相关性分析得到进一步证实。这些结果表明M3和M5代表了可能调控GC中上皮行为和转移倾向的关键功能单元。
3.3 LNM候选基因的筛选及预后相关性
为了识别LNM的预后驱动因素，我们首先对来自M3和M5模块的127个枢纽基因进行单变量Cox回归，识别出22个与总生存期显著相关的基因（P < 0.05）。为了增强稳健性，我们应用了三种机器学习算法：LASSO回归选择了9个预后基因，随机森林对前10个特征
