ArchR官网教程学习笔记4:ArchR的降维

系列回顾：
ArchR官网教程学习笔记1：Getting Started with ArchR
ArchR官网教程学习笔记2:基于ArchR推测Doublet
ArchR官网教程学习笔记3:创建ArchRProject

这一章作者用了大量的篇幅来解释ArchR是如何进行降维的。其中有很多专业术语，我看着也困难，翻译出来感觉也不是人话。。。个人觉得看得懂更好，看不懂也不用太纠结，主要是会运行代码，知道有这么个步骤就行了。

由于数据的稀疏性，scATAC-seq的降维非常具有挑战性。在scATAC-seq中，一个特定的位点的可接近性可以在一个等位基因上，两个等位基因上，或者没有等位基因。即使在高质量的scATAC-seq数据中，大多数可接近的区域没有调换（transposed），导致许多loci的可接近的等位基因是0。此外，当我们看到(例如)在单个细胞中三个Tn5插入一个peak里，数据的稀疏性使我们怀疑在这个位点上是否真的有3倍高的可接近性（相比于其他细胞）。因此，许多分析策略都在二进制的scATAC-seq数据矩阵上工作。这个二进制的矩阵仍然大部分是0。然而，需要注意的是，scATAC-seq中的0可能意味着“不可接近”或“未被采样”，而且从生物学角度来看，这两种推断是非常不同的。正因为如此，1代表有信息而0没有。这种低信息量使得scATAC-seq数据变得稀疏。

如果你要在这个稀疏的矩阵上执行一个标准的降维操作，如主成分分析（PCA），并绘制前两个主成分，你不会得到想要的结果，因为稀疏性导致在所有的0位置上的细胞间高度相似。为了解决这个问题，我们使用分层降维方法。首先，我们使用Latent Semantic Indexing(LSI)，一种来自处理自然语言的方法，最初设计用于基于字数相似度来评估文档的相似性。这个解决方案是为自然语言处理创建的，因为数据是稀疏的和noisy的(许多不同的单词和许多低频率的单词)。Cusanovich等人Cusanovich et al. (Science 2015)首次将LSI引入scATAC-seq中。在scATAC-seq中，不同的样本是documents，不同的区域/峰是words。首先，我们通过每个单细胞的深度归一化来计算频率。然后，这些值按照反向文档频率进行规范化（文档频率根据特征出现的频率对特征进行衡量，以识别更“特异的”而非commonly可访问的特征）。由此得到TF-IDF矩阵，这个矩阵反映了一个word(即区域/峰)对document(即样品)的重要性。然后，通过奇异值分解(SVD)技术，识别样本间最有价值的信息，并在较低维空间中表示出来。LSI允许你将稀疏矩阵的维度从几千减少到几十或几百。在此基础上，采用UMAP或t-SNE等传统的降维技术来实现数据的可视化。在ArchR中，这些可视化方法称为嵌入（embeddings）。

（一）ArchR的LSI的实现

ArchR实现了一些不同的LSI方法，并且我们已经在多个不同的测试数据集上对其中的许多方法进行了基准测试。ArchR的默认的LSI实现与Timothy Stuart在Signac中引入的方法有关，该方法使用一个term频率，该频率被深度归一化为常数(10,000)，然后使用逆文档频率归一化，再对结果矩阵进行对数变换(又名log(TF-IDF))。

LSI降维的关键输入之一是初始矩阵。到目前为止，scATAC-seq的两种主要策略是(1)使用峰区或(2)使用全基因组。然而，在峰区域使用LSI本质上是具有挑战性的，因为在降维之前，我们没有进行聚类（cluster）或集群特定的峰（cluster-specific peaks）。此外，在clustering之前call peaks会模糊细胞特异类型的峰。而且，在实验中加入新样本时，任何union峰set都会发生变化，使得该策略不太稳定。第二种策略是使用全基因组分片，通过使用一致且无偏差的特征集(全基因组分片)来解决这些问题。然而，对所有细胞的全基因组分片会使矩阵过大。由于这个原因，大多数LSI实现使用size大于或等于5kb的“片”。这大大降低了该方法的分辨率，因为大多数可接近区域只有几百个bp长度。

由于Arrow files的设计方式，ArchR能够使用全基因组500bp分片非常快速地执行LSI。这解决了分辨率的问题，并允许在call peaks之前鉴定clusters。问题是500bp的bin会生成大约600万个特征。虽然ArchR能够通过将矩阵分组将大量数据读取到R中，但我们也实现了一种“estimated LSI”方法，该方法对总细胞的一个子集执行初始维数缩减。这种estimated LSI方法有两个主要的用途- (i)加速降维，(ii)随着细胞数的减少，将减少数据的粒度（granularity）。这种粒度的减少可以用于减少数据中的批次效应。然而，它也可能掩盖真实的生物学效应，所以estimated LSI方法应该在人工监督下使用。

（二）迭代LSI

