Seurat识别细胞类群的原理（FindNeighbors和FindClusters）

Posted 2023-05-06

参考技术A

众所周知，seurat在降维之后主要依据两个函数来进行细胞分类，这里我们来深入了解一下seurat如何进行细胞分类的。
首先我们来看有关分类的两个函数

我们来一一解决其中的问题

Jaccard index , 又称为Jaccard相似系数（Jaccard similarity coefficient）用于比较有限样本集之间的相似性与差异性。Jaccard系数值越大，样本相似度越高。
给定两个 集合 A,B，Jaccard 系数定义为A与B交集的大小与A与B并集的大小的比值，定义如下：

首先来说annoy ，annoy全称“Approximate Nearest Neighbors Oh Yeah”，是一种适合实际应用的快速相似查找算法。Annoy 同样通过建立一个二叉树来使得每个点查找时间复杂度是O(log n)，和kd树不同的是，annoy没有对k维特征进行切分。annoy的每一次空间划分，可以看作聚类数为2的KMeans过程。收敛后在产生的两个聚类中心连线之间建立一条垂线（图中的黑线），把数据空间划分为两部分。在划分的子空间内不停的递归迭代继续划分，直到每个子空间最多只剩下K个数据节点，划分结束。最终生成的二叉树具有如下类似结构，二叉树底层是叶子节点记录原始数据节点，其他中间节点记录的是分割超平面的信息。

接下来查看第二个函数

里面有很多算法方面的内容，对于数学而言， 不会就是不会 ，没什么炫耀的,但是要知道其原理。
，这里分享给大家一些链接，大家可以多多研究研究。
Kronecker delta函数
leiden
Louvain
SLM
请保持愤怒，让王多鱼倾家荡产

单细胞数据处理小细节汇总

参考技术A

1. Seurat对象查看当前的Assay

在进行了SCTransform操作后，矩阵默认会变成SCT矩阵，如果不加设置，后续的PCA等操作都是基于SCT矩阵。
修改DefaultAssay：

2. Seurat使用FindVariables找到高变基因后增删高变基因

3. 不同运行步骤中的文件来源和储存位置⚠️

⚠️：PCA的值是可以被覆盖的，使用三步法对矩阵进行标准化后进行PCA后再使用SCT矩阵进行标准化，PCA的矩阵变成了SCT的PCA矩阵，原有的PCA矩阵不会保留。后续的TSNE和UMAP降维图也和三步法不一样。

4. @data标准化矩阵 和 @scale.data 归一化矩阵 的区别
单细胞RNA 测序数据中，文库之间测序覆盖率的系统差异通常是由细胞间的cDNA 捕获或PCR 扩增效率方面的技术差异引起的，这归因于用最少的起始材料难以实现一致的文库制备。标准化旨在消除这些差异（例如长度，GC 含量），以使它们不干扰细胞之间表达谱的比较。这样可以确保在细胞群体中观察到的任何异质性或差异表达都是由生物学而不是技术偏倚引起的（批次矫正仅在批次之间发生，并且必须同时考虑技术偏差和生物学差异，标准化只需考虑技术差异）。
软件Seurat 提供了三种标准化的方法，分别为LogNormalize、CLR、RC，通常情况下我们采用LogNormalize 的方式进行标准化，计算公式为：log1p((Feature counts/total counts) ∗ scale. factor)

归一化的目的则是使特征具有相同的度量尺度

参考： Seurat的normalization和scaling

5. 关于有些细胞属于同一个cluster但是在umap或者tsne图上相聚较远的问题：UMAP和TSNE是各自的算法在PCA降维的基础上再进行非线形降维，在二维图上把其各自算法认为相近的细胞聚在一起。但是FindClusters输入的不是UMAP或TSNE降维的数据，而是FindNeighbors的数据，而FindNeighbors输入的数据是PCA降维数据，是用另外一种算法计算的细胞之间的距离。因此会出现有些细胞被认为是同一个cluster，但是在umap或者tsne图上相聚较远（尤其是一些散在的，脱离主群的细胞）

6. marker基因鉴定，查看marker基因的表达是使用RNA矩阵还是sct矩阵？
这是一个争议性问题，两个都可以，目前建议最好使用RNA矩阵。
sct的到的count并不是真实的基因表达值，而是通过scaledata倒推出来的，它是一个回归，运算之后的残差。

7. 关于FindAllMarks找到的基因
如下图，先看cluster0的Marker基因：cluster0的差异基因是cluster0的细胞和剩下的所有的cluster合在一起的细胞做对比得到的。pct.1是这个基因在cluster0中的表达比例（S100A8在cluster0的细胞中的表达比例是100%），pct.2是这个基因在除了cluster0以外的所有细胞中的表达比例（S100A8在除cluster0以外的细胞中的表达比例是51.2%）。avg_log2FC是表达差异倍数，p_val_adj是校正后的p值。

8. 在提取Marker基因时比较好的办法：因为单细胞矩阵算出来的结果，p_val_adj==0的有很多，所以可以先把p_val_adj==0先提出来。再把p_val_adj<0.01的按差异倍数靠前的20/30/500...（按需要）个基因提出来，然后把两个矩阵合在一起（取交集）用来做细胞鉴定。（结合p值和fc来做筛选效果更好）

⚠️：提取没有核糖体和线粒体的marker基因更好。（这些基因对鉴定没有帮助）

