ArchR实现了LSI降维算法

ArchR_Logo_Integrated

数据的稀疏性使得对scATAC-seq数据降维十分困难。在二倍体的scATAC-seq数据中,一个位点有三种测量状态:一个allele被测量、两个allele被测量、两个allele都没有被测量。即使是在高质量的scATAC-seq数据里,大部分的开放区域也没有被转座酶转座,这就导致了大量位点的检测丰度为0,即两个allele都没有被测量到。此外,例如当我们在某个单细胞的某个peak中观察到3个转座位点时,由于数据的稀疏性,我们并不能保证说这个peak的开放程度就比那些其他细胞中只观察到一个转座位点的相同位置的peak的开放程度大。因此基于这个原因,大量的scATAC-seq数据的分析方法都是基于二值化后的数据矩阵进行的,即将原始的计数矩阵转化为01矩阵。这个二值化后的矩阵的大量元素仍然是零,因为我们测量到的开放区域相对于细胞实际开放的区域来说要少得多。然而这里我们需要明白,二值化矩阵中的零既可能表示这个区域本身没有开放,或者实验没有检测到它开放,而这两者本身就很难进行区别。显然,第一种情况下的零不包含任何开放信息,而第二种情况下的零包含了开放信息,但是实验却没有捕获到。正是这种只包含了占开放区域总体极少比例的特性造成了scATAC-seq数据的稀疏性。

如果你对这个稀疏矩阵进行常规的降维操作例如PCA,并画出PC1和PC2的关系图,你不会得到想要的结果,因为大量的零使得所有细胞间的相似度都很高。为了解决这个问题,ArchR使用了分层降维的方法。首先,ArchR使用了潜在语义搜索(Latent Semantic Indexing, LSI),这是一种在自然语言处理领域广泛应用的算法,基于单词的计数评估文档的相似性。自然语言中的文档跟ATAC数据一样,也具有很强的稀疏性和噪音。一个文档中往往会出现大量的单词和大量的低频单词。LSI最早是由Cusanovich et al. (Science 2015)引入scATAC-seq数据处理中来的。在scATAC-seq数据中,我们将不同的细胞看作不同的文档,将不用的染色质区域或者peak看作自然文档中的单词。首先我们在单个细胞内对peak的fragments覆盖数目进行正则化以计算每个peak的频率。这些频率值通过反向文档频率(Inverse Document Frequency, IDF)再次进行正则化。IDF通过特征出现的次数估计特征对于不同文档识别的重要性,显然那些出现频率较低的单词在某个文档中特异性出现的机率很大,而这样的单词对于不同文档间的识别就十分重要,因为那些出现次数很高的单词大概率在每个文档中出现的次数都很高,它们对不同类文档间的识别帮助较小。最终计算出的特征频率-反转文档频率(TF-IDF)矩阵能够反应出这个单词(或者叫region/peak)对文档(或者叫细胞)的重要程度。然后,通过奇异值分解(SVD)能够识别细胞间有价值的信息,并在低维空间中加以展示。LSI能够极大的降低稀疏的(二值化的)计数矩阵的维度,然后使用TSNE或者UMAP进行可视化。ArchR将这些可视化方法叫做embeddings。

ArchR对LSI的实现

ArchR实现了几种不同的LSI实现,并在多个不同的测试数据集上进行了基准测试。ArchR默认的LSI与Signac使用的LSI相似:use a term frequency that has been depth normalized to a constant (10,000) followed by normalization with the inverse document frequency and then log-transforming the resultant matrix (aka log(TF-IDF))。

LSI降维只需要一个初始矩阵。到目前为止,scATAC-seq数据降维的两个主要策略是:使用peak矩阵或者由全基因组切分而成的tile矩阵(SnapATAC中叫做bin矩阵)。然而,对peak矩阵使用LSI降维本身就很有挑战性,因为peak往往是针对单个cluster富集出来的,而聚类是在降维之后,所以在对peak矩阵降维之前,我们并没有相关的聚类信息。如果我们不进行聚类而直接使用整个数据集上富集出来的peak,这些peak不能反映出不同cluster之间的差异性,所以从逻辑上来说,使用peak矩阵进行降维几乎是行不通的。另一方面,peak的富集严重依赖于样本,当数据集的样本发生变化时,peak的集合也可能发生变化,这导致peak矩阵极不稳定。相比之下,第二种策略可行性要强得多。对全基因组进行切片(genome-wide tiles),即将全基因组分成相邻等距的窗口,计算每个窗口中的fragments数量,ArchR中称为tile矩阵,SnapATAC中叫做bin矩阵。另一个显而易见的事实是,这样会得到大量的窗口,而其中很大一部分窗口是没有意义的,例如计数为0或者常数,位于blacklist区域等等,如果不加过滤会使得初始矩阵格外的庞大。解决这个问题一般有两个途径:控制窗口大小(一般不小于5kbp)和对窗口进行过滤。显然前者是控制窗口数量的有效措施,后者更著重于提高窗口的质量。切片法 也有其缺点,因为大多数染色质开放区域实际上只有几百碱基对的长度,而我们的窗口大小都是几千碱基对长度,这使得tile矩阵对开放区域并不是特别敏感,使用tile矩阵降维的方法分辨率会更低。当然解决办法也是有的——减小窗口大小,但是会耗费更长的计算时间和更多的内存。

ArchR由于设计了箭头文件,能够较快的计算较小窗口的tile矩阵,因此ArchR中窗口的默认大小是500bp。较小的窗口大小能够缓解tile矩阵对开放区域不敏感的缺陷。伴随而来的另一个问题是,500bp大小的窗口会使得切分出的窗口数量多达600万个。ArchR实现了数据的分块读取,同时也实现了另一种 预估式的LSI(estimated LSI) 方法对读取的子集进行预降维。这种预估式的LSI主要有两个用途:

  1. 加速降维
  2. 当减少用于(单次)降维的细胞数量时,相当于减少了数据的 粒度(granularity)

原文中使用了 粒度 这个词,大意应该是对数据的细分程度。对数据总体分的块越多,意味着每次读取的细胞数量越少,粒度越小。粒度代表了单位IO中的样本数量。减少数据的粒度在一定程度上可以缓解 批次效应。但是,过小的粒度也可能会消解掉真实的生物学差异,所以使用这种方法应该十分小心。

----- For reprint please indicate the source -----
0%