scTenifoldXct：基于流形对齐与基因调控网络的细胞通讯分析新方法

1. 项目概述：为什么我们需要更聪明的细胞通讯分析工具？

在单细胞RNA测序（scRNA-seq）技术普及的今天，我们得以窥见组织中成千上万个细胞的基因表达图谱。然而，理解这些细胞如何“交谈”——即细胞通讯——才是揭示组织功能、发育过程乃至疾病机制的关键。细胞通讯的核心是配体-受体（Ligand-Receptor, LR）相互作用：一个细胞（发送者）分泌配体蛋白，与另一个细胞（接收者）表面的受体蛋白结合，从而触发下游信号传导。传统上，识别这些关键的LR对主要依赖于一个简单的假设：相互作用的配体和受体基因应该同时高表达。因此，许多工具（如CellChat、NATMI）将两个基因表达量的乘积作为相互作用的强度指标。

但现实往往更复杂。我在分析多个炎症和肿瘤数据集时反复遇到一个困境：一些已知在生物学上至关重要的信号通路，其对应的LR基因在测序数据中表达量并不高，甚至很低。例如，在皮肤炎症中起关键作用的S100A8-ITGB2对，其表达水平远低于许多其他LR对。如果仅凭表达量筛选，这些“低调”但功能强大的信号分子就会被淹没在背景噪音中，导致我们错过真正的生物学故事。这就像在一个嘈杂的派对上，只注意那些喊得最大声的人，而忽略了角落里低声交谈却决定派对走向的关键人物。

这正是scTenifoldXct试图解决的问题。它不再单纯依赖表达量的绝对值，而是引入了一个更精巧的框架：通过流形对齐（Manifold Alignment）技术，在考虑细胞内基因调控网络（Gene Regulatory Network, GRN）背景的前提下，量化细胞间基因的相似性。简单来说，它不仅要看两个基因（配体和受体）自己“喊”得多响，还要看它们在各自细胞内部的“社交网络”（即基因调控网络）中处于什么位置，以及这两个网络结构是否“对齐”。这种方法的核心洞见在于，功能上有关联的基因，即使在表达量上不突出，它们在调控网络中的相对位置也可能是相似的。这种方法显著提升了对低表达但功能重要LR对的检测灵敏度，并能发现现有知识库（如OmniPath）中未记录的潜在新相互作用。

2. scTenifoldXct核心原理深度拆解：从网络构建到流形对齐

要理解scTenifoldXct的强大之处，我们需要深入其三个核心步骤。这不仅仅是流程介绍，更是理解其为何能突破传统方法局限的关键。

2.1 第一步：构建联合相似性矩阵——搭建跨细胞对话的桥梁

这一步的目标是为后续的流形对齐准备一个全面的“关系图谱”，它包含细胞内和细胞间两个维度的信息。

2.1.1 构建单细胞基因回归网络（scGRNs）

对于你感兴趣的两种细胞类型（例如A：成纤维细胞，B：树突状细胞），scTenifoldXct会分别为它们构建基因调控网络。这里采用的方法是主成分回归（Principal Component Regression）。

为什么用PC回归而不是更复杂的方法？在单细胞数据中，基因数量（n）通常远大于细胞数量（m），这导致了经典的“高维诅咒”问题。直接进行多元线性回归会面临严重的共线性和过拟合。PC回归通过先对自变量（所有其他基因的表达矩阵）进行主成分分析（PCA），提取出主要的方向（主成分，PCs），然后用这些互不相关的主成分作为新的自变量来回归目标基因。这有效地降低了维度，过滤了噪音，并且计算上相对高效稳定。最终，对于细胞类型A，我们得到一个邻接矩阵 (W_A)，其中每个元素 (W_A(i, j)) 表示基因j对基因i的调控系数（可正可负，代表激活或抑制）。(W_B) 同理。

注意：这里构建的是“基因回归网络”，它近似但不等同于真实的基因调控网络。它反映了基因表达层面的统计关联，是数据驱动的，不需要预先知道转录因子（TF）与靶基因的对应关系。用户也可以直接输入自己构建的高质量网络（如通过SCENIC、GENIE3等方法推断的）来替代这一步。

