单细胞monocle2分析

seurat单细胞数据到monocle2的分析流程 单细胞个性化分析 单细胞拟时序分析

conda activate snapatacv1
conda install -c bioconda bioconductor-monocle
#上面的方法 会报错
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("monocle")
#R4.0可能不行，如果后续很多包报错以至于无法正常运行程序，则必须安装在R3.6版本才能用。

数据读入

思维导图

library(dplyr)
library(Seurat) 
library(patchwork)
library(hdf5r)
library(SeuratDisk)
library(monocle)
packageVersion("monocle")
seurat_obj <- LoadH5Seurat('./data/cs8_meso_new.h5seurat')
`%!in%`=Negate(`%in%`)
table(seurat_obj@meta.data$zonghe_new)
sc_qc=subset(seurat_obj,subset = zonghe_new%!in%c('Endo','NOTO',15,18,5))
sc_meso <- subset(seurat_obj,subset = zonghe_new %in% c('EXE.Meso.Proge', 'Inter.Meso','LP.Meso'))
sc_qc = subset(x = seurat_obj, zonghe_new = c(1, 2))
#importHSMM(sc_qc, import_all = TRUE) #因seurat4数据对象的改变，这个函数已经失效
#手动读取seurat对象
data <- as(as.matrix(sc_qc@assays$RNA@data), 'sparseMatrix')
pd <- new('AnnotatedDataFrame', data = sc_qc@meta.data)
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
fd <- new('AnnotatedDataFrame', data = fData)
#Construct monocle HSMM 这是在构建稀疏矩阵数据集，这种数据集更适合单细胞大量数据的处理
HSMM <- newCellDataSet(data,
                         phenoData = pd,
                         featureData = fd,
                         lowerDetectionLimit = 0.5,
                         expressionFamily = negbinomial.size())
HSMM

质控

主要的质控和分群应该在seurat里就完成，不要再monocle做分群和质控

#估计尺寸因子和分散度
HSMM <- estimateSizeFactors(HSMM)
HSMM <- estimateDispersions(HSMM)
#过滤低质量数据
HSMM <- detectGenes(HSMM, min_expr = 0.1)
print(head(fData(HSMM)))
expressed_genes <- row.names(subset(fData(HSMM),
    num_cells_expressed >= 10))
print(head(pData(HSMM)))
#差异基因筛选
#也可输入seurat筛选出的高变基因：expressed_genes_s <- VariableFeatures(sc_qc) 
HSMM <- setOrderingFilter(HSMM, expressed_genes)
plot_ordering_genes(HSMM)
HSMM <- reduceDimension(HSMM, max_components = 2,
    method = 'DDRTree') #降维，如果不先执行这一步后面可能会报错
diff <- differentialGeneTest(HSMM[expressed_genes],fullModelFormulaStr="~zonghe_new",cores=8) 
#~后面是表示对谁做差异分析的变量，理论上可以为p_data的任意列名
head(diff)

##差异表达基因作为轨迹构建的基因,差异基因的选择标准是qval<0.01,decreasing=F表示按数值增加排序
deg <- subset(diff, qval < 0.01) #基因数目在2000上下比较合适
deg <- deg[order(deg$qval,decreasing=F),]
head(deg)
ordergene <- row.names(deg)[order(deg$qval)][1:800] #如果gene太多了可以选出来头部基因

##差异基因的结果文件保存
write.table(deg,file="train.monocle.DEG.xls",col.names=T,row.names=F,sep="\t",quote=F)

轨迹分析

## 轨迹构建基因可视化
#ordergene <- rownames(deg) 
HSMM <- setOrderingFilter(HSMM, ordergene)  
#这一步是很重要的，在我们得到想要的基因列表后，我们需要使用setOrderingFilter将它嵌入HSMM对象，后续的一系列操作都要依赖于这个list。
#setOrderingFilter之后，这些基因被储存在HSMM@featureData@data[["use_for_ordering"]]，可以通过table(HSMM@featureData@data[["use_for_ordering"]])查看
pdf("train.ordergenes.pdf")
plot_ordering_genes(HSMM)
dev.off()
#出的图黑色的点表示用来构建轨迹的差异基因，灰色表示背景基因。红色的线是根据第2步计算的基因表达大小和离散度分布的趋势(可以看到，找到的基因属于离散度比较高的基因)

