?
前面一篇文章介紹了monocle2使碾，今天介紹一下monocle3。其實拘鞋，Monocle3和Monocle2并沒有本質(zhì)上的區(qū)別矢门，只是把降維圖從DDRTree改成了UMAP。原因可能是包的作者認(rèn)為UMAP比DDRTree降維更能反映高維空間的數(shù)據(jù)隔躲。
monocle3的官網(wǎng)地址：https://cole-trapnell-lab.github.io/monocle3/
其中宣旱，monocle3的三個主要功能：
1. 分群叛薯、計數(shù)細胞
2. 構(gòu)建細胞軌跡
3. 差異表達分析
工作流程如下圖：
===安裝===
if(!requireNamespace("BiocManager", quietly = TRUE))
? ? ? ? ? ?install.packages("BiocManager")
BiocManager::install(version ="3.10")
BiocManager::install(c('BiocGenerics','DelayedArray', 'DelayedMatrixStats',
?????????????????????? 'limma', 'S4Vectors','SingleCellExperiment',
?????????????????????? 'SummarizedExperiment','batchelor', 'Matrix.utils'))
install.packages("devtools")
devtools::install_github('cole-trapnell-lab/leidenbase')
devtools::install_github('cole-trapnell-lab/monocle3')
===數(shù)據(jù)準(zhǔn)備====
我們還使用pbmc3k的數(shù)據(jù)集案训，由于pbmc都是分化成熟的免疫細胞强霎，理論上并不存在直接的分化關(guān)系蓉冈，因此不適合用來做擬時軌跡分析轩触，這里仍然僅作為學(xué)習(xí)演示。
library(Seurat)
library(monocle3)
library(tidyverse)
library(patchwork)
rm(list=ls())
dir.create("Monocle3")
setwd("Monocle3")
?
//創(chuàng)建CDS對象并預(yù)處理數(shù)據(jù)
pbmc <- readRDS("pbmc.rds")
?
data <- GetAssayData(pbmc, assay ='RNA', slot = 'counts')
cell_metadata <- pbmc@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)
===預(yù)處理====
1. 標(biāo)準(zhǔn)化和PCA降維
?（RNA-seq是使用PCA家夺，如果是處理ATAC-seq的數(shù)據(jù)用Latent Semantic Indexing)
//preprocess_cds函數(shù)相當(dāng)于seurat中NormalizeData+ScaleData+RunPCA
cds <- preprocess_cds(cds, num_dim = 50)
plot_pc_variance_explained(cds)
注：圖中橫軸的50個PC就是preprocess_cds()中num_dim 設(shè)置的50個PC脱柱，如果這里設(shè)置的100，圖的橫軸也會展示100個PC拉馋。
2.?可視化
?
umap降維
cds <- reduce_dimension(cds,preprocess_method = "PCA")?
plot_cells(cds)
color_cells_by參數(shù)設(shè)置umap圖的顏色榨为，可以是colData(cds)中的任何一列。
cds <- cluster_cells(cds)?
colnames(colData(cds))
//以之前的Seurat分群來添加顏色煌茴，和原有的Seurat分群對比
p1 <- plot_cells(cds, reduction_method="UMAP", color_cells_by="seurat_clusters") + ggtitle('cds.umap')
//從seurat導(dǎo)入整合過的umap坐標(biāo)
cds.embed <- cds@int_colData$reducedDims$UMAP
int.embed <- Embeddings(pbmc, 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="seurat_clusters") + ggtitle('int.umap')
p = p1|p2
注：左圖是monocle3重新降維的結(jié)果，右圖是之前seurat降維的結(jié)果蔓腐。
如果細胞數(shù)目特別多（>10,000細胞或更多）矩乐，可以設(shè)置一些參數(shù)來加快UMAP運行速度。在reduce_dimension()函數(shù)中設(shè)置umap.fast_sgd=TRUE可以使用隨機梯度下降方法（fast stochastic gradient descent method）加速運行回论。還可以使用cores參數(shù)設(shè)置多線程運算散罕。
可視化指定基因：
ciliated_genes <- c("CD4","CD52","JUN")
plot_cells(cds,genes=ciliated_genes,label_cell_groups=FALSE,show_trajectory_graph=FALSE)
也可以使用tSNE降維：
cds <- reduce_dimension(cds, reduction_method="tSNE",preprocess_method = "PCA")
plot_cells(cds, reduction_method="tSNE", color_cells_by="seurat_clusters")
隨后也可使用Monocle3分cluster，鑒定每個cluster的marker基因并進行細胞注釋等等傀蓉。由于在Seurat的操作中已經(jīng)對數(shù)據(jù)進行了注釋欧漱，就不再使用Monocle3進行這些操作。
====構(gòu)建細胞軌跡===
軌跡學(xué)習(xí)Learn the trajectory graph（使用learn_graph()函數(shù)）
## 識別軌跡
cds <- learn_graph(cds)
plot_cells(cds, color_cells_by = "cell_type", label_groups_by_cluster=FALSE,label_leaves=FALSE, label_branch_points=FALSE)
上面這個圖將被用于許多下游分析葬燎，比如分支分析和差異表達分析误甚。
plot_cells(cds, color_cells_by = "cell_type", label_groups_by_cluster=FALSE,label_leaves=TRUE, label_branch_points=TRUE,graph_label_size=3)
注：黑色的線顯示的是graph的結(jié)構(gòu)。數(shù)字帶白色圓圈表示不同的結(jié)局萨蚕，也就是葉子靶草。數(shù)字帶黑色圓圈代表分叉點，從這個點開始岳遥，細胞可以有多個結(jié)局奕翔。這些數(shù)字可以通過label_leaves和label_branch_points參數(shù)設(shè)置。
2. 細胞按擬時排序
?
在學(xué)習(xí)了graph之后浩蓉，我們就可以根據(jù)學(xué)習(xí)的發(fā)育軌跡（擬時序）排列細胞派继。
為了對細胞進行排序，我們首先需要告訴Monocle哪里是這個過程的起始點捻艳。也就是需要指定軌跡的'roots'驾窟。
?
手動選擇root
cds <- order_cells(cds)
運行上面的代碼后，會跳出這個窗口认轨∩鹇纾可以手動在圖上選擇一個位置，然后點擊Done。（比如圖上我選擇的是右邊開始的位置的點）可以選擇多個位置恩急。
plot_cells(cds, color_cells_by = "pseudotime", label_cell_groups = FALSE, label_leaves = FALSE,? label_branch_points = FALSE)
saveRDS(cds, file = "cds.rds")
===差異表達分析=====
?尋找擬時軌跡差異基因
//graph_test分析最重要的結(jié)果是莫蘭指數(shù)（morans_I）杉畜，其值在-1至1之間，0代表此基因沒有
//空間共表達效應(yīng)衷恭，1代表此基因在空間距離相近的細胞中表達值高度相似此叠。
Track_genes <- graph_test(cds,neighbor_graph="principal_graph", cores=6)
Track_genes <-Track_genes[,c(5,2,3,4,1,6)] %>% filter(q_value < 1e-3)
write.csv(Track_genes,"Trajectory_genes.csv", row.names = F)
2. 挑選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="seurat_clusters",min_expr=0.5, ncol= 2)
FeaturePlot圖
p <- plot_cells(cds,genes=Track_genes_sig, show_trajectory_graph=FALSE,
??????????????? label_cell_groups=FALSE,? label_leaves=FALSE)
p$facet$params$ncol <- 5
ggsave("Genes_Featureplot.pdf",plot = p, width = 20, height = 8)
尋找共表達基因模塊
Track_genes <-read.csv("Trajectory_genes.csv")
genelist <- pull(Track_genes,gene_short_name) %>% as.character()
gene_module <-find_gene_modules(cds[genelist,], resolution=1e-1, cores = 6)
write.csv(gene_module,"Genes_Module.csv", row.names = F)
cell_group <-tibble::tibble(cell=row.names(colData(cds)),
????????????????????????????cell_group=colData(cds)$seurat_clusters)
agg_mat <-aggregate_gene_expression(cds, gene_module, cell_group)
row.names(agg_mat) <-stringr::str_c("Module ", row.names(agg_mat))
p <- pheatmap::pheatmap(agg_mat,scale="column", clustering_method="ward.D2")
ggsave("Genes_Module.pdf", plot =p, width = 8, height = 8)
提取擬時分析結(jié)果返回seurat對象
pseudotime <- pseudotime(cds,reduction_method = 'UMAP')
pseudotime <-pseudotime[rownames(pbmc@meta.data)]
pbmc$pseudotime <- pseudotime
p = FeaturePlot(pbmc, reduction ="umap", features = "pseudotime")
ggsave("Pseudotime_Seurat.pdf",plot = p, width = 8, height = 6)
saveRDS(pbmc, file ="sco_pseudotime.rds")
本文使用 文章同步助手 同步
