search

ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
108 stars 30 forks source link

R 包精读 | 使用 Harmony 快速、敏感且准确地整合单细胞数据 #5573

Closed ixxmu closed 1 week ago

ixxmu commented 1 week ago

https://mp.weixin.qq.com/s/jpugbjhVvOQIESTCZcbW-A

ixxmu commented 1 week ago

R 包精读 | 使用 Harmony 快速、敏感且准确地整合单细胞数据 by BioJournal Link

Basic Information

  • 英文标题： Fast, sensitive and accurate integration of single-cell data with Harmony
  • 中文标题：快速、敏感且准确地整合单细胞数据与Harmony
  • 发表日期：18 November 2019
  • 文章类型：Article
  • 所属期刊：Nature Methods
  • 文章作者：Ilya Korsunsky | Soumya Raychaudhuri
  • 文章链接：https://www.nature.com/articles/s41592-019-0619-0

Abstract

  1. 单细胞RNA测序数据集的日益多样化使得能够在广泛的生物学和临床条件下全面表征细胞类型。
  2. 然而，将它们一起分析是具有挑战性的，特别是当数据集使用不同技术进行检测时，因为生物学和技术差异是交织在一起的。
  3. 我们介绍了Harmony（https://github.com/immunogenomics/harmony），一种将细胞投影到共享嵌入中的算法，在这种嵌入中，细胞按细胞类型而不是数据集特定条件进行分组。
  4. Harmony同时考虑了多个实验和生物学因素。
  5. 在六次分析中，我们展示了Harmony相较于先前发表的算法的优越性能，同时需要更少的计算资源。
  6. Harmony能够在个人计算机上整合约10^6个细胞。
  7. 我们将Harmony应用于来自具有较大实验差异的数据集的外周血单核细胞、五项胰腺胰岛细胞研究、小鼠胚胎发生数据集以及scRNA-seq与空间转录组数据的整合。

Main

  1. 最近的技术进步使得在一次实验中无偏地单细胞转录组分析数千个细胞成为可能。
  2. 诸如人类细胞图谱（HCA）和加速药物合作项目等项目，展示了日益增长的原始人类组织参考数据集。
  3. 虽然单个实验逐步扩展了我们对细胞类型的理解，但一个全面的健康和疾病细胞目录将需要整合来自不同捐赠者、研究和技术平台的多数据集的能力。
  4. 此外，在转化研究中，跨组织和临床条件的联合分析对于识别疾病扩展的细胞群体至关重要。
  5. 由于来自不同研究的单细胞RNA测序数据集中的有意义生物学变异通常被数据来源混淆，研究人员开发了无监督的多数据集整合算法。
  6. 这些方法将来自不同实验条件和生物学背景的细胞嵌入到一个共同的低维嵌入中，以实现跨数据集的共享细胞类型识别。
  7. 在这里，我们介绍Harmony，一种用于健壮、可扩展和灵活的多数据集集成的算法，以应对无监督单细胞RNA测序联合嵌入的四个关键挑战：扩展到大数据集，识别广泛的细胞群体和细粒度的亚群体，适应复杂实验设计的灵活性，以及跨模态集成的能力。
  8. 我们将Harmony应用于多种示例，包括细胞系、使用不同技术检测的外周血单核细胞（PBMCs）、来自多个供体和研究的胰腺胰岛细胞的荟萃分析、小鼠胚胎发育的纵向样本，以及单细胞RNA测序数据与空间转录组学数据的跨模态集成。
  9. Harmony作为一个R包在github上可用（https://github.com/immunogenomics/harmony），包含用于独立分析和Seurat7管道分析的功能。

Results

Harmony iteratively learns a cell-specific linear correction function

Harmony 通过迭代学习一个特定细胞的线性校正函数

  1. 和谐算法首先从低维嵌入细胞开始，例如主成分分析（PCA），（补充说明1和方法）。
  2. 利用这种嵌入，和谐算法首先将细胞分组到多数据集簇中（图1a）。
  3. 我们使用软聚类将细胞分配到潜在的多个簇中，以考虑到细胞状态之间的平滑过渡。
  4. 这些簇作为代理变量，而不是用于识别离散的细胞类型。
  5. 我们开发了一种新的软k-means聚类算法，该算法倾向于包含来自多个数据集的细胞的簇（方法）。
  6. 包含来自少数数据子集的细胞的簇会通过信息论指标进行惩罚。
  7. 和谐算法允许多种不同的惩罚，以适应多种技术或生物学因素，例如不同的批次和组织来源。
  8. 软聚类保留了离散和连续的拓扑结构，同时避免了因在多个数据集中过快地最大化表示而可能导致的局部最小值。
  9. 聚类后，每个数据集都有特定的簇质心（图1b），用于计算簇特定的线性校正因子（图1c）。
  10. 由于簇对应于细胞类型和状态，簇特定的校正因子对应于单个细胞类型和细胞状态特定的校正因子。
  11. 通过这种方式，和谐算法学习了一个简单的线性调整函数，该函数对内在的细胞表型敏感。
  12. 最后，每个细胞被分配这些项的簇加权平均值，并通过其特定的线性因子进行校正（图1d）。
  13. 由于每个细胞可能处于多个簇中，每个细胞都有一个潜在的独特校正因子。
  14. 和谐算法迭代这四个步骤，直到细胞簇分配稳定为止。

Fig. 1: Overview of Harmony algorithm.

  • PCA 将细胞嵌入到降维的空间中。
  • Harmony 接受这个降维空间中的细胞坐标，并运行一个迭代算法来调整数据集特定的效应。
  • a, Harmony 使用模糊聚类将每个细胞分配到多个簇中，同时一个惩罚项确保每个簇内数据集的多样性最大化。
  • b, Harmony 为每个簇计算一个全局质心，以及每个簇的数据集特定质心。
  • c, 在每个簇内，Harmony 基于质心为每个数据集计算一个校正因子。
  • d, 最后，Harmony 用一个细胞特定因子校正每个细胞：这是步骤 a 中通过细胞的软簇分配加权的数据集校正因子的线性组合。
  • Harmony 重复步骤 a 到 d，直到收敛。
  • 随着每一轮的进行，簇分配和数据集之间的依赖性逐渐减弱。
  • 数据集用颜色表示，细胞类型用不同的形状表示。

Quantifying performance in cell-line data