2.1.2 计算基因间串扰分数

这是连接两个细胞类型的关键。对于细胞类型A中的基因i和细胞类型B中的基因j，scTenifoldXct计算一个串扰分数（Crosstalk Score）：(score(i_A, j_B))。

其计算公式并非简单的表达量乘积，而是： [ S_{ij} = (\mu_i \mu_j) \times ( \frac{\sigma_i}{\mu_i} + \frac{\sigma_j}{\mu_j} ) ] 其中，(\mu) 和 (\sigma) 分别代表基因在所有细胞中的平均表达量和标准差。

这个设计的精妙之处在于：

1. 均值项 ((\mu_i \mu_j)): 保留了传统思想，相互作用的基因对倾向于共表达。
2. 变异项 (( \frac{\sigma}{\mu} )，即变异系数): 引入了表达变异性的考量。一个基因在细胞群体中表达波动大（高变异系数），可能意味着它受到精细调控，或在特定细胞亚群中高度特异表达，这种异质性可能恰恰是细胞通讯特异性的来源。将两个基因的变异系数相加，放大了那些表达模式“活跃”或“特异”的基因对的分数。

因此，这个分数同时考虑了表达强度和表达特异性，使得那些在特定功能亚群中高表达、而在其他细胞中低表达的LR对也能获得较高的分数。所有基因对的串扰分数构成了矩阵 (S)。

2.1.3 组装联合相似性矩阵

最后，将上述两部分组合成一个大的联合相似性矩阵(W)： [ W = \begin{bmatrix} W_A & S \ S^T & W_B \end{bmatrix} ]

这个矩阵是scTenifoldXct所有分析的基石。它的对角线块 (W_A) 和 (W_B) 描述了细胞内的调控关系（基因与基因如何相互影响），而非对角线块 (S) 描述了细胞间的潜在相互作用关系（一个细胞的基因与另一个细胞的基因如何关联）。将两者放在同一个矩阵中，意味着后续的流形对齐过程会同时考虑“细胞内环境”和“细胞间吸引力”，从而更准确地判断两个基因是否真的在功能上配对。

2.2 第二步：流形对齐——在统一的空间中寻找“门当户对”

流形对齐是scTenifoldXct的灵魂。什么是流形？你可以想象为一个弯曲的、低维的空间，你的高维数据（这里是基因）分布在这个空间里。(W_A) 和 (W_B) 分别定义了细胞类型A和B的基因在其各自高维表达空间中的“关系网”（即流形结构）。

流形对齐的目标是找到一个共同的低维潜在空间，将来自 (W_A) 和 (W_B) 的基因都映射到这个空间里。在这个共同空间中，如果来自细胞A的配体基因L和来自细胞B的受体基因R的距离非常近，那就意味着，综合考虑了它们各自的细胞内调控背景和彼此间的串扰潜力后，它们在功能上是高度匹配的，因此更可能发生真实的相互作用。

scTenifoldXct采用神经网络来学习这个对齐映射，而不是传统的基于特征值分解的方法。这样做的好处是：

• 非线性能力：基因调控关系往往是非线性的，神经网络能更好地捕捉这种复杂模式。
• 计算效率与稳定性：对于大规模单细胞数据（数万个基因），基于神经网络的优化比求解大规模矩阵特征值问题更稳定、更高效。
• 灵活性：易于扩展和集成其他类型的约束或信息。

这个过程可以类比为“相亲大会”。(W_A) 和 (W_B) 分别是两位嘉宾（基因）的“个人资料”（兴趣爱好、社交圈）。流形对齐就像是一个智能匹配系��，它不只看两个人的表面条件（表达量），而是深入分析他们的资料（调控网络），将他们投射到一个“三观匹配度”空间。在这个空间里距离近的两个人，说明他们的内在特质更契合，更有可能成功配对（发生LR相互作用）。

2.3 第三步：识别显著的LR对——从距离到生物学意义

经过流形对齐，每个基因（包括所有配体和受体）在这个低维潜在空间中都有一个坐标。对于任何一对潜在的LR对（配体基因L来自细胞A，受体基因R来自细胞B），我们可以计算它们之间的欧几里得距离(d(L, R))。

