通常情况下,单细胞转录组拿到亚群后会进行更细致的分群,或者看不同样本不同组别的内部的细胞亚群的比例变化。
这就是个性化分析阶段,这个阶段取决于自己的单细胞转录组项目课题设计情况,我们了解到的各式各样的分析点,并不是通用的。比如如果要比较细胞亚群比例,就必须要有多个样本,如果是单个样本,其分析的内容及方法也会不尽相同,类似的问题以及涉及到的分析很多。
今天小编带大家来学习下“单细胞个性化分析之细胞亚群继续分群”,相信大多数单细胞转录组的小伙伴都能用到。
细胞亚群进一步划分原则
理论上细胞亚群是可以无限划分的,因为世界上没有两个一模一样的细胞,关键是要把握一个度,什么样的差异可以判定为不同细胞亚群,什么样的差异是可以容忍的细胞类群内部异质性。有一个策略就是找出主要因素和次要因素。主要因素划分为主要亚群,比如外周血里面的T,B细胞当然是不同亚群,但是T细胞里面还可以继续划分:CD4或者CD8的T细胞,甚至继续划分, 如下图所示:
亚群细化分析实操
加载R包
1.{
2.library(Seurat)
3.library(dplyr)
4.library(reticulate)
5.library(sctransform)
6.library(cowplot)
7.library(ggplot2)
8.library(viridis)
9.library(tidyr)
10.library(magrittr)
11.library(reshape2)
12.library(readxl)
13.library(stringr)
14.library(cowplot)
15.library(scales)
16.library(readr)
17.library(progeny)
18.library(gplots)
19.library(tibble)
20.library(grid)
21.library(rlang)
22.
23.theme_set(theme_cowplot())
24.use_colors <- c(
25. Alveolar_Macrophages = "#6bAEd6",
26. `Monocyte-derived macrophages`= "#fff500",
27. Monocytes= "#FA9FB5",
28. `Myeloid dendritic cells`= "#DD3497",
29. Plasmacytoid_dendritic_cells= "#7A0177",
30. `Conventional_T_cells`= "#c2e699",
31. `Regulatory T cells`= "#006837",
32. `CD8+_T_cells`= "#bcbddc",
33. NK_cells= "#4a1486",
34. B_cells= "#969696",
35. Plasma_cells= "#636363")
36.}
37.imm_anno <- readRDS("imm_anno.RDS")
38.imm_anno@meta.data$cell_type_imm <-
39. ordered(
40. imm_anno@meta.data$cell_type_imm, levels = c("Alveolar_Macrophages",
41. "Monocyte-derived macrophages",
42. "Monocytes","Myeloid dendritic cells",
43. "Plasmacytoid_dendritic_cells",
44. "Conventional_T_cells",
45. "Regulatory T cells",
46. "CD8+_T_cells",
47. "NK_cells",
48. "B_cells",
49. "Plasma_cells")
50. )
51.DimPlot(imm_anno, reduction = 'umap', split.by = 'tissue_type', label = T, cols = use_colors)
52.##产生聚类点图
53.
细化分析的开端是发现(感觉)某群细胞有多个亚群,如下面的图发现单核来源巨噬有多群。
1.### 实际上同一群细胞内部小亚群比较更有意思
2.###导入M1
3.m1m2_pws <- read_lines("./data/CLASSICAL_M1_VS_ALTERNATIVE_M2_MACROPHAGE_UP.gmt") %>%
4. lapply(str_split, "\\t") %>%
5. unlist(recursive = F) %>%
6. lapply(function(x) setNames(list(x[-c(1:2)]), x[1])) %>%
7. unlist(recursive = F)
8.
9..### 导入M2
10.m1m2_pws <- append(m1m2_pws, read_lines("./data/CLASSICAL_M1_VS_ALTERNATIVE_M2_MACROPHAGE_DN.gmt") %>%
11. lapply(str_split, "\\t") %>%
12. unlist(recursive = F) %>%
13. lapply(function(x) setNames(list(x[-c(1:2)]), x[1])) %>%
14. unlist(recursive = F))
15.
16.### 运用addmodulescore函数,加入富集分数到metadata中
17.imm_anno <- AddModuleScore(object = imm_anno, features = m1m2_pws, name = c("m1up", "m1dn"), nbin = 12)### 提出单核来源巨噬来分析macrophage(mdm)
18.scRNA_mdm <- subset(imm_anno, subset = cell_type_imm %in% c("Monocyte-derived macrophages"))
19.scRNA_mdm=ScaleData(scRNA_mdm)
20.
21.### 重新降维聚类
22.scRNA_mdm <- RunPCA(scRNA_mdm)
23.ElbowPlot(scRNA_mdm, ndims = 50)
24.
25.library(ggplot2)
26.scRNA_mdm<- RunUMAP(scRNA_mdm, dims = 1:20)
27.scRNA_mdm <- FindNeighbors(scRNA_mdm, dims = 1:20)
28.for (i in c(0.01,0.2, 0.3, 0.4, 0.5, 1)) {
29. scRNA_mdm <- FindClusters(scRNA_mdm, resolution = i)
30. print(DimPlot(scRNA_mdm, reduction = "umap") + labs(title = paste0("resolution: ", i)))
31.}
32.### cluster数自己决定,不要太多即可,否则相似的群多
33.table(scRNA_mdm$SCT_snn_res.0.2)
34.Idents(scRNA_mdm) <-scRNA_mdm@meta.data$SCT_snn_res.0.2
35.
36.DimPlot(scRNA_mdm,reduction = 'umap',group.by = 'SCT_snn_res.0.2',label = T,split.by = 'tissue_type')
37.ggsave2(filename = 'Macrophage亚群.pdf',path = './result_3v3/analysis',width = 8,height = 4)
1.### 不同mdm小亚群比较通路
2.VlnPlot(scRNA_mdm, features = c("m1up1", "m1dn2"),
3. group.by = "SCT_snn_res.0.2", pt.size = 0,
4. cols = use_colors)
5.### 第四群有一个显著的M1活化
6.ggsave2(filename = 'M1M2巨噬细胞细分.pdf',path = './result_3v3/analysis/',width = 8,height =4)
1.### 单核细胞来源的巨噬细胞间具有较大异质性,所以考虑再分别注释
2.
3.mdm_markers <- FindAllMarkers(scRNA_mdm, only.pos = T, min.pct = 0.25, min.diff.pct = 0.25)
4.
5.# group_by和top_n组合
6.top_mdm_markers <- mdm_markers %>% group_by(cluster) %>% top_n(10, wt = avg_log2FC)
7.
8.DoHeatmap(scRNA_mdm,features=top_mdm_markers$gene,group.by = 'SCT_snn_res.0.2')+scale_fill_viridis()
9.ggsave2(filename = '热图细分.pdf',path = './result_3v3/analysis/',width = 10,height =9)
1.#### 第4群CCL5+,其实还有CD8A+,大家认为,这是一群新的巨噬,还是由于细胞污染呢~
2.FeaturePlot(scRNA_mdm,features = 'CCL5',cols = viridis(10),label = T)
3.ggsave2(filename = 'CCL5_mdm.pdf',path = './result_3v3/analysis/',width = 5,height =4.5)
4.
5.### 据此我们可以继续定义细胞亚群,例如cluster1为C1Q+的巨噬细胞!
亚群可视化
亚群细分好后,通常我们期望能直观展示分析结果,下面给大家展示2种常见的可视化方式。
某些文章里面会把主要和次要细胞亚群同一个tSNE图展现,实际上,细胞二维散点图,是没办法写全部细胞亚群的生物学功能定义的, 通常也是把主要细胞亚群标记上去。然后每个细胞亚群再进行继续分群。只不过是在同一个散点图上面展示(如下图所示)。
但上面这种展示形式展示的细胞数量太多,一般人很难看出来不同细胞亚群内部具体如何划分子亚群,以及不同亚群到底以什么标记基因来进行区分。有一个策略是,把标记基因的表达量在所有亚群及子亚群里面展现,以热图形式,如下图: