9 Seurat
Seurat最初是作为scRNA-seq数据的聚类工具，但是在过去的几年中，该软件包包含越来越多的功能。目前Seurat是一个流行的R软件包，可以对scRNA-seq数进行质量控制，分析和探索，涵盖本课程的许多任务。
注意 超过\(5000\)个细胞的数据集推荐使用Seurat。对于较小的数据集，可以使用SC3。
We recommend using Seurat for datasets with more than \(5000\) cells. For smaller dataset a good alternative will be SC3.
注意 本章使用是 该教程.
9.1 构建Seurat对象
我们将分析10X Genomics免费提供的Peripheral Blood Mononuclear Cells (PBMC)数据集。用Illumina NextSeq 500对2,700个细胞进行了测序。原始数据在 这里。
从读入数据开始，Seurat中的所有功能都已配置为与稀疏矩阵配合使用，从而为Drop-seq/inDrop/10x数据节省大量内存和速度。
library(Seurat)
library(dplyr)
library(cowplot)
# Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "data/pbmc3k_filtered_gene_bc_matrices/hg19/")
# Examine the memory savings between regular and sparse matrices
dense.size <- object.size(x = as.matrix(x = pbmc.data))
dense.size
sparse.size <- object.size(x = pbmc.data)
sparse.size
dense.size/sparse.size
# Initialize the Seurat object with the raw (non-normalized data). Keep all
# genes expressed in >= 3 cells (~0.1% of the data). Keep all cells with at
# least 200 detected genes
pbmc <- CreateSeuratObject(counts = pbmc.data, min.cells = 3, min.features = 200, project = "10X_PBMC", assay = "RNA")9.2 标准预处理流程
以下步骤包括Seurat中scRNA-seq数据的标准预处理流程。包括创建Seurat对象，基于QC指标的细胞选择和过滤，数据标准化和缩放，以及高可变基因的检测。
9.2.1 质控和选择细胞
虽然CreateSeuratObject使用基本的minimum gene-cutoff，但可以根据技术或生物参数在此阶段过滤掉细胞。 Seurat允许用户定义的标准探索QC指标和过滤细胞。在下面的例子中，对基因和分子counts进行可视化，绘制它们之间的关系，并排除multiplets具有异常数量的基因的细胞。当然，这不是排除细胞doublets的保证方法，但将其作为用户定义过滤异常细胞的示例。 我们还根据存在的线粒体基因的百分比过滤细胞。
# The number of genes and UMIs (nGene and nUMI) are automatically calculated
# for every object by Seurat. For non-UMI data, nUMI represents the sum of
# the non-normalized values within a cell We calculate the percentage of
# mitochondrial genes here and store it in percent.mito using AddMetaData.
# We use object@raw.data since this represents non-transformed and
# non-log-normalized counts The % of UMI mapping to MT-genes is a common
# scRNA-seq QC metric.
mito.genes <- grep(pattern = "^MT-", x = rownames(pbmc@assays[["RNA"]]), value = TRUE)
percent.mito <- Matrix::colSums(pbmc@assays[["RNA"]][mito.genes, ])/Matrix::colSums(pbmc@assays[["RNA"]])
# AddMetaData adds columns to object@meta.data, and is a great place to
# stash QC stats
#Seurat v2 function, but shows compatibility in Seurat v3
pbmc <- AddMetaData(object = pbmc, metadata = percent.mito, col.name = "percent.mito")
#in case the above function does not work simply do:
pbmc$percent.mito <- percent.mito
VlnPlot(object = pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"), ncol = 3)
# GenePlot is typically used to visualize gene-gene relationships, but can
# be used for anything calculated by the object, i.e. columns in
# object@meta.data, PC scores etc. Since there is a rare subset of cells
# with an outlier level of high mitochondrial percentage and also low UMI
# content, we filter these as well
par(mfrow = c(1, 2))
FeatureScatter(object = pbmc, feature1 = "nCount_RNA", feature2 = "percent.mito")
# We filter out cells that have unique gene counts (nFeature_RNA) over 2,500 or less than
# 200 Note that > and < are used to define a'gate'.
#-Inf and Inf should be used if you don't want a lower or upper threshold.
pbmc <- subset(x = pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mito > -Inf & percent.mito < 0.05 )9.3 数据标准化
从数据集中删除不需要的细胞后，下一步是数据标准化。 默认情况下，采用全局缩放归一化方法“LogNormalize”，通过总表达对每个细胞的基因表达值进行标准化，乘以比例因子（默认为10,000），并对结果进行对数转换。
9.4 检测高可变基因
Seurat计算高度可变的基因，并关注下游分析。FindVariableGenes计算每个基因的表达均值和离散程度，将这些基因放入bin中，然后计算每个bin离散程度的z-score。这有助于控制方差和表达均值之间的关系。(Macosko et al.)没有更改该函数，但很快就会出现用于HVG鉴定的新方法。 建议用户设置这些参数以在dispersion plot标记异常值，但具体的参数设置可能会根据数据类型，样本中的异质性和标准化策略而有所不同。 这里的参数识别~2000个可变基因，并且代表UMI数据的典型参数设置，其归一化为总共1e4个分子。
查看FindVariableFeatures输出的输出。基因不存储在对象中，但可以通过这种方式访问。
head(x = HVFInfo(object = pbmc))
## mean variance variance.standardized
## AL627309.1 0.003411676 0.003401325 0.9330441
## AP006222.2 0.001137225 0.001136363 0.9924937
## RP11-206L10.2 0.001895375 0.001892500 0.9627290
## RP11-206L10.9 0.001137225 0.001136363 0.9924937
## LINC00115 0.006823351 0.006779363 0.9062135
## NOC2L 0.107278241 0.159514698 0.78493099.5 缩放数据并删除不需要的变异来源
单细胞数据集可能包含“不感兴趣”的变异来源。这不仅包括技术噪音，还包括批次效应，甚至包括生物变异来源（细胞周期）。 正如 Buettner et al, NBT, 2015所建议的那样，从分析中回归这些信号可以改善下游降维和聚类结果。 为了减轻这些信号的影响，Seurat基于用户定义的变量构建线性模型预测基因表达。这些模型的缩放z-score残差存储在scale.data中，用于降维和聚类。
我们可以建立由批次效应驱动的基因表达中的细胞-细胞变异，细胞比对率（由Drop-seq数据的Drop-seq工具提供），检测到的分子数量和线粒体基因表达回归模型。对于循环细胞，我们还可以学习‘cell-cycle’ score(示例) 并对其进行回归。 在这个有丝分裂后血细胞的简单例子中，我们对每个细胞检测到的分子数量以及线粒体基因含量百分比建立回归模型。
Seurat v2.0将此回归实现为数据扩展过程的一部分。RegressOut函数已经被弃用，并将其替换为ScaleData中的vars.to.regress参数。
9.6 线性降维
–> Seurat v2: 接下来对缩放后的数据执行PCA。 默认情况下，object@var.genes中的基因用作输入，但可以使用pc.genes定义。 通常对高度可变基因运行降维可以提高性能。对于UMI数据，特别是在消除技术因素之后，PCA在更大的基因集（包括整个转录组）上运行时返回相似（尽管更慢）的结果。
–> refered to Seurat v3 (latest): 高可变基因可通过函数HVFInfo(object)访问。尽管RunPCA有一个features参数，用于指定计算PCA的features，但测试修改该参数值，PCA结果相同，表明features参数中提供的基因并不完全是用于计算PCA。不知道该函数是直接获取HVG还是不考虑它们。
–> refered to Seurat v2: Seurat提供了几种可视化PCA中细胞和基因的方法，包括PrintPCA, VizPCA, PCAPlot, 和 PCHeatmap
–> Seurat v3 (latest): Seurat v3 提供以下可视化的方法: - PCA - 定量特征着色的PCA图 - 单个细胞散点图 - 各基因的散点图 - 小提琴和ridge图 - Violin and Ridge plots - 热图
# Dimensional reduction plot, with cells colored by a quantitative feature
FeaturePlot(object = pbmc, features = "MS4A1")
# Scatter plot across single cells, replaces GenePlot
FeatureScatter(object = pbmc, feature1 = "MS4A1", feature2 = "PC_1")
# Scatter plot across individual features, repleaces CellPlot
CellScatter(object = pbmc, cell1 = "AGTCTACTAGGGTG", cell2 = "CACAGATGGTTTCT")
DimHeatmap允许探索数据集中异质性的主要来源，并且确定下游分析包含的PCs。细胞和基因根据其PCA分数排序。将cells.use设置为数字可以绘制两端的“extreme”细胞，从而大大加快了大型数据集的绘图速度。虽然一种监督分析，但是为探索相关基因集的有用工具。
Seurat 3.0中不再支持ProjectPCA函数。
9.7 统计显著主成分
为克服scRNA-seq数据广泛技术噪音，Seurat根据其PCA评分对细胞进行聚类，每个PC基本上代表相关基因集中信息组合的“metagene”。 因此，确定下游分析使用多少PC是重要的一步。
Macosko et al受jackStraw程序启发实现重采样检验。 随机扰动数据的一个子集（默认为1％）并重新运行PCA，构建基因得分的“零分布”，并重复此过程。富集了低P值基因集的PC为显著PC。
# NOTE: This process can take a long time for big datasets, comment out for
# expediency. More approximate techniques such as those implemented in
# PCElbowPlot() can be used to reduce computation time
pbmc <- JackStraw(object = pbmc, reduction = "pca", dims = 20, num.replicate = 100, prop.freq = 0.1, verbose = FALSE)JackStrawPlot函数提供了一个可视化工具，用于比较PC的p值分布和均匀分布（虚线）。显著PC将显示强烈富集低p值的基因（虚线上方的实线）。 在这种情况下，PC 1-10都很重要。
确定使用哪些PC的更特别的方法是查看主成分的标准偏差图，并在图中拐点地方绘制cutoff。这可以通过ElbowPlot来完成。 在这个例子中，看起来拐点落在PC5周围。
选择PC,识别数据集的真实维度，是Seurat的重要一步，但对用户来说可能具有挑战性/不确定性。因此，我们建议考虑这三种方法。第一个是更多的监督分析，探索PC以确定异质性的相关来源，并且可以与GSEA一起使用。第二个基于随机空模型的统计检验，但对于大型数据集来说耗时的，并且可能不会返回明确的PC cutoff。第三种是常用的启发式算法，可以立即计算。在这个例子中，所有三种方法产生了类似的结果，可以在PC 7-10之间选择cutoff。这里采用jackStraw结果，看到PCHeatmap在这些PC中返回可解释的信号（包括canonical dendritic cell markers）。虽然结果只是受到cutoff的微小变化的影响（可以在下面进行测试），但强烈建议您始终探索选择进行下游分析的PCs。
9.8 聚类
Seurat现在包括基于图形的聚类方法，虽然聚类分析的距离度量（基于先前识别的PC）保持不变，但将细胞距离矩阵分成cluster的方法已经大大改进。我们的方法受到最近scRNA-seq数据基于图聚类方法 SNN-Cliq, Xu and Su, Bioinformatics, 2015和CyTOF data PhenoGraph, Levine et al., Cell, 2015的启发。简而言之，这些方法将细胞嵌入到图结构中 - 例如K-最近邻（KNN），相邻的细胞具有相似的基因表达模式，然后将该图划分为高度互连的“quasi-cliques”或“quasi-cliques”。与PhenoGraph一样，我们首先根据PCA空间中的欧氏距离构建KNN图，并根据其局部邻域中的重叠部分（Jaccard相似性）优化任意两个细胞之间的边权。为了聚类细胞，应用模块化优化技术，例如Louvain算法（默认）或SLM SLM, Blondel et al., Journal of Statistical Mechanics以迭代方式将细胞分组在一起来优化标准模块化函数。
FindClusters函数实现上述过程，并包含一个分辨率参数，用于设置下游聚类的“granularity”，值越大cluster越多。 将此参数设置在0.6-1.2之间对包含3k细胞的单细胞数据集返回不错的记过。较大的数据集通常会增加最佳分辨率。 最新的聚类结果将存储在seurat_clusters的元数据中。
首先计算k-最近邻并构造SNN图(FindNeighbors)，然后运行FindClusters。
9.9 非线性降维
Seurat继续使用tSNE作为可视化和探索这些数据集的强大工具。虽然不再建议直接在tSNE上进行聚类，上述基于图聚类cluster中的细胞在tSNE图中应该定位在一起。这是因为tSNE旨在将高维空间中具有相似局部邻域的细胞在低维空间放在一起。tSNE的输入，建议使用与聚类分析相同的PC，尽管genes.use参数也支持基于缩放基因表达计算tSNE。
pbmc <- RunTSNE(object = pbmc, dims.use = 1:10, do.fast = TRUE)
# note that you can set do.label=T to help label individual clusters
DimPlot(object = pbmc, reduction = "tsne")
9.10 运行UMAP
为了并排显示这两个条件，可以使用split.by参数来显示由cluster着色的每个条件。 To visualize the two conditions side-by-side, we can use the split.by argument to show each condition colored by cluster.
pbmc <- RunUMAP(pbmc, reduction = "pca", dims = 1:20)
DimPlot(pbmc, reduction = "umap", split.by = "seurat_clusters")
您可以在此时保存对象，以便可以轻松地将其加载或者与协作者共享，而无需重新运行上面执行的计算密集型步骤。
9.11 差异表达分析
Seurat可以通过差异表达找到cluster的makers。默认情况下，通过与所有其他细胞进行比较，识别单个cluster的（在ident.1中指定）的正负marker。 FindAllMarkers为所有cluster自动执行此过程，但也可以检验cluster之间或所有细胞间。
min.pct参数要求在两组细胞中的任一组中最小检测到基因的百分比，并且thresh.test参数要求基因在两组之间以一定量差异表达（平均）。 可以将这两者设置为0，但随着时间的显著增加 - 因为这将检验大量不太可能具有高度差异基因。 加速这些计算的另一种选择，可以设置max.cells.per.ident。这将对每个类进行下采样，使其不多于设置的细胞。 虽然这会损失效力，但显著增加速度，并且最高度差异表达的基因可能仍会上升到顶部。
# find all markers of cluster 1
cluster1.markers <- FindMarkers(object = pbmc, ident.1 = 1, min.pct = 0.25)
print(x = head(x = cluster1.markers, n = 5))# find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(object = pbmc, ident.1 = 2, ident.2 = c(0, 3), min.pct = 0.25)
print(x = head(x = cluster5.markers, n = 5))# find markers for every cluster compared to all remaining cells, report
# only the positive ones
pbmc.markers <- FindAllMarkers(object = pbmc, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25)
pbmc.markers %>% group_by(cluster) %>% top_n(2, avg_logFC)Seurat有多种差异表达检验，可以使用test.use参数设置（详细信息见 DE vignette）。 例如，ROC检验返回marker的“classification power”（范围从0 - random, to 1 - perfect）。
Seurat提供了几种用于可视化marker表达的工具。
* VlnPlot (显示cluster见的表达概率分布),
* FeaturePlot (tSNE或PCA可视化基因表达，也是最常用的),
建议尝试以下可视化方法：
* RidgePlot,
* CellPlot
* DotPlot
DoHeatmap为给定的细胞和基因生成表达热图。在这种情况下，为每个cluster绘制前20个markers（如果小于20，绘制所有marker）。
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(10, avg_logFC)
# setting slim.col.label to TRUE will print just the cluster IDS instead of
# every cell name
DoHeatmap(object = pbmc, features = top10$gene, label = TRUE)
9.12 确定cluster细胞类型
对于该数据集，使用经典marker确定cluster细胞类型。 Fortunately in the case of this dataset, we can use canonical markers to easily match the unbiased clustering to known cell types.
current.cluster.ids <- c(0, 1, 2, 3, 4, 5, 6, 7)
new.cluster.ids <- c("CD4 T cells", "CD14+ Monocytes", "B cells", "CD8 T cells", "FCGR3A+ Monocytes", "NK cells", "Dendritic cells", "Megakaryocytes")
pbmc@active.ident <- plyr::mapvalues(x = pbmc@active.ident, from = current.cluster.ids, to = new.cluster.ids)
DimPlot(object = pbmc, reduction = "tsne", do.label = TRUE, pt.size = 0.5)
9.13 细分细胞类型
如果扰乱我们上面的一些参数选择（例如，设置resolution=0.8或更改PC的数量），可能会看到CD4 T细胞细分为两个群体。 您可以探索此细分以查找区分两个T细胞marker基因。但是，在重新聚类(覆盖object@ident)之前，可以存储我们重命名的identities，以便以后恢复。
# First lets stash our identities for later
pbmc <- StashIdent(object = pbmc, save.name = "ClusterNames_0.6")
# Note that if you set save.snn=T above, you don't need to recalculate the
# SNN, and can simply put: pbmc <- FindClusters(pbmc,resolution = 0.8)
pbmc <- FindClusters(object = pbmc, reduction.type = "pca", dims.use = 1:10, resolution = 0.8, print.output = FALSE)# Demonstration of how to plot two tSNE plots side by side, and how to color
# points based on different criteria
plot1 <- DimPlot(object = pbmc, reduction = "tsne", do.return = TRUE, no.legend = TRUE, do.label = TRUE)
plot2 <- DimPlot(object = pbmc, reduction = "tsne", do.return = TRUE, group.by = "ClusterNames_0.6", no.legend = TRUE, do.label = TRUE)
plot_grid(plot1, plot2)
# Find discriminating markers
tcell.markers <- FindMarkers(object = pbmc, ident.1 = 0, ident.2 = 1)
# Most of the markers tend to be expressed in C1 (i.e. S100A4). However, we
# can see that CCR7 is upregulated in C0, strongly indicating that we can
# differentiate memory from naive CD4 cells. cols.use demarcates the color
# palette from low to high expression
FeaturePlot(object = pbmc, features = c("S100A4", "CCR7"), cols = c("green", "blue"))
memory/naive划分有点弱，可以查看更多细胞看看这是否更有说服力。与此同时，我们可以恢复原来cluster identities以进行下游处理。
9.14 sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 7 x64 (build 7601) Service Pack 1
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Chinese (Simplified)_People's Republic of China.936
## [2] LC_CTYPE=Chinese (Simplified)_People's Republic of China.936
## [3] LC_MONETARY=Chinese (Simplified)_People's Republic of China.936
## [4] LC_NUMERIC=C
## [5] LC_TIME=Chinese (Simplified)_People's Republic of China.936
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] cowplot_1.0.0 dplyr_0.8.3 Seurat_3.0.2 knitr_1.23
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.0 tidyr_0.8.3 viridisLite_0.3.0
## [4] jsonlite_1.6 splines_3.6.0 lsei_1.2-0
## [7] R.utils_2.9.0 gtools_3.8.1 Rdpack_0.11-0
## [10] assertthat_0.2.1 yaml_2.2.0 ggrepel_0.8.1
## [13] globals_0.12.4 pillar_1.4.2 lattice_0.20-38
## [16] reticulate_1.12 glue_1.3.1 digest_0.6.20
## [19] RColorBrewer_1.1-2 SDMTools_1.1-221.1 colorspace_1.4-1
## [22] htmltools_0.3.6 Matrix_1.2-17 R.oo_1.22.0
## [25] plyr_1.8.4 pkgconfig_2.0.2 bibtex_0.4.2
## [28] tsne_0.1-3 listenv_0.7.0 bookdown_0.12
## [31] purrr_0.3.2 scales_1.0.0 RANN_2.6.1
## [34] gdata_2.18.0 Rtsne_0.15 tibble_2.1.3
## [37] ggplot2_3.2.0 ROCR_1.0-7 pbapply_1.4-0
## [40] lazyeval_0.2.2 survival_2.44-1.1 magrittr_1.5
## [43] crayon_1.3.4 evaluate_0.14 R.methodsS3_1.7.1
## [46] future_1.14.0 nlme_3.1-140 MASS_7.3-51.4
## [49] gplots_3.0.1.1 ica_1.0-2 tools_3.6.0
## [52] fitdistrplus_1.0-14 data.table_1.12.2 gbRd_0.4-11
## [55] stringr_1.4.0 plotly_4.9.0 munsell_0.5.0
## [58] cluster_2.1.0 irlba_2.3.3 compiler_3.6.0
## [61] rsvd_1.0.1 caTools_1.17.1.2 rlang_0.4.0
## [64] grid_3.6.0 ggridges_0.5.1 htmlwidgets_1.3
## [67] igraph_1.2.4.1 bitops_1.0-6 rmarkdown_1.14
## [70] npsurv_0.4-0 gtable_0.3.0 codetools_0.2-16
## [73] reshape2_1.4.3 R6_2.4.0 gridExtra_2.3
## [76] zoo_1.8-6 future.apply_1.3.0 KernSmooth_2.23-15
## [79] metap_1.1 ape_5.3 stringi_1.4.3
## [82] parallel_3.6.0 Rcpp_1.0.1 sctransform_0.2.0
## [85] png_0.1-7 tidyselect_0.2.5 xfun_0.8
## [88] lmtest_0.9-37