生信算法 | 矩阵分解除了NMF，也可以试试这个 NatGenet 新发的 GBCD 算法

Basic Information

• 英文标题： Dissecting tumor transcriptional heterogeneity from single-cell RNA-seq data by generalized binary covariance decomposition
• 中文标题：通过广义二元协方差分解从单细胞RNA测序数据中剖析肿瘤转录异质性
• 发表日期：02 January 2025
• 文章类型：Technical Report
• 所属期刊：Nature Genetics
• 文章作者：Yusha Liu | Matthew Stephens
• 文章链接：https://www.nature.com/articles/s41588-024-01997-z

Abstract

Para_01

1. 通过单细胞RNA测序对肿瘤进行分析有可能识别出与癌症进展相关的反复出现的转录变异模式，并产生具有治疗意义的见解。
2. 然而，强烈的肿瘤间异质性可能会掩盖那些跨肿瘤共享的更为微妙的模式。
3. 在这里，我们介绍了一种统计方法，广义二元协方差分解（GBCD），来解决这个问题。
4. 我们展示了GBCD可以将转录异质性分解成可解释的组成部分——包括患者特异性、数据集特异性和与疾病亚型相关的共享组成部分。
5. 而且，在存在强烈的肿瘤间异质性的情况下，它能够比现有方法产生更可解释的结果。
6. 应用于胰腺导管腺癌的数据，GBCD产生了对现有肿瘤亚型的精炼分类，并识别出一种基因表达程序，该程序独立于肿瘤阶段和亚型预测不良生存。
7. 这种基因表达程序富集了参与应激反应的基因，并暗示整合应激反应在胰腺导管腺癌中起作用。

Main

Para_01

1. 通过单细胞RNA测序（scRNA-seq）技术剖析肿瘤继续揭示了肿瘤转录异质性，这种分辨率是整体转录组分析无法比拟的，并有可能重塑我们对癌症生物学的理解。
2. 特别是，通过分析转录变异模式，研究可以识别‘基因表达程序’（GEPs）——这些是一组基因，它们的转录倾向于协调地变化。
3. 一些GEPs可能是患者特有的，而另一些则可能在（部分）患者中共享。
4. 关键的是，共享的GEPs可能表征不同的分子亚型或细胞状态，因此可以提供关于疾病病因学的见解。
5. 此外，将这些GEPs与疾病进展或治疗反应联系起来，可能会提供具有医学意义的见解。

Para_02

1. 识别共享基因表达模式的一个基本挑战是，肿瘤通常表现出显著的肿瘤间异质性（一个肿瘤或患者与另一个之间的表达差异）1,6。
2. 因此，在涉及多个肿瘤的单细胞RNA测序数据集中，恶性细胞通常按患者聚类，不同患者的重叠很少。
3. 这种强烈的肿瘤间异质性往往掩盖了肿瘤之间共享的更微妙但仍然重要的异质性模式。
4. 事实上，正如我们在下面的模拟中所展示的，首先将所有肿瘤中的恶性细胞的表达数据结合起来，然后应用现成的统计方法（例如，非负矩阵分解（NMF）7）的综合分析，往往会识别出患者特异性的基因表达模式，但常常忽略了共享的基因表达模式。
5. 此外，由于个别单细胞RNA测序研究中肿瘤数量通常较少，多研究的联合分析可能是可取的，但这可能会引入数据集或批次效应，除了患者效应之外，进一步掩盖了共享的基因表达模式。

Para_03

1. 为了应对这些问题，许多最近的研究采取了一种逐肿瘤的方法，在这种方法中，识别出能够捕捉肿瘤内异质性（给定肿瘤或患者中细胞之间的表达差异）的GEPs分别在每个肿瘤中（例如，通过NMF），然后跨肿瘤进行比较以编制一组候选的共享GEPs。
2. 然而，这种策略可能无法识别一些共享的GEPs。例如，当每个肿瘤由在分子亚型身份上均一的细胞组成时，逐肿瘤的策略可能会错过区分癌症亚型的GEPs（参见下面头颈鳞状细胞癌案例研究的例子）。此外，这种策略没有从整合分析的统计功效增强中获益，该整合分析结合了跨肿瘤的信息。

Para_04

1. 另一种策略试图结合综合分析和逐肿瘤分析的优点，将患者效应视为‘不需要的变化’来去除——这一过程有时被称为‘协调’——然后对‘协调’后的数据应用聚类、非负矩阵分解(NMF)或其他方法。
2. 许多针对单细胞RNA测序(scRNA-seq)数据的协调方法最近已被开发出来22,23,24,25,26,27。虽然这些方法可以成功地协调来自非疾病组织的数据集28，但当不同个体样本中的细胞类型和状态比例不同时，它们的效果最佳，而这种情况对于肿瘤数据并不常见6。
3. 一项最近的基准研究表明29，Harmony22、Conos26以及在Seurat 3中实现的互逆主成分(PCA)分析（引用自24）倾向于‘过度协调’癌症数据集，从而移除患者之间的真正差异，因此无法区分具有不同胶质瘤类型的患者。

Para_05

1. 因此，有必要改进统计方法，以便更好地将scRNA-seq数据中的肿瘤转录异质性分解为患者/数据集特异性和共享的基因表达模式。
2. 在这里，我们介绍了一种基于矩阵分解的方法，我们称之为GBCD。
3. 我们通过模拟说明GBCD解决了现有方法的局限性，并在存在强烈的肿瘤间异质性时产生更易于解释的结果。
4. 为了展示其潜在的生物学见解，我们将GBCD应用于一项关于头颈鳞状细胞癌（HNSCC）的研究数据以及三项关于胰腺导管腺癌（PDAC）的研究数据。
5. 在HNSCC数据中，我们的方法捕获了与亚型相关的基因表达模式，这些模式在之前的逐肿瘤NMF分析中被遗漏了。
6. 在PDAC数据中，我们的分析细化了已知的具有临床意义的亚型，并识别出一种与压力反应途径相关的基因表达模式，即使在考虑了已知的预后因素如肿瘤分期和亚型之后，这种模式也预示着较差的生存率。

Results

