Monocle3

一、monocle3介绍

monocle3可以执行的3个主要的功能

对细胞进行聚类、分类和计数。Monocle 3 可识别新的（可能是罕见的）亚型细胞

构建单细胞轨迹。通过拟时序分析帮助大家解析生物体发育、疾病等过程中细胞发生的变化。这是最主要的功能。

差异表达分析。Monocle 3 包括一个复杂但易于使用的差分表达系统，可以表征新的细胞类型和状态始于与其他更好理解的细胞进行比较。

 主要流程

1.读取数据，创建cell_data_set对象。

2.数据预处理：标准化，去除批次效应

3.降维

4.聚类

5.进行差异基因表达分析

6.拟时序分析

二、关于不同数据的读取办法

（一）Bioconductor的ExpressionSet对象：

monocle3读取的数据要包含3个部分：

• expression_matrix：表达式值的数字矩阵，其中行是gene，列是cell
• cell_metadata:数据框，其中行是cell，列是细胞属性（如细胞类型、培养条件、获取日期等）
• gene_metadata :数据框，其中行是特征（例如基因），列是基因属性，例如生物型，GC含量等。

输入前应该确保几个“等式”

•  expression_matrix 列数= cell_metadata 行数 并且两者要相匹配
•  expression_matrix 行数= gene_metadata 列数 并且两者要相匹配
•  gene_metadata的其中一列应为“gene_short_name”,表示每个基因的基因符号或简要名称

创建对象示例 

# Load the data



expression_matrix <- readRDS(url(":/trapnell-lab/software/monocle3/celegans/data/cao_l2_expression.rds"))




cell_metadata <- readRDS(url(":/trapnell-lab/software/monocle3/celegans/data/cao_l2_colData.rds"))




gene_annotation <- readRDS(url(":/trapnell-lab/software/monocle3/celegans/data/cao_l2_rowData.rds"))




 



 



# Make the CDS object



cds <- new_cell_data_set(expression_matrix,




                         cell_metadata = cell_metadata,




                         gene_metadata = gene_annotation)

（二）从 10X 输出生成cell_data_set

要找到正确的文件，必须提供包含未修改的Cell Ranger“outs”文件夹的文件夹的路径。文件结构应如下所示：10x_data/outs/filtered_feature_bc_matrix/，其中filtered_feature_bc_matrix包含 files 、 和 。（还可以处理Cell Ranger V2数据，其中“feature”被替换为“gene”，并且文件未被gz压缩。(这个读取起来很容易报错)

读取示例:注意使用函数是load_cellranger_data

# Provide the path to the Cell Ranger output.



cds <- load_cellranger_data("~/Downloads/10x_data")

或者如果有以下3个文件，也可直接读取：注意使用函数是load_mm_data

,                    和

cds <- load_mm_data(mat_path = "~/Downloads/", 



                    feature_anno_path = "~/Downloads/", 



                    cell_anno_path = "~/Downloads/")

（三）处理大量数据

当细胞量或基因量比较大的话，可以使用稀疏矩阵的方法读取 

cds <- new_cell_data_set(as(umi_matrix, "sparseMatrix"),




cell_metadata = cell_metadata,




gene_metadata = gene_metadata)

（四）多个cds对象组合

# make a fake second cds object for demonstration



cds2 <- cds[1:100,]




 



big_cds <- combine_cds(list(cds, cds2))

二、数据预处理

（一）数据标准化

cds <- preprocess_cds(cds, num_dim = 100)#主要是完成数据的标准化＋PCA




plot_pc_variance_explained(cds)#查看的是每个维度的变异分数




cds <- align_cds(cds, alignment_group = "plate")#去除批次效应

查看每个PC的变异分数：plot_pc_variance_explained() 

 （二）降维分析并可视化

   Monocle 3默认使用UMAP降维，降维方法有两种可以选择：PCA和LSI，RNA-seq数据选择         PCA，如果是ATAC-seq则使用LSI

cds <- reduce_dimension(cds,preprocess_method = 'PCA')

• 降维方法method = c("PCA", "LSI")，默认是PCA
• 默认维度num_dim是50
• 初步降维前归一化norm_method 默认是log处理
• 设置降维方法是PCA时，自动进行标准化处理

可视化：

plot_cells(cds)

 对分群进行着色

