Monocle3 基本流程复现

参考：单细胞之轨迹分析-3：monocle3 - 简书 (jianshu.com)，教程中的方法是使用new_cell_data_set()函数创建了一个新的 cds对象，这个函数常用于从头开始创建一个新的rds对象。而本方法使是用as.cell_data_set()函数，直接读入经过Seurat上游处理（SCTransform、RunPCA、RunUMAP、FindNeighbors、FindClusters）后的obj创建rds对象。

• new_cell_data_set ( )：从头创建CDS对象并预处理数据，创建了一个全新的对象，这样很繁琐，还要再做一次降维聚类分群。因为我们导入了seurat对象里的表达矩阵，meta信息和genelist，所以这个cds是没有进行降维聚类等操作，导致后面的拟时序分析是做不了的，官网上说拟时序分析是基于低维度，所以必须对cds进行降维聚类分群。
• as.cell_data_set ( ): 直接读入经过Seurat上游处理（SCTransform、RunPCA、RunUMAP、FindNeighbors、FindClusters）后的obj创建rds对象，不用再进行预处理（preprocess_cds）和降维（reduce_dimension）步骤。

数据：PBMC 3K：Seurat - Guided Clustering Tutorial • Seurat (satijalab.org)

# 安装monocle3
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install(version = "3.10")
BiocManager::install(c('BiocGenerics', 'DelayedArray', 'DelayedMatrixStats',
                       'limma', 'S4Vectors', 'SingleCellExperiment',
                       'SummarizedExperiment', 'batchelor', 'Matrix.utils'))
install.packages("devtools")
devtools::install_github('cole-trapnell-lab/leidenbase')
devtools::install_github('cole-trapnell-lab/monocle3')
-------------------------------------------------------------------------------

library(Seurat)
library(monocle3)
library(tidyverse)
library(patchwork)
library(ggplot2)
library(dplyr)
library(SeuratWrappers)

dir.create("C:/Users/Bolin/Desktop/Monocle3")
setwd("C:/Users/Bolin/Desktop/Monocle3")

pbmc.data <- Read10X("C:/Users/Bolin/Desktop/hg19")
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
obj <- pbmc %>%
  SCTransform(verbose = FALSE, vars.to.regress = 'nCount_RNA') %>% 
  RunPCA(verbose = FALSE, npcs = 20) %>% 
  RunUMAP(verbose = FALSE,dims = 1:6) %>% 
  FindNeighbors(verbose = FALSE) %>% 
  FindClusters(verbose = FALSE)

# new_cell_data_set
-------------------------------------------------------------------------------
data <- GetAssayData(obj, assay = 'RNA', slot = 'counts') # data <- obj@assays$RNA@counts
cell_metadata <- obj@meta.data
gene_annotation <- data.frame(gene_short_name = rownames(data))
rownames(gene_annotation) <- rownames(data)
cds <- new_cell_data_set(data,
                          cell_metadata = cell_metadata,
                          gene_metadata = gene_annotation)
# 预处理：标准化和PCA降维                          
cds <- preprocess_cds(cds, num_dim = 50) # preprocess_cds函数相当于seurat中NormalizeData + ScaleData + RunPCA
plot_pc_variance_explained(cds)

# 可视化
cds <- reduce_dimension(cds,preprocess_method = "PCA") #preprocess_method默认是PCA
plot_cells(cds)                                    
-------------------------------------------------------------------------------

# as.cell_data_set

cds <- as.cell_data_set(obj) # SeuratWrappers这个包提供这个函数，需要从git下载
# 之前已经使用 Seurat 对数据进行了UMAP降维，并使用as.cell_data_set()函数将数据转换为cds 对象。cluster_cells() 函数基于 UMAP 降维结果，在 cds 对象中执行聚类操作
cds <- cluster_cells(cds, reduction_method = "UMAP")
cds <- learn_graph(cds, use_partition = TRUE) #轨迹学习Learn the trajectory graph

p1 <- plot_cells(cds, reduction_method="UMAP", color_cells_by="seurat_clusters") + ggtitle('cds.umap')


-----------------------------------------------------------------------------------
# 一、手动选择root
    cds <- order_cells(cds)

# 二、定义root
    ## 定义函数 get_earliest_principal_node ：根据指定的细胞类型（celltype）名称，找到Seurat对象中对应细胞的主成分图（UMAP）上的最近顶点
    get_earliest_principal_node <- function(cds, time_bin=c('DEGC')){
      # 首先找到指定celltype的 ID
      cell_ids <- which(colData(cds)[, "celltype"] == time_bin)
      # 获取主成分图（UMAP）中细胞投影到最近分支结点的信息
      closest_vertex <- cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
      closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
      # 找到在指定的celltype出现次数最多的分支节点
      root_pr_nodes <- igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names(which.max(table(closest_vertex[cell_ids,]))))] # igraph::V()函数获取主成分图中的所有顶点信息
      root_pr_nodes
    }