#HSMM <- reduceDimension(HSMM, max_components = 2, method = 'DDRTree')
HSMM <- orderCells(HSMM)
#⚠️使用root_state参数可以设置拟时间轴的根，如下面的拟时间着色图中可以看出，左边是根。根据state图可以看出，根是State1，若要想把另一端设为根，可以按如下操作
#HSMM <- orderCells(HSMM, root_state = 5) #把State5设成拟时间轴的起始点
#自己配置时序图的颜色
col_new<-c("#6f87a0","#a4a0af","#1c81b7","#7CCD7C","#76EEC6","#5cbebf","#87CEFA","#cce087")
col_new<-c("#828F2A","#006B69","#4C8744","#3178A8","#7477B4","#AB72AC","#117B5C","#00748C")
###new_label### APS  connect Epi EXE.Meso.Progenitor     Inter.Meso    LP.Meso   NMP     PSM 
#正确方向是EPi-APS-EXE
#library(ggsci)
#HSMM@phenoData@data$State
#HSMM <- orderCells(HSMM, root_state = 5) #, root_state = 3
p1=plot_cell_trajectory(HSMM,color_by="Pseudotime", size=1,show_backbone=TRUE) 
p2=plot_cell_trajectory(HSMM, color_by = "zonghe_new") + facet_wrap("~zonghe_new", nrow = 1) + scale_color_manual(values = col_new)
p3=plot_cell_trajectory(HSMM, color_by = "zonghe_new")  + scale_color_manual(values = col_new)
p4 <- plot_cell_trajectory(HSMM, color_by = "State")
pdf("merge.pdf",width = 15,height = 8)
p1|p3|p4
dev.off()
pdf("main.pdf",width = 10,height = 8)
p3
dev.off()
pdf("explit.pdf",width = 12,height = 6)
p2
dev.off()
p5 <- plot_complex_cell_trajectory(HSMM, x = 1, y = 2,
                                   color_by = "zonghe_new")+
scale_color_manual(values = col_new) +
  theme(legend.title = element_blank()) 
pdf("main_2.pdf",width = 10,height = 8)
p5
dev.off()
df <- pData(HSMM)
p6 <- ggplot(df, aes(Pseudotime, colour = zonghe_new, fill=zonghe_new)) +
  geom_density(bw=0.5,size=1,alpha = 0.5)+ scale_fill_manual(name = "", values = col_new)+scale_color_manual(name = "", values = col_new)
pdf("main_3.pdf",width = 12,height = 8)
p6
dev.off()

一些细节

改名

HSMM@phenoData@data$zonghe_new[HSMM@phenoData@data$zonghe_new == "AD_meso"] <- "EXE.Meso.Progenitor" HSMM@phenoData@data$zonghe_new[HSMM@phenoData@data$zonghe_new == "inter_meso"] <- "Inter.Meso" HSMM@phenoData@data$zonghe_new[HSMM@phenoData@data$zonghe_new == "LP_meso"] <- "LP.Meso" HSMM@phenoData@data$zonghe_new[HSMM@phenoData@data$zonghe_new == "NOTO"] <- "Node"
HSMM@phenoData@data$zonghe_new[HSMM@phenoData@data$zonghe_new == "PS"] <- "Connecting Stalk"

seurat_obj$zonghe_new[which(seurat_obj$zonghe_new == "AD_meso")] <- "EXE.Meso.Progenitor"
seurat_obj$zonghe_new[which(seurat_obj$zonghe_new == "inter_meso")] <- "Inter.Meso"
seurat_obj$zonghe_new[which(seurat_obj$zonghe_new == "LP_meso")] <- "LP.Meso"
seurat_obj$zonghe_new[which(seurat_obj$zonghe_new == "NOTO")] <- "Node"
seurat_obj$zonghe_new[which(seurat_obj$zonghe_new == "PS")] <- "Connecting Stalk"

参考网址

单片眼镜 (cole-trapnell-lab.github.io)

单细胞之轨迹分析-2：monocle2 原理解读+实操 - 简书 (jianshu.com)

安装软件推荐使用github的社区修改版本来本地安装，避免新版本的R4与monocle2兼容性的bug

what is the “condition” in orderCells() function · Issue #434 · cole-trapnell-lab/monocle-release (github.com)

以及…

还想要更多代码？不知道怎么搭建？想让别人帮忙生信分析？欢迎来我的闲鱼咨询！价格绝对全网最低