生信程序 | NatGenet | 使用潜在嵌入多变量回归分析多条件单细胞数据

Basic Information

• 英文标题：Analysis of multi-condition single-cell data with latent embedding multivariate regression
• 中文标题：使用潜在嵌入多变量回归分析多条件单细胞数据
• 发表日期：03 January 2025
• 文章类型：Article
• 所属期刊：Nature Genetics
• 文章作者：Constantin Ahlmann-Eltze | Wolfgang Huber
• 文章链接：https://www.nature.com/articles/s41588-024-01996-0

Abstract

Para_01

1. 识别不同条件下异质组织中的基因表达差异是一项基本的生物学任务，这可以通过多条件单细胞RNA测序（RNA-seq）来实现。
2. 目前的数据分析方法将组成细胞分为聚类，旨在代表细胞类型，但这种离散分类往往是对潜在生物学的一种不满意的模型。
3. 在这里，我们引入了一种称为潜在嵌入多元回归（LEMUR）的模型，该模型在或在做出离散分类承诺之前进行操作。
4. LEMUR（1）整合来自不同条件的数据，（2）预测每个细胞的基因表达变化作为条件和其在潜在空间中的位置的函数，以及（3）对于每个基因，识别一个具有一致差异表达的紧凑细胞邻域。
5. 我们将LEMUR应用于癌症、斑马鱼发育和阿尔茨海默病中的空间梯度，展示了其广泛的应用性。

Main

Para_01

1. 过早离散化连续变量会导致数据分析中的伪影和功效损失；然而，在多条件单细胞数据中处理细胞类型和状态的多样性时，这是主要的方法。
2. Lähnemann等人1描述了克服在下游分析前依赖于聚类或（离散）细胞类型分配的问题，这是单细胞数据分析中的一个重大挑战之一。

Para_02

1. 单细胞RNA测序可以用于研究实验干预或观察条件对一组异质细胞的影响，例如来自组织活检或类器官的细胞。
2. 通常，同一份样本中的细胞共享相同的条件，但来自多种细胞类型和状态（例如，在分化路径或细胞衰老路径中的位置、细胞周期、代谢）。
3. 多条件单细胞RNA测序相对于批量测序的优势在于能够区分不同条件下相应细胞之间的表达变化，以及不同细胞类型或状态之间的表达变化。

Para_03

1. 我们提出了一种生成模型和推理方法来解决多条件单细胞数据分析中的三个任务：（1）将数据整合到一个共同的潜在空间中，（2）对于每个细胞，预测它在任何一种条件下的表达情况，以及（3）发现有趣且具有统计显著性的差异表达模式。

Para_04

1. 对于第一个任务，存在许多方法，这些方法或多或少明确地旨在使在公共潜在空间中仍然代表细胞类型或细胞状态而不是外部条件的剩余变化。

Para_05

1. 然而，只有部分整合方法解决了第二个任务的反事实预测。
2. Harmony2、Seurat3和MNN4没有将集成嵌入中的位置映射回基因表达空间的标准方法。
3. 相比之下，scVI5、scGen6、CPA7和CellOT8使用自动编码器，在其中编码器将基因表达空间映射到潜在空间，解码器则反向映射回来。
4. 潜在嵌入多元回归（LEMUR）不分别学习编码器和解码器，而是拟合解析可逆函数，从集成嵌入生成条件特异性基因表达值。

Para_06

1. 对于第三项任务，跨条件差异表达分析，最先进的方法是采用一个综合嵌入，将细胞分配到簇中，并使用从bulk RNA测序分析中已知的方法分别对每个给定的簇找到差异表达的基因（‘伪批量’）。
2. 在这里，我们将这一过程颠倒过来。我们利用LEMUR反事实预测来计算每个细胞和每个基因的差异表达统计量，然后选择具有一致差异表达的连接细胞集。

Para_07

1. LEMUR模型将数据中的变化分解为四个来源：
2. 明确已知的条件，
3. 未明确已知但假定可以用低维流形表示的细胞类型或状态，
4. 二者之间的相互作用，
5. 未解释的残余变异性。

