安裝
#安裝monocle2包
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("monocle")
#安裝nomocle依賴包
library(monocle)
install.packages("devtools")
devtools::install_github("cole-trapnell-lab/monocle-release@develop")
創(chuàng)建文件夾存儲(chǔ)結(jié)果
#創(chuàng)建一個(gè)文件夾用于寫分析結(jié)果
result.name <- "cca_monocle_result"
if(!dir.exists(result.name)){
dir.create(result.name)
}
用CCA后的數(shù)據(jù)進(jìn)行擬時(shí)序分析
# 該代碼用于進(jìn)行擬時(shí)序分析
#加載分析使用的包
library(Seurat)
library(monocle)
library(ggplot2)
library(cowplot)
library(Matrix)
library(dplyr)
##使用monocle2進(jìn)行擬時(shí)序分析
#構(gòu)造表達(dá)及注釋數(shù)據(jù)氢卡，提取CCA之后的數(shù)據(jù),intergrated是批次校正后的數(shù)據(jù)
exp.matrix<-as(as.matrix(exp.seurat@assays$integrated@data), 'sparseMatrix')
feature_ann<-data.frame(gene_id=rownames(exp.matrix),gene_short_name=rownames(exp.matrix))
rownames(feature_ann)<-rownames(exp.matrix)
exp_fd<-new("AnnotatedDataFrame", data = feature_ann)
sample_ann<-exp.seurat@meta.data
rownames(sample_ann)<-colnames(exp.matrix)
exp_pd<-new("AnnotatedDataFrame", data =sample_ann)
#生成monocle對(duì)象
exp.monocle<-newCellDataSet(exp.matrix,phenoData =exp_pd,featureData =exp_fd,expressionFamily=negbinomial.size())
head(pData(exp.monocle))
head(fData(exp.monocle))
#計(jì)算sizefactor译秦，類似于標(biāo)準(zhǔn)化
exp.monocle <- estimateSizeFactors(exp.monocle)
exp.monocle <- estimateDispersions(exp.monocle)
#根據(jù)seurat cluster計(jì)算差異表達(dá)基因并挑選用于構(gòu)建擬時(shí)序軌跡的基因
diff_test_res<-differentialGeneTest(exp.monocle,fullModelFormulaStr = "~seurat_clusters")
ordering_genes<-row.names (subset(diff_test_res, qval < 0.01))
exp.monocle<-setOrderingFilter(exp.monocle, ordering_genes)
plot_ordering_genes(exp.monocle)
#DDRTree的方法降維并構(gòu)建擬時(shí)序
exp.monocle<-reduceDimension(exp.monocle, max_components = 2, reduction_method = "DDRTree")
exp.monocle<-orderCells(exp.monocle)
colnames(pData(exp.monocle))
#修改monocle對(duì)象中的列名示例
names(pData(exp.monocle))[names(pData(exp.monocle))=="seurat_clusters"]="Cluster"
#將不同分組情況的擬時(shí)序軌跡圖畫到一起
plot1<-plot_cell_trajectory(exp.monocle, color_by = "Cluster",cell_size=1)
plot2<-plot_cell_trajectory(exp.monocle, color_by = "sample",cell_size=1)
plot3<-plot_cell_trajectory(exp.monocle, color_by = "batch",cell_size=1)
plot4<-plot_cell_trajectory(exp.monocle, color_by = "State",cell_size=1)
plot5<-plot_cell_trajectory(exp.monocle, color_by = "Pseudotime",cell_size=1)
pdf(paste0("./",result.name,"/trajectory_plot.pdf"),width = 16,height = 10)
CombinePlots(plots = list(plot1, plot2,plot3,plot4,plot5),legend = NULL)
dev.off()
rm(plot1,plot2,plot3,plot4,plot5)
save.image("./monocle2_analysis_pipeline.RData")
結(jié)果展示
image.png
