Basic Information
英文标题： Method of moments framework for differential expression analysis of single-cell RNA sequencing data 中文标题：单细胞 RNA 测序数据差异表达分析的方法框架 发表日期：24 October 2024 文章类型：Resource 所属期刊：Cell 文章作者：Min Cheol Kim | Chun Jimmie Ye 文章链接：https://www.sciencedirect.com/science/article/pii/S0092867424011449
Highlights
Para_01
用于单细胞 RNA 测序的统计模型将测量噪声和表达噪声区分开来 高效的重采样允许进行校准良好的假设检验 Memento 使研究基因对扰动的协调表达成为可能 Memento 映射与基因表达均值、变异性和相关性相关的位点
Summary
Para_01
单细胞 RNA 测序（scRNA-seq）数据的差异表达分析对于表征实验因素如何影响基因表达分布至关重要。 然而，区分细胞间变异的生物学和技术来源以及评估细胞组之间定量比较的统计显著性仍然具有挑战性。 我们引入了 Memento，这是一种用于从 scRNA-seq 数据中进行稳健且高效的均值表达、变异性和基因相关性差异分析的工具，可扩展到数百万个细胞和数千个样本。 我们将 Memento 应用于 70,000 个气管上皮细胞，以鉴定干扰素响应基因；应用于 160,000 个 CRISPR-Cas9 扰动的 T 细胞，以重建基因调控网络；应用于 120 万个外周血单核细胞（PBMCs），以绘制细胞类型特异性的数量性状位点（QTLs）；并应用于包含 5000 万个细胞的 CELLxGENE Discover 资源库，以比较任意细胞组。 在所有情况下，Memento 相比现有方法识别出更多显著且可重复的均值表达差异。 它还识别出变异性和基因相关性的差异，这表明扰动赋予了不同的转录调控机制。
Graphical abstract
Keywords
scRNA-seq analysis; Differential expression; Bootstrap method; Gene-regulatory networks; Spatial genomics; Single-cell transcriptomics; CELLxGENE Discover; CD4 T cells; Memento; Gene expression variance
Introduction
Para_01
基因表达由细胞的遗传组成和环境相互作用内在决定，可能会因内在噪声（源自mRNA转录和降解）和与细胞特定状态相关的外在噪声而波动。 虽然遗传和环境历史对细胞群体中的表达变异性有显著贡献，但随机转录噪声也会影响细胞对扰动的反应以及细胞的发育和分化。 表征确定性和随机因素如何共同影响基因表达分布是理解转录控制如何建立、维持以及可能被打破的核心。 这些见解可以阐明那些基因型-表型关系未完全解释的现象背后的机制，例如不稳定、不完全外显率和可变表现度。
Para_02
细胞群体中的基因表达分布主要由其均值和方差以及相关衍生指标来描述。 组成型表达的管家基因以恒定的速率进行转录和降解，预计符合泊松分布。 然而，大多数基因表现出过分散，方差高于预期，且同一生物途径内的基因通常在转录上是相关的。 这些观察结果与一个模型一致，该模型认为相关基因的表达受相似的顺式调控元件调控，这些元件与一组共同的转录因子相互作用，这些转录因子在"开启"和"关闭"状态之间循环。 直到最近，研究基因表达的分布，特别是多个基因的联合分布，在技术上一直具有挑战性，并且主要在可以进行遗传改造的模式生物中进行。
Para_03
单细胞 RNA 测序（scRNA-seq）已成为一种系统且高效的方法，用于在实验因素（包括细胞外刺激、基因扰动和自然遗传变异）下分析细胞的转录组。 理论上，scRNA-seq 数据的分析可以揭示确定性和随机性因素如何共同塑造基因表达的分布。 然而，仍然需要能够比较细胞组之间分布参数（包括均值、变异性以及基因相关性）差异的差异分析方法。 为了评估均值表达的差异，通常的做法是对伪批量谱型进行差异表达分析，这些伪批量谱型是通过将定义为聚类的细胞组的转录计数汇总生成的。 尽管伪批量方法没有充分利用单细胞作为重复测量的优势，但它们出人意料地优于显式建模观察到的 scRNA-seq 数据分布的方法。 此外，用于评估基因表达变异性及基因对之间相关性差异的方法非常少。
Para_04
单细胞 RNA 测序（scRNA-seq）数据的广义差异表达分析仍然是一个严峻的挑战，主要是由于两个关键的统计限制。 首先，将观察到的细胞间变异分解为其组成成分——生物学变异和测量噪声——是一个重大障碍。 这一困难源于基因转录和 scRNA-seq 采样过程中涉及的分子数量较少（图 1A）。 大多数现有方法实施了参数化模型，旨在解释观察到的稀疏转录计数的高于预期的方差。 然而，这些模型并未明确建模测量噪声，这是 scRNA-seq 工作流程固有的欠采样特性的副产品。 重要的是，准确估计生物学变异对于有效建模基因对之间的相关性至关重要。 其次，确定特定细胞群之间均值、变异或基因相关性的统计显著性仍然是一个未解决的问题。 许多现有方法利用渐近理论来确定比较均值的假设检验的显著性，通常会产生未校准的 p 值。 这在需要进行数千次比较的研究中尤其成问题，因为未校准的 p 值违反了多重检验校正的假设。 此外，大多数现有方法需要精确指定参数模型，并且在有效结合层次结构和连续协变量方面缺乏灵活性。 因此，它们并未明确考虑从多路复用工作流程中生成的生物学和技术重复，这些工作流程适应了越来越多的个体或条件。 像 DESCEND 这样的方法是值得注意的例外，它们利用灵活定义的广义线性模型，在理论上能够有效解决这一问题。 然而，当建模 scRNA-seq 数据固有的复杂层次结构时，这些模型经常遇到显著的计算障碍，并且仅限于特定的细胞间变异模型。 确实，最近的研究报告称，与伪批量方法相比，scRNA-seq 方法在测试均值差异时表现不佳。
图1. Memento 工作流程用于差异均值、变异性和基因相关性测试 (A) 单细胞 RNA 测序 (scRNA-seq) 样本在文库制备和测序过程中每个细胞内的 RNA 转录本。在 scRNA-seq 采样后，观察到的转录本计数的基因表达的均值、变异性和相关性的模式不再类似于实际分布。 (B) Memento 将 scRNA-seq 建模为超几何抽样过程，使用矩估计法估计表达分布参数（均值、残差方差和相关性），通过有效的自助法估计置信区间 (CIs)，并在两组细胞之间测试表达参数的差异。 (C) Memento 在四个应用中表征了约 70,000 个人类气管上皮细胞对外源性细胞因子的反应，从约 170,000 个人类 CD4+ T 细胞中重建基因调控网络，这些细胞受到 CRISPR-Cas9 的扰动，绘制了来自 162 名系统性红斑狼疮 (SLE) 患者和 99 名健康对照者的 1.2 百万外周血单核细胞 (PBMCs) 的基因表达遗传决定因素，并在 CELLxGENE Discover 数据库中比较了任意组细胞。另见图 S1。
Para_04
为了应对这些统计和方法上的挑战，我们提出了 Memento，这是一种端到端的方法，用于从单细胞 RNA 测序（scRNA-seq）数据中估计均值、残差方差和基因相关性，并为这些参数的假设检验提供了统计框架（图 1B）。 Memento 采用多变量超几何抽样过程，并利用 scRNA-seq 数据的稀疏性，实现了一种高效的自举策略，用于在不同细胞群之间对估计参数进行统计比较。 通过模拟和实际数据分析，我们证明了 Memento 在广泛的基因表达分布和采样效率下能够产生准确的参数估计，计算出适合多重检验校正的校准良好的检验统计量，并实现亚线性运行时间。 我们在四个应用中展示了 Memento 的广泛适用性，旨在阐明实验和遗传因素如何影响人类细胞中的基因表达分布（图 1C）。 首先，我们对 70,000 个气管上皮细胞进行了 scRNA-seq 分析，这些细胞受到细胞外干扰素（IFNs）的刺激，并研究了刺激如何随时间调制响应基因的变异性和相关性。 其次，我们对 170,000 个 T 细胞进行了 Perturb-seq 分析，绘制了定义广泛 T 细胞激活方面的基因调控网络。 第三，我们重新分析了从 250 名个体收集的 120 万个细胞，以识别特定细胞类型中与均值、变异性和基因相关性相关的遗传变异。 最后，我们利用 Chan Zuckerberg Initiative (CZI) CELLxGENE Discover Census API 实现了一种近似自举策略，便于在 5000 万个细胞的 CELLxGENE 数据库中对任意细胞群进行近乎实时的比较。 在这多样化的应用中，Memento 与现有方法相比，一致地识别出实验组之间更多的显著且可重复的均值表达差异。它还识别出了表达变异性和基因相关性的差异，从而揭示了由扰动引起的不同的转录调控模式。 Memento 使用 Python 实现，兼容 scanpy，并可以从 https://github.com/yelabucsf/scrna-parameter-estimation 下载。
Results
Statistical model of scRNA-seq
单细胞 RNA 测序的统计模型
Para_01
自从诞生以来，单细胞 RNA 测序（scRNA-seq）尽管在分子生物学方面不断进步，但仍然产生了稀疏数据，即使在相同环境中暴露的遗传上相同的细胞之间也表现出高度的细胞间变异（图1A）。 将这种变异分解为生物变异和测量噪声的成分对于单细胞 RNA 测序数据的差异表达分析至关重要。
待补充
图 S1. 验证超几何抽样模型，与图 1 相关 (A) 单步超几何抽样很好地近似了从捕获和测序过程中的复合抽样。 (B) 表征将两步抽样近似为单步超几何抽样的影响。顶部：不同捕获效率和测序饱和度下分布之间的 Wasserstein 距离的热图。 底部：在单一捕获效率和测序饱和度下两种不同抽样过程的直方图（A 中的一个单元格）。 (C) 每个点代表模拟中的一个基因。y 轴是方差估计误差（|ln(v_true) - ln(v_estimated)|），x 轴是基因的真实模拟均值。
Estimating distributional parameters of gene expression from scRNA-seq
从单细胞 RNA 测序数据估计基因表达的分布参数
Para_01
据我们所知，这是首次使用超几何抽样过程对单细胞 RNA 测序数据进行建模，这可能是由于通过最大似然估计分布参数的复杂性所致。 在这里，我们在假设超几何抽样的情况下，推导了给定 Y_c 的 X_c 的一阶（均值）、二阶（方差）和混合（协方差）矩的方法矩估计器（详见 STAR 方法部分推导和详细信息）： μ̂_g,memento = 1 / (∑_c N_c) * r_g^，其中 r_g^ 是基因 g 的 Good-Turing 校正计数； σ̂_g,memento^2 = 1 / n_cells * c * (Y_cg^2 - Y_cg * (1 - q)) / (N_c^2 * q^2) - μ̂_g^2； σ̂_g_i,g_j,memento = 1 / n_cells * ∑_c (Y_cg_i * Y_cg_j) / (N_c^2 * q^2) - μ̂_g_i * μ̂_g_j
待补充
Para_03
我们进行了广泛的模拟，将 Memento 的超几何估计器与 scHOT 使用的朴素插件估计器、Zhang 等人引入的泊松近似下的经验贝叶斯估计器（Memento 估计器的一个特例，设置 q = 0）以及从 BASiCS 得出的估计值进行了比较（见 STAR 方法中的朴素和泊松估计器的形式）。 在一系列 q 值下，Memento 的超几何估计器产生了准确的均值估计（林氏一致性相关系数 ρc > 0.8，10 个细胞；ρc > 0.98，100 个细胞），残差方差（ρc > 0.98，100 个细胞）和基因相关性（ρc > 0.98，100 个细胞）（图 2A）。 此外，Memento 在不同 q 值下产生的残差方差和基因相关性估计值稳定，优于其他估计器，适用于低效率和高效率的单细胞 RNA 测序工作流程。 尽管所有估计器对高表达基因的准确性更高，但 Memento 即使对于低表达基因也优于其他方法（图 S1C）。 这些模拟基于单步采样方法，如上所述，该方法有效地近似了建模 RT 和测序的两步采样过程。
图2. Memento 在模拟和真实数据中的性能 (A) 使用10个细胞估计均值（左）、使用100个细胞估计变异（中）以及使用100个细胞估计基因相关性（右）与模拟的真实值（y轴）的林氏一致性，不同总体转录捕获效率（x轴）。阴影区域表示标准误差。 (B) Memento 从 DropSeq 数据估计的值与 smFISH 估计的同一黑色素瘤细胞群体的值（y轴）的皮尔逊相关性，对于均值（左）、变异（中）和基因相关性（右），使用不同数量的 DropSeq 细胞（x轴）。阴影区域表示标准误差。 (C) 功效（y轴）与假发现率（FDR）（x轴）比较现有方法与 Memento 进行 DM（左）、DV（中）和 DC（右）分析。 (D) 单细胞 DM 分析（绿色）与伪批量 DM 分析（红色）的一致性曲线下面积（AUC）（x轴），使用 Squair 等人和 Perez 等人的数据集。 (E) 三种方法在 DM 和 DV 分析中随细胞数量（x轴）的运行时间（y轴）。另见图 S2 和 S3。
Para_03
为了进一步验证 Memento 参数估计的准确性，我们重新分析了一组包含基于液滴的单细胞 RNA 测序（DropSeq）和单分子荧光原位杂交（smFISH）数据的数据集。 该数据集之前使用 SAVER 进行了分析，SAVER 是一种通过借用相似基因和细胞的信息来估计基因相关性的插补方法，已被证明优于其他方法（图 2B）。 对于同时使用 DropSeq 和 smFISH 测定的基因，当使用非常少的细胞时，Memento 的平均估计值比其他方法使用的简单估计器略有改进（考虑了 21 个基因；使用 100 个细胞时，ρ = 0.58 和 ρ = 0.54）。 对于残差方差，Memento 的估计值与 smFISH 获得的估计值显著更相关（考虑了 14 个基因；ρ = 0.71），而简单估计器（ρ = 0.56）和 BASiCS（ρ = 0.61）使用所有可用的 8,498 个细胞时则不然。 最后，对于基因相关性，Memento（ρ = 0.53）也显著优于简单估计器（ρ = 0.29）、SAVER（ρ = 0.38）和 scVI（ρ = 0.23），使用所有细胞时也是如此。 重要的是，Memento 在不使用插补方法（如 SAVER）和变分推断方法（如 scVI）所需的额外基因的情况下，产生了更好的基因相关性估计值。 这一优势不仅体现在估计的计算效率上（Memento，17 秒 vs. SAVER，30 分钟，用于 14 个基因对），还可能产生更适合特定下游分析（如遗传映射）的估计值，因为在这些分析中，插补可能会无意中引入混淆效应。 这些结果强调了 Memento 参数估计的准确性，这通过模拟和与基准 smFISH 数据的比较分析得到了证明。
Hypothesis testing using highly efficient bootstrapping
使用高效自助法进行假设检验
Para_01
假设检验的目标是确定细胞组之间估计参数（如均值、变异性和基因相关性）的观察差异是否在统计上显著，与零假设相比。 当测试成千上万的基因时，这是单细胞 RNA 测序实验中常见的问题，多重检验问题是主要关注点：在预测预期验证数量的同时，提名一组可行的候选基因进行实验跟进。 因此，在多重检验校正下，适当校准零假设下的检验统计量变得至关重要。 尽管使用矩估计法提供了计算效率和建模灵活性，但建立估计参数的统计显著性需要通过数据的自助法计算置信区间（CIs）。 使用标准方案对大量细胞进行自助法抽样（有放回抽样）将需要大量的计算资源，这在时间和内存上都是不可行的，尤其是在处理大型数据集时。
Para_02
Memento 实现了一种利用 scRNA-seq 数据稀疏性的洗牌方案，以实现快速、内存高效且高度并行化的自助法。 我们的方案基于一个关键观察，即观察到的独特转录本计数的数量远小于细胞数量（图 S2A），即使对于独特观察到的计数对也是如此（图 S3B），尽管程度较小。 因此，每个自助法迭代只需要重新采样每个基因的 K 个独特转录本计数，这些计数从多项分布中按观察频率比例重新采样（图 S2C），而不是从包含 N 个元素（细胞）的多项分布中重新采样单个细胞的计数（多项分布（N, 1/N, …, 1/N））。 这种方法最终使得每次重新采样迭代拟合一个显著较小的加权数据集（K << N）。 为了适应多重实验，我们使用元回归框架扩展了自助法策略，将每个重复视为数据的单独子组，从而实现分层重新采样。 这种方法使我们能够在尊重数据生成过程的同时量化不确定性，例如从不同个体中采样的细胞。 在模拟中，Memento 的自助法策略在广泛的基因范围内产生了与朴素自助法重新采样相当的均值、残差方差和基因相关性的零分布的高精度估计（图 S2D 和 S2E）。 通过使用自助法量化参数估计的置信区间，Memento 计算了适合多重检验校正的 DM、DV 和差异相关性（DC）的校准经验 p 值（图 S2F）。
图 S2. 单细胞 RNA 测序数据的离散性质使得准确、快速重采样成为可能，与图 2 相关 (A) 在 IFN-B 数据集中随机选择的基因对的唯一基因对数量。 (B) 在 IFN-B 数据集中随机选择的基因对的唯一基因对数量。 (C) 重新采样频率而非表达量的概念图。 (D) Memento 与全样本重采样的高效对比。代表性基因和基因对的对数均值、对数残差方差和相关性的样本分布。 (E) Memento 重采样与朴素重采样的全面比较。y 轴是通过比较 Memento 和朴素重采样的分布得到的 KS 检验 p 值，分别针对均值（顶部）、残差方差（中部）和相关性（底部）。x 轴是基因的平均表达量分位数。 (F) 从 Memento 计算出的 p 值的代表性示例。
图 S3. 进一步验证 Memento 的差异均值方法，与图 2 相关 (A) 低表达和高表达基因在模拟中的错误发现率。(B) 在各种方法中 IFN-A 和 IFN-B 之间差异表达基因的一致性。
Para_02
为了证明 Memento 在保持高统计功效的同时能够准确估计假阳性，我们模拟了一个包含两个不同细胞群体的数据集。 为了与实际数据保持相关性，我们使用了从重组干扰素-β（rIFNB）刺激前后的 CD4+ T 细胞的真实数据集中提取的参数。 我们展示了对于 DM、DV 和 DC，Memento 在指定显著性阈值下估计了预期的假阳性数量，同时在检测真实参数差异方面实现了最高的功效（图 2C）。 此外，我们观察到现有的单细胞 RNA 测序 DM 方法过于宽松（t 检验、Wilcoxon 秩和检验），而伪批量 DM 方法则过于保守（edgeR、DESeq2），这与 Squair 等人的结果一致（图 S3A）。 Squair 等人此前将这一结果归因于大多数单细胞 RNA 测序数据集中存在的复制水平异质性，并建议使用伪批量方法来简化层次结构。 通过直接考虑层次结构，Memento 即使在存在不同程度的异质效应时，也能在每个显著性阈值下产生预期的假阳性率（FPR）。 除了模拟之外，我们还使用配对的单细胞和批量 RNA 测序样本对 Memento 进行了基准测试，使用了 Squair 等人使用的数据集（图 2D，左）和一个来自系统性红斑狼疮（SLE）患者的额外数据集（图 2D，右）。 在这两个数据集中，Memento 从单细胞 RNA 测序数据中产生的 DM 结果与批量 RNA 测序分析结果最一致。 最后，Memento 在 IFN-α 和 IFN-β 刺激的纤毛细胞中识别出最多的一致差异表达基因，这两种都是已知的 I 型干扰素受体的配体（图 S3B）。
Para_03
与现有的 DM、DV 和 DC 方法相比，Memento 在计算速度上快了几个数量级，允许扩展到数百万个细胞（图 2E）。 在一个规模与新兴 scRNA-seq 数据集相当的模拟中（每组包含 10^6 个细胞），对 1,000 个基因进行 DM 和 DV 分析，每个基因使用 10,000 次引导迭代仅需 13 分钟使用单个 CPU。 Memento 的多核实现促进了多个基因的并行化，进一步将运行时间减少到使用 6 个 CPU 时的 2-3 分钟。 特别是在 DV 和 DC 分析中，Memento 使用等效计算资源比现有方法实现了高达 1,000 倍的计算速度提升。 这些结果证实了 Memento 的引导策略在高计算效率下产生了准确的效果大小置信区间估计。 这最终导致了校准良好的检验统计量，使得 scRNA-seq 数据的假设检验可以扩展到包含数百万个细胞的组（有关重采样策略和假设检验的详细描述，请参见 STAR 方法）。
图3. 使用 Memento 映射 HTECs 对细胞外干扰素的转录反应 (A) 整个 HTEC 数据集的 UMAP 图，按识别的细胞类型着色（左），放大后的纤毛细胞按刺激着色（中），以及时间标签（右）。 (B) 在 6 小时后对 IFN-α（x 轴）与 IFN-β（左）、IFN-γ（中）和 IFN-λ（右）的平均表达量的对数倍变化（LFC）。 (C) 四种类型干扰素（每个热图中的列）在 5 个时间点上的 LFC 的分层聚类热图。突出显示了 I 型（绿色）和 II 型特异性（蓝色）反应。 (D) 随时间变化的基因共表达网络，其中青色节点表示典型 ISGs，洋红色节点表示非典型 ISGs。高相关性（Memento ρ > 0.6）的基因对相连。 (E) 纤毛细胞的基线表达变异性（y 轴）与均值（x 轴）的关系。 (F) 典型和非典型 ISGs 的持续敏感性（y 轴）与均值（x 轴）的关系。***p < 0.001。 (G) 对 IFN-β（左）和 IFN-γ（右）反应的变异性变化（y 轴）与均值变化（x 轴）的关系。蓝色点代表典型 ISGs。另见图 S4。
Differential variability and gene correlation in response to exogenous IFN
对外源性IFN的反应中的差异变异性和基因相关性
Para_01
虽然干扰素（IFNs）是促进抗病毒免疫的强大细胞因子，它们也在炎症和自身免疫性疾病的发生中起作用。 它们的作用——通过自分泌和旁分泌信号诱导基因表达——已经得到充分记录；然而，刺激细胞中的转录组反应的异质性仍 largely 未被探索。 使用 Memento，我们研究了干扰素刺激对人支气管上皮细胞（HTECs）基因表达分布的影响。 我们使用多重单细胞 RNA 测序（mux-seq）分析了来自两名健康供体的 69,958 个 HTECs，探索了包括未刺激对照和各种干扰素刺激的条件：I 型（IFN-α 和 IFN-β）、II 型（IFN-γ）和 III 型（IFN-λ）。 分析在多个刺激后时间点进行：3、6、9、24 和 48 小时。 降维、最近邻识别和 Leiden 聚类产生了 7 种可识别的细胞类型，使用均匀流形近似和投影（UMAP）可视化：神经内分泌细胞、离子细胞、毛细胞、基底细胞、基底/俱乐部细胞、杯状细胞和纤毛细胞（图 3A）。 我们随后的分析仅集中在纤毛细胞上，已知纤毛细胞是病毒性感染的主要靶标，包括 SARS-CoV-2，并且以其强大的干扰素反应而闻名。
Para_02
我们确定了5018个基因表现出差异平均表达（DMGs，错误发现率 [FDR] < 0.01），在未刺激的纤毛细胞和用四种干扰素中的任何一种刺激6小时后的纤毛细胞之间。 比较分析显示，干扰素-α在DMGs中的倍数变化（FCs）与干扰素-β和干扰素-λ相似（ρ = 0.96）。 相比之下，与干扰素-γ相比，整体倍数变化的相关性较低（ρ = 0.70；图3B），这是由于存在I型和II型干扰素特异性的DMGs。 在此，我们将响应任何干扰素上调的DMGs定义为干扰素刺激基因（ISGs）。 ISGs在时间点上的层次聚类揭示了跨干扰素的动态转录组反应，包括主要组织相容性复合体（MHC）II类基因的早期诱导和一个独特的基因簇，包括PLAAT2、BTN3A1和DUOX2（图3C）。 我们还识别了每种干扰素特有的模式，例如，一组典型的ISGs（IFIT2、IFITM2和ISG15）在响应干扰素-λ时表现出晚期诱导，而在响应I型干扰素时在整个时间过程中持续诱导（图3C）。 有趣的是，某些由一种干扰素更多诱导的基因（例如，由干扰素-γ诱导的MHC II类基因）在其他干扰素中也表现出相似的时间行为，这表明存在独特和共享的调控机制。
Para_03
虽然 DM 分析揭示了经典和非经典 ISGs 的诱导，但它未能阐明这些基因是否受到相同的转录调控控制。 为了绘制 IFN 基因相关网络及其子组件，我们使用 Memento 识别 ISG 对在不同刺激和时间点之间的差异相关性（图 3D）。 对所得基因相关矩阵进行聚类分析，揭示了对 IFN-β 反应的不同 ISG 亚群，形成了未刺激细胞、刺激细胞或两者的簇——这些区别仅通过 DM 分析是无法识别的。 例如，包括 MX1、OAS1 和 IFI6 在内的经典 ISGs 即使在没有外源 IFN 存在的情况下也保持高相关性（图 3D，青色节点）。 在 IFN-β 刺激下，最初由经典 ISGs 组成的相关网络扩展到包括非经典 ISGs，如 MHC I 类分子和其他与抗原呈递相关的基因，这些基因在未刺激细胞中不相关（图 3D，洋红色节点）。 与聚类分析一致，Memento 识别出更多非经典 ISGs 中的差异相关基因对（DCGs，FDR < 0.1）（860 个 DCGs，占总对数的 34%）而不是经典 ISGs（421 个 DCGs，占总对数的 16%）。 值得注意的是，当考虑所有基因对以及仅考虑平均表达量显著变化的对时，基因对之间相关性的增加不能用其平均表达量的增加来解释（图 S4A）。
图 S4. HTEC 数据中基因共表达和变异模式，与图 3 相关 (A) 两个基因之间的相关性变化（y 轴）与平均值变化的乘积（x 轴）的关系（左图），以及筛选出平均值变化较大的基因对（右图）。 （B）每种类型基因在不同细胞类型中的表达变异性（y 轴）。
Para_03
我们假设，在未刺激的细胞中，经典 ISGs 由于感知到持续的 IFN 而相互关联，并且在特定细胞群内协调诱导 ISGs。 持续的 IFN 已被描述为在细胞间自然地形成 ISG 表达梯度，并在病毒防御、免疫细胞稳态和自身免疫中发挥重要作用。 在我们的数据集中，与非经典 ISGs 相比，经典 ISGs 在未刺激的细胞中表现出更大的变异性（图 3E），这与先前记录的细胞因子和非细胞因子之间表达变异性的差异一致（图 S4B）。 在使用 Memento 识别出的 761 个差异可变基因（DVGs，FDR < 0.1）中，未刺激的纤毛细胞与用四种 IFNs 刺激 6 小时后的细胞相比，有 394 个在未刺激的细胞中高度可变（FDR < 0.005），并且富集了经典 ISGs（GSEA IFN-α/IFN-β 信号调整 p = 3.35 × 10−12），包括 IFIT1、IFIT3 和 MX1。
Para_04
我们接下来评估了每个 ISG 对稳态 IFN 的敏感性，这通过比较 Ifnar 基因敲除小鼠和野生型（WT）小鼠在没有外源 IFN 条件下巨噬细胞基因表达的倍数变化来估计。 这一分析表明，经典 ISG 对稳态 IFN 的敏感性显著高于非经典 ISG（p < 2.73 × 10^-10；图 3F）。 值得注意的是，在用 IFN-β 刺激后（以及在较小程度上，用 IFN-γ 刺激后），相当比例的经典 ISG 的变异性降低（分别为 78% 和 39%；图 3G，FDR < 0.1），这意味着外源刺激使细胞环境均一化，消除了对稳态 IFN 反应的异质性影响。
Para_05
我们的研究结果突显了 Memento 在分析基因表达分布和揭示受 IFN 信号影响的转录调控网络方面的强大能力。 通过利用 Memento 来解析基因表达的均值、方差和相关性的影响，我们阐明了在存在和不存在 IFN 的情况下决定细胞行为的复杂调控相互作用，为细胞如何调节其对环境信号的转录组响应提供了新的视角。
Differential expression analysis of perturbed CD4+ T cells maps gene-regulatory networks in T cell activation
扰动的CD4+ T细胞差异表达分析绘制了T细胞活化中的基因调控网络
Para_01
将 CRISPR-Cas9 介导的基因组扰动与单细胞 RNA 测序（scRNA-seq）结合，为在多种体外系统中进行正向遗传筛选创造了新的机会。 利用 Memento，我们分析了约 173,000 个经过 CRISPR-Cas9 扰动的人类 CD4+ T 细胞，以绘制调节其激活和极化的转录调控网络。 细胞通过单导向 RNA (sgRNA) 慢病毒感染与 Cas9 蛋白电穿孔（SLICE）进行扰动，随后进行多路测序（mux-seq）。 使用一组 280 个 sgRNA，我们针对了 140 个转录调控因子（TRs），这些因子因其高表达（在批量 RNA 测序中位于前四分之一）或其结合位点的差异可及性（通过批量转座酶可及染色质测序 [ATAC-seq] 检测到）而被选择在活化的 CD4+ T 细胞中。 在 Cas9 电穿孔和多轮选择和增殖后，来自 9 名供体的活化 CD4+ T 细胞使用多路测序进行了分析。
图4. 使用 Perturb-seq 和 Memento 重建 T 细胞激活的基因调控网络 (A) 本研究中扰动调节因子的选择标准，基于表达（上）和结合位点富集（下）。 (B) 每个基因（行）在被相应 sgRNA 扰动的细胞（列）中的平均基因表达热图。 (C) 左：每个 sgRNA 扰动的细胞（列）中 DMGs（行）的平均基因表达热图。右：从 WT 细胞估计的相同 DMGs 的基因-基因相关矩阵。 (D) 每个调节因子与其下游靶标在 WT 细胞中的相关性。 (E) 未考虑调节因子之间相互作用的双部分基因调控网络，由 Perturb-seq 数据的 DM 分析构建。 (F) 包括调节因子之间遗传相互作用的基因调控网络，利用 DM 和 DC 分析构建。 (G) 在 TSS 周围不同窗口内，相互作用或非相互作用调节因子对的结合位点数。误差条表示标准误差。 (H) LGALS3BP 的染色体位置以及预测通过 DM 和 DC 分析相互作用的 IRF1 和 PRDM1 的结合位点。另见图 S5。
Para_01
为了评估每个 sgRNA 的切割效率，我们对来自每个供体的 sgRNA 池和编辑细胞的 DNA 中的 280 个位点中的 268 个进行了靶向扩增测序。 268 个 sgRNA 的平均切割效率定义为编辑细胞在目标位点的测序覆盖率与相应 sgRNA 在池中的测序覆盖率之比，估计为 21%，标准差为 15%（图 S5A）。 14 个 sgRNA 的切割效率低于 2.0%（标准差 1.7%；Z 分数，p < 0.05），被指定为未切割阴性对照（WT）。 我们的筛选的稳健性和有效性通过两项质量控制分析得到了证实。 首先，我们使用 Memento 确认了转导相应 sgRNA 的细胞中目标基因显著下调（图 4B）。 其次，观察到 WT 细胞之间的平均基因表达相关性更高（ρ = 0.50）或转导相同基因靶标的 sgRNA 的细胞之间的相关性更高（ρ = 0.44），而转导不同基因靶标的 sgRNA 的细胞之间的相关性则为 0（KS 测试 p < 2.2 × 10^−16；图 S5B）。
图 S5. 与图 4 相关的基因扰动和自然遗传变异分析的附加实验细节和结果 (A) 我们估计了每个 sgRNA 的平均 sgRNA 敲除效率（x 轴）（y 轴）。每个点代表平均敲除效率，误差线表示不同供体之间的标准偏差。(B) WT 引导、针对同一基因的引导和随机引导对的转录组相关性分布。(C) 细胞类型特异性 ATAC 峰中的 eQTL 富集（欧洲人）。(D) 使用通过 Memento 和伪批量方法找到的 eGenes 进行疾病 LDSC 分数回归富集。
Para_01
利用 Memento，我们在比较野生型细胞与至少一个 sgRNA 扰动的细胞时，鉴定出 7,641 个具有差异甲基化基因（DMGs）的基因（FDR < 0.05）。 对 sgRNA 之间的 DMGs 平均基因表达进行层次聚类，揭示了施加类似转录组效应的 sgRNA 聚类以及对这些扰动反应类似的基因聚类（图 4C）。 我们鉴定出五个与核糖体（FDR < 5.35 × 10^-24）、细胞毒性（FDR < 0.014）、抗原呈递（FDR < 0.0011）和增殖（FDR < 0.001）显著相关的 DMGs 聚类。 此外，使用 Memento 计算的 DMGs 的成对相关矩阵显示，在初始五个 DMG 聚类中的每个聚类内存在额外的亚聚类，并且在野生型和扰动细胞中均持续存在（图 4C）。 有趣的是，虽然抗原处理基因的平均表达受一组共享的 sgRNA 调控，但 MHC II 类基因——特别是 HLA-DPA1、HLA-DRA、HLA-DRB1 和 HLA-DPB1——表现出强烈的关联，表明它们的协调表达可能受额外的转调控因子控制。
Para_02
在探索 Memento 在检测基因相关性变化中的实用性时，我们假设可以在不进行组合扰动的情况下识别 TRs 之间的遗传相互作用。 为了验证这一假设，我们进行了遗传相互作用分析，重点关注由 DMGs 及其 TRs 组成的对，称为 TR-DMGs（见 STAR 方法）。 具体来说，我们关注那些在敲除后导致 DMGs 表达减少的调节因子。 与我们的预期一致，TR-DMGs 通常在 WT 细胞内表现出正相关（二项检验，p < 0.00668；图 4D）。
Para_03
在没有基因相互作用的情况下，两个转录调节因子（R1 和 R2）可以独立调控目标基因（G）；因此，敲除一个调节因子不应明显影响另一个调节因子的功能（图4E）。 相比之下，在存在相互作用的情况下，敲除一个调节因子（例如，R1）可能会影响 R2 对 G 的调控能力。 这种效应可以通过检测当 R1 被扰动时 R2 和 G 之间基因相关性的变化来识别（图4F）。 采用这一策略，我们确定了 564 个基因相互作用，涉及 432 个独特的调节因子对（FDR < 0.1，图4F）。 验证这些相互作用，结合 ENCODE 的染色质免疫沉淀测序（ChIP-seq）数据分析显示，相互作用的转录调节因子对更可能在目标基因转录起始位点（TSS）附近具有共定位的结合位点，而非相互作用的对则不然（图4G）。
Para_04
例如，我们确定 IRF1 调控 LGALS3PB（从 DM 表达分析中明显看出），并且在 WT 细胞中与 LGALS3PB 保持强烈的相关性（ρ = 0.28）。 PRDM1 的敲除导致 IRF1 和 LGALS3PB 之间的相关性显著下降（Δρ = −0.38），这表明 PRDM1 和 IRF1 在调控 LGALS3PB 方面可能存在相互作用。 与这些观察结果一致，LGALS3BP 在其 TSS 附近有 IRF1 和 PRDMB1 的结合位点（图 4H）。
Para_05
这些结果证明了 Memento 在分析正向遗传 Perturb-seq 筛选中的能力。 我们强调了 DC 分析在划分共享调控元件的基因集方面的潜力——尽管这些基因集参与了不同的途径——以及重建调控 T 细胞激活的转调控因子的遗传相互作用。
Genetic analysis of population-scale scRNA-seq
群体规模单细胞 RNA 测序的遗传分析
Para_01
大规模人群单细胞 RNA 测序（scRNA-seq）数据的日益可用，为绘制特定细胞类型中与邻近基因（顺式）表达分布变化相关的遗传变异铺平了道路。 现有的研究主要利用伪批量方法，如矩阵表达数量性状位点（eQTLs），来识别影响平均表达的顺式 eQTLs（cis-eQTLs）。 尽管线性混合模型最近已被应用于 scRNA-seq 数据中的顺式 eQTLs 映射，但它们受到计算效率低、仅关注均值比较以及基础参数模型易受误指定的影响。 我们认为，与伪批量方法相比，Memento 的参数估计精度更高，并且在推断过程中能够考虑个体内部和个体之间的变异，这将提高检测顺式 eQTLs 的能力，并发现新的变异性和相关性 QTLs（vQTLs 和 cQTLs，分别）。 此外，高效分层自助策略的实施有望适用于大规模人群 scRNA-seq 数据集，而这些数据集对于参数化线性混合模型来说可能是计算上无法克服的。 为了证明这一点，我们将 Memento 应用于重新分析一个现有的 scRNA-seq 数据集，该数据集包括来自 162 名系统性红斑狼疮（SLE）患者和 99 名健康捐赠者的 120 万个外周血单核细胞（PBMCs）。
Para_02
数据分别针对每种报告的细胞类型进行了分析：CD4+ T 细胞（T4）、CD8+ T 细胞（T8）、自然杀伤（NK）细胞、经典单核细胞（cMs）和非经典单核细胞（ncMs）。 东亚血统和欧洲血统的人群分别进行了分析，并通过后续比较在这两个人群之间进行了复制分析。 对于每种不同的细胞类型和血统群体，Memento 映射了顺式遗传变异——特别是那些距离转录起始位点（TSS）100 kb 以内的变异——与平均表达量、表达变异性以及基因相关性相关联，生成了校准良好的 p 值（图 5A）。
图5. 使用Memento映射平均QTL（eQTL）、vQTL和cQTL。(A) Memento计算的预期p值（y轴）与理论p值（x轴）的分位数-分位数（QQ）图，用于eQTL、vQTL和cQTL。对于eQTL，叠加了伪批量方法（矩阵eQTL）的p值QQ图。 (B) Memento和基于伪批量的矩阵eQTL在从一个更大队列（OneK1K）中识别出的eQTL恢复中的接收者操作特征（ROC）曲线。 (C) Memento和矩阵eQTL在不同个体数量下的eQTL恢复能力（y轴）。分析在CD4+ T细胞（T4）、B细胞（B）、经典单核细胞（cMs）和自然杀伤细胞（NKs）上进行。 (D) 细胞类型特异性eQTL在细胞类型特异性ATAC峰中的富集。每个条目代表在一个细胞类型（列）中检测到的eQTL在另一个细胞类型（行）中检测到的ATAC峰中的富集程度。强度为−log10(p值)。 (E) 在每个细胞类型中检测到的eQTL在相同细胞类型中检测到的细胞类型特异性ATAC峰中的富集。误差线表示标准误差。 (F) 一个vQTL的例子。每个个体在chr6:31326612处基因型不同的表达变异（y轴）。 (G) 显示每个基因型代表性个体HLA-C表达分布的直方图。 (H) 一个cQTL的例子。不同基因型个体在chr12:69688073处的JUNB-LYZ基因相关性（y轴）。 (I) 所有供体（灰色）和一个代表性个体（黑色）的单细胞中LYZ表达（y轴）与JUNB表达（x轴）的散点图。
Para_02
Memento 和 Matrix eQTL 在检测顺式 eQTL 方面的效能和假阳性率进行了比较分析，基准由 OneK1K 研究提供的 1,000 名非重叠个体组成。 值得注意的是，在东亚和欧洲队列中，Memento 在识别顺式 eQTL 方面表现出更高的效能（AUC = 0.85），超过了 Matrix eQTL（AUC = 0.81），同时保持了相同的假阳性率（图 5A 和 5B）。 总体而言，Memento 在两个群体中均优于 Matrix eQTL，分别在东亚人中跨细胞类型复制了 1,606 个 vs. 855 个顺式 eQTL，同样在欧洲人中复制了 1,778 个 vs. 958 个顺式 eQTL。 此外，在常见的多路测序实验队列规模范围内，Memento 在 80 个个体的情况下平均效能提高了 15%——这一指标在 50 个个体的情况下增加到 32%，每个个体平均有 440 个细胞（图 5C）。
Para_03
我们随后探讨了 Memento 检测到的顺式 eQTL 数量增加是否也提高了开放染色质区域内的富集以及与疾病的关联。 在东亚队列中，Memento 在特定细胞类型中鉴定出的顺式 eQTL 在由另一项对批量分选免疫细胞进行 ATAC-seq 的研究注释的细胞类型特异性开放染色质区域中更为丰富（匹配细胞类型的 p 值：B 细胞，9.0 × 10−9 对比 0.04；T4 细胞，9.3 × 10−4 对比 0.11；T8 细胞，0.03 对比 0.58；NK 细胞，6.67 × 10−8 对比 0.03；cM 细胞，2.1 × 10−11 对比 0.67；ncM 细胞，1.0 × 10−6 对比 0.46；图 5D 和 5E）。 在欧洲队列中也观察到了类似的富集增益（图 S5C）。 连锁不平衡（LD）评分回归（LDSC）分析发现，Memento 鉴定出的顺式 eQTL 也更丰富地与全基因组关联研究（GWAS）中的免疫介导疾病相关联，这表明细映射性能有所提高（图 S5D）。
Para_04
除了映射顺式 eQTL 外，Memento 还能够识别与表达变异和基因相关性相关的遗传变异，提供了遗传变异可能影响基因表达的替代机制的见解。 利用 Memento，我们确定了 10,607 个影响所有细胞类型中 733 个基因的表达 vQTL。 例如，HLA-C 表达的变异性在 chr6:31326612 的不同基因型之间有所不同（图 5F），A 等位基因放大了 HLA-C 的表达变异性，但对其平均值没有显著影响（图 5G）。 对于映射 cQTL，我们专注于测试至少具有一个显著顺式 eQTL 的基因与已知转录因子之间的相关性，从而具体测试遗传变异可能调节转录因子对基因表达的影响的假设。 我们在所有细胞类型中映射了 3,726 个 cQTL，涉及 238 个基因对。 例如，chr12:69688073 的 SNP 不仅与 LYZ 的平均表达有关，还与 JUNB 和 LYZ 之间的相关性有关。 有趣的是，SNP 周围 1 kbp 内存在一个 JUNB 结合位点，这表明 JUNB 可能作为 LYZ 的反式调控因子，其调控强度受该 SNP 位点的基因型影响。
Para_05
这些发现强调了 Memento 作为一种可扩展的方法，适用于大规模人群单细胞 RNA 测序数据的遗传分析，提供了更高的统计能力来识别顺式表达数量性状位点（cis-eQTLs），并引入了映射变异数量性状位点（vQTLs）和细胞类型特异性数量性状位点（cQTLs）的能力。 这些进展不仅提高了疾病关联的精细定位，还揭示了遗传变异可能调节基因表达的新机制。
Census-scale differential expression analysis across cell types, individuals, and disease states
跨细胞类型、个体和疾病状态的人口规模差异表达分析
Para_01
上述应用展示了 Memento 在不同数据集上进行广义差异表达分析的广泛适用性，包括对 IFNs 刺激的气管上皮细胞时间响应的分析、从 CD4+ T 细胞的 Perturb-seq 数据映射基因调控网络，以及对人群队列收集的单细胞 RNA 测序数据的大规模遗传分析。 这些应用和模拟表明，Memento 一贯优于现有方法，提供了一套独特的功能来比较方差和相关性，而不仅仅是均值，并且极其高效，能够扩展到数百万个细胞和数万个重复样本。
Para_02
全球范围内大规模的单细胞 RNA 测序（scRNA-seq）数据存储库的出现，对能够高效比较数据集并在确保适当校准的统计行为的计算技术提出了新的需求。 截至 2023 年 11 月，CELLxGENE Discover 包含了来自 1,102 个数据集和数千名个体的 5000 万个独特细胞，其 Census API 提供了对其中大多数数据的访问。 与由单一研究项目生成的、具有明确假设的 scRNA-seq 数据集不同，CELLxGENE Discover 的用户访问这一资源时，心中有着多种多样的比较分析目标。 例如，一个用户可能对不同器官系统中相同细胞类型的表达差异感兴趣。 另一个用户可能对不同疾病状态个体之间相同细胞类型的表达差异感兴趣。 在任何带有标记细胞类型的 scRNA-seq 数据集中，细胞组之间的可能比较数量都非常大（图 6A 和 6B）。 此外，多个数据集可以合并，以提高跨数据集存在的相同细胞组之间的比较能力。
图6. 扩展 Memento 用于 CZI CELLxGENE Discover 内的近实时差异表达分析 (A) 在 CELLxGENE 内的 SLE PBMC 数据集的 UMAP 图。 (B) 可以在细胞组内和组间进行的不同比较的数量。 (C 和 D) 预计算模式和完整模式之间的 (C) 差异均值和 (D) 差异变异性的显著性（p 值）比较。 (E) 查询时进行的比较次数与运行时间的关系（不包括预计算）。 (F) 使用 CELLxGENE 分析多个数据集以识别 pDCs 和 cDCs 之间的差异表达基因的示意图。 (G) 比较 pDCs 和 cDCs 时，结合多个数据集（青色）和单独使用每个数据集（灰色）的 p 值的 QQ 图。参见补充图 S6。
Para_02
差异表达方法需要在普查中高效地执行用户定义的细胞组之间的准确、校准良好的比较，以接近实时的速度交付结果，以便于网络门户集成。 尽管 Memento 在细胞数量增加时表现出色的可扩展性，如图2F所示，但其实时结果交付受到每次比较都需要进行引导操作的限制，当子集包含多个生物学和技术重复时，这一限制变得更加明显。 为了扩展 Memento 的广泛适用性，我们与 CZI 合作，利用 CELLxGENE Discover Census API 执行引导操作，并量化整个语料库中预定义细胞组的不确定性（见 STAR 方法）。 这种扩展允许预先计算标准误差，然后利用这些标准误差通过加权最小二乘法实现接近实时的差异表达分析。 因此，从这种预先计算模式得出的标准误差为全模式中使用的引导方法提供了有效的近似，简化了分析过程。
Para_03
为了评估 Memento 在预计算模式和完整模式之间的协议，我们对狼疮数据集中（图5所示）单个供体的 CD4+ T 细胞和 cMs 进行了差异表达分析，该数据集也包含在 CELLxGENE Discover 中。 鉴于分析涉及相同的基础数据，我们预期结果将非常相似。 主要差异将归因于两个 Memento 版本，预计算模式利用了整个 CELLxGENE Discover 数据集的估计细胞大小。 我们的预期通过观察完整模式和近似预计算模式之间效应大小估计值的强相关性（图 S6）得到了证实。 在显著性水平上也观察到类似的强相关性，由 −log10(p 值) 表示（图 6C 和 6D）。 重要的是，与在完整模式下执行 Memento 进行各种细胞组比较相比，确定效应大小和 p 值的计算时间显著减少（图 6E）。
图 S6. 使用完整版和预计算版 Memento 计算的效果大小比较，与图 6 相关
Para_03
Memento 在大规模普查数据上的独特应用在于其增强的能力来比较细胞群体，特别是对于在单个数据集中较为罕见的细胞群体尤为有益。 为了说明这一点，我们利用 Memento 的预计算模式来识别传统树突状细胞（cDC）和浆细胞样树突状细胞（pDC）之间的差异甲基化基因（DM 基因）。 这些细胞类型分别占 CELLxGENE Discover 中免疫细胞 scRNA-seq 数据集的 5.8% 和 4.0%（图 6F）。 在分析 CELLxGENE Discover 中的 23 个独立数据集，涵盖 362,619 个细胞时，我们发现跨这些数据集的联合分析显著提高了统计功效，与任何单个数据集的分析相比（图 6G）。 这些结果突显了 Memento 矩估计器的效率及其引导方法的适应性，使其能够有效地应用于广泛的普查存储库。
Discussion
Para_01
随着可扩展工作流程的发展，出现了单细胞 RNA 测序（scRNA-seq）数据集，其中在细胞组之间定量比较基因表达分布是一项关键任务。 这些任务包括在实验条件下比较单细胞表达谱、由基因组编辑引起的不同的遗传扰动以及继承不同等位基因的个体之间的比较。 最初的观察表明，实验和遗传扰动主要导致基因表达的微妙变化而非明确的细胞状态，这突显了需要能够比较基因表达分布的方法。 然而，能够促进对大量细胞和广泛协变量（例如数百种体外扰动或数百万遗传多态性）进行假设检验的可扩展计算方法仍然很少。 此外，目前能够测试基因表达变异性和基因相关性差异的方法更少，而这些参数是 scRNA-seq 独特捕获的。
待补充
Para_03
MoMs 估计器在分层模型中的实现普遍面临通过重采样建立置信区间（CIs）的挑战，因为将抽样过程纳入解析置信区间和 p 值的推导过程中，在没有进一步假设的情况下可能会变得极其复杂。 尽管重采样在计算上可能非常昂贵，尤其是在细胞数量较大时，我们采用的近似引导法重新采样的是唯一计数的数量，而不是单个细胞的数量。 由于我们的假设检验框架利用了近似引导法，理论上它应该与现有的参数模型和其他类型的估计器兼容，从而能够为各种单细胞测序分析方法提供更好的经验 p 值估计。 例如，可以设计一个估计器用于实验中，其中 mRNA 抽样过程不能被近似为单一步骤，而需要更深入的处理。 此外，我们还展示了 Memento 的原理可以扩展到以极高效的方式结合多个数据集进行差异表达分析，通过预先加载昂贵的计算，为研究人员提供更好的工具来与庞大的资源（如 CELLxGENE Discover）互动。
Para_04
通过在四个概念验证设置中的应用，我们证明了 Memento 在检测不同研究中差异表达基因方面具有更高的能力。 我们展示了我们的均值估计器在较低细胞数量时特别准确，且我们的推断与批量 RNA-seq 实验的结果更加一致。 此外，我们证明了差异变异性和相关性分析可以识别出使用 DM 分析无法检测到的新基因调控关系。 在多种数据集中得到验证，Memento 成为一种高度适应性和可扩展的方法，适用于包含许多重复和实验条件的大型单细胞 RNA-seq 数据集的定量分析。
Limitations of the study
研究的局限性
Para_05
当前版本的 Memento 存在一些限制，为后续版本提供了改进的机会。 首先，由于采用了引导方法，Memento 目前不支持包含细胞水平协变量或连续样本协变量。 在这两种情况下，一种方法可以是根据细胞水平协变量对细胞进行分组，或对连续协变量进行离散化。 其次，我们的方法在联合基因分析方面的能力仅限于估计和比较基因相关性。 鉴于许多生物途径在低维流形内运行，未来的增强功能应能够进行全面的联合分析，涉及两个以上的基因。 此外，Memento 假设检验框架的灵活性应能无缝促进这些适应。
Resource availability
Lead contact
主要联系人
Para_01
进一步的信息和资源、试剂的请求应直接联系主要联系人，Chun Jimmie Ye（jimmie.ye@ucsf.edu）。
Materials availability
材料可用性
Para_01
本研究未产生新的独特试剂。
Data and code availability
数据和代码的可用性
单细胞 RNA 测序数据已存入 GEO，并于出版日期起公开可用。登录号列在关键资源表中。
所有原始代码均已存放在 GitHub 上，并于出版之日公开。相关网址列在关键资源表中。
重新分析本文报告的数据所需的所有其他信息，可根据请求从主要联系人处获得。
Acknowledgments
Para_01
我们承认在撰写和编辑本手稿的过程中使用了人工智能工具。 这些工具帮助组织和精炼内容，提高了最终文档的清晰度和连贯性。 作者对所呈现研究的知识产权内容和准确性负全责。 C.J.Y. 受到 NIH 授予的 R01HG011239、R01AI136972 和 R01AI045073 以及 Chan Zuckerberg Initiative 和 Cancer Research Institute 的支持。 他还是 Chan Zuckerberg Biohub 和 Arc Institute 的研究员，并且是 Parker Institute for Cancer Immunotherapy (PICI) 的成员。 本文所呈现的研究还受益于 JDRF 和 NIH 的风湿病精准医学资源中心的支持。
Author contributions
Para_01
M.C.K. 和 C.J.Y. 构思了项目并撰写了论文。 M.C.K. 开发了 Memento 并进行了所有分析。 D.S.L. 和 E.G. 生成了 HTEC 数据。 A.L. 和 R.G. 生成了 Perturb-seq 数据。 A.T. 和 P.E.G.-N. 贡献了 Memento 在 CELLxGENE Discover 中的实现。 V.N.、E.S. 和 A.M. 提供了宝贵的统计和实验反馈。
Declaration of interests
Para_01
C.J.Y. 是 DropPrint Genomics（现为 ImmunAI）和 Survey Genomics 的创始人，并持有这两家公司的股权；是 Related Sciences 和 ImmunAI 的科学顾问委员会成员，并持有这两家公司的股权；同时是 Maze Therapeutics 的顾问，并持有该公司的股权。 C.J.Y. 获得了来自 Chan Zuckerberg Initiative、Chan Zuckerberg Biohub、Arc Institute、Parker Institute for Cancer Immunotherapy、Genentech、BioLegend、ScaleBio 和 Illumina 的研究支持。 A.M. 是 Site Tx、Arsenal Biosciences、Spotlight Therapeutics 和 Survey Genomics 的联合创始人；担任 Site Tx、Spotlight Therapeutics 和 Survey Genomics 的董事会成员；并担任 Site Tx、Arsenal Biosciences、Cellanome、Spotlight Therapeutics、Survey Genomics、NewLimit、Amgen 和 Tenaya 的科学顾问委员会成员。 A.M. 持有 Arsenal Biosciences、Site Tx、Cellanome、Spotlight Therapeutics、NewLimit、Survey Genomics、Tenaya 和 Lightcast 的股票，并从 Site Tx、Arsenal Biosciences、Cellanome、Spotlight Therapeutics、NewLimit、Gilead、Pfizer、23andMe、PACT Pharma、Juno Therapeutics、Tenaya、Lightcast、Trizell、Vertex、Merck、Amgen、Genentech、GLG、ClearView Healthcare、AlphaSights、Rupert Case Management、Bernstein 和 ALDA 收取了费用。 A.M. 是 Offline Ventures 的投资者和非正式顾问，并且是 EPIQ 的客户。 Marson 实验室获得了来自 Parker Institute for Cancer Immunotherapy、Emerson Collective、Arc Institute、Juno Therapeutics、Epinomics、Sanofi、GlaxoSmithKline、Gilead 和 Anthem 的研究支持，并从 Genscript 和 Illumina 获得了试剂。
STAR★Methods
Key resources table
关键资源表
Experimental model and study participant details
实验模型和研究参与者详情
Human tracheal epithelial cells
人支气管上皮细胞
Para_01
人气管上皮细胞根据既定协议从已故器官捐赠者处收集。 冷冻细胞悬液被重新激活，并在上皮生长培养基（EGM）中培养，该培养基由F-12营养混合物（Gibco）和Dulbecco改良Eagle培养基（Invitrogen）按3:1（体积比）混合，含有5%胎牛血清（Gibco）、0.4 μg/mL氢化可的松（Sigma-Aldrich）、5 μg/mL胰岛素（Sigma-Aldrich）、8.4 ng/mL霍乱毒素（Sigma-Aldrich）、10 ng/mL表皮生长因子（Invitrogen）、24 μg/mL腺嘌呤（Sigma-Aldrich）和10 μM Y-27632（Enzo Life Sciences）。 EGM每周更换三次，直到培养皿长满，此时用0.25%胰蛋白酶处理30分钟进行传代。 对于气液界面培养，扩增的基底细胞以每6.5 mm Transwell插件（Corning 3470）50,000个细胞的密度接种到涂有人胎盘胶原蛋白（Sigma-Aldrich）的插件上，并按照制造商的说明使用Pneumacult ALI（StemCell）培养21-28天。
Study subjects and genotyping for Perturb-seq
Perturb-seq 的研究对象和基因型分析
Para_01
我们的样本参加了 PhenoGenetic 研究（年龄 18 至 56 岁，平均 29.9 岁），作为 Immvar 队列的一部分，这些样本在大波士顿地区招募。 每位捐赠者均签署书面同意书参与研究，并且健康，没有任何炎症疾病、自身免疫疾病、慢性代谢障碍或慢性感染障碍的历史。 我们在 OmniExpressExome 芯片上对 56 名高加索样本进行了基因分型，并排除了 2080 个呼叫率 < 90% 的 SNP（占总数的 0.22%），1521 个 Hardy Weinberg P < 0.0001 的 SNP（占 0.16%）和 259,860 个 MAF < 0.1 的 SNP（占 27.04%），共计 960,919 个 SNP 中的这些 SNP。 使用 Michigan Imputation Server 和 Haplotype Reference Consortium Panel Version r1.1 对这些基因型进行推算后，共获得了 5,324,560 个 SNP，然后从中选择了我们九位捐赠者的子集。
Method details
方法详情
Inteferon stimulation of HTECs
干扰素对HTECs的刺激
Para_01
从第27天开始，在收获前0、24、39、42和45小时分别加入干扰素刺激（IFN-β: 10 ng/ml, IFN-α2: 10 ng/ml, IFN-γ: 10 ng/ml, IFN-λ2: 10 ng/ml），最终时间点分别为3、6、9、24和48小时。 在收获当天，吸出基础培养基，并用PBS冲洗基底室和顶室两次。 冲洗两次后，向基底室和顶室中加入胰蛋白酶-EDTA（0.25% Fisher 编号 25200072）（基底室300 μl，顶室100 μl），并在37°C下孵育30分钟，每10分钟混匀一次。 用300 μl维持培养基终止胰蛋白酶化，转移到1.5 ml离心管（Eppendorf 编号 022431021），并在4°C下以350xg离心5分钟。 将细胞重悬于94 μl细胞染色缓冲液（Biolegend 编号 420201）中，并在冰上用5 μl TruStain FcX（Biolegend 编号 422302）封闭10分钟。 封闭后的细胞在冰上用1 μl Biolegend Totalseq-B 哈希标签（Biolegend Totalseq-B 哈希标签 1-11）染色30分钟。 用1 ml细胞染色缓冲液终止染色，并在4°C下以300xg离心5分钟，然后再用1 ml细胞染色缓冲液洗涤两次。 将细胞重悬于100 μl 0.05% BSA的PBS中，并通过Countess II（Fisher 编号 A27977）计数。 计数后的细胞等量分为两池，并在4°C下以300xg离心5分钟。 细胞通过100 μm滤器（Corning 编号 431752）过滤，进行最终计数，每池加载到两个10x 3'v3通道上。 按照10x 3'v3用户指南准备文库。 样品在NovaSeq S4的三个通道上测序。
Regulator target identification and CROP-seq library generation
调节因子靶标识别和CROP-seq文库生成
Para_01
我们的文库包含140个靶向调节因子（转录因子和RNA结合蛋白），每个调节因子有2个sgRNA。 每个调节因子都是根据95和105名健康供体激活的CD4+ T细胞的基因表达和可及性数据无偏选择的。 为了使用RNA-seq数据获取高表达的调节因子，我们进行了TMM归一化，并取了高表达基因的上四分位数，然后筛选出那些是调节因子的基因。 为了使用ATAC-seq数据获取具有高可及性结合位点的调节因子，我们在激活的可及染色质区域中富集了HOMER数据库中的所有结合位点。 我们取了高表达调节因子和高可及性结合位点的并集，总共得到140个调节因子（图1B）。
Para_02
用于克隆 CROP-Seq 文库的主干质粒是 CROPseq-Guide-Puro，购自 Addgene（Addgene. 质粒 #86708）。 我们为选定的 140 个调控因子中的每一个使用了来自 Brunello 文库的两个 sgRNA 寡核苷酸序列。 sgRNA 文库的寡核苷酸从 Integrated DNA Technologies (IDT) 购买，并使用 Datlinger 等人描述的方法克隆到 CROPseq 质粒主干中。 慢病毒使用 UCSF ViraCore 生产。
SLICE experiment and sequencing
SLICE 实验和测序
Para_01
原代人 CD4+ T 细胞通过使用 EasySep 人 CD4+ T 细胞分离试剂盒（STEMCELL，货号 #17952）从外周血单核细胞（PBMCs）中通过磁性负选择分离。 细胞在 X-Vivo 培养基中培养，该培养基由 X-Vivo15 培养基（Lonza，货号 #04-418Q）组成，含有 5% 胎牛血清、50 mM 2-巯基乙醇和 10 mM N-乙酰-L-半胱氨酸。 在分离当天（第 1 天），细胞在没有刺激的情况下在培养基中休息 24 小时。 在分离后的第二天（第 2 天），细胞用 ImmunoCult 人 CD3/CD28 T 细胞激活剂（STEMCELL，货号 #10971）和 50 U/mL 的 IL-2 刺激。 刺激后 24 小时（第 3 天），1 μL 的慢病毒直接加入到培养的 T 细胞中，并轻轻混合。 24 小时后（第 4 天），细胞被收集、离心，并用 PBS 洗涤两次。 然后，细胞被重悬于 Lonza 电转缓冲液 P3（Lonza，货号 #V4XP-3032）中。 将 Cas9 蛋白（MacroLab, Berkeley，40 mM 储备液）以 1:10 的体积比加入到细胞悬浮液中。 细胞被转移到 96 孔电转杯板（Lonza，货号 #VVPA-1002）中，使用 Lonza Nucleofector 96 孔穿梭系统和脉冲代码 EH115（Lonza，货号 #VVPA-1002）进行核转染。 电转后立即向每个电转孔中加入预热的培养基，并将 96 孔板置于 37 度下 20 分钟。 然后，细胞被转移到含有 50 U/mL IL-2 的 X-Vivo 培养基中，以 1 × 10^6 细胞/mL 的密度在适当的组织培养容器中培养。 两天后，在培养基中加入 1.5 μg/mL 的嘌呤霉素进行选择。 每两天扩增一次细胞，添加新鲜的含 50 U/mL IL-2 的培养基，保持细胞密度为 1 × 10^6 细胞/mL。 在实验的最后一天（第 13 天），使用 Vi-CELL XR 计数来自九位供体的细胞，并等量混合以获得最终的 180,000 个细胞，体积为 60 μL 的 PBS。 混合后的细胞由 UCSF 人类遗传学研究所（IHG）基因组核心实验室使用 10X Chromium 单细胞 v2（PN-120237）的 16 个孔进行处理，每个孔单独索引，按照制造商的协议操作。 最终文库在 Nova-seq 上进行了测序，总读数为 6.7B 个。 为了最大化检测细胞中 sgRNA 的概率，我们进一步扩增并测序了 10X cDNA 文库中的 sgRNA 转录本，直至接近饱和，如前所述（98%）。
Quantification and statistical analysis
量化和统计分析
Modeling scRNA-seq as a hypergeometric sampling process
将单细胞 RNA 测序建模为超几何抽样过程
Para_01
单细胞 RNA 测序（scRNA-seq）中的测量噪声可以归因于几乎所有工作流程中共有的至少三个分子生物学过程的低效率：1）仅捕获部分表达转录本进行逆转录（RT）为 cDNA，2）每个聚合酶链反应（PCR）循环中仅扩增部分 cDNA 分子，3）仅对部分扩增的 cDNA 进行测序。 尽管独特分子标识符（UMIs）的发展在很大程度上消除了建模 PCR 引入的噪声的需求，但由不完全转录本捕获用于 RT 和测序过程中不完全 cDNA 采样引起的噪声仍然存在，最终导致观察到的计数分布减弱。
Para_02
我们使用一个灵活的分层模型来建模从单细胞 RNA 测序 (scRNA-seq) 获得的计数数据，该模型明确考虑了表达转录本计数的生成过程以及通过大规模并行 scRNA-seq 方法对 mRNA 分子的采样。 正如正文中所述，我们对 scRNA-seq 采样过程的完整模型可以概括如下： 细胞 c 中的表达转录本计数 Z_c ∼ P_Z, 细胞 c 的总转录本计数 N_c = 1^T Z_c, 细胞 c 中的归一化转录本计数 X_c = Z_c / N_c, 细胞 c 中观察到的转录本计数 Y_c ∼ MultiHG(Z_c, N_c, q N_c) = MultiHG(X_c N_c, N_c, q N_c),
Para_03
q 是表示在观察到的单细胞 RNA 测序实验中最终被计数为唯一分子标识符（UMIs）的表达转录本比例的随机变量。 在我们上面讨论的大多数单细胞 RNA 测序工作流程中的噪声来源时，它既考虑了逆转录采样效率，也考虑了测序过程中转录本的采样。 在极端情况下，如果文库被测序至饱和，则 q 减少为逆转录采样效率；然而，在大多数实验中，文库并未测序至饱和，而是测序到已知比例的独特分子。 通过广泛的模拟，我们证明这种复合噪声过程可以通过使用一个多重超几何过程来很好地近似，该过程使用 E[q] 的值，该值是逆转录采样效率（针对特定实验技术可用）和测序采样效率（从预处理管道如 CellRanger 可用）的乘积（图 S1A 和 S1B）。
Para_04
然后我们用多元超几何分布来建模mRNA捕获过程。 给定(Z₁, Z₂, Z₃, …, Zₖ)成分（即基因），总计数N = ∑(i=1到c) Zᵢ，和样本数量n ∈ {0, 1, …, N"}，多元超几何分布的概率质量函数(PMF)由下式给出：(公式1) pₘᵤₗₜᵢₕᵍ(Y; Z₁, Z₂, …, Zₖ, N, n) = ∏(i=1到G) (Zᵢ 选择 Yᵢ) / (N 选择 n)
Para_05
在先前的研究中，完整的超几何处理通过一系列近似简化，从超几何模型到泊松模型： Y_c ∼ MultiHG(Z_c, N_c, qN_c)，观察到的转录本计数 Y_c ∼ Multinomial(\frac{Z_c"}{N_c"}, qN_c)，观察到的转录本计数 Y_{cg"} ∼ Bn(\frac{Z_{cg"}}{N_c"}, qN_c)，观察到的转录本计数 Y_{cg"} ∼ Poi(\frac{Z_{cg"}}{N_c"}qN_c)，观察到的转录本计数 Y_{cg"} ∼ Poi(qZ_{cg"})，观察到的转录本计数
Para_06
Y_cg 是向量 Y_c 中的一个元素，因为泊松模型认为每个基因的采样是独立的。 正如我们在接下来的部分中讨论的那样，当 q 非常小（接近 0）时，完整的超几何处理和泊松简化得到的估计器非常相似，但随着 q 值的增加，两者之间的差异会变得更大，因为单细胞 RNA 测序实验流程的改进。
Method of moments estimation of expressed transcript counts
表达转录本计数的矩估计方法
Para_01
我们将从回顾 Zhang 等人在确定单细胞 RNA 测序实验最佳测序深度时首次提出的泊松估计器的推导开始。 首先，回顾之前介绍的单细胞 RNA 测序的泊松抽样模型，其中 Nc 表示每个细胞的总表达转录本，q 是总体抽样效率，Xcg 是真实的相对 mRNA 表达量，Ycg ~ Poi(qNcXcg)。
Para_02
对于泊松变量 A ∼ Poi(λ)，A 的矩为 E[A] = λ 和 E[A^2] = λ^2 + λ。 同样地，对于我们的模型，我们可以写出给定其他变量 q、N_c 和 X_c 的 Y_cg 的矩的方程。(方程 2) E[Y_cg | X_cg, N, q] = X_cg N_c q E[Y_cg^2 | X_cg, N, q] = X_cg^2 N_c^2 q_c^2 + X_cg q_c N_c E[Y_cgi Y_cj | X_cgi, X_cj, N, q] = E[X_cgi X_cj N_c^2 q^2 | X_cgi, X_cj, N, q] = X_cgi X_cj N_c^2 q^2
Para_03
将第一个矩方程代入第二个，我们得到：E[Y_cg^2 - Y_cg | X_cg, N, q] = X_cg^2 N_c^2 q^2
Para_04
这些方程导致了对 μ̂g,Poi、σ̂g,Poi² 和 σ̂gigj,Poi 的估计器，分别是 Xcg 的均值、方差和协方差，通过将所有单元格的矩平均化得到：(公式 4)μ̂g,Poi = Ê[Xcg] = 1/ncells ∑c Ycg / (Nc q)，σ̂g,Poi² = Ê[Xcg²] − Ê[Xcg]² = 1/ncells ∑c (Ycg² − Ycg) / (Nc² q²) − (1/ncells ∑c Ycg / (Nc q))²，σ̂gigj,Poi = Ê[Xcgi Xcgj] − Ê[Xcgi] Ê[Xcgj] = 1/ncells ∑c (Ycgi Ycgj / (Nc² q²)) − (1/ncells ∑c Ycgi / (Nc q))(1/ncells ∑c Ycgj / (Nc q))
Para_05
现在，让我们考虑完整的多变量超几何模型，Yc∼=MultiHG(XcNc,Nc,qNc)。 对于随机向量 A∼MultiHG(K,N,n)，A 的矩为： E[Ai]=nKiN E[Ai2]=n(N−n)/(N−1) * Ki/N * (1−Ki/N) + n2Ki2/N2 E[AiAj]=−n(N−n)/(N−1) * KiKj/N2 + n2KiKj/N2
Para_06
我们可以再次写出矩方程，这次是多元超几何模型的。 E[Ycg|Xcg,Nc,q] = qNc * (Xcg / Nc) = Xcg * Nc * q E[Ycg^2|Xcg,Nc,q] = qNc * ((Nc - qNc) / (Nc - 1)) * (Xcg / Nc) * (1 - (Xcg / Nc)) + q^2 * Nc^2 * (Xcg^2 / Nc^2) ≈ qNc * (1 - q) * Xcg * (1 - Xcg) + q^2 * Nc^2 * Xcg^2 = Xcg^2 * (Nc^2 * q^2 - Nc * q * (1 - q)) + Xcg * Nc * q * (1 - q) E[YcgiYcgj|Xcgi,Xcgj,N,q] = -qNc * ((Nc - qNc) / (Nc - 1)) * (Xcgi * Xcgj / Nc^2) + q^2 * Nc^2 * (Xcgi * Xcgj / Nc^2) ≈ q^2 * Nc^2 * Xcgi * Xcgj - q * (1 - q) * Nc * Xcgi * Xcgj = Xcgi * Xcgj * (Nc^2 * q^2 - Nc * q * (1 - q))
Para_07
将第一个矩方程代入第二个矩方程，我们得到：E[Ycg2−(1−q)Ycg|Xcg,N,qc]=Xcg2(Nc2q2−Ncq(1−q)) E[Ycg2−(1−q)Ycg|Xcg,N,qc] = Xcg2(Nc2q2−Ncq(1−q))
Para_08
用于推导第二和第一对数矩的近似假设是 Nc 远大于 1。 对于大多数哺乳动物细胞，表达的转录本数量约为 10^5，这些近似应该成立。 类似于基于泊松模型的估计器，我们可以从完整的多变量超几何模型中推导出基于这些矩方程的估计器：（公式 8）μ̂g,HG = Ē[Xcg] = (1/n_cells) ∑c Ycg / (Nc q)，σ̂g,HG² = Ē[Xcg²] - Ē[Xcg]² = (1/n_cells) ∑c (Ycg² - Ycg (1-q)) / (Nc² q² - Nc q (1-q)) - ((1/n_cells) ∑c Ycg / (Nc q))² ≈ (1/n_cells) ∑c (Ycg² - Ycg (1-q)) / (Nc² q²) - ((1/n_cells) ∑c Ycg / (Nc q))²，σ̂gigj,HG = Ē[Xcgi Xcgj] - Ē[Xcgi] Ē[Xcgj] = (1/n_cells) ∑c Ycgi Ycgj / (Nc² q² - Nc q (1-q)) - ((1/n_cells) ∑c Ycgi / (Nc q)) ((1/n_cells) ∑c Ycgj / (Nc q)) ≈ (1/n_cells) ∑c Ycgi Ycgj / (Nc² q²) - ((1/n_cells) ∑c Ycgi / (Nc q)) ((1/n_cells) ∑c Ycgj / (Nc q))
Para_09
最后，为了完整性，我们写出均值、方差和协方差的朴素估计器。
（公式 9）μ̂g,naive = Ê[Xcg] = 1/ncells ∑c Ycg / (Ncq)
σ̂g,naive^2 = Ê[Xcg^2] − Ê[Xcg]^2 = 1/ncells ∑c Ycg^2 / (Nc^2 q^2) − (1/ncells ∑c Ycg / (Ncq))^2
σ̂gigj,naive = Ê[Xcgi Xcgj] − Ê[Xcgi] Ê[Xcgj] = 1/ncells ∑c Ycgi Ycgj / (Nc^2 q^2) − (1/ncells ∑c Ycgi / (Ncq))(1/ncells ∑c Ycgj / (Ncq))
到目前为止，均值和协方差的估计器在朴素、泊松和 HG 估计器之间是相同的。然而，方差的估计器，它对残差方差和相关性的测量有所贡献，是这三组估计器之间的关键差异。
重要的是，可以明显看出 HG 方差估计器包括了朴素和泊松估计器：
当 q 趋近于 0 时，σ̂²_g,HG = σ̂²_g,Poi；当 q 趋近于 1 时，σ̂²_g,HG = σ̂²_g,naive。
这些结果表明，当 q 较小时，总体采样效率较低，HG 估计器的行为与泊松估计器非常相似。
当 q 接近 1 时，假设单细胞 RNA 测序工作流程完美且我们捕获了所有表达的转录本，HG 估计器收敛到朴素估计器，因为没有噪声过程。
随着单细胞 RNA 测序工作流程的改进和 q 变大，HG 估计器成为 Zhang 等人提出的估计器的推广，以适应具有不同 q 值的不同类型实验工作流程。
我们还讨论了 q 在细胞之间不是常数的情况。
在推导上述估计器时，假设 q 是一个已知的常数，我们不需要为每个细胞单独估计它。
然而，对于某些单细胞 RNA 测序技术，在测序未饱和的情况下，q 实际上可能是围绕其均值 E[q] 的分布。
实验上，我们可以通过使用 spike-in RNA 控制来实际测量每个细胞的 q 值来缓解这一问题。
因为 q 不出现在泊松估计器中，即使可以测量每个细胞的 q 值，也无法显式地考虑 q 的变异性。
通过这里推导的超几何估计器，我们可以简单地用每个细胞测量的 qc 值替代上述的 q。
Improving mean estimation with Good-Turing method
使用 Good-Turing 方法改进均值估计
Para_01
我们在这一节中的推导展示了不同的朴素、泊松和超几何估计器来估计方差和相关性，但均值估计器是相同的。 我们试图进一步改进均值估计，特别是在细胞数量较少的情况下，这可能出现在许多生物样本和扰动组合的实验中。 保持我们的超几何模型，我们从 Good-Turing 频率估计中获得灵感，该方法可以用于估计以前未见过的物种的频率。 Good-Turing 估计指出，如果基因 i 的转录本在包含总共 N 个转录本的池中出现了 r 次，并且出现 r 次的基因数量为 n_r，我们应该将基因 A 的频率估计为：(公式 11) 1/N * (r+1) * n_{r+1} / n_r
Para_02
我们将此方程应用于单细胞数据的生物样本层面，最终得到平均估计值：(公式 12) rg = ∑c Ycg, 基因 g 在样本中的计数；nr = ∑g 1(rg = r)，计数为 r 的基因数量；rg* = (rg + 1)nr+1 / nrg；μ̂g,memento = Ē[Xcg] = (1 / ∑c Nc) rg*
Estimating cell sizes by trimming variable genes
通过修剪可变基因来估算细胞大小
Para_01
上述 HG 估计方程中出现的 Ncqc 值指的是细胞大小，它作为每个细胞的归一化因子。 这些常数确保即使捕获的转录本比例在细胞之间有所不同，估计也不会受到这种技术噪声的影响。 我们可以将 Ncqc 分解为两个部分：一个常数 numi 和 γc，使得 Ncqc = numi * γc。 最简单的估计 γc 的方法是首先计算 numi = (1 / ncells) * Σc (1^T * Yc)，然后设置 γc = (1 / numi) * (1^T * Yc)，进行总计数归一化。
Para_02
这是张等人工作中提出的泊松估计器如何估计细胞大小的方法。 在 Memento 中，我们通过首先计算数据集中所有细胞的残差方差，并去除表现出高变异性的基因，提供了一种替代方法。 这种方法假设数据集中的大多数基因不应该差异表达，最不变化的基因适合作为归一化的基准。 使用非差异表达基因的想法已经在其他方法中被采用，例如其他方法。 默认情况下，Memento 使用最不变化的 10% 的基因。 在去除变异基因后形成基因集 G* 后，我们使用以下公式计算 γc：γc = (δ + ∑(g ∈ G*) Ycg) / (δ + 1/ncells ∑c ∑(g ∈ G*) Ycg)
Para_03
这里的 δ 值作为估计细胞大小的正则化因子；当这个值较高时，表示数据集不需要大小因子归一化（采样在细胞间确实恒定，例如在测序饱和时）。 默认情况下，Memento 使用细胞 c 的 median(∑g∈G*Ycg) 作为 δ 值。
Para_04
需要注意的是，文献中存在更复杂的归一化方法。 Memento 可以轻松地将这些计算细胞大小的替代方法整合到其流程中。
Computing the residual variance
计算残差方差
Para_01
单细胞 RNA 测序（scRNA-seq）数据中的均值和方差通常高度相关，测量表达变异时必须考虑这种相关性。 BASiCs 通过在拟合均值和过度离散参数之间进行多成分非线性回归来考虑这种依赖关系。 我们没有先拟合负二项分布然后从过度离散参数中回归出均值，而是直接取估计的真实均值和方差，并拟合一个简单的多项式回归。 我们对给定细胞群的所有基因使用单个拟合多项式（默认次数为 2），该细胞群由细胞类型、实验条件或批次定义。 我们发现，即使这个简单的回归也能够很大程度上消除 scRNA-seq 数据中存在的均值-方差依赖性。
Efficient bootstrapping by exploiting data sparsity
通过利用数据稀疏性实现高效的引导
Para_01
通常，生成置信区间和计算假设检验的 p 值需要对数据分布和估计器本身做出某些假设。 例如，为了计算线性回归模型系数的 p 值，我们通常假设数据是正态分布的，并且系数的抽样分布也是正态的。 在单细胞 RNA 测序（scRNA-seq）的背景下，我们的估计器可以在不对表达转录本计数的分布做出任何假设的情况下，测量平均值、变异性和基因相关性。 然而，在不对数据本身和估计值的抽样分布做出任何假设的情况下，很难计算出我们估计器的解析置信区间。
Para_02
自助法是一种估计任意统计量抽样分布的方法，无需对数据生成过程做出大量假设。 在 Memento 中，我们提出了一种策略，以极其高效的方式在单细胞 RNA 测序（scRNA-seq）数据中执行自助法。 具体来说，在一个包含 N 个细胞 x1, x2, x3, …, xN 的单基因数据集中，我们可以将每个观察值的出现次数建模为多项分布 Multinomial(N, 1/N, …, 1/N)。 如果有 K 个独特的计数，每个计数有 nk 个细胞，我们可以将重新采样的分布重写为 Multinomial(N, n1/N, …, nK/N)。
Para_03
在考虑标准化转录本丰度时，我们必须考虑到每个细胞中的总转录本数量（Nc）。 虽然这从技术上会在每个细胞中创建不同的 Nc，使得我们的方案变得不那么有用，但将不同细胞的 Nc 分类到少量离散区间内的策略可以很好地近似参数的真实自助分布。 通过模拟，我们展示了随着区间数量的增加，真实自助分布和近似自助分布几乎相同（图 S2D 和 S2E）。
Hypothesis testing and extension to account for replicates in multiplexed scRNA-seq experiments
假设检验及其在多重 scRNA-seq 实验中考虑重复样本的扩展
Para_01
考虑一个包含两组细胞 A 和 B 的场景，我们计算了每组的感兴趣参数 t，并计算了它们的差异 Δt。 t 取决于我们想要执行的测试类型；我们将计算平均值、残差方差和相关性，分别用于测试平均值、变异性和共表达的差异。 然后，我们通过 B 次迭代进行自助法抽样，以生成测试统计量 Δt 的抽样分布，从 Δt1 到 ΔtB。 如果我们希望检验备择假设 H1: Δt ≠ 0 与零假设 H0: Δt = 0，我们首先通过从 Δt1, …, ΔtB 中减去 Δt 来生成零分布，形成 Δt1*, …, ΔtB*，类似于 Efron 和 Tibshirani 提出的策略。 然后，我们可以计算该检验的实现显著性水平（ASL），公式如下：ASL = {2/B * Σ(i=1 到 B) 1(Δt > Δti*) 如果 Δt ≥ 0；2/B * Σ(i=1 到 B) 1(Δt < Δti*) 如果 Δt < 0
Para_02
生成带有重复样本（例如不同个体）的单细胞 RNA 测序数据的趋势日益增加，特别是在多路复用工作流程中。 考虑一个具有两个条件和 n 个重复样本的实验。 然后，我们提出一个元分析框架，首先将细胞分为 2n 组，并使用 2n 个观察值进行元回归： [lnμ1, lnμ2, ⋮, lnμ2n-1, lnμ2n], [lnσ˜1, lnσ˜2, ⋮, lnσ˜2n-1, lnσ˜2n], [ρ1, ρ2, ⋮, ρ2n-1, ρ2n] ~ β[W1, W2, ⋮, W2n-1, W2n] + α[1, 1, ⋮, 1, 1] 其中 μi、σ˜i 和 ρi 分别表示第 i 个重复样本中估计的均值、残差方差和相关性，Wi 表示条件。 然后，我们可以对回归系数进行 B 次自助法抽样，以得出原始统计量 βˆ 和自助法统计量 βˆ1, …, βˆB。 类似于非重复样本的情况，我们可以通过从 βˆ1, …, βˆB 中减去 βˆ 来生成零分布 βˆ1*, …, βˆB*。 我们还可以进一步计算 ASL： ASL = {2/B ∑i=1B 1(βˆ > βˆi*) if βˆ ≥ 0, 2/B ∑i=1B 1(βˆ < βˆi*) if βˆ < 0}
Para_03
或者，零分布 βˆi∗ 可以近似为正态分布 N(0, σ∗²)，显著性水平可以计算为 2(1−Φ(|βˆ|/σ∗))。
Para_04
该框架可以轻松扩展到包含许多协变量，包括批次变量和感兴趣变量之间的交互作用，通过在方程13中的模型中引入额外的协变量来实现，即除了处理变量W之外提供额外的列。 细胞群组层面的信息，如年龄、性别、基因型，可以类似于将其纳入广义线性模型的方式进行纳入。 目前，Memento 还不处理单细胞层面的协变量。
Para_05
作为一个技术性的补充，我们注意到计算 ASL 的过程假设了感兴趣的检验统计量的抽样分布是平移不变的。 通过广泛的模拟，我们确认对于我们在 Memento 中考虑的检验统计量，此过程在零假设下产生了校准良好的结果（图 S2F）。 如果使用自定义检验统计量，检查假设检验结果的校准非常重要。 Memento 还提供了选项，假设效应大小的抽样分布是具有未知方差的正态分布，该方差使用自助法进行估计，这有助于加快假设检验的速度。 在这项工作中，这种近似仅用于分析自然变异的影响（图 5）。
Pre-processing the rIFNB1 PBMC dataset
预处理 rIFNB1 PBMC 数据集
Para_01
我们使用了原始聚类和 rIFNB1 数据集的 tSNE 可视化，该数据集已存放在 Gene Expression Omnibus 中，访问编号为 GSE96583。 关于此数据集预处理的更多细节可以在原始论文中找到。
Para_02
对于所有分析，我们选择了平均观察表达量 E[Ycg] = 0.07 的基因，这是该实验的可靠性限制。 关于可靠性限制的更多细节可以在 Zhang 等人的研究中找到。 该值是根据 10X Chromium V1 报告的 UMI 捕获效率以及该实验的测序饱和度（约为 90%）计算得出的。
Extracting mean and variance from scRNA-seq data for simulation
从单细胞 RNA 测序数据中提取均值和方差用于模拟
Para_01
我们使用了 PBMC 数据集作为模拟的基础。 我们希望模拟单细胞 RNA 轮廓，使其均值和方差的分布处于 scRNA-seq 的现实范围内。 为了实现这一点，我们使用 Memento 估计器从 CD4+ T 细胞中估计了 5000 个最高表达基因的均值、方差和相关性。 这些值随后被设定为模拟的真值参数。 我们使用了两组真值参数——一组是从未刺激的细胞中估计的（munstim, vunstim），另一组是从用 IFN-B 刺激的细胞中估计的（mstim, vstim）。
Simulating transcriptomes with given means, variances, and gene-gene correlations
模拟具有给定均值、方差和基因-基因相关性的转录组
使用真实相关参数生成均值为零、方差为一的高斯样本
通过取每个点的逆高斯累积分布函数来计算copula
通过使用先前计算出的指定均值和离散向量评估负二项分布的百分位点函数来生成边缘分布。
Para_01
该过程使用高斯Copula方法从具有指定相关矩阵和负二项边缘分布的联合分布中生成多变量样本。 请注意，Memento 不对基础分布做任何假设，这里使用负二项分布是为了与过去模拟单细胞RNA测序数据的策略保持一致。
Para_02
在模拟出"真实"的转录组后，我们使用超几何分布以总体捕获效率 q（结合文库制备和测序的采样）对转录本进行采样。
Comparing Memento, BASiCS, and scHOT for estimation
比较 Memento、BASiCS 和 scHOT 用于估计
Para_01
为了生成图2A，我们使用了从实际单细胞RNA测序数据集中估计的munstim和vunstim，以及从scikit-learn包中的make_psd函数采样的相关矩阵C，同时改变了总体捕获效率qreal。 我们使用Memento（超几何分布）、Memento（q=0）和朴素估计器来估计均值和相关性。 我们使用Memento（超几何分布）、Memento（q=0）、朴素和BASiCS估计器来估计方差。 我们通过使用BASiCS输出的离散度估计，根据负二项分布的均值-方差关系计算方差。 由于无法直接计算smFISH数据中的残差方差，我们在这次分析中用变异系数代替残差方差。 对于每个总体捕获效率值，这个过程重复了20次以生成置信区间。 我们使用了10个细胞的模拟来估计均值，100个细胞的模拟来估计变异性及基因相关性。
Para_02
对于均值、朴素、Memento（q=0）（泊松估计器）和BASiCS估计是相同的。BASiCS不适用于基因相关性。
Para_03
此外，为了研究基因的平均表达如何影响方差估计的准确性，我们将真实的模拟均值与方差估计的平均误差进行了比较。
Simulating genes with differential mean, variability, and coexpression
模拟具有不同均值、变异性和共表达的基因
Para_01
为了在现实和复杂的实验环境中比较 Memento 在差异表达方面的性能与其他方法，我们使用了具有层次生成的层次模拟。
Para_02
在这个模拟中，保留了150个基因的均值、变异性和相关性，其余基因则去除了这些差异（见STAR方法）。 对于均值差异的模拟，我们创建了一个包含生物学重复的数据集，以模拟多重实验设计。
Para_03
对于差异均值，我们首先计算 Δm = log(mstim) − log(munstim)。 为了指定真实差异表达基因，我们将 Δm 中小于 0.1 的任何元素设置为 0。 我们模拟了包含 2 个重复的数据，创建了 4 组细胞：未刺激重复 1，刺激重复 1，未刺激重复 2，未刺激重复 2。 我们生成了四组均值向量：m1,unstim = N(munstim, 0.25)，m2,unstim = N(munstim, 0.25)，m1,stim = m1,unstim + Δm + N(0, 0.25)，m2,stim = m2,unstim + Δm + N(0.25)。
Para_04
这些均值向量表示跨重复样本（如个体）存在的基线变异和异质治疗效应（来自不同重复样本的细胞可能不会以相同的方式响应）。 然后，我们模拟了不同数量的细胞（1000、1000、1100、1100），以模拟每个重复样本的不同样本量，使用在《具有给定均值、方差和基因-基因相关性的转录组模拟》中描述的程序。 对于差异均值模拟，我们将所有方差设置为 vunstim，并且不诱导基因之间的相关性。
Para_05
对于差异变异性，我们首先计算了 Δv = log(vstim) - log(vunstim)。 为了指定真实差异变异基因，我们将 Δv 中小于 0.1 的任何元素设置为 0。 我们模拟了包含 2 个重复的数据，创建了 4 个细胞组：未刺激重复 1、刺激重复 1、未刺激重复 2、未刺激重复 2。 我们生成了四组方差向量：v1,unstim = N(vunstim, 0.25)，v2,unstim = N(vunstim, 0.25)，v1,stim = v1,unstim + Δv + N(0, 0.25)，v2,stim = v2,unstim + Δv + N(0.25)。
Para_06
这些方差向量代表了重复样本（如个体）之间存在的基线变异和异质性治疗效应（来自不同重复样本的细胞可能不会以相同的方式响应）。 然后，我们模拟了不同数量的细胞（500、500、700、700），以模拟每个重复样本的不同样本量，使用在‘具有给定均值、方差和基因-基因相关性的转录组模拟’中描述的程序。 对于差异变异性模拟，我们将均值向量设置为与模拟差异均值相同的方式。
Para_07
对于差异相关性，我们采用了与差异均值和变异相似的方法，但通过从 scikit-learn 模型的 make_psd 函数生成的两个随机相关矩阵中相减来生成 Δcorr。
Para_08
类似于用于比较估计性能的模拟，我们使用超几何分布对"真实的"转录组中的转录本进行采样，总体捕获效率为 q（结合了文库制备和测序的采样）。
Comparing DE methods, BASiCS, and scHOT for differential mean, variability, and correlation
比较DE方法、BASiCS和scHOT在差异均值、变异性和相关性方面的表现
Para_01
为了将 Memento 与已建立的差异表达方法进行比较，我们使用了 Squair 等人（2022）使用的相同参数。 对于 BASiCS，我们在不加内标模式和回归模式下运行了该方法。 对于 scHOT，我们使用默认参数，并使用二进制掩码来计算"伪时间"上的参数。
Para_02
我们在不同细胞数量下进行了20次重复模拟，以生成图2C。
Para_03
对于均值、朴素、q=0（泊松估计器）和 BASiCS 估计是相同的。BASiCS 不适用于基因相关性，SCVI/SAVER 不提供基因变异性的直接估计。
Clustering the HTEC transcriptomes
聚类 HTEC 转录组
Para_01
我们使用 scanpy 工具套件的默认值进行了过滤、归一化和聚类。 细胞类型是根据已知的 HTECs 标志基因手动识别的。
Para_02
类似于 rIFNB1 数据集，我们选择了平均观察表达量 E[Ycg]=0.07 的基因，这是该实验的可靠性限制。
Clustering the correlation matrices for genes with differential mean expression
对具有差异均值表达的基因的相关矩阵进行聚类
Para_01
通过将每个刺激和时间点与未刺激对照组进行比较，使用 Memento 在纤毛细胞中鉴定了 DMGs。 使用 Memento 计算了 IFN-β 刺激条件下每个时间点的 DMGs 之间的相关性。 然后使用 sklearn Python 包中的 AgglomerativeClustering 函数对 6 小时时的时间点的相关矩阵进行了聚类。 选择了基因数量最多的前 4 个聚类进行绘图。
Identifying highly variable genes at baseline
识别基线水平下高变异基因
Para_01
我们使用 Memento 的单样本模式计算了转录组中每个基因的供体平均表达均值和变异性，这些基因的平均 UMI 计数大于 0.07。 然后，我们使用 EnrichR 进行基因集富集分析，以获得显著富集的基因集。
Estimation of cutting efficiency
切削效率的估算
Para_01
切割效率被估计为含有特定插入缺失（通过读取计数）的总 DNA 比例，该比例经过实验中具有特定 sgRNA 的细胞相对比例的归一化处理。
Para_02
需要注意的是，虽然分母和分子都不能是一个大于1的数字，但我们对用于归一化的具有特定引导序列的细胞比例的估计存在误差，导致少数引导序列的切割效率 > 1。
Visualizing gene regulatory networks
可视化基因调控网络
Para_01
为了生成图4E中的基因调控网络（GRNs），我们首先使用了一组调控因子与其差异均值表达基因的配对列表来定义一个二分图，然后在Cytoscape中进行了可视化。 接着，我们在同一个先前可视化的网络中（图4F）添加了由差异相关基因（DCGs）发现的相互作用调控因子之间的连接。
Identifying candidate interactions for differential correlation analysis
识别用于差异相关性分析的候选相互作用
Para_01
对于转录调节因子 TR，我们首先确定了所有 TR 作为转录激活子作用的 DMGs，在整个 KO 中 DM 系数小于 0。然后，我们在 WT 细胞中计算了每个 TR-DMG 对的相关性，并通过选择那些在 WT 中具有显著相关性（ρ > 0.1）的对来构建最终的 TR-DMG 对集合。对于这些 TR-DMG 对中的每一个，我们测试了不同 sgRNAs 靶向除 TR 以外的转录调节因子时的相关性差异。最终的相互作用集是通过过滤 FDR < 0.1 获得的。
Counting genes with shared TFBS for pairs of transcription factors
计算转录因子对之间具有共享转录因子结合位点的基因数量
Para_01
对于转录因子 TF1 和 TF2，我们首先在 ENCODE 数据集中的 ChIP-seq 数据中确定了它们的转录结合位点（TFBSs）。 然后，我们取已知基因转录起始位点（TSSs）的位置，并测量每个 TSS 附近最近的 TFBS 距离。 接下来，我们统计了在 TSS 附近一系列窗口大小内同时具有 TF1 和 TF2 的 TFBSs 的基因数量，窗口大小从 10 个碱基对到 100,000 个碱基对不等。 我们对随机选择的转录因子对以及通过差异相关性分析鉴定为相互作用的转录因子对进行了这一过程。
Assessing the tonic sensitivity of ISGs
评估ISGs的稳态敏感性
Para_01
我们使用了 Gough 等人的持续敏感性测量结果，其中作者比较了 IFNAR1-KO 和 WT 巨噬细胞中 ISGs 的表达。 这两组之间的倍数变化被定义为持续敏感性，这是我们用于图 3D 的数值。
eQTL discovery using pseudobulk approach and Memento
使用伪批量方法和Memento进行eQTL发现
Para_01
我们使用了由 Perez 等人生成的单细胞数据集，该数据集对系统性红斑狼疮（SLE）患者和健康对照组的外周血单核细胞进行了表征。 我们保留了该研究中使用的相同细胞类型分类。
Para_02
为了使用伪批量方法识别 eQTLs，我们首先通过将每个细胞的表达量归一化到每个细胞的总 UMI 数，计算所有个体中每个基因的平均值，并对每个均值进行 log(x+1) 转换来创建细胞类型和个体水平的伪批量。 我们过滤掉了单细胞数据集中 UMI 计数低于 0.01 的基因。
Para_03
我们分别对亚洲和欧洲人群进行了 Matrix eQTL 分析，使用了与 Perez 等人相同的基因型和协变量。 对于 Memento，我们也分别对这两个群体进行了测试，使用了相同的基因型和协变量。 我们使用了 Memento 的分层重采样模式。
Enrichment of eQTLs in ATAC peaks
eQTL在ATAC峰中的富集
Para_01
我们使用了 Perez 等人使用的同一组 ATAC 峰。 对于每个 SNP，我们标记了该细胞类型特异性 ATAC 峰是否覆盖了 SNP 的位置。 然后，我们使用 Wilcoxon 秩和检验比较了细胞类型峰内的 eQTL 候选基因的 p 值与细胞类型峰外的 eQTL 候选基因的 p 值。
Comparison of eQTLs with OneK1K cohort
与OneK1K队列的eQTL比较
Para_01
为了计算 ROC 曲线并在 5 中进行功效分析，我们将使用两种方法发现的 eQTL 与 Yazar 等人报道的 eQTL 进行了比较。 我们使用这个更大的数据集作为金标准来比较应用于 SLE 数据集的方法学。 具体来说，我们计算功效为能够在较小数据集中复制的 Onek1k 命中的比例。 我们通过打乱个体的基因型（同时保持个体-细胞分配不变）并计算 P < 0.05 的 SNP-基因对的比例来计算假阳性率。
Precomputation of estimates and standard errors in the CELLxGENE Discover database
CELLxGENE Discover 数据库中估计值和标准误差的预计算
Para_01
CELLxGENE Census 数据库提供了 RNA 表达计数，作为一个包含 M 个细胞和 N 个基因的 MxN 矩阵。 每个细胞都标注了元数据，指定了其细胞类型、数据集、组织、检测类型、供体、疾病、性别、发育阶段、种族和悬浮类型。 Census 数据是稀疏的，即如果某个细胞中某个基因的测量表达量不是正数，则不会显式存储。
Para_02
从这些人口普查数据中，我们根据细胞类型、数据集、组织、检测类型、供体、疾病、性别、发育阶段、种族和悬浮类型对细胞进行了分组。 然后，对于每个细胞组中的每种表达基因，我们使用上述相同的估计器和重采样策略计算了均值的对数、均值的标准误差、方差和方差的标准误差。
Para_03
这些预计算的值随后被保存，以便在无需重新计算这些估计器的情况下高效地执行重复的差异表达分析。 包含 3000 万个细胞的普查数据被缩减为 14 万个细胞组，因此其大小减少了两个数量级。
Hypothesis testing using precomputed standard errors
使用预先计算的标准误差进行假设检验
Para_01
为了计算两种不同细胞群体之间的差异表达，这两种细胞群体因特定处理而有所不同，可以通过使用两个不同的细胞注释值来过滤预计算数据，从而获取两组子集。 在计算差异表达时，所有剩余的细胞注释都被视为协变量。
Para_02
对于本文中呈现的细胞类型比较，我们使用了基本加权最小二乘法（WLS）来结合预先计算的均值和残差方差估计（作为响应变量）以及它们的标准误差（作为权重）。 为了在图6中的所有23个数据集中进行差异表达分析，我们使用了供体协变量作为独热编码变量，以及它们与细胞类型独热编码变量的交互项。
