简介

  Seurat是一个r包，被设计用于单细胞rna-seq数据的细胞质控和分析。Seurat旨在使用户能够识别和解释单细胞转录组数据中的异质性来源，同时提供整合不同类型的单细胞数据的函数。目前Seurat软件版本已更新到V3。
  上篇介绍了单细胞rna-seq分析的第一个步骤——数据质控和定量（10X单细胞测序之cellranger介绍）。接下来将介绍下第二个步骤——细胞质控、聚类及差异分析。这些分析由Seurat软件全部可以完成，同时还提供了可视化结果，功能非常强大。

---

分析

  Seurat可以分析多个方向内容，如下图所示，像单个样品的标准分析、多个样品合并分析、空间转录组分析、scRNA-seq 和scATAC-seq的联合分析，同时还提供了被用户经常请求的可视化图，说实话，这软件真的很亲民，太为用户着想了。

image.png
  由于这段时间使用Seurat只进行了单个样品和多个样品的分析，因此本次只介绍这块的分析内容，像其他空间转录组、ATAC联合分析的内容待后续学习后更新。

单样品分析

创建Seurat对象

  在我们使用cellranger分析完成后，得到了每个样品中的细胞表达矩阵，如下所示：

├── filtered_feature_bc_matrix #过滤后矩阵信息
│   ├── barcodes.tsv.gz  #过滤（细胞中表达基因的数目在阈值内）后的细胞总数文件
│   ├── features.tsv.gz #所有细胞表达基因的并集
│   └── matrix.mtx.gz #坐标文件，第一列是表达基因数的坐标，第二轮是count值，第三列是细胞

  我们接下来直接可以使用这些数据创建Seurat对象，创建方法如下：

library(dplyr)
library(Seurat)

# 使用Read10X函数读取矩阵数据，得到一个稀疏矩阵
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_feature_bc_matrix/")
# 使用CreateSeuratObject函数创建Seurat对象，这里要求细胞中基因数目不能小于200，且基因至少在三个细胞中有表达，否则过滤掉
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc
## An object of class Seurat 
## 13714 features across 2700 samples within 1 assay 
## Active assay: RNA (13714 features)

  Seurat对象就是一个S4类，里面装着单细胞数据集，如count矩阵、细胞矩阵、聚类信息都存储在这个容器中。我们可以通过查看这个对象来得到我们想要的信息。如我们想看一下前5个细胞的前5个基因的表达矩阵，可以这样使用：

pbmc@assays$RNA@counts[1:5,1:5]

细胞选择

  由于低质量细胞或空液滴通常只有很少的基因、双细胞可能表现出异常高的基因数目、低质量/濒死细胞通常表现出广泛的线粒体污染三个原因，我们需要进行细胞数目的选择。使用percentagefeatureset函数计算线粒体qc指标，该函数计算源自一组基因的计数百分比。

#使用从MT-作为线粒体基因集
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
#qc指标可视化
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

image.png

  我们从上图的可视化结果进行细胞选择。我们过滤具有超过2500或少于200个基因计数的细胞，同时过滤掉线粒体比例超过5%的细胞。

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

数据标准化

  从数据集中移除不需要的细胞后，下一步是标准化数据。默认情况下，我们使用“lognormalize”全局缩放的归一化方法，通过总表达值对每个细胞的基因表达值归一化，并将其乘以缩放因子(默认为10,000)，最后对结果进行对数变换，标准化数据存储在pbmc@assays$RNA@data中。

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

高变异基因选择

  计算数据集中表现出高细胞间差异的基因子集(即，它们在一些细胞中高表达，而在另一些细胞中低表达)。在下游分析中关注这些基因有助于突出单细胞数据集中的生物信号。默认情况下，每个数据集选择2000个高变异基因用于下游分析。

pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# 前10的高变异基因
top10 <- head(VariableFeatures(pbmc), 10)
# 将前10的高变异基因在图中标注出来
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
CombinePlots(plots = list(plot1, plot2))

image.png

缩放数据

  使用线性变换(“scaling”)处理数据，过程：改变每个基因的表达，使细胞间的平均表达为0；缩放每个基因的表达，使细胞间的差异为1，这样高表达基因就不会影响后续分析；这步是pca等降维技术之前的标准预处理步骤，缩放的数据储存在pbmc@assays$RNA@scale.data中。

all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

