轨迹分析系列：

• 单细胞之轨迹分析-1：RNA velocity

拟时间序列分析（Pseudotime分析）的字面意思是通过构建细胞间的变化轨迹来重塑细胞随着时间的变化过程。从具体的分类分析和复杂程度来说，可以分为细胞轨迹分析和细胞谱系分析。
其原理此前已经介绍过，参考单细胞测序的轨迹推断。
细胞轨迹分析指的是简单模型的细胞变化轨迹分析，通常指的是细胞沿着某个过程有特定化的变化终点，轨迹具有简单树状结构，一端是“根”，另一端是“叶”；细胞谱系分析通常指的是某类祖源细胞，在特定条件下，有多个发育轨迹和命运，变化过程类似复杂树状分支变化过程。因此，简单细胞轨迹分析和细胞谱系分析原理上类似，复杂程度有所区别，当然，基于此的分析手法和方式也会有所不同。

伪时间是衡量单个细胞在细胞分化等过程中取得了多大进展的指标。在许多生物学过程中，细胞并不是完全同步的。在细胞分化等过程的单细胞表达研究中，捕获的细胞在分化方面可能分布广泛。也就是说，在同一时间捕获的细胞群中，有些细胞可能已经很长时间了，而有些细胞甚至还没有开始这个过程。当想要了解在细胞从一种状态转换到另一种状态时所发生的调节更改的顺序时，这种异步性会产生主要问题。跟踪同时捕获的细胞间的表达可以产生对基因动力学一个大致的认识，该基因表达的明显变异性将非常高。Monocle根据每个cell在学习轨迹上的进展对其进行排序，从而缓解了由于异步而产生的问题。Monocle不是跟踪表达式随时间变化的函数，而是跟踪沿轨迹变化的函数，我们称之为伪时间。伪时间是一个抽象的分化单位：它只是一个cell到轨迹起点的距离，沿着最短路径测量。轨迹的总长度是由细胞从起始状态移动到结束状态所经历的总转录变化量来定义的。

生信技能树此前推送过一篇文章：拟时序分析就是差异分析的细节剖析，指出拟时序分析就是差异分析的细节剖析，所有的大样本量差异分析都可以转为拟时序分析，可以加深我们对拟时序分析的理解。

Monocle2是做单细胞拟时分析最有名的R包。
相较还在持续开发中的Monocle3来说，Monocle2更稳定且更倾向于半监督的分析模式，更适合针对感兴趣的细胞亚群做个性化分析。

monocle2 文章链接：https://www.nature.com/articles/nmeth.4402

文章中的核心理论为：每个细胞都可以表示为高维空间中的一个点，在高维空间中，每个维对应着一个有序基因的表达水平。高维数据首先通过几种降维方法，如PCA(默认)、扩散映射等，投射到低维空间。Monocle 2然后在自动选择的一组数据质心上构造一棵生成树（DDRTree算法）。然后，该算法将细胞移动到它们最近的树的顶点，更新顶点的位置以适应细胞，学习新的生成树，并迭代地继续这个过程，直到树和细胞的位置已经收敛。在这个过程中，Monocle 2保持了高维空间和低维空间之间的可逆映射，从而既学习了轨迹，又降低了数据的维数。一旦Monocle 2学会了树，用户就会选择一个tip作为根。计算每个单元的伪时间作为其沿树到根的测地线距离，并根据主图自动分配其分枝。因为monocle2学习树结构，与其他方法相比，分支结构自动出现。当它更新细胞位置并细化树时，monocle2简化了轨迹的结构，修剪了小的分支，这样最终的轨迹只保留了描述细胞状态显著差异的分支。

monocle2官网：http://cole-trapnell-lab.github.io/monocle-release/docs/

摘要

