monocle3系列之三：聚类

这里使用官方提供的示例数据进行继续讲解，首先载入所需要的数据：

library(monocle3)
library(dplyr)

expression_matrix <- readRDS(url("https://depts.washington.edu:/trapnell-lab/software/monocle3/celegans/data/cao_l2_expression.rds"))
cell_metadata <- readRDS(url("https://depts.washington.edu:/trapnell-lab/software/monocle3/celegans/data/cao_l2_colData.rds"))
gene_annotation <- readRDS(url("https://depts.washington.edu:/trapnell-lab/software/monocle3/celegans/data/cao_l2_rowData.rds"))

cds <- new_cell_data_set(expression_matrix,
                         cell_metadata = cell_metadata,
                         gene_metadata = gene_annotation)

1）数据预处理

这一步的目的主要是对数据标准化和去除批次效应。

## num_dim是指定主成分的数目
cds <- preprocess_cds(cds, num_dim = 100)

如果你不确定num_dim填多少，可以先画个图看看，保证所使用的主成分数目足以捕捉数据集中所有细胞基因表达的主要变异：

plot_pc_variance_explained(cds)

2）降维

可以使用tSNE或者UMAP（默认）

cds <- reduce_dimension(cds)
plot_cells(cds)
## plot_cells(cds, color_cells_by="cao_cell_type")
## plot_cells(cds, genes=c("cpna-2", "egl-21", "ram-2", "inos-1"))

如果想用tSNE降维，需要明确指定：

cds <- reduce_dimension(cds, reduction_method="tSNE")
plot_cells(cds, reduction_method="tSNE", color_cells_by="cao_cell_type")

3）去除批次效应

如果你的数据中包含多个批次，需要先去除批次效应：

cds <- align_cds(cds, num_dim = 100, alignment_group = "plate")
cds <- reduce_dimension(cds)
plot_cells(cds, color_cells_by="plate", label_cell_groups=FALSE)

4）聚类

cds <- cluster_cells(cds, resolution=1e-5)
plot_cells(cds)
plot_cells(cds, color_cells_by="partition", group_cells_by="partition") ## partition是 PAGA algorithm中比cluster更大一些的分组

5）找每个cluster的差异基因

We could group the cells according to cluster, partition, or any categorical variable in colData(cds).

marker_test_res <- top_markers(cds, group_cells_by="partition", 
                               reference_cells=1000, cores=8)

通过pseudo_R2筛选特异性marker。

top_specific_markers <- marker_test_res %>%
                            filter(fraction_expressing >= 0.10) %>%
                            group_by(cell_group) %>%
                            top_n(1, pseudo_R2)

top_specific_marker_ids <- unique(top_specific_markers %>% pull(gene_id))

将上述特异性marker可视化：

plot_genes_by_group(cds,
                    top_specific_marker_ids,
                    group_cells_by="partition",
                    ordering_type="maximal_on_diag",
                    max.size=3)