Para_08

1. LEMUR 在一个名为 lemur 的 Bioconductor R 包和一个名为 pyLemur 的 Python 包中实现。

Results

Para_01

1. 图1a概述了LEMUR的工作流程。该方法以基因×细胞大小的数据矩阵作为输入。
2. 它假设已经进行了适当的预处理，包括规模因子归一化和方差稳定变换11。
3. 此外，它还期望有两个元数据表用于细胞：一个分类变量，对于每个细胞，指明其来源的样本（独立实验单元，例如组织活检或类器官），以及一个设计矩阵12。

Fig. 1: Conceptual overview of LEMUR.

• a, 四步工作流程。
• b, LEMUR 核心的矩阵分解。
• c, 关于 a 中的每一步详细信息：步骤 1，为每个条件分别拟合一个线性子空间。不同条件下的子空间通过仿射变换相互关联，这些仿射变换由协变量参数化。为了这个可视化，绘制了三维基因空间中的两个二维子空间；实际维度更高。
• 步骤 2，在对照组和处理组预测值之间计算差异表达统计量 Δ 作为差值。该可视化展示了步骤 1 的俯视图。
• 步骤 3，对于每个基因，将具有相似 Δ 值且彼此靠近的细胞分组为邻域。
• 步骤 4，在每个邻域内对细胞应用伪批量差异表达检验。
• Ctrl.，对照；DE，差异表达；dim，维度；trt，处理。

Para_02

1. 设计矩阵编码一个或多个表示实验处理或观察条件的协变量。
2. 它类似于差异表达工具（如limma13和DESeq2（参考文献14））中的设计矩阵，并且可以解释完全一般的实验或研究设计。
3. 设计矩阵可以包括不需要的变化来源（例如，实验批次）和我们感兴趣的变化来源（例如，处理状态）。
4. 我们将协变量值的独特组合称为条件。
5. 在两个条件比较的简单情况下，设计矩阵是一个两列矩阵，其中第一列（截距）的所有元素都是1，第二列的元素是0或1，表示每个细胞来自哪个条件。

Matrix factorization with known covariate information

具有已知协变量信息的矩阵分解

Para_01

1. 考虑每个细胞作为一个由其测量的基因表达谱定义的高维空间中的点。
2. 流形假设认为数据集中在该高维空间内的一个（未知的）低维流形上。
3. 经验上已经发现，通过主成分分析（PCA）的前几十个分量所构成的线性向量空间可以对那个低维流形做出有用的近似。

Para_02

1. LEMUR 的核心思想是使用 PCA 的多条件扩展来表示多条件单细胞 RNA 测序数据（图 1b）。
2. 给定一个数据矩阵 Y，PCA 可以用来近似 Y ≈ RZ，其中包含两个较小的矩阵：第一个矩阵 R 包含主成分，这些主成分作为低维线性子空间的基础；第二个矩阵 Z 包含每个细胞相对于该基础的坐标。
3. 在这种基本形式下，没有地方可以明确编码已知的实验或研究协变量。
4. LEMUR 通过加入回归分析组件来实现这一功能。

Para_03

1. 我们不是使用一个单一子空间，而是为每个条件找到一个单独的子空间。
2. 为此，我们让表示子空间的矩阵 R(X) 依赖于设计矩阵 X 中提供的协变量（图 1c，步骤 1）。
3. 通过这一假设，我们解决了上述分解任务：已知的变化源被编码在 X 中；细胞类型和状态变化由细胞坐标 Z 表示。
4. 虽然 X 是明确已知的，但 Z 是潜在的；也就是说，它是从数据中估计出来的。
5. 这种构建方式允许建模交互作用，也就是说，不同细胞类型和状态下跨条件的基因表达变化是不同的。
6. 这是 LEMUR 的主要特征。

Para_04