量化细胞系数据中的性能

  1. 我们首先使用三个精心控制的数据库来评估Harmony，以评估其在整合（数据集混合）和准确性（细胞类型不混合）方面的性能。
  2. 通过混合所有细胞，无论其细胞身份，可以实现完美的整合。
  3. 同样，通过将细胞划分为广泛的簇，而在小邻域内不混合数据集，可以实现高准确性。
  4. 在这种情况下，定义了广泛的细胞状态，但细粒度的细胞亚状态和亚型因原始数据集而混淆。
  5. 为了量化这种嵌入的整合和准确性，我们定义了一个客观指标：每个细胞局部邻域的局部逆辛普森指数（LISI，见方法）。
  6. 为了评估整合，我们采用‘整合LISI’（iLISI，见图2a），它定义了邻域中有效数据集的数量。
  7. 仅由单个数据集表示的邻域获得iLISI为1，而由两个数据集的等量细胞组成的邻域获得iLISI为2。
  8. 请注意，即使在理想混合下，如果数据集的细胞数量不同，iLISI也会小于2。
  9. 为了评估准确性，我们使用‘细胞类型LISI’（cLISI，见图2b），这是相同的数学度量，但应用于细胞类型而非数据集标签。
  10. 准确的整合应保持cLISI为1，反映嵌入中独特细胞类型的分离。
  11. 错误的嵌入将包括cLISI为2的邻域，表明邻居有两种不同类型的细胞。

Fig. 2: Quantitative assessment of dataset mixing and cell-type accuracy with cell-line datasets.

  • a, iLISI 衡量嵌入中数据集之间的混合程度，从未混合空间的 1 到完全混合空间中的 B，即分析中的数据集数量。
  • b, cLISI 使用相同的公式衡量整合准确性，但计算基于细胞类型标签。
  • 一个准确的嵌入在每个邻域中的 cLISI 接近 1，反映了不同细胞类型的分离。
  • 分析了来自纯细胞系数据集（紫色和黄色）的 n = 3,255 个 Jurkat 细胞和 n = 2,859 个人胚胎肾 293T（HEK293T）细胞，以及来自混合（绿色）细胞系数据集的 n = 1,799 个 Jurkat 细胞和 n = 1,565 个 HEK293T 细胞。
  • c,d, 在 Harmony 整合之前，细胞按数据集（c）和已知细胞类型（d）分组。
  • 为每个细胞的邻域计算了 iLISI（c）和 cLISI（d），并用密度图进行了总结。
  • e,f, 在 Harmony 整合之后，混合数据集中的细胞与其他数据集混合（e），通过将 Jurkat 细胞与 Jurkat 细胞混合，HEK293T 细胞与 HEK293T 细胞混合实现（f）。
  • 在 Harmony 嵌入中重新计算了 iLISI（e）和 cLISI（f）。
  1. 我们从两个细胞系的三个数据集开始：（1）纯Jurkat，（2）纯293T和（3）50/50混合数据集。这些数据集非常适合说明和评估，因为每个细胞都可以明确标记为Jurkat或293T（补充图1a）。
  2. 彻底的整合会将混合数据集中的1,799个Jurkat细胞与纯Jurkat数据集中的3,255个细胞混合，并将混合数据集中的1,565个293T细胞与纯293T数据集中的2,859个细胞混合。
  3. 因此，我们预计平均iLISI值范围从1（反映无整合）到1.8（= 1/[(1,799/(1,799 + 2,859))^2 + (3,255/(1,799 + 3,255))^2]）对于Jurkat细胞，以及1.5（= 1/[(1,565/(1,565 + 2,859))^2 + (2,859/(1,565 + 2,859))^2]）对于293T细胞，反映最大准确整合。
  4. 应用标准PCA流程后进行UMAP嵌入显示，细胞广泛地按数据集和细胞类型分组。这一点既直观可见，又通过高准确度量化（图2c,d），低cLISI（中位iLISI 1.00，95%置信区间（CI）[1.00, 1.00]）反映了这一点。
  5. 然而，iLISI（中位iLISI 1.01，95% CI [1.00, 1.61]）也较低，反映了不完全整合以及每种细胞类型内部丰富的结构，反映了原始数据集。
  6. 经过Harmony处理后，50/50数据集的细胞适当地混合到纯数据集中（图2e）。增加的iLISI（中位iLISI 1.59，95% CI [1.27, 1.97]）反映了数据集的混合，而低cLISI（图2f，中位iLISI 1.00，95% CI [1.00, 1.02]）反映了Jurkat和293T细胞的准确分离。
  7. iLISI和cLISI提供了一种定量比较多种算法整合和准确性的方法。我们使用MNN Correct、BBKNN、MultiCCA和Scanorama重复了整合和LISI分析，并显示它们产生的嵌入在统计上整合效果较差（补充结果，图1b和表1）。
  8. 这个基准测试展示了评估混合和准确性的两个关键指标，并表明Harmony在这两个指标上表现良好，在对细胞系数据集的严格控制分析中。
  9. LISI的一个潜在陷阱是它对大小差异极大的数据集敏感。
  10. 在这种情况下，大多数邻域可能被单个数据集主导，LISI值变得难以解释（补充说明2）。

Harmony scales to large data

Harmony 适用于大规模数据

  1. 我们评估了Harmony的计算性能，测量了总运行时间和最大内存使用量。
  2. 为了展示Harmony与其他方法的可扩展性对比，我们对HCA数据12（来自16位捐赠者和两种组织的528,688个细胞）进行了降采样，创建了包含500,000、250,000、125,000、60,000和30,000个细胞的五个基准数据集。
  3. 我们报告了所有基准数据集的运行时间和内存（补充表2和3）。
  4. Harmony的运行时间在所有数据集中表现良好（图3a），从30,000个细胞的4分钟到500,000个细胞的68分钟，比MultiCCA和MNN Correct快30到200倍。
  5. 对于最多包含125,000个细胞的数据集，Harmony、BBKNN和Scanorama的运行时间相当。
  6. 与其他算法相比，Harmony所需的内存显著减少（图3b）：在30,000个细胞上仅需0.9千兆字节（GB），在500,000个细胞上需7.2 GB。
  7. 在125,000个细胞时，Harmony所需的内存比Scanorama、MNN Correct和Seurat MultiCCA少30到50倍；这些其他方法无法扩展到超过125,000个细胞。
  8. Harmony返回的集成嵌入（图3c）比其他竞争算法多得多（补充结果），允许识别跨组织和捐赠者的共享细胞类型（图3d）。
  9. 这些结果表明，Harmony在计算上高效，能够甚至在个人计算机上分析大型数据集（105到106个细胞）。

Fig. 3: Computational efficiency benchmarks. BBKNN, Scanorama, MNN Correct and MultiCCA are compared on five downsampled HCA datasets of increasing sizes.

  • a,b, 显示了分析每个数据集所需的总运行时间（a）和最大内存（b）。Scanorama、MultiCCA和MNN Correct在250,000和500,000个细胞数据集上因内存请求过多而被终止。
  • c, 这里可视化了500,000个细胞分析中Harmony嵌入的组织混合情况。这包括来自8个独立供体的n = 239,794个脐带血细胞和来自8个独立供体的n = 260,206个骨髓细胞。
  • d, 在Harmony嵌入中，通过典型标记对聚集的细胞群体进行标注：前T细胞、CD4幼稚T细胞、CD4记忆T细胞、T调节细胞、CD8幼稚T细胞、CD8效应T细胞、自然杀伤细胞、前B细胞、幼稚B细胞、记忆B细胞、浆细胞、浆细胞样树突状细胞（pDC）、常规树突状细胞（DC）、粒细胞巨噬细胞祖细胞（GMP）、CD16−单核细胞（CD14单核）、CD16+单核细胞（CD16单核）、同时阳性表达巨核细胞标记的单核细胞群体（PPBP单核）、巨核细胞（Mk）、红细胞祖细胞（Eryth）以及一群造血干细胞和多能祖细胞（HSC/MPP）。源数据