核心逻辑：距离越近，相互作用可能性越大。

接下来是统计检验。scTenifoldXct对所有可能的LR对（通常是数万甚至数百万对）计算距离，然后通过一种非参数检验方法（如置换检验）来评估观察到的距离是否显著小于随机预期。最终，通过控制错误发现率（FDR，例如<0.05），筛选出一组显著的LR相互作用对。

与双样本分析的扩展：scTenifoldXct的强大之处还在于其双样本比较框架。比如，比较疾病组织与正常组织。它并不是独立分析两个样本然后比较结果，而是构建一个更大的耦合联合矩阵，将两个样本的四个网络（样本1的A、B细胞网络，样本2的A、B细胞网络）和两个串扰矩阵一起进行流形对齐。在这个统一的潜在空间中，分别计算同一LR对在两个样本中的距离 (d_{disease}) 和 (d_{normal})，其差值 (\Delta d = d_{disease} - d_{normal}) 的大小和方向（正负）就反映了相互作用的强度在疾病状态下是增强还是减弱。再通过卡方检验等统计方法识别出差异显著的LR对。

3. 实战演练：从数据到生物学洞见

理论讲得再多，不如亲手操作一遍。下面我将以一篇经典的炎症性皮肤病scRNA-seq数据为例，带你走通scTenifoldXct的完整分析流程，并解读结果。

3.1 环境准备与数据加载

scTenifoldXct主要是一个R包。确保你的R版本在4.0以上。

# 1. 安装必要的R包 if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("scTenifoldXct") # 安装并加载其他辅助包 install.packages(c("Matrix", "igraph", "ggplot2", "pheatmap")) library(scTenifoldXct) library(Matrix) library(ggplot2) # 2. 加载示例数据 # 假设我们有一个SingleCellExperiment对象`sce`，其中包含基因表达矩阵和细胞类型注释 # 列：cells， 行：genes # 细胞类型信息存储在`sce$celltype`中 # 我们关注两种细胞类型：'Fibroblast' 和 'Dendritic_cell' # 提取目标细胞类型的表达矩阵和标签 cell_types_of_interest <- c('Fibroblast', 'Dendritic_cell') idx <- sce$celltype %in% cell_types_of_interest sub_sce <- sce[, idx] count_matrix <- counts(sub_sce) # 获取原始计数矩阵 cell_labels <- sub_sce$celltype

3.2 运行scTenifoldXct核心分析

现在，我们将运行单样本分析，寻找成纤维细胞和树突状细胞之间显著的LR相互作用。

# 3. 运行scTenifoldXct单样本分析 # 主要函数：scTenifoldXct() # 参数说明： # - X: 基因 x 细胞的表达矩阵（建议使用log归一化后的数据，如log1p(CPM)） # - cellCluster: 长度等于细胞数的向量，指定每个细胞的类型 # - cellTypeA: 发送者细胞类型名称 # - cellTypeB: 接收者细胞类型名称 # - pc_ratio: 用于PC回归的主成分比例，默认0.2（即使用前20%的PCs） # - nNet: 用于构建稳健网络的子采样网络数量（默认10） # - nCom: 流形对齐后潜在空间的维度（默认50） # - nCore: 并行计算使用的核心数 # 对计数矩阵进行简单的log转换（模拟标准化数据） log_norm_matrix <- log1p(t(t(count_matrix) / colSums(count_matrix) * 1e4)) # CPM + log1p set.seed(123) # 设置随机种子保证结果可重复 result <- scTenifoldXct( X = log_norm_matrix, cellCluster = cell_labels, cellTypeA = 'Fibroblast', cellTypeB = 'Dendritic_cell', pc_ratio = 0.2, nNet = 10, nCom = 50, nCore = 4 # 根据你的电脑配置调整 ) # 4. 查看结果 # result对象是一个列表，包含多个组件 names(result) # 最重要的组件是`result$results`，它是一个数据框，包含所有LR对的详细信息 lr_results <- result$results head(lr_results[order(lr_results$p.value), ]) # 按p值排序查看最显著的LR对

输出结果数据框通常包含以下列：

