之前的课程我们已经对monocle2做了一个完备的教程，很多粉丝也一直在催monocle3的教程，这里我们带领大家先简单的认识一下monocle3的特点与功能。原本是想简单介绍一下，由于monocle3实在是漏洞百出，没想到边吐槽边漫谈居然录了五十分钟。视频中提到的Rmakdown还没有制作完，大家不要着急，我们先把这节课的测试文件与代码放在这个链接中，monocle3系列课程制作完毕会再进行更新，请大家持续关注：

你也可以直接前往monocle3官网教程进行学习：

https://cole-trapnell-lab.github.io/monocle3/docs/differential/

视频总是上传失败，大家直接去B站看吧：
https://www.bilibili.com/video/BV1br4y1x7Hf?p=8

一、初识monocle

1.monocle3与monocle2的主要区别

1.1.monocle3的主要功能

Clustering, classifying, and counting cells：主要任务是帮助大家找到一些具有特殊功能的亚群。
Constructing single-cell trajectories：这个功能大家需要monocle的最根本原因，通过拟时序分析帮助大家解析生物体发育、疾病等过程中细胞发生的变化。
Differential expression analysis：差异计算monocle2也有，但我们实测的感觉是不如Seurat的，这里作者表明monocle3在这一块进行了优化升级，我们后面具体看看用起来效果如何。

1.2.先谈缺点

目前monocle3已经更新到β版本了，作者在官网也承认了缺点，monocle3 α已经是不推荐使用的，可能会存在bug，但是monocle3 β仍然处于搭建中的状态，也就是说monocle3仍然是可能存在bug的，并且我们之前讲绪论的时候说到monocle1、2都发表在了Nature系列之上，但是monocle3迟迟没有发表，并且目前发表的文章还是使用monocle2的比较多，monocle3的不稳定性可能是重要原因。不过我相信新发表的论文很快就会由monocle2转向monocle3。
官方表述如下

image.png

1.3.再看看优点
咱们抛开可能存在的bug不谈，monocle3相较于monocle2具有以下几点优势：(1)最大的优点就是计算量变大了，可以处理百万级别的单细胞数据集，也就是说整个器官、甚至整个胚胎的矩阵交给monocle3处理完全没压力。
(2)代码结构性优化：这点我要吐槽一下，monocle系列的语法我一直觉得很奇怪，默认参数也很不人性化
(3)支持UMAP算法的降维，这个也非常Nice，速度比tSNE快的不是一星半点。
(4)支持多谱系的拓扑结构：换句话说拟时序的轨迹可以做的很复杂
(5)相较于原来的RGE算法，新的approximate graph abstraction能够计算不连续的、平行的拓扑结构
(6)新的基因表达量计算及差异分析方法被引入，也就是说原来的differentialGeneTest()和BEAM()可以被替代。
(7)可以像Seurat的多样本整合一样对拟时序对象进行整合：这个功能可以说是刚需了，换句话说，如果你有合适的、已构建好拟时序的参考数据集，可以直接把自己的数据跟参考数据集进行投影、比对。
(8)数据整合时也可同时加上注释：这有点类似于Seurat中的TransferData，可以利用已经注释好的参考数据给现行数据添加注释。
(9)对monocle对象的读取、加载、转换做了一定的优化，我们后面可以看看效果如何
(10)优化了负二项分布模型：也就是说对处理count的优化
(11)可视化提供了3D展示功能

2.安装Monocle 3
保证自己的环境符合这些版本要求，不然报错很麻烦 R >= 4.1.0 , Bioconductor = 3.14,这样装上的monocle3 >= 1.2.7

2.2.R语言中的安装

setwd('~/analysis/monocle3/')
if (!requireNamespace("BiocManager", quietly = TRUE))  install.packages("BiocManager")package.version
('BiocManager')BiocManager::install(version = "3.14")#先把依赖包都装上
if(!require(monocle3))BiocManager::install(c('BiocGenerics', 'DelayedArray', 'DelayedMatrixStats',                                             'limma', 'lme4', 'S4Vectors', 'SingleCellExperiment',                                             'SummarizedExperiment', 'batchelor', 'Matrix.utils',                                             'HDF5Array', 'terra', 'ggrastr'))if(!require(devtools))install.packages("devtools")
if(!require(monocle3))devtools::install_github('cole-trapnell-lab/monocle3')#这是稳定版
if(!require(monocle3))devtools::install_github('cole-trapnell-lab/monocle3', ref="develop")#这是测试版

2.3.通过conda安装
即使在R4.1.0中还是很容易有很多包装不上

