10X单细胞（10X空间转录组）多样本批次效应去除分析之RCA2

Posted 2023-03-21

参考技术A

Data processing can also be carried out with Seurat. Here is an example how you can combine a RCA analysis with data preprocessed in Seurat.

Using the same 10x data as before, we generate a Seurat object and perform an initial analysis:

To run RCA, no further processing steps would be needed. However, we want to also compare the RCA result to the Seurat based clustering, therefore we first go on with a Seurat based analysis:

Based on the Elbowplot (not shown here), we use 20 PCs for further analysis.

We generate a UMAP of the data stored in the Seurat object using the umap R package:

We use the RCA function createRCAObject to generate a RCA object from the raw and optionally also the normalized data stored in our Seurat object.

Next, we can compute the projection, cluster the data, and estimate the most likely cell type for each cell as above:

Using the RCA cell type labels, RCA and Seurat clusters, we generate two new UMAPs whose coordinates are based on the PCs derived from HVGs and that are colored according to RCA clusters and cell type labels.

The RCA clusters show a high concordance to the Seurat clusters shown in the previous UMAP.

For greater convenience the results of RCA can be saved within the Seurat object for further analysis.

Also, a UMAP reduction based on the projection space can be added to the Seurat object:

RNA velocity describes the rate of gene expression change for an individual gene at a given time point based on the ratio of its spliced and unspliced messenger RNA (mRNA). Here, we describe how one can use the scvelo package, in Python, to visualize RNA velocity on the RCA generated result.

To transfer spliced RNA counts to scvelo, first transpose the raw RCA data matrix to get a cells x genes matrix, and export it to a CSV file.

In addition, export the RCA projection and UMAP embeddings to respective CSV files too.

Create an iPython notebook in the same folder and import the required packages as below.

Then, create a Scanpy object using the raw counts from the CSV file.

Populate the PCA slot in the Scanpy object as the projection data from RCA.

Populate the UMAP slot in the Scanpy object as the umap coordinates from RCA.

Load the unspliced loom object generated by velocyto .

Then, merge the spliced and unspliced objects together as described below:

As recommended by the scvelo tutorial, perform the following steps to compute RNA velocity:

It is possible that not all barcodes had sufficient quality of both spliced and unspliced reads, and thus some cells may have been discarded during the merging process. To ensure your cell type labels are still maintained, export the merged data observations from the merged scvelo object to a CSV file.

In R, load this CSV file in and extract the RCA labels and filter only those which were considered in the merged data by scvelo.

Note: If your cell names have underscores in them, scanpy will automatically split the cell name into barcode and sample_batch.

In this case, replace the last line of the above block of code with the following:

Now export these cluster labels to a CSV file.

Back in the scvelo iPynb, load this RCA cluster annotation table and set it as the observation slot of your merged data.

10X单细胞（10X空间转录组）降维分析之UMAP

参考技术A

UMAP ，全称uniform manifold approximation and projection，统一流形逼近与投影，是基于黎曼几何和代数拓扑的理论框架结构构建的。在处理大数据集时，UMAP优势明显，运行速度更快，内存占用小。Etienne Becht等人2019年在Nature Biotechnology上发表一篇文章将其应用在生物学数据上并阐述了UMAP在处理单细胞数据方面的应用和优势。

如果你不知道tSNE是什么，它是如何工作的，也没有读过2008年的革命性的van der Maaten & Hinton原稿，可以参考我的那文章 10X单细胞（10X空间转录组）降维分析之tSNE（算法基础知识） 。尽管tSNE对一般的单细胞基因组学和数据科学产生了巨大的影响，但人们普遍认为它有一些缺点，这些缺点很快将得到解决。（ tSNE的缺点在上次分享的文章中也做过详细的介绍 ）。

看看上面的图，我想说的是 t分布应该提供全局距离信息，因为它们将高维空间中相距较远的点推到低维空间中更远的点。

