细胞亚群的生物学命名
作者 | 单细胞天地小编 刘小泽
课程链接在:http://jm.grazy.cn/index/mulitcourse/detail.html?cid=55
第二单元第8讲:细胞亚群的生物学命名
上一次我们好不容易得到了这个1.9G的RData,也就是作者自己做出来的Seurat对象,那么我们要怎么利用它去进一步探索呢?
加载上次运行的结果
首先还是三步走
1rm(list = ls())
2options(warn=-1)
3suppressMessages(library(Seurat))
然后加载进来之前保存的5.1G PBMC对象
1start_time <- Sys.time()
2load('./patient1.PBMC.output.Rdata')
3end_time <- Sys.time()
4end_time - start_time
5# Time difference of 12.6741 secs
重温上次的聚类结果
1# 使用Seurat 2.3.4版本
2colP<-c('green4',
3 'pink',
4 '#FF7F00',
5 'orchid',
6 '#99c9fb',
7 'dodgerblue2',
8 'grey30',
9 'yellow',
10 'grey60',
11 'grey',
12 'red',
13 '#FB9A99',
14 'black'
15)
16TSNEPlot(PBMC,
17 colors.use = colP,
18 do.label = T)
开始新的探索
看看作者整合的数据对批次的处理
也就是把四个时间点映射到上面的tsne坐标中,并且理论上应该是:每群细胞都覆盖到四个时间点
1TSNEPlot(PBMC,group.by = "TimePoints")
再用table
对比一下
1> table(PBMC@meta.data$TimePoints,PBMC@ident)
2
3 0 1 2 3 4 5 6 7 8 9 10 11 12
4 PBMC_ARD614 665 726 572 559 420 302 457 313 283 123 17 11 68
5 PBMC_EarlyD27 43 173 245 85 120 110 543 59 91 29 7 3 84
6 PBMC_Pre 369 527 197 93 146 393 4 76 17 48 25 187 0
7 PBMC_RespD376 800 433 555 677 636 516 119 324 204 200 170 11 39
可视化一些marker基因
这些marker基因也是来源于文章的Supp Fig.7,基于他们对免疫知识的了解
1allGenes = row.names(PBMC@raw.data)
2markerGenes <- c(
3 "CD3D",
4 "CD3E",
5 "TRAC",
6 "IL7R",
7 "GZMA",
8 "FCGR3A",
9 "CD14",
10 "MS4A1",
11 "FCER1A"
12)
13# 判断这些marker是不是存在表达矩阵中
14> markerGenes %in% allGenes
15[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
16# 选择性保存pdf格式
17# pdf('patient1_pBMC_marker_FeaturePlot.pdf', width=10, height=15)
18FeaturePlot(object = PBMC,
19 features.plot =markerGenes,
20 cols.use = c("grey", "blue"),
21 reduction.use = "tsne")
利用上面marker基因在不同细胞群的特殊表达,其实就能得到每个群的细胞命名,不过这些marker基因的获得以及和细胞名称的对应,是需要一段时间的研究才能得到的。我们这里只是展示如何操作
赋予每个cluster细胞类型
根据我们得到的cluster和原文的命名
就可以对应得到一个细胞名称列表:
1cat >celltype-patient1-PBMC.txt
20 "B cells"
31 "CD4+ T cells"
42 "Naive memory T cells"
53 "Classical monicytes"
64 "CD8+ effector T cells"
75 "NK cells"
86 "rm1"
97 "Non-classical monocytes"
108 "Dendritic cells"
119 "rm2"
1210 "CD8+ cytotoxic T cells"
1311 "Myeloid cells"
1412 "rm3"
假如我们是先有了这个列表,核心就是要将分群的clsuter数字与细胞名称和颜色对应起来
利用Seurat V3
版本3的优势是可以直接提取出seurat对象的cluster数字,并且提供了RenameIdents
函数,直接进行转换
1a=read.table('celltype-patient1-PBMC.txt')
2new.cluster.ids <- as.character(a[,2])
3names(new.cluster.ids) <- levels(PBMC_V3)
4PBMC_V3 <- RenameIdents(PBMC_V3, new.cluster.ids)
5
6DimPlot(PBMC_V3, reduction = "tsne", label = TRUE, pt.size = 0.5, cols = colP) + NoLegend()
利用Seurat V2
版本2就需要自己去对应:首先要了解它的分群信息存储在PBMC@ident
中,然后要将全部12874个细胞的分群编号与celltype-patient1-PBMC.txt
表中的第一列编号对应(只有这样,才能和表中的第二列对应上)
用到对应关系时,首先思考
match
能不能做到;如果要做,需要准备什么;对应关系搞清楚
1# 先与表中第一列对应
2match(as.numeric(as.character(PBMC@ident)),a[,1])
3
4# 第一点需要注意的是:为什么先用as.character后用as.numeric,而不是直接用as.numeric?
5# 原因就是PBMC@ident存储的是因子型变量,直接取只会得到它们的位置信息,而不是真实的分群信息
6> head(as.numeric(PBMC@ident))
7[1] 2 1 2 2 2 6
8> head(as.character(PBMC@ident))
9[1] "1" "0" "1" "1" "1" "5"
10> head(as.numeric(as.character(PBMC@ident)))
11[1] 1 0 1 1 1 5
12
13# 第二点需要注意的是:match函数的规则是,A要在B中找到对应位置,那么就是 match(A,B)
接下来就可以得到对应的第二列,也就是细胞名称
1labels=a[match(as.numeric(as.character(PBMC@ident)),a[,1]),2]
2
3# 检查一下,主要看数量的对应
4> table(labels)
5labels
6 B cells CD4+ T cells CD8+ cytotoxic T cells
7 1877 1859 219
8 CD8+ effector T cells Classical monicytes Dendritic cells
9 1322 1414 595
10 Myeloid cells NK cells Naive memory T cells
11 212 1321 1569
12Non-classical monocytes rm1 rm2
13 772 1123 400
14 rm3
15 191
16> table(PBMC@ident)
17
18 0 1 2 3 4 5 6 7 8 9 10 11 12
191877 1859 1569 1414 1322 1321 1123 772 595 400 219 212 191
添加到metadata中,方便后面使用
1PBMC@meta.data$labels=labels
2
3TSNEPlot(PBMC, group.by = 'labels',
4 colors.use = colP,
5 do.label = T)
但很明显,颜色标记和原文不同。因为之前只是对应了分群的编号和细胞名称,此外还需要修改颜色的顺序
1# 修改颜色顺序,思想就是:将现在的细胞名与表中的细胞名对应一下,然后这个顺序就是颜色出现的顺序
2colP=colP[match(levels(as.factor(labels)),a[,2])]
3TSNEPlot(PBMC, group.by = 'labels',
4 colors.use = colP,
5 do.label = T)
再按时间拆分分群结果
做出文章的这张图
首先得到我们的四个时间点
1> TimePoints = PBMC@meta.data$TimePoints
2> table(TimePoints)
3TimePoints
4 PBMC_ARD614 PBMC_EarlyD27 PBMC_Pre PBMC_RespD376
5 4516 1592 2082 4684
然后用一个函数`SubsetData`
举一个例子,绘制Pre
时期的分群图
1# SubsetData函数其实也是利用TRUE/FALSE取子集
2PBMC_Pre = SubsetData(PBMC,TimePoints =='PBMC_Pre')
3TSNEPlot(PBMC_Pre,
4 colors.use = c('green4', 'pink', '#FF7F00', 'orchid', '#99c9fb', 'dodgerblue2', 'grey30', 'yellow', 'grey60', 'grey', 'red', '#FB9A99', 'black'),
5 do.label = F)
6# ggsave('PBMC_Pre_tSNE.pdf')