1. 这一构造的第二个特征是跨条件整合：来自不同条件下观察到的细胞的数据被映射到一个共同的潜在空间。
2. 默认情况下，这种整合基于相应子空间的对齐，但可以通过表明某些在不同条件下观察到的细胞彼此对应的指示信息来选择性地改进。
3. 此类信息可以以显式地标的形式出现，即表达独特标记基因的细胞，也可以通过诸如Harmony2等方法利用的细胞统计特性来获得。
4. 我们通过向模型中添加仿射变换S(X)来考虑这种对应关系。

Para_05

1. 我们模型的第三个特征是，对于每个细胞，它能够预测在任何条件下该细胞的基因表达谱会是什么样子，即使这种条件只在一个条件下被观察到。
2. 事实上，这些预测不仅适用于那些细胞被观察到的位置，还适用于所有位置，也就是说，也可以用于假设性地内插或外推细胞类型和状态。
3. 我们利用这些预测来发现跨潜在空间区域（即相同或相似的细胞类型和状态）的基因表达变化。

Para_06

1. 图2展示了LEMUR方法的风格化插图。
2. LEMUR为每个条件拟合一个一维子空间（线），每个条件由施加在共同基空间上的旋转参数化。
3. 参数模型为每个条件中的每个细胞产生预测的表达值，我们寻找潜在空间中的区域（在这里，我们有两个主要区域，左和右）在这些区域中预测值始终为正或负。
4. 在方法部分，我们提供了更正式的数学规范。

Fig. 2: Stylized example with two genes observed in two groups of cells in two conditions.

• 散点图展示了基因空间和特定条件的一维子空间。详情见方法部分。
• 预测了每个细胞中的两个基因的差异表达。基因1仅在‘左’细胞类型中相对于对照组在处理组中表达较少；基因2仅在‘右’细胞类型中在处理组中上调表达。
• 数据来源

Cluster-free differential expression analysis

无聚类差异表达分析

Para_01

1. 我们可以使用这个参数模型预测细胞在任何条件下的表达。
2. 两种条件之间的差异表达（图1c，步骤2）仅仅是它们预测值之间的差，并且可以针对每个细胞进行计算，即使对于任一细胞的数据仅在一个条件的一个重复中被观察到（图2b）。

Para_02

1. 由此产生的差异表达估计矩阵Δ有两个用途：首先，我们可以将每个基因的差异表达值可视化为潜在空间的函数。潜在空间的维度通常是十到一百，仅用于可视化，我们使用进一步的非线性降维到二维散点图，例如均匀流形近似和投影（UMAP）15。示例见图4和图6。其次，我们使用Δ来指导差异表达区域的识别，即一致显示某个特定基因差异表达的一组细胞（图1c，步骤3）。这些区域的目的是连接、凸出且最大，也就是说，如果扩展这些区域，差异表达模式将变得不成比例地不一致且不显著。
2. 由此产生的差异表达估计矩阵Δ有两个用途：首先，我们可以将每个基因的差异表达值可视化为潜在空间的函数。潜在空间的维度通常是十到一百，仅用于可视化，我们使用进一步的非线性降维到二维散点图，例如均匀流形近似和投影（UMAP）15。示例见图4和图6。其次，我们使用Δ来指导差异表达区域的识别，即一致显示某个特定基因差异表达的一组细胞（图1c，步骤3）。这些区域的目的是连接、凸出且最大，也就是说，如果扩展这些区域，差异表达模式将变得不成比例地不一致且不显著。

Para_03

1. 为了排名和评估我们对发现的社区的信心水平，我们不试图测量预测 Δ 的统计不确定性。
2. 相反，我们使用原始数据的伪批量聚集：如果有可用的，则使用原始计数，否则使用对数归一化的值。
3. 对于每个样本，我们将原始计数相加或取该社区中细胞的对数归一化值的平均值，以获得特定于社区的基因 × 样本表（图 1c，步骤 4），然后进行差异表达检验，使用 glmGamPoi16、edgeR17 或 limma13。
4. 这里没有单独的句子，但为了完整性，将上一句视为包含所有信息的句子。

Outputs

Outputs

Para_01

1. LEMUR产生以下输出：

