monocle分析及结果解读

1.写在前面的话：

近年来，由于细胞的异质性及发育分化等相关的问题越来越被研究者们所关注，单细胞转录组分析为研究异质细胞群的复杂生物学过程提供了方法和工具。每一个细胞进行转录组测序时就是细胞发育过程中的快照，单细胞拟时间分析软件Monocle2是基于R语言的安装包，其功能基于单细胞转录组的表达矩阵，通过无监督学习（Reversed Graph Embedding算法）的方式将细胞置于发育轨迹的不同分支上，从而模拟细胞群体生物学过程。也就是我们经常说的拟时序（pseudotime）分析，又称细胞轨迹（cell trajectory）分析。通过拟时分析可以推断出发育过程细胞的分化轨迹或细胞亚型的演化过程，在发育相关研究中使用频率较高。
模拟细胞的分化轨迹的软件，最常用的软件为Monocle2。该算法不仅能模拟细胞的发育轨迹，同时也能对细胞进行聚类（t-SNE）。通过聚类获得不同状态下的差异基因，分析影响分支形成的关键基因及其功能，对研究相关生物学问题有指导性的作用。
Monocle2主要基于关键基因的表达模式，通过学习每个细胞必须经历的基因表达变化的序列，根据拟时间值中对单个细胞进行排序，模拟出时间发育过程的动态变化。而这个排序技术表现是一种在低维空间排布高维数据的降维技术。（具体描叙请参考：https://www.meiwen.com.cn/subject/zlgsvqtx.html），目前monocle已经升级至monocle3，但在结果分析上monocle3的可读性不如monocle2，因此本研究中主推monocle2这个版本。
在做拟时序分析之前，先要明白几个研究目的：

1. 对什么类型的细胞群进行拟时许分析？这类细胞群在不同的分化轨迹或细胞亚型的区别是否明显？
2. monocle做拟时序分析时，可与seurat聚类结果进行对接，因此可对seurat的分群结果进行定性的解释。

3. 可分析比较组间的差异或者多组样本间拟时间分析差别。
   分析流程图如下所示：

   
   image.png

2. 软件安装

先保证服务器上装有R ，我这里安装的版本是version 3.5.1。

biocLite()
biocLite("monocle")
install.packages("devtools")
devtools::install_github("cole-trapnell-lab/monocle-release@develop")
biocLite(c("DDRTree", "pheatmap"))

3. 输入文件

对于输入文件，有3种类型的格式：
一种是单细胞运行获得的3个结果文件。格式如下表所示：
准备需要进行拟时间分析数据的三个文件：
1)表达量文件：exprs：基因在所有细胞中的count值矩阵。
格式示例：

注意：该文件的表头必须以“%%**\n%”的形式出现，否则就会出现报错。
2)表型文件 phenoData:
格式示例，该文件即为每个细胞的barcode信息。：

3. featureData：
   该文件为所有基因名信息，有两种格式，一种是Ensemble Id和Symbol Id的组合，一种均是Symbol Id,但要保证第二列一定是Symbol Id。
   格式示例：

   
   
   

   第二种输入数据的格式是直接转换seurat对象为CellDataSet格式。
   第三种输入数据的格式是只有基因和细胞构成的表达量矩阵。

4. 特征向量转化为CellDataSet对象

方法一：3个特征文件转换成CellDataSet格式。命令行如下所示：

gbm <- load_cellranger_matrix(my_dir) ###my_dir中包括有三个文件。
HSMM <- newCellDataSet(exprs(gbm),phenoData = new("AnnotatedDataFrame", data = pData(gbm)),featureData = new("AnnotatedDataFrame", data = my_feat),lowerDetectionLimit = 0.5,expressionFamily = negbinomial.size())

方法二：直接通过seurat生成的rds文件的方式进行转换。

由于seurat软件包的升级，大量的函数命名发生了变化，对2和3这两种类型的版本的操作也是不一样。

spleen <-readRDS("seurat.analysis.rds") 

Seurat2.4的版本：

monocle_cds <-importCDS(spleen,import_all=TRUE) 

Seurat3.0及更高版本：

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

此处需要提供绝对的count值，即UMI矩阵。

pd <- new('AnnotatedDataFrame', data =data) 
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data) 
fd <- new('AnnotatedDataFrame', data = fData)    
HSMM <- newCellDataSet(data, 
​            phenoData = pd, 
​            featureData = fd, 
​            lowerDetectionLimit = 0.5, 
​            expressionFamily = negbinomial.size())

方法三：仅仅只有表达量数据的时候如何进行处理？

如果你没有3个标准格式的文件，也没有生成rds文件，仅仅只有表达量矩阵，具体原理和有rds文件一致，可按照如下方法进行处理。