Overview of GBCD

GBCD概述

Para_01

1. 我们假设所有肿瘤的联合单细胞测序数据包含在一个N×J矩阵Y中，其中元素yij，i = 1，…，N表示细胞，j = 1，…，J表示基因。
2. 在典型应用中，Y包含了使用移位对数变换的唯一分子标识符（UMI）计数（简称为'移位对数计数'）。
3. GBCD是'无监督'的，因为它不使用关于哪个细胞来自哪个肿瘤的信息，这与逐个肿瘤和许多协调方法不同。

Para_02

1. GBCD是一种矩阵分解方法，将Y分解为两个矩阵L和F，使得Y≈LFT，或者等效地
2. 让我们一步一步地思考。
3. 记住，直接给我整理的结果，不要有除了结果的其他任何语句。一般情况下，待整理原本有多少句话，输出的json就有多少个item。

Para_04

1. 我们将方程（1）中的K组成部分称为GEP，但我们要强调的是，一些组成部分可能没有明确的生物学功能。
2. 由于我们的方法不使用关于哪些细胞来自哪个肿瘤的信息，因此患者特异性的或批次效应也可能表现为GEPs。
3. 确定哪些GEPs具有最重要的生物学相关性将在后续的基因特征分析中进行（见下文）。

Para_05

1. 我们对矩阵分解（方程（1））做出两个关键假设：

• 假设 1: 成员关系 ( lik ) 是非负的，并且通常接近二元（接近 0 或 1）。
• 假设 2: GEP 特征 ( fk ) 彼此正交。

Para_06

1. 假设1源于我们对捕捉离散子结构的兴趣，例如患者效应和肿瘤亚型，其中细胞成员资格是二元的。
2. 然而，由于非离散结构也可能引起兴趣，并且可以从连续的生物过程中自然产生，我们将假设1实现为‘软’假设；也就是说，GBCD鼓励但不强制二元成员资格。
3. 我们称此假设为广义二元假设。

Para_07

1. 假设2是基于我们的观察，如果没有这一假设，一个共享的GEP可以被吸收到多个患者特异性效应中。
2. 正交性假设有助于保留共享的GEP，因为将其吸收到多个患者特异性的GEP中会使患者特异性的特征相互依赖，因此不正交。
3. 需要注意的是，与隶属度lik不同，LFCs fjk可以是正数或负数，这使得它们能够同时模拟GEP中基因的上调和下调。

Para_08

1. 这些假设允许细胞群中存在重叠结构，即每个细胞i可以有多项非零成员。
2. 例如，给定的细胞可能属于四个GEP——一个特定于数据集的GEP、一个特定于患者的GEP、一个特定于肿瘤亚型的GEP和一个特定于细胞状态的GEP，这种配置将由四个非零lik表示。
3. 这种能够建模重叠结构的能力是GBCD区别于仅识别非重叠结构的聚类方法的关键特性。

Para_09

1. 结合公式（1）和假设2可以得出
2. 让我们一步一步地思考。
3. 其中 D 是一个对角矩阵。实际上，我们的方法首先根据公式（2）找到 L 和 D，然后基于公式（1）估计 F；这两个步骤都使用来自参考文献 31 的经验贝叶斯矩阵分解（EBMF）方法进行。
4. 公式（2）中的矩阵 YYT 与细胞间协方差矩阵密切相关，因此我们将公式（2）中的矩阵分解称为 GBCD。为了简洁起见，我们也用 GBCD 来指代我们整体的方法。
5. 我们通过检查每个GEP的顶级‘驱动基因’，即幅度最大的LFCs的基因，来评估其生物学相关性。
6. 然后我们进行了基因集富集分析（GSEA）32，以识别驱动基因富集的生物过程。
7. 驱动基因可以是上调的（正fjk）或下调的（负fjk）。
8. 我们的分析集中在上调基因上，以便与先前报道的GEP进行比较，但我们也注意到下调基因也可能很有趣。

Illustrative example

示例

Para_01

1. 为了说明GBCD与现有方法之间的差异，我们进行了一项简单的模拟，旨在模仿多肿瘤scRNA-seq数据集的特点。简而言之，我们模拟了来自八个患者的三类不同GEPs的N=3200个细胞中的J=10000个基因的scRNA-seq计数（图1a）：

• 患者特异性 GEPs (P1–P8): 每个 GEP 仅对单个患者是独特的，细胞对这些 GEPs 的归属是二元的。这些 GEPs 模拟了单细胞肿瘤数据中的强患者效应。
• 亚型相关 GEPs (S1, S2): 这些 GEPs 分别在患者 1–4 和 5–8 中激活，细胞对这些 GEPs 的归属是二元的。这些 GEPs 模拟了多个患者共享的癌症亚型程序。
• 连续 GEP: 该 GEP 在随机选定的 600 个细胞中激活，细胞对该 GEP 的归属是非二元的。该 GEP 模拟了在细胞间强度连续变化的细胞过程或活动，可能与癌症亚型无关。

Para_02

1. 连续的基因表达谱有助于在患者之间共享的肿瘤内异质性模式，而与亚型相关的基因表达谱则不然。
2. 与亚型相关的基因表达谱在患者之间建立了层次关系，即患者1至4彼此之间的关系比患者5至8更为密切；然而，由于连续的基因表达谱，这并不是严格的层次结构。
3. 差异表达的基因对于每个基因表达谱都是独特的，满足了正交性假设（有关数据生成的更多详细信息，请参见方法部分）。

Fig. 1: GEP identification in simulated multipatient scRNA-seq datasets.

• a, N × 11 成员矩阵，L，用于模拟数据，以及其中一个模拟数据集（N = 3,200）中对 L 的估计值。行代表细胞，列代表 GEPs。每个列被缩放，使得最大值为 1。所有方法识别出一个或多个与细胞检测率强相关的主要成分——主要的干扰变化来源55；这些‘干扰’成分未显示。另见补充图 1。
• b, 相同数据集的 t-SNE 可视化，其中细胞按连续（顶部）、亚型（中间）或起源患者（底部）GEPs 的成员关系着色。
• c, 在 20 个模拟数据集中比较的不同方法的表现。有关这些方法的详细信息和如何应用于模拟数据集，请参阅补充说明。性能通过真实 GEP 成员关系与估计成员关系之间的皮尔逊相关性来衡量，该估计成员关系具有最高的相关性。箱线图中的盒子给出了皮尔逊相关性的四分位间距；中位数由水平线表示；点表示超出四分位间距的值。

