查看: 1893|回复: 3

[其他] 单细胞分组数据分析学习笔记分享

[复制链接]

迅猛龙

Rank: 8Rank: 8

主题
166
注册时间
2020.6.16
在线时间
99 小时

发表于 2020.11.2 11:26:29 | 显示全部楼层 |阅读模式
本帖最后由 基迪奥-Jt桃 于 2020.11.2 11:25 编辑

对于分组数据来说,由于处理效应的存在,处理组和对照组的基因表达水平是存在差异的,因此直接合并(merge)两组细胞后进行细胞分群可能会存在问题。但Seurat提供了一套整合数据的方法,可“消除”处理效应对细胞分群的影响,具体算法可阅读参考资料中相应的文献,原理如下:

(Cell,2019)

本教程的主要内容包括如何使用整合后的数据进行细胞分群,找到对照组和处理组中保守的细胞标记和差异基因。

数据准备

[AppleScript] 纯文本查看 复制代码
#载入Seurat包;
library(Seurat)
#读入表达矩阵;
ctrl.data <- read.table(file = "../data/immune_control_expression_matrix.txt", sep = "\t")
stim.data <- read.table(file = "../data/immune_stimulated_expression_matrix.txt", sep = "\t")
#创建对照组数据Seurat对象;
ctrl <- CreateSeuratObject(counts = ctrl.data, project = "IMMUNE_CTRL", min.cells = 5)
ctrl$stim <- "CTRL"
#创建处理组数据Seurat对象;
stim <- CreateSeuratObject(counts = stim.data, project = "IMMUNE_STIM", min.cells = 5)
stim$stim <- "STIM"
#将两个对象合并到同一个list中;
ifnb.list = list(ctrl, stim)


以上是从头创建数据对象的方法,在最新的Seurat官方教程中,范例数据是储存在SeuratData包中的,也就是说,我们可以直接安装SeuratData包,并调用其中的范例数据。

Tips:

以上创建Seurat对象的方法仅作为本文的补充,接下来主要用SeuratData包中对应的数据做演示。

[AppleScript] 纯文本查看 复制代码
#SeuratData包主流安装方法如下:
devtools::install_github('satijalab/seurat-data')


但是,这种安装方法太慢了,很难安装成功;可尝试将安装包使用迅雷下载到本地,复制到R软件安装目录(C:\Program Files\R\R-4.0.2\bin\x64),对于Windows用户可打开cmd,通过命令进行本地安装SeuratData包:

rcmd INSTALL ifnb.SeuratData_3.1.0.tar.gz

[AppleScript] 纯文本查看 复制代码
#载入R包;
library(SeuratData)
#查看SeuratDat包含的数据集;
AvailableData()
#载入本次的范例数据;
data(ifnb)
#查看范例数据说明信息;
?ifnb
数据对象的结构大致如下:


[AppleScript] 纯文本查看 复制代码
#在本范例数据中,将PBMCs分成受干扰素β刺激组和对照组,这里按照分组信息(stim或orig.ident)将案例数据merge后的对象拆分成两个对象,储存到list对象中;
ifnb.list <- SplitObject(ifnb, split.by = "stim")

#自定义函数func,用于对数据进行标准化并计算高变基因;
func <- function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
}

#对列表中的两个对象使用自定义函数进行标准化和高变基因计算;
ifnb.list <- lapply(X = ifnb.list, FUN = func)


Tips:

lapply()函数:对列表中的对象逐个应用函数,返回结果仍为列表;lapply()可视作代替循环的一种方法。

整合数据

[AppleScript] 纯文本查看 复制代码
#以Seurat objects作为input,用FindIntegrationAnchors /ˈæŋkər/ 函数识别锚点(耗时2~5min),并使用这些锚点将两个数据集整合到在一起,这里用到的是IntegrateData函数。
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, dims = 1:20)
immune.combined <- IntegrateData(anchorset = immune.anchors, dims = 1:20)
head(immune.anchors@anchors)