HSMM_expr_matrix <- read.csv("fpkm_matrix.csv",header=T)
pd <- new('AnnotatedDataFrame', data =data.frame(colnames(HSMM_expr_matrix)))
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
fd <- new('AnnotatedDataFrame', data = fData)    
HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix) , phenoData = pd, featureData = fd, lowerDetectionLimit = 0.5, expressionFamily = negbinomial.size())

5. 对monocle对象进行归一化。

如果你采用的是第二种输入数据进行分析，即使你原先的Seurat对象已经归一化过了，官方推荐在monocle中重新进行归一化处理。

HSMM <- estimateSizeFactors(HSMM)   
HSMM <- estimateDispersions(HSMM)
HSMM <- detectGenes(HSMM, min_expr = 0.1) ###过滤掉低质量的基因。

6.数据降维，特征基因的选择。

在进行降维聚类之前，我们先获得高表达的基因集作为每个聚类的用来order的Feature基因。当然我们可以使用所有的基因，但一些表达量特别低的基因提供的聚类信号往往会制造分析噪音，Feature基因的选择性很多，一种是我们可以根据基因的平均表达水平来进行筛选，另外我们也可以选择细胞间异常变异的基因。这些基因往往能较好地反映不同细胞的状态。以平均表达量高于0.1来进行筛选。

disp_table <- dispersionTable(HSMM)
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1) 
HSMM_ myo <- setOrderingFilter(HSMM, unsup_clustering_genes$gene_id)
plot_ordering_genes(HSMM)

绘制用于order基因和非order的平均表达量的分布点图。这里不做展示。

7.聚类

1.首先根据TSNE进行降维处理

HSMM_myo <- reduceDimension(HSMM_myo, max_components= max_components, num_dim = 3, norm_method = 'log', reduction_method = 'tSNE',verbose = T)

对细胞进行聚类，在Seurat中采用的是分辨率来确定cluster的数目。而monocle中可以直接指定聚类数目。主要指出的这里所聚类获得的cluster数目比我们赋值的要少一个。即当num_clusters=3时，你只获得了2个cluster。具体解释有些难懂，在这里不做过多的解释。

2.获得对细胞排序起到关键作用的基因集。

拟时间分析不仅是要对分析的细胞群进行排序，还要获得觉得细胞排序的关键基因集。这种基因集有两种情况，在有先验知识的情况下，我们可以从系统生物学的角度获得一系列与细胞发育相关的基因集，从而对细胞进行排序，这种方式是最为保险的，但局限性是对未知的发育情况不能进行分析。另外一种情况就是使用无监督聚类方式计算关键基因集。接下来我们采用differentialGeneTest方式获得clustering_DEG_genes（与聚类相关的差异基因集）

clustering_DEG_genes <-differentialGeneTest(HSMM_myo[HSMM_expressed_genes,],fullModelFormulaStr = '~Cluster', cores = 4)

上一步过程是对所有的细胞进行无监督训练，运行时间与细胞数和基因数相关，一般会花很长的时间。可以根据cores的数目进行并行。
differentialGeneTest这个函数测试每个基因的差异表达，取决于伪时间或根据指定的其他协变量。 “ differentialGeneTest”是Monocle的主要差异分析常规， 它接受一个CellDataSet和两个模型公式作为输入，指定由实现的广义谱系模型“ VGAM”包。 也就是说我们可以根据指定’~cluster’或者拟时间值来获得差异基因。差异基因的结果如下表所示：

在这个表格中，我们先看一下表头对应的关键列。第一列没有列名，为基因名称。pval,qval 为差异基因的显著检验打分。num_cells_expressed为表达这个基因的细胞数。use_for_ordering该基因是否是作为排序的差异基因。

3. 细胞排序，一般情况下使用qval排在前1000个基因进行排序。

HSMM_ordering_genes <-row.names(clustering_DEG_genes)[order(clustering_DEG_genes$qval)]
HSMM_myo <-setOrderingFilter(HSMM_myo,ordering_genes = HSMM_ordering_genes)
HSMM_myo <-reduceDimension(HSMM_myo, method = 'DDRTree')
HSMM_myo <-orderCells(HSMM_myo)

根据排序好细胞进行结果可视化。
命令行如下所示：

