seurat涉及的数据分析包括很多步骤。
之前只顾着干活儿，也没有系统的整理过分析中的具体内容。
这里就参照网上大神们分享的帖子，来梳理一下。
1、读入数据。
这个函数就是根据输入矩阵/数据框，创建Seurat对象的。重要步骤是 设置 ident 和添加 meta.data。
dat <- Read10X(data.dir = "your/work/path")
organoids <- CreateSeuratObject(counts = dat, project = "Organoids", min.cells = 3, min.features = 200)
*min.cells 表示一个基因至少要在3个细胞中被检测到，否则不要。
*min.features 参数指定每个细胞需要检测的最小基因数量。此参数将过滤掉质量较差的细胞，这些细胞可能只是封装了随机barcodes，而没有任何真实的细胞。通常，检测到的基因少于100或者200个的细胞不会被考虑进行分析。
这里还是设计一个知识点就是R里面的S3类和S4类。
list一般情况下被认为是S3类，S4类是指使用slots存储数据的格式。（如果说的不对欢迎中纠错）
这里读进去的数值是三个文本文件创建的稀疏矩阵。
什么是稀疏矩阵？
在[矩阵]中，若数值为0的元素数目远远多于非0元素的数目，并且非0元素分布没有规律时，则称该矩阵为稀疏矩阵；与之相反，若非0元素数目占大多数时，则称该矩阵为稠密矩阵。
如果自己手上有单细胞数据，那个matrix文件里面包含很多0。
2、NormalizeData
因为在测序之前会对抓取到的RNA进行PCR扩增，所以需要考虑文库深度的对测序的影响，所以需要对上一步得到的稀疏矩阵进行Normalize。
Normalize的方式：每个细胞每个基因的特征计数除以该细胞的特征总计数，再乘以scale.factor（默认10000），然后使用log1p进行对数转换。（log1p=log(n+1)）
Normalize之后的数据储存在seurat[['RNA']]@data这里。
combined.data <- NormalizeData(combined.data, normalization.method = "LogNormalize", scale.factor = 10000)
2.1、单个样本Normalize和多个样本合并之后Normalize的结果是否存在差异。
首先我们合并数据的时候一般直接用的是merge，所以不同样本的细胞的数值不会发生变化。
2.2、scale factor为什么是10000？
这里截取whitebird的解释。
scale factor--from whitebird
2.3、为什么要进行log转化？
做过传统转录组分析的家人们都明白，用转录组数据的FPKM和TPM绘制热图等等的时候，因为数值的变化范围太过巨大，都需要进行一个log转换，让数据压缩在一个区间。
其次，也是最重要的改变数据分布：测序数值本身不符合正态分布，log转换能让数据趋近于正态分布，方便后续的进一步分析。
3、FindVariableFeatures()
高变异基因就是highly variable features（HVGs），就是在细胞与细胞间进行比较，选择表达量差别最大的基因，Seurat使用FindVariableFeatures函数鉴定高可变基因，这些基因在不同细胞之间的表达量差异很大（在一些细胞中高表达，在另一些细胞中低表达）。
默认情况下，会返回2,000个高可变基因用于下游的分析，如PCA等。
3.1、vst的基本流程
算法实现在 FindVariableFeatures.default() 中。
目的是在var~mean曲线中，不同mean值区域都能挑选var较大的基因。
- 使用loess（局部加权回归）拟合平滑曲线模型
- 获取模型计算的值作为y=var.exp值
- var.standarlized = get variance after feature standardization: (每个基因 - mean)/sd 后 取var(). 注意sd=sqrt(var.exp)
- 按照 var.standarlized 降序排列，取前n(比如2000)个基因作为高变基因。
combined.data <- FindVariableFeatures(combined.data, selection.method = "vst", nfeatures = 2000)
4、ScaleData（）归一化
单细胞基因表达counts矩阵数据经过NormalizeData()处理后，还需要进行scale。
[ScaleData()]函数将基因的表达转换为Z分数（值以 0 为中心，方差为 1）。
它存储在 seurat_obj[['RNA']]@scale.data，用于下游的PCA降维。
默认是仅在高可变基因上运行标准化。
combined.data <- ScaleData(combined.data)
最开始分析单细胞的时候，这里有点疑惑。
为什么前面已经Normalize，这里还要scale一下？
"Scaling与Normalization的区别"
scale改变的是数据的范围，normalize改变的是数据的分布。
scale是将数据的分布限定在一个范围内，这样子方便比较。normalize却是将偏态分布转换成趋近于正态分布。
这里引用简书“whitebird”所写的内容。
R语言的Z score计算是通过[scale()]函数求得，Seurat包中ScaleData()函数也基本参照了scale()函数的功能。
scale方法中的两个参数：center和scale
- center和scale默认为真，即T或者TRUE
- center为真表示数据中心化：数据集中的各项数据减去数据集的均值
如有数据集1, 2, 3, 6, 3，其均值为3
那么中心化之后的数据集为1-3,2-3,3-3,6-3,3-3，即：-2,-1,0,3,0 - scale为真表示数据标准化：指中心化之后的数据在除以数据集的标准差，即数据集中的各项数据减去数据集的均值再除以数据集的标准差。
如有数据集1, 2, 3, 6, 3，其均值为3，其标准差为1.87
那么标准化之后的数据集为(1-3)/1.87,(2-3)/1.87,(3-3)/1.87,(6-3)/1.87,(3-3)/1.87，即：-1.069,-0.535,0,1.604,0
Z score的概念是指原始数据距离均值有多少个标准差。当以标准差为单位进行测量时，Z score衡量的是一个数值偏离总体均值以上或以下多少个标准差。如果原始数值高于均值，则 Z score得分为正，如果低于均值，则Z score为负。
Z score其实是标准正态分布（Standard Normal Distribution），即平均值μ=0，标准差 σ=1 的正态分布。SND标准正态分布的直方图如下所示：
Z score分布--from whitebird
5、RunPCA()降维处理
Seurat使用RunPCA函数对标准化后的表达矩阵进行PCA降维处理。
默认情况下，只对前面选出的2000个高可变基因进行线性降维，也可以通过feature参数指定想要降维的数据集。
RunPCA之后会返回5个PC和对应的一大堆基因。
每个PC对应60个基因，分为Positive和Negative，各30个。
Positive和Negative就是PC轴的正负映射关系，正值为Positive，负值为Negative。
返回的是正值和负值绝对值最大的top30，可以理解为对所有细胞区分度最大的基因。
combined.data <- RunPCA(combined.data)
View(pbmc@reductions[["pca"]]@feature.loadings)
View(pbmc@reductions[["pca"]]@cell.embeddings)
这是PCA之后得到的两个结果。
第一个是每个PC对应的基因。
第二个是每个细胞对应的PC上的坐标。
得到的PCA的结果还可以用以下两种方式来看。
DimPlot(pbmc,reduction = 'pca')
p1=FeaturePlot(pbmc,features = "PC_1", order = T)
p2=FeaturePlot(pbmc,features = "PC_2", order = T)
p1|p2
6、FindNeighbors()
首先计算每个细胞的KNN，也就是计算每个细胞之间的相互距离，依据细胞之间邻居的overlap来构建snn graph。
计算给定数据集的k.param最近邻。也可以选择(通过compute.SNN)，通过计算每个细胞最近邻之间的邻域重叠(Jaccard索引)和其邻近的k.param来构造SNN。
combined.data <- FindNeighbors(combined.data, dims = 1:30)
具体在计算细胞之间的距离的时候呢，用得到的KNN算法，即邻近算法。
但是，这个算法我也不太懂，但是其中有个事情还蛮有趣。
就是判定邻近的模块有两个：annoy、rann
其中这个annoy全称又叫做Approximate Nearest Neighbors Oh Yeah，名字还蛮可爱的。
下面这位博主对于这个算法有非常详细的解释，有兴趣的家人们自行前往观看。
FindNeighbors {Seurat} - 简书 (jianshu.com)
7、FindClusters()
就是在已经计算完细胞之间的距离之后，对这些细胞进行分类。
可以指定分为几类细胞。
但是很多参考资料里面最重要强调的都只是一个参数：resolution。
resolution这个参数设置的大小决定了细胞类型的多少，值越大细胞类型越多。
具体分析的时候很多时候就会问：到底多少合理？到底应该分为几群？
其实对于测序的人来讲，很多时候确实也不太清楚到底有多少种细胞。
那就有两种解决办法。
1）直接看tSNE的图，物理距离就是判断的一种方法。当物理距离很近的一群细胞被拆开了，那就说明可能没拆开之前是合理的。但是，这种方法呢就简单粗暴一些。
2）有另外一个包clustree，可以对你的分群数据进行判断。
如下图中，当分为4、5和6群的时候，细胞之间没有太多的交叉，都是在进一步细分。
但是再往下，那些半透明的箭头就表示，从一个细胞群跑到了另外一个细胞群，说明就不太合适。
clustree
所以，参照上述两种判断方法，可以得到结果。
8、RunUMAP()、RunTSNE()
上述两种方法其实都是对数据进行降维。
为什么需要降维？
如果将单个细胞看作一个数据点，那么检测的基因数就是其对应的变量数，也就是我们所说的维数。一般一个人类scrna-seq 能检测到～25k基因，于是基因表达数据也被称为高维数据。我们知道， 并不是所有的基因都会表达，并且由于单细胞测序技术的限制，不是所有的转录分子都能被成功捕获，再加上测序深度的差异， 每个细胞中约能检测到10%～50%的转录分子，这导致了许多基因计数为0. 虽然通过数据的预处理我们能过滤掉零计数基因，剩下的基因维数仍然可以高达～15k。并且基因之间的高度相关性也使得表达数据包含了许多冗余信息， 从而掩盖了真正的生物学差异。因此，我们需要通过降维抽取数据的概要（component），减少数据噪音，并提高下游分析速度。降维的另外一个好处就是可以在低维中实现数据的可视化。
常见的降维方法有PCA (principle component analysis)、MDS (multi-dimensional scaling)、Sammon mapping、Isomap 、t-SNE （t-distributed stochastic neighbor embedding）和UMAP（Uniform Approximation and Projection method）。
其中PCA， t-SNE和UMAP在scRNA-seq中使用非常普遍。
之前的分析中已经用PCA降维了，为什么这里还要降维？
PCA输入的是@scale.data，是属于线性降维
RunUMAP()和RunTSNE()输入的数据是PCA之后的数据（PCA@cell.embedding），属于非线性降维。输出的结果保存在各自的slots里面。
这两种方法有什么区别？
摘抄自：（跟着小鱼头学单细胞测序-scRNA-seq数据的降维和可视化 - 云+社区 - 腾讯云 (tencent.com)
）
t-SNE和UMAP是另外两种非线性的降维方法，由于其漂亮的可视化效果，这两种方法在单细胞数据教程中非常受欢迎。与PCA不同，这两种方法不仅限于线性变换，也不受制于准确表示远距离种群之间的距离，因此使得它们在低位空间如何对细胞进行排列具有更大的自由度，从而在其二维的可视化效果图中能够将复杂的细胞群分成许多不同的簇，使得效果图比PCA更加直观和容易解释。
t-SNE是以牺牲整体结构（global structure）为代价，着重于捕获局部结构（local similarly），因此它可能会夸大细胞群体之间的差异而忽略群体之间的潜在联系。如果简单的通过t-SNE结果图来解读细胞簇之间的关系，产生的结论可能会被图中的集群大小和位置所误导。在效果图中，t-SNE倾向于将密集的簇膨胀，并且压缩稀疏的簇，因此我们不能简单的通过图上的大小来衡量细胞群的差异。并且由于t-sne无法保证能保留距离较远的簇的相对位置，我们也不能简单通过图中的位置来确定远距离细胞簇之间的关系。相比之下，UMAP与t-SNE类似，同时在低维空间保留了高维空间细胞间的关系，因此UMAP更好的保留并反映了细胞群潜在的拓扑结构，对于细胞轨迹推断（trajectory inference）分析来说更实用 [4]。
在下图中，我们可以看到UMAP的可视化更趋向于一种紧凑的视觉效果，群簇之间的空间更大，也保留了更多的global structure，因此在选择可视化图中大家可以根据具体的需要来选择。根据小编的实战经验，如果对轨迹推断没有要求，只是看细胞簇群，UMAP的效果会更干净，但是在细胞数目非常大情况下，由于UMAP最大限度的保留了全局结构，这也使得每个簇群的分辨率降低，可能会使得簇群重叠，从而遮盖一些小的簇群，而t-SNE通常能将所有簇群尽可能的“铺开”，所以这种情况建议大家两种都画，然后比较一下。从速度方面来说，同一个数据UMAP的速度要比t-SNE快，这也是UMAP变得更受欢迎的重要原因。
9、降维之后的可视化DimPlot、featureplot等等。
处理完上述数据之后，可视化就很简单了。
具体的参数自行查阅。
以上。
欢迎查漏补缺。
网友评论