Seurat简介
Seurat---几乎是当前单细胞RNA-seq分析领域的不可或缺的工具，特别是基于10X公司的cellrange流程得出的结果，可以方便的对接到Seurat工具中进行后续处理，简直是带给迷茫在单细胞数据荒漠中小白的一眼清泉，相对全面的功能，简洁的操作命令，如丝般顺滑。
更新记录
Seurat目前最新版为V3，第一版是为处理空间转录组而设计，第二版针对基因droplet的单细胞技术而开发的用于单细胞质控，降维，聚类，鉴定细胞类型marker等多重功能，目前更新的第三版可以整合多组学，多批次，多实验方法的数据，把握了后续数据分析需要多组学整合的热点策略，符合下一步分析思路更全面，生物学解释更透彻的科研需求。
以下为Seurat更新的记录：
2019年4月16日发布3.0版
2018年11月2日版本3.0 alpha发布
- 发布的预印本描述了在单细胞数据集中识别“锚点”的新方法
2018年3月23日 2.3版本发布 - 提高速度和内存效率
- 用于分析来自Microwell-seq Mouse Cell Atlas数据集的~250,000个细胞的新视图
2018年1月10日： 2.2版本发布 - 支持多数据集对齐
- 评估对齐性能的新方法
2017年10月16日： 2.1版发布 - 支持多模式单细胞数据
- 支持用于差异表达式测试的MAST和DESeq2包
2017年7月26日： 2.0版本发布 - 发布预印本用于scRNA-seq数据集的综合分析
- 数据集集成，可视化和探索的新方法
- 对代码库进行重大重组，以强调清晰度和清晰的文档
2016年10月4日： 1.4版本发布 - 为UMI计数数据添加了负二项回归和差分表达式测试的方法
- 合并和下采样Seurat对象的新方法
2016年8月22日： 1.3版本发布 - 改进的聚类方法 - 有关详细信息，请参阅FAQ
- 所有函数都支持稀疏矩阵
- 消除不需要的变异来源的方法
- 一致的函数名称
- 更新的可视化
2015年5月21日： Drop-Seq手稿发布。版本1.2发布 - 增加了对光谱t-SNE（非线性降维）和密度聚类的支持
- 新的可视化 - 包括pcHeatmap，dot.plot和feature.plot
- 扩展包文档，减少导入包负担
- Seurat代码现在托管在GitHub上，可以通过devtools包轻松安装
- 小错误修复
2015年4月13日： 空间映射手稿发布。版本1.1发布
Seurat安装
直接在CRAN加载安装
# Enter commands in R (or R studio, if installed)
install.packages('Seurat')
library(Seurat)
同时SeuratV3支持使用UMAP算法进行降维聚类可视化，需要配置对应的pythonUMAP环境，方可执行UMAP对应功能。
Seurat分析单细胞RNA-seq数据
官网提供了基于10X公司的PBMC样本的数据作为示范；
首先读入实例数据，利用Read10X函数读取cellranger流程的输出文件，返回唯一的分子识别（UMI）计数矩阵。该矩阵中的值表示在每个细胞（列）中检测到的每个特征（即基因;行）的分子数。然后使用计数矩阵来创建一个Seurat对象。该对象充当容器，其包含单细胞数据集的数据（如计数矩阵）和分析（如PCA或聚类结果）
library(dplyr)
library(Seurat)
setwd("/Users/liuyang/Downloads/03project/02scRNA/10X_Seurat")
rm(list = ls())
# Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "/Users/liuyang/Downloads/03project/02scRNA/10X_Seurat/filtered_gene_bc_matrices/hg19/")
# Initialize the Seurat object with the raw (non-normalized data).
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc
标准的预处理工作流程基于QC指标，数据标准化和缩放以及高度可变特征的检测来选择和过滤细胞。
QC并选择细胞进行进一步分析
Seurat允许用户根据自己定义的标准轻松探索QC指标和过滤单元格。常用的一些QC指标包括
- 在每个细胞中检测到的独特基因的数量。
- 低质量细胞或空液滴通常只有很少的基因
- 细胞双峰或多重峰可能表现出异常高的基因计数
- 同样，细胞内检测到的分子总数（与独特基因强烈相关）
- 映射到线粒体基因组的读数百分比（对于已经组装到染色体水平或有线粒体基因组的样本）
- 低质量/垂死细胞通常表现出广泛的线粒体污染
- 我们使用
PercentageFeatureSet函数计算线粒体QC指标，该函数计算源自一组特征的计数百分比 - 我们使用所有基因
MT-集作为一组线粒体基因
pbmc[["percent.mt"]] <- PercentageFeatureSet(object = pbmc, pattern = "^MT-")
VlnPlot(object = pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
# FeatureScatter is typically used to visualize feature-feature relationships, but can be used for anything calculated by
# the object, i.e. columns in object metadata, PC scores etc.
plot1 <- FeatureScatter(object = pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(object = pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))
标准化数据
从数据集中删除不需要的单元格后，下一步是标准化数据。默认情况下，采用全局缩放规范化方法LogNormalize，通过总表达式对每个单元格的要素表达式度量进行标准化，将其乘以比例因子（默认为10,000），并对结果进行对数转换。归一化值存储在pbmc[["RNA"]]@data。
pbmc <- NormalizeData(object = pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
识别高度可变的特征（特征选择）
然后计算在数据集中表现出高细胞间差异的特征子集（即，它们在一些细胞中高度表达，而在其他细胞中低表达）。作者和其他人研究者发现，在下游分析中关注这些基因有助于突出单细胞数据集中的生物信号。
这里详细描述了SeuratV3中的程序，并通过直接建模单细胞数据中固有的均值 - 方差关系对先前版本进行了改进，并在FindVariableFeatures函数中实现。默认情况下，为每个数据集返回2,000个要素。这些将用于下游分析，如PCA。
pbmc <- FindVariableFeatures(object = pbmc, selection.method = "vst", nfeatures = 2000)
# Identify the 10 most highly variable genes
top10 <- head(x = VariableFeatures(object = pbmc), 10)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(object = pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
CombinePlots(plots = list(plot1, plot2))
缩放数据
接下来，应用线性变换（'缩放'），这是在PCA等降维技术之前的标准预处理步骤。该ScaleData函数：
改变每个基因的表达，使得跨细胞的平均表达为0
缩放每个基因的表达，以便跨细胞的方差为1
该步骤在下游分析中给予相同的权重，因此高表达的基因不占优势结果存储在 pbmc[["RNA"]]@scale.data
执行线性降维聚类
对缩放数据执行PCA。默认情况下，只有先前确定的变量要素用作输入，但features如果要选择其他子集，则可以使用参数定义。
瑟拉提供定义可视化细胞和功能的几种有用的方式PCA，包括VizDimReduction，DimPlot，和DimHeatmap
VizDimLoadings(object = pbmc, dims = 1:2, reduction = "pca")
DimPlot(object = pbmc, reduction = "pca")
特别是DimHeatmap允许容易地探索数据集中的异质性的主要来源，并且在试图决定包括哪些PC以用于进一步的下游分析时可以是有用的。细胞和特征都根据其PCA分数排序。设置cells为数字可以绘制光谱两端的“极端”单元格，从而大大加快了大型数据集的绘图速度。虽然显然是一种监督分析，但我们发现这是探索相关特征集的有用工具。
DimHeatmap(object = pbmc, dims = 1, cells = 500, balanced = TRUE)
确定数据集的“维度”
为了克服scRNA-seq数据的任何单一特征中的广泛技术噪声，Seurat基于其PCA分数对细胞进行聚类，每个PC基本上代表在相关特征集中组合信息的“元特征”。因此，顶部主要组件表示数据集的强大压缩。但是，我们应该选择包含多少个组件？10？20？100？在Macosko 等人中，实施了一项受JackStraw程序启发的重采样测试。我们随机置换数据的子集（默认为1％）并重新运行PCA，构建特征分数的“空分布”，并重复此过程。我们将“重要”PC识别为具有低p值特征的强大丰富的PC。
# NOTE: This process can take a long time for big datasets, comment out for expediency. More approximate techniques such
# as those implemented in ElbowPlot() can be used to reduce computation time
pbmc <- JackStraw(object = pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(object = pbmc, dims = 1:20)
该JackStrawPlot功能提供了一种可视化工具，用于比较每台PC的p值分布和均匀分布（虚线）。“重要的”PC将显示具有低p值的特征的强烈丰富（虚线上方的实线）。在这种情况下，看起来在前10-12PC之后，重要性急剧下降。
JackStrawPlot(object = pbmc, dims = 1:15)
另一种启发式方法生成“肘图”：基于每个（ElbowPlot函数）解释的方差百分比对主成分进行排序。在这个例子中，我们可以观察PC9-10周围的“肘部”，这表明大多数真实信号都是在前10台PC中捕获的。
ElbowPlot(object = pbmc)
识别数据集的真实维度 - 对用户来说可能具有挑战性/不确定性。因此，我们建议考虑这三种方法。第一个是更多的监督，探索PC以确定异质性的相关来源，并且可以与GSEA一起使用。第二个实现基于随机空模型的统计测试，但对于大型数据集来说是耗时的，并且可能不会返回明确的PC截止。第三种是常用的启发式算法，可以立即计算。在这个例子中，所有三种方法产生了类似的结果，但我们可能有理由在PC 7-12之间选择任何作为截止值的东西。
我们在这里选择了10，但鼓励用户考虑以下内容：
树突细胞和NK爱好者可以认识到与PC12和13强烈相关的基因定义了罕见的免疫亚群（即MZB1是浆细胞样DC的标记）。然而，这些组是如此罕见，如果没有先验知识，很难将这种大小的数据集与背景噪声区分开来。
我们鼓励用户使用不同数量的PC（10,15甚至50！）重复下游分析。正如您将观察到的，结果通常没有显着差异。
在选择此参数时，我们建议用户在较高的一侧犯错。例如，仅使用5台PC进行下游分析确实会对结果造成严重影响。
聚类细胞
Seurat v3采用基于图形的聚类方法，建立在（Macosko 等人）的初始策略之上。重要的是，驱动聚类分析的距离度量（基于先前识别的PC）保持不变。然而，我们将细胞距离矩阵分成簇的方法已经大大改进。我们的方法受到最近手稿的启发，这些手稿将基于图形的聚类方法应用于scRNA-seq数据[SNN-Cliq，Xu和Su，Bioinformatics，2015]和CyTOF数据[PhenoGraph，Levine 等，Cell，2015]。简而言之，这些方法将单元格嵌入图形结构中 - 例如K-最近邻（KNN）图形，在具有相似特征表达模式的单元格之间绘制边缘，然后尝试将此图形划分为高度互连的“准集团”或“社区”。
与PhenoGraph一样，我们首先根据PCA空间中的欧氏距离构建KNN图，并根据其局部邻域中的共享重叠（Jaccard相似性）细化任意两个单元之间的边权重。使用该FindNeighbors函数执行该步骤，并将先前定义的数据集维度（前10个PC）作为输入。
为了聚类细胞，我们接下来应用模块化优化技术，例如Louvain算法（默认）或SLM [SLM，Blondel 等，Journal of Statistical Mechanics]，将细胞迭代分组在一起，目的是优化标准模块化功能。该FindClusters函数实现此过程，并包含一个分辨率参数，用于设置下游群集的“粒度”，增加的值会导致更多的群集。我们发现将此参数设置在0.4-1.2之间通常会返回大约3K单元格的单细胞数据集的良好结果。较大的数据集通常会增加最佳分辨率。可以使用该Idents功能找到群集。
pbmc <- FindNeighbors(object = pbmc, dims = 1:10)
pbmc <- FindClusters(object = pbmc, resolution = 0.5)
运行非线性降维（UMAP / tSNE）
Seurat提供了几种非线性降维技术，如tSNE和UMAP，可视化和探索这些数据集。这些算法的目标是学习数据的基础流形，以便在低维空间中将相似的单元放在一起。上面确定的基于图的聚类内的单元应该共同定位这些降维图。作为UMAP和tSNE的输入，我们建议使用相同的PC作为聚类分析的输入。
# If you haven't installed UMAP, you can do so via reticulate::py_install(packages = 'umap-learn')
pbmc <- RunUMAP(object = pbmc, dims = 1:10)
# note that you can set `label = TRUE` or use the LabelClusters function to help label individual clusters
DimPlot(object = pbmc, reduction = "umap")
寻找差异表达的特征（聚类生物标志物）
Seurat可以帮助您找到通过差异表达定义聚类的标记。默认情况下，它标识单个群集（指定ident.1）的正负标记，与所有其他单元格进行比较。FindAllMarkers为所有群集自动执行此过程，但您也可以测试群集彼此之间或所有单元格。
该min.pct论证要求在两组单元中的任一组中以最小百分比检测特征，并且thresh.test参数要求在两组之间以一定量差异地表示（平均）特征。您可以将这两个设置为0，但时间会急剧增加 - 因为这将测试大量不太可能具有高度歧视性的功能。作为加速这些计算的另一种选择，max.cells.per.ident可以设置。这将对每个标识类进行下采样，使其不再具有比设置的更多的单元格。虽然通常会出现功率损失，但速度增加可能是显着的，并且最高度差异表达的功能可能仍然会升至顶部。
cluster1.markers <- FindMarkers(object = pbmc, ident.1 = 1, min.pct = 0.25)
head(x = cluster1.markers, n = 5)
Seurat有几个差分表达式测试，可以使用test.use参数设置
cluster1.markers <- FindMarkers(object = pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
我们提供了几种用于可视化标记表达的工具。VlnPlot（显示跨群集的表达概率分布），以及FeaturePlot（在tSNE或PCA图上可视化特征表达）是我们最常用的可视化。我们还建议您探索RidgePlot，CellScatter以及DotPlot查看数据集的其他方法。
VlnPlot(object = pbmc, features = c("MS4A1", "CD79A"))
FeaturePlot(object = pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
DoHeatmap为给定的单元格和要素生成表达式热图。在这种情况下，我们正在为每个群集绘制前20个标记（或所有标记，如果小于20）。
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(object = pbmc, features = top10$gene) + NoLegend()
将单元类型标识分配给集群
根据每个类群中高表达的特异基因，我们可以利用已知细胞类型特定的标记基因，对每个类群确定细胞类型。
new.cluster.ids <- c("Memory CD4 T", "CD14+ Mono", "Naive CD4 T", "B", "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Mk")
names(x = new.cluster.ids) <- levels(x = pbmc)
pbmc <- RenameIdents(object = pbmc, new.cluster.ids)
DimPlot(object = pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