Identification of broad and fine-grained PBMC subpopulations

PBMC亚群的广泛和细粒度识别

  1. 为了评估Harmony在更具挑战性的场景下的表现，我们收集了三组人类PBMCs数据集，每组数据都在Chromium 10X平台上进行检测，但使用不同的制备协议：3′端v1（3pV1）、3′端v2（3pV2）和5′端（5p）化学方法。
  2. 将所有细胞混合后，我们进行了联合分析。
  3. 在整合之前，细胞主要按数据集分组（图4a，中位iLISI 1.00，95% CI [1.00, 1.00]）。

Fig. 4: Fine-grained subpopulation identification in PBMCs across technologies.

  • 三个PBMC数据集使用10X进行了检测，采用了不同的文库构建协议：5′（橙色，n = 7,697个细胞），3′ V1（紫色，n = 4,809个细胞）和3′ V2（绿色，n = 8,380个细胞）。
  • a，在整合之前，细胞按数据集分组。
  • b，经过Harmony整合后，数据集混合在一起。
  • c,d，使用典型标记物（c），我们识别出（d）五种共享的T细胞亚型和两种共享的B细胞亚型。
  1. Harmony 将细胞聚类成生物学上连贯的组（补充图 6a），并消除了每个簇内的数据集特定变异。
  2. 在最终的整合空间中，数据集混合得很好（图 4b，中位 iLISI 1.96，95% CI [1.36, 2.56]），比其他方法更好（中位 iLISI 1.02，95% CI [1.00, 2.11]）。
  3. 我们确认超过 83% 的细胞在 Harmony 中的 iLISI 值显著（假发现率 (FDR) < 5%）高于任何其他算法（补充图 6b,c）。
  4. 为了评估准确性，在每个数据集中，我们分别用主要预期群体的典型标记（方法）对广泛的细胞簇进行注释（补充图 6d）：单核细胞（CD14+ 或 CD16+）、树突状细胞（FCER1A+）、B 细胞（CD20+）、T 细胞（CD3+）、巨核细胞（PPBP+）和自然杀伤（NK）细胞（CD3−/GNLY+）在聚类之前。
  5. 我们观察到 Harmony（中位 cLISI 1.00，95% CI [1.00, 1.15]）保留了细胞类型之间的差异，类似于其他整合方法（中位 cLISI 1.00，95% CI [1.00, 1.28]）。
  6. 与其他算法相比，更大的数据集整合提供了一个独特的机会来识别细粒度的细胞亚型（补充图 6e）。
  7. 使用典型标记（图 4e），我们识别了共享的细胞亚群（图 4f），包括幼稚 CD4 T 细胞（CD4+/CCR7+）、效应记忆 CD4 T 细胞（CD4+/CCR7−）、Treg 细胞（CD4+/FOXP3+）、记忆 CD8 细胞（CD8+/GZMK−）、效应 CD8 T 细胞（CD8+/GZMK+）、幼稚 B 细胞（CD20+/CD27−）和记忆 B 细胞（CD20+/CD27+）。
  8. 在其他算法产生的嵌入中，中位 iLISI 未超过 1.1（补充表 5）。
  9. 因此，上述识别的亚型存在于数据集特定的，而不是数据集混合的簇中（补充图 7）。
  10. 即使我们对整个细胞群体进行下采样，创建具有非重叠细胞类型的不平衡数据集，我们也能保持高质量的结果（补充结果和补充图 8–10）。
  11. 在完整和下采样的 PBMC 数据集中，Harmony 对参数选择具有鲁棒性，尤其是多样性惩罚（补充结果和补充图 11）。
  12. 这些结果表明，Harmony 成功地考虑了不同协议之间的技术变异，并整合了许多不同的细胞类型，同时保留了数据中的大尺度和细粒度结构，即使数据集之间有非重叠的群体。

Simultaneous integration across donors and technologiesidentifies rare pancreas islet subtypes

跨供体和技术的同步整合识别出罕见的胰腺胰岛亚型

  1. 我们考虑了一个更复杂的实验设计，其中必须同时对多个变异源进行整合。
  2. 我们从五个独立的研究中收集了人类胰腺胰岛细胞13,14,15,16,17，每个研究都使用了不同的技术平台。
  3. 跨这些平台的整合尤其具有挑战性，因为每个数据集中的细胞数量和每个细胞中测量的独特基因数量差异很大（补充图12）。
  4. 此外，在两个数据集中13,14，作者注意到了显著的供体特异性效应。
  5. 成功整合这些研究必须考虑到不同技术的影响以及这些研究中使用的36名供体的影响。
  6. 这些效应可能会以不同方式影响细胞类型。
  7. Harmony是目前唯一专门为单细胞数据设计的整合算法，能够明确地整合多个变量。
  8. 与之前一样，我们使用在每个数据集中独立识别的典型细胞类型来评估细胞类型准确性cLISI（补充图13a）：α细胞（GCG+）、β细胞（MAFA+）、γ细胞（PPY+）、δ细胞（SST+）、腺泡细胞（PRSS1+）、导管细胞（KRT19+）、内皮细胞（CDH5+）、星状细胞（COL1A2+）和免疫细胞（PTPRC+）。
  9. 由于有两个整合变量，我们评估了供体iLISI和技术iLISI。
  10. 在整合之前，PCA按技术将细胞分开（图5a，中位iLISI 1.00，95% CI [1.00, 1.06]），按供体分开（图5b，中位iLISI 1.42，95% CI [1.00, 5.50]），以及按细胞类型分开（中位cLISI 1.00，95% CI [1.00, 1.48]）。
  11. 供体iLISI的广泛范围反映了在CEL-seq、CEL-seq2和Fluidigm C1数据集中，许多供体在整合前已经很好地混合了。
  12. Harmony通过技术（图5c，中位iLISI 2.17，95% CI [1.02, 3.91]）和供体（图5d，中位iLISI 5.05，95% CI [1.24, 10.05]）整合细胞（补充图13b），同时根据先前定义的类型正确合并细胞（准确性>98%，补充图13c）。
  13. 我们还使用BBKNN、Scanorama、MNN Correct和MultiCCA对这些数据进行了基准测试，以及使用limma18，后者可以在多个层面上进行校正。
  14. 只有Harmony能够显著混合供体和技术（补充结果和补充图14）。