然而，这种良好的意愿被成本函数(KL-divergence)的选择所扼杀，我们将在后面看到其原因。

(1),可以显著降低计算时间高维图像由于求和或集成是一个代价高昂的计算过程。想想马尔可夫链蒙特卡罗(MCMC)它基本上是试图近似地计算在贝叶斯规则的分母上的积分（UMAP使用最近邻居的数量而不是perplexity）

(2)定义perplexity, UMAP则定义了没有log2函数的最近邻居k的个数，即:

UMAP使用稍微不同的高维概率对称

symmterization是必要的因为UMAP融合在一起的点与本地不同的指标(通过参数ρ),它可能发生图A和B节点之间的重量不等于B之间的权重和节点。为什么UMAP使用这种对称而不是tSNE使用的对称还不清楚。我将在下一篇文章(从头开始编写UMAP)中展示我对不同的对称化规则的实验，这并没有使我相信这是如此重要的一步，因为它对最终的低维嵌入式产生了很小的影响。
UMAP使用曲线族1 / (1+a*y^(2b))在低维中建模距离概率，不是完全的学生t分布，但非常非常相似，请注意再次没有应用标准化:

其中，对于默认UMAP超参数a≈1.93,b≈0.79(实际上，对于min_dist = 0.001)。在实践中，UMAP从非线性最小二乘拟合到带有min_dist超参数的分段函数中找到a和b:

为了更好地理解曲线族1 / (1+a*y^(2b))的行为，让我们画出不同a和b的曲线:

我们可以看到曲线族对参数b非常敏感，在大的参数b处，在小的参数y处，曲线族形成了一种高峰。这意味着在UMAP超参数min_dist之下，所有的数据点都是同样紧密相连的。由于Q(Y)函数的行为几乎像一个Heaviside阶跃函数，这意味着UMAP为所有在低维空间中相互靠近的点分配了几乎相同的低维坐标。min_dist正是导致在UMAP维数降低图中经常观察到的超密集集群的原因。

为了演示如何准确地找到a和b参数，让我们展示一个简单的分段函数(其中高峰部分是通过min_dist参数定义的)，并使用函数族1 / (1+a y^(2b))通过优化来拟合它。curve_fit来自Scipy Python库。作为拟合的结果，我们得到了函数1 / (1+a y^(2b))的初值a和初值b。

由于我们需要知道交叉熵的梯度，以便以后实现梯度下降，让我们快速计算它。忽略只包含p(X)的常数项，我们可以将交叉熵重新写一下，并将其微分如下:

图拉普拉斯、谱聚类、拉普拉斯Eignemaps、扩散图、谱嵌入等，实际上是指将矩阵分解和邻接图方法结合起来解决降维问题的同一种有趣的方法。在这种方法中，我们首先构造一个图(或knn图)，然后通过构造拉普拉斯矩阵用矩阵代数(邻接矩阵和度矩阵)将其形式化，最后分解拉普拉斯矩阵，即求解特征值分解问题。

我们可以使用scikit-learn Python库，并使用spectralembedded函数在演示数据集(即与癌症相关的成纤维细胞(CAFs) scRNAseq数据)上轻松地显示初始的低维坐标:

最后，UMAP使用随机梯度下降(SGD)代替常规梯度下降(GD)，如tSNE / FItSNE，这既加快了计算速度，又减少了内存消耗。

现在让我们简要地讨论一下为什么他们说tSNE只保留数据的局部结构。可以从不同的角度来理解tSNE的局部性。首先,我们有σ参数Eq。(1)本地数据点集这样互相“感觉”。因为成对欧几里得距离衰减指数的概率,在小的σ值,它基本上是零遥远的点(大型X)和快速增长仅为最近的邻居(小X)。相比之下,在大的σ,遥远而近点的概率成为限制可比和σ→∞,概率就等于1为所有任何一对点之间的距离,即成为等距点。

