主要通过MonocleR包，使用反向图形嵌入(Reversed Graph Embedding)的机器学习算法，来学习描述细胞如何从一种状态过渡到另外一种状态的轨迹。其中分析的前提需要一张展现细胞转录特征相似性关系的图，Monocle2使用DDTree降维图，Monocle3使用UMAP降维图。 轨迹分析的前提是待分析的细胞有紧密的发育关系，PBMC细胞不是很好的的示例数据，我们选择T细胞群体演示一下。Monocle建议导入原始表达矩阵，由它完成数据标准化和其他预处理。

数据的导入与处理

rm(list = ls())
library(monocle)
library(Seurat)
library(tidyverse)
scRNAsub <- readRDS("scRNAsub.rds")  #scRNAsub是上一节保存的T细胞子集seurat对象
#生成mycds对象
data <- as(as.matrix(scRNAsub@assays$RNA@counts), 'sparseMatrix')
#count矩阵
pd <- new('AnnotatedDataFrame', data = scRNAsub@meta.data)
# meta表转成特定格式
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
fd <- new('AnnotatedDataFrame', data = fData)
# 基因名表转成特定格式
mycds <- newCellDataSet(data,
                        phenoData = pd,
                        featureData = fd,
                        expressionFamily = negbinomial.size())
#expressionFamily参数用于指定表达矩阵的数据类型，有几个选项可以选择：
#稀疏矩阵用negbinomial.size()，FPKM值用tobit()，logFPKM值用gaussianff()
#mycds是Monocle为我们的数据生成的对象，相当于我们在seurat使用的scRNA对象。数据导入后需要进行标准化和其他预处理：
mycds <- estimateSizeFactors(mycds)
mycds <- estimateDispersions(mycds, cores=1, relative_expr = TRUE)

与seurat把标准化后的表达矩阵保存在对象中不同，monocle只保存一些中间结果在对象中，需要用时再用这些中间结果转化。经过上面三个函数的计算，mycds对象中多了SizeFactors、Dipersions、num_cells_expressed和num_genes_expressed等信息。

选择代表基因

完成数据导入和预处理后，就可以考虑选择哪些基因代表细胞的发育特征，Monocle官网教程提供了4个选择方法：

• 选择发育差异表达基因
• 选择clusters差异表达基因
• 选择离散程度高的基因
• 自定义发育marker基因

前三种都是无监督分析方法，细胞发育轨迹生成完全不受人工干预；最后一种是半监督分析方法，可以使用先验知识辅助分析。第一种方法要求实验设计有不同的时间点，对起点和终点的样本做基因表达差异分析，挑选显著差异的基因进行后续分析。对于没有时序设计的实验样本，可以使用第2、3种方法挑选基因。第2种方法要先对细胞降维聚类，然后用clusters之间差异表达的基因开展后续分析。Monocle有一套自己的降维聚类方法，与seurat的方法大同小异，很多教程直接使用seurat的差异分析结果。第3种方法使用离散程度高的基因开展分析，seurat有挑选高变基因的方法，monocle也有自己选择的算法。本案例数据不具备使用第1、4种方法的条件，因此这里只演示2、3种方法的使用。

##使用clusters差异表达基因
diff.wilcox = FindAllMarkers(scRNAsub)
all.markers = diff.wilcox %>% select(gene, everything()) %>% subset(p_val<0.05)
diff.genes <- subset(all.markers,p_val_adj<0.01)$gene
mycds <- setOrderingFilter(mycds, diff.genes)
p1 <- plot_ordering_genes(mycds)

##使用seurat选择的高变基因
var.genes <- VariableFeatures(scRNAsub)
mycds <- setOrderingFilter(mycds, var.genes)
p2 <- plot_ordering_genes(mycds)
##使用monocle选择的高变基因
disp_table <- dispersionTable(mycds)
disp.genes <- subset(disp_table, mean_expression >= 0.1 & dispersion_empirical >= 1 * dispersion_fit)$gene_id
mycds <- setOrderingFilter(mycds, disp.genes)
p3 <- plot_ordering_genes(mycds)
##结果对比
p1|p2|p3
#选择不同的基因集，拟时分析的结果不同，实践中可以几种方法都试一下。

降维及细胞排序

使用disp.genes开展后续分析

#降维及细胞排序
使用disp.genes开展后续分析
#降维
mycds <- reduceDimension(mycds, max_components = 2, method = 'DDRTree')
#排序
mycds <- orderCells(mycds)
#State轨迹分布图
p1 <- plot1 <- plot_cell_trajectory(mycds, color_by = "State")
#分面展示
p2 <- plot_cell_trajectory(mycds, color_by = "State") + facet_wrap(~State, nrow = 1)
p1|p2