线性降维

  接下来，我们对缩放后的数据执行pca降维，默认降至50维。PCA(Principal Component Analysis)，即主成分分析方法，是一种使用最广泛的数据降维算法。PCA的主要思想是将n维特征映射到k维上，在尽可能保留原始数据信息的同时降低数据维度来加速数据分析。过程就是从原始高维的空间中按顺序地找一组相互正交的坐标轴系统，新的坐标轴的选择与数据本身是密切相关的。其中，第一个新坐标轴选择是原始数据中方差最大的方向，第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的，第三个轴是与第1,2个轴正交的平面中方差最大的。依次类推，可以得到n个这样的坐标轴。大部分方差都包含在前面k个坐标轴中，后面的坐标轴所含的方差几乎为0。于是，我们可以忽略余下的坐标轴，只保留前面k个含有绝大部分方差的坐标轴，实现对数据降维。

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#选择前5个维度进行查看
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)

## PC_ 1 
## Positive:  CST3, TYROBP, LST1, AIF1, FTL 
## Negative:  MALAT1, LTB, IL32, IL7R, CD2 
## PC_ 2 
## Positive:  CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1 
## Negative:  NKG7, PRF1, CST7, GZMB, GZMA 
## PC_ 3 
## Positive:  HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1 
## Negative:  PPBP, PF4, SDPR, SPARC, GNG11 
## PC_ 4 
## Positive:  HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1 
## Negative:  VIM, IL7R, S100A6, IL32, S100A8 
## PC_ 5 
## Positive:  GZMB, NKG7, S100A8, FGFBP2, GNLY 
## Negative:  LTB, IL7R, CKB, VIM, MS4A7

  使用DimPlot可视化函数查看样品中的细胞在pca中的分布情况：

DimPlot(pbmc, reduction = "pca")

image.png

  使用DimHeatmap可视化函数查看样品中前500个细胞在前6个PCA中的热图：

DimHeatmap(pbmc, dims = 1:6, cells = 500, balanced = TRUE)

image.png

确定数据集的维度

  每个pc本质上代表一个“元特征”，它将相关特征集中的信息组合在一起。因此，越在顶部的主成分越可能代表数据集。然而，我们应该选择多少个主成分才认为我们选择的数据包含了绝大部分的原始数据信息呢？我们可以使用JackStraw函数来选择PC数目：

#将数据的1%打乱重新运行pca，构建特征得分的“零分布”，这个过程重复100次
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
plot1<-JackStrawPlot(pbmc, dims = 1:15)
plot2<-ElbowPlot(pbmc)
CombinePlots(plots = list(plot1, plot2))

image.png

  左图是前15个pc的p值分布和均匀分布，越显著的PC，其p值越低，从这图我们可以看出在前10-12的pc之后，分布有一个急剧下降趋势。右图是一个肘型图，我们可以观察到pc9-10周围有一个“肘部”,这表明大部分真实信号是在前10个pc中捕获的，这里选择10个pc进行后续分析。

细胞聚类

  Seurat3使用基于图聚类算法（graph-based clustering approach）进行细胞聚类分析。这部分聚类按我的理解是在pca空间中分布着各个细胞，首先在这个空间中的计算离每个细胞最近的k个细胞的欧氏距离来构造一个个knn图，并基于任意两个细胞在其局部邻域中的共享重叠来细化它们之间的边缘权重，尝试将该图划分为高度互连的社区，最后应用模块化优化技术（默认Louvain 算法）将细胞迭代分组在一起。该步骤使用FindNeighbors和FindClusters两个函数完成：

#计算最邻近距离
pbmc <- FindNeighbors(pbmc, dims = 1:10)
#聚类，包含设置下游聚类的“间隔尺度”的分辨率参数resolution ，增加值会导致更多的聚类。
pbmc <- FindClusters(pbmc, resolution = 0.5)

  可以使用idents函数找到聚类情况：

# 查看前5个细胞的聚类id
head(Idents(pbmc), 5)
## AAACATACAACCAC AAACATTGAGCTAC AAACATTGATCAGC AAACCGTGCTTCCG AAACCGTGTATGCG 
##              1              3              1              2              6 
## Levels: 0 1 2 3 4 5 6 7 8

