Monocle3 r语言 - 锐客网

一、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("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"))# Make the CDS object cds <- new_cell_data_set(expression_matrix, cell_metadata = https://www.it610.com/article/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 features.tsv.gz、barcodes.tsv.gz 和 matrix.mtx.gz。（还可以处理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
matrix.mtx,feature.tsv和barcodes.tsv

cds <- load_mm_data(mat_path = "~/Downloads/matrix.mtx", feature_anno_path = "~/Downloads/features.tsv", cell_anno_path = "~/Downloads/barcodes.tsv")

（三）处理大量数据当细胞量或基因量比较大的话，可以使用稀疏矩阵的方法读取

cds <- new_cell_data_set(as(umi_matrix, "sparseMatrix"), cell_metadata = https://www.it610.com/article/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")

文章图片

【Monocle3】
查看某些基因的表达分布

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 <- as.matrix(closest_vertex[colnames(cds), ]) root_pr_nodes <- igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names (which.max(table(closest_vertex[cell_ids,]))))]root_pr_nodes }#选择想要定义为起点的细胞类型 cds <- order_cells(cds, root_pr_nodes=get_earliest_principal_node(cds))#定义起点 #另一种方法 myselect <- function(cds,select.classify,my_select){ cell_ids <- which(colData(cds)[,select.classify] == my_select) closest_vertex <- cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex closest_vertex <- as.matrix(closest_vertex[colnames(cds), ]) root_pr_nodes <- igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names (which.max(table(closest_vertex[cell_ids,]))))]root_pr_nodes } cds <- order_cells(cds, root_pr_nodes=myselect(cds,select.classify = '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,select.classify = '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) %>% as.character()#基因表达趋势图，绘图,画不出来，一直报错。 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) %>% as.character() gene_module <- find_gene_modules(cds[genelist,], resolution=1e-2, cores = 5) cell_group <- tibble::tibble(cell=row.names(colData(cds)), cell_group=colData(cds)$celltype) agg_mat <- aggregate_gene_expression(cds, gene_module, cell_group) row.names(agg_mat) <- stringr::str_c("Module ", row.names(agg_mat)) pheatmap::pheatmap(agg_mat, scale="column", clustering_method="ward.D2")