安装所需要的包

BiocManager::install("monocle")
library(monocle)
library(Seurat)

载入数据集

Neuron <- readRDS("Neuron.rds")
mono_data <- subset(Neuron, idents = c('PEP',"MAAC"))

准备CellDataSets对象所需的数据文件

表达矩阵

matrix <- as(as.matrix(mono_data@assays$RNA@counts), 'sparseMatrix')
matrix <- as(as.matrix(mono_data@assays$RNA@counts), 'dgCMatrix')
matrix <- as.matrix(mono_data@assays$RNA@data)

phenoData

sample<- mono_data@meta.data
sample$cluster <- mono_data@active.ident

featureData

geneAnn<- data.frame(gene_short_name = row.names(matrix), 
                   row.names = row.names(matrix))

构建CDS对象

pd <- new("AnnotatedDataFrame", data = sample)
fd <- new("AnnotatedDataFrame", data = geneAnn)
cds <- newCellDataSet(matrix, phenoData = pd, featureData = fd,
                      lowerDetectionLimit = 0.1,
                      expressionFamily = negbinomial.size()
                      )

expressionFamily参数指定表达矩阵的数据类型，稀疏矩阵为negbinomial.size()或negbinomial()（很慢）；TPM或FPKM为tobit；log转换FPKM/TPMs为gaussianff()。

数据预处理

cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)

estimateSizeFactors：对细胞间mRNA的差异进行归一化，去除由于样本批次带来的影响（去批次效应）；estimateDispersions：帮助以后进行差异表达分析
注：只有在构建CellDataSet时，使用negbinomial()或negbinomial.size()参数，estimateSizeFactors()和estimateDispersions()时才会工作，并且是必需的。

进行轨迹分析

进行差异基因选择

DEG_genes <- differentialGeneTest(cds,fullModelFormulaStr = '~cluster')
ordering_genes <- row.names(expressed_genes)[order(expressed_genes$pval)][1:153]

可视化选择的基因

cds <- setOrderingFilter(cds, ordering_genes)
plot_ordering_genes(cds)

step 2：基于高可变基因的降维分析

cds <- reduceDimension(cds, max_components = 2,
                       method = 'DDRTree')

step 3：沿着时间轨迹，排序细胞

cds <- orderCells(cds)
plot_cell_trajectory(cds)

plot_cell_trajectory

plot_cell_trajectory(cds, color_by = "Pseudotime")+theme_void()

trajectory

参考资料：

CNS单细胞数据的拟时序分析（monocle2算法的实现）和绘图 (qq.com)