• 一个通用的低维潜在空间表示所有细胞（Z），
• 参数化变换 R 和 S，用于将条件特定的潜在空间相互映射，
• 每个基因和细胞在任意条件下的预测表达值 (\hat{{\bf{Y}}})，
• 从设计矩阵构建的任意对比中，每个基因和细胞的预测差异表达 Δ，
• 对于每个基因和对比，细胞的邻域和显著性统计测量（P 值）。

Para_02

1. 这些中的最后一个通常是最终用户最关心的；其余的对于诊断、质量评估、可视化以及数据的进一步建模使用是有用的。

Para_03

1. 显式的变换R和S的参数化意味着它们可以很容易地在观测数据集之外进行插值或外推，并且可以从低维嵌入反向转换回数据空间。
2. 从这个意义上说，LEMUR是一个生成模型。

Performance assessment

绩效评估

Para_01

1. 我们使用了13个公开可用的多条件单细胞数据集（在数据可用性中列出）评估了LEMUR的表现。
2. 我们遵循Ahlmann-Eltze和Huber11的方法一致地预处理了每个数据集。

Para_02

1. 首先，我们考虑了整合性能：细胞的联合低维表示在保留潜空间中编码的生物信号的同时，如何去除已知协变量（如批次和处理效应）的痕迹？
2. 图3a展示了Kang等人18研究中狼疮患者样本使用干扰素β或载体对照治疗的数据集上的这一情况。
3. 我们通过计算每个细胞的k=20个邻居中有多少来自相同条件（k最近邻混合）来衡量协变量的去除。
4. 对于两个条件平衡的数据集，理想的方法k-最近邻混合值为k/2=10。
5. 我们通过分别比较嵌入数据的聚类与原始数据的聚类来衡量生物信号的保留，用调整后的兰德指数（ARI）的平均值来衡量。
6. 理想的方法得分接近ARI=1。

Fig. 3: Performance assessment.

• UMAPs of the latent spaces for the data from Kang et al.18.
• The k-NN mixing coefficient and the ARI are defined in the main text.
• b, Density plots of the bootstrapped mean performance for ARI and k-NN mixing across 13 datasets.
• To adjust for dataset-dependent variation, we divided the k-NN and ARI scores by the average per dataset.
• c, Scatterplots of predicted (pred.) expression under treatment (y axis) against observed expression (x axis) for 500 genes in each of eight cell types (same data as in a).
• d, Prediction error as in c, across the same 13 datasets as in b (gray points).
• Red points show the mean.
• e, Simulation setup.
• Left, for one of the implanted genes, a UMAP is shown of the LEMUR latent space, where the color indicates whether an expression change in this gene was simulated for that cell.
• Center and right, simulated expression values.
• f, Left, predicted log fold change for the gene from e (Δ=Ŷ^B-Ŷ^A).
• Center and right, the set of cells inside or outside the inferred neighborhood.
• g, Comparison of observed false discovery proportion and true positive rate (TPR) for 11 datasets, with ten replicates and the overall mean shown as a large point.
• The nominal FDR was fixed to 10%.
• FP, false positive; IFN, interferon; m, months; NK, natural killer; TP, true positive.

Para_03

1. 在13个数据集中，LEMur在这两项指标上的表现与Harmony相似（图3b）。
2. 其他方法在这两项指标之间做出了不同的权衡，没有哪种方法明显优于其他方法。
3. 扩展数据图1显示，在另外七个附加指标中，结果保持一致。

Para_04

1. 运行 LEMUR 的计算成本处于预期范围内的低端，并且与其他近似PCA方法相当。例如，在Goldfarb-Muren数据上计算前50个潜在维度，该数据包含24,178个细胞和20,953个基因（占用4GB的RAM），使用近似的irlba算法19耗时35秒和24GB的RAM。
2. 相比之下，不进行整合地拟合 LEMUR 需要103秒并且需要33GB的RAM。
3. 使用地标或Harmony对细胞进行对齐分别增加了2秒和95秒（扩展数据图2）。

Para_05