单细胞水平的研究使人们可以描述复杂生理过程和高度异质性细胞群体的转录调控。这些研究有助于发现识别特定细胞亚型的基因、标记生物过程中间状态的基因，以及在两种不同的细胞命运之间过渡态的基因。在许多单细胞研究中，单个细胞以不同步的方式执行基因表达过程。实际上，每个细胞都是正在研究的转录过程的一个瞬间。Monocle包是分析单细胞测序的工具。
Monocle引入了在伪时间(拟时间)内对单个细胞排序的策略，利用单个细胞的非同步进程，将它们置于与细胞分化等生物学过程相对应的轨迹上。Monocle利用先进的机器学习技术(反向图嵌入)从单细胞数据中学习显式的主图（展现细胞转录特征相似性关系的图，Monocle2使用DDTree降维图，Monocle3使用UMAP降维图）来对细胞进行排序，Monocle的机器学习算法可以依据上述降维图形，学习描述细胞如何从一种状态过渡到另一种状态的轨迹。Monocle假设轨迹是树状结构，一端是“根”，另一端是“叶”。一个细胞在生物过程的开始，从根开始沿着主干进行，直到它到达第一个分支。然后，该细胞必须选择一条路径，并沿着树移动越来越远，直到它到达一片叶子。一个细胞的假时间值是它返回根所需的距离。降维方面monocle与seurat的过程大同小异，首先进行数据标准化，其次选择部分基因代表细胞转录特征 ，最后选用适当的算法降维。这可以强大而准确地解决复杂的生物过程。
Monocle也可以进行聚类(即使用t-SNE和密度峰值聚类)和差异基因表达测试，使人们能够识别在不同状态下差异表达的基因，沿着生物过程以及不同的细胞命运时基因表达的变化。Monocle是专为单细胞RNA-Seq研究设计的，但也可以用于其他分析。

Introduction

Monocle 2包括新的和改进的算法用于细胞分类和计数，执行细胞亚群之间的差异表达分析，以及细胞轨迹重建。

Monocle主要可以进行以下三种分析：

• 细胞的聚类、分类和计数。 Single-cell RNA-Seq experiments allow you to discover new (and possibly rare) subtypes of cells. Monocle helps you identify them.
• 重建单细胞轨迹。 In development, disease, and throughout life, cells transition from one state to another. Monocle helps you discover these transitions.
• 差异表达分析。 Characterizing new cell types and states begins with comparing them to other, better understood cells. Monocle includes a sophisticated but easy to use system for differential expression.

首先，Monocle 2使用一种简单的、无偏的和高度可扩展的统计程序来选择具有轨迹进展特征的基因。然后，它采用了一类流形学习算法，旨在在高维单细胞RNA-seq数据中嵌入一个主图。以前的方法是通过启发式分析细胞之间的成对距离来推断分支结构，而Monocle 2可以使用这张图来直接识别发育的命运决定。我们已经通过广泛的基准测试证明，Monocle 2优于其他工具，如Wishbone，而不需要用户指定轨迹的结构。

Monocle的使用

分析流程示意图：

1. 安装

install.packages("devtools")
devtools::install_github("cole-trapnell-lab/monocle-release@develop")

或者

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("monocle")

2. 创建CellDataSet

2.1 将Seurat object中数据提取来创建⚠️
数据准备：pbmc3k数据下载（轨迹分析的前提是待分析的细胞有紧密的发育关系，PBMC细胞不是很好的的示例数据，在此仅作为演示。）
由于该数据集没有对细胞类型进行注释，因此我们参考seurat标准流程对这个数据集进行注释。

library(dplyr)
library(Seurat)
library(patchwork) #用来做拼图的包，后面的p1|p1|p3在一张图上展示三个图就是这个包的功劳

# Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
# Initialize the Seurat object with the raw (non-normalized data).
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
pbmc <- RunUMAP(pbmc, dims = 1:10)
new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono", 
    "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
pbmc[['cell_type']] <- pbmc@active.ident #将注释结果添加到metadata
levels(pbmc)
saveRDS(pbmc, file = "pbmc.rds")

运行完上述代码，得到注释好的数据集。

从Seurat对象中提取构建CDS对象所需要的3个输入文件：表达矩阵信息、基因信息和表型信息
⚠️monocle输入的是count矩阵（不建议使用data矩阵）

library(monocle)
pbmc <- readRDS("pbmc.rds") #导入注释好的seurat对象（已注释）

##提取表型信息--细胞信息(建议载入细胞的聚类或者细胞类型鉴定信息、实验条件等信息)
expr_matrix <- as(as.matrix(pbmc@assays$RNA@counts), 'sparseMatrix')
##提取表型信息到p_data(phenotype_data)里面 
p_data <- pbmc@meta.data 
p_data$celltype <- pbmc@active.ident  ##整合每个细胞的细胞鉴定信息到p_data里面。如果已经添加则不必重复添加
##提取基因信息 如生物类型、gc含量等
f_data <- data.frame(gene_short_name = row.names(pbmc),row.names = row.names(pbmc))
##expr_matrix的行数与f_data的行数相同(gene number), expr_matrix的列数与p_data的行数相同(cell number)

