其他
单细胞测序数据进阶分析—《拟时序分析》4.初识monocle3
之前的课程我们已经对monocle2做了一个完备的教程，很多粉丝也一直在催monocle3的教程，这里我们带领大家先简单的认识一下monocle3的特点与功能。原本是想简单介绍一下，由于monocle3实在是漏洞百出，没想到边吐槽边漫谈居然录了五十分钟。视频中提到的Rmakdown还没有制作完，大家不要着急，我们先把这节课的测试文件与代码放在这个链接中，monocle3系列课程制作完毕会再进行更新，请大家持续关注：
代码及测试文件链接：https://pan.baidu.com/s/16lPDInnenycR9sTcABolSg?pwd=tdp7你也可以直接前往monocle3官网教程进行学习：
https://cole-trapnell-lab.github.io/monocle3/docs/differential/
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。
官方表述如下
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语言中的安装
Hide
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中还是很容易有很多包装不上
Hide
#如果上面这些包安装时出现了这种报错：#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.版本信息
总之，无论你通过哪种方式安装，要保证我们以下软件的版本信息均一致，避免不必要的麻烦
Hide
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.03.开始学习
3.1.monocle3对象的构建
Hide
#还记得我们在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"Hide
#缺点就是这个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 50Hide
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_matricesHide
#依然是报错#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"Hide
#总体使用起来还是感觉版本之间有些不兼容3.2.初步学习一下
Hide
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.4091Hide
names(cds@colData)#引号里的内容在这找## [1] "plate" "cao_cluster" "cao_cell_type" "cao_tissue"
## [5] "Size_Factor"Hide
## 降维，默认是"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'Hide
#可以支持多线程、点赞，如果你下游参数都选择默认的话，这里要选择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写的，那就意味着在终端操作这一步就不可能了
Hide
cds <- order_cells(cds)#速度倒是挺快，即使是7GB的对象，几分钟也就跑完了选择过程如下：
Hide
plot_cells(cds)#看看效果如何Hide
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)Hide
pr_test_res <- graph_test(cds, neighbor_graph="principal_graph", cores=1,verbose = F)#这个函数很讨厌，在Rmarkdown中会输出很长的进度条Hide
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"最后，取子集与合并
Hide
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)拟时序系列教程
B站连续播放起来比较方便：
https://www.bilibili.com/video/BV1br4y1x7Hf?p=1往期推送单细胞测序数据进阶分析—《拟时序分析》1.概论
单细胞测序数据进阶分析—《拟时序分析》2.monocle概论单细胞测序数据进阶分析—《拟时序分析》3.monocle2实操：完整版
往期回顾
单细胞数据分析系列教程：
B站视频，先看一遍视频再去看推送操作，建议至少看三遍：https://www.bilibili.com/video/BV1S44y1b76Z/单细胞测序基础数据分析保姆级教程，代码部分整理在往期推送之中：手把手教你做单细胞测序数据分析（一）——绪论
手把手教你做单细胞测序数据分析（二）——各类输入文件读取
手把手教你做单细胞测序数据分析（三）——单样本分析
手把手教你做单细胞测序数据分析（四）——多样本整合手把手教你做单细胞测序数据分析（五）——细胞类型注释
手把手教你做单细胞测序数据分析（六）——组间差异分析及可视化手把手教你做单细胞测序数据分析（七）——基因集富集分析
上游fastq文件处理：
单细胞分析的最上游——处理Fastq文件：cellranger
单细胞分析的最上游——处理Fastq文件：dropseqRunner有奖提问,这个smart-Seq2数据实战的比对率为何如此低？
细胞通讯B站连续播放起来比较方便：
https://www.bilibili.com/video/BV1Ab4y1W7qx?p=1往期推送单细胞测序数据进阶分析—《细胞通讯》1.概论
单细胞测序数据进阶分析—《细胞通讯》2.1CellChat基础分析教程
单细胞测序数据进阶分析—《细胞通讯》2.2CellChat多组别分析
给你们申请到了免费的128GB云服务器
拟时序分析B站连续播放起来比较方便：
https://www.bilibili.com/video/BV1br4y1x7Hf?p=1往期推送单细胞测序数据进阶分析—《拟时序分析》1.概论
单细胞测序数据进阶分析—《拟时序分析》2.monocle概论单细胞测序数据进阶分析—《拟时序分析》3.monocle2实操：完整版
其他单细胞相关技术贴也在这里：细胞的数量由誰决定？单细胞中应该如何做GSVA？答读者问（三）：单细胞测序前景答读者问（四）：如何分析细胞亚群答读者问（六）：Seurat中如何让细胞听你指挥答读者问（八）：为什么Read10X也会报错?
答读者问（十）整合后的表达矩阵，如何拆分出分组信息？
答读者问（十一）如何一次性读取一个目录下的cellranger输出文件？
给你安排一个懂生信的工具人（十）：不学编程 零代码完成单细胞测序数据分析：Loupe Browser
什么？不做单细胞也能分析细胞类群和免疫浸润？
答读者问 （十三）查看Seurat对象时的ERROR：type='text'
各类单细胞对象（数据格式）转换大全（一）批量整理好GEO中下载的单细胞数据
答读者问 （十四）Seurat中分类变量处理技巧
答读者问 （十五）稀疏矩阵转matrix, as.matrix函数是下下策
答读者问 （十六）做单细胞测序到底需要多少内存
单细胞文献阅读：文献阅读（三）、单细胞测序解析糖尿病肾病中肾小球的动态变化
文献阅读（四）、单细胞测序技术解析健康人与T2D患者的胰岛差异文献阅读（六）、小鼠全肾单细胞测序开篇之作文献阅读（七）、一篇不花钱就能白嫖的文章文献阅读（八）、不会吧不会吧，Nature都能白嫖？文献阅读（十一）、高氧下小鼠肺发育损伤的ScRNA图谱文献阅读（十二）、IgAN & STRT-Seq
文献阅读（十三）、老树开新花——EGFR、肿瘤、免疫+scRNA-Seq
文献阅读（十五）、癌前基质细胞驱动BRCA1肿瘤发生
非技术帖：关于单细胞的事 谈谈后面的计划趁机预告一波有关提高"Biomamba 生信基地"运行效率事宜
独享服务器2.0即将上线，128GB免费试用
如何联系我们
大家可以阅读完这几篇之后添加笑一笑也就算了有关提高"Biomamba 生信基地"运行效率事宜