拟时序分析指的是根据不同细胞亚群基因表达量随时间的变化情况，构建细胞谱系发育。这个时间并不是真的时间，而是一个虚拟的时序列，是指细胞与细胞之间的转化和演替的顺序和轨迹。
monocle是我们经常用的拟时序分析工具，可以使用bioconductor安装。目前已经更新到第3版本，monocle3安装有点麻烦，可以参考官网https://cole-trapnell-lab.github.io/monocle3/。
Monocle使用一种算法来学习每个细胞必须经历的基因表达变化序列，作为动态生物学过程的一部分。一旦它了解了基因表达变化的整体“轨迹”，Monocle就可以将每个细胞放置在轨迹中的适当位置。Monocle依靠一种叫做反向图嵌入的机器学习技术来构建单细胞轨迹https://www.nature.com/articles/nmeth.4402
seurat对象载入
单细胞测序一般使用seurat分析，monocle分析seurat对象是需要进行一些处理。
library(monocle)
class(PBMC_CD4Tcell)[[1]]
[1] "Seurat"
data = as(as.matrix(PBMC_CD4Tcell@assays$RNA@data), "sparseMatrix")
pd = new("AnnotatedDataFrame", data = PBMC.combined.CD4Tcell@meta.data)
fData = data.frame(gene_short_name = row.names(data), row.names = row.names(data))
fd = new("AnnotatedDataFrame", data = fData)
clustered_pbmc_monocle = newCellDataSet(data, phenoData = pd, featureData = fd,
lowerDetectionLimit = 0, expressionFamily = negbinomial.size())
拟时序分析
monocle2拟时序分析主要基于以下三个步骤
1、基因筛选：Monocle可选取感兴趣的基因（差异基因），并利用这些基因来构造数据。
2、降低维度：选择用于细胞排序的基因，对数据进行降维处理。
3、pseudotime对细胞排序：通过将表达数据投影到较低维空间，构建细胞间的分化轨迹。
#计算SizeFactor和Dispersions
clustered_pbmc_monocle <- estimateSizeFactors(clustered_pbmc_monocle)
clustered_pbmc_monocle <- estimateDispersions(clustered_pbmc_monocle, parallel = T)
#选择定义细胞发展的基因
disp_table <- dispersionTable(clustered_pbmc_monocle)
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)
#数据降维
Mono.cds <- setOrderingFilter(clustered_pbmc_monocle, unsup_clustering_genes$gene_id)
Mono.cds <- reduceDimension(
Mono.cds,
max_components = 2,
method = 'DDRTree')
#计算psudotime值
Mono.cds <- orderCells(Mono.cds)
plot_cell_trajectory(Mono.cds,cell_size = 1,color_by = "State")
结果如下图，各个点代表一个细胞，具有相似状态的细胞被聚在一起。每个分支点代表一个可能的生物学过程的决策点。
ordering_genes = unsup_clustering_genes$gene_id
plot_pseudotime = plot_pseudotime_heatmap(Mono.cds.all[ordering_genes, ], num_clusters = 3,
cores = 10, return_heatmap = T, show_rownames = F)
细胞间状态转换相关基因表达热图如下图,热图显示大体分为两类基因，一类从低到高表达，一类从高到底表达，还有一类中间高表达（也可多分几类）。类似基因共表达模块，具有相同表达模式的基因聚在一起。如下图第3个cluster基因，在分化的起始高表达，随着拟时间的增加，表达量逐渐降低。那么这些模块的基因可能就存在这某个生物学的功能或者通路，接下来就可以进行进一步的挖掘。
按聚类数3划分基因聚类关系，可以提取这三部分基因进行go、kegg注释。
clusters <- cutree(plot_pseudotime$tree_row, k = 3)
clustering <- data.frame(clusters)
clustering[,1] <- as.character(clustering[,1])
colnames(clustering) <- "Gene_Clusters"
#接下来就可以用rownames(clustering)提取基因，进行GO KEGG分析了
head(clustering)
# Gene_Clusters
#AAK1 1
#ABCA7 1
#ABHD14B 2
#ABI2 2
#ABLIM1 2
#AC004057.1 2
其他方法
TSCAN
文献 10.1093/nar/gkw430
http://www.bioconductor.org/packages/release/bioc/html/TSCAN.html
