单细胞分析之monocle3

2014/07/06来源：单细胞天地评论1,555

monocle3简介

monocel3的优势

从UMAP图识别发育轨迹，可以继承Seurat的质控、批次校正和降维分析结果，实现“一张图”展现细胞的聚类、鉴定和轨迹分析结果。
自动对UMAP图分区（partition），可以选择多个起点，轨迹分析算法的逻辑更符合生物学现实。

除了轨迹分析的主要功能，monocle3差异分析方法也有其独到之处，可以做一些与seurat不好实现的分析。

monocel3的安装

先安装一些依赖包，大家安装前可以查看一下这些包是否已经安装过了。

BiocManager::install(c('BiocGenerics', 'DelayedArray', 'DelayedMatrixStats',
                       'limma', 'S4Vectors', 'SingleCellExperiment',
                       'SummarizedExperiment', 'batchelor', 'Matrix.utils'))

然后安装实现umap图分区的包leidenbase，最后安装monole3

install.packages("devtools")
devtools::install_github('cole-trapnell-lab/leidenbase')
devtools::install_github('cole-trapnell-lab/monocle3')

安装有困难的朋友可以使用我的镜像kinesin/rstudio:1.2，下载链接见《kinesin_rstudio的日常升级二》，使用方法见《华为云配置单细胞分析环境及报错处理》。

monocle3分析实践

数据来源

数据来自Immune Landscape of Viral- and Carcinogen-Driven Head and Neck Cancer，数据集GEO编号：GSE139324。建议大家自己下载。原数据集有63个scRNA的数据，都是分选的CD45 免疫细胞。考虑到计算资源问题，挑选了10个样本用于此次演示。

最后用于轨迹分析的细胞是提取的T细胞亚群。

创建seurat对象

library(Seurat)
library(monocle3)
library(tidyverse)
library(patchwork)
rm(list=ls())
 
##==准备seurat对象列表==##
dir <- dir("GSE139324/")           #GSE139324是存放数据的目录
dir <- paste0("GSE139324/",dir)
sample_name <- c('HNC01PBMC', 'HNC01TIL', 'HNC10PBMC', 'HNC10TIL', 'HNC20PBMC',
                 'HNC20TIL',  'PBMC1', 'PBMC2', 'Tonsil1', 'Tonsil2')
scRNAlist <- list()
for(i in 1:length(dir)){
counts <- Read10X(data.dir = dir[i])
scRNAlist[[i]] <- CreateSeuratObject(counts, project=sample_name[i], min.cells=3, min.features = 200)
scRNAlist[[i]][["percent.mt"]] <- PercentageFeatureSet(scRNAlist[[i]], pattern = "^MT-")
}   
scRNA <- merge(scRNAlist[[1]], scRNAlist[2:10])
##提取表达矩阵用于细胞类型鉴定
count <- GetAssayData(scRNA, assay = "RNA", slot = "counts")
saveRDS(count, "count.rds")

利用azimuth鉴定细胞类型

Seurat官网近期推出了在线细胞类型鉴定服务，可以准确鉴定pbmc细胞，建议大家尝试一下。

点击箭头所指的Browse按钮，将上一步保存的count.rds文件上传到网站，网址：http://azimuth.satijalab.org/app/azimuth

可以按自己的需要设置质控标准，网页会同步展示小提琴图

质控参数确定后，点击Map cells to reference就可以鉴定细胞类型了。

比对成功后，先点击左边红框标注的Download Results，然后下载预测结果

提取T细胞亚群

## 读取azimuth鉴定的结果
pred <- read.delim("azimuth_pred.tsv")
head(pred, 2)
#                          cell predicted.id predicted.score mapping.score
#1 HNC01PBMC_AAACCTGAGGAGCGTT-1    CD14 Mono       1.0000000     0.8842671
#2 HNC01PBMC_AAACCTGCACGGACAA-1      CD8 TEM       0.3651781     0.8838794
## 提取只含T细胞的子集
pred.T <- subset(pred, pred$predicted.id %in% c('CD4 CTL',
                                                'CD4 Naive', 
                                                'CD4 Proliferating',
                                                'CD4 TCM', 
                                                'CD4 TEM', 
                                                'CD8 Naive', 
                                                'CD8 TCM',
                                                'CD8 Proliferating', 
                                                'CD8 TEM',
                                                'MAIT', 
                                                'Treg',
                                                'dnT', 
                                                'gdT'))
scRNAsub <- scRNA[,as.character(pred.T$cell)]
## T细胞子集去除批次效应
scRNAsub <- SCTransform(scRNAsub) %>% RunPCA() 
scRNAsub <- RunHarmony(scRNAsub, group.by.vars = "orig.ident",
                       assay.use = "SCT", max.iter.harmony = 10)  