在scRNA-seq中，识别variable基因是计算降维的常用方法(如PCA)。之所以这样做，是因为这些高度可变的基因在生物学上有可能是很重要的，这就减少了实验噪音。在 scATAC-seq 中，数据是二进制的，因此你不能识别variable峰来进行降维。我们尝试使用最易接近的特性作为LSI的输入。然而，当运行多个样本时，结果显示噪音高，重现性低。为了解决这个问题，我们引入了“迭代LSI”方法（Satpathy, Granja et al. Nature Biotechnology 2019和Granja, Klemm and McGinnis* et al. Nature Biotechnology 2019）。这种方法在最易接近的分片上计算初始LSI转换，并识别不存在批处理混淆的低分辨率clusters。例如，当对外周血单核细胞进行检测时，将识别出与主要细胞类型(T细胞、B细胞和单核细胞)相对应的簇。然后ArchR计算每个cluster所有特征的平均可接近性。之后，ArchR识别出这些cluster中变化最高的峰，并再次将这些特性用于LSI。在第二次迭代中，变化最高的峰更类似于在scRNA-seq LSI实现中使用的可变基因。用户可以设置应该执行多少次LSI迭代。我们发现这种方法可以使观察到的批次效应最小化，并允许在一个合理大小的特征矩阵上进行降维操作。

为了在ArchR中执行迭代LSI，我们使用addIterativeLSI()函数。默认参数应该涵盖大多数情况，但我们鼓励你探索可用的参数，以及它们如何影响你的数据集。使用?addIterativeLSI查看更详细的说明。最常见的调整参数是iterations、varFeatures和resolution。重要的是要注意到LSI是不确定性的。这意味着即使你以完全相同的方式使用完全相同的参数运行LSI，你也不会得到完全相同的结果。当然，它们高度相似，但不完全相同。因此，一旦确定了理想的降维，请确保保存你的ArchRProject或相关的LSI信息。

为了练习使用，我们要创建一个reducedDims对象，称为“IterativeLSI”：

> projHeme2 <- addIterativeLSI(
  ArchRProj = projHeme2,
  useMatrix = "TileMatrix", 
  name = "IterativeLSI", 
  iterations = 2, 
  clusterParams = list( #See Seurat::FindClusters
    resolution = c(0.2), 
    sampleCells = 10000, 
    n.start = 10
  ), 
  varFeatures = 25000, 
  dimsToUse = 1:30
)

Checking Inputs...
ArchR logging to : ArchRLogs\ArchR-addIterativeLSI-28e836691ea-Date-2020-11-20_Time-00-32-15.log
If there is an issue, please report to github with logFile!
2020-11-20 00:32:17 : Computing Total Across All Features, 0.004 mins elapsed.
......
2020-11-20 00:35:38 : Finished Running IterativeLSI, 3.357 mins elapsed.

如果您看到下游有细微的批次效应，可做的是添加更多的LSI迭代，并从较低的初始clustering分辨率开始，如下所示。此外，可变特征的数量可以降低，以增加对更多的可变特征的关注。

出于举例的目的，我们将把这个reducedDims对象命名为“IterativeLSI2”，但是我们不会在下游分析里使用它：

> projHeme2 <- addIterativeLSI(
    ArchRProj = projHeme2,
    useMatrix = "TileMatrix", 
    name = "IterativeLSI2", 
    iterations = 4, 
    clusterParams = list( #See Seurat::FindClusters
        resolution = c(0.1, 0.2, 0.4), 
        sampleCells = 10000, 
        n.start = 10
    ), 
    varFeatures = 15000, 
    dimsToUse = 1:30
)

（三）Estimated LSI

对于非常大的scATAC-seq数据集，ArchR可以利用LSI投射来估计LSI降维。这个过程类似于迭代LSI的工作流，但是LSI过程有所不同。首先，使用随机选择的“landmark”细胞的子集进行LSI降维。其次，使用从“标志”细胞确定的反向文档频率对其余细胞进行TF-IDF规范化。第三，将这些归一化的细胞投影到由“标志”细胞定义的SVD子空间中。这导致了基于小细胞群的LSI转换为剩余细胞的“标志”。这种estimated LSI过程在ArchR中是有效的，因为当将新的细胞投射到LSI的“标志”细胞时，ArchR迭代地从每个样本中读取细胞，而LSI不将它们全部存储在内存中。这种优化使得内存消耗最小，并进一步提高了超大数据集的可伸缩性。重要的是，所需的“标志”细胞集大小取决于数据集中不同细胞的比例。

ArchR通过addIterativeLSI()函数利用设置sampleCellsFinal和projectCellsPre参数来访问Estimated LSI。samplesCellsFinal指定“标志”细胞群的大小，而projectCellsPre告诉ArchR使用该“标志”细胞群对其余细胞进行投影。

（四）使用Harmony进行批次效应矫正

有时迭代LSI方法不足以校正强的批次效应。因此，ArchR实现了一个常用的批次效应修正工具Harmony，它最初是为scRNA-seq设计的。我们提供了一个封装器，它将降维对象从ArchR直接传递给HarmonyMatrix()函数。其他参数可以通过附加参数(…)直接传递给函数中的HarmonyMatrix()。更多细节请参见?addHarmony()。

> projHeme2 <- addHarmony(
  ArchRProj = projHeme2,
  reducedDims = "IterativeLSI",
  name = "Harmony",
  groupBy = "Sample"
)

Harmony 1/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony converged after 3 iterations

这个过程会在我们的projHeme2对象里创建一个新的reducedDims对象，叫“Harmony”。