seuratV3简介及实操

1字数 3862阅读 1442

Seurat简介

  Seurat---几乎是当前单细胞RNA-seq分析领域的不可或缺的工具,特别是基于10X公司的cellrange流程得出的结果,可以方便的对接到Seurat工具中进行后续处理,简直是带给迷茫在单细胞数据荒漠中小白的一眼清泉,相对全面的功能,简洁的操作命令,如丝般顺滑。

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))
feature-feature relationships
标准化数据

从数据集中删除不需要的单元格后,下一步是标准化数据。默认情况下,采用全局缩放规范化方法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")      
PC

特别是DimHeatmap允许容易地探索数据集中的异质性的主要来源,并且在试图决定包括哪些PC以用于进一步的下游分析时可以是有用的。细胞和特征都根据其PCA分数排序。设置cells为数字可以绘制光谱两端的“极端”单元格,从而大大加快了大型数据集的绘图速度。虽然显然是一种监督分析,但我们发现这是探索相关特征集的有用工具。

DimHeatmap(object = pbmc, dims = 1, cells = 500, balanced = TRUE)
PC热图
确定数据集的“维度”

为了克服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)
JackStrawPlot

另一种启发式方法生成“肘图”:基于每个(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()
定细胞类型
参考材料:

https://satijalab.org/seurat/

公众号首发平台