• ligand: 配体基因名
• receptor: 受体基因名
• distance: 在潜在空间中的欧氏距离
• p.value: 原始p值
• p.adj: 经过多重检验校正后的p值（如FDR）
• score: 可能包含的额外评分（如串扰分数）

3.3 结果筛选与注释

我们需要根据校正后的p值筛选显著结果，并与已知数据库进行比对。

# 5. 筛选显著LR对 (例如 FDR < 0.05) significant_lr <- subset(lr_results, p.adj < 0.05) nrow(significant_lr) # 查看显著对的数量 # 6. 与已知数据库（如OmniPath）进行比对，区分已知和新颖相互作用 # scTenifoldXct可能内置了查询功能，或者我们可以使用其他R包如`OmnipathR` if (!require("OmnipathR", quietly = TRUE)) BiocManager::install("OmnipathR") library(OmnipathR) # 获取OmniPath中的配体-受体互作信息 op_interactions <- import_ligrecextra_interactions() # 简化信息 known_pairs <- unique(paste(op_interactions$source_genesymbol, op_interactions$target_genesymbol, sep = "-")) # 标记我们的结果 significant_lr$pair <- paste(significant_lr$ligand, significant_lr$receptor, sep = "-") significant_lr$is_known <- significant_lr$pair %in% known_pairs # 查看已知和新颖相互作用的比例 table(significant_lr$is_known) # 7. 功能富集分析（以受体基因为例） if (!require("clusterProfiler", quietly = TRUE)) BiocManager::install("clusterProfiler") if (!require("org.Hs.eg.db", quietly = TRUE)) BiocManager::install("org.Hs.eg.db") library(clusterProfiler) library(org.Hs.eg.db) # 提取所有显著LR对中的受体基因 receptor_genes <- unique(significant_lr$receptor) # 将基因符号转换为Entrez ID gene_ids <- bitr(receptor_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) # 进行GO生物过程富集分析 ego <- enrichGO(gene = gene_ids$ENTREZID, OrgDb = org.Hs.eg.db, keyType = "ENTREZID", ont = "BP", # 生物过程 pAdjustMethod = "BH", qvalueCutoff = 0.05, readable = TRUE) # 可视化前10个富集结果 dotplot(ego, showCategory=10, title="GO Enrichment Analysis of Receptors")

3.4 网络可视化

理解LR对所在的调控上下文至关重要。我们可以利用结果中的网络信息进行可视化。

# 8. 提取并可视化围绕特定LR对的调控网络 # 假设我们对'CCL19-CCR7'这对感兴趣 target_ligand <- 'CCL19' target_receptor <- 'CCR7' # 从结果对象中提取网络数据（具体函数可能因包版本而异，请参考文档） # 通常，result对象中会包含对齐后的基因坐标或网络邻接矩阵 # 这里假设我们有一个函数可以提取以特定基因为中心的子网络 # 例如，提取与CCL19和CCR7最相关的前N个基因（基于回归网络中的权重） plot_integrated_network <- function(result, lig, rec, top_n = 15) { # 此函数为示例，实际实现需根据scTenifoldXct包的输出结构调整 # 1. 从result$network_A (Fibroblast网络) 中找到与配体lig最相关的top_n个基因（TF或其他基因） # 2. 从result$network_B (Dendritic网络) 中找到与受体rec最相关的top_n个基因 # 3. 使用igraph包绘制二分网络图 # 边颜色：正调控为红色，负调控为蓝色 # 边粗细：与权重绝对值成正比 # 突出显示LR对本身（如用绿色高亮） library(igraph) # ... (具体的网络提取和绘图代码) # 返回一个igraph图对象或直接绘图 } # 调用函数（需根据实际数据结构实现） # net_graph <- plot_integrated_network(result, target_ligand, target_receptor) # plot(net_graph, layout = layout_with_fr, ...)

实操心得：在运行scTenifoldXct时，pc_ratio和nCom是两个关键的超参数。pc_ratio控制PC回归中保留的信息量，对于噪声较大的数据，可以适当降低（如0.1）以增强鲁棒性；对于细胞数很多、信号强的数据，可以增加（如0.3）以保留更多细节。nCom是潜在空间的维度，通常设置在20-100之间。维度太低会丢失信息，太高则会引入噪音并增加计算负担。建议通过检查结果中LR对距离的分布稳定性来调整，或者对已知的阳性对照LR对进行测试，看其排名是否合理。