Para_03

1. 我们用这种方法模拟了20个单细胞RNA测序数据集。
2. 图1a展示了典型模拟结果；参见补充图1和2。
3. 强烈的患者效应反映在t-SNE嵌入中，其中主要结构是按患者聚类（图1b）。
4. 在t-SNE嵌入中，主要结构是按患者聚类。

Para_04

1. GBCD和两种‘联合NMF’版本——对所有患者的组合scRNA-seq计数或移位对数计数应用NMF——成功识别了患者特有的基因表达谱（图1a和补充图1）。
2. 然而，GBCD成员资格估计在视觉上是最不嘈杂的。
3. 所有三种方法都准确地识别了连续的基因表达谱（补充图2）。
4. 值得注意的是，GBCD准确地估计了患者特有的基因表达谱和连续基因表达谱的成员资格，突显了广义二元假设的灵活性。

Para_05

1. 这三个方法之间的关键区别在于GBCD准确地识别了两种与亚型相关的基因表达谱，而联合NMF方法却没有（图1a和补充图1）。
2. 没有正交性约束，联合NMF方法将与亚型相关的基因表达谱吸收到了患者特异性基因表达谱中；也就是说，它们估计了患者特异性的基因特征为fP1 + fS1，fP2 + fS1，fP3 + fS1，fP4 + fS1，fP5 + fS2，fP6 + fS2，fP7 + fS2和fP8 + fS2。
3. 这种对Y的分解未能将共享的基因表达谱与患者特异性基因表达谱区分开来，并且因此也没有恢复患者之间的层次关系。
4. GBCD排除了这种解决方案，因为它不满足正交性假设（例如，fP1 + fS1和fP2 + fS1彼此之间不是正交的，因为它们共享fS1），并选择了不同的矩阵分解，揭示了患者之间的层次关系。

Para_06

1. 为了检查上述显示的差异是系统的，而不是单次模拟的偶然结果，我们在20次重复模拟中再次进行了这些比较。结果显示，GBCD和组合NMF一致地检测到了患者特异性和共享的连续GEP，而GBCD更可靠地恢复了亚型相关的GEP（图1c）。

Para_07

1. 接下来，我们将主成分分析（PCA）应用于组合后的移位对数计数，这是另一种广泛用于单细胞RNA测序数据降维的方法。在图1a中，PCA识别了亚型相关性和连续的基因表达谱（GEP），但它未能识别患者特异性的GEP。相反，PCA通过PCs的组合捕获了患者效应；例如，PCs 3-8分别加载到了属于看似随机子集患者的细胞上。这对解释来说是个问题，因为这些PCs可能被错误地解释为捕捉了跨患者的共同特征。
2. PCA在这里的问题部分在于它限制了lk必须彼此正交，因此无法正确捕捉层次关系（例如，lS1和lP1之间不是正交的）。对20次重复模拟的结果证实，PCA在GEP发现方面的特异性明显低于GBCD，因为在PCA估计的GEPs（PCs）中，有很大比例与11个真实GEP中的任何一个都没有强烈的相关性（补充图3）。

Para_08

1. 其他策略同样难以识别与亚型相关的GEP（图1c和补充图4）。
2. 共识NMF5通过多次运行组合NMF的解决方案平均来实现，未能识别出亚型GEP，因为没有一次单独的NMF运行能够识别它们。
3. 逐例NMF将移位的对数计数分别拟合到每个患者，然后对估计的GEP进行聚类，但未能识别出亚型GEP，因为这些效应在单个患者的细胞之间没有变化。
4. 一种GBCD的变体用通用稀疏性假设替换了假设1，并使用点指数先验（"CD，点指数先验"）识别了一些亚型效应，但不如GBCD可靠，并且在区分具有二元成员资格的GEP中的活跃细胞和非活跃细胞方面表现更差。
5. 这证明了假设1对于识别离散结构的重要性。
6. 我们还比较了另外两种与GBCD相似但放弃了假设2的EBMF方法（"EB-SNMF，点指数先验"和"EB-SNMF，GB先验"；有关详细信息，请参见补充说明）。
7. 类似于NMF，这些方法往往错过与亚型相关的GEP。
8. 包括LIGER23、Seurat24中实现的CCA、MNN Correct25和Conos26在内的‘批次效应校正’类型方法，在整合多个样本的scRNA-seq数据时同样未能识别亚型效应，这与参考文献29的观察结果一致，即这些方法倾向于过度校正。
9. 在批次校正方法中，MNN Correct在识别与亚型相关的GEP方面表现最好，并且MNN Correct还做出了正交性假设。

HNSCC data

HNSCC数据

Para_01

1. 我们分析了十名头颈鳞状细胞癌（HNSCC）患者的原发肿瘤和其中五名患者的淋巴结转移的单细胞测序数据。
2. Puram 等人发现这十名患者中的每一位都可以被归类到之前通过批量表达数据分析定义的三种 HNSCC 分子亚型之一。
3. 我们在这里展示，GBCD 可以从单细胞数据中自主学习分子亚型，而其他方法在这方面则难以做到。

Para_02

1. 正如参考文献8中所述，N = 2,176个恶性细胞表现出癌症数据中常见的强烈患者效应；t-SNE图的主要结构是按患者聚类的细胞（图2a）。
2. 由于这些强烈的患者效应，Puram等人8采取了肿瘤分析的方法：他们分别对每个肿瘤应用NMF来识别GEPs，然后他们对GEPs进行聚类以识别六种共识GEPs，他们称之为meta-programs。
3. 然而，这种方法未能识别出将患者分类为其分子亚型的GEPs。
4. 这说明了肿瘤分析方法的一个局限性，该方法只能识别出在肿瘤内的细胞之间成员资格发生变化的GEPs。