#构建CDS对象
pd <- new('AnnotatedDataFrame', data = p_data) 
fd <- new('AnnotatedDataFrame', data = f_data)
#将p_data和f_data从data.frame转换AnnotatedDataFrame对象。
cds <- newCellDataSet(expr_matrix,
                      phenoData = pd,
                      featureData = fd,
                      lowerDetectionLimit = 0.5,
                      expressionFamily = negbinomial.size())

FPKM/TPM值通常是对数正态分布的，而UMIs或读计数使用负二项更好地建模。要处理计数数据，需要将负二项分布指定为newCellDataSet的expressionFamily参数：

negbinomial.size()和negbinomial()：输入的表达矩阵为UMI,一般适用于10x的数据；negbinomial()的结果更准确，但是计算更耗时；一般建议采用negbinomial.size()。
稀疏矩阵用negbinomial.size()
tobit():适用于输入的表达矩阵为FPKM或者TPM, 构建monocle2的class时会自动进行log化计算
gaussianff():输入为log化后的FPKM或者TPM。(目前在单细胞数据中，FPKM已不多用，smart-seq2平台数据一般采用TPM)

2.2 直接读取表达矩阵来创建

library(data.table)
##读取数据
data <- fread("fpkm.txt",data.table = F,header = T)
pd <-  fread("metadata.txt",data.table = F,header = T)
fd <- fread("gene_annotations.txt",data.table = F,header = T)
##创建
pd <- new("AnnotatedDataFrame", data = pd)
fd <- new("AnnotatedDataFrame", data = fd)
HSMM <- newCellDataSet(as.matrix(data),
                       phenoData = pd, featureData = fd,
                       expressionFamily = tobit())
###如果数据量大，建议转化为稀疏矩阵
HSMM <- newCellDataSet(as(as.matrix(data), "sparseMatrix"),
                       phenoData = pd, 
                       featureData = fd,
                       expressionFamily = tobit())

2.3 将Seurat对象直接转化为CellDataSet对象

importCDS(pbmc)

如果我们想要从Seurat对象或SCESet中导入所有的插槽，我们可以设置参数'import_all'为TRUE。#(默认为FALSE或只保留最小数据集)

3. 估计size factor和离散度

size facotr帮助我们标准化细胞之间的mRNA的差异。
离散度值可以帮助我们进行后续的差异分析。
（类似于seurat的数据归一化处理）

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

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

4. 过滤低质量的细胞

大多数单细胞工作流程至少会包含一些由死细胞或空孔组成的库。同样重要的是要删除doublets:由两个或多个细胞意外生成的库。这些细胞可以破坏下游步骤，如伪时间排序或聚类。要知道一个特定的基因有多少个表达，或者一个给定的细胞有多少个基因表达，通常是很方便的。Monocle提供了一个简单的函数来计算这些统计数据。
因为Seurat已经完成细胞过滤，此步可省略
但由于Seurat是通过基因表达量来对细胞进行过滤。因此，我们可以通过下面的代码用表达某基因的细胞的数目对基因进行过滤，从而得到后续操作需要的基因。

cds <- detectGenes(cds, min_expr = 0.1) #这一操作会在fData(cds)中添加一列num_cells_expressed
print(head(fData(cds)))#此时有13714个基因
expressed_genes <- row.names(subset(fData(cds),
    num_cells_expressed >= 10)) #过滤掉在小于10个细胞中表达的基因，还剩11095个基因。

5. 细胞分类

Monocle官网教程提供了4个分类方法：
Classifying cells by type
Clustering cells without marker genes
Clustering cells using marker genes
Imputing cell type
建议先将细胞注释好再进行Monocle分析，不建议使用monocle做细胞分类。

6. 轨迹定义基因选择及可视化和构建轨迹

The ordering workflow
Step 1: choosing genes that define progress
Step 2: reducing the dimensionality of the data
Step 3: ordering the cells in pseudotime

