這里使用官方提供的示例數(shù)據(jù)進行繼續(xù)講解，首先載入所需要的數(shù)據(jù)：

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）數(shù)據(jù)預處理

這一步的目的主要是對數(shù)據(jù)標準化和去除批次效應。

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

如果你不確定num_dim填多少，可以先畫個圖看看，保證所使用的主成分數(shù)目足以捕捉數(shù)據(jù)集中所有細胞基因表達的主要變異：

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）去除批次效應

如果你的數(shù)據(jù)中包含多個批次，需要先去除批次效應：

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)