#如果上面这些包安装时出现了这种报错：
#configure: error: gdal-config not found or not executable.
#那么你需要在Linux中运行：
#sudo apt-get -y update && apt-get install -y  libudunits2-dev libgdal-dev libgeos-dev libproj-dev
#需要有root权限
#当然我还是推荐用conda安装，作者要求装什么版本咱就装什么版本，以下注释代码在Linux中完成
#conda activate monocle3
#conda insyall r-base==4.1.0
#conda install gcc#conda install seurat
#conda install r-seuratobject
#conda install r-terra==1.5.12
#conda install r-units
#conda install r-raster
# conda install r-spdata
# conda install r-sf
# conda install r-spdep
#  conda install r-ragg
#  conda install r-ggrastr
#  conda install r-devtoolsif(!require(monocle3))devtools::install_github('cole-trapnell-lab/monocle3')
#安装好上述依赖包后可运行
#conda search monocle3#conda install r-monocle3
#依赖的gdal以及python也会自动装好
#这个时候装好monocle3的R的路径应该在~/miniconda3/envs/monocle3/bin/R
#安装好了都加载一下试试library(monocle3)package.version("monocle3")
## [1] "1.0.0"library(Seurat)library(SeuratObject)library(tidyselect)library(dplyr)

2.4.版本信息
总之，无论你通过哪种方式安装，要保证我们以下软件的版本信息均一致，避免不必要的麻烦

sessionInfo()

sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /home/biomamba/miniconda3/envs/newmonocle3/lib/libopenblasp-r0.3.20.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] tidyselect_1.1.2            sp_1.5-0                   
##  [3] SeuratObject_4.1.0          Seurat_4.1.1               
##  [5] monocle3_1.0.0              SingleCellExperiment_1.16.0
##  [7] SummarizedExperiment_1.24.0 GenomicRanges_1.46.1       
##  [9] GenomeInfoDb_1.30.1         IRanges_2.28.0             
## [11] S4Vectors_0.32.3            MatrixGenerics_1.6.0       
## [13] matrixStats_0.62.0          Biobase_2.54.0             
## [15] BiocGenerics_0.40.0         EBImage_4.36.0             
## [17] openintro_2.3.0             usdata_0.2.0               
## [19] cherryblossom_0.1.0         airports_0.1.0             
## [21] forcats_0.5.1               stringr_1.4.0              
## [23] dplyr_1.0.9                 purrr_0.3.4                
## [25] readr_2.1.2                 tidyr_1.2.0                
## [27] tibble_3.1.7                ggplot2_3.3.6              
## [29] tidyverse_1.3.1            
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.4.0           backports_1.4.1        plyr_1.8.7            
##   [4] igraph_1.3.0           lazyeval_0.2.2         splines_4.1.0         
##   [7] listenv_0.8.0          scattermore_0.8        digest_0.6.29         
##  [10] htmltools_0.5.2        viridis_0.6.2          tiff_0.1-11           
##  [13] fansi_1.0.3            magrittr_2.0.3         tensor_1.5            
##  [16] cluster_2.1.3          ROCR_1.0-11            tzdb_0.3.0            
##  [19] globals_0.15.1         modelr_0.1.8           spatstat.sparse_2.1-1 
##  [22] jpeg_0.1-9             colorspace_2.0-3       rvest_1.0.2           
##  [25] ggrepel_0.9.1          haven_2.5.0            xfun_0.31             
##  [28] crayon_1.5.1           RCurl_1.98-1.7         jsonlite_1.8.0        
##  [31] spatstat.data_2.2-0    progressr_0.10.1       survival_3.3-1        
##  [34] zoo_1.8-10             glue_1.6.2             polyclip_1.10-0       
##  [37] gtable_0.3.0           zlibbioc_1.40.0        XVector_0.34.0        
##  [40] leiden_0.4.2           DelayedArray_0.20.0    future.apply_1.9.0    
##  [43] abind_1.4-5            scales_1.2.0           DBI_1.1.3             
##  [46] spatstat.random_2.2-0  miniUI_0.1.1.1         Rcpp_1.0.8.3          
##  [49] xtable_1.8-4           viridisLite_0.4.0      spatstat.core_2.4-4   
##  [52] reticulate_1.25        htmlwidgets_1.5.4      httr_1.4.3            
##  [55] RColorBrewer_1.1-3     ellipsis_0.3.2         ica_1.0-2             
##  [58] pkgconfig_2.0.3        uwot_0.1.11            deldir_1.0-6          
##  [61] sass_0.4.1             dbplyr_2.2.0           locfit_1.5-9.5        
##  [64] utf8_1.2.2             rlang_1.0.2            reshape2_1.4.4        
##  [67] later_1.2.0            munsell_0.5.0          cellranger_1.1.0      
##  [70] tools_4.1.0            cli_3.3.0              generics_0.1.2        
##  [73] broom_0.8.0            ggridges_0.5.3         evaluate_0.15         
##  [76] fastmap_1.1.0          fftwtools_0.9-11       goftest_1.2-3         
##  [79] yaml_2.3.5             knitr_1.39             fs_1.5.2              
##  [82] fitdistrplus_1.1-8     RANN_2.6.1             nlme_3.1-158          
##  [85] pbapply_1.5-0          future_1.26.1          mime_0.12             
##  [88] xml2_1.3.3             compiler_4.1.0         rstudioapi_0.13       
##  [91] plotly_4.10.0          png_0.1-7              spatstat.utils_2.3-1  
##  [94] reprex_2.0.1           bslib_0.3.1            stringi_1.7.6         
##  [97] highr_0.9              rgeos_0.5-9            lattice_0.20-45       
## [100] Matrix_1.4-1           vctrs_0.4.1            pillar_1.7.0          
## [103] lifecycle_1.0.1        spatstat.geom_2.4-0    lmtest_0.9-40         
## [106] jquerylib_0.1.4        RcppAnnoy_0.0.19       data.table_1.14.2     
## [109] cowplot_1.1.1          bitops_1.0-7           irlba_2.3.5           
## [112] patchwork_1.1.1        httpuv_1.6.5           R6_2.5.1              
## [115] promises_1.2.0.1       KernSmooth_2.23-20     gridExtra_2.3         
## [118] parallelly_1.32.0      codetools_0.2-18       MASS_7.3-57           
## [121] assertthat_0.2.1       withr_2.5.0            sctransform_0.3.3     
## [124] GenomeInfoDbData_1.2.7 mgcv_1.8-40            parallel_4.1.0        
## [127] hms_1.1.1              rpart_4.1.16           grid_4.1.0            
## [130] rmarkdown_2.14         Rtsne_0.16             shiny_1.7.1           
## [133] lubridate_1.8.0