colnames(colData(cds))#可以选任意一个进行着色




#[1] "plate"         "cao_cluster"   "cao_cell_type" "cao_tissue"   



#[5] "Size_Factor" 



plot_cells(cds, color_cells_by="cao_cell_type")

 查看某些基因的表达分布

ciliated_genes <- c("che-1",




                    "hlh-17",




                    "nhr-6",




                    "dmd-6",




                    "ceh-36",




                    "ham-1")




 



plot_cells(cds,




           genes=ciliated_genes,




           label_cell_groups=FALSE,




           show_trajectory_graph=FALSE)

tsne降维

可以看到效果不是很好

#也可以用tsne的方法进行降维，但是效果不好



cds <- reduce_dimension(cds, reduction_method="tSNE")




plot_cells(cds, color_cells_by="cao_cell_type", reduction_method="tSNE")




plot_cells(cds, reduction_method="tSNE", 



           color_cells_by="cao_cell_type",




           cell_size = 2,




           group_label_size = 3)

三、细胞聚类（cluster your cells）

这是细胞分群过程中非常重要的一步。Monocle利用Louvain community detection 这个非监督聚类方法（它是phenoGraph算法的一部分），它和我们常见的Kmeans、hclust层次聚类等还不太一样，这个方法更倾向于找到一个网络中联系紧密的部分，而经常忽略节点的特性；而我们常见的聚类呢，它一般会忽视整体中各个部分的联系，而通过计算两个节点(目标)之间的距离（如欧式距离、曼哈顿距离、余弦相似度等） 找到相似的特性

（一）聚类

这里的聚类已经不再是单单的聚类，按照作者的说法，是把细胞进行分区，后去的轨迹分析都会针对这个分区来进行。partition是一些小的cluster的聚合，如果根据cluster进行着色的话会显得比较乱。

cds <- cluster_cells(cds)#聚类之后，每一个cluster会自成一个”拟时序轨迹“




plot_cells(cds, color_cells_by = "partition")

（二）利用learn graph在每个partition中寻找主路径 

识别轨迹：对每个cluster进行主成分分析，这步以后拟时序图谱就已经初步产生了

cds <- learn_graph(cds)#对每个cluster进行主成分分析，这步以后拟时序图谱就已经初步产生了




plot_cells(cds,




           color_cells_by = "cao_cell_type",




           label_groups_by_cluster=FALSE,




           label_leaves=FALSE,




           label_branch_points=FALSE)




plot_cells(cds,




           color_cells_by = "cao_cell_type",




           label_cell_groups=FALSE,




           label_leaves=TRUE,#展示分支




           label_branch_points=TRUE,#展示分支节点




           graph_label_size=1.5)

 

• 图中的黑线就是整个架构（注意到整个图并非完全连接的，毕竟是按照partition聚类，每个partition中的细胞差异有点大）；
• 浅灰色的圆圈表示”叶片“表示发育轨迹中的不同结局（可以认为拟时序分析的图是一个”根-茎-叶“结构），用参数label_leaves控制；
• 黑色的圆圈表示分支节点，预示着其中的细胞会有不同发展方向，用参数label_branch_points控制；
• 圈中数字大小表示出现时间的先后

四、拟时序分析

（一）首先告诉monocle3，哪个是轨迹分析的起点

1.手动点击添加

cds <- order_cells(cds)




plot_cells(cds,




           color_cells_by = "pseudotime",




           label_cell_groups=FALSE,




           label_leaves=FALSE,




           label_branch_points=FALSE,




           graph_label_size=1.5)

 展示某些基因在时间轨迹中的表达情况

#某个基因的展示



plot_genes_in_pseudotime(cds[1:5,],




                         color_cells_by="cao_cell_type",




                         min_expr=0.5)

 

 2.代码添加，get_earliest_principal_node()函数可以指定起点的名称，在这里指定Am/PH sheath cells这个细胞类型附近作为起点

#方法一



get_earliest_principal_node <- function(cds, my_select="Am/PH sheath cells"){




  cell_ids <- which(colData(cds)[, "cao_cell_type"] == my_select)




  closest_vertex <-




    cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex



  closest_vertex <- (closest_vertex[colnames(cds), ])




  root_pr_nodes <-




    igraph::V(principal_graph(cds)[["UMAP"]])$name[(names




                                                              ((table(closest_vertex[cell_ids,]))))]




  



  root_pr_nodes




}#选择想要定义为起点的细胞类型