4. 方法对比与性能评估：scTenifoldXct强在哪里？

在原文的炎症皮肤数据案例中，作者将scTenifoldXct与五种主流方法（NATMI, SingleCellSignalR, Connectome, iTALK, CellChat）进行了系统比较。我们可以从中总结出scTenifoldXct的几个突出优势，并在自己的分析中加以验证。

4.1 检测灵敏度：发现“沉默的”关键信号

传统方法严重依赖表达量。在炎症皮肤数据中，像S100A8-ITGB2和CCL26-CCR6这样的LR对，其基因表达水平很低，因此被所有其他五种方法遗漏。然而，大量的生物学证据表明，S100A8-ITGB2在特应性皮炎的发病机制中通过调节细胞因子和皮肤屏障蛋白发挥关键作用；CCL26-CCR6则与炎症性树突状细胞和记忆T细胞的招募密切相关。

scTenifoldXct成功检测到了这些对。根本原因在于其流形对齐机制。即使S100A8和ITGB2表达量不高，但它们在各自细胞（如角质形成细胞和免疫细胞）的基因调控网络中，可能处于与炎症反应通路相关的关键“枢纽”位置，或者它们的表达模式（高变异系数）显示出细胞亚群特异性。当两个网络在潜在空间中对齐时，这些在调控背景上“门当户对”的基因就会被拉近，从而被识别为显著的相互作用。

验证建议：在你自己的数据分析中，可以特别关注那些被scTenifoldXct识别为显著、但表达乘积排名很靠后的LR对。通过查阅文献或进行功能富集分析，验证它们是否属于已知但低表达的通路，这可能是新发现的起点。

4.2 综合性能：高精度与高召回率的平衡

作者使用了一个“金标准”集（被其他五种方法共同检测到的40个LR对）来绘制ROC曲线和精确率-召回率曲线。scTenifoldXct取得了最高的AUROC（0.89）和平均精度（AP，0.87）。这表明scTenifoldXct在识别已知真实相互作用方面，既准确（高精度）又全面（高召回率）。

这得益于其整合策略：传统方法是“先独立分析，后比较”，而scTenifoldXct的双样本分析是“联合对齐，直接比较差异”。后者在数学上更优雅，能更好地捕捉样本间细微但一致的差异，避免了独立分析引入的批次效应或噪声放大问题。

4.3 提供机制性解释：从相互作用到调控网络

这是scTenifoldXct区别于绝大多数工具的独特优势。它不仅告诉你“谁和谁说话”，还能通过其内置的基因回归网络，展示这个LR对在细胞内是如何被调控的。

例如，在分析CCL19-CCR7这对相互作用时，scTenifoldXct构建的整合网络显示，在树突状细胞中，转录因子JUN（AP-1复合体的亚基）正向调控CCR7的表达。这与已知生物学知识——AP-1蛋白在皮肤炎症中发挥作用——完全吻合。同时，它还发现NF-κB亚基REL与CCR7有强关联，这与CCR7能激活NF-κB通路的实验证据一致。

给你的分析带来的价值：当你得到一个显著的LR对列表后，不要止步于此。利用scTenifoldXct生成的网络数据（通常是回归系数矩阵），去探索：

• 上游调控：哪些转录因子最可能调控这个配体或受体基因？
• 协同通路：这个LR基因是否与细胞内其他关键信号分子（如细胞因子、激酶）紧密相连？
• 跨细胞协调：发送细胞和接收细胞内部的调控网络是否存在镜像或互补关系？

这能将你的分析从“现象描述”提升到“机制假设”，为后续的实验验证提供强有力的线索。

5. 常见问题、排查技巧与高级应用

在实际使用中，你可能会遇到各种问题。以下是我在多次实践中总结的经验和解决方案。

5.1 数据预处理与输入要求

5.2 结果解读与验证陷阱

