我们知道monocle2 最重要的一个分析是对分支进行差异分析（如下图）
image.png
经典的教程总是：
BEAM_res <- BEAM(cds1, branch_point = 1, cores = 4)
BEAM_res <- BEAM_res[order(BEAM_res$qval),]
my_branched_heatmap <- plot_genes_branched_heatmap(cds1[row.names(subset(BEAM_res, qval < 1e-4)),],branch_point= 1,
num_clusters = 4,cores = 4,use_gene_short_name = TRUE,show_rownames = FALSE,return_heatmap = TRUE)
那如果我的分支并非只有一个而是多个，我怎么批量化分析呢？（如下图）
image.png
也就是说如何在不看图的情况下获得分支的个数呢？很简单，以下命令，敲重点：
mst_branch_nodes <- cds@auxOrderingData[[cds@dim_reduce_type]]$branch_points
image.png
这样你就可以循环对分支进行分析，代码如下：
for(i in 1:length(mst_branch_nodes)){
tryCatch({
print(Sys.time())
print("run BEAM")
BEAM_res <- BEAM(cds1, branch_point = i, cores = 4)
BEAM_res <- BEAM_res[order(BEAM_res$qval),]
BEAM_res <- BEAM_res[BEAM_res$qval<0.05,]
BEAM_res <- BEAM_res[,c("gene_short_name", "pval", "qval")]
write.csv(BEAM_res,paste(outdir,paste0(prefix,'_branch_point',i,'_genes.csv'),sep="/"),quote=F)
# check out the results
head(BEAM_res)
print(Sys.time())
graph <- paste(outdir,paste0(prefix,'_branch',i,'_heatmap.pdf'),sep='/')
print(graph)
pdf(graph, w=8, h=12)
my_branched_heatmap <- plot_genes_branched_heatmap(cds1[row.names(subset(BEAM_res, qval < 1e-4)),],branch_point = i,num_clusters = 4,cores = 4,use_gene_short_name = TRUE,show_rownames = FALSE,return_heatmap = TRUE) ###show_rownames = TRUE
dev.off()
my_row <- my_branched_heatmap$annotation_row
my_row <- data.frame(cluster = my_row$Cluster,gene = row.names(my_row),stringsAsFactors = FALSE)
head(my_row[my_row$cluster == 3,'gene'])
my_gene <- row.names(subset(fData(cds1),gene_short_name %in% head(my_row[my_row$cluster == 1,'gene'])))
graph <- paste(outdir,paste0(prefix,'genes_branch',i,'_pseudotime.pdf'),sep='/')
pdf(graph, w=6, h=12)
plot_genes_branched_pseudotime(cds1[my_gene,],branch_point = i,ncol = 2)
dev.off()
}, error = function(e) {
cat("monocle2 BEAM",conditionMessage(e),"\n\n")
})yu){
tryCatch({
print(Sys.time())
print("run BEAM")
BEAM_res <- BEAM(cds1, branch_point = i, cores = 4)
BEAM_res <- BEAM_res[order(BEAM_res$qval),]
BEAM_res <- BEAM_res[BEAM_res$qval<0.05,]
BEAM_res <- BEAM_res[,c("gene_short_name", "pval", "qval")]
write.csv(BEAM_res,paste(outdir,paste0(prefix,'_branch_point',i,'_genes.csv'),sep="/"),quote=F)
# check out the results
head(BEAM_res)
print(Sys.time())
graph <- paste(outdir,paste0(prefix,'_branch',i,'_heatmap.pdf'),sep='/')
print(graph)
pdf(graph, w=8, h=12)
my_branched_heatmap <- plot_genes_branched_heatmap(cds1[row.names(subset(BEAM_res, qval < 1e-4)),],branch_point = i,num_clusters = 4,cores = 4,use_gene_short_name = TRUE,show_rownames = FALSE,return_heatmap = TRUE) ###show_rownames = TRUE
dev.off()
my_row <- my_branched_heatmap$annotation_row
my_row <- data.frame(cluster = my_row$Cluster,gene = row.names(my_row),stringsAsFactors = FALSE)
head(my_row[my_row$cluster == 3,'gene'])
my_gene <- row.names(subset(fData(cds1),gene_short_name %in% head(my_row[my_row$cluster == 1,'gene'])))
graph <- paste(outdir,paste0(prefix,'genes_branch',i,'_pseudotime.pdf'),sep='/')
pdf(graph, w=6, h=12)
plot_genes_branched_pseudotime(cds1[my_gene,],branch_point = i,ncol = 2)
dev.off()
}, error = function(e) {
cat("monocle2 BEAM",conditionMessage(e),"\n\n")
})
网友评论