Seurat学习与使用(一)

简介

  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")
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 151,511评论 1 330
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 64,495评论 1 273
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 101,595评论 0 225
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 42,558评论 0 190
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 50,715评论 3 270
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 39,672评论 1 192
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,112评论 2 291
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 29,837评论 0 181
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 33,417评论 0 228
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 29,928评论 2 232
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,316评论 1 242
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 27,773评论 2 234
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,253评论 3 220
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 25,827评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,440评论 0 180
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 34,523评论 2 249
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 34,583评论 2 249

推荐阅读更多精彩内容