Monocle 3 | 太牛了！单细胞必学R包！~（五）（差异分析之聚类比较与模块鉴定）

1写在前面

准备出去玩耍了，今天就不废话了，直接上主题吧。🥳

monocle3做差异分析也是牛的一米！~🌾

2用到的包

rm(list = ls())
library(tidyverse)
library(monocle3)

3示例数据

我们还是载入之前用过的一个数据集吧。😘

expression_matrix <- readRDS("./cao_l2_expression.rds")
cell_metadata <- readRDS("./cao_l2_colData.rds")
gene_annotation <- readRDS("./cao_l2_rowData.rds")

还是前面一套老操作，具体的就不讲述了，不清楚的翻看之前的教程吧。🥳

cds <- new_cell_data_set(expression_matrix,
                         cell_metadata = cell_metadata,
                         gene_metadata = gene_annotation)
cds <- preprocess_cds(cds, num_dim = 100)
cds <- reduce_dimension(cds)
cds <- cluster_cells(cds, resolution=1e-5)

colData(cds)$assigned_cell_type <- as.character(partitions(cds))
colData(cds)$assigned_cell_type <- dplyr::recode(colData(cds)$assigned_cell_type,
                                                 "1"="Body wall muscle",
                                                 "2"="Germline",
                                                 "3"="Motor neurons",
                                                 "4"="Seam cells",
                                                 "5"="Sex myoblasts",
                                                 "6"="Socket cells",
                                                 "7"="Marginal_cell",
                                                 "8"="Coelomocyte",
                                                 "9"="Am/PH sheath cells",
                                                 "10"="Ciliated neurons",
                                                 "11"="Intestinal/rectal muscle",
                                                 "12"="Excretory gland",
                                                 "13"="Chemosensory neurons",
                                                 "14"="Interneurons",
                                                 "15"="Unclassified eurons",
                                                 "16"="Ciliated neurons",
                                                 "17"="Pharyngeal gland cells",
                                                 "18"="Unclassified neurons",
                                                 "19"="Chemosensory neurons",
                                                 "20"="Ciliated neurons",
                                                 "21"="Ciliated neurons",
                                                 "22"="Inner labial neuron",
                                                 "23"="Ciliated neurons",
                                                 "24"="Ciliated neurons",
                                                 "25"="Ciliated neurons",
                                                 "26"="Hypodermal cells",
                                                 "27"="Mesodermal cells",
                                                 "28"="Motor neurons",
                                                 "29"="Pharyngeal gland cells",
                                                 "30"="Ciliated neurons",
                                                 "31"="Excretory cells",
                                                 "32"="Amphid neuron",
                                                 "33"="Pharyngeal muscle")
cds

4寻找亚型间的差异基因

首先我们提取一下神经元数据集。😬

提取一下neurons数据集吧。🧐

neurons_cds <- cds[,grepl("neurons", colData(cds)$assigned_cell_type, ignore.case=TRUE)]
plot_cells(neurons_cds, color_cells_by="partition")

这里可以看到神经元有许多亚型，所以也许不同的神经元簇对应着不同的亚型。⚙️

为了研究哪些基因在不同的簇中表达不同，我们可以使用之前介绍的回归分析工具。🧰

当然在这里，Monocle提供了另一种方法来寻找UMAP中不同细胞群之间的差异基因。🧬

函数graph _ test ()使用了一个来自空间自相关分析的统计数据，称为Moran’s I。😘

pr_graph_test_res <- graph_test(neurons_cds, neighbor_graph="knn", cores=8)

pr_deg_ids <- row.names(subset(pr_graph_test_res, q_value < 0.05))

head(pr_deg_ids)

5寻找共调控基因模块

非常简单，调用find_gene_modules即可。🥸

gene_module_df <- find_gene_modules(neurons_cds[pr_deg_ids,], resolution=1e-2)

1️⃣ 可视化方案一：👇

这里需要计算一个module内的基因表达的总和，需要用到aggregate_gene_expression。🥰

cell_group_df <- tibble::tibble(cell=row.names(colData(neurons_cds)), 
                                cell_group=partitions(cds)[colnames(neurons_cds)])
agg_mat <- aggregate_gene_expression(neurons_cds, gene_module_df, cell_group_df)
row.names(agg_mat) <- stringr::str_c("Module ", row.names(agg_mat))
colnames(agg_mat) <- stringr::str_c("Partition ", colnames(agg_mat))

pheatmap::pheatmap(agg_mat, 
                   cluster_rows=T, 
                   cluster_cols=T,
                   scale="column", 
                   clustering_method="ward.D2",
                   fontsize=6)

2️⃣ 可视化方案二：👇

第二种方案就是使用plot_cells进行可视化。🥰

如果有很多模块，很难看出每个module在哪里表达，所以我们只看它们的一个子集。🤒

plot_cells(neurons_cds, 
           genes=gene_module_df %>% filter(module %in% c(8, 28, 33, 37)),
           group_cells_by="partition",
           color_cells_by="partition",
           show_trajectory_graph=F)