3.开始学习
3.1.monocle3对象的构建

#还记得我们在monocle2中教过大家的吗，monocle对象的构建
# data <- as(as.matrix(data4mono@assays$RNA@counts), 'sparseMatrix')
# pd <- new('AnnotatedDataFrame', data = data4mono@meta.data)
# fData <- data.frame(gene_short_name = row.names(data4mono), row.names = row.names(data4mono))
# fd <- new('AnnotatedDataFrame', data = fData)
# mycds <- newCellDataSet(data,#                         phenoData = pd,
#                         featureData = fd,
#                         expressionFamily = negbinomial.size())
#monocle2的作者比momocle2良心多了，不仅提供了测试数据的链接
# 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"))
#我已经替容易被墙的同学下到本地了：
expression_matrix <- readRDS('author.pro/expression_matrix.rds')
cell_metadata <- readRDS('author.pro/cell_metadata.rds')
gene_annotation <- readRDS('author.pro/gene_annotation.rds')
#还把蹩脚的AnnotatedDataFrame输入格式改成了dataframeclass(cell_metadata)

## [1] "data.frame"

#缺点就是这个cds文件好像有点大。。。可能是为了展示monocle3对于大型数据的计算能力
#构建cds对象，cds <- new_cell_data_set(expression_matrix,                         cell_metadata = cell_metadata,                         gene_metadata = gene_annotation)
#monocle2的函数名是，newCellDataSet，一定看清楚函数名，不要 function not found 了还不知道哪里出错了
#取个子集加快演示速度mysubset <- c()
for(runplate in unique(cds@colData$plate)){  mytemp<- cds@colData %>% as.data.frame() %>% filter(plate==runplate) %>% 
rownames() %>% sample(50)  
mysubset <- cbind(mysubset,mytemp)}#mysubset <- sample(1:20271,2000)
cds <- cds[,mysubset]table(cds@colData$plate)

## 
## 001 002 003 004 005 006 007 008 009 010 
##  50  50  50  50  50  50  50  50  50  50

saveRDS(cds,'test.data/simple.rds')#另外cds对象现在也可以view了
try(cds <- load_cellranger_data("test.data/filtered_gene_bc_matrices/hg19/"))

## Error in load_cellranger_data("test.data/filtered_gene_bc_matrices/hg19/") : 
##   Could not find directory: test.data/filtered_gene_bc_matrices/hg19//outs/filtered_gene_bc_matrices

