摘要：本笔记展示了 Monocle3 拟时序分析的详细过程。

# 简介

在发展过程中、对刺激做出反应以及在整个生命周期中，细胞从一种功能 “状态” 过渡到另一种状态。处于不同状态的细胞表达不同的基因组，产生动态的蛋白质和代谢物组合来执行它们的工作。当细胞在不同状态间移动时，它们经历一个转录重组的过程，其中一些基因被沉默，而其他基因被新激活。这些瞬态状态往往难以表征，因为在更稳定的终点状态之间纯化细胞可能是困难或不可能的。单细胞 RNA 测序可以让你在无需纯化的情况下观察这些状态。然而，要做到这一点，我们必须确定每个细胞在可能状态范围内的位置。

Monocle 首次引入了使用 RNA 测序进行单细胞轨迹分析的策略。 Monocle 不是通过实验将细胞纯化成离散状态，而是使用算法来学习细胞作为动态生物过程的一部分必须经历的基因表达变化序列。一旦学习了基因表达变化的总体 “轨迹”， Monocle 就可以将每个细胞放置在轨迹中的适当位置。然后，你可以使用 Monocle 的差异分析工具包来找到在轨迹过程中受调控的基因，如 “伪时间功能变化中的基因寻找” 部分所述。如果这个过程有多个结果， Monocle 将重建一个 “分支” 轨迹。这些分支对应于细胞的 “决策”， Monocle 提供了强大的工具来识别受它们影响并参与决策制定的基因。你可以在 “单细胞轨迹中分析分支” 部分看到如何分析分支。

# 分析流程
## 加载包