cds <- order_cells(cds, root_pr_nodes=get_earliest_principal_node(cds))#定义起点




#另一种方法



myselect <- function(cds,,my_select){




  cell_ids <- which(colData(cds)[,] == my_select)




  closest_vertex <-




    cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex



  closest_vertex <- (closest_vertex[colnames(cds), ])




  root_pr_nodes <-




    igraph::V(principal_graph(cds)[["UMAP"]])$name[(names




                                                              ((table(closest_vertex[cell_ids,]))))]




  



  root_pr_nodes



}



cds <- order_cells(cds, 



                   root_pr_nodes=myselect(cds, = 'cao_cell_type',




                                          my_select = "Body wall muscle")




#绘图



plot_cells(cds,




           color_cells_by = "pseudotime",




           label_cell_groups=FALSE,




           label_leaves=FALSE,




           label_branch_points=FALSE,




           graph_label_size=1.5)|plot_cells(cds,




           color_cells_by = "cao_cell_type",




           label_cell_groups=FALSE,




           label_leaves=FALSE,




           label_branch_points=FALSE,




           graph_label_size=1.5)#两种着色方式

可以看到以pseudotime为依据着色的话，就是以Body wall muscle为起点的紫色区域，序号1就是这个细胞所在区域。

（二）3D版拟时序分析

3D轨迹实际上就是降维时选前3个主成分 => max_components = 3，后续都和2D保持类似

cds_3d <- reduce_dimension(cds, max_components = 3)#降维




cds_3d <- cluster_cells(cds_3d)#聚类




cds_3d <- learn_graph(cds_3d)#轨迹学习




cds <- order_cells(cds, 



                   root_pr_nodes=myselect(cds, = 'cao_cell_type',




                                          my_select = "Body wall muscle")





)#选择拟时序分析的起点，和前面2D一样，用代码选择起点




cds_3d_plot_obj <- plot_cells_3d(cds_3d, color_cells_by="cao_cell_type")#绘图

五、差异基因分析

差异表达分析，进行基因的轨迹分析，找出哪一些基因在这个轨迹中变化最大，graph_test()这个函数可以进行一个叫做Moran's I 的spatial autocorrelation分析方法

subset_pr_test_res <- graph_test(cds, neighbor_graph="principal_graph", cores=4)#cores是指定核心数

运行好之后会得到一个表格，

Moran's I这个指标叫：莫兰指数，澳大利亚统计学家莫兰于1950年提出。它是用来度量空间相关性的一个指标，数据经过方差归一化之后，会落在-1.0和1.0之间。这个值大于0表示空间正相关，值越大空间相关性越强；等于0时空间是随机性分布；小于0时表示空间负相关，值越小空间差异越大
morans_I（空间共表达）,其数值越靠近1代表这个基因在空间距离相近的细胞中表达值越相似，0则代表没有空间共表达效应 

基因特征图

#我们这里就按照morans_I的大小来选择前十个基因



top10 <- subset_pr_test_res %>% top_n(n=10, morans_I) %>%




  pull(gene_short_name) %>% ()




 



#基因表达趋势图，绘图,画不出来，一直报错。



p <- plot_genes_in_pseudotime(cds[top10,],




                              color_cells_by="cao_cell_type", 



                              min_expr=0.5, 



                              ncol = 4)




#基因特征图



plot_cells(cds, genes=top10, show_trajectory_graph=FALSE,




           label_cell_groups=FALSE,  label_leaves=FALSE)

 寻找共表达基因模块

这个是寻找module与细胞/cluster/partition的相关性，module是指协同表达的基因模块，每个模块的基因可能有相似的表达趋势。

##寻找共表达模块



genelist <- pull(subset_pr_test_res, gene_short_name) %>% ()




gene_module <- find_gene_modules(cds[genelist,], resolution=1e-2, cores = 5)




cell_group <- tibble::tibble(cell=(colData(cds)), 



                             cell_group=colData(cds)$celltype)




agg_mat <- aggregate_gene_expression(cds, gene_module, cell_group)




(agg_mat) <- stringr::str_c("Module ", (agg_mat))




pheatmap::pheatmap(agg_mat, scale="column", clustering_method="ward.D2")