## 降维聚类
ElbowPlot(scRNA, ndims = 50)
pc.num=1:30
scRNAsub <- RunUMAP(scRNAsub, reduction="harmony", dims=pc.num) %>%
            FindNeighbors(reduction="harmony", dims=pc.num) %>% 
            FindClusters(resolution=0.8)   
pred.T <- data.frame(pred.T[, c(2,3,4)], row.names = pred.T$cell)
scRNAsub <- AddMetaData(scRNAsub, metadata = pred.T)           
DimPlot(scRNAsub, reduction = "umap", group.by = "predicted.id", label = T)   NoLegend()

monocle3轨迹分析

##创建CDS对象并预处理数据
data <- GetAssayData(scRNAsub, assay = 'RNA', slot = 'counts')
cell_metadata <- scRNAsub@meta.data
gene_annotation <- data.frame(gene_short_name = rownames(data))
rownames(gene_annotation) <- rownames(data)
cds <- new_cell_data_set(data,
                         cell_metadata = cell_metadata,
                         gene_metadata = gene_annotation)
#preprocess_cds函数相当于seurat中NormalizeData ScaleData RunPCA
cds <- preprocess_cds(cds, num_dim = 50)
#umap降维
cds <- reduce_dimension(cds, preprocess_method = "PCA")
p1 <- plot_cells(cds, reduction_method="UMAP", color_cells_by="celltype")   ggtitle('cds.umap')
##从seurat导入整合过的umap坐标
cds.embed <- cds@int_colData$reducedDims$UMAP
int.embed <- Embeddings(scRNAsub, reduction = "umap")
int.embed <- int.embed[rownames(cds.embed),]
cds@int_colData$reducedDims$UMAP <- int.embed
p2 <- plot_cells(cds, reduction_method="UMAP", color_cells_by="celltype")   ggtitle('int.umap')

monocle降维的结果

导入整合过的umap坐标后作图

## Monocle3聚类分区
cds <- cluster_cells(cds)
p1 <- plot_cells(cds, show_trajectory_graph = FALSE)   ggtitle("label by clusterID")
p2 <- plot_cells(cds, color_cells_by = "partition", show_trajectory_graph = FALSE)   
        ggtitle("label by partitionID")
p = wrap_plots(p1, p2)

单细胞分析之monocle3-图片1

## 识别轨迹
cds <- learn_graph(cds)
p = plot_cells(cds, label_groups_by_cluster = FALSE, label_leaves = FALSE, 
               label_branch_points = FALSE)

单细胞分析之monocle3-图片2

##细胞按拟时排序
# cds <- order_cells(cds) 存在bug，使用辅助线选择root细胞
p   geom_vline(xintercept = seq(-7,-6,0.25))   geom_hline(yintercept = seq(0,1,0.25))
embed <- data.frame(Embeddings(scRNAsub, reduction = "umap"))
embed <- subset(embed, UMAP_1 > -6.75 & UMAP_1 < -6.5 & UMAP_2 > 0.24 & UMAP_2 < 0.25)
root.cell <- rownames(embed)
cds <- order_cells(cds, root_cells = root.cell)
plot_cells(cds, color_cells_by = "pseudotime", label_cell_groups = FALSE, 
           label_leaves = FALSE,  label_branch_points = FALSE)

画几条辅助线用于找到root细胞

完成拟时分析的细胞排序结果

monocle3差异分析

##寻找拟时轨迹差异基因
#graph_test分析最重要的结果是莫兰指数（morans_I），其值在-1至1之间，0代表此基因没有
#空间共表达效应，1代表此基因在空间距离相近的细胞中表达值高度相似。
Track_genes <- graph_test(cds, neighbor_graph="principal_graph", cores=10)
#挑选top10画图展示
Track_genes_sig <- Track_genes %>% top_n(n=10, morans_I) %>%
                   pull(gene_short_name) %>% as.character()
#基因表达趋势图
plot_genes_in_pseudotime(cds[Track_genes_sig,], color_cells_by="predicted.id", 
                              min_expr=0.5, ncol = 2)
#FeaturePlot图
plot_cells(cds, genes=Track_genes_sig, show_trajectory_graph=FALSE,
               label_cell_groups=FALSE,  label_leaves=FALSE)
##寻找共表达模块
genelist <- pull(Track_genes, gene_short_name) %>% as.character()
gene_module <- find_gene_modules(cds[genelist,], resolution=1e-2, cores = 10)
cell_group <- tibble::tibble(cell=row.names(colData(cds)), 
                             cell_group=colData(cds)$predicted.id)
agg_mat <- aggregate_gene_expression(cds, gene_module, cell_group)
row.names(agg_mat) <- stringr::str_c("Module ", row.names(agg_mat))
pheatmap::pheatmap(agg_mat, scale="column", clustering_method="ward.D2")

基因表达趋势图

FeaturePlot图

共表达模块热图