运行非线性降维

  seurat提供了几种非线性降维技术，例如tsne和umap，来可视化和探索这些数据集。我们将上一步PC维度中的数据进一步降至二维，实现数据的可视化。之前pca是线性降维，能很好处理高维度数据，tsne和umap是非线性降维方式，对低维度数据处理速度快，但难以处理高维度数据。相比tsne，umap能更好地反映原始数据结构，有着更好的连续性。

pbmc <- RunUMAP(pbmc, dims = 1:10)
pbmc <- RunTSNE(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap")
DimPlot(pbmc, reduction = "tsne")

image.png

差异分析

  seurat可以通过差异表达找到每个聚类的markgene，差异分析可以有多种形式，如找到所有聚类的markene（如cluster1中所有的markgene是指cluster1相对于其余所有cluster是差异的）、两个cluster之间的差异分析、某个cluster中两个样品之间差异分析等。

#找到cluster1中的markgene
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
#找到cluster5和cluster0、3之间的markgene
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
#找到每个cluster相比于其余cluster的markgene，只报道阳性的markgene
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

添加细胞注释信息

  如果你对每个cluster已经鉴定出细胞类型，就可以在可视化中将细胞类型标注在图中。

new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono", 
    "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

image.png

可视化

  如果你想画一下可视化图，可以参考官网，里面有很多被用户反馈需求多的图及作图方法，这里列举几个常见的。

#查看特定基因在聚类图中的表达量分布情况
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", 
    "CD8A"))

image.png

#绘制每个cluster按logfc排名前10的markgene的热图
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()

image.png

#查看特定基因在cluster中的表达山脊图
features <- c("LYZ", "CCL5")
RidgePlot(pbmc, features = features, ncol = 2)

image.png

#查看特定基因在cluster中的分布点图，点的大小代表表达该基因的细胞比例，颜色代表平均表达水平
DotPlot(pbmc, features = features) + RotatedAxis()

image.png

---

多个样品分析

  上面介绍的是单个样品Seurat分析的基本过程，但往往我们分析的数据会有多个样品，如果是多个样品，就需要先进行样品数据整合分析。在单个样品的细胞质控完成后，可使用如下过程进行整合：

#合并数据
Merge_data <- merge(J016_12h, y=c(J016_72h, J016_5d, J016_45d)] 
                add.cell.ids = c('J016_12h', 'J016_72h', 'J016_5d', 'J016_45d')], project = "four_merged")
split.list <- SplitObject(Merge_data , split.by = "orig.ident")
#标准化
pancreas.list<-lapply(X=split.list,FUN=function(x){
  x <- NormalizeData(x)
  x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})
#选择参考数据集，如果选择所有样品当作参考数据集，则可以直接使用把pancreas.list当作reference.list
reference.list <- pancreas.list[c('J016_12h', 'J016_72h', 'J016_45d')]
#寻找数据集之间的锚点，利用两个或多个独立单细胞测序数据集中存在的少部分属于同一细胞类型的数据点对整个数据集进行“锚定”进而整合
pancreas.anchors <- FindIntegrationAnchors(object.list = reference.list, dims = 1:30)
#利用锚点整合数据
pancreas.integrated <- IntegrateData(anchorset = pancreas.anchors, dims = 1:30)

  运行Integratedata后，seurat对象将包含一个带有整合表达矩阵的新Assay。注意，原始值(未校正值)仍然存储在“RNA”分析中的对象中，因此您可以来回切换。

# 切换至IntegrateData
DefaultAssay(pancreas.integrated) <- "integrated"
# 运行标准分析
pancreas.integrated <- ScaleData(pancreas.integrated, verbose = FALSE)
pancreas.integrated <- RunPCA(pancreas.integrated, npcs = 30, verbose = FALSE)
pancreas.integrated <- RunTSNE(pancreas.integrated, reduction = "pca", dims = 1:30)
DimPlot(pancreas.integrated, reduction = "tsne", group.by = "orig.ident")

image.png

差异分析

  分析cluster中的markgene需要使用的数据是RNA对象，需要进行Assay切换。

DefaultAssay(pancreas.integrated) <- "RNA"
pancreas.integrated<-NormalizeData(pancreas.integrated, normalization.method = "LogNormalize", scale.factor = 10000)
all.markers <- FindAllMarkers(pancreas.integrated, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

保存与读取数据

#保存中间数据
saveRDS(pancreas.integrated, file = "pancreas.integrated.rds")
#读取中间数据
pancreas.integrated<-readRDS("pancreas.integrated.rds")