nodes_vec <- c(get_earliest_principal_node(cds,"Endothelial cell"),get_earliest_principal_node(cds,"Muscle cell"))
cds = order_cells(cds, root_pr_nodes=nodes_vec,reduction_method = "UMAP")

-----------------------------------------------------------------------------------

p2 = plot_cells(
                cds = cds,
                color_cells_by = "pseudotime",
                show_trajectory_graph = TRUE,
                group_cells_by = 'cluster',cell_size = 1,label_leaves=FALSE) + scale_color_viridis(option = "D")

ggsave(plot=p2,file=paste0('monocle3.pdf'),width=5,height=4)

p2 = plot_cells(cds,
                color_cells_by = "celltype",
                label_groups_by_cluster=FALSE,
                label_leaves=FALSE,
                label_branch_points=FALSE,cell_size = 1, label_cell_groups=FALSE) + ggplot2::scale_color_manual(values=cor)

ggsave(plot=p2,file=paste0('Celltype.pdf'),width=5.5,height=4)



# 差异表达分析：寻找拟时轨迹差异基因
## graph_test分析最重要的结果是莫兰指数（morans_I），其值在-1至1之间，0代表此基因没有
## 空间共表达效应：1代表此基因在空间距离相近的细胞中表达值高度相似。
Track_genes <- graph_test(cds, neighbor_graph="principal_graph", cores=1)
Track_genes <- Track_genes[,c(2,3,4,5,6)] %>% filter(q_value < 1e-3)
write.csv(Track_genes, "Trajectory_genes.csv", row.names = F)


## top10 基因表达趋势图
gene_short_name = rownames(cds)
Track_genes_sig <- Track_genes %>% top_n(n=10, morans_I) %>% pull(gene_short_name) %>% as.character()
p3 <- plot_genes_in_pseudotime(cds[Track_genes_sig,], color_cells_by="seurat_clusters",min_expr=0.5, ncol = 2)
ggsave("Genes_Jitterplot.pdf", plot = p3, width = 8, height = 6)


## FeaturePlot图
p4 <- plot_cells(cds, genes=Track_genes_sig, show_trajectory_graph=T,label_cell_groups=T,label_leaves=F)
p4$facet$params$ncol <- 5
ggsave("Genes_Featureplot.pdf", plot = p4, width = 10, height = 8)


## 寻找共表达基因模块热图
Track_genes <- read.csv("Trajectory_genes.csv")
genelist <- pull(Track_genes, gene_short_name) %>% as.character()
gene_module <- find_gene_modules(cds[genelist,], resolution=1e-1, cores = 1)
write.csv(gene_module, "Genes_Module.csv", row.names = F)
cell_group <- tibble::tibble(cell=row.names(colData(cds)), cell_group=colData(cds)$seurat_clusters)
agg_mat <- aggregate_gene_expression(cds, gene_module, cell_group)
row.names(agg_mat) <- stringr::str_c("Module ", row.names(agg_mat))
p5 <- pheatmap::pheatmap(agg_mat, scale="column", clustering_method="ward.D2")
ggsave("Genes_Module.pdf", plot = p5, width = 10, height = 10)


## 提取拟时分析结果返回seurat对象
pseudotime <- pseudotime(cds, reduction_method = 'UMAP')
pseudotime <- pseudotime[rownames(obj@meta.data)]
obj$pseudotime <- pseudotime
p6 = FeaturePlot(obj, reduction = "umap", features = "pseudotime")

ggsave("Pseudotime_Seurat.pdf", plot = p6, width = 8, height = 6)
saveRDS(pbmc, file = "sco_pseudotime.rds")

pdf("Combined_Plots.pdf")
   print(p5)
   print(p1)
   print(p2)
   print(p3)
   print(p4)
   print(p6)
dev.off()

 

1.寻找拟时轨迹差异基因：

拟时轨迹是指在发育过程中细胞状态随时间（或类似指标）的变化。通过寻找拟时轨迹差异基因，我们可以找到在不同发育阶段或状态中表达水平有显著变化的基因。这些差异基因可能在细胞发育和分化过程中扮演重要的角色，例如调控特定发育途径或细胞类型的标志物。拟时轨迹差异基因的鉴定有助于理解细胞发育过程中的基因调控和细胞状态转换。

2.空间共表达效应：

空间共表达效应是指在空间中距离相近的细胞之间基因表达的相似性。在组织或器官中，相邻的细胞可能具有类似的功能和特性，因此它们的基因表达也可能有相似之处。通过研究空间共表达效应，我们可以揭示细胞在空间上的组织结构和相互作用模式。这有助于了解细胞在组织中的空间分布、细胞-细胞相互作用以及细胞特异性的功能。

通过寻找拟时轨迹差异基因和研究空间共表达效应，我们可以深入了解单细胞数据中的动态细胞状态和空间结构，为细胞发育和组织功能的研究提供重要线索。