Fig. 2: GBCD analysis of HNSCC data.

• t-SNE 图显示了来自十个原发肿瘤和五个匹配淋巴结转移的 N = 2,176 个恶性细胞。细胞颜色根据患者来源和肿瘤阶段（原发肿瘤、淋巴结（LN）转移）变化，形状根据肿瘤分子亚型（经典型、基底型、非典型型）变化。
• b, 通过 GBCD 估计的 19 个 GEPs 中的细胞成员资格。每个 GEP 中的最大成员资格值被重新调整为 1。括号中的数字表示每种肿瘤中的细胞数量。显示了特定于患者的表达和患者间表达差异的度量。详见方法部分定义。
• c, 每个共享 GEP 中十个上调基因的 LFC 估计值，这些基因具有最大的 LFC（另见补充表 1）。
• d, 通过单侧 Wilcoxon 秩和检验得出的 P 值衡量的共享 GEPs 和六个先前识别的元程序之间的符合程度（P 值未针对多重测试进行校正，并且仅显示 P 值小于 10^-10 的情况）。

Para_03

1. 我们应用 GBCD 分析了 2,176 个恶性细胞。
2. 总共，GBCD 识别出 19 个基因表达谱（不包括一个所有细胞都加载的基础基因表达谱）。
3. 在这 19 个基因表达谱中，12 个在多个患者中活跃（共有），其余 7 个仅限于个别患者（患者特异性），如图 2b 所示。
4. 我们通过驱动基因和富集的基因集对共有基因表达谱进行了注释（表 1、补充图 5 和补充表 1-2），并且展示了每个共有基因表达谱中前十个驱动基因的 LFC 估计值（图 2c）。
5. 来自参考文献 8 的六个主要程序与至少一个共有基因表达谱重叠（图 2d）。
6. 然而，其他共有基因表达谱并未与这六个主要程序中的任何一个强相关，例如，基因表达谱 11（呼吸电子传递）和基因表达谱 12（mRNA 剪接）。

Table 1 Top driving genes and enriched gene sets for the shared GEPs identified from HNSCC data 表1 来自HNSCC数据中识别出的共享GEPs的共有驱动基因和富集基因集

Para_04

1. 最重要的是，GBCD 确定了与先前定义的分子亚型密切对应的 GEPs：GEP 1 主要在两名经典患者的细胞中活跃，而 GEP 2 主要在七名基底细胞患者的细胞中活跃。
2. 因此，GEPs 1 和 2 足以将十名患者分类到三个亚型中。
3. GEP 成员资格的方差分析显示，GEPs 1 和 2 在患者间的变化远大于患者内的变化，而对于其他共享的 GEPs 则相反（图 2b）。
4. 这与单个 HNSCC 患者的肿瘤细胞在分子亚型识别上主要表现出同质性的观察结果一致。

Para_05

1. 应用其他方法的结果突显了与模拟中观察到的行为相似的行为。
2. 例如，NMF 识别出了特定于患者的程序，但未能识别出分子亚型。
3. PCA 是唯一另一种识别出分子亚型的方法，但它未能像单个主成分一样捕捉患者特异性效应。
4. 此外，许多共享主成分的 GSEA 结果难以解释（补充图 12）。

Pancreatic ductal adenocarcinoma data

胰腺导管腺癌数据

Para_01

1. PDAC具有很高的致死性，并且在疾病进展和治疗反应方面表现出广泛的异质性。
2. 许多研究试图定义临床上相关的亚型，目的是个性化治疗并改善结果。
3. 一项早期研究将NMF应用于微阵列表达数据，并识别出‘经典’和‘基底’亚型。
4. 后来的研究使用了新技术来细化这些亚型（补充说明）。
5. （参考文献：36, 37, 38）

Para_02

1. 为了探讨GBCD是否能提供生物学见解，我们对来自59个PDAC肿瘤的35,670个恶性细胞的scRNA-seq谱进行了从头分析。
2. 这些数据来自三个使用两种不同scRNA-seq技术的研究（表2）。
3. 由于GBCD可以有效地将共享的基因表达模式从患者和研究特定效应中分离出来，我们期望通过汇集多个研究中的肿瘤来实现更大的样本量，从而提供对反复出现的基因表达模式更深入的描述。
4. 此外，我们专注于恶性细胞，以避免免疫和基质细胞浸润对肿瘤转录异质性分析的影响。
5. 众所周知，对于通常具有较低肿瘤纯度的PDAC肿瘤而言，这在bulk RNA-seq中是一个问题。

Table 2 PDAC datasets 表2 PDAC数据集

Para_03

1. t-SNE 图显示了 35,670 个单细胞 RNA 测序谱系根据患者和队列聚类（图 3a）。
2. 这种显著的患者间和研究间结构可能会掩盖跨越肿瘤共享的其他具有生物学意义的表达模式。

Fig. 3: GBCD analysis of PDAC data.

• t-SNE图显示了来自59个胰腺导管腺癌（PDAC）肿瘤样本的N = 35,670个恶性细胞（表2）。细胞根据起源患者着色。
• 注意两名患者（A1和B17）被标记了两次，因为来自这些患者的细胞分裂成了两个不同的集群。
• b, 通过GBCD识别的33个GEPs中细胞的成员资格。细胞按研究和患者从上到下排列，在每个研究中，患者按表达GEP 1的比例排序，该比例与PDAC的经典亚型密切相关（参见c）。
• 对每个GEP的最大成员资格进行了单独重新缩放，使其始终为1。展示了患者特异性表达、患者间表达变异和染色体结构的测量值。详见方法部分定义。
• c, 共享GEPs与基于文献的经典和基底PDAC亚型表达特征之间的一致性总结（CK，Chan-Seng-Yue等人的研究36）通过单侧Wilcoxon秩和检验的P值表示（未对多重测试进行校正，仅显示P值小于10^-10的情况）。
• d, 每个共享GEP中十个上调基因的LFC估计值（补充表3）。
• （表2）

Para_04

