加载需要的包
library(monocle)
library(ggplot2)
library(Seurat)
加载函数
order.gene参数可选项：
- variable
- difftest
CreatMonocleTrajectory <- function(seurat.object,order.gene,out.path){
exp_matrix <- GetAssayData(object = seurat.object[["integrated"]],slot = "data")
exp_matrix <- 2^exp_matrix
pheno_Data <- seurat.object@meta.data
feature_Data <- data.frame(gene_short_name = row.names(exp_matrix),row.names = row.names(exp_matrix))
pd <- new("AnnotatedDataFrame",data = pheno_Data)
fd <- new("AnnotatedDataFrame",data = feature_Data)
cds <- newCellDataSet(as(as.matrix(exp_matrix), "sparseMatrix"), phenoData = pd, featureData = fd,expressionFamily = tobit())
#expressionFamily参数，稀疏矩阵用negbinomial.size()，FPKM值用tobit()，logFPKM值用uninormal()
cds <- estimateSizeFactors(cds)
print("estimateSizeFactors OK")
#cds <- estimateDispersions(cds)
#print("estimateDispersions OK")
#cds <- detectGenes(cds,min_expr=0.1) seurat已经过滤过无需再过滤
print("cds have been constructed")
if (order.gene=="variable"){
variable_genes <- VariableFeatures(seurat.object)
cds <-setOrderingFilter(cds,variable_genes)
print("setOrder OK")
}
if (order.gene=="difftest"){
clustering_DEG_genes <-differentialGeneTest(cds,fullModelFormulaStr = '~seurat_clusters', cores = 4)
ordering_genes <-row.names(subset(clustering_DEG_genes,qval<0.01))
cds <-setOrderingFilter(cds,ordering_genes)
print("setOrder OK")
}
#plot_4 <- plot_ordering_genes(cds)
#ggsave('ordering_variable_genes.pdf',plot_4,path=out.path)
print("start reduceDimension")
cds <- reduceDimension(cds,method = 'DDRTree')
print("start orderCells")
cds <- orderCells(cds)
plot_5 <- plot_cell_trajectory(cds, color_by='seurat_clusters')
ggsave('monocle_trajectory.pdf',plot_5,path=out.path)
plot_6 <- plot_cell_trajectory(cds, color_by='orig.ident')
ggsave('monocle_trajectory_sample.pdf',plot_6,path=out.path)
plot_7 <- plot_cell_trajectory(cds, color_by='seurat_clusters') + facet_wrap(~seurat_clusters, nrow=3)
ggsave('monocle_trajectory_cluster_sep.pdf',plot_7,path=out.path)
plot_8 <- plot_cell_trajectory(cds, color_by='orig.ident') + facet_wrap(~orig.ident, nrow=3)
ggsave('monocle_trajectory_sample_sep.pdf',plot_8,path=out.path)
for (i in highlight_genes){
key_gene_plot <- plot_cell_trajectory(cds, markers=i,use_color_gradient=TRUE,,cell_size=0.5)+scale_colour_gradient(low ="white",high ="#3333FF")
ggsave(paste0(i,"_gene_plot.pdf"),key_gene_plot,path=out.path)
}
plot_9 <- plot_cell_trajectory(cds, color_by='State')
ggsave('monocle_trajectory_by_state.pdf',plot_9,path=out.path)
saveRDS(cds,file=paste(out.path,'monocle_cds.rds',sep='/'))
return(cds)
}
加载数据
path_out_genomics_debach_FAP <- "/data/zhaoxueya/project/gut/result/10Xgenomics/debach/subset/FAP"
path_out_genomics_debach_FAP_monocle <- "/data/zhaoxueya/project/gut/result/10Xgenomics/debach/subset/FAP/monocle"
seurat_object_genomics_FAP <- readRDS(file=paste(path_out_genomics_debach_FAP,"seurat_object_subset.rds",sep="/"))
CreatMonocleTrajectory(seurat.object=seurat_object_genomics_FAP,order.gene="variable",out.path=path_out_genomics_debach_FAP_monocle)