有些基因比如Foxp3，对细胞鉴定很重要，但常常在筛选Marker基因的时候筛选不出来。 非负矩阵分解 可能更好。
参考： 过滤线粒体核糖体基因

9. 提取亚群

⚠️ 新提取的亚群需要重新进行降维聚类 （和大群相比，标记基因发生了变化），并重新寻找marker基因，重新分群，注释。❗️subset提取子集后，不同样本间批次校正的信息也被去除，需要重新进行批次矫正

参考： Seurat取子集时会用到的函数和方法 ⚠️⚠️⚠️

10. 取子集后如何去除空子集（还存在这个level，但里面包含的细胞为0，如何去除） as.factor(as.character())

11. 双细胞的预测和去除如DoubletFinder建议单样本进行，不建议双样本一起预测。除此之外，其他步骤都可以多样本一起做，质控的时候也可以多样本一起做，但是建议每个样本都单独看一看。

12. 单细胞多样本整合：merge()；多样本拆分：SplitObject()

13. 在做多组数据整合，每个组又有多个样本的时候，最好把单独的每个样本当成批次，而不是把不同的组当成批次。

14. 多核运算

参考： 单细胞数据分析中future包的使用

15. pbmc3k.final@commands$FindClusters 可以查看FindClusters运行时间和记录。Seurat是记录其分析过程的，也可查看command下其他操作

16. 关于质控标准：同一组织的最好用同样的标准，不同组织的可以不一样。不同组织线粒体含量等可能存在差异。

17. 可视化的方法总结
参考： https://www.jianshu.com/p/0d1e2c7d21a4

18. circos图绘制

19. 单细胞数据思维导图，有利于查看单细胞数据格式。
https://www.jianshu.com/p/7560f4fd0d77

20. 对于旧版本Seurat对象的更新
scRNA <- UpdateSeuratObject(scRNA)
UpdateSeuratObject SeuratObject ：Update old Seurat object to accommodate new features

21. 对有些操作需要用到python设置的情况

22. 单细胞数据做pooling的好处：可以尽量的降低dropout的问题。（dropout就是矩阵中的zero，这些zero实际上并不是0，而是每个液滴里面起始反应量太低了。而一般的反转录效率只能到30%左右，70%的转录本实际上在反转录那一步是被丢掉的，这是单细胞测序一个比较大的问题）。
但是一旦做了pooling，你必须要证明pooling对结果是没有影响的（下图的右面三个图）。

23. Seurat的VlnPlot中的combine参数，在如下画三个基因的情况下，设置成T就画一张图，设置成False，会将三个基因各画一张图。

24. rev()这一步是将横坐标的基因反过来排序

这两个画出来的图横坐标基因的顺序是相反的（见NicheNet）

25. 堆叠小提琴图的绘制
完成这个需求有以下几种实现方法：
1. Seurat包直接就可以实现（stack = T）
2. 通过Seuart->scanpy来实现，第一张是Seurat包VlnPlot函数画的图，第二张是scanpy中stacked_violin函数画的图，那么现在问题就变成为Seurat对象到scanpy对象的转换
3. 用R原生函数实现StackedVlnPlot的方法
4. 使用基于scanpy包衍生的scanyuan包来说实现
5. 使用R包MySeuratWrappers来实现

最简单的方法1如下：

如果不设置level，会按字母顺序排列，case会自动排在con前面。

使用Seurat的 RenameIdents 函数也可以

x: table
margin: a vector giving the margins to split by. E.g., for a matrix 1 indicates rows, 2 indicates columns, c(1, 2) indicates rows and columns. When x has named dimnames, it can be a character vector selecting dimension names.

得到的HC_1样本的orig.ident默认是样本名中第一个_号的前一部分。所以要保证矩阵的列名是 样本名_细胞barcode 这样的格式。
如果有多个分组，例如两个样本矩阵中细胞分别命名为HC_1_barcode，HC_2_barcode，在直接通过如下方法得到两个Seurat对象，再对其进行merge之后，两个样本会被合并成一个。也就是样本信息只保留了第一个_号之前的HC，没有保留_号之后的1和2。

为了避免这种情况，可以在构建Seurat对象时通过参数进行设置

⚠️PC数的选择：Seurat官网提供的三种方法只能给出PC数的粗略范围，选择不同PC数目，细胞聚类效果差别较大，因此，需要一个更具体的PC数目。作者提出一个确定PC阈值的三个标准：

一般先选默认分辨率(0.8)，大概可能会分出十几个群。因为最终都是要注释到每一个barcode，所以首先可以看大类marker的分布(不受分辨率影响)，可以根据marker基因的分布来调整分辨率。是否需要精细的分群得看精细的分群对研究有没有决定作用，还有很重要的一点是 看分出的各个cluster在Findallmarkers给出的结果中marker的热图是不是能明显分开 。精细划分的细胞本来就很类似，如果有部分小群的热图明显分不开或者非常类似，就可以考虑把分辨率调小。

这实际上是没有必要的必须保持一致的。下游的都是用pca之后的，pca是为了压缩数据。
umap和tsne是为了可视化(仅仅是可视化)，但是FindNeighbor是计算细胞间距离矩阵。找类群数目和可视化可以说没有关系。

map函数：
R语言循环第三境界：purrr包map函数！
浅析R语言中map(映射)与reduce(规约)

参考： monocle2

查看不同细胞群的中位基因也是一样

查看不同样品的中位基因也是一样

或者也可以