很多人在交流群问到:自己的单细胞数据分析里面的拟时序环节使用了monocle2,会不会投稿时候被审稿人disss,毕竟monocle2的作者都在强力推荐monocle3这样的更新版本了。
其实大家简单的搜索就能发现 trapnell 实验室虽然出了 monocle3 ,而且写的很清楚:Monocle 2 and Monocle 3 alpha are deprecated, Please use our new package, Monocle 3 ,如下所示的链接 :
但是 monocle3 目前仍然是在 github上面,并不是在cran或者Bioconductor,它的安装方式是:
# install monocle3 through the cole-trapnell-lab GitHub, execute:
install.packages("devtools")
devtools::install_github('cole-trapnell-lab/monocle3')
# If you wish to install the develop branch of monocle3, execute:
devtools::install_github('cole-trapnell-lab/monocle3', ref="develop")
对中国大陆地区的小伙伴来说,并不是很友好。而且它 依赖于一个过期的包,https://cran.r-project.org/src/contrib/Archive/Matrix.utils/ ,这样的话很多人(代码能力不够)压根就安装不上啊。
反倒是 Monocle2就很友好啦,安装超级简单, 在bioconductor就有啊 :https://bioconductor.org/packages/release/bioc/html/monocle.html
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("monocle")
目前我看到的这个monocle是 2.26.0 版本 ,如果大家的R版本跟这个冲突,似乎也是让初学者头疼的事情。
不过,使用起来也是标准代码即可,而且结果的解读都是我们反复分享过的。这个拟时序分析流程,我们早就在前面的教程:拟时序分析就是差异分析的细节剖析,我们展现了一个表达量矩阵如何去走Monocle2分析,通常我们的表达量矩阵在seurat对象里面, 首先导出,然后构建Monocle2对象,过滤细胞,选择基因,然后降维的时候选择默认DDRTree算法即可。
完全不惧往往是门外汉的审稿人啊,比如张泽民的最新nature单细胞文章:《Liver tumour immune microenvironment subtypes and neutrophil heterogeneity》,就是3个算法一起使用:
3个算法一起使用
最后的结果如下所示:
e, Monocle trajectories of neutrophils coloured by tissues (left), cluster identities (middle) and CytoTRACE scores (right).
可视化
Each dot represents a single cell. Cell orders are inferred based on the expression of the most variable genes across neutrophil clusters.
# 一般来说,我们不会对seurat对象里面的全部的单细胞亚群进行拟时序
#
# seurat=subset(seurat,idents = c(1,3))
DimPlot(seurat)
#选择需要构建细胞分化轨迹的细胞类型(subset提取感兴趣的名字,上面已经准备好名字了)
expr_matrix=seurat@assays$RNA@counts#使用counts表达值
sample_sheet<-seurat@meta.data#将实验信息赋值新变量
gene_annotation=data.frame(gene_short_name=rownames(seurat))#构建一个含有基因名字的数据框
rownames(gene_annotation)=rownames(seurat)#将上述数据框的行名赋值基因名字
pd <- new("AnnotatedDataFrame", data = sample_sheet)#将实验信息变量转化为monocel可以接收的对象
fd <- new("AnnotatedDataFrame", data = gene_annotation)#将基因注释变量转化为monocle可以接收的对象
cds <- newCellDataSet(expr_matrix, phenoData = pd, featureData = fd,
expressionFamily=negbinomial.size())#创建一个monocle的对象
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
# 并不是所有的基因都有作用,所以先进行挑选,合适的基因用来进行聚类。
disp_table <- dispersionTable(cds)
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)
cds <- setOrderingFilter(cds, unsup_clustering_genes$gene_id)
plot_ordering_genes(cds)
plot_pc_variance_explained(cds, return_all = F) # norm_method='log'
# 其中 num_dim 参数选择基于上面的PCA图
cds <- reduceDimension(cds, max_components = 2, num_dim = 6,
reduction_method = 'tSNE', verbose = T)
cds <- clusterCells(cds, num_clusters = 6)
plot_cell_clusters(cds, 1, 2 )
table(pData(cds)$Cluster)
colnames(pData(cds))
# seurat$RNA_snn_res.0.2
table(pData(cds)$Cluster,pData(cds)$RNA_snn_res.0.2)
plot_cell_clusters(cds, 1, 2 )
save(cds,file = 'input_cds.Rdata')
# 接下来很重要,到底是看哪个性状的轨迹
colnames(pData(cds))
table(pData(cds)$Cluster)
table(pData(cds)$Cluster,pData(cds)$RNA_snn_res.0.2)
plot_cell_clusters(cds, 1, 2 )
pData(cds)$Cluster = pData(cds)$RNA_snn_res.0.2
Sys.time()
diff_test_res <- differentialGeneTest(cds,
fullModelFormulaStr = "~Cluster")
Sys.time()
# Select genes that are significant at an FDR < 10%
sig_genes <- subset(diff_test_res, qval < 0.1)
sig_genes=sig_genes[order(sig_genes$pval),]
head(sig_genes[,c("gene_short_name", "pval", "qval")] )
cg=as.character(head(sig_genes$gene_short_name))
# 前面是找差异基因,后面是做拟时序分析
# 第一步: 挑选合适的基因. 有多个方法,例如提供已知的基因集,这里选取统计学显著的差异基因列表
ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
ordering_genes
cds <- setOrderingFilter(cds, ordering_genes)
plot_ordering_genes(cds)
cds <- reduceDimension(cds, max_components = 2,
method = 'DDRTree')
# 第二步: 降维。降维的目的是为了更好的展示数据。函数里提供了很多种方法,
# 不同方法的最后展示的图都不太一样, 其中“DDRTree”是Monocle2使用的默认方法
cds <- reduceDimension(cds, max_components = 2,
method = 'DDRTree')
# 第三步: 对细胞进行排序
cds <- orderCells(cds)
# 最后两个可视化函数,对结果进行可视化
plot_cell_trajectory(cds, color_by = "Cluster")
ggsave('monocle_cell_trajectory_for_Cluster.pdf')
length(cg)
plot_genes_in_pseudotime(cds[cg,],
color_by = "Cluster")
ggsave('monocle_plot_genes_in_pseudotime_for_Cluster.pdf')
# https://davetang.org/muse/2017/10/01/getting-started-monocle/
my_cds_subset=cds
# pseudotime is now a column in the phenotypic data as well as the cell state
head(pData(my_cds_subset))
# 这个differentialGeneTest会比较耗费时间
my_pseudotime_de <- differentialGeneTest(my_cds_subset,
fullModelFormulaStr = "~sm.ns(Pseudotime)",
cores = 1)
# 不知道为什么无法开启并行计算了
head(my_pseudotime_de)
save( my_cds_subset,my_pseudotime_de,file = 'output_of_monocle.Rdata')