``` r
library(monocle3)

载入需要的程辑包：Biobase

载入需要的程辑包：BiocGenerics


载入程辑包：'BiocGenerics'

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    anyDuplicated, append, as.data.frame, basename, cbind, colnames,
    dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
    grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
    rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which.max, which.min

Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

载入需要的程辑包：SingleCellExperiment

载入需要的程辑包：SummarizedExperiment

载入需要的程辑包：MatrixGenerics

载入需要的程辑包：matrixStats


载入程辑包：'matrixStats'

The following objects are masked from 'package:Biobase':

    anyMissing, rowMedians


载入程辑包：'MatrixGenerics'

The following objects are masked from 'package:matrixStats':

    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
    colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
    colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
    colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
    colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
    colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
    colWeightedMeans, colWeightedMedians, colWeightedSds,
    colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
    rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
    rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
    rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
    rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
    rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
    rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
    rowWeightedSds, rowWeightedVars

The following object is masked from 'package:Biobase':

    rowMedians

载入需要的程辑包：GenomicRanges

载入需要的程辑包：stats4

载入需要的程辑包：S4Vectors


载入程辑包：'S4Vectors'

The following objects are masked from 'package:base':

    expand.grid, I, unname

载入需要的程辑包：IRanges

载入需要的程辑包：GenomeInfoDb


载入程辑包：'monocle3'

The following objects are masked from 'package:Biobase':

    exprs, fData, fData<-, pData, pData<-

library(dplyr)

载入程辑包：'dplyr'

The following objects are masked from 'package:GenomicRanges':

    intersect, setdiff, union

The following object is masked from 'package:GenomeInfoDb':

    intersect

The following objects are masked from 'package:IRanges':

    collapse, desc, intersect, setdiff, slice, union

The following objects are masked from 'package:S4Vectors':

    first, intersect, rename, setdiff, setequal, union

The following object is masked from 'package:matrixStats':

    count

The following object is masked from 'package:Biobase':

    combine

The following objects are masked from 'package:BiocGenerics':

    combine, intersect, setdiff, union

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(Seurat)

载入需要的程辑包：SeuratObject

载入需要的程辑包：sp

The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
which was just loaded, were retired in October 2023.
Please refer to R-spatial evolution reports for details, especially
https://r-spatial.org/r/2023/05/15/evolution4.html.
It may be desirable to make the sf package available;
package maintainers should consider adding sf to Suggests:.


载入程辑包：'sp'

The following object is masked from 'package:IRanges':

    %over%


载入程辑包：'SeuratObject'

The following object is masked from 'package:SummarizedExperiment':

    Assays

The following object is masked from 'package:GenomicRanges':

    intersect

The following object is masked from 'package:GenomeInfoDb':

    intersect

The following object is masked from 'package:IRanges':

    intersect

The following object is masked from 'package:S4Vectors':

    intersect

The following object is masked from 'package:BiocGenerics':

    intersect

The following object is masked from 'package:base':

    intersect


载入程辑包：'Seurat'

The following object is masked from 'package:SummarizedExperiment':

    Assays

# 数据预处理

重建轨迹的工作流程与聚类的工作流程非常相似，但它有一些额外的步骤。为了说明这个工作流程，我们将使用另一个秀丽隐杆线虫数据集，这个数据集来自 Packer & Zhu 等人的研究。他们的研究包括对整个发育中胚胎的时间序列分析。我们将检查数据的一个小子集，其中包括大部分神经元。我们将像处理 L2 数据一样加载它：

expression_matrix <- readRDS(url("https://depts.washington.edu:/trapnell-lab/software/monocle3/celegans/data/packer_embryo_expression.rds"))

cell_metadata <- readRDS(url("https://depts.washington.edu:/trapnell-lab/software/monocle3/celegans/data/packer_embryo_colData.rds"))

gene_annotation <- readRDS(url("https://depts.washington.edu:/trapnell-lab/software/monocle3/celegans/data/packer_embryo_rowData.rds"))

cds <- new_cell_data_set(expression_matrix,

cell_metadata = cell_metadata,

gene_metadata = gene_annotation)

数据预处理的工作方式与聚类分析中的完全相同。这次，我们将使用一种不同的批次校正策略，包括 Packer & Zhu 等人在他们原始分析中所做的：

注意：您的数据不会有这里演示的加载批次信息，您将使用您自己的批次信息进行批次校正。

cds <- preprocess_cds(cds, num_dim = 50)

cds <- align_cds(cds, alignment_group = "batch", residual_model_formula_str = "~ bg.300.loading + bg.400.loading + bg.500.1.loading + bg.500.2.loading + bg.r17.loading + bg.b01.loading + bg.b02.loading")

除了使用 align_cds() 的 alignment_group 参数对细胞组（即批次）进行对齐外，我们还使用 residual_model_formula_str。这个参数用于减去连续效应。您可以使用它来控制诸如每个细胞中线粒体读数的比例之类的事情，这有时被用作每个细胞的质量控制指标。在这个实验中（就像在许多单细胞 RNA 测序实验中一样），一些细胞会自发裂解，释放它们的 mRNA 到单细胞文库制备前的细胞悬液中。这种 “上清 RNA” 在一定程度上污染了每个细胞的转录组轮廓。幸运的是，估算每批细胞中背景污染的水平并减去它是相当直接的，这正是 Packer 等人在原始研究中所做的。每一列 bg.300.loading , bg.400.loading 等对应于一个细胞可能被污染的背景信号。将这些列作为 residual_model_formula_str 中的项传递，告诉 align_cds () 在降维、聚类和轨迹推断之前减去这些信号。注意，您可以使用 alignment_group 、 residual_model_formula 或两者调用 align_cds() 。

# 降维和可视化结果

接下来，我们降低数据的维度。然而，与使用 UMAP 和 t-SNE 都能很好工作的聚类不同，在这里我们强烈建议您使用 UMAP，这是默认方法：

cds <- reduce_dimension(cds)

plot_cells(cds, label_groups_by_cluster=FALSE,  color_cells_by = "cell.type")

如您所见，尽管我们只查看了这个数据集的一小部分， Monocle 重建了一个带有众多分支的轨迹。将手动注释覆盖在 UMAP 上显示，这些分支主要被一种细胞类型占据。

与聚类分析一样，您可以使用 plot_cells() 来可视化单个基因沿轨迹的变化情况。让我们看看在纤毛神经元中表达模式有趣的一些基因：

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)

这样，我们可以看到这些特定基因在纤毛神经元中的表达模式，以及它们在细胞轨迹中如何变化。

# 聚类你的细胞

尽管细胞可能会连续地从一个状态过渡到下一个状态，没有它们之间的明确边界，但 Monocle 并不假设数据集中的所有细胞都来自一个共同的转录 “祖先”。在许多实验中，实际上可能存在多个不同的轨迹。例如，在对感染做出反应的组织中，组织驻留免疫细胞和基质细胞将具有非常不同的初始转录组，并且对感染的反应也大不相同，因此它们应该是同一轨迹的一部分。

Monocle 能够通过其聚类程序学习何时应该将细胞放在同一轨迹中，而不是分开的轨迹中。回想一下，我们运行 cluster_cells() 时，每个细胞不仅被分配到一个簇中，还被分配到一个分区中。当您正在学习轨迹时，每个分区最终将成为一个单独的轨迹。我们像以前一样运行 cluster_cells() 。

cds <- cluster_cells(cds)

plot_cells(cds, color_cells_by = "partition")

这样，我们可以看到每个细胞被分配到哪个分区中，并且这些分区如何对应于不同的轨迹。

# 学习轨迹图

接下来，我们将使用 learn_graph() 函数在每个分区内拟合一个主要图：

cds <- learn_graph(cds)

plot_cells(cds,

color_cells_by = "cell.type",

label_groups_by_cluster=FALSE,

label_leaves=FALSE,

label_branch_points=FALSE)

这个图将被用于许多后续步骤，如分支分析和差异表达。通过这种方式，我们可以更加深入地理解不同细胞类型在生物过程中的动态变化和轨迹。

这张图将使用于下游分析。

# 按拟时序排序细胞

一旦我们学习了一个图，我们就准备好根据细胞在发育程序中的进展来排序细胞。Monocle 用拟时序来衡量这一进展。下面的框定义了拟时序。

为了对细胞进行排序，我们需要告诉 Monocle 生物过程的 “起点” 在哪里。我们通过选择我们标记为轨迹 “根” 的图中的区域来实现这一点。在时间序列实验中，这通常可以通过找到 UMAP 空间中被早期时间点的细胞占据的位置来完成：

plot_cells(cds,

color_cells_by = "embryo.time.bin",

label_cell_groups=FALSE,

label_leaves=TRUE,

label_branch_points=TRUE,

graph_label_size=1.5)

黑线显示了图的结构。注意图并不是完全连接的：不同分区的细胞位于图的不同组成部分中。图中带数字的圆圈表示图中的特殊点。每个叶子（由浅灰色圆圈表示）对应于轨迹的不同结果（即细胞命运）。黑色圆圈表示分支节点，细胞可以走向多个结果。您可以使用 plot_cells 的 label_leaves 和 label_branch_points 参数控制这些是否在图中显示。请注意，圆圈中的数字仅供参考。

现在我们知道早期细胞所在的位置，我们可以调用 order_cells() ，它将计算每个细胞在拟时序中的位置。为此， order_cells() 需要您指定轨迹图的根节点。如果您没有作为参数提供它们，它将启动一个图形用户界面来选择一个或多个根节点。

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)

注意，有些细胞是灰色的。这意味着它们有无限的拟时序，因为它们无法从所选的根节点到达。通常，任何没有根节点的分区上的细胞都会被分配一个无限的拟时序。通常，您应该为每个分区选择至少一个根。

通常希望以编程方式指定轨迹的根，而不是手动选择它。下面的函数通过首先将细胞根据它们最接近的轨迹图节点进行分组，然后计算每个节点上来自最早时间点的细胞的比例，然后选择最多早期细胞占据的节点，并将其作为根返回。

# 一个辅助函数来识别根主要点：

get_earliest_principal_node <- function(cds, time_bin="130-170"){

    cell_ids <- which(colData(cds)[, "embryo.time.bin"] == time_bin)

    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))

通过编程选择的根节点传递给 order_cells() 会产生以下结果：

plot_cells(cds,

color_cells_by = "pseudotime",

label_cell_groups=FALSE,

label_leaves=FALSE,

label_branch_points=FALSE,

graph_label_size=1.5)

注意，我们可以通过首先使用 partitions() 函数按分区对细胞进行分组，然后在每个分区的基础上这样做。这将导致所有细胞被分配一个有限的拟时序。

# 根据分支子集细胞

根据轨迹中的分支对细胞进行子集化通常很有用。函数 choose_graph_segments 允许您交互式地这样做。

cds_sub <- choose_graph_segments(cds)

# 使用 3D 轨迹可视化

cds_3d <- reduce_dimension(cds, max_components = 3)

cds_3d <- cluster_cells(cds_3d)

cds_3d <- learn_graph(cds_3d)

cds_3d <- order_cells(cds_3d, root_pr_nodes=get_earliest_principal_node(cds))

cds_3d_plot_obj <- plot_cells_3d(cds_3d, color_cells_by="partition")

# 寻找拟时序相关基因

我们如何找到在轨迹不同路径上差异表达的基因？我们如何找到那些仅限于轨迹开始的基因？或者从轨迹中排除的基因？

再次，我们使用 graph_test() ，这次传递参数 neighbor_graph="principal_graph" ，这告诉它测试轨迹上相似位置的细胞是否有相关的表达：

ciliated_cds_pr_test_res <- graph_test(cds, neighbor_graph="principal_graph", cores=4)

pr_deg_ids <- row.names(subset(ciliated_cds_pr_test_res, q_value < 0.05))

这里有一些根据 graph_test() 得分非常显著的有趣基因：

plot_cells(cds, genes=c("hlh-4", "gcy-8", "dac-1", "oig-8"),

show_trajectory_graph=FALSE,

label_cell_groups=FALSE,

label_leaves=FALSE)

像以前一样，我们可以将轨迹变量基因收集到模块中：

gene_module_df <- find_gene_modules(cds[pr_deg_ids,], resolution=c(10^seq(-6,-1)))

这里我们绘制了 Packer & Zhu 等人注释的每组细胞类型的聚合模块得分：

cell_group_df <- tibble::tibble(cell=row.names(colData(cds)),

cell_group=colData(cds)$cell.type)

agg_mat <- aggregate_gene_expression(cds, gene_module_df, cell_group_df)

row.names(agg_mat) <- stringr::str_c("Module ", row.names(agg_mat))

pheatmap::pheatmap(agg_mat,

scale="column", clustering_method="ward.D2")

、

我们也可以将 gene_module_df 传递给 plot_cells() ，就像我们在上面比较 L2 数据中的簇时所做的。

plot_cells(cds,

genes=gene_module_df %>% filter(module %in% c(27, 10, 7, 30)),

label_cell_groups=FALSE,

show_trajectory_graph=FALSE)

Monocle 提供了另一个绘图函数，有时可以更清晰地显示沿单一路径的基因动态。您可以使用 choose_cells() 选择路径，或者通过按簇、细胞类型或其他仅限于路径的注释来子集化细胞数据集。让我们选择这样一个路径，即 AFD 细胞：

AFD_genes <- c("gcy-8", "dac-1", "oig-8")

AFD_lineage_cds <- cds[rowData(cds)$gene_short_name %in% AFD_genes,

colData(cds)$cell.type %in% c("AFD")]

AFD_lineage_cds <- order_cells(AFD_lineage_cds)

函数 plot_genes_in_pseudotime() 获取一小组基因，并根据拟时序显示它们的动态：

plot_genes_in_pseudotime(AFD_lineage_cds,

color_cells_by="embryo.time.bin",

min_expr=0.5)

您可以看到 dac-1 在其他两个基因之前被激活。

# 工作环境

sessionInfo()

R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Linux Mint 21.2

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=zh_CN.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=zh_CN.UTF-8        LC_COLLATE=zh_CN.UTF-8    
 [5] LC_MONETARY=zh_CN.UTF-8    LC_MESSAGES=zh_CN.UTF-8   
 [7] LC_PAPER=zh_CN.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=zh_CN.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] Seurat_5.0.1                SeuratObject_5.0.1         
 [3] sp_2.1-0                    dplyr_1.1.3                
 [5] monocle3_1.3.4              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.4            MatrixGenerics_1.6.0       
[13] matrixStats_1.0.0           Biobase_2.54.0             
[15] BiocGenerics_0.40.0        

loaded via a namespace (and not attached):
  [1] Rtsne_0.16             minqa_1.2.6            colorspace_2.1-0      
  [4] deldir_1.0-9           ellipsis_0.3.2         ggridges_0.5.4        
  [7] XVector_0.34.0         RcppHNSW_0.5.0         spatstat.data_3.0-1   
 [10] rstudioapi_0.15.0      leiden_0.4.3           listenv_0.9.0         
 [13] ggrepel_0.9.3          RSpectra_0.16-1        fansi_1.0.4           
 [16] codetools_0.2-18       splines_4.1.2          knitr_1.45            
 [19] polyclip_1.10-6        spam_2.9-1             jsonlite_1.8.7        
 [22] nloptr_2.0.3           ica_1.0-3              cluster_2.1.2         
 [25] png_0.1-8              uwot_0.1.16            spatstat.sparse_3.0-2 
 [28] sctransform_0.4.1      shiny_1.7.5            compiler_4.1.2        
 [31] httr_1.4.7             lazyeval_0.2.2         Matrix_1.6-4          
 [34] fastmap_1.1.1          cli_3.6.1              later_1.3.1           
 [37] htmltools_0.5.6        tools_4.1.2            igraph_1.5.1          
 [40] dotCall64_1.0-2        gtable_0.3.4           glue_1.6.2            
 [43] GenomeInfoDbData_1.2.7 reshape2_1.4.4         RANN_2.6.1            
 [46] Rcpp_1.0.11            scattermore_1.2        vctrs_0.6.3           
 [49] spatstat.explore_3.2-3 nlme_3.1-155           progressr_0.14.0      
 [52] lmtest_0.9-40          spatstat.random_3.1-6  stringr_1.5.0         
 [55] xfun_0.40              globals_0.16.2         lme4_1.1-34           
 [58] mime_0.12              miniUI_0.1.1.1         lifecycle_1.0.3       
 [61] irlba_2.3.5.1          goftest_1.2-3          future_1.33.0         
 [64] terra_1.7-46           zlibbioc_1.40.0        MASS_7.3-55           
 [67] zoo_1.8-12             scales_1.2.1           spatstat.utils_3.0-3  
 [70] promises_1.2.1         parallel_4.1.2         RColorBrewer_1.1-3    
 [73] yaml_2.3.7             gridExtra_2.3          pbapply_1.7-2         
 [76] reticulate_1.34.0      ggplot2_3.4.4          stringi_1.7.12        
 [79] fastDummies_1.7.3      boot_1.3-28            rlang_1.1.1           
 [82] pkgconfig_2.0.3        bitops_1.0-7           evaluate_0.22         
 [85] lattice_0.20-45        tensor_1.5             purrr_1.0.2           
 [88] ROCR_1.0-11            htmlwidgets_1.6.2      patchwork_1.1.3       
 [91] cowplot_1.1.1          tidyselect_1.2.0       parallelly_1.36.0     
 [94] RcppAnnoy_0.0.21       plyr_1.8.9             magrittr_2.0.3        
 [97] R6_2.5.1               generics_0.1.3         DelayedArray_0.20.0   
[100] pillar_1.9.0           fitdistrplus_1.1-11    abind_1.4-5           
[103] survival_3.2-13        RCurl_1.98-1.12        tibble_3.2.1          
[106] future.apply_1.11.0    KernSmooth_2.23-20     utf8_1.2.3            
[109] spatstat.geom_3.2-5    plotly_4.10.2          rmarkdown_2.25        
[112] grid_4.1.2             data.table_1.14.10     digest_0.6.33         
[115] xtable_1.8-4           tidyr_1.3.0            httpuv_1.6.11         
[118] munsell_0.5.0          viridisLite_0.4.2