提供了编程的基础技术教程

zl程序教程

您现在的位置是：首页 >  .Net

当前栏目

马上都2023了，但是CNS级别单细胞文章仍然是使用monocle2

2023-02-18 16:30:35 时间

很多人在交流群问到：自己的单细胞数据分析里面的拟时序环节使用了monocle2，会不会投稿时候被审稿人disss，毕竟monocle2的作者都在强力推荐monocle3这样的更新版本了。

其实大家简单的搜索就能发现 trapnell 实验室虽然出了 monocle3 ，而且写的很清楚：Monocle 2 and Monocle 3 alpha are deprecated， Please use our new package, Monocle 3 ，如下所示的链接 ：

  • http://cole-trapnell-lab.github.io/monocle-release/docs/#installing-monocle
  • https://cole-trapnell-lab.github.io/monocle3/docs/installation/

但是 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算法即可。

而且有大量CNS文章给你背书

完全不惧往往是门外汉的审稿人啊，比如张泽民的最新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对象直接走monocle2流程

# 一般来说，我们不会对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')

相关文章