Step 1: 选择定义过程的基因
推断单细胞轨迹是一个机器学习问题。在单细胞RNA-Seq中，低水平表达的基因通常非常嘈杂，但有些基因包含有关细胞状态的重要信息。因此，轨迹推断的第一步就是选择Monocle将用作机器学习方法输入的基因。这叫做特征选择，它对轨迹的形状有很大的影响。
Monocle主要基于关键基因的表达模式，通过学习每个细胞必须经历的基因表达变化的序列，根据拟时间值中对单个细胞进行排序，模拟出时间发育过程的动态变化。而这个排序技术表现是一种在低维空间排布高维数据的降维技术。
Monocle提供了多种工具来选择基因，这些基因将产生一个健壮、准确和具有生物学意义的轨迹。你可以使用这些工具来执行一个完全“无监督”的分析，在这个分析中，Monocle不知道你认为哪个基因是重要的。或者，你可以利用marker gene知识来定义生物学进展，从而形成Monocle的轨迹。我们认为这种模式是“半监督”的分析。

• Monocle官网教程提供了4个选择方法：
  选择发育差异表达基因
  选择clusters差异表达基因
  选择离散程度高的基因
  自定义发育marker基因

前三种都是无监督分析方法，细胞发育轨迹生成完全不受人工干预；最后一种是半监督分析方法，可以使用先验知识辅助分析。

##使用seurat选择的高变基因⚠️
express_genes <- VariableFeatures(pbmc)
cds <- setOrderingFilter(cds, express_genes)
plot_ordering_genes(cds)
##使用clusters差异表达基因
deg.cluster <- FindAllMarkers(pbmc)
express_genes <- subset(deg.cluster,p_val_adj<0.05)$gene
cds <- setOrderingFilter(cds, express_genes)
plot_ordering_genes(cds)
##使用monocle选择的高变基因⚠️
disp_table <- dispersionTable(cds)
disp.genes <- subset(disp_table, mean_expression >= 0.1 & dispersion_empirical >= 1 * dispersion_fit)$gene_id
cds <- setOrderingFilter(cds, disp.genes)
plot_ordering_genes(cds)

理想情况下，我们希望尽可能少地使用正在研究的系统生物学的先验知识。我们希望从数据中发现重要的排序基因，而不是依赖于文献和教科书，因为这可能会在排序中引入偏见。我们将从一种更简单的方法开始，但是我们通常推荐一种更复杂的方法，称为“dpFeature”。

#这一步输入的expressed_genes来自于步骤4。
#⚠️⚠️后续分析使用的是该方法
#也可输入seurat筛选出的高变基因：expressed_genes <- VariableFeatures(pbmc) 
diff <- differentialGeneTest(cds[expressed_genes,],fullModelFormulaStr="~cell_type",cores=1) 
#~后面是表示对谁做差异分析的变量，理论上可以为p_data的任意列名
head(diff)

##差异表达基因作为轨迹构建的基因,差异基因的选择标准是qval<0.01,decreasing=F表示按数值增加排序
deg <- subset(diff, qval < 0.01) #选出2829个基因
deg <- deg[order(deg$qval,decreasing=F),]
head(deg)

##差异基因的结果文件保存
write.table(deg,file="train.monocle.DEG.xls",col.names=T,row.names=F,sep="\t",quote=F)

## 轨迹构建基因可视化
ordergene <- rownames(deg) 
cds <- setOrderingFilter(cds, ordergene)  
#这一步是很重要的，在我们得到想要的基因列表后，我们需要使用setOrderingFilter将它嵌入cds对象，后续的一系列操作都要依赖于这个list。
#setOrderingFilter之后，这些基因被储存在cds@featureData@data[["use_for_ordering"]]，可以通过table(cds@featureData@data[["use_for_ordering"]])查看
pdf("train.ordergenes.pdf")
plot_ordering_genes(cds)
dev.off()
#出的图黑色的点表示用来构建轨迹的差异基因，灰色表示背景基因。红色的线是根据第2步计算的基因表达大小和离散度分布的趋势(可以看到，找到的基因属于离散度比较高的基因)

⚠️选择的用于排序的基因数目一般在2000左右比较合适

gene数太多的话也可以选择top基因

ordergene <- row.names(deg)[order(deg$qval)][1:400]

这样用的文献：肿瘤浸润髓系细胞的pan-cancer单细胞转录图谱

Step 2: 降维
一旦细胞有序排列，我们就可以在降维空间中可视化轨迹。所以首先选择用于细胞排序的基因，然后使用反向图嵌入(DDRTree)算法对数据进行降维。

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

Step 3: 拟时间轴轨迹构建和在拟时间内排列细胞
将表达数据投射到更低的维度空间，通过机器学习描述细胞如何从一种状态过渡到另一种状态的轨迹。假设轨迹具有树状结构，一端是“根”，另一端是“叶”。尽可能地将最佳树与数据匹配起来。这项任务被称为“歧管学习”，在生物过程的开始阶段，细胞从根部开始，沿着主干前进，直到到达第一个分支（如果有的话）。然后，细胞必须选择一条路径，沿着树走得越来越远，直到到达一片叶子。一个细胞的伪时间值是它回到根的距离。

根据order gene的表达趋势，将细胞排序并完成轨迹构建

cds <- orderCells(cds)
#⚠️使用root_state参数可以设置拟时间轴的根，如下面的拟时间着色图中可以看出，左边是根。根据state图可以看出，根是State1，若要想把另一端设为根，可以按如下操作
#cds <- orderCells(cds, root_state = 5) #把State5设成拟时间轴的起始点

可视化：根据cds@phenoData@data中的表型信息(metadata)对细胞上色

1. 以pseudotime值上色 (Pseudotime是monocle2基于细胞基因表达信息计算的概率，表示时间的先后。)

pdf("train.monocle.pseudotime.pdf",width = 7,height = 7)
plot_cell_trajectory(cds,color_by="Pseudotime", size=1,show_backbone=TRUE) 
dev.off()

图中123是分叉点，不代表特别意义，也不代表时间先后。这张图按照sudotime，时间先后是从左往右，左边是起点（root)。起点的设置可以使用orderCells的root_state参数，将右侧设置为起点。

2. 以细胞类型上色

pdf("train.monocle.celltype.pdf",width = 7,height = 7)
plot_cell_trajectory(cds,color_by="cell_type", size=1,show_backbone=TRUE)
dev.off()

3. 以细胞状态上色

pdf("train.monocle.state.pdf",width = 7,height = 7)
plot_cell_trajectory(cds, color_by = "State",size=1,show_backbone=TRUE)
dev.off()

state的多少是monocle算出来的，不能调整，与输入的用于轨迹学习的基因有关。分叉和顶点之间或者顶点和顶点之间为一个state，与发育轨迹时间先后没有关系，与细胞类型也不完全相关。

4. 按照seurat分群排序细胞

pdf("seurat.clusters.pdf",width = 7,height = 7)
plot_cell_trajectory(cds, color_by = "seurat_clusters")
dev.off()

5. 以细胞状态上色（拆分）“分面”轨迹图，以便更容易地查看每个状态的位置。

pdf("train.monocle.state.faceted.pdf",width = 10,height = 7)
plot_cell_trajectory(cds, color_by = "State") + facet_wrap("~State", nrow = 1)
dev.off()

⚠️这一部分的配色可以直接对接ggsci，也可以使用scale_color_manual()自己设置颜色

library(ggsci)
p1=plot_cell_trajectory(cds, color_by = "cell_type")  + scale_color_npg() 
p2=plot_cell_trajectory(cds, color_by = "State")  + scale_color_nejm()
colour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB","#98FB98","#F08080")
p3=plot_cell_trajectory(cds, color_by = "State")  + scale_color_manual(values = colour)
p1|p2|p3

觉得这种轨迹不太直观的也可以画成树形图

p1 <- plot_cell_trajectory(cds, x = 1, y = 2, color_by = "celltype") + 
  theme(legend.position='none',panel.border = element_blank()) + #去掉第一个的legend
  scale_color_manual(values = colour) 
p2 <- plot_complex_cell_trajectory(cds, x = 1, y = 2,
                                   color_by = "celltype")+
  scale_color_manual(values = colour) +
  theme(legend.title = element_blank()) 
p1|p2

还可以画沿时间轴的细胞密度图

library(ggpubr)
df <- pData(cds) 
## pData(cds)取出的是cds对象中cds@phenoData@data的内容
View(df)
ggplot(df, aes(Pseudotime, colour = cell_type, fill=cell_type)) +
  geom_density(bw=0.5,size=1,alpha = 0.5)+theme_classic2()

从上面的矩阵中也能提取不同state的细胞，分别画图

手动设置颜色注意

ClusterName_color_panel <- c(
  "Naive CD4 T" = "#DC143C", "Memory CD4 T" = "#0000FF", "CD14+ Mono" = "#20B2AA",
  "B" = "#FFA500", "CD8 T" = "#9370DB", "FCGR3A+ Mono" = "#98FB98",
  "NK" = "#F08080", "DC" = "#0000FF", "Platelet" = "#20B2AA"
)
ggplot(df, aes(Pseudotime, colour = cell_type, fill=cell_type)) +
  geom_density(bw=0.5,size=1,alpha = 0.5)+theme_classic2()+ scale_fill_manual(name = "", values = ClusterName_color_panel)+scale_color_manual(name = "", values = ClusterName_color_panel)

提取感兴趣的细胞（进行后续分析）

#比如对State7的细胞感兴趣
pdata <- Biobase::pData(cds)
s.cells <- subset(pdata, State=="7") %>% rownames()
save(s.cells, file = "Monocle_state7.rda")

保存结果

write.csv(pData(cds), "pseudotime.csv")
save(cds, file = "cds.rda")

7. 指定基因的可视化

直接查看一些基因随细胞状态等的表达变化

##选择前4个top基因并将其对象取出
keygenes <- head(ordergene,4)
cds_subset <- cds[keygenes,]
##可视化：以state/celltype/pseudotime进行
p1 <- plot_genes_in_pseudotime(cds_subset, color_by = "State")
p2 <- plot_genes_in_pseudotime(cds_subset, color_by = "cell_type")
p3 <- plot_genes_in_pseudotime(cds_subset, color_by = "Pseudotime")
plotc <- p1|p2|p3
ggsave("Genes_pseudotimeplot.pdf", plot = plotc, width = 16, height = 8)

#指定基因
s.genes <- c("SELL","CCR7","IL7R", "CD84","CCL5","S100A4")
p1 <- plot_genes_jitter(cds[s.genes,], grouping = "State", color_by = "State")
p2 <- plot_genes_violin(cds[s.genes,], grouping = "State", color_by = "State")
p3 <- plot_genes_in_pseudotime(cds[s.genes,], color_by = "State")
plotc <- p1|p2|p3
ggsave("Genes_Jitterplot.pdf", plot = plotc, width = 16, height = 8)

由于state1-7并不代表分化时间先后，因此这张图还要结合其他的图来看

⚠️拟时序展示单个基因表达量

colnames(pData(cds))
pData(cds)$CCL5 = log2( exprs(cds)['CCL5',]+1)
p1=plot_cell_trajectory(cds, color_by = "CCL5")  + scale_color_gsea()
pData(cds)$S100A4 = log2(exprs(cds)['S100A4',]+1)
p2=plot_cell_trajectory(cds, color_by = "S100A4")    + scale_color_gsea()
library(patchwork)
p1+p2

8. 寻找拟时相关的基因（拟时差异基因）

使用回归算法
注意：不要使用多核运算，经常会出现警告
Monocle的主要工作是通过生物过程（如细胞分化）将细胞按顺序排列，而不知道要提前查看哪些基因。一旦这样做了，你就可以分析细胞，找到随着细胞进展而变化的基因。

官方给出的差异分析有三大方法：
1、Basic Differential Analysis
2、Finding Genes that Distinguish Cell Type or State
3、Finding Genes that Change as a Function of Pseudotime
我们重点关注第三个：根据伪时间功能寻找差异基因

sm.ns函数指出Monocle应该通过表达式值拟合自然样条曲线，以帮助它将表达式的变化描述为进程的函数。

• 寻找拟时差异基因（qvalue体现基因与拟时的密切程度）绘制热图

#这里是把排序基因（ordergene）提取出来做回归分析，来找它们是否跟拟时间有显著的关系
#如果不设置，就会用所有基因来做它们与拟时间的相关性
Time_diff <- differentialGeneTest(cds[ordergene,], cores = 1, 
                                  fullModelFormulaStr = "~sm.ns(Pseudotime)")
Time_diff <- Time_diff[,c(5,2,3,4,1,6,7)] #把gene放前面，也可以不改
write.csv(Time_diff, "Time_diff_all.csv", row.names = F)
Time_genes <- Time_diff %>% pull(gene_short_name) %>% as.character()
p=plot_pseudotime_heatmap(cds[Time_genes,], num_clusters=4, show_rownames=T, return_heatmap=T)
ggsave("Time_heatmapAll.pdf", p, width = 5, height = 10)

⚠️ ⚠️⚠️前面通过设置num_clusters将热图聚成了四个cluster，如过想要把每个cluster的基因单独提出来做分析

p$tree_row
# Call:
# hclust(d = d, method = method)
# Cluster method   : ward.D2 
# Number of objects: 2829 
clusters <- cutree(p$tree_row, k = 4)
clustering <- data.frame(clusters)
clustering[,1] <- as.character(clustering[,1])
colnames(clustering) <- "Gene_Clusters"
table(clustering)
# 1    2    3    4 
# 570 1031  506  722 
write.csv(clustering, "Time_clustering_all.csv", row.names = F)

这样就把每个基因属于哪个cluster提取出来了，后续可以做每个cluster的富集分析。参考：拟时序分析的热图提取基因问题

• 拟时差异基因热图绘制（提取了前100个）

Time_genes <- top_n(Time_diff, n = 100, desc(qval)) %>% pull(gene_short_name) %>% as.character()
p = plot_pseudotime_heatmap(cds[Time_genes,], num_clusters=4, show_rownames=T, return_heatmap=T)
ggsave("Time_heatmapTop100.pdf", p, width = 5, height = 10)

这个图的横轴是拟时间，cluster1的基因是在拟时排序起点高表达的基因，cluster2的基因则是在拟时排序的重点高表达的。cluster数的多少是由plot_pseudotime_heatmap函数中的num_clusters参数定义的

注：plot_pseudotime_heatmap函数可以来可视化所有monocle的假时间依赖性基因。plot_pseudotime_heatmap采用CellDataSet对象（通常只包含重要基因的子集），并生成平滑的表达曲线，非常类似于plot_genes_in_pseudotime。然后它对这些基因进行聚类，并使用pheatmap软件包进行绘图。绘出的热图可以让我们观测到假时间依赖性基因中的不同基因模块在不同的时间内共同变化，能比较好的回答时间序列基因表达中“哪些基因遵循相似的动力学趋势”这一常见问题。

• 显著差异基因按热图结果排序并保存

hp.genes <- p$tree_row$labels[p$tree_row$order]
Time_diff_sig <- Time_diff[hp.genes, c("gene_short_name", "pval", "qval")]
write.csv(Time_diff_sig, "Time_diff_sig.csv", row.names = F)

另外：也可手动选择基因来绘制热图，查看其表达模式

marker_genes <- row.names(subset(fData(cds),
                                 gene_short_name %in% c("MEF2C", "MEF2D", "MYF5",
                                                        "ANPEP", "PDGFRA","MYOG",
                                                        "TPM1",  "TPM2",  "MYH2",
                                                        "MYH3",  "NCAM1", "TNNT1",
                                                        "TNNT2", "TNNC1", "CDK1",
                                                        "CDK2",  "CCNB1", "CCNB2",
                                                        "CCND1", "CCNA1", "ID1")))

diff_test_res <- differentialGeneTest(cds[marker_genes,],
                                      fullModelFormulaStr = "~sm.ns(Pseudotime)")
sig_gene_names <- row.names(subset(diff_test_res, qval < 0.1))
plot_pseudotime_heatmap(cds[sig_gene_names,],
                        num_clusters = 6,
                        cores = 1,
                        show_rownames = T)

9. 单细胞轨迹的“分支”分析

⚠️上一步寻找拟时相关的基因是全局的，找拟时起点和终点相关的基因，这一步则是寻找和分叉点相关的基因。

单细胞轨迹常常包括分支。这些分支的产生是因为细胞执行不同的基因表达程序。在发育过程中，当细胞做出命运选择时，分支出现在轨迹中：一个发育谱系沿着一条路径前进，而另一个谱系产生第二条路径。Monocle包含分析这些分支事件的广泛功能。Monocle提供了一个特殊的统计测试：分支表达式分析建模，或BEAM。
BEAM(Branched expression analysis modeling)是一种统计方法，用于寻找以依赖于分支的方式调控的基因。

plot_cell_trajectory(cds, color_by = "State")

结合前面color_by="Pseudotime"的图，这个图默认的是最左边是root(pseudotime 0)，所以这个图是status1是最早的，在分叉2处分为status1和status7，status2又随后在分叉3处分为status3和status4，status4随后又在分叉1处分出status6和status5。

BEAM进行统计分析

BEAM_res <- BEAM(cds[ordergene,], branch_point = 1, cores = 2) 
#这里用的是ordergene，也就是第六步dpFeature找出来的基因。如果前面用的是seurat的marker基因，记得改成express_genes
#BEAM_res <- BEAM(cds, branch_point = 1, cores = 2) #对2829个基因进行排序，运行慢
BEAM_res <- BEAM_res[order(BEAM_res$qval),]
BEAM_res <- BEAM_res[,c("gene_short_name", "pval", "qval")]
head(BEAM_res)
#           gene_short_name         pval         qval
# CD79A               CD79A 2.757676e-73 7.782161e-70
# TCL1A               TCL1A 1.574889e-65 2.222168e-62
# IGLL5               IGLL5 2.356778e-64 2.216942e-61
# S100A9             S100A9 1.504319e-58 1.061297e-55
# S100A8             S100A8 6.028175e-57 3.402302e-54
# LINC00926       LINC00926 3.180527e-55 1.495908e-52
write.csv(BEAM_res, "BEAM_res.csv", row.names = F)
plot_genes_branched_heatmap(cds[row.names(subset(BEAM_res,
                                                  qval < 1e-4)),],
                            branch_point = 1, #绘制的是哪个分支
                            num_clusters = 4, #分成几个cluster，根据需要调整
                            cores = 1,
                            use_gene_short_name = T,
                            show_rownames = T)#有632个gene，太多了

该热图显示的是同一时间点两个谱系的变化，热图的列是伪时间的点，行是基因。这张图最上面的条条，灰色的代表分叉前，左边红色代表左边这个cell fate，右边蓝色代表右边这个cell fate，从热图中间往右读，是伪时间的一个谱系，往左是另一个谱系。基因是被按照等级聚类的，需要结合生物学知识来进行解读。

⚠️⚠️因为前面设置的branch_point = 1，根据按照State绘制的trajectory图，在1这个分枝上分出的是5和6这两个status，所以上图比较的是state5和state6这两个branch，而且上图中的pre-branch使用的是4, 2, 7, 1这四个status（BEAM tries to traverse backward from the cell on the branch point all the way back to the root cell (the cell with pseudotime 0) and use all those cells as the the pre-branch）。而且，上图左边的Cell fate 1指的是state5，Cell fate 2指的是state6（Cell fate 1 corresponds to the state with small id while cell fate 2 corresponds to sate with bigger id ）。

#选前100个基因可视化
BEAM_genes <- top_n(BEAM_res, n = 100, desc(qval)) %>% pull(gene_short_name) %>% as.character()
p <- plot_genes_branched_heatmap(cds[BEAM_genes,],  branch_point = 1, 
                                 num_clusters = 3, show_rownames = T, return_heatmap = T)
ggsave("BEAM_heatmap.pdf", p$ph_res, width = 6.5, height = 10)

这个图的解读：
· How to map cell fate to branches?
· monocle2 拟时间分支点分析结果解读

#显著差异基因(top100)按热图结果排序并保存
##如果要所有的差异基因，就把前面所632个基因的热图存为p
hp.genes <- p$ph_res$tree_row$labels[p$ph_res$tree_row$order]
BEAM_sig <- BEAM_res[hp.genes, c("gene_short_name", "pval", "qval")]
write.csv(BEAM_sig, "BEAM_sig.csv", row.names = F)

选择上面热图中或显著差异基因中感兴趣的基因进行可视化

head(BEAM_sig)
#        gene_short_name         pval         qval
# STX11            STX11 2.928177e-20 1.180474e-18
# CEBPD            CEBPD 2.137741e-22 1.058369e-20
# TYROBP          TYROBP 3.157539e-54 1.272939e-51
# FCER1G          FCER1G 1.094272e-45 2.573364e-43
# SRGN              SRGN 7.564014e-24 4.356255e-22
# CTSD              CTSD 2.666673e-17 7.601364e-16
genes <- row.names(subset(fData(cds),
                          gene_short_name %in% c( "SITX11", "CEBPD", "TYROBP")))

plot_genes_branched_pseudotime(cds[genes,],
                               branch_point = 1,
                               color_by = "State",
                               ncol = 1)

可以看到，在pre-branch细胞向 5/6这两个分叉分化时，这两个基因都是在state5中高表达的，也就是state5的分支相关基因。

⚠️注意：不同分组间的细胞尽量不要放在一起做轨迹分析，同一组的生物学重复可以一起分析。明显没有生物学相关性的细胞也不要放在一起做轨迹分析。

10. 文献里的monocle

有时候细胞群不会完美的按照分叉排列。如下，只要一类细胞占某一个分叉细胞的大部分也可以。

Single-cell reconstruction of the adult human heart during heart failure and recovery reveals the cellular landscape underlying cardiac function, Extended Data Fig. 8a-e

分析过程中使用的主要函数汇总