###Cluster轨迹分布图
p1 <- plot_cell_trajectory(mycds, color_by = "seurat_clusters")
p2 <- plot_cell_trajectory(mycds, color_by = "seurat_clusters") + facet_wrap(~seurat_clusters, nrow = 1)
p1/p2

#Pseudotime轨迹图
plot3 <- plot_cell_trajectory(mycds, color_by = "Pseudotime")
plot3

#Monocle基因可视化
s.genes <- c("ITGB1","CCR7","KLRB1","GNLY")
# 点图（抖动）
p1 <- plot_genes_jitter(mycds[s.genes,], grouping = "State", color_by = "State")
# 小提琴图
p2 <- plot_genes_violin(mycds[s.genes,], grouping = "State", color_by = "State")
# 伪时间图
p3 <- plot_genes_in_pseudotime(mycds[s.genes,], color_by = "State")
plotc <- p1|p2|p3
plotc

拟时相关基因聚类热图

Monocle中differentialGeneTest()函数可以按条件进行差异分析，将相关参数设为fullModelFormulaStr = "~sm.ns(Pseudotime)"时，可以找到与拟时先关的差异基因。我们可以按一定的条件筛选基因后进行差异分析，全部基因都输入会耗费比较长的时间。建议使用cluster差异基因或高变基因输入函数计算。分析结果主要依据qval区分差异的显著性，筛选之后可以用plot_pseudotime_heatmap函数绘制成热图。

拟时相关基因聚类热图

##cluster差异基因
diff.wilcox = FindAllMarkers(scRNAsub)
all.markers = diff.wilcox %>% select(gene, everything()) %>% subset(p_val<0.05)
diff.genes <- subset(all.markers,p_val_adj<0.01)$gene
sig_diff.genes <- subset(diff.genes,p_val_adj<0.0001&abs(avg_logFC)>0.75)$gene
sig_diff.genes <- unique(as.character(sig_diff.genes))
diff_test <- differentialGeneTest(mycds[sig_diff.genes,], cores = 1, 
                                  fullModelFormulaStr = "~sm.ns(Pseudotime)")
sig_gene_names <- row.names(subset(diff_test, qval < 0.01))
p1 = plot_pseudotime_heatmap(mycds[sig_gene_names,], num_clusters=3,
                             show_rownames=T, return_heatmap=T)
#这一步报错了。。暂时不知道怎么解决。。。

解决方法·
diff.genes <- read.csv('subcluster/diff_genes_wilcox.csv')
sig_diff.genes <- subset(diff.genes,p_val_adj<0.0001&abs(avg_log2FC)>0.75)$gene
sig_diff.genes <- unique(as.character(sig_diff.genes))
diff_test <- differentialGeneTest(mycds[sig_diff.genes,], cores = 1, 
                                  fullModelFormulaStr = "~sm.ns(Pseudotime)")
sig_gene_names <- row.names(subset(diff_test, qval < 0.01))
p1 = plot_pseudotime_heatmap(mycds[sig_gene_names,], num_clusters=3,
                             show_rownames=T, return_heatmap=T)


#高变基因
disp_table <- dispersionTable(mycds)
disp.genes <- subset(disp_table, mean_expression >= 0.5&dispersion_empirical >= 1*dispersion_fit)
disp.genes <- as.character(disp.genes$gene_id)
diff_test <- differentialGeneTest(mycds[disp.genes,], cores = 4, 
                                  fullModelFormulaStr = "~sm.ns(Pseudotime)")
sig_gene_names <- row.names(subset(diff_test, qval < 1e-04))#以qval为指标，挑选显著
p2 = plot_pseudotime_heatmap(mycds[sig_gene_names,], num_clusters=5,
                             show_rownames=T, return_heatmap=T)

BEAM分析

单细胞轨迹中通常包括分支，它们的出现是因为细胞的表达模式不同。当细胞做出命运选择时，或者遗传、化学或环境扰动时，就会表现出不同的基因表达模式。BEAM(Branched expression analysis modeling)是一种统计方法，用于寻找以依赖于分支的方式调控的基因。

#BEAM分析
disp_table <- dispersionTable(mycds)
disp.genes <- subset(disp_table, mean_expression >= 0.5&dispersion_empirical >= 1*dispersion_fit)
disp.genes <- as.character(disp.genes$gene_id)
mycds_sub <- mycds[disp.genes,]
plot_cell_trajectory(mycds_sub, color_by = "State")
beam_res <- BEAM(mycds_sub, branch_point = 1, cores = 1)
beam_res <- beam_res[order(beam_res$qval),]
beam_res <- beam_res[,c("gene_short_name", "pval", "qval")]
mycds_sub_beam <- mycds_sub[row.names(subset(beam_res, qval < 1e-4)),]
plot_genes_branched_heatmap(mycds_sub_beam,  branch_point = 1, num_clusters = 3, show_rownames = T)
saveRDS(mycds, file="mycds.rds")