Fig. 5: Integration of pancreatic islet cells by both donor and technology.

  • 来自36位捐赠者的14,746个人类胰腺胰岛细胞在五种不同的技术平台上进行了检测：inDrop（4位捐赠者，8,569个细胞），Fluidigm C1（13位捐赠者，638个细胞），Smart-Seq2（10位捐赠者，2,355个细胞），CEL-seq（5位捐赠者，946个细胞）和CEL-seq2（4位捐赠者，2,238个细胞）。
  • a,b，细胞最初按（a）技术分组，用不同颜色表示，以及（b）捐赠者，用颜色的深浅表示。
  • c,d，Harmony同时在（c）技术和（d）捐赠者之间整合细胞。
  • e，在Harmony嵌入中的聚类。
  • f，不同细胞类型的比例。
  • g,h，我们将β内质网应激细胞群体（绿色，n = 306个细胞）中的内质网应激基因（g）和β内分泌功能基因（h）的log表达量（标准化为每十万计数的CP10K）与β细胞群体（橙色，n = 3,374个细胞）中的表达量进行了比较。
  • i,j，同样，我们将α内质网应激细胞群体（紫色，392个细胞）中的内质网应激基因（i）和α内分泌功能基因（j）的表达量与α细胞群体（绿色，n = 3,978个细胞）进行了比较。
  • 差异表达分析采用双尾调节t检验进行。FDR通过Benjamini–Hochberg程序计算。
  • k，两种内质网应激群体的丰度在n = 36位捐赠者中进行了相关性分析。相关性通过Spearman’s rho计算，名义P值通过算法AS 89估计。只进行了一次测试，因此未应用多重假设校正。
  1. Harmony 能够识别出五个数据集中罕见的（<2%）细胞亚型（图5e）。
  2. 我们使用典型标记物对先前描述的亚型进行标记：活化星形细胞（PDGFRA+）、静止星形细胞（RGS5+）、肥大细胞（BTK+）、巨噬细胞（C1QC+）以及内质网（ER）应激下的β细胞（图5g）。
  3. β ER 应激细胞可能代表一个功能失调的群体，其特征为 ER 应激基因 DDIT3（FDR = 2.1 × 10−187）、ATF3（FDR = 3.2 × 10−79）、ATF4（FDR = 1.9 × 10−76）和 PPP1R15A（FDR = 1.8 × 10−145）。
  4. 这一簇细胞显著低表达对β细胞身份和功能关键的基因：PDX1（FDR = 1.1 × 10−24）、MAFA（FDR = 4.0 × 10−30）、INSM1（FDR = 3.2 × 10−4）和 NEUROD1（FDR = 9.6 × 10−42）（图5h）。
  5. 此外，Sachdeva 等人提出，PDX1 缺乏会使β细胞功能减弱，并使其暴露于 ER 应激诱导的凋亡中。
  6. 我们还观察到了一个据我们所知尚未被先前描述的α细胞亚群。
  7. 这个簇也富含与内质网应激相关的基因（图5i，DDIT3（FDR = 2.6 × 10−29），ATF3（FDR = 5.5 × 10−55），ATF4（FDR = 2.4 × 10−74）和PPP1R15A（FDR = 6.7 × 10−113））。
  8. 与β内质网应激群体类似，这些α细胞也显著低表达维持正常功能所必需的基因：GCG（FDR = 4.2 × 10−13），ISL1（FDR = 1.0 × 10−22），ARX（FDR = 3.8 × 10−26）和MAFB（FDR = 1.0× 10−71）（图5j）。
  9. 最近的一项研究报道了小鼠α细胞中的内质网应激，并将这种应激与功能失调的胰高血糖素分泌联系起来。
  10. 此外，我们发现α和β内质网应激细胞的比例在所有数据集的捐赠者之间存在显著相关性（Spearman r = 0.46，P = 0.004，图5k）。
  11. 这些结果表明，在人类糖尿病期间，α细胞损伤可能与β细胞功能障碍并行发生。
  12. 这个稀有群体在原始研究中几乎无法识别，因为它要么太稀少，要么被捐赠者特异性变异所掩盖（补充图15）。

Harmony integrates time course developmental trajectories

和声整合了时间进程发展轨迹

  1. 我们接下来评估了Harmony整合具有平滑轨迹数据集的能力，而不是离散的细胞类型。
  2. 我们分析了来自小鼠造血系统八个时间点的15,875个细胞，从E6.75到E8.5以及混合囊胚形成期。
  3. 这个数据集对Harmony提出了新的挑战：每个样本都是在不同的发育时间点采集的。
  4. 因此，没有任何一个样本包含所有细胞类型（补充图16a）。
  5. 细胞最初按样本ID分组，而不是按细胞类型分组（补充图16b）。
  6. 经过Harmony整合后，样本混合良好，主要按先前定义的细胞类型分组（补充图16c）。
  7. Harmony保留了发育细胞的连续结构，而不是错误地将细胞聚类成离散群体。
  8. Harmony的软聚类对于保持平滑性至关重要：它将过渡群体中的细胞分配到更多聚类中，而将更终端群体中的细胞分配到更少聚类中（补充图16d）。
  9. 使用校正后的Harmony嵌入，我们使用DDRTree进行了轨迹分析（补充图16e）。
  10. 我们恢复了一个分支轨迹结构，正确捕捉了从共同中胚层和血液内皮祖群体到分化内皮和红细胞群体的进展。
  11. Harmony保留了两个血液祖群体之间的分离，但在三个红细胞群体中没有（补充图16e）。
  12. 这表明这些不同的红细胞群体可能是批次效应而非生物学现象。
  13. 在基准测试中（补充结果），我们发现除了MNN Correct外，其他算法未能保持平滑性或错误地混合了不同的祖群体（补充图17）。

Harmony integrates dissociated scRNAseq with spatially resolved datasets

Harmony 将分离的单细胞RNA测序与空间解析数据集整合

  1. 现在我们考虑整合用不同模态测量的细胞，测量细胞生物学的不同方面。这种类型的整合特别强大，因为它允许推断在一种模态中测量但在另一种模态中未测量的特征。
  2. Harmony整合了一个空间解析的多重错误稳健荧光原位杂交（MERFISH）数据集和一个来自小鼠下丘脑前区（Bregma +0.5到−0.6 mm）的解离单细胞RNA测序（scRNAseq）数据集，使用154个共享基因（图6a）。
  3. 从嵌入中，我们推断未测量基因的空间定位，以识别神经元转录因子中先前未报告的空间模式。
  4. 我们下载了Moffit等人发表的两组数据：来自六只小鼠的30,370个细胞，使用解离scRNAseq（10×）测量，以及来自一只小鼠的64,373个细胞，使用MERFISH测量。
  5. 空间解析数据集是135个使用MERFISH检测的基因和20个使用双色单分子FISH（smFISH）检测的基因的组合。
  6. 解离的10×数据提供了关于完整转录组的信息（质量控制后为22,067个基因），但没有任何空间信息。
  7. 相比之下，MERFISH在该数据集中仅应用于有限数量的目标基因。
  8. 这项分析中的主要挑战是重叠特征的有限数量。成功的整合将合并相似的细胞类型，并允许推断剩余21,913个基因的空间模式。