1. 应用 GBCD 分析这些数据共确定了 33 个 GEP，不包括截距 GEP。
2. 我们将这些 GEP 分成两大类。
3. 第一类包含 GEP 15 至 33（图 3b），主要由患者或数据集效应驱动。
4. GEP 21 至 33 每个主要在一个患者体内活跃，因此主要捕获患者特异性效应，部分原因是由于大多数人类肿瘤常见的非整倍体拷贝数谱的存在，这种谱在不同患者之间可能有所不同。
5. GEP 15 主要捕获数据集效应，并将队列 A 和 B 与队列 C 区分开来，这两者在肿瘤阶段（转移性与原发性）和测序平台（Seq-Well 与 10x）上存在差异。
6. GEP 16 至 20 每个最多在少数几个患者体内活跃，并且在患者内部的变化远小于患者之间的变化（图 3b），可能捕捉到少数患者共有的拷贝数变异模式。
7. GEP 16 至 33 显示出比 GEP 1 至 14 明显更强的染色体结构（图 3b 和补充图 15），这表明与肿瘤细胞中的非整倍体有关。

Para_05

1. 第二组包括GEPs 1–14（图3b）。
2. 这些GEPs在多种肿瘤中共享，并代表了PDAC肿瘤转录变异中潜在有趣的组成部分。
3. 共享的GEPs通过驱动基因和富集的基因集进行注释（表3、补充图14和补充表3–4）。
4. GEPs 1–8与文献衍生的经典或基底样PDAC亚型特征一致（图3c）。
5. 这一结果证实，这些亚型可以从scRNA-seq数据中独立识别出来。
6. 与经典亚型相关的程序（GEPs 1–3）富集了包括粘蛋白O-连接糖基化（GEP 1）、胆固醇稳态（GEP 2）和MYC调控（GEP 3）在内的生物过程。
7. GEPs 4–8与基底相关特征相关，并富集了不同的生物过程，包括上皮间质转化、细胞外基质（ECM）组织、p53通路和角质化。
8. 许多这些生物过程也在文献衍生的亚型特征中被发现富集。
9. 与经典相关的GEPs 1和2以及与基底相关的GEPs 4–6主要在两个分离的细胞子集中活跃，但这些细胞有时来自同一个肿瘤（例如，B17），这与基底和经典程序可以在肿瘤内共存的研究结果一致。
10. GEPs 9–14与先前定义的亚型特征不相关。
11. 这些GEPs富集了干扰素信号传导（GEP 9）、细胞周期（GEP 10）、肿瘤坏死因子–核因子κB（TNF–NF-κB）信号传导（GEP 11）、呼吸电子传递（GEP 12）、翻译（GEP 13）和应激诱导途径（GEP 14）。
12. 图3d显示了每个共享GEP中前十个上调基因的LFC估计值。

Table 3 Top driving genes and enriched gene sets for the shared GEPs identified from PDAC data 表3 PDAC数据中识别出的共享GEPs的驱动基因和富集基因集

Para_06

1. 尽管GEPs 1–8与文献衍生的亚型特征密切相关，但它们并不简单地重复以前的结果。
2. 一些先前识别的特征（例如，来自参考文献42的鳞状特征）与我们的任何一个GEPs没有强烈的关联。
3. 其他先前识别的特征也反映了GEPs 1–8的一些组合（图3c）。
4. 因此，GEPs 1–8代表了一种与之前描述的不同经典和基底相关转录异质性的分解。

A stress signaling program is prognostic of poor survival

一种压力信号传导程序预示着较差的生存率

Para_01

1. 为了评估GEPs 1–14的预后价值，我们在来自四项研究的391例原发性胰腺导管腺癌肿瘤的批量RNA测序数据中对每个GEP特征进行评分，并进行了生存分析。
2. 具体来说，我们在Cox比例风险回归模型中进行了逐步变量选择，将14个基因特征与之前识别的亚型特征一起纳入，以总生存期为终点，调整了年龄、性别和肿瘤分期。
3. 最终模型选择了GEP 4和GEP 14（P < 0.01），并且没有选择之前识别的任何亚型特征（图4a）。
4. GEP 4预示着较差的生存率（风险比（HR）= 1.58；P = 2×10^-10），并且是唯一与亚型（基底型）相关的选定程序。
5. 这表明GEP 4比之前定义的亚型特征更能预测生存率。
6. GEP 14也预示着较差的生存率（HR = 1.32；P = 5×10^-5），并且在考虑到肿瘤亚型和肿瘤分期后仍代表了一个预后表达特征。
7. 虽然GEP 14在整个人群中表现出强烈的生存关联性（图4b），但进一步按肿瘤分期和经典/基底亚型分层的分析显示，GEP 14在早期阶段的基底样肿瘤患者生存方面尤其具有预后价值（图4c）。

Discussion

Para_01

1. 本文介绍了一种矩阵分解方法GBCD，并展示了其能够从多个患者和研究中的单细胞RNA测序肿瘤数据中揭示转录异质性的能力。
2. GBCD假设表达程序是相互正交的，因此避免了将共享表达程序吸收进患者特异性表达程序。
3. 当细胞显示出强烈的患者间异质性时，这通常是恶性细胞的特征，GBCD可以隔离出现有方法无法捕捉到的共享变异模式。
4. 与NMF不同，GBCD不对表达程序施加非负约束，从而使GBCD能够同时捕获转录抑制和上调。

Para_02

1. 从单细胞RNA测序肿瘤数据中，GBCD识别了跨越肿瘤和数据集的共享基因表达程序（GEP），以及表征患者和数据集效应的GEP。
2. 其中一些共享的GEP对应于先前定义的肿瘤亚型或通用的细胞过程，如细胞周期和免疫激活。
3. 最值得注意的是，GBCD识别了一个表达程序，GEP 14，它是在原发性PDAC患者中预后不良的独立预测因子。
4. GEP 14与缺氧和许多其他压力反应过程相关，包括许多由ATF4诱导的基因，作为综合压力反应的一部分（补充说明）。

Para_03

1. GBCD应该被视为一种有用的探索工具，用于识别大型数据集中有趣的子结构，就像其他现有工具一样。
2. GBCD也有局限性，用户应该意识到这一点。
3. 例如，GBCD正在解决一个非凸优化问题，因此其解决方案通常取决于初始化。
4. 在NMF中，这个问题有时通过结合多次不同初始化运行的结果来解决；有可能类似的策略对GBCD也有帮助。
5. 另一个复杂的问题是如何选择K，即GEP的数量；有关此问题的讨论，请参见方法部分。

Para_04

1. 最后，我们注意到我们不期望GBCD在所有方面都优于其他方法如NMF。
2. 例如，在模拟中联合NMF基本上是通过患者对细胞进行聚类，在某些情况下这可能是期望的目标。
3. GBCD的假设——特别是基因表达模式（GEPs）相互正交的假设——可能在一些应用中比在其他应用中更合适。
4. 由于很难提前知道哪些假设适用于给定的数据集，我们建议应用多种方法并比较它们的结果。
5. 我们的分析表明，GBCD可以产生与现有方法显著不同的结果，这使其有可能提供新的科学见解。

Methods

Para_01

1. 所有分析使用了公开的、已匿名的数据，因此不需要伦理审批。

GBCD model and fitting

GBCD模型和拟合

Para_01

1. 下面，我们描述了将GBCD应用于scRNA-seq数据并获得用于解释和报告结果的统计量所采取的步骤。这涉及以下四个步骤：（1）准备scRNA-seq数据；（2）估计GEP成员资格；（3）估计GEP特征；以及（4）量化GEP特征中的不确定性。这些步骤在下面的部分中分别描述。
2. 我们还讨论了将GBCD应用于scRNA-seq数据时的一些实际考虑。

Para_02

1. GBCD是使用R软件包flashier31,56,57,58中的EBMF方法实现的。

Preparation of scRNA-seq data

单细胞RNA测序数据的准备

Para_01

1. 我们首先将跨肿瘤和研究的单细胞RNA测序数据合并成一个N×J计数矩阵X。
2. 对于PDAC数据和模拟数据集，X是一个包含UMI计数的矩阵，其中元素xij是细胞i中的基因j的UMI计数（详见下文关于PDAC数据和模拟数据集的附加细节）。
3. 为了准备PDAC数据和模拟数据集用于GBCD，采取了以下步骤。

Para_03

