【原】cytofWorkflow之聚类分群（四）

健明 2021-07-14

展开全文

前面我们公布了《cytof数据资源介绍（文末有交流群）》，现在就开始正式手把手教学。

上一讲我们构造好了SingleCellExperiment对象，后续全部的分析都会以这个SingleCellExperiment对象为准，大家务必熟悉SingleCellExperiment对象的各种结构，教程见：cytofWorkflow之构建SingleCellExperiment对象（二）。有了这个SingleCellExperiment对象，而且经过了合理的质量控制，接下来就可以进行聚类分群拉！

第一次聚类分群需要指定自己需要多少个亚群

通常建议第一次分群，细一点要好，后面可以人工合并。所以我这里先指定20个亚群，代码如下：


pro='basic_cluster_k20'

set.seed(1234)
# 
# We call ConsensusClusterPlus() with maximum number of clusters maxK = 20.
sce <- cluster(sce, features = "type",
               xdim = 10, ydim = 10, maxK = 20, seed = 1234)

pdf(paste0(pro,'_cluster_plotExprHeatmap_row_clust_F.pdf'))
plotExprHeatmap(sce, features = "type", 
                by = "cluster_id", k = "meta20", 
                row_clust = F,
                bars = TRUE, perc = TRUE)
dev.off()

pdf(paste0(pro,'_plotClusterExprs.pdf'))
plotClusterExprs(sce, k = "meta20", features = "type")
dev.off()

pdf(paste0(pro,'_cluster_plotMultiHeatmap.pdf'))
plotMultiHeatmap(sce, 
                 hm1 = "type",   k = "meta20", 
                 row_anno = FALSE, bars = TRUE, perc = TRUE)
dev.off()

我们直接看最后一个图，这20个细胞亚群如下所示：

初步分成了20个群

理论上不同的亚群肯定是有不一样的抗体信号，比如上图的第3和4亚群，都是高表达CD45,CD3,CD4以及CD7，毫无疑问都是CD4T细胞啦，但是它们也有细微差异，就是CD123啦。假如我们并不关心这么细的分群，就可以把它们都合并为CD4T细胞即可。

然后看看不同细胞亚群的空间分布情况

我们这里采用t-SNE/UMAP 两个方法来可视化，代码如下：


# 节约计算资源
# run t-SNE/UMAP on at most 500/1000 cells per sample
set.seed(1234)
sce <- runDR(sce, "TSNE", cells = 1e3, features = "type")
sce <- runDR(sce, "UMAP", cells = 1e3, features = "type")
# plotDR(sce, "UMAP", color_by = "CD4")

library(ggplot2)
p1 <- plotDR(sce, "TSNE", color_by = "meta20") + 
  theme(legend.position = "none")
p2 <- plotDR(sce, "UMAP", color_by = "meta20")
lgd <- get_legend(p2)
p2 <- p2 + theme(legend.position = "none")
plot_grid(p1, p2, lgd, nrow = 1, rel_widths = c(5, 5, 2))
ggsave2(filename = paste0(pro,'_umap_vs_tSNE.pdf'))

# facet by sample
plotDR(sce, "TSNE", color_by = "meta20", facet_by = "sample_id")
ggsave2(filename = paste0(pro,'_TSNE_by_samples.pdf'))
# facet by condition
plotDR(sce, "TSNE", color_by = "meta20", facet_by = "condition")
ggsave2(filename = paste0(pro,'_TSNE_by_condition.pdf'))


# facet by sample
plotDR(sce, "UMAP", color_by = "meta20", facet_by = "sample_id")
ggsave2(filename = paste0(pro,'_umap_by_samples.pdf'))
# facet by condition
plotDR(sce, "UMAP", color_by = "meta20", facet_by = "condition")
ggsave2(filename = paste0(pro,'_umap_by_condition.pdf'))

plotCodes(sce, k = "meta20")

pdf(paste0(pro,'_cluster_som100_plotMultiHeatmap.pdf'))
plotMultiHeatmap(sce, 
                 hm1 = "type",  k = "som100", m = "meta20", 
                 row_anno = FALSE, col_anno = FALSE, bars = TRUE, perc = TRUE)
dev.off()

sce@metadata
plot_grid(labels = c("A", "B"),
          plotDR(sce, "UMAP", color_by = "meta20"),
          plotDR(sce, "UMAP", color_by = "meta8"))

plotAbundances(sce, k = "meta20", by = "sample_id")
ggsave2(filename = paste0(pro,'_plotAbundances_barplot.pdf'))
plotAbundances(sce, k = "meta20", by = "cluster_id", shape_by = "patient_id")
ggsave2(filename = paste0(pro,'_plotAbundances_boxplot.pdf'))

save(sce,file = 'K20_output_of_cytofWorkflow.Rdata')