单细胞转录组测序分析--初探Seurat

Posted 2023-05-06

参考技术A 时代发展的步伐总是毫不留情的将你甩在身后，连车尾灯都看不见。当你还在沉迷于普通转录组数据挖掘时，已经有人悄悄的搞上单细胞了。单细胞转录组测序，顾名思义就是在单个细胞的分辨率基础上去研究细胞内的基因表达等，其主要目的是为了研究不同细胞类型的基因表达异质性，从而解决相关生物学问题。谈到单细胞就不得不提一下当下火爆的10x Genomics服务商了，具体参见 10x Genomics 。本篇文章暂时不介绍10x，主要介绍单细胞转录组数据分析软件Seurat。
Seurat软件是一个R包，可以说是单细胞转录组测序分析的明星软件，很多单细胞测序文章都会引用该软件，引用次数也是杠杠的，而且也有详细的 在线教程 。本文也主要是根据其教程介绍一下使用Seurat软件分析一个样本的单细胞转录组数据的步骤及注意事项，供大家讨论。

导入分析需要的包

Seurat软件提供了很友好的函数可以直接读取10x Genomics的输出结果

导入文件后便可以创建Seurat对象

创建完Seurat对象后，Seurat将数据保存在不同的slot中，如[e-mail], filter_10x_objectt@data, [e-mail], filter_10x_object@ident，其中raw.data存放的是每个细胞中每个gene的原始UMI数据，data存放的是gene的表达量，meta.data存放的是每个细胞的统计数据如UMI数目，gene数目等，ident此时存放的是project信息。

由于技术原因，一个GEM中可能会包含2个或多个细胞，也可能不包含细胞，这时候可以通过观察每个barcode中的基因数目或UMI数目来判断。

上图展示的是每个barcode中的基因数目和UMI数目的关系，一般二者都成正相关关系，有个别barcode的基因数目和UMI数目过高，有可能就是包含2个细胞的GEM，可以考虑在后续分析中将其过滤掉。
我们不仅仅可以观察每个barcode的基因数目，还可以计算每个barcode中的线粒体基因含量等，从而更加仔细的观察数据的质量。

这张图片展示了每个barcode中基因数目、UMI数目以及线粒体基因含量的分布情况，根据上述2张图片就可以大致确定是否需要过滤哪些数据进行后续分析。
Seurat提供了一个很好用的数据过滤函数：

以上就是数据的预处理过程了，接下来就进入正式的分析阶段，包括数据的标准化、归一化、数据降维以及聚类分析等。

FindVariableGenes算法：首先计算基因的平均表达量，然后计算基因的离散度；接下来根据平均表达值将基因分成20块并计算每块的离散度的Z值。
如上图：横坐标代表基因的平均表达量，纵坐标代表基因的离散度的Z值，标有基因名的点就是由函数中的cutoff值决定的，改变cutoff值，这些标记也会随之改变。
数据的线性回归、中心化和比例化：对数据进行线性回归分析，去除不想要的变异源。
中心化：首先计算基因A在所有细胞中的平均表达量，然后分别将每个细胞中基因A的表达值减去平均值。
比例化：在中心化的基础上，首先计算基因A在所有细胞中的中心化值后的标准差，然后分别将每个细胞中基因A的中心化值除以标准差。这些步骤都在一个函数中完成。

单细胞转录组测序产生的数据是数万个基因在数万个细胞中的表达情况，属于典型的高维数据。如果把1个基因视为1个坐标轴的话，那么一个细胞的空间位置就是在数万个坐标轴中的定位，这样的话相同细胞类型的细胞就应该挨在一起，我们就可以根据细胞的空间位置判断细胞亚群了。可是我最多也就认识三维坐标啊，咋办，能不能把这些高维数据投影到二维坐标呢，那就交给PCA和t-SNE吧。PCA和tSNE都是数据降维分析方法，PCA属于线性降维，tSNE属于非线性降维。我们先执行PCA分析，使高维数据的信息最大程度保留在低维数据中，PCA分析利用的是保存在scale.data的值。

执行完PCA分析后，就要根据PCA得分来进行聚类分析了，但是在进行聚类分析之前，需要选择使用对少个主成分进行计算。每个主成分实际上代表的是相关基因集的信息，因此确定多少个主成分是一个重要的步骤，我们可以根据PCElbowPlot函数来判断。

从上图可以看到，拐点出现在10-15之间，我们可以选择15来进行聚类分析。Seurat采用的是基于图形的聚类方法，即利用PCA空间中的欧几里德距离构造一个KNN图（数学好的可以留下来帮忙讲讲）。

好了，到此我们就知道了我们的数据中有多少种细胞亚群了，怎么可以少得了图片展示呢。超棒的可视化方法tSNE要上场了。tSNE的目标是将在高维空间中具有相似局部邻域的细胞，在低维空间中放在一起。

既然我们知道了有多少种细胞亚群，那么是不是就要分析一下这些亚群间的差异性呢，交给FindAllMarkers吧。FindAllMarkers能够同时计算所有亚群的差异性（分别计算每个亚群与剩下的所有细胞的差异性）。

得到差异表达基因后，当然要进行展示了。

好了，剩下的就是进行生物学知识挖掘了，例如根据这些差异基因推断细胞类型啊之类的。
关于单个样本的单细胞转录组数据分析就介绍到这儿了，那多个样本的分析会有什么不同呢，我们下次再说吧。

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

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

天行健，君子以自强不息