scATAC分析神器ArchR初探-使用ArchR降维（4）

4-使用ArchR降维

由于稀疏性，使用scATAC-seq降低尺寸具有挑战性的数据。在scATAC-seq中，可以在一个等位基因（两个等位基因或一个等位基因）上访问特定位点。即使在高质量的scATAC-seq数据中，大多数可访问区域也不会转座，这导致许多具有0个可访问等位基因的基因座。此外，当我们在单个单元格的单个峰区域内看到（例如）三个Tn5插入片段时，数据的稀疏性使我们无法自信地确定该单元格中的该位点实际上比只有一个单元格的另一个单元格高出三倍在同一站点插入一个。因此，许多分析策略都对二进制化的scATAC-seq数据矩阵起作用。由于转置很少，因此该二值化矩阵最终仍大部分为0。然而，重要的是要注意，scATAC-seq中的0可能表示“不可访问”或“未采样”，并且从生物学的角度来看，这两个推论有很大不同。因此，1具有信息，而0没有。如此低的信息量是我们的scATAC-seq数据的来源稀疏的。

如果要在此稀疏插入计数矩阵上执行标准降维（例如主成分分析）并绘制前两个主成分，则将无法获得所需的结果，因为稀疏性导致所有0都具有较高的小区间相似度职位。为了解决这个问题，我们使用分层降维方法。首先，我们使用潜在语义索引（LSI），这是自然语言处理中的一种方法，最初旨在根据字数评估文档的相似性。该解决方案是为自然语言处理而创建的，因为数据稀疏且嘈杂（许多不同的单词和许多低频单词）。LSI由Cusanovich等人首次针对scATAC-seq引入。（2015年科学）。对于scATAC-seq，不同的样本是文档，不同的区域/峰是单词。首先，我们通过每个单元的深度归一化来计算频率项。然后，通过反文档频率对这些值进行归一化，反文档频率通过对特征进行加权的频率来对特征进行加权，以识别更“特定”而不是通常可访问的特征。最终的词频-反文档频率（TF-IDF）矩阵反映了单词（即区域/峰值）对文档（即样本）的重要性。然后，通过一种叫做奇异值分解（SVD），最有价值跨样本的信息被识别并在较低维度的空间中表示。LSI使您可以将稀疏插入计数矩阵的维数从数千减少到数十或数百。然后，可以使用更常规的降维技术，例如均匀流形近似和投影（UMAP）或t分布随机邻居嵌入（t-SNE），将数据可视化。在ArchR中，这些可视化方法称为嵌入。

4.1 ArchR的LSI实现

ArchR实现了几种不同的LSI实现，我们已经在多个不同的测试数据集中对许多方法进行了基准测试。ArchR的默认LSI实现与Timacy Stuart在Signac中引入的方法有关，该方法使用的术语频率已被深度归一化为常数（10,000），然后使用反文档频率进行归一化，然后对结果矩阵进行对数转换（aka log(TF-IDF)）。

降低LSI尺寸的关键输入之一是起始矩阵。到目前为止，scATAC-seq中的两个主要策略是（1）使用峰区域或（2）全基因组图块。但是，将峰区域用于LSI本身就具有挑战性，因为在降维之前我们没有簇或簇特定峰。此外，在聚类之前在聚集的细胞上调用峰会掩盖特定于细胞类型的峰。而且，将新样品添加到实验中时，任何联合峰集都会改变，从而使该策略的稳定性降低。第二种策略是使用全基因组切片，通过使用一致且无偏的特征集（全基因组切片）来缓解这些问题。但是，所有区域中所有细胞的全基因组图块矩阵可能会变得过大。为此原因，大多数实现都使用大于或等于5 KB的小块。因为大多数可访问区域只有几百个碱基对长，所以这大大降低了方法的分辨率。

由于Arrow文件的设计方式，ArchR能够使用全基因组范围的500 bp切片快速执行LSI。这解决了分辨率问题，并允许在调用峰之前识别簇。挑战在于，500 bp的条带生成约600万个特征，并按图块矩阵包含在单元中。虽然ArchR可以通过分块相关矩阵将大量数据读取到R中，但我们还实现了“估计LSI”方法，该方法对所有单元的子集执行初始降维。这种估计的LSI方法具有两个主要用途-（i）加快降维速度；（ii）减少初始降维时使用的单元数，这会降低数据的粒度。粒度的减少可用于您的优势，以减少数据中的批量影响。但是，它也可能掩盖真实的生物学信息，因此应在密切的人工监督下使用估计的LSI方法。

4.2迭代潜在语义索引（LSI）

在scRNA-seq中，识别可变基因是计算降维的常用方法（例如PCA）。这样做是因为这些高度可变的基因在生物学上更可能具有重要意义，并且可以减少实验噪音。在scATAC-seq中，数据是二进制的，因此您无法识别可变峰以降低维数。我们没有确定最大的峰，而是尝试使用最易访问的功能作为LSI的输入。但是，运行多个样品时的结果显示出较高的噪声水平和低重现性。为了解决这个问题，我们引入了``迭代LSI''方法（Satpathy *，Granja *等人，Nature Biotechnology 2019和Granja *，Klemm *和McGinnis *等人，Nature Biotechnology 2019）。此方法在最易访问的图块上计算初始LSI转换，并标识没有批次混淆的较低分辨率的群集。例如，当对外周血单核细胞进行检测时，这将鉴定出与主要细胞类型（T细胞，B细胞和单核细胞）相对应的簇。然后，ArchR计算所有要素中所有这些集群的平均可访问性。然后，ArchR识别这些群集中变化最大的峰，并将这些功能再次用于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
)

如果您在下游看到具有细微的批处理效果，则另一个选择是添加更多的LSI迭代，并从较低的初始群集分辨率开始，如下所示。另外，可变特征的数量可以减少以增加对更多可变特征的关注。

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
)

4.3估计的LSI

对于非常大的scATAC-seq数据集，ArchR可以估计带有LSI投影的LSI维数减少。此过程与迭代LSI工作流程相似，但是LSI过程不同。首先，将随机选择的“地标”单元的子集用于LSI降维。其次，使用从界标单元确定的反文档频率对其余单元进行TF-IDF归一化。第三，将这些归一化的单元投影到由界标单元定义的SVD子空间中。这导致了基于一小部分单元的LSI转换，这些单元被用作剩余单元的投影的界标。这种估计的LSI程序对于ArchR是有效的，因为在将新单元投影到界标单元LSI中时 ArchR迭代地从每个样本读取单元，而LSI将其投影而不将其全部存储在内存中。这种优化导致最小的内存使用，并进一步提高了超大型数据集的可伸缩性。重要的是，所需的界标集大小取决于数据集中不同像元的比例。

addIterativeLSI()通过设置sampleCellsFinal和projectCellsPre参数，可以通过功能在ArchR中访问估计的LSI 。samplesCellsFinal指定界标单元格子集的大小，并projectCellsPre告诉ArchR使用此界标单元格子集进行其余单元格的投影。

4.4通过谐波校正批效应

有时，迭代LSI方法不足以解决强批处理差异。因此，ArchR实现了一种常用的批处理校正工具Harmony，该工具最初是为scRNA-seq设计的。我们提供了一个包装器，该包装器会将降维对象从ArchR直接传递给HarmonyMatrix()函数。附加参数可以HarmonyMatrix()通过附加参数（...）直接传递到函数中。请参阅?addHarmony()以获取更多详细信息。用户应了解针对其特定应用的批处理纠正的注意事项。

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

参考材料：

https://www.archrproject.com/