拟时分析或者叫轨迹分析，是单细胞转录组分析中常用的分析方法。主要用于推断出发育过程细胞的分化或者基因表达变化。原理一句话概括是对细胞进行排序，模拟出细胞发育过程。
目前Monocle存在三个版本：Monocle2、Monocle3、Monocle-alpha
注：Alpha版: 此版本表示该软件在此阶段主要是以实现软件功能为主，通常只在软件开发者内部交流，一般而言，该版本软件的Bug较多，需要继续修改。Beta版: 该版本相对于α版已有了很大的改进，消除了严重的错误，但还是存在着一些缺陷，需要经过多次测试来进一步消除，此版本主要的修改对像是软件的UI。RC版: 该版本已经相当成熟了，基本上不存在导致错误的BUG，与即将发行的正式版相差无几。Release版: 该版本意味“最终版本”，在前面版本的一系列测试版之后，终归会有一个正式版本，是最终交付用户使用的一个版本。该版本有时也称为标准版。一般情况下，Release不会以单词形式出现在软件封面上，取而代之的是符号(R)。
因为目前Monocle3还没有更新到Release版，同时Seurat包也提出为什么Seurat对象目前不能直接转到Monocle3,主要原因就是Monocle3目前不稳定。（其实是服务器一直安装不上Monocle3…尴尬)
排序原理
排序的基本要素：
- 对什么排序
- 如何判断先后顺序
- 如何寻找分支点（如果分支的话） 拟时分析：（1）关键基因选择 （2）拟时间（排序空间）（3）排序
Monocle2学习
功能包括(将Seurat基础功能基本都包含)：
- 聚类、降维、计数
- 构建轨迹
- 差异表达分析
具体分析流程代码：
# 安装加载monocle2
BiocManager::install("monocle")
library(monocle)
setwd("F:\\github\\newfile\\monocle学习")
# if(FALSE)代表可选步骤
# 构建CellDataSet对象
# 基于ExpressionSet对象而来(之前多用于芯片数据存储)
# 三个输入文件：
# exprs：表达矩阵文件，行为基因名，列为细胞
# phenoData：注释对象，细胞注释，行为细胞，列为细胞属性(如细胞类型，细胞状态等)
# featureData：注释对象，基因注释，行为基因，列为基因属性(如GC含量等)
# phenoData的行数以及行名需要与exprs的列数列名对应
# featureData的行数以及行名需要与exprs的行数行名对应
# featureData其中必须需要一列命名为“gene_short_name"
expr_matrix <- read.table("Root_single_cell_wt_datamatrix.csv",header = T,row.names = 1,sep = ",")
sample_sheet <- read.table("wt_sample.csv",header = T,sep = ",")
rownames(sample_sheet)=colnames(expr_matrix)
gene_annotation <- read.table("wt_gene.txt",header = T)
rownames(gene_annotation)=rownames(expr_matrix)
pd <- new("AnnotatedDataFrame", data = sample_sheet)
fd <- new("AnnotatedDataFrame", data = gene_annotation)
ARA <- newCellDataSet(as.matrix(expr_matrix),
phenoData = pd, featureData = fd,expressionFamily = negbinomial.size()) # expressionFamil参数设置分布
# FPKM/TPM---log-mormally distributed;UMIs---negative binomial
# 预处理
ARA <- estimateSizeFactors(ARA) # 用于评估每个细胞的sizefactor, 用于后续归一化；
ARA <- estimateDispersions(ARA) # 计算基因的方差
head(pData(ARA))
head(fData(ARA))
# 筛选基因(没有固定标准)
disp_table <- dispersionTable(ARA)
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)
ARA <- setOrderingFilter(ARA, unsup_clustering_genes$gene_id)
# 降维分析
ARA <- reduceDimension(ARA, max_components = 2, method = 'DDRTree') # 计算量大，运行巨耗内存，建议服务器运行
# psudotime分析
ARA <- orderCells(ARA)
plot_cell_trajectory(ARA)
# 对于10X数据导入monocle
# 通过Seurat导入
## 切换到工作目录，读入三个文本文件.
pbmc.data <- Read10X(data.dir = ".")
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "10XPBMC")
## seurat2 可以直接读入
sce <- importCDS(pbmc)
## 由seurat3 导入函数newimport 替代上面的importCDS功能.
newimport <- function(otherCDS, import_all = FALSE) {
if(class(otherCDS)[1] == 'Seurat') {
requireNamespace("Seurat")
data <- otherCDS@assays$RNA@counts
if(class(data) == "data.frame") {
data <- as(as.matrix(data), "sparseMatrix")
}
pd <- tryCatch( {
pd <- new("AnnotatedDataFrame", data = otherCDS@meta.data)
pd
},
#warning = function(w) { },
error = function(e) {
pData <- data.frame(cell_id = colnames(data), row.names = colnames(data))
pd <- new("AnnotatedDataFrame", data = pData)
message("This Seurat object doesn't provide any meta data");
pd
})
# remove filtered cells from Seurat
if(length(setdiff(colnames(data), rownames(pd))) > 0) {
data <- data[, rownames(pd)]
}
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
fd <- new("AnnotatedDataFrame", data = fData)
lowerDetectionLimit <- 0
if(all(data == floor(data))) {
expressionFamily <- negbinomial.size()
} else if(any(data < 0)){
expressionFamily <- uninormal()
} else {
expressionFamily <- tobit()
}
valid_data <- data[, row.names(pd)]
monocle_cds <- newCellDataSet(data,
phenoData = pd,
featureData = fd,
lowerDetectionLimit=lowerDetectionLimit,
expressionFamily=expressionFamily)
if(import_all) {
if("Monocle" %in% names(otherCDS@misc)) {
otherCDS@misc$Monocle@auxClusteringData$seurat <- NULL
otherCDS@misc$Monocle@auxClusteringData$scran <- NULL
monocle_cds <- otherCDS@misc$Monocle
mist_list <- otherCDS
} else {
# mist_list <- list(ident = ident)
mist_list <- otherCDS
}
} else {
mist_list <- list()
}
if(1==1) {
var.genes <- setOrderingFilter(monocle_cds, otherCDS@assays$RNA@var.features)
}
monocle_cds@auxClusteringData$seurat <- mist_list
} else if (class(otherCDS)[1] == 'SCESet') {
requireNamespace("scater")
message('Converting the exprs data in log scale back to original scale ...')
data <- 2^otherCDS@assayData$exprs - otherCDS@logExprsOffset
fd <- otherCDS@featureData
pd <- otherCDS@phenoData
experimentData = otherCDS@experimentData
if("is.expr" %in% slotNames(otherCDS))
lowerDetectionLimit <- otherCDS@is.expr
else
lowerDetectionLimit <- 1
if(all(data == floor(data))) {
expressionFamily <- negbinomial.size()
} else if(any(data < 0)){
expressionFamily <- uninormal()
} else {
expressionFamily <- tobit()
}
if(import_all) {
# mist_list <- list(iotherCDS@sc3,
# otherCDS@reducedDimension)
mist_list <- otherCDS
} else {
mist_list <- list()
}
monocle_cds <- newCellDataSet(data,
phenoData = pd,
featureData = fd,
lowerDetectionLimit=lowerDetectionLimit,
expressionFamily=expressionFamily)
# monocle_cds@auxClusteringData$sc3 <- otherCDS@sc3
# monocle_cds@auxOrderingData$scran <- mist_list
monocle_cds@auxOrderingData$scran <- mist_list
} else {
stop('the object type you want to export to is not supported yet')
}
return(monocle_cds)
}
sce <- newimport(pbmc)