有趣的是,如果我们扩大成对欧几里得距离的概率高维度成泰勒级数在σ→∞,我们会在第二近似幂律:

关于两两欧几里得距离的幂律类似于多维定标法(MDS)的成本函数，MDS是通过保存每对点之间的距离来保存全局距离，而不管它们是相距很远还是很近。一个可以解释这个大的σtSNE远程数据点之间的相互作用,所以是不完全正确的说tSNE只能处理当地的距离。然而，我们通常会受到perplexity有限值的限制，Laurens van der Maaten建议perplexity的取值范围在5到50之间，尽管在局部信息和全局信息之间可能会有一个很好的折衷，那就是使用平方根≈N^(1/2)来选择perplexity，其中N为样本量。相反的极限,σ→0,我们最终的极端“局部性”高维概率的行为类似于狄拉克δ函数的行为。

另一种理解tSNE“局部性”的方法是考虑KL-divergence函数。假设X是高维空间中点之间的距离Y是低维空间中点之间的距离

根据kl -散度的定义：

方程(9)的第一项对于X的大小都是趋近于0的，对于X的大小也是趋近于0的，因为指数趋近于1，而log(1)=0。对于大X，这一项仍然趋近于0因为指数前因子趋近于0的速度快于对数趋近于负无穷。因此，为了直观地理解kl散度，只考虑第二项就足够了:

这是一个看起来很奇怪的函数，让我们画出KL(X, Y)

这个函数的形状非常不对称。如果点在高维度X之间的距离很小,指数因子变成1和对数项行为日志(1 + Y ^ 2)这意味着如果Y是在低维空间的距离大,将会有一个大的惩罚,因此tSNE试图减少Y在小X为了减少罚款。相反，对于高维空间中的长距离X, Y基本上可以是0到∞之间的任何值，因为指数项趋于0，并且总是胜过对数项。因此，在高维空间中相距遥远的点，在低维空间中可能会相互靠近。因此，换句话说，tSNE并不能保证高维空间中相距较远的点在低维空间中会保持较远的距离。然而，它确实保证了在高维空间中相邻的点在低维空间中保持相邻。所以tSNE不是很擅长远距离投射至低维,所以它只保留本地数据结构提供了σ不去∞。

与tSNE不同，UMAP使用交叉熵(CE)作为成本函数，而不是KL-divergence

这导致了地方-全球结构保护平衡的巨大变化。在X的小值处，我们得到了与tSNE相同的极限，因为第二项由于前因子和对数函数比多项式函数慢的事实而消失:

因此，为了使惩罚规则最小化，Y坐标必须非常小，即Y→0。这与tSNE的行为完全一样。但是，在大X的相反极限，即X→∞时，第一项消失，第二项的前因子为1，得到:

这里，如果Y很小，我们会得到一个很大的惩罚，因为Y在对数的分母上，因此，我们鼓励Y很大，这样，对数下的比率就变成了1，我们得到零惩罚。因此，我们在X→∞处得到Y→∞，所以从高维空间到低维空间的整体距离保持不变，这正是我们想要的。为了说明这一点，让我们绘制UMAP CE成本函数:

在这里，我们可以看到图的“右”部分看起来与上面的kl散度曲面非常相似。这意味着在X低的时候，为了减少损失，我们仍然想要Y低。然而，当X很大时，Y的距离也要很大，因为如果它很小，CE (X, Y)的损失将是巨大的。记住，之前，对于KL (X, Y)曲面，在X很大的情况下，我们在高Y值和低Y值之间没有差别，这就是为什么CE (X, Y)代价函数能够保持全局距离和局部距离。

我们知道UMAP是速度比tSNE担忧)时大量的数据点,b)嵌入维数大于2或3,c)大量环境维度的数据集。在这里,让我们试着了解UMAP要优于tSNE来自于数学和算法实现。

tSNE和UMAP基本上都包含两个步骤：

然而，我注意到UMAP的第一步比tSNE快得多。这有两个原因：

