偷偷问一下，关注了吗？
内容获取
1、购买打包合集（《KS科研分享与服务》付费内容打包集合），价格感人，可以加入微信VIP群（答疑交流群，甚至有小伙伴觉得群比代码更好），可以获取建号以来所有内容，群成员专享视频教程，提前更新，其他更多福利！
2、《KS科研分享与服务》公众号有QQ群，进入门槛是20元（完全是为了防止白嫖党，请理解），请考虑清楚。群里有免费推文的注释代码和示例数据（终身拥有），没有付费内容，群成员福利是购买单个付费内容半价！
需要者详情请联系作者（非需要者勿扰，处理太费时间）：
Monocle2个性化作图
library(Seurat)library(monocle)library(viridis)#========================================================================================setwd("D:\\KS项目\\公众号文章\\视频教程---monocle2分析测试")#load data-加载之前做了cytroTRACE2的那个obj#可以结合辅助看看cell起点load("D:/KS项目/公众号文章/CytroTRACE2：单细胞转录组数据推断细胞发育潜能/sce_trace.RData")DimPlot(cytotrace2_sce, label = T, pt.size = 1)#这个数据是大群里面直接提取了目标细胞进行分析#========================================================================================#monocle pieline function#pay attention---Seuratobject version 5.0.0ks_run_Monocle2 <- function(object, #seurat obj or expression matrix (建议数据格式转为matrix,如果数据量大转化为稀疏矩阵as(as.matrix(data), "sparseMatrix"))layer, #used when object is a seurat objassay, #used when object is a seurat objlowerDetectionLimit = 0.1,VARgenesM=c("dispersionTable","seurat","differentialGeneTest"),cellAnno=NULL,define_root=F,root_state,reverse=NULL){if(class(object)[1] == 'Seurat') {data <- GetAssayData(object=object, layer=layer, assay=assay)#get expression matrix datapd <- new("AnnotatedDataFrame", data = object@meta.data)fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))fd <- new("AnnotatedDataFrame", data = fData)if(all(data == floor(data))) {expressionFamily <- negbinomial.size()} else if(any(data < 0)){expressionFamily <- uninormal()} else {expressionFamily <- tobit()}#Creates a new CellDateSet object.monocle_cds <- newCellDataSet(data,phenoData = pd,featureData = fd,lowerDetectionLimit=0.1,expressionFamily=expressionFamily)}else{print("This fucntions only apply for a seurat obj")}#Estimate size factors and dispersions#数据处理monocle_cds <- estimateSizeFactors(monocle_cds)#size facotr标准化细胞之间的mRNA差异monocle_cds <- estimateDispersions(monocle_cds)#质量控制-filter cellsmonocle_cds <- detectGenes(monocle_cds, min_expr = 0.1)# print(head(fData(monocle_cds)))# print(head(pData(monocle_cds)))# expressed_genes <- row.names(subset(fData(mouse_monocle), num_cells_expressed >= 10))monocle_cds <- monocle_cds[fData(monocle_cds)$num_cells_expressed >= 10, ]#select methods for VariableFeaturesif(VARgenesM=="dispersionTable"){disp_table <- dispersionTable(monocle_cds)ordering_genes <- subset(disp_table,mean_expression >= 0.1 &dispersion_empirical >= 1.5* dispersion_fit)$gene_id}if(VARgenesM=="seurat"){ordering_genes <- VariableFeatures(FindVariableFeatures(object, assay = "RNA"), assay = "RNA")}if(VARgenesM=="differentialGeneTest"){diff_test_res <- differentialGeneTest(monocle_cds,fullModelFormulaStr = paste0("~",cellAnno))##~后面是表示对谁做差异分析的变量diff_test_res_sig <- diff_test_res[order(diff_test_res$qval,decreasing=F),]ordering_sce <- diff_test_res_sig[diff_test_res_sig$qval< 0.01,]if(nrow(ordering_sce)>3000){ordering_genes <- ordering_sce$gene_short_name[1:3000]}else{ordering_genes <- rdering_sce$gene_short_name}}#Marks genes for clusteringmonocle_cds <- setOrderingFilter(monocle_cds, ordering_genes)plot_ordering_genes(monocle_cds)#clustermonocle_cds <- reduceDimension(monocle_cds, max_components = 2,reduction_method = 'DDRTree')#order cellsmonocle_cds <- orderCells(monocle_cds, reverse=reverse)if(define_root){monocle_cds <- monocle_cds <- orderCells(monocle_cds,root_state = root_state)}return(monocle_cds)}
sce_CDS <- ks_run_Monocle2(object = cytotrace2_sce,layer = 'counts',assay = "RNA",VARgenesM="dispersionTable",cellAnno = "cluster")
#可视化拟时plot_cell_trajectory(sce_CDS,show_branch_points = F,color_by = "Pseudotime",cell_size = 3)+scale_color_gradientn(colours = colorRampPalette(c('blue','green','yellow','red'))(100))
#按照celltype着色plot_cell_trajectory(sce_CDS,show_branch_points = F,color_by = "cluster",cell_size = 3)+theme_dr(arrow = grid::arrow(length = unit(0, "inches")))+theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+scale_color_manual(values = c("#E69253", "#EDB931", "#E4502E", "#4378A0"))+geom_curve(aes(x=9,y=2,yend=0,xend=5),arrow=arrow(length=unit(0.03,"npc")),size=1,curvature=0,color="black")+annotate("text", x=6.5, y=1.5, label = "bold(YSMP)", parse = TRUE, angle=45, size=4)+geom_curve(aes(x=-4,y=-1,yend=3,xend=-6),arrow=arrow(length=unit(0.03,"npc")),size=1,curvature=0,color="black")+annotate("text", x=-4, y=1, label = "bold(Monocyte~Fate)", parse = TRUE, angle=-75, size=4)+geom_curve(aes(x=-5,y=-3,yend=-5,xend=-6),arrow=arrow(length=unit(0.03,"npc")),size=1,curvature=0,color="black")+annotate("text", x=-6.5, y=-4, label = "bold(Neu~Fate)", parse = TRUE, angle=75, size=4)
#结合可视化cytroTRACE2结果plot_cell_trajectory(sce_CDS,show_branch_points = F,color_by = "CytoTRACE2_Relative",cell_size = 3)+scale_colour_gradientn(colours = (c( "#000004FF", "#3B0F70FF", "#8C2981FF", "#DE4968FF", "#FE9F6DFF", "#FCFDBFFF")),na.value = "transparent",limits=c(0,1),breaks = c(0,1),labels=c("(More diff.)", "(Less diff.)"),name = "Relative\norder \n" ,guide = guide_colorbar(frame.colour = "black", ticks.colour = "black"))+theme_dr(arrow = grid::arrow(length = unit(0, "inches")))+#坐标轴主题修改theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())
#=======================================================================================#拟时差异基因expressed_genes=row.names(subset(fData(sce_CDS),use_for_ordering=="TRUE"))peu_gene <- differentialGeneTest(sce_CDS[expressed_genes,],fullModelFormulaStr = "~sm.ns(Pseudotime)",cores = 1)# write.csv(peu_gene,file='peu_gene.csv')#保存好文件，这个分析过程挺费时间peu_gene <- peu_gene[which(peu_gene$qval<0.01 & peu_gene$num_cells_expressed>100),]#筛选显著是的# peu_gene %>% arrange(qval) -> peu_gene#按照qval排个序# peu_gene <- peu_gene[1:100,] #这里我们取前100个基因演示#改造热图函数，展示需要的基因source('./new_heatmap.R')library(pheatmap)library(grid)p1 <-plot_pseudotime_heatmap(sce_CDS[peu_gene$gene_short_name,],num_clusters = 4,cores = 2,show_rownames = T,return_heatmap =T,hmcols = viridis(256),use_gene_short_name = T,clustercolor = c("#d2981a", "#a53e1f", "#457277", "#8f657d"))p1#只展示感兴趣的基因genes <- c("C1QB","C1QC","C1QA","MRC1","LGMN","MS4A7","MAF","FOLR2","HLA-DPA1","CLEC10A","IL10RA","CD163","KCTD12","CLEC7A","MS4A6A","CD14","ITM2A","CYTL1","MDK","SELP","CD24","S100A8","S100A9","S100A12")source('./add.flag.R')add.flag(p1, kept.labels = genes, repel.degree = 0.2)
#提取module基因，进行富集分析，结合热图展示library(clusterProfiler)library(ggplot2)source('./Monocle2_gene_enrichment.R')GOannlysis <- Monocle2_gene_enrichment(p,knum=4,species='org.Hs.eg.db',pvalueCutoff=0.05,qvalueCutoff=0.05)table(GOannlysis$Cluster)#module1基因最少，没有显著性富集# 1 2 3 4# 0 76 43 77write.csv(GOannlysis, file = 'GOannlysis.csv')
library(RColorBrewer)= T,color_by = "State",cell_size = 3)#BEAM#查看在分叉处的差异基因BEAM_res <- BEAM(sce_CDS, branch_point = 3, cores=2)BEAM_res <- BEAM_res[order(BEAM_res$qval),]BEAM_res <- BEAM_res[,c("gene_short_name", "pval", "qval")]# write.csv(BEAM_res,file = "BEAM_res.csv")branch_p <- plot_genes_branched_heatmap(sce_CDS[row.names(subset(BEAM_res, qval < 1e-8)),],branch_point = 3,num_clusters = 4,cores = 2,hmcols = colorRampPalette(c("steelblue", "white","darkorange2", "orangered4"))(100),use_gene_short_name = T,show_rownames = T,return_heatmap = T)
觉得我们分享有些用的，点个赞再走呗！
关注我们获取精彩内容：
关注不迷路：扫描下面二维码关注公众号！
B站视频号链接：https://space.bilibili.com/471040659?spm_id_from=333.1007.0.0
关注 KS科研分享与服务，
认清正版优质内容和服务！
优质内容持续输出，物超所值！
合作联系：ks_account@163.com