1. 接下来，我们评估了LEMUR预测跨条件基因表达的能力。
2. 我们使用它来预测在对照条件下观察到的细胞在治疗条件下的基因表达，并将这些预测与实际上已经接受处理的细胞的数据进行比较。
3. 为了避免过拟合，我们在‘保留’的细胞上评估了预测结果，这些细胞的数据未用于训练。
4. 由于在两种条件下观察到的单个细胞之间没有直接对应关系，我们考虑了注释细胞类型的平均值。
5. 图3c显示了Kang等人18号数据集中八种细胞类型中500个最易变基因的预测-观察对比散点图，涉及四种方法：CPA、scVI、LEMUR和无变化（恒等）的简单预测。
6. 在13个数据集上，LEMUR通过观测值和预测值之间的L2距离测量显示出最小的预测误差（图3d）。
7. 扩展数据图3显示了结果在六个附加指标上的一致性。

Para_06

1. 在第三组比较中，我们测试了 LEMUR 识别具有一致差异表达的细胞群的能力。
2. 我们取用了 Kang 等人18数据中的对照条件下的所有细胞，将它们随机分配到条件A或B，并在部分细胞中植入差异表达的基因。
3. 图3e展示了一个例子。
4. LEMUR 准确地识别了表达变化，并推断出一个与模拟真实情况重叠良好的细胞邻域（图3f）。

Para_07

1. 我们扩展了这一分析，增加了十个更多的半合成数据集，每个数据集中植入了200个不同的差异表达基因，以评估LEMUR差异表达检验中的I型错误控制情况。
2. LEMUR平均控制了假发现率（FDR）(图3g，顶部)。
3. 此外，它比对所有细胞进行伪批量测试或对细胞类型或聚类进行子集独立测试更具有功效，这些测试方法类似于Crowell等人9中的方法，或如miloDE10中的邻域测试方法(图3g，底部)。
4. 扩展数据图4评估了LEMUR方法的其他变体的FDR控制和功效。
5. 结果显示，考虑从邻域推断中产生的选择后偏差非常重要，并且通过数据分割成功解决了这一问题。

Para_08

1. 总体而言，这些基准测试表明 LEMUR（1）成功地整合了来自不同条件的单细胞数据，（2）在没有访问先前的细胞聚类或分类的情况下检测细胞类型和状态特异性差异表达模式，以及（3）提供了准确的统计I型错误控制和良好的功效。接下来，我们将LEMUR应用于不同生物数据集的分析。
2. 在下文中，我们应用 LEMUR 来分析不同的生物学数据集。

Treatment–control data: panobinostat in glioblastoma

治疗对照数据：panobinostat在胶质母细胞瘤中的应用

Para_01

1. 赵等人报告了胶质母细胞瘤活检的单细胞RNA测序数据。
2. 五名患者的样本在两种条件下进行了检测：对照组和帕比诺司他，一种组蛋白去乙酰化酶（HDAC）抑制剂（图4a）。
3. 经过质量控制，数据包含了65,955个细胞（扩展数据表1）。

Fig. 4: Analysis of five glioblastoma tumor biopsies.

• a, 实验设计。
• b, 对数转换数据的UMAP（第一列）和LEMUR产生的潜在嵌入（第二和第三列），按处理和细胞类型着色（第一行和第二行分别）。
• c, 肿瘤细胞内LMO2的差异表达分析。内部和外部邻域内的细胞的分面UMAP。中间的散点图显示了供体和条件下的伪批量表达值。右侧，比较邻域内外的差异（红色箭头）。P值基于负二项计数模型的双侧似然比检验。
• d, 在控制条件下，子群体内外细胞的比较火山图。注释为涉及翻译活性的基因用红色突出显示。位于水平线以上的基因集FDR为10%。
• Nbh.表示邻域；panob.表示panobinostat；diff. in diff.表示差异中的差异；expr.表示表达；cytopl. transl.表示细胞质翻译；scRNA-seq表示单细胞RNA测序。数据来源

Para_02

