2022-01-20 拟时序分析-monocle2包

安装

#安装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")

创建文件夹存储结果

#创建一个文件夹用于写分析结果
result.name <- "cca_monocle_result"
if(!dir.exists(result.name)){
  dir.create(result.name)
}

用CCA后的数据进行拟时序分析

# 该代码用于进行拟时序分析

#加载分析使用的包
library(Seurat)
library(monocle)
library(ggplot2)
library(cowplot)
library(Matrix)
library(dplyr)

##使用monocle2进行拟时序分析
#构造表达及注释数据，提取CCA之后的数据,intergrated是批次校正后的数据
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对象
exp.monocle<-newCellDataSet(exp.matrix,phenoData =exp_pd,featureData =exp_fd,expressionFamily=negbinomial.size())
head(pData(exp.monocle))
head(fData(exp.monocle))

#计算sizefactor，类似于标准化
exp.monocle <- estimateSizeFactors(exp.monocle)
exp.monocle <- estimateDispersions(exp.monocle)

#根据seurat cluster计算差异表达基因并挑选用于构建拟时序轨迹的基因
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的方法降维并构建拟时序
exp.monocle<-reduceDimension(exp.monocle, max_components = 2, reduction_method = "DDRTree")
exp.monocle<-orderCells(exp.monocle)
colnames(pData(exp.monocle))

#修改monocle对象中的列名示例
names(pData(exp.monocle))[names(pData(exp.monocle))=="seurat_clusters"]="Cluster"

#将不同分组情况的拟时序轨迹图画到一起
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")

结果展示