Fig. 6: Harmony integrates spatially resolved transcriptomic with dissociated scRNAseq datasets.

  • 小鼠脑下丘脑前区细胞采用两种技术并行检测。对来自6只动物的30,370个解离细胞的完整转录组使用10X技术进行 profiling。对来自1只动物完整组织中的64,373个细胞，使用MERFISH技术在原位 profiling 了155个基因。比例尺，2毫米。
  • Harmony将两种模态的细胞整合到一个共享的嵌入中，正确合并了先前识别的12种细胞类型。源数据
  1. 在Harmony之前，细胞主要按模态分组（中位数iLISI 1.00，95% CI [1.00, 1.04]，见补充图18a）。
  2. Harmony混合了两种模态，10×和MERFISH（中位数iLISI 1.15，95% CI [1.00, 1.99]，见图6b和补充图18b），并根据它们之前的细胞类型标签合并细胞，覆盖了12个主要的神经元和胶质细胞群体（整合后预测准确率为93.7%，整合前为94%，见补充图18c）。
  3. 作为独立验证，我们展示了Harmony嵌入与原始论文中的细粒度聚类分析一致（补充结果和补充图18d,e）。
  4. 我们接下来使用核k近邻（kNN）插补方法，基于其最近的10X数据集邻居来预测MERFISH数据集细胞中的基因表达（方法部分）。
  5. 我们使用五折交叉验证来评估性能，发现我们准确预测了154个基因中的150个，对高表达基因的预测表现更为一致（补充结果和补充图18f）。
  6. 原始论文识别了区分空间分布的神经元亚型的关键转录因子。
  7. 应用Harmony后，我们可以轻松识别空间自相关的转录因子。
  8. 我们测量了JASPAR29中记录的所有已知脊椎动物转录因子的空间自相关性（Moran’s I）（补充表6和7）。
  9. 如预期的那样，所有与神经元亚型相关的因子在空间上显著（FDR < 5%）自相关。
  10. 此外，我们发现50个转录因子在兴奋性神经元中显著（FDR < 5%，Moran’s I > 0.1）自相关，而在抑制性神经元中有37个。
  11. 我们确认这些新因子至少是某一神经元亚型的阳性标记。
  12. 事实上，空间相关性与其最佳亚型的曲线下面积显著相关（Spearman’s r = 0.71，P = 2.6 × 10−62）（补充图18g）。
  13. 我们接下来对Satb1进行了跟进研究，这是一种与小鼠皮层中间神经元存活相关的染色质组织者，因为在抑制性神经元中它表现出强烈的自相关性（Moran’s I = 0.44）。
  14. 预测的Satb1表达定位在前部切片中（补充图18h）。
  15. 为了验证这一定位，我们将预测的表达与艾伦脑图谱中匹配的Satb1表达图像进行了比较。
  16. 这些图像选择了相似的前后位置（介于Bregma +0.8至−0.2毫米之间）和室形态（以绿色显示）。
  17. 测量的Satb1表达在定性上与预测表达的模式相似。
  18. 这154个基因是由作者精心挑选的，与下丘脑前区细胞类型在生物学上相关。
  19. 使用相关性指导的降采样方法，我们仅用60个代表性的锚基因就能获得相似的结果（补充结果和补充图18i,j）。

Discussion

  1. Harmony 接受一个细胞矩阵和每个细胞的协变量标签作为输入。
  2. 其输出是一个调整后的坐标矩阵，与输入矩阵具有相同的维度。
  3. 因此，Harmony 应该在应用完整的分析流程之前首先使用。
  4. 下游分析，如聚类、轨迹分析和可视化，然后可以使用调整后的 Harmony 嵌入。
  5. Harmony 不改变单个基因的表达值以解释数据集特定的差异。
  6. 我们建议使用一种批次感知的方法，例如混合效应线性模型，来进行差异表达分析。
  7. Harmony 使用聚类方法将细胞所在的大规模非线性嵌入空间划分为更小的线性区域。
  8. 采用离散聚类时，可能会过度离散化空间，从而丢失种群间的平滑过渡（例如，参见补充说明3）。
  9. Harmony 通过使用软聚类避免了过度离散化。
  10. 在这种策略下，在鼠胚胎发育数据集中，Harmony 有效模拟了终端种群和过渡状态，保留了从共同祖细胞到内皮和造血谱系的平滑过渡和分支事件。
  11. 随着高通量单细胞技术的出现，现在可以检测细胞生物学的不同相互作用方面。
  12. 有效整合跨多种模态的数据集可能有助于揭示这些相互作用。
  13. Harmony 将空间分辨的 MERFISH 细胞与分离的 10× 细胞嵌入，尽管这些技术之间基因重叠有限且检测效率差异显著。
  14. 因此，我们能够分析那些从未用 MERFISH 测量的基因的空间表达模式。
  15. Harmony 整合 10X 和 MERFISH 并估算未测量基因表达的能力依赖于重叠锚基因集的质量。
  16. 选择最佳的锚基因是一个重要的开放性问题。
  17. 通过我们的交叉验证和最小基因集分析，我们发现了两个可能有助于指导空间探针集设计的趋势。
  18. 首先，我们可以从锚基因中去除相关基因，而不会显著降低整合质量。
  19. 其次，较少丰度的基因更难预测。
  20. 因此，包含低丰度且最大程度非冗余基因的空间探针集应能通过 Harmony 与适当匹配的 scRNAseq 数据集估算大量基因。
  21. 在应用单细胞整合算法之前，通常会对基因进行批次敏感的预处理步骤。
  22. 特别是，一些研究人员会在将细胞合并到一个矩阵之前，分别对数据集中的基因表达值进行缩放。
  23. 虽然这种策略可能使得整合某些数据集（补充结果和补充图19a,b）变得更容易，在这些数据集中所有细胞群体都存在于所有数据集中，但当数据集由重叠但不完全相同的群体组成时（补充图19c），它可能会增加误差。
  24. 因此，在本文中我们不使用这种缩放策略。
  25. 在我们的分析流程中，我们避免所有批次敏感的预处理步骤。
  26. 我们只是简单地合并数据，并在组合的数据集上执行PCA。
  27. Harmony 模型能够识别并消除已知变异来源的影响。
  28. 然而，不希望的变异也可能来自未知来源。
  29. 像 SVA32 和 PEER33 这样的方法通过线性模型在大量转录组学中推断并消除潜在的变异来源。
  30. 在未来的工作中，我们计划扩展 Harmony 以识别和消除不希望的潜在影响。
  31. 随着自动化分类算法的兴起，用户在整合之前可能能够分配有关细胞类型的信息。
  32. 例如，过去我们的团队曾使用分类器来定义软性或硬性身份，以通过线性判别分析对单个细胞进行概率分类。
  33. Harmony有一个高级选项，可以使用这种概率细胞类型分配来初始化软聚类分配（参见HarmonyMatrix）。
  34. 如果之前的分配是有帮助的，Harmony将更快地收敛到准确的解决方案。
  35. 如果之前的分配不准确，Harmony可以在聚类步骤中适当地重新分配细胞。
  36. Harmony框架为一些令人兴奋的未来应用奠定了基础。
  37. 首先，它可以扩展以准确建模基因计数。
  38. 这将允许用户对需要完整基因表达谱的方法应用Harmony预处理，而不是低维细胞嵌入，例如RNA速度。
  39. 接下来，我们设想将Harmony专门化，以快速将细胞映射到一个覆盖全面组织、生物体和临床条件的数十亿细胞参考库。
  40. 这一应用将使单个实验中的细胞与更大的参考库进行快速比较，并在过程中几秒钟内注释已知和新的细胞类型和状态。