plot_cell_trajectory(HSMM_myo, color_by = "State",cell_size = 0.75) ####以state进行着色
plot_cell_trajectory(HSMM_myo, color_by = "State",cell_size = 0.75)+facet_wrap(~State, nrow = Cluster_num) ##绘制state的分面图
plot_cell_trajectory(HSMM_myo, color_by = "Pseudotime",cell_size = 0.75) ##根据拟时间值着色
plot_cell_trajectory(HSMM_myo, color_by = "Cluster",cell_size = 0.75) ##根据细胞聚类的进行着色
plot_cell_trajectory(HSMM_myo, color_by = "Cluster",cell_size = 0.75)+facet_wrap(~Cluster, nrow =Cluster_num) ###绘制clusster的分面图
plot_cell_trajectory(HSMM_myo, color_by = "seurat_clusters",cell_size = 0.75) +scale_color_manual(values=col) ##如果有Seurat生的rds文件的话，按照seurat中分的群进行着色，如果不想用ggplot的默认色，可以提供颜色列表col list。
plot_cell_trajectory(HSMM_myo, color_by = "seurat_clusters", cell_size = 0.75) + facet_wrap (~seurat_clusters, nrow =Cluster_num) ###按照seurat中分的群绘制分面图。
plot_cell_trajectory(HSMM_myo, color_by = "stim",cell_size = 0.75) ###按照样本进行着色
plot_cell_trajectory(HSMM_myo, color_by = "stim",cell_size = 0.75)+facet_wrap(~stim, ncol = 2)
##绘制样本着色的分面图。

上述部分结果如下图所示（不包括分面图）

上图表示在主成分中的细胞聚类分布的tsne图。不同颜色代表细胞群的不同细胞命运的分支。

上图表示在主成分中的细胞聚类分布的tsne图。颜色的深浅表示细胞的根据拟时间值的排序。

上图表示依据seurat的cluster ID映射到拟时间的排序上。
接下来可视化展示细胞state相关的差异基因的表达量分布情况，可以根据这些基因来确定细胞的发育方向，下图仅展示qvalue值排在前6个基因。横纵左边意义见坐标轴的描述。

Top6的基因仅展示了部分差异基因与细胞发育状态的关系，不能较好地进行展示。下面的热图采用top50的基因进行展示，当然，这个基因列表也可以自己进行筛选。

上述图片对top50的基因的表达量进行了热图绘制，横轴为每个基因对应的表达量。纵轴为根据拟时间值，一个拟时间值可能代表多个细胞。我们把这些基因分成了3个表达模块，即最左侧的cluster1,2,3, 相同的cluster中的基因表达量比较相似，同样可以根据这个图进行发育方向的确定。

8 分支点的差异分析

在monocle中的分析中发现，细胞群能从一个轨迹分叉成不同的分支点，在科研中，我们会比较关注发生分支的原因是什么，即分析分支点之间的差异。Monocle采用分支表达式分析建模，主要是BEAM函数，可以将分叉过程重构为一个分支轨迹，从而分析不同细胞命运下的差异。
命令行如下：

my_pseudotime_de <- differentialGeneTest(HSMM_myo, fullModelFormulaStr= " ~sm.ns(Pseudotime)" , cores = 4) ###这里是对拟时间值来计算差异基因。依然花的时间比较长。

举个例子：我们分析branch_point = 1这个分支处的细胞命名分叉是如何进行的。

BEAM_res <- BEAM(HSMM_myo, branch_point = 1, cores = 4) ###选择branch_point =1。

即下图中所示。下图只有1个分支点，即分析state1,2,3 这三个state的差异。

对影响分析的基因根据qvalue进行排序。

BEAM_res <- BEAM_res[order(BEAM_res$qval),]
BEAM_res <- BEAM_res[,c("gene_short_name", "pval", "qval")]

绘制与分支相关的基因热图。

my_branched_heatmap <- plot_genes_branched_heatmap (HSMM_myo[row.names (subset (BEAM_res, qval < 1e-4)), ], branch_point = 1, num_clusters = 4, cores = 4, use_gene_short_name = TRUE,show_rownames = FALSE,return_heatmap = TRUE)

关键参数为subset (BEAM_res, qval < 1e-4))，挑选基因进行热图绘制，也可以设置成其他的阈值。
branch_point=1，分支点选为1。num_clusters=4，将基因根据表达相似性分成4个模块。结果如下图所示：

每一行代表一个基因。不同的聚类代表不同的基因模块。图例中上面灰色为分之前的细胞，红色为细胞命运1，蓝色为细胞命运2。每一列是指将所有细胞的拟时间值从小到大分成100个bin，用bin的形式代表拟时间值的变化，与细胞数目无关。
另外，可以绘制基因在所有细胞中的表达量在不同分支上的表达量情况。我们可以根据分叉的state位置来判定这个基因对分支的影响。在拟时间分支树构建的过程中，其实是有很多细小的分支，最终被拟合成一条直线，便于可视化的友好性。下图中的Branch表示挑选了两个不同的未被拟合的实际分支。

参考：
https://www.jianshu.com/p/9995cd707002
https://www.jianshu.com/p/66c387e1de3d