1. 图4b的左列显示了输入数据Y的二维UMAP。大部分可见的变化与已知的协变量相关：供体和处理条件。
2. 还有一些变化与活检中的不同细胞类型有关。我们使用LEMUR吸收供体和处理效应到R中，将潜在空间维度设置为P=60。
3. 图4b的中间列显示，在将S(x)固定为单位矩阵后，每个细胞的潜在坐标矩阵Z的UMAP。
4. 来自不同样本的细胞更加混合在一起，样本内的细胞异质性变得更加明显。
5. 在使用S通过Harmony的最大多样性聚类来匹配样本间的细胞亚群后（图4b，右列），这种图像变得更加清晰。
6. 在这里，大量的肿瘤细胞群体和三个非肿瘤亚群体变得显而易见。

Para_03

1. 连续改进的潜在空间表示从左到右进一步展示在图4b的底部一行，其中点根据我们从选定标记基因的表达和已知的染色体非整倍体现象（扩展数据图5a，b）获得的细胞类型分配进行着色。

Para_04

1. LEMUR 的线性潜在空间易于解释。
2. 这在扩展数据图 5c 中得到了体现，该图将 PCA21 的双变量图概念扩展到了多条件设置。
3. 我们可以通过绘制基因的载荷向量相对于 Z 坐标系的位置来探索任何基因的高表达或低表达如何影响细胞在潜在空间 Z 中的位置。

Para_05

1. 使用 LEMUR 的差异表达检测，我们发现帕比司他导致了 25% 的所有基因（10,000 个基因中的 2,498 个）在 10% 的假发现率下发生了细胞亚群特异性的表达变化。
2. 扩展数据图 6 显示了七个基因的差异表达和推断的邻域。

Para_06

1. 我们专注于肿瘤细胞，识别出对帕比诺司他治疗反应不同的亚群（图4c）。
2. 在来自所有五名患者的肿瘤细胞的一个亚群中（共9,430个肿瘤细胞），帕比诺司他治疗导致LMO2下调，
3. 而在大多数肿瘤细胞（n = 36,535）中，LMO2的表达保持不变（图4c，右侧）。
4. LMO2的产物与转录因子TAL1、TCF3和GATA形成蛋白复合物；
5. 它对于血管生成至关重要，并且最初是在T细胞急性淋巴细胞性白血病中被发现的一种致癌基因22。
6. Kim等人23研究了LMO2在胶质母细胞瘤中的作用，发现更高的表达与患者生存率降低相关。
7. 他们得出结论，LMO2可能是一个具有临床意义的药物靶点。
8. 为了进一步表征对帕比诺司他治疗通过下调LMO2作出响应的肿瘤细胞亚群，
9. 我们将这些细胞在对照条件下的整体基因表达谱与其它肿瘤细胞进行了比较。
10. 我们发现核糖体基因的表达较低，这与翻译活动减少一致（图4d）。

Time course data: zebrafish embryo development

时间序列数据：斑马鱼胚胎发育

Para_01

1. Saunders等人报告了斑马鱼胚胎发育图谱，其中包括在受精后18小时到48小时之间每隔2小时收集的16个时间点上的967个胚胎和838,036个细胞的数据（图5a）。
2. 他们利用单细胞RNA测序数据将每个细胞分配到一个通用的细胞类型分类方案，并研究了沿发育时间出现和消失的细胞类型的时态动态。
3. 我们问基因表达变化是否能揭示更多的生物学现象。
4. 因此，我们寻找在共享相同细胞类型注释的细胞间系统性不同的基因表达的时间模式。
5. 为此，我们利用LEMUR预测任何基因在潜在空间中的任何时间点的表达。
6. 为了用较少的参数表示每个基因的时间依赖性，我们使用了具有三个自由度的自然三次样条，遵循Smyth等人的方法。
7. 为了建模潜在空间的变化，我们在潜在空间Z中对来自中枢神经系统的两个细胞（红色和绿色）进行线性插值，并且类似地，在来自表皮（一种覆盖发育胚胎早期阶段的外层上皮层）的两个细胞（紫色和青绿色）之间进行插值（图5b和扩展数据图7a）。