Methods

Harmony

和谐

Overview

概述

  1. Harmony算法输入细胞的PCA嵌入（Z）及其批次分配（ϕ），并返回一个批次校正后的嵌入（(\hat Z)）。
  2. 该算法，如下面的算法1所总结，在两个互补阶段之间迭代：最大多样性聚类（算法2）和基于混合模型的线性批次校正（算法3）。
  3. 聚类步骤使用批次校正后的嵌入(\hat Z)来计算细胞到簇的软分配，编码在矩阵R中。
  4. 校正步骤使用这些软簇来从原始嵌入计算一个新的校正嵌入。
  5. Harmony的高效实现，包括聚类和校正子程序，可作为R包的一部分在https://github.com/immunogenomics/harmony上获取。
  6. 请注意，校正程序使用的是Z，而不是(\hat Z)来消除混杂因素的影响。
  7. 通过这种方式，我们将校正限制在原始嵌入的线性模型中。
  8. 另一种方法会使用最后一次迭代的输出(\hat Z)作为校正程序的输入。
  9. 因此，最终的(\hat Z)将是原始嵌入的一系列线性校正的结果。
  10. 虽然这允许更具有表现力的转换，但我们发现，在实践中，这可能会导致数据过度校正。
  11. 我们选择限制转换反映了引言中的观点。
  12. 也就是说，如果我们在校正前对细胞类型有完美了解，我们会在每个细胞类型内线性回归消除批次效应。

Algorithm 1 Harmony

算法1 和谐

Glossary

