❝❞ Monocle2之前写过，已经有一年了快，后来有人看到monocle官网消息，monocle2已经废弃，还埋怨我给的代码过时 ，我也极力解释monocle2分析还是可以做的，所以我出了Monocle3，算是与时俱进！ 这些不管了吧，至少2022最新的高分文章中，Monocle2还是依然被使用的。况且有些人他就喜欢用monocle2，你说气人不气人。这里我更新下之前的monocle2，一方面流程上更加清晰，另一方面，这个可恶的monocle2因为某些包更新后导致出错，这里借鉴了一些修改方法。 首先，我原先的monocle用的2.24.1的版本，但是居然第一步构建monocle对象的时候就出错了。后来网上有人提出使用低版本的monocle，例如2.20.0，但是降低版本会导致其他一些更新的包不能用。所以一不做，二不休，我直接更新了monocle2，直接到2.26.0，然而并没有什么卵用，依然有错。不过从seurat到monocle对象，参考了一个函数，可以一步构建，再也不用提取martix，gene name，防止出错！
第二个错误就是在ordercell那里会出错，BEAM函数那里还会出错。修改方法参考:(https://github.com/cole-trapnell-lab/monocle-release/issues/434）总之就是修改这两个函数的部分内容，然后直接将整个monocle文件夹粘贴到你装R包的路径，直接加载就可以使用了。
setwd("D:/KS项目/公众号文章/monocle2更新版")library(Seurat)devtools::load_all("C:/Users/AppData/Local/R/win-library/4.2/monocle")library(dplyr)
mouse_monocle_fun <- seurat_to_monocle(mouse_data, assay = "RNA", slot = "counts")
#===========================================================================#Estimate size factors and dispersionsmouse_monocle <- estimateSizeFactors(mouse_monocle_fun)mouse_monocle <- estimateDispersions(mouse_monocle)#质量控制mouse_monocle <- detectGenes(mouse_monocle, min_expr = 0.1)print(head(fData(mouse_monocle)))print(head(pData(mouse_monocle)))expressed_genes <- row.names(subset(fData(mouse_monocle), num_cells_expressed >= 10))pData(mouse_monocle)$Total_mRNAs <- Matrix::colSums(exprs(mouse_monocle))mouse_monocle <- mouse_monocle[,pData(mouse_monocle)$Total_mRNAs < 1e6]
#====================method1====================================================#monocle计算高变基因cds_monocle <- mouse_monocledisp_table <- dispersionTable(cds_monocle)unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)head(unsup_clustering_genes)cds_monocle <- setOrderingFilter(cds_monocle, unsup_clustering_genes$gene_id)save(cds_monocle, file = 'cds_monocle.RData')#====================method2====================================================#使用seurat计算cds_seurat <- mouse_monocleseurat_variable_genes <- VariableFeatures(FindVariableFeatures(mouse_data, assay = "RNA"), assay = "RNA")cds_seurat <- setOrderingFilter(cds_seurat, seurat_variable_genes)plot_ordering_genes(cds_DGT)save(cds_seurat, file = 'cds_seurat.RData')
#====================method3====================================================#differentialGeneTestcds_DGT <- mouse_monoclediff_test_res <- differentialGeneTest(cds_DGT,fullModelFormulaStr = "~celltype")ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))# ordering_genes <- row.names(diff_test_res)[order(diff_test_res$qval)][1:2000]#或者选择前多少的基因cds_DGT <- setOrderingFilter(cds_DGT, ordering_genes)plot_ordering_genes(cds_DGT)save(cds_DGT, file = "cds_DGT.RData")
##################################################################################接下来我们以第三种方法进行#降维,排序cds_DGT <- reduceDimension(cds_DGT, max_components = 2,reduction_method = 'DDRTree')cds_DGT <- orderCells(cds_DGT)cds_DGT <- orderCells(cds_DGT, root_state = NULL, num_paths = NULL, reverse = T) ##################################################################################初步可视化plot_cell_trajectory(cds_DGT, cell_size = 2.2, color_by = "celltype") + facet_wrap(~celltype, nrow = 2)plot_cell_trajectory(cds_DGT, color_by = "Pseudotime")plot_cell_trajectory(cds_DGT, color_by = "orig.ident",cell_link_size = 1.5)plot_cell_trajectory(cds_DGT, color_by = "State")
#基因表达拟时图plot_cell_trajectory(cds_DGT, markers = "Il1b",use_color_gradient=T,cell_size = 1,cell_link_size = 1.5)#各个分组拟时图plot_cell_trajectory(cds_DGT, color_by = "orig.ident") + facet_wrap(~orig.ident, nrow = 2)
##################################################################################查看随着pseudotime变化的基因cds_DGT_pseudotimegenes <- differentialGeneTest(cds_DGT,fullModelFormulaStr = "~sm.ns(Pseudotime)")cds_DGT_pseudotimegenes_sig <- subset(cds_DGT_pseudotimegenes, qval < 0.01)write.csv(cds_DGT_pseudotimegenes_sig,"cds_DGT_pseudotimegenes_sig.csv")
#筛选一些基因去做图，这里完全是为了作图，实际展示基因按照自己按要求cds_DGT_pseudotimegenes_sig <- subset(cds_DGT_pseudotimegenes, qval < 1e-320)Time_genes <- cds_DGT_pseudotimegenes_sig %>% pull(gene_short_name) %>% as.character()p <- plot_pseudotime_heatmap(cds_DGT[Time_genes,], num_cluster = 4, show_rownames = T, return_heatmap = T)
就本例中，我们的细胞没有多余的State，所以也用不到branch分析了。但是当我用BEAM函数的时候，居然还会出错，具体的错误如下：
BEAM_res <- BEAM(cds_DGT, branch_point = 1, cores = 1)#Error in if (isSparseMatrix(exprs(X))) { : the condition has length > 1#In addition: Warning message:#In class(cellData) != "matrix" && isSparseMatrix(cellData) == FALSE :#解决办法就是将cores提高，只要不是1就行
BEAM_res <- BEAM(cds_DGT, branch_point = 1, cores = 2)#Error in checkForRemoteErrors(val) : # 2 nodes produced errors; first error: could not find function "calculate_NB_dispersion_hint"#In addition: Warning message:
结果依然出错，找不到这个calculate_NB_dispersion_hint函数，将其添加到BEAM.R中依然没有起色，果断放弃，不做这一步也可以。不晓得到底哪里的问题，不过有些小伙伴不会遇到这样的问题，听天由命吧，毕竟没有monocle2，我们还有monocle3！
|