整合后的数据在assays 的slots中添加了integrated数据对象,如下图。


[AppleScript] 纯文本查看 复制代码
#保存整合后的临时对象,便于之后反复调用;
save(immune.anchors,file = "immune.anchors.Rdata")
save(immune.combined,file = "immune.combined.Rdata")


细胞分群

接下来可以对所有的细胞进行常规的分群流程啦!

[AppleScript] 纯文本查看 复制代码
#将默认的assay(active.assay)设为"integrated";
DefaultAssay(immune.combined) <- "integrated"

#对数据进行归一化和PCA降维;
immune.combined <- ScaleData(immune.combined, verbose = FALSE)
immune.combined <- RunPCA(immune.combined, npcs = 30, verbose = FALSE)
#对细胞进行聚类、分群、UMAP可视化;
immune.combined <- RunUMAP(immune.combined, reduction = "pca", dims = 1:20)
immune.combined <- FindNeighbors(immune.combined, reduction = "pca", dims = 1:20)
immune.combined <- FindClusters(immune.combined, resolution = 0.5)

library(cowplot)
library(patchwork)
p1 <- DimPlot(immune.combined, reduction = "umap", group.by = "stim")
p2 <- DimPlot(immune.combined, reduction = "umap", label = TRUE)
plot_grid(p1, p2)



[AppleScript] 纯文本查看 复制代码
#也可以对两个分组单独绘制;
DimPlot(immune.combined, reduction = "umap", split.by = "stim")



鉴定保守细胞标记

[AppleScript] 纯文本查看 复制代码
#将默认的Assay(active.assay)设为"RNA";
DefaultAssay(immune.combined) <- "RNA"
#差异分析(即查找保守makergene,耗时1~2min);这里计算cluster 6 (NK细胞)中与刺激条件无关的保守标记基因。
nk.markers <- FindConservedMarkers(immune.combined, ident.1 = 6, grouping.var = "stim", verbose = FALSE)
head(nk.markers)



[AppleScript] 纯文本查看 复制代码
#绘制保守marker genes的umap映射图;
#多图比较耗时(1~2min);
FeaturePlot(immune.combined, features = c("CD3D", "SELL", "CREM", "CD8A"), min.cutoff = "q9")



[AppleScript] 纯文本查看 复制代码
#对细胞亚群名称(active.ident)进行重命名;
immune.combined <- RenameIdents(immune.combined,
`0` = "CD14 Mono", `1` = "CD4 Naive T", `2` = "CD4 Memory T",
`3` = "CD16 Mono", `4` = "B", `5` = "CD8 T", `6` = "NK", `7` = "T activated",
`8` = "DC", `9` = "B Activated", `10` = "Mk", `11` = "pDC", `12` = "Eryth")
#重绘umap分群图;
DimPlot(immune.combined, label = TRUE,pt.size=0.5,label.size = 5)



组间差异分析

接下来鉴定组间差异基因,这里用处理组和对照的naive T cells 、CD14 单核细胞的平均表达量绘制散点图,突出显示对干扰素刺激表现出显著应答的基因。

[AppleScript] 纯文本查看 复制代码
theme_set(theme_cowplot())
#提取CD4 Naive T细胞并计算平均表达量;
t.cells <- subset(immune.combined, idents = "CD4 Naive T")
#将stim列作为默认id(active.ident)列;
Idents(t.cells) <- "stim"
avg.t.cells <- log1p(AverageExpression(t.cells, verbose = FALSE)$RNA)
avg.t.cells$gene <- rownames(avg.t.cells)
#同样,提取CCD14 Mono细胞并计算平均表达量;
cd14.mono <- subset(immune.combined, idents = "CD14 Mono")
Idents(cd14.mono) <- "stim"
avg.cd14.mono <- log1p(AverageExpression(cd14.mono, verbose = FALSE)$RNA)
avg.cd14.mono$gene <- rownames(avg.cd14.mono)