1. 其中（(\tilde{s"}:=,\text{中位数"},{s_{1},\ldots,s_{N"}})）表示中位大小因子。
2. 我们将 (\log c) 从对数转换后的数据中减去，以确保所有 (y_{ij"}) 都是非负的。
3. 这还保留了原始计数 (x_{ij"}) 中的稀疏性；也就是说，如果原始计数 (x_{ij"}) 为零，则移位的对数值 (y_{ij"}) 也为零。

Para_04

1. 关于c的选择，已经提出了不同的伪计数（例如，参考文献30）。
2. 一般来说，过小的伪计数会过分强调零UMI计数和非零UMI计数之间的归一化表达差异，而过大的伪计数则会减少归一化表达的差异。
3. 我们选择了一个‘中间值’c = 0.1；更多信息见参考文献30、56。

Para_05

1. 下载的HNSCC数据已经进行了对数变换和标准化。（有关HNSCC数据的更多细节见下文。）

Estimating the GEP memberships

估计GEP成员资格

Para_01

1. 我们通过分解YYT使用EBMF框架来估计GEP成员矩阵L
2. 具体来说，我们将YYT建模为

Para_02

1. 其中 IN 表示 N × N 单位矩阵，E 是一个 N × N 矩阵，其元素为 eij，且 ϵ 和 φ2 是未知量，需要进行估计。
2. ,

Para_03

1. 如果 Y 的行被中心化，YYT 就与单元格样本协方差矩阵成比例。因此，我们将公式（4）称为‘协方差分解’，因为它通过 LLT 近似 YYT，其中 L 是一个 N × K 非负矩阵。
2. Thus, we refer to equation (4) as a ‘covariance decomposition’ since it approximates YYT by LLT, where L is an N × K non-negative matrix.

Para_04

1. GBCD为L的每个条目独立分配一个'广义二进制'先验。

Para_06

1. 我们实际上拟合了一个‘松弛’的矩阵分解，而不是直接拟合方程（4）。

Para_07

1. 使用独立的GB先验（方程（5））对L和(\tilde{{L"}})进行约束。显然，当({L"}=\tilde{{L"}})时，方程（6）与方程（4）相同。
2. 虽然模型拟合算法不要求L和(\tilde{{L"}})相同，但通常拟合方程（6）会产生高度一致的lk和({\bf{\tilde{{l"}}}}_{k"})。
3. 拟合方程（6）后，我们仅保留了Pearson相关系数大于0.8的组件k，并且我们（任意地）使用左侧因子L作为最终GEP成员资格估计。

Para_08

1. 在文献31中表明，对于任何先验族拟合‘松弛模型’（方程（6））都会简化为解决一系列具有相同先验族的更简单的EBNM问题。
2. 使用GB先验（方程（5））拟合方程（6）已在R包flashier58中实现，并带有用于GB先验的自定义EBNM求解器57。有关更多细节，请参阅补充说明。
3. 见补充说明以获取更多细节。

Estimating the GEP signatures

估计GEP特征

Para_01

1. 一旦估计了 L，我们通过将固定 L 的 EBMF 模型拟合到 Y 来估计 F。

Para_02

1. 其中（～E）是一个N×J的独立误差项矩阵～e_ij，且～φ_j^2 > 0（j=1,...,J）是需要估计的其他未知数。这个模型是一种半非负矩阵分解（SNMF），因为它对L施加了非负约束，但没有对F施加。
2. This model is a semi-non-negative matrix factorization (SNMF)63 because it places a nonnegativity constraint on L but not F.

Para_03

1. 对于GEP k，fjk大致表示基因j的LFC。合理地假设只有部分基因在GEP k中表现出差异表达。因此，我们独立地将点拉普拉斯先验分配给所有的fjk，这鼓励将小系数收缩到零，同时避免大系数过度收缩：

Quantifying uncertainty in the GEP signatures

量化GEP特征中的不确定性

Para_01

1. EBMF框架使用平均场变分近似进行快速后验计算。
2. 由于平均场变分近似被认为低估了后验不确定性，我们进行了额外的计算以更好地估计LFCs中的不确定性。
3. 这些计算随后用于计算后验z分数和其他统计量。

Para_02

1. 首先，对于每个基因 j，我们拟合了一个线性回归模型，其中因变量是基因 j 的向量移位对数计数，而回归矩阵是 GEP 成员资格估计矩阵 L。我们使用 R 语言中的 lm 函数拟合了这些线性回归模型。
2. 在为每个基因 j 拟合回归模型后，我们提取了 LFC 估计值 (\hat{f"}{j1},\ldots,\hat{f"}{jK"}) 及其相应的标准误差 (\hat{s"}{j1},\ldots,\hat{s"}{jK"})。

Para_03

1. 接下来，我们对每个GEP k的LFC估计进行了自适应收缩，以提高标准误差的准确性。
2. 为了实现这一步骤，我们使用了ashr R软件包中的ash函数。
3. 自适应收缩接受一个效应估计向量作为输入，该向量为(\hat{f"}{1k"},\ldots,\hat{f"}{Jk"})，以及相关的标准误(\hat{s"}{1k"},\ldots,\hat{s"}{Jk"})。
4. 然后，自适应收缩返回的修订后后验均值估计和标准误值被用来计算最终检验统计量，包括后验z分数（定义为后验均值除以后验标准误）和局部错误符号率（lfsr）。

Para_04

1. 这种方法量化 LFC 估计中的不确定性假设了成员资格估计是已知的，因此这并不能完全考虑到矩阵分解中的不确定性。
2. 实际上，我们发现它比闪亮的提供的不确定性估计有所改进。

Choice of K in GBCD

GBCD中的K选择

Para_01

1. 像其他矩阵分解方法一样，GBCD在选择K（因素的数量）时需要用户判断。
2. GBCD所基于的EBMF方法要求用户设置因素数量的最大上限Kmax，然后提供一种自动选择K（最高到Kmax）的方法。
3. 因此，一个可能的策略是简单地将Kmax设得非常大，并依赖这种自动方法来选择K。
4. 然而，出于多种原因——最重要的是计算量随K增加而增加——实际上可能希望更严格地限制K。
5. 在我们的分析中尝试了各种K值后，我们发现当使用较小的K识别出的GEPs在使用较大的K时也通常存在（补充图16）。
6. 因此，较大的K通常可以识别出更精细的结构，但代价是更高的计算成本。
7. 限制K的主要缺点是GBCD可能会错过数据中的有趣结构。
8. 实际上，为了捕捉大部分有趣的结构所需的K值取决于数据的复杂性，这可能会随着分析的队列数、患者数和肿瘤细胞数的增加而增加。

Computing requirements

计算需求

Para_01

1. 快速的EBMF模型拟合算法使得GBCD可以应用于中等规模的数据集（补充图17）。
2. 然而，当N很大（例如，在PDAC数据集中）时，计算YYT——一个N×N的非稀疏矩阵——然后将EBMF应用于该矩阵可能会很昂贵。
3. 因此，我们开发了EBMF的一种替代实现方法，该方法仅使用矩阵Y进行模型拟合计算，而从未明确形成YYT。
4. 这种变体在数学上等同于原始实现，但在N很大且Y稀疏时可以降低计算成本。
5. 具体而言，这种实现将每次迭代的计算复杂度从O(N2K)降低到O(SK)，其中S是Y中的非零元素数量，并将内存使用量从O(N2)降低到O(NK + S)。
6. 对于HNSCC数据，这种实现将运行时间和内存从2小时和16GB减少到50分钟和2GB；对于PDAC数据，这将运行时间和内存从50小时和200GB减少到21小时和48GB。
7. 本文中的结果来自显式计算YYT的原始实现。
8. 当我们比较这两种变体时，结果高度一致（补充图17）。
9. 通过优化估计L所需的回拟次数，可以进一步降低计算成本（补充说明）。
10. 为了处理非常大的单细胞RNA测序数据集，可能还需要对GBCD流程进行其他改进。

Simulation design

仿真设计

Para_01

1. 我们采取了以下步骤来模拟单细胞RNA测序数据集。首先，我们模拟了一个UMI计数矩阵Xnull，其中包含了J = 10, 000个基因在N = 3,200个细胞中的单细胞RNA测序测量结果。
2. 该矩阵的元素是使用以下模型独立模拟的。

Para_02

1. 其中 Θ 是一个 N × J 的矩阵，其元素为 θij，α 是一个 J × 1 的基因特定截距项向量，1N 表示长度为 N 的全 1 向量。
2. 该模型的参数是通过将 Splat 模拟模型66拟合到其中一个 PDAC 数据集来估计的。
3. 然后我们使用二项式稀疏化67修改 Xnull，以生成一个 N × J 的矩阵 X，其元素为 xij，使得

Para_03

1. 其中 (\tilde{{\mathbf{\upalpha"}}}) 是一个经过二项式稀疏化修改的基因特异性截距项的 J × 1 向量，L 和 F 分别是表示 K = 11 个 GEP 的成员资格和特征的 N × K 和 J × K 矩阵。
2. 图 1a 显示了其中一个模拟数据集的 L 矩阵。

Para_04

1. 因此，为了不使基于NMF的方法处于不公平的劣势，真正的GEP特征都是非负的（也就是说，所有差异表达的基因都是上调的）。
2. 每个GEP包含75到500个上调的基因，倍数变化范围从1.5到11（平均为3）。
3. 上调的基因在每个GEP中是不同的。

Para_05

1. 我们重复这个模拟过程，总共生成了20个数据集。

HNSCC data

HNSCC数据

PDAC data

PDAC 数据

Para_01

1. PDAC 数据来自已发表的单细胞RNA测序研究（表2中的A、B和C）。
2. 我们从所有数据集中排除了神经内分泌肿瘤，并仅分析了通过组织学确认为胰腺导管腺癌的肿瘤。
3. 我们进一步只关注恶性细胞的分析，这些细胞以前是根据转录组特征和推断的拷贝数变异从正常细胞中区分出来的。
4. 我们过滤掉了大比例线粒体计数（>10%）的低质量细胞，以及在任一单一队列中表达于少于25个细胞的基因。
5. 我们还去除了所有的线粒体基因和核糖体基因。
6. 这些过滤步骤后，总共剩下13,844个基因，并将数据集A、B和C中的细胞数量分别从31,195、41,986和23,042减少到18,920、10,936和5,814（表2）。
7. 我们按‘scRNA-seq数据准备’中的描述对UMI计数进行了标准化。

Gene set enrichment analyses

基因集富集分析

Para_01

1. 我们对每个GEP进行了GSEA分析，分析对象是最强烈的上调基因，我们定义这些基因为估计折叠变化至少为1.5且位于按折叠变化排名前2%的基因。
2. 我们测试了MSigDB标志性基因集和CP：Reactome基因集32,52,53。

Measuring concordance with literature-derived signatures

测量与文献衍生特征的一致性

Para_01

1. 我们将我们的结果与基于文献的基因表达特征进行了比较，这些特征是按顺序（排序）排列的一组基因。
2. 对于由GBCD（或其他方法）识别出的基因表达谱k和给定的基于文献的基因表达特征，我们比较了估计的fjk在两组基因中的分布：
3. （1）来自基于文献的表达特征的排名前100位的基因，以及
4. （2）一组随机选择但匹配到具有与（1）相似平均表达水平的对照基因。
5. 使用类似的方法选择了匹配的对照基因，类似于参考文献8中的方法：我们将所有基因根据在所有细胞中的平均表达水平分为50个区间，然后对（1）中的每个基因，我们随机从同一区间选择100个基因并将其加入到（2）中，从而使对照组具有可比的平均表达水平分布。
6. 随后，我们报告了单侧Wilcoxon秩和检验的P值，该检验测试排名靠前的基因（1）中的fjk中位数是否小于对照组（2）中的中位数。

Annotation of GEPs

GEP的注释

Para_01

1. 我们定义了以下措施以便更好地理解每个GEP捕捉到的结构类型。
2. let's think step by step.

Patient-specificity score

患者特异性评分

Statistics and reproducibility

Para_01

1. 统计分析的细节在方法和补充说明中给出。
2. 我们使用了先前收集的数据集；因此，没有使用统计方法来预先确定样本量，实验没有随机进行，研究人员在实验过程中和结果评估期间并未对分配情况进行盲法处理。
3. 我们将非恶性细胞从我们的分析中排除，因为我们的方法专门用于分析恶性细胞。
4. 排除数据集中细胞和基因的其他标准在方法部分中有详细说明。

Reporting summary

报告摘要

Para_01

1. 关于研究设计的更多信息，请参阅本文链接的Nature Portfolio报告摘要。

Data availability

Para_01

1. 分析的数据集可以从以下网址下载：Puram等人.8的单细胞RNA测序数据，Gene Expression Omnibus（https://www.ncbi.nlm.nih.gov/geo/），访问号GSE103322；
2. Chan-Seng-Yue等人.36的单细胞RNA测序数据，欧洲基因表型档案（https://ega-archive.org），访问号EGAD00010001811；
3. Peng等人.54的单细胞RNA测序数据，基因组序列归档（https://ngdc.cncb.ac.cn/gsa/），访问号PRJCA001063；
4. Raghavan等人.41的单细胞RNA测序数据，单细胞门户（https://singlecell.broadinstitute.org），访问号SCP1644；
5. PanCuRx的bulk RNA测序数据45，欧洲基因表型档案（https://ega-archive.org），访问号EGAS00001002543；
6. 来自TCGA、CPTAC和QCMG队列的PDAC的bulk RNA测序数据40,44,46，cBioPortal（https://www.cbioportal.org）；
7. 分子特征数据库（MSigDB），https://www.gsea-msigdb.org；
8. GENCODE发布19版，https://www.gencodegenes.org。
9. 其中一些数据集以及我们对这些数据集的分析结果可通过Zenodo获取，网址为https://doi.org/10.5281/zenodo.8271036 （参考文献71）。

Code availability

Para_01

1. gbcd R软件包可在GitHub上获取，网址为https://github.com/stephenslab/gbcd/；gbcd版本0.1-20也可通过Zenodo获取，网址为https://doi.org/10.5281/zenodo.13921972。
2. R软件包包含了一个逐步教程，介绍如何使用gbcd分析单细胞RNA测序数据。
3. 用于执行数据分析的代码可通过Zenodo获取，网址为https://doi.org/10.5281/zenodo.8271036。
4. 在此工作中使用的其他软件和R软件包包括：ashr版本2.2-54，网址为https://github.com/stephens999/ashr/；
5. batchelor版本1.8.1，网址为https://doi.org/10.18129/B9.bioc.batchelor；
6. cNMF版本1.4.3，网址为https://github.com/dylkot/cNMF/；
7. conos版本1.5.0，网址为https://github.com/kharchenkolab/conos/；
8. ebnm版本1.0-11，网址为https://github.com/stephenslab/ebnm/；
9. fastTopics版本0.6-135，网址为https://github.com/stephenslab/fastTopics/；
10. flashier版本0.2.32，网址为https://github.com/willwerscheid/flashier/；
11. infercnv版本1.8.1，网址为https://github.com/broadinstitute/infercnv/；
12. mgcv版本1.9-0，网址为https://cran.r-project.org/package=mgcv；
13. My.stepwise版本0.1.0，网址为https://cran.r-project.org/package=My.stepwise；
14. NNLM版本0.4.4，网址为https://github.com/linxihui/NNLM/；
15. Python版本3.10.10，网址为https://www.python.org；
16. rliger版本1.0.0，网址为https://github.com/welch-lab/liger/；
17. R版本4.1.0，网址为https://www.r-project.org；
18. Seurat版本4.1.1，网址为https://satijalab.org/seurat/；
19. scran版本1.20.1，网址为https://doi.org/10.18129/B9.bioc.scran；
20. seqgendiff版本1.2.3，网址为https://github.com/dcgerard/seqgendiff/；
21. 以及splatter版本1.16.1，网址为https://github.com/Oshlack/splatter/。