Monocle2流程
Monocle是用于单细胞拟时序分析的著名R包，通常Monocle2最为常用，而Monocle3仍在开发中。接下来随小编一起来看一看使用Monocle2包进行拟时序分析的简单攻略吧！
01 单细胞数据文件
在考虑进行拟时序分析前，您应该已经完成了单细胞数据的基础分析，包括降维、聚类、细胞类型注释等步骤，相应的Seurat_Object我们在以下称为soj。这里我们使用万乘基因自主研发的10K Genomics-Perseus™平台产生的3` 转录组pbmc数据作为展示。
往期推文：实测数据 | 万乘基因10K平台单细胞核实测数据展示！ (qq.com)；实测案例丨万乘基因10K平台单细胞转录组测序性能测评展示！ (qq.com)。
library(Seurat)soj=readRDS('10k_pbmc.RDS')DimPlot(soj, reduction = "umap", label = T,pt.size = 0.03,label.size = 4,repel = TRUE)
图1. 此处使用万乘基因自主研发的10K Genomics-Perseus™平台产生的3` 转录组pbmc数据作为展示
02 安装相应R包
在Rstudio中安装相应的R包，具体可以查看官网https://cole-trapnell-lab.github.io/monocle-release/docs/，之后我们就可以开始正式的分析了!
03 拟时序轨迹分析
3.1 提取目标细胞群
对于已有的Seurat_Object，我们可以选取感兴趣的几种细胞群，将数据提取出来，用于接下来的分析。这里我们截取T细胞亚群作为展示样例：
library(monocle)#提取目标细胞群soj.downsample=subset(soj,idents=c("NK","Naive_T_cell","CD8+_T"))#此处200指每种细胞群提取200个细胞soj.downsample=subset(soj.downsample,downsample=200)
3.2 构建CellDataSet
首先需要从Seurat_Object中提取我们需要的数据
#提取表达矩阵data <- as(as.matrix(soj.downsample@assays$RNA@counts), 'sparseMatrix')#提取表型信息pd <- new('AnnotatedDataFrame', data = soj.downsample@meta.data)fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))fd <- new('AnnotatedDataFrame', data = fData)
然后构建CDS对象
3.4 降维 使用反向图嵌入DDRTree算法进行数据降维； 然后就可以尝试作图查看结果了！ 图2. 按“拟时间”着色展示的细胞轨迹图；此处的“拟时间”为0的起点为计算程序自动选择，实际使用时可以按实际需求，结合生物学知识，手动选择起点。 图3. 按“细胞类型”着色展示的细胞轨迹图 04 拟时序差异基因分析 4.1 差异基因计算预热图 对先前用作轨迹构建的特征基因进行回归分析，计算它们是否与拟时间有显著关系。 图6. 通过热图展示与拟时间有显著关系的特征基因，可以直观地发现哪些基因有随着拟时间的改变而发生表达量的显著增高或降低。 热图的横坐标为拟时间的顺序，纵坐标一格为一个基因。通过热图可以发现一些基因有随着拟时间增高逐渐表达变高或变低的倾向，这或许可以成为我们未来的研究对象。 参考文献： Trapnell, C., Cacchiarelli, D., Grimsby, J. et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol 32, 381–386 (2014). Qiu, X., Hill, A., Packer, J. et al. Single-cell mRNA quantification and differential analysis with Census. Nat Methods 14, 309–315 (2017). monocle_cds <- newCellDataSet(data, phenoData = pd, featureData = fd, lowerDetectionLimit = 0.5, expressionFamily = negbinomial.size())
monocle_cds <- estimateSizeFactors(monocle_cds)monocle_cds <- estimateDispersions(monocle_cds)
#使用FindVariableFeatures得到高变基因soj.downsample=FindVariableFeatures(soj.downsample,nfeatures = 2000)ordering_genes=VariableFeatures(soj.downsample)#将得到的基因嵌入CDSmonocle_cds <- setOrderingFilter(monocle_cds, ordering_genes)
monocle_cds <- reduceDimension(monocle_cds, max_components = 2, method = 'DDRTree')
monocle_cds <- orderCells(monocle_cds)
png("Pseudotime_trajectory_soj.downsample.png")plot_cell_trajectory(monocle_cds, color_by = "Pseudotime")dev.off()png("Annotation_trajectory_soj.downsample.png")plot_cell_trajectory(monocle_cds, color_by = "annotation",cell_size = 0.5)dev.off()png("State_trajectory_soj.downsample.png")plot_cell_trajectory(monocle_cds, color_by = "State")dev.off()#选取State2作为起点monocle_cds <- orderCells(monocle_cds,root_state = 2)
png("Annotation_trajectory_soj.downsample_faceted.png",width = 1500)plot_cell_trajectory(monocle_cds, color_by = "annotation")+facet_wrap("~annotation",nrow = 1)dev.off()#计算差异基因Time_diff <- differentialGeneTest(monocle_cds[ordering_genes,], cores = 1, fullModelFormulaStr = "~sm.ns(Pseudotime)")#保存为csv文件write.csv(Time_diff, "Time_diff.csv")#使用100个基因作热图。num_clusters=3可选的将热图中基因聚为3个clusterTime_genes <- Time_diff %>% pull(gene_short_name) %>% as.character()Time_genes <- top_n(Time_diff, n = 100, desc(qval)) %>% pull(gene_short_name) %>% as.character()p = plot_pseudotime_heatmap(monocle_cds[Time_genes,], num_clusters=3, show_rownames=T, return_heatmap=T)ggsave("Time_heatmap_Top100.png", p, width = 5, height = 10)4.2 按Cluster提取基因
如果需要提取在上一步的热图中特定Cluster的基因进行分析，可以进行以下操作。
clusters <- cutree(p$tree_row, k = 3)clustering <- data.frame(clusters)clustering[,1] <- as.character(clustering[,1])colnames(clustering) <- "Gene_Clusters"#查看各Cluster基因数table(clustering)#存放在csv中write.csv(clustering, "Time_clustering.csv", row.names = T)
4.3 按热图顺序提取基因
hp.genes <- p$tree_row$labels[p$tree_row$order]Time_heat_sig <- Time_diff[hp.genes, c("gene_short_name", "pval", "qval")]#通过p值筛选p<0.05的显著基因Time_diff_sig <- filter(Time_diff,Time_diff$pval<0.05)write.csv(Time_diff_sig, "Time_diff_sig.csv", row.names = T)