library(ggplot2)
genes.to.label = c("ISG15", "LY6E", "IFI6", "ISG20", "MX1", "IFIT2", "IFIT1", "CXCL10", "CCL8")
p1 <- ggplot(avg.t.cells, aes(CTRL, STIM)) + geom_point(color="#a8df65") + ggtitle("CD4 Naive T Cells")
p1 <- LabelPoints(plot = p1, points = genes.to.label, repel = TRUE)
p1



[AppleScript] 纯文本查看 复制代码
p2 <- ggplot(avg.cd14.mono, aes(CTRL, STIM)) + geom_point( color="#a8df65") + ggtitle("CD14 Monocytes")
p2 <- LabelPoints(plot = p2, points = genes.to.label, repel = TRUE)
p2



从上图可知,许多基因在这两种细胞类型中都上调表达,可能存在一种保守的干扰素应答通路。接下来,我们可以分析在相同类型的细胞中,在干扰素处理下会有哪些基因的表达发生变化。

[AppleScript] 纯文本查看 复制代码
#首先,在immune.combined对象的meta.data数据框中添加两列:celltype.stim和celltype;
immune.combined$celltype.stim <- paste(Idents(immune.combined), immune.combined$stim, sep = "_")
immune.combined$celltype <- Idents(immune.combined)
#更改immune.combined对象的id(active.ident)为celltype.stim列;
Idents(immune.combined) <- "celltype.stim"
#然后,我们使用FindMarkers函数查找B细胞中处理组和对照组之间的差异基因;
b.interferon.response <- FindMarkers(immune.combined, ident.1 = "B_STIM", ident.2 = "B_CTRL", verbose = FALSE)
#查看前15个差异基因;
head(b.interferon.response, n = 15)



[AppleScript] 纯文本查看 复制代码
#接下来使用FeaturePlot或VlnPlot函数的split.by参数,按分组对特定差异基因进行可视化;
FeaturePlot(immune.combined, features = c("CD3D", "GNLY", "IFI6"),
split.by = "stim", max.cutoff = 3, cols = c("grey70", "red"))



从上面的结果可以看出,因为CD3D和GNLY基因是典型的细胞类型标记(T细胞和NK/CD8 T细胞),实际上不受干扰素的影响,在对照组和处理组中显示相似的基因表达模式。而IFI6是核心干扰素应答基因,几乎在所有细胞类型中都上调表达。

[AppleScript] 纯文本查看 复制代码
#使用感兴趣的差异基因绘制小提琴图;
plots <- VlnPlot(immune.combined, features = c("LYZ", "ISG15", "CXCL10"), split.by = "stim",group.by = "celltype", pt.size = 0, combine = FALSE)
wrap_plots(plots = plots, ncol = 1)



从上图的结果可以看出,在干扰素刺激后,CXCL10在单核细胞和B细胞中表现出明显的上调,而基因ISG15几乎在所有细胞类型中都上调表达。

好啦,今天的单细胞数据分析学习笔记就分享到这里啦!

参考资料:
Stuart T, Butler A, Hoffman P, et al. Comprehensive integration of single-cell data[J]. Cell, 2019, 177(7): 1888-1902. e21.
https://satijalab.org/seurat/vignettes.html

本文作者:基迪奥-莫北         


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
新的一天加油!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
1
注册时间
2019.11.26
在线时间
2 小时

发表于 2020.11.2 19:08:38 | 显示全部楼层
哈哈哈哈哈哈哈哈哈
回复 支持 反对

使用道具 举报

帝王蝶

Rank: 4

主题
0
注册时间
2020.2.12
在线时间
35 小时

发表于 2020.11.2 19:11:36 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
1
注册时间
2016.9.4
在线时间
138 小时

发表于 2020.11.3 08:26:33 | 显示全部楼层
今日已签到
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表