寫在前面：今天突然想寫這個(gè)教程馆类，花了一下午看資料歼培，其實(shí)我前面是有一些基礎(chǔ)的，但是看的還是會很亂季率，所以想把這個(gè)包的脈絡(luò)寫一下野瘦，希望能夠可以跟大家討論一下并且自己以后使用這個(gè)包的時(shí)候也相當(dāng)于翻閱說明書吧！ｌ骸鞭光！
首先講一段話：
我的一個(gè)使用情景是，當(dāng)我跑完seruat的標(biāo)準(zhǔn)流程后泞遗，不可能只看細(xì)胞數(shù)量之間的差異吧惰许，所以想看細(xì)胞之間的發(fā)育軌跡，就想到了monocle史辙，Attention：我是已經(jīng)完成了seurat分析Ｐ诼颉佩伤！
所以大家能從上面這段話得到什么信息呢？信息：就是我只想用monocle對我的細(xì)胞進(jìn)行軌跡分析卦睹，然而monocle這個(gè)包功能其實(shí)不單單是進(jìn)行軌跡分析畦戒，它其實(shí)一直想要代替seruat這個(gè)包方库，所以它幾乎包括seurat所有的分析流程
假如你是沒有進(jìn)行seurat分析的話结序，你可以嘗試一下用monocle進(jìn)行分析
#加載包
library(Seurat)
library(monocle3)
library(tidyverse)
library(patchwork)
#導(dǎo)入數(shù)據(jù)
#在這里你需要導(dǎo)入三個(gè)數(shù)據(jù)：細(xì)胞的表達(dá)矩陣，meta信息纵潦，gene_annotation
data <- 細(xì)胞的表達(dá)矩陣
cell_metadata <- meta信息
gene_annotation <- data.frame(gene_short_name = rownames(data))
rownames(gene_annotation) <- rownames(data)
#創(chuàng)建CDS對象
cds <- new_cell_data_set(data,
cell_metadata = cell_metadata,
gene_metadata = gene_annotation) #這樣就獲得了monocle的對象了徐鹤，可類比于創(chuàng)建seurat對象
#思考一下我們在創(chuàng)建了seurat對象后要干嘛呢？降維聚類分群邀层，那么monocle也不例外
#NormalizeData+ScaleData+RunPCA
cds <- preprocess_cds(cds, num_dim = 50) #preprocess_cds函數(shù)相當(dāng)于seurat中NormalizeData+ScaleData+RunPCA
plot_pc_variance_explained(cds) #像seurat一樣展示pc數(shù)
#umap降維
cds <- reduce_dimension(cds,preprocess_method = "PCA") #preprocess_method默認(rèn)是PCA
#tSNE降維
cds <- reduce_dimension(cds, reduction_method="tSNE")
#聚類
cds <- cluster_cells(cds)
#那么我們到這里就完成了分群了返敬，后面就是對于每一個(gè)群的細(xì)胞注釋和細(xì)胞類型的結(jié)果展示了，這個(gè)大家可以去看周運(yùn)來老師的簡書寥院，而我的目的不是這個(gè)
還記得我提出來的background嗎劲赠，我已經(jīng)用seurat進(jìn)行降維聚類分群和細(xì)胞注釋，所以我沒有必要monocle包給我做這些秸谢，對于使用了seurat分析流程的我們只想用monocle來做一個(gè)擬時(shí)序分析該怎么做呢凛澎？
#加載包
library(Seurat)
library(monocle3)
library(tidyverse)
library(patchwork)
rm(list=ls())
#導(dǎo)入我們的用seurat包分析得到的S4對象
pbmc <- readRDS("scRNAsub.rds")
#容易弄混的地方來了
#首先我們要得到seurat對象的表達(dá)矩陣
data <- GetAssayData(pbmc, assay = 'RNA', slot = 'counts')
#其次獲得meta信息
cell_metadata <- pbmc@meta.data
#然后要獲得gene_annotation
gene_annotation <- data.frame(gene_short_name = rownames(data))
rownames(gene_annotation) <- rownames(data)
#創(chuàng)建CDS對象
cds <- new_cell_data_set(data,
cell_metadata = cell_metadata,
gene_metadata = gene_annotation)
#為啥要這樣呢，這怎么解釋呢估蹄，你用monocle的包塑煎，所以要聽它的，還是要走一遍某些流程
#然后就是降維聚類分群
#NormalizeData+ScaleData+RunPCA
cds <- preprocess_cds(cds, num_dim = 50) #preprocess_cds函數(shù)相當(dāng)于seurat中NormalizeData+ScaleData+RunPCA
plot_pc_variance_explained(cds) #像seurat一樣展示pc數(shù)
#umap降維
cds <- reduce_dimension(cds,preprocess_method = "PCA") #preprocess_method默認(rèn)是PCA
#tSNE降維
cds <- reduce_dimension(cds, reduction_method="tSNE")
#聚類
cds <- cluster_cells(cds)
這里必須強(qiáng)調(diào)一下：剛開始我會糾結(jié)臭蚁，monocle的降維聚類分群的結(jié)果會不會跟seurat不一樣呢最铁，答案是肯定的，算法都不同垮兑，能一樣嗎冷尉，那豈不是改變了我們的初衷，只想做擬時(shí)序分析
那么解決辦法來了系枪，我們可以用seurat分析得到umap結(jié)果覆蓋掉monocle分析得到的umap結(jié)果
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 #因此我們以后只要是畫umap的點(diǎn)圖就跟seurat的圖片形狀是一樣的
#畫圖看一下
plot_cells(cds, reduction_method="UMAP", color_cells_by="seurat_clusters")
在這里我要講兩點(diǎn)：
第一：為什么一定要做降維聚類分群雀哨，因?yàn)槲覀兪怯胢onocle創(chuàng)建了一個(gè)全新的對象，我們只是導(dǎo)入了seurat對象里的表達(dá)矩陣嗤无，meta信息和genelist震束，所以這個(gè)cds是沒有進(jìn)行降維聚類等操作，導(dǎo)致后面的擬時(shí)序分析是做不了的当犯，假如你看了monocle的官網(wǎng)的話啊垢村，你會知道的，官網(wǎng)上說-擬時(shí)序分析是基于低維度嚎卫，所以必須對cds進(jìn)行降維聚類分群嘉栓。
第二：已經(jīng)用seurat分析過的我們宏榕，可以在monocle流程里省去哪一些步驟呢，1.我們省去了找
marker基因侵佃，2.我們省去了細(xì)胞注釋麻昼，但是我們加進(jìn)去的meta信息里有細(xì)胞類型和其他的一些信息，所以其實(shí)這一步可以因?yàn)槲覀冇昧藄eurat分析并輸入了它的meta信息而代替馋辈，所以我們還是相當(dāng)于做了完整的monocle流程抚芦，只是有一些是用seurat做的結(jié)果罷了（歡迎一起討論這一點(diǎn)）
下面開始才是我們想要的：
擬時(shí)序分析
第一步:軌跡推斷
cds <- learn_graph(cds) #是的，就這么簡單迈螟，用 learn_graph即可
#畫圖 我們可以根據(jù)meta信息來畫出相關(guān)的圖片
head(colData(cds))
plot_cells(cds,
color_cells_by = "celltype",
label_groups_by_cluster=FALSE,
label_leaves=FALSE,
label_branch_points=TRUE,
group_label_size=4,
cell_size=1.5) #這個(gè)圖就可以看出細(xì)胞直接的分化軌跡了
#黑色的線顯示的是graph的結(jié)構(gòu)叉抡。
#數(shù)字帶白色圓圈表示不同的結(jié)局，也就是葉子答毫。
#數(shù)字帶黑色圓圈代表分叉點(diǎn)褥民，從這個(gè)點(diǎn)開始，細(xì)胞可以有多個(gè)結(jié)局洗搂。
#這些數(shù)字可以通過label_leaves和label_branch_points參數(shù)設(shè)置消返。
第二步：選擇根 就是選擇細(xì)胞的起源 這也是monocle3與monocle2最大的區(qū)別，手動選擇(自定義)根 耘拇，使用其中一種就可以Ｄ旒铡！驼鞭！
#第一種：
cds = order_cells(cds) #在前陣子有bug 無法展示秦驯，現(xiàn)在好像又可以了
plot_cells(cds,
color_cells_by = "pseudotime",
label_cell_groups=FALSE,
label_leaves=TRUE,
label_branch_points=TRUE,
graph_label_size=1.5,
group_label_size=4,cell_size=1.5)
#第二種：用輔助線
p + geom_vline(xintercept = seq(-7,-6,0.25)) + geom_hline(yintercept = seq(0,1,0.25))
embed <- data.frame(Embeddings(pbmc, 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)
#第三種：定義一個(gè)函數(shù)
# a helper function to identify the root principal points:
get_earliest_principal_node <- function(cds, time_bin="epithelial cells"){
cell_ids <- which(colData(cds)[, "celltype"] == 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)) #很多人這里一直在哭訴error，那是都沒有理解這一步在干嘛挣棕，很多無腦運(yùn)行別人的代碼译隘，別人代碼選擇XX細(xì)胞作為根，你的數(shù)據(jù)集里又沒有洛心，所以報(bào)錯(cuò)說沒有node啦
第一種方法：
可以看出我選擇紅點(diǎn)的地方為發(fā)育起始點(diǎn)固耘，我用的是第一種方法
第二種方法：
上面的發(fā)育軌跡是在一個(gè)二維的坐標(biāo)上看的，monocle3搞了一個(gè)三維的發(fā)育軌跡词身，所以我們可以在三維空間內(nèi)看發(fā)育軌跡(其實(shí)我覺得也那樣厅目，無非是加多一個(gè)維度罷了） 但是假如你想畫三維的圖，你要定義一個(gè)function法严，就像第三種找root一樣损敷，定義一個(gè)能找root的函數(shù)，因?yàn)榫退闶侨S深啤，你也要給別人一個(gè)根呀拗馒，否則怎么給你畫順序呢？
#第一步：定義一個(gè)function
get_earliest_principal_node <- function(cds, time_bin="epithelial cells"){
cell_ids <- which(colData(cds)[, "celltype"] == 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
}
#第二步：降維到三維溯街，并plot
cds_3d = reduce_dimension(cds, max_components = 3) #這里就降到三維嘛
cds_3d = cluster_cells(cds_3d)
cds_3d = learn_graph(cds_3d) #為啥這慢呢诱桂，因?yàn)樗窃谌S的基礎(chǔ)上去進(jìn)行l(wèi)earn
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="celltype")
cds_3d_plot_obj #這個(gè)會由于網(wǎng)絡(luò)的原因畫不出來
到這里大家應(yīng)該能體會到了為什么前面不管你是不是seurat處理過的洋丐，都要對cds這個(gè)對象進(jìn)行降維了吧，不然法畫圖呀
接下來是差異分析挥等，就是找出哪一些基因在這個(gè)軌跡中變化最大
#計(jì)算基因按照軌跡的變化
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)
#挑選top10畫圖展示
Track_genes_sig <- Track_genes %>% top_n(n=10, morans_I) %>%
pull(gene_short_name) %>% as.character()
#基因表達(dá)趨勢圖
plot_genes_in_pseudotime(cds[Track_genes_sig,], color_cells_by="celltype",
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)
#在這里運(yùn)用到文章的話友绝，應(yīng)該可以這樣用，你可以用這些變化top基因來講一些東西肝劲，或者你可以show你自己fouce on的gene
然后就是還可以做一個(gè)尋找共表達(dá)模塊迁客，怎么應(yīng)用呢，就可以說是這個(gè)模塊導(dǎo)致他們狀態(tài)發(fā)生變化涡相，然后查看一個(gè)模塊里的基因哲泊，再講一些東西
##尋找共表達(dá)模塊
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)$celltype)
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")
最后剩蟀，我們可以把monocle的擬時(shí)序結(jié)果導(dǎo)入seurat對象里了
#提取擬時(shí)分析結(jié)果返回seurat對象
pseudotime <- pseudotime(cds, reduction_method = 'UMAP')
pseudotime <- pseudotime[rownames(pbmc@meta.data)]
pbmc$pseudotime <- pseudotime
p = FeaturePlot(pbmc, reduction = "umap", features = "pseudotime")
# pseudotime中有無限值催蝗，無法繪圖。
ggsave("Pseudotime_Seurat.pdf", plot = p, width = 8, height = 6)
saveRDS(pbmc, file = "sco_pseudotime.rds")
References:
單細(xì)胞分析十八般武藝5：monocle3 - 云+社區(qū) - 騰訊云 (tencent.com)
scRNA-seq數(shù)據(jù)分析 || Monocle3 - 簡書 (jianshu.com)
單細(xì)胞之軌跡分析-3：monocle3 - 簡書 (jianshu.com)
假如有講錯(cuò)或理解不到位的地方歡迎大家指出哈!!!