• 如何判断结果是否可靠？

  1. 阳性对照：检查你的系统中已知的、经典的LR对（如皮肤炎症中的CCL19-CCR7）是否出现在显著结果列表中，且排名靠前。
  2. 富集分析一致性：对显著LR对的配体或受体基因进行GO或KEGG富集分析。结果是否富集到与你生物学问题相关的通路（如“免疫应答”、“细胞趋化”、“炎症反应”）？
  3. 表达模式验证：虽然scTenifoldXct不依赖高表达，但一个完全不在任何细胞中表达的基因被预测为配体/受体是可疑的。用点图或小提琴图检查关键LR对的表达分布，确保它们至少在部分细胞中有表达信号。

• 发现了大量“新颖”相互作用，怎么办？这是scTenifoldXct的优势，也可能是挑战。首先，不要盲目相信所有新预测。

  1. 数据库局限性：OmniPath等数据库并不完整。通过PubMed或STRING数据库进行文献挖掘，验证这些新颖对是否有间接证据支持。
  2. 跨物种谨慎性：如果使用小鼠数据，而用人类数据库注释，需考虑基因同源性的问题。
  3. 功能一致性：将这个新颖LR对的上下游基因（从其调控网络中找）一起做富集分析。如果它们都富集到同一个生物学过程，那么这个新颖相互作用的可信度就大大增加。

5.3 高级应用：双样本比较分析

双样本比较（如疾病vs正常）是scTenifoldXct的亮点。操作流程与单样本类似，但需要准备两个样本的数据。

# 假设有normal_sce和disease_sce两个SingleCellExperiment对象 # 分别提取相同两种细胞类型的数据并合并，同时保留样本来源信息 # 构建一个包含所有细胞的表达矩阵X_all # 构建一个细胞标签向量cell_labels_all，格式可以是“样本_细胞类型”，如“Normal_Fibroblast” # 构建一个样本来源向量sample_id_all # 使用scTenifoldXct的双样本比较函数（具体函数名请查阅最新文档，可能为`scTenifoldXct2`或通过参数指定） # 以下为概念性代码 diff_result <- scTenifoldXct_compare( X = X_all, cellCluster = cell_labels_all, sampleID = sample_id_all, # 指定每个细胞属于哪个样本 cellTypeA = 'Fibroblast', cellTypeB = 'Dendritic_cell', sample1 = 'Normal', sample2 = 'Disease' ) # 结果中会包含每个LR对在两个样本中的距离，以及距离差和显著性p值 # 筛选在疾病中显著增强或减弱的相互作用 enhanced_pairs <- subset(diff_result, diff_distance < 0 & p.adj < 0.05) # 距离减小表示增强 weakened_pairs <- subset(diff_result, diff_distance > 0 & p.adj < 0.05) # 距离增大表示减弱

解读差异结果时需注意：一个LR对相互作用的“增强”，在scTenifoldXct中表现为潜在空间距离的“减小”。这可能是由于配体/受体表达上调，也可能是其细胞内调控网络发生了重编程，使得两个基因在功能上更“匹配”。因此，结合差异表达分析（DE）一起看，能提供更完整的图景：是表达变化主导，还是网络重布线主导？

5.4 与下游分析的衔接

scTenifoldXct的结果可以无缝对接多种下游分析：

• 网络可视化与模块挖掘：将显著的LR对及其相关的上游调控基因（从scGRN中提取）导入Cytoscape进行高级网络分析和模块识别。
• 轨迹分析与细胞通讯动态：如果你有单细胞拟时序数据，可以将细胞按拟时序排序，然后滑动窗口运行scTenifoldXct，观察关键LR相互作用强度如何随时间变化。
• 整合空间转录组数据：虽然scTenifoldXct目前不直接处理空间信息，但其预测的LR对可以与空间共定位信息结合。例如，使用ST数据验证预测的配体和受体是否在空间相邻的斑点中共同表达。

最后，记住任何计算工具都是辅助。scTenifoldXct提供了一个强大的、基于网络的假设生成框架。它给出的每一个显著LR对，尤其是那些新颖的、低表达的预测，都是一个等待实验生物学验证的、充满潜力的科学故事起点。将计算预测与领域知识、实验验证紧密结合，才是解锁细胞通讯奥秘的正确方式。