Single Cell——轨迹分析（Monocle3）

monocle3

Monocle 3包提供了分析单细胞基因表达实验的工具包。

Monocle3和Monocle2并没有本质上的区别，只是把降维图从DDRTree改成了UMAP。原因可能是包的作者认为UMAP比DDRTree降维更能反映高维空间的数据。

monocle3的三个主要功能：

• 分群、计数细胞
• 构建细胞轨迹
• 差异表达分析

细胞聚类及鉴定亚群

安装包

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
# 一些依赖包
bioc_pkgs <- c('BiocGenerics', 'DelayedArray', 'DelayedMatrixStats',
                       'limma', 'S4Vectors', 'SingleCellExperiment',
                       'SummarizedExperiment')
for (bpkg in bioc_pkgs){
   
  if (! require(bpkg,character.only=T) ) {
   
    BiocManager::install(bpkg,ask = F,update = F)
    require(bpkg,character.only=T)
    }
}
 
# (可选)如果要在细胞聚类时设定分辨率参数(resolution)，就要安装一个python包
install.packages("reticulate")
reticulate::py_install("louvain")
# 现在安装3版本，还是要通过github
devtools::install_github('cole-trapnell-lab/monocle3')
# 重启Rstudio，检测是否成功
library(Seurat)
library(monocle3)
library(tidyverse)
library(patchwork)
library(ggplot2)
library(dplyr)

加载数据

• expression_matrix 是一个行为基因，列为细胞样本的表达矩阵
• cell_metadata 是一个数据框，行为细胞，列是细胞的属性（比如细胞类型、培养环境、培养时间等）
• gene_metadata是一个数据框，行为feature信息（比如基因），列是基因属性（比如GC含量）

必须要满足：

1. expression_matrix的行与gene_metadata的行对应
2. expression_matrix的列与cell_metadata的行对应

# 下载数据(如果已经有了，可以跳过)
expression_matrix <- readRDS(url("http://staff.washington.edu/hpliner/data/cao_l2_expression.rds"))
cell_metadata <- readRDS(url("http://staff.washington.edu/hpliner/data/cao_l2_colData.rds"))
gene_annotation <- readRDS(url("http://staff.washington.edu/hpliner/data/cao_l2_rowData.rds"))
 
# 加载、探索数据
expression_matrix <- readRDS(file = 'cao_l2_expression.rds')
cell_metadata <- readRDS(file = 'cao_l2_colData.rds')
gene_annotation <- readRDS(file = 'cao_l2_rowData.rds')
 
dim(expression_matrix)
expression_matrix[1:3,1:3]
cell_metadata[1:3,1:3]
head(gene_annotation)

将数据存储在cell_data_set对象中

# 创建CDS(Cell Data Set)对象
cds <- new_cell_data_set(expression_matrix,
                                   cell_metadata = cell_metadata,
                                   gene_metadata = gene_annotation)
cds 

数据预处理

# 对数据进行规范化和预处理
# 对数据进行规范化和预处理
cds <- preprocess_cds(cds, num_dim = 100) #大约2 mins
plot_pc_variance_explained(cds)

我们可以看到，使用超过100台PC只会捕获少量额外的变化，并且每增加一台PC都会使Monocle的下游步骤变慢。

降低维度并可视化细胞

可以使用t-SNE，这在单细胞rna测序中非常流行，或者使用UMAP。Monocle 3默认使用UMAP，它更适合于RNA-seq中的聚类和轨迹分析。调用reduce_dimension()将数据的维数降低到X， Y平面，以便我们可以轻松地绘制它，请：

# umap降维
cds <- reduce_dimension(cds, preprocess_method = "PCA",reduction_method = c("UMAP"))

plot_cells(cds)

图中的每个点都代表cds对象中的不同细胞，可以看到不同的细胞分成了不同的群，有的群有几千个细胞，有的群只有少数几个，我们这里直接使用测试数据中做好的细胞标记

# 数据中一共提供了30类细胞(细胞类型及数量见：)
table(colData(cds)$"cao_cell_type")
# 根据这个上色
plot_cells(cds, color_cells_by="cao_cell_type")
# (另外除了用cao_cell_type细胞类型，还能用下面的任意一个)
colnames(colData(cds))

可以看到图中许多细胞类型再UMAP结果中靠的很近

根据基因表达量进行映射：

#根据基因表达量进行映射
plot_cells(cds, genes=c("cpna-2", "egl-21", "ram-2", "inos-1"))

使用tSNE方法降维的话：

#进行tSNE降维
cds <- reduce_dimension(cds, reduction_method="tSNE")
# 然后对tSNE结果可视化
plot_cells(cds, reduction_method="tSNE", color_cells_by="cao_cell_type")

其实可以看到，tSNE结果的各个细胞亚群内部的关联性不如UMAP

检查、移除批次效应 => residual_model_formula_str参数

为了更加真实地反映差异基因的来源是生物学因素，而不是技术因素（比如细胞板批次、细胞培养差别、测序差别等），就首先要检查有没有所谓的”技术批次“对数据差异产生影响。

当进行降维分析时，应该经常检查批次效应，最好在colData(cds) 添加一列，其中注明细胞的批次信息（比如来自哪个细胞板），然后就能根据这个对细胞的批次进行可视化，结果最好是各个细胞亚群的不同批次混在一起，这样我们就能认为：降维所采用的”特异性“的特征并不来自批次

#对细胞的批次进行可视化
plot_cells(cds, color_cells_by="plate", label_cell_groups=FALSE)

如果真的担心不同批次会产生其他的影响，可以在降维之前对批次进行校正，保证后面的各个细胞群只来自一个批次，那么它们之间的分群差异就更可能是由于生物因素导致

# 假入要对批次进行校正
cds = preprocess_cds(cds, num_dim = 100, residual_model_formula_str = "~ plate")
# 校正后再降维
cds = reduce_dimension(cds)
# 再次检查批次效应
plot_cells(cds, color_cells_by="plate", label_cell_groups=FALSE)

细胞聚类

这是细胞分群过程中非常重要的一步。Monocle利用Louvain community detection 这个非监督聚类方法（它是phenoGraph算法的一部分），它和我们常见的Kmeans、hclust层次聚类等还不太一样，这个方法更倾向于找到一个网络中联系紧密的部分，而经常忽略节点的特性；而我们常见的聚类呢，它一般会忽视整体中各个部分的联系，而通过计算两个节点(目标)之间的距离（如欧式距离、曼哈顿距离、余弦相似度等） 找到相似的特性

按照cluster可视化

cds = cluster_cells(cds, resolution=c(10^seq(-6,-1)))
plot_cells(cds)
# 这个可视化函数如果不加任何参数，它默认对不同的cluster上色
# 用下面函数查看有哪些cluster：
cds@clusters$UMAP$clusters #或者
clusters(cds) 

按照partition进行可视化

不过这样细分看的太混杂，于是 Alex Wolf et al 利用他们开发的 PAGA 算法将一些小的cluster聚合成一个大的partition

plot_cells(cds, color_cells_by="partition", group_cells_by="partition")
# 检查有多少partition
cds@clusters$UMAP$partitions

按照每个cluster对应的细胞类型

除了默认按照cluster编号进行标记不同的亚群，还可以按照每个cluster对应的细胞类型进行可视化：

#按照每个cluster对应的细胞类型进行可视化
plot_cells(cds, color_cells_by=