接下来，UMAP实际上在第二步中也变得更快了。这种改善也有几个原因:

在这篇文章中，我们了解到尽管tSNE多年来一直服务于单细胞研究领域，但它有太多的缺点，如速度快、缺乏全球距离保存。UMAP总体上遵循了tSNE的哲学，但是引入了一些改进，例如另一个成本函数和缺少高维和低维概率的标准化。

除了运行速度快，内存占用小等特点，UMAP在处理细胞学数据时还有一个大的优势，就是可以反映细胞群体之间分化的连续性和组织性。下面将通过文献中的数据【2】来为大家详细讲解。

对同一组数据分别进行tSNE和UMAP降维，该数据为多达30万个从8种不同组织富集得到的T细胞和NK细胞的样本，并使用Phenograph聚类把细胞分为6大类，每种颜色代表一种细胞。从图中可以看出，UMAP和tSNE都可以较好地把不同类别的细胞分开。但tSNE倾向于把相同细胞群划分为更多的群，如图显示，黑色圈中CD8 T细胞，在tSNE结果中，群数更多，距离更远。

同样这组数据用组织来源对UMAP和t-SNE图上细胞的进行颜色区分，可以观察到一个有意思的现象。与UMAP相比，t-SNE更加倾向于根据它们的来源来分离总体细胞。而 UMAP则会兼顾细胞群的类别和来源来排列，如图中在CD4 T细胞和CD8 T细胞群内，细胞的排列与来源也会有一定的规律性，都是大致从脐带血（CB）和外周血单核细胞（PBMC），到肝脏（Liver）和脾脏（Spleen），最后到一端的扁桃或另一端的皮肤（Skin）、肠道（Gut）和肺（Lung）。

通过驻留记忆T细胞标志物CD69/CD103、记忆T细胞标志物CD45 RO和naïve T细胞标志物CCR7表达群的分布，可以观察到UMAP可以展示出T细胞连续的分化阶段。而tSNE结果中，这些群之间也是连续的，但是却没有非常明显的沿轴结构。同样的现象也在造血细胞系统中被观察到。由此可见， UMAP在大数据集的处理时可以展现细胞集群的连续性。

对三组数据（Samusik、Wong、Han_400k）分别进行数据随机降低至100-200,000之间不同的数量级，形成小数据集。纵轴为小数据集与原始数据集的相关性，代表降维方法在不同数据量上的可重复性。UMAP表现最好，数据集越大，优势越明显。

下图是UMAP和t-SNE对一套784维Fashion MNIST高维数据集降维到3维的效果的比较。

虽然这两种算法都表现出强大的局部聚类并将相似的类别分组在一起，但UMAP还将这些相似类别的分组彼此分开。另外，UMAP降维用了4分钟，而多核t-SNE用了27分钟。

UMAP的两个最常用的参数：n_neighbors 和 min_dist，它们可有效地用于控制最终结果中局部结构和全局结构之间的平衡。

最重要的参数是 n_neighbors ，近似最近邻居数。它有效地控制了UMAP局部结构与全局结构的平衡，数据较小时，UMAP会更加关注局部结构，数据较大时，UMAP会趋向于代表大图结构，丢掉一些细节。

第二个参数是 min_dist，点之间的最小距离。此参数控制UMAP聚集在一起的紧密程度，数据较小时，会更紧密。较大的值会更松散，而将重点放在保留广泛的拓扑结构上。

t-SNE和UMAP大部分的表现非常相似，但以下示例明显例外：宽而稀疏的cluster中有密集的cluster(如下图所示)。UMAP无法分离两个嵌套的群集，尤其是在维数较高时。

UMAP在初始图形构造中局部距离的使用可以解释该算法无法处理情况的原因。由于高维点之间的距离趋于非常相似(维数的诅咒)，所以可能会因此将其混合在一起。

算法很难，所以懂的人才显得牛

天行健，君子以自强不息