词汇表

  1. 为了参考，我们定义了所有在Harmony函数中使用的数据结构。对于每一个数据结构，我们定义了其维度和可能的值，以及对其在上下文中含义的直观描述。维度以d表示，即嵌入的维度（例如，主成分的数量）；B表示批次的数量；N表示样本的数量；Nb表示批次b中的样本数量；K表示簇的数量。
  2. Z 属于 R 的 d×N 维空间，表示输入嵌入，将在 Harmony 中进行修正。这通常是细胞的 PCA 嵌入。
  3. （\hat Z \in {\Bbb R"}^{d\times N"}）由 Harmony 输出的集成嵌入。
  4. （R ∈ [0,1]^{K × N"}）细胞（列）到簇（行）的软聚类分配矩阵。每一列是一个概率分布，因此总和为1。错误！！！ - 待补充
  5. （Pr_b ∈ [0,1]^B）批次的频率。
  6. O属于[0,1]的K×B矩阵，表示在簇（行）和批次（列）中观察到的细胞共现矩阵。
  7. （E ∈ [0,1]^{K × B"}）在假设簇和批次分配独立的情况下，细胞在簇和批次中的预期共现矩阵。
  8. （Y属于[0,1]的d×K维矩阵）k-means聚类算法中的簇质心位置。

Assumptions about input data

关于输入数据的假设

  1. 在本文中，我们始终在细胞的低维嵌入上使用Harmony。
  2. 为了明确这个低维空间的性质，我们明确提出了三个假设：
  3. 细胞被嵌入到低维空间中，作为 PCA 的结果。PCA 嵌入捕捉了基因表达的变异，形成紧凑的正交空间。因此，Harmony 的默认输入现在是一个经过库大小标准化的基因表达数据矩阵。我们首先对这个高维矩阵进行 PCA，然后使用特征值缩放的特征向量作为低维嵌入输入到 Harmony 中。
  4. 基因表达已针对库大小进行了标准化。在 RNA-seq 中，每个细胞的测序深度不同，这导致每个细胞的库大小也不同。因此，在进行 PCA 之前，考虑这一技术变异来源是最佳实践。在本文中，我们使用了标准的每10,000个计数的对数转换（CP10K），具体方法在方法部分有描述。经过深度转换后，表达值被转换为每个细胞内的相对频率。因此，不可能在一组细胞中每个基因都表现出上调。
  5. 由欧几里得距离引导的低维最近邻结构应与常见的相似性度量（如余弦相似性和相关性）保持一致。这可以通过计算稀疏最近邻图并比较邻接矩阵轻松检查。邻居重叠少于 20% 的细胞可以被视为异常值。违反这一假设的常见方式是模拟围绕原点的细胞（即所有嵌入均为 0）。我们在实际的单细胞 RNA 测序数据中没有发现这种情况。该假设在使用余弦距离比较细胞的整合方法中是常见的，例如 MNN 校正和 Scanorama。

Maximum diversity clustering

最大多样性聚类

  1. 我们开发了一种聚类算法，以最大化簇内批次之间的多样性。
  2. 我们如下介绍这种方法。
  3. 首先，我们回顾了一个先前发表的用于软k-means聚类的目标函数。
  4. 然后，我们在此目标函数中添加了一个最大化多样性的正则化项，并推导出这个正则化项作为对两个随机变量之间统计依赖性的惩罚：批次成员资格和簇分配。
  5. 接着，我们推导并展示了用于优化目标函数的算法的伪代码。
  6. 最后，我们解释了实现的关键细节。

Background: entropy regularization for soft k-means

背景：软K-均值算法的熵正则化

  1. 经典K均值聚类的基本目标函数，其中每个单元格恰好属于一个簇，由单元格到其分配的质心的距离定义。
  2. 在上文中，Z是数据的某个特征空间，由质心Y共享。Rk,i可以取值0或1，表示细胞i属于簇k的成员资格。
  3. 为了将其转化为软聚类目标，我们遵循参考文献36的方向，并在R上添加一个由超参数σ加权的熵正则化项。
  4. 现在，Rki可以取值在0到1之间，只要对于给定的细胞i，簇成员资格的和Σk Rki等于1。
  5. 也就是说，Ri必须是一个在支持[1,K]上的正确概率分布。
  6. 当σ趋近于0时，这种惩罚趋近于硬聚类。当σ趋近于无穷大时，R的熵超过了数据质心距离。在这种情况下，每个数据点被平等地分配到所有簇中。
  7. 让我们一步步思考。

Objective function for maximum diversity clustering

最大多样性聚类的目标函数

  1. Harmony聚类算法的完整目标函数基于前一部分内容。除了软分配正则化之外，下面的函数还对所有定义的批次变量中批次多样性较低的簇进行惩罚。这种惩罚在下文中推导，取决于簇和批次的身份，Ω(R, ϕi) = Σi,k Rki log(Oki/Eki)ϕi。
  2. 对于每个批次变量，我们添加一个新的参数θ。
  3. 这里θ决定了批次成员资格和聚类分配之间依赖关系的惩罚程度。
  4. 当∀fθ = 0时，问题回归到(2)，对依赖关系没有惩罚。
  5. 随着θ的增加，目标函数更倾向于批次f和聚类分配之间的独立性。
  6. 当θ接近无穷大时，它将产生一个退化解。
  7. 在这种情况下，每个聚类在批次f中具有等效的分布。
  8. 然而，细胞与质心之间的距离可能较大。
  9. 最后，为了在梯度计算中便于表示，添加了σ这一项。
  10. 我们发现，当我们计算 Z 和 Y 之间的余弦距离，而不是欧几里得距离时，这种聚类效果最佳。
  11. Haghverdi 等人8 显示，当向量进行 log2 归一化时，平方欧几里得距离等同于余弦距离。
  12. 因此，假设所有 Zi 和 Yk 都具有单位 log2 范数，上述平方欧几里得距离可以重写为点积。

Cluster diversity score

集群多样性得分

  1. 在这里，我们讨论并推导了前一部分定义的多样性惩罚项 Ω(.)。
  2. 为简便起见，我们讨论单个批次变量的多样性，因为多个批次惩罚项在目标函数中是可加的。
  3. Ω(.) 的目标是惩罚批次身份与簇分配之间的统计依赖性。
  4. 在统计学中，两个离散随机变量之间的依赖性通常用 χ2 统计量来衡量。
  5. 这个测试考虑了两个随机变量的不同值一起出现的频率。
  6. 观察到的共现次数 (O) 与独立情况下预期的次数 (E) 进行比较。
  7. 出于实际原因，我们不直接使用 χ2 统计量。
  8. 相反，我们使用 Kullback Leibler 散度 (DKL)，这是两个分布之间的信息论距离。
  9. 在本节中，我们在概率簇分配矩阵 R 的背景下定义了 O 和 E 分布，以及 DKL 惩罚。
  10. 接下来，我们定义基于 R 的 KL 散度。请注意，O 和 E 都依赖于 R。
  11. 然而，在下面的推导中，我们展开了一个 O 项。
  12. 这在后面的优化过程中具有功能性的作用。
  13. 直观地讲，在单个细胞的 R 更新步骤中，我们计算所有其他细胞上的 O 和 E。
  14. 通过这种方式，我们决定如何根据当前批次在聚类中的分布情况，将单个细胞分配到聚类中。

Optimization

优化

  1. 方程（4）的优化允许采用期望最大化框架，迭代进行簇分配（R）和质心（Y）估计。

Cluster assignment R

簇分配 R

  1. 使用与参考文献36相同的策略，我们求解每个细胞i的最优分配Ri。首先，我们设置带有对偶参数λ的拉格朗日函数，并求解关于每个簇k的偏导数。
  2. 让我们一步步思考。
  3. 上面的分母项确保 Ri 的总和为1。在实际操作中（算法2），我们计算分子并将其除以总和。

Centroid estimation Y

质心估计 Y

Algorithm 2 Maximum Diversity Clustering

算法2 最大多样性聚类

Implementation details

实现细节

  1. 上述推导出的 R 和 Y 的更新步骤构成了最大多样性聚类的核心，如算法 2 所示。
  2. 本节解释该伪代码的其他实现细节。
  3. 同样，为简化起见，我们独立讨论与每个单一批次变量相关的多样性惩罚项 θ、ϕ、O 和 E 的细节。

Block updates of R

阻止R的更新

  1. 与常规的k-means不同，上述针对R的优化过程无法忠实地并行化，因为O和E的值会随着R的变化而变化。
  2. 因此，精确解依赖于在线过程。
  3. 为了提高速度，我们可以对这个过程进行粗粒度处理，并以小块（例如，5%的数据）更新R。
  4. 同时，O和E在保留的数据上计算。
  5. 在实践中，这种方法在块大小足够小的情况下成功最小化了目标函数。
  6. 在算法中，这些块被包含在for循环中的更新块中。

Centroid initialization

质心初始化

  1. 我们使用基础 R 语言中的 k-means 函数实现常规 k-means 聚类来初始化簇质心。
  2. 我们进行十次随机重启，并保留最佳结果。
  3. 然后，我们对质心进行 log2 归一化，以准备在算法 2 中进行球形 k-means 聚类，即最大多样性聚类。

Regularization for smoother penalty

用于平滑惩罚的正则化

θ discounting

θ 折扣

  1. 多样性惩罚，通过θ加权，确保来自一个批次的细胞在所有簇中均匀混合。这个假设在细胞数量较少的批次中更容易被打破。
  2. 批次越小，通过抽样论证，某些细胞类型在批次中未出现的可能性就越大。
  3. 将这样的批次分散到所有簇中会导致错误的聚类。
  4. 为了防止这种情况，我们为每个批次分配其自己的θb项，并根据批次中的细胞数量进行缩放。
  5. 上述θmax是未折现的θ值，适用于足够大的批次。
  6. 乘数因子(\left[ {1 - \mathrm{exp"}\left( { - N_b/K\tau } \right)^2} \right])的范围从0到1。
  7. 该因子在小批次大小Nb时呈指数级变化，而在足够大的Nb时趋于平稳。
  8. 超参数τ可以解释为从单个批次分配到每个簇的最小细胞数。
  9. 默认情况下，我们使用τ = 5到τ = 20之间的值。

K, the number of clusters

K，簇的数量

  1. 在和谐软聚类中使用的簇数 K 应设置为反映数据集大小和复杂性的值。
  2. 簇数太少将无法捕捉到生物上不同的细胞类型和状态的数量。
  3. 簇数太多将过分强调批次特定的异常值，并阻碍有效整合。
  4. 作为一个启发式方法，我们假设数据集最多有 100 种不同的细胞类型，并且每个簇至少应有 30 个细胞。
  5. 我们将默认的簇数 K 设置为介于这两个极端之间，针对 N 个细胞。

Linear mixture model correction

线性混合模型校正

  1. 在本节中，我们将所有需要从原始嵌入中整合出去的效果称为批次效应。
  2. 这并不意味着这些效应纯粹是技术性的。
  3. 这种术语只是为了方便使用。

Mixture of experts model

专家混合模型

  1. 一旦我们定义了批次多样化的簇，我们希望从每个簇中去除批次特定的变异。
  2. 我们通过改进Jordan和Jacobs38提出的原始专家混合模型来实现这一点。
  3. 在Harmony的背景下，每个细胞被概率性地分配到一个小的专家集合中。
  4. 这种分配在之前的聚类步骤中已经计算完成。
  5. 在给定簇/专家的情况下，专家混合模型假设响应变量和自变量之间存在线性关系。
  6. 因此，我们在簇/专家k的条件下，为响应变量定义了一个高斯概率分布。
  7. 在这里，均值是独立变量的函数（μk = WkTϕ*i），而协方差则不是（σk2I）。
  8. 请注意，上述设计矩阵（ϕ*）与聚类步骤中使用的设计矩阵（ϕ）不同。
  9. 我们对原始设计矩阵ϕ进行扩展，以包括一个截距项：ϕ* = 1 || ϕ。
  10. 这些截距项捕捉了每个簇或专家中批次独立（即细胞类型）的变异。
  11. 我们还可以通过修改ϕ*来实现更复杂的行为，例如参考映射（Reference mapping）。
  12. 通过这种生成式公式，我们可以独立地为每个簇或专家求解线性模型的参数（Wk）。
  13. 上述 diag(Rk) 是簇 k 的成员项的对角矩阵。Z 是原始 PCA 嵌入的矩阵
  14. Wk 中的每一列对应一个主成分维度，从 PC1 到 PCd。Wk 的第一行对应批次独立的截距项。随后的行（1 到 B）对应于原始设计矩阵 ϕ 中的一热编码批次分配。
  15. 为了获得细胞特定的校正值，我们根据簇分配概率分布 R 对 Wk 进行期望计算。具体来说，每个细胞由一个批次独立的截距建模，该截距表示在 Wk 的第一行（Wk[0,.]）中，而其批次依赖的项则由剩余行表示：Wk[1:B,.]。
  16. 剩下的是细胞类型特定的截距 Wk[0,.] 和细胞特定的残差 ϵi。
  17. 不幸的是，对于设计矩阵ϕ*，公式(13)中没有解。
  18. 这是因为ϕ不是满秩的，因此ϕ diag(Rk) ϕ*T不可逆。
  19. 这种奇异性源于一个独热编码的分类变量之和等于截距的事实。
  20. 为了解决这种共线性，我们对Wk中的非截距项施加log2范数惩罚，类似于岭回归。
  21. 这会使Wk项收缩到0。
  22. 就像在岭回归中一样，我们不是求逆ϕ* diag(Rk) ϕT，而是求逆ϕ diag(Rk) ϕ*T + λI，这是非奇异的。
  23. 现在Wk的解变为

Algorithm 3 Mixture of Experts Correct

算法3 专家混合校正

  1. function correct (Z, R, ϕ)
  2. let's think step by step.
  3. （\hat Z \leftarrow Z）
  4. （φ^ * ← 1||φ）

Reference mapping

参考映射

Caveat

警告

Performance and benchmarking

性能和基准测试

LISI metric

LISI 指标

  1. 评估批量校正和数据集整合过程中的混合程度是一个开放性问题。
  2. 几个研究组已经提出了方法来量化嵌入空间内由kNN图定义的局部邻域中的批次多样性。
  3. Buttner等人提供了一个统计测试来评估混合程度，而Azizi等人报告了这些分布的熵。
  4. 我们的局部多样性度量与这些方法相关，因为我们从kNN图开始。
  5. 然而，我们的方法考虑了这些方法未涉及的两个问题。
  6. 首先，该度量标准应对局部距离更加敏感。例如，一个包含100个细胞的邻域可以在四个批次中均匀混合。
  7. 然而，在邻域内，细胞可能按批次聚类。
  8. 第二个问题是解释问题。
  9. kBET提供了一个统计测试来评估混合的显著性，但不清楚当数据集具有截然不同的细胞类型比例时，所有邻域是否都应显著混合。
  10. Azizi等人使用熵作为多样性的度量，但不清楚如何解释编码邻域分布所需的比特数。
  11. 我们的多样性得分，LISI，解决了这两个问题。
  12. 为了对局部多样性敏感，我们构建了基于高斯核的邻域分布。
  13. 这为邻域中的细胞提供了基于距离的权重。
  14. 当前实现使用固定的困惑度（默认30）来计算这些局部分布，这已被证明比固定邻居数量更为平滑的函数。
  15. 我们通过使用逆辛普森指数（1/Σb=1Bp(b)）来解决第二个解释问题。
  16. 这里的概率是指上述局部邻域分布中的批次概率。
  17. 该指数是在抽取到两个来自同一批次的细胞之前，预期需要抽取的细胞数量。
  18. 如果邻域仅由一个批次组成，则只需要一次抽取。
  19. 如果它是两个批次的等量混合，则平均需要两次抽取。
  20. 因此，该指数报告了局部邻域中的有效批次数量。
  21. 我们的多样性得分，LISI，结合了这两个特征：基于困惑度的邻域构建和逆辛普森指数。
  22. LISI为每个细胞分配一个多样性得分。
  23. 这个得分是该细胞邻域中的有效批次数量。
  24. 计算LISI的代码可在https://github.com/immunogenomics/LISI获取。
  25. LISI的主要缺点在于处理大小差异极大的数据集的情况。
  26. 我们在补充说明2中展示了这种情况下的结果。

Significance between embeddings

嵌入之间的意义

Time and memory

时间和记忆

  1. 我们对所有分析进行了执行时间和最大内存使用量的基准测试。
  2. 所有任务都在Linux服务器上运行，分配了六个核心和120 GB的内存。
  3. 这些机器配备了Intel Xeon E5-2690 v.3处理器。
  4. 为了评估执行时间和最大内存使用量，我们使用了Linux的time工具（在我们的系统上是/usr/bin/time）并带有-v标志来记录内存使用情况。
  5. 执行时间是从{Elapsed time"}字段记录的。
  6. 最大内存使用量是从{Maximum resident set size"}字段记录的。

Cell-type prediction accuracy

细胞类型预测准确度

  1. 除了cLISI之外，我们还通过细胞类型预测来测量嵌入的准确性。
  2. 简要地说，我们基于每个细胞的邻近细胞来预测其类型，并将准确性计算为正确预测的频率。
  3. 具体而言，我们基于余弦距离计算一个细胞的30个最近邻。
  4. 然后，我们使用宽度为σ = 0.1的径向基函数核对其进行加权。
  5. 最后，我们对邻近细胞的细胞类型进行加权求和，并返回最可能的细胞类型。
  6. 形式上，对于细胞i，其30个最近邻细胞（NN(i)），嵌入矩阵Z和细胞类型标签向量T，计算最可能的细胞类型标签t：

Data availability

  1. 本文分析的所有数据均可通过在线资源公开获取。
  2. 我们在补充表8中提供了所有数据源的链接。

Code availability

  1. Harmony 和 LISI 作为 R 包可在 https://github.com/immunogenomics/harmony 和 https://github.com/immunogenomics/lisi 上获取。
  2. 用于重现主要分析结果的脚本将可在 https://github.com/immunogenomics/harmony2019 上获取。
  3. 此外，附注作为补充说明包含在内。
  4. 补充说明 1 提供了 Harmony 的详细操作指南，将理论算法组件与其代码实现相连接。
  5. 补充说明 2 展示了 LISI 指标及其统计显著性的评估方法。
  6. 补充说明 1 使用了模拟数据集来演示 Harmony。

Change history

[ul]- 26 August 2020 In the supplementary information originally posted for this article, the Supplementary Results and Supplementary Notes 1–3 were missing. The error has been corrected online.