#依然是报错#10X文件可以直接这么读取
from10X <- load_mm_data(mat_path = "test.data/filtered_gene_bc_matrices/hg19/outs/matrix.mtx",                         feature_anno_path = "test.data/filtered_gene_bc_matrices/hg19/outs/genes.tsv",                         cell_anno_path = "test.data/filtered_gene_bc_matrices/hg19/outs/barcodes.tsv")
#这样貌似是没问题的#如果两种方法你都运行不了，咱们还是曲线救国，自己想办法吧
pbmc <- Read10X('test.data/filtered_gene_bc_matrices/hg19/outs/')
pbmc <- CreateSeuratObject(counts = pbmc)
seurat2cds <- new_cell_data_set(as(as.matrix(pbmc@assays$RNA@counts), 'sparseMatrix'),                                cell_metadata = as.data.frame(pbmc@meta.data),                                gene_metadata = data.frame(gene_short_name = row.names(pbmc),                                                            row.names = row.names(pbmc)))
#成功，看一下数据类型class(seurat2cds)#没问题，是cds

## [1] "cell_data_set"
## attr(,"package")
## [1] "monocle3"

#总体使用起来还是感觉版本之间有些不兼容

3.2.初步学习一下

cds <- preprocess_cds(cds, num_dim = 100)#抽平、预处理，这一步默认调用最大线程cds <- align_cds(cds, alignment_group = "plate")#去除批次效应，看起来要比Seurat中方便的多

## Aligning cells from different batches using Batchelor. 
## Please remember to cite:
##   Haghverdi L, Lun ATL, Morgan MD, Marioni JC (2018). 'Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors.' Nat. Biotechnol., 36(5), 421-427. doi: 10.1038/nbt.4091

names(cds@colData)#引号里的内容在这找

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

## 降维，默认是"Aligned"方式cds <- reduce_dimension(cds,cores=5)

## No preprocess_method specified, and aligned coordinates have been computed previously. Using preprocess_method = 'Aligned'

## Note: reduce_dimension will produce slightly different output each time you run it unless you set 'umap.fast_sgd = FALSE' and 'cores = 1'

#可以支持多线程、点赞，如果你下游参数都选择默认的话，这里要选择UMAP
#每次降维可能效果均有细微差别，但是'umap.fast_sgd = FALSE' and 'cores = 1'可以保证结果相同
#试试设置随机数是否可以解决这个问题
#cds <- reduce_dimension(cds,reduction_method = 'tSNE',cores=5)#如果你还是想用tSNE
## 聚类分群
cds <- cluster_cells(cds)#不能手动调用多线程，但是其实这步运算起来很快## 拟时序
cds <- learn_graph(cds)#这一步就有点慢了，即使我们的数据很小，这步主要是通过降维数据执行reversed graph embedding(RGE)算法

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%

这一步整活了，可以交互式地选择起点，看了一下这个功能应该是依赖shiny写的，那就意味着在终端操作这一步就不可能了

cds <- order_cells(cds)#速度倒是挺快，即使是7GB的对象，几分钟也就跑完了

选择过程如下：

image.png

plot_cells(cds)#看看效果如何

image.png

names(cds@colData)
## [1] "plate"         "cao_cluster"   "cao_cell_type" "cao_tissue"   ## [5] "Size_Factor"#寻找差异基因，这里相当于Seurat::FindAllMarkers()
gene_fits <- fit_models(cds, model_formula_str = "~plate")
fit_coefs <- coefficient_table(gene_fits)
emb_time_terms <- fit_coefs %>% filter(term == "plate")
emb_time_terms <- emb_time_terms %>% mutate(q_value = p.adjust(p_value))sig_genes <- emb_time_terms %>% filter (q_value < 0.05) %>% pull(gene_short_name)

pr_test_res <- graph_test(cds,                            neighbor_graph="principal_graph",                           cores=1,verbose = F)#这个函数很讨厌，在Rmarkdown中会输出很长的进度条

pr_deg_ids <- row.names(subset(pr_test_res, q_value < 0.05))
head(pr_deg_ids)#看一下最终的基因列表

## [1] "WBGene00000001" "WBGene00000006" "WBGene00000010" "WBGene00000020"
## [5] "WBGene00000024" "WBGene00000030"

最后，取子集与合并

cds1 <- cds[1:100,]

cds2 <- cds[101:200,]big_cds <- combine_cds(list(cds, cds2))

## Warning in estimate_size_factors(cds): Your CDS object contains cells with zero
## reads. This causes size factor calculation to fail. Please remove the zero read
## cells using cds <- cds[,Matrix::colSums(exprs(cds)) != 0] and then run cds <-
## estimate_size_factors(cds)

欢迎关注同名公众号~