Monocle拟时间分析

Monocle可以进行下面三种类型的分析，可以与seurat无缝连接：

1 对细胞进行聚类、分组和计算；
2.进行拟时间分析；
3.差异表达分析；

2.Constructing Single Cell Trajectories

Trajectory step 1: choose genes that define a cell's progress

理论上希望找到一组能够随着研究过程的进展而增加(或减少)表达的基因。同时尽可能少地使用正在研究的系统生物学的先验知识。我们希望从数据中发现重要的排序基因，而不是依赖于文献和教科书，因为这可能会在排序中引入偏见。将从一种简单的方式开始，但是推荐使用另外一种更有效的方法“dpFeature”(后面会提到)。
首先，找到一组基因的有效方法之一是简单地比较过程开始时和结束时收集的细胞（前提是你有这一组信息），并找到差异表达的基因。

#这一种情况适用于有一个明确的时间节点（这里就是换培养基）
diff_test_res <- differentialGeneTest(HSMM_myo[expressed_genes,],
              fullModelFormulaStr = "~Media")
ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))

#确定ordering_genes之后，将它整合到HSMM对象中
HSMM_myo <- setOrderingFilter(HSMM_myo, ordering_genes)
plot_ordering_genes(HSMM_myo)

Trajectory step 2: reduce data dimensionality

#降维处理
HSMM_myo <- reduceDimension(HSMM_myo, max_components = 2,method = 'DDRTree')

Trajectory step 3: order cells along the trajectory

HSMM_myo <- orderCells(HSMM_myo)
plot_cell_trajectory(HSMM_myo, color_by = "Hours")

根据上面的拟时间轨迹，可以看出0时收集的细胞聚集在树的一端，其他时期分布在数的两条分支上。Monocle是不知道树的哪边是起始端，因此往往需要再次调用orderCells函数，使用root_state参数去指定起始端。

#下面就是进行的再次调用orderCells函数的步骤
#首先按照树的分段进行颜色区分（state参数）
plot_cell_trajectory(HSMM_myo, color_by = "State")

可以看到，一共有5种状态，就是有五个颜色

#下面这个函数用来识别哪一段树枝包含时间0的细胞最多?这里不是太理解
GM_state <- function(cds){
  if (length(unique(pData(cds)$State)) > 1){
    T0_counts <- table(pData(cds)$State, pData(cds)$Hours)[,"0"]
    return(as.numeric(names(T0_counts)[which
          (T0_counts == max(T0_counts))]))
  } else {
    return (1)
  }}
HSMM_myo <- orderCells(HSMM_myo, root_state = GM_state(HSMM_myo))
plot_cell_trajectory(HSMM_myo, color_by = "Pseudotime")

#如果你细胞数太多，可以分开看每个细胞亚群所处的位置
plot_cell_trajectory(HSMM_myo, color_by = "State") +
    facet_wrap(~State, nrow = 1)

如果没有时间序列，就需要根据生物学知识找到root。就这个实验而言，一个高度增殖的祖细胞群产生两种有丝分裂的细胞。因此，root就应该为高表达增殖marker的那一群。

blast_genes <- row.names(subset(fData(HSMM_myo),gene_short_name %in% c("CCNB2", "MYOD1", "MYOG")))
plot_genes_jitter(HSMM_myo[blast_genes,],grouping = "State",min_expr = 0.1)

#为了确认顺序是正确的，我们可以选择几个肌源性进展过程中的标记
HSMM_expressed_genes <-  row.names(subset(fData(HSMM_myo),num_cells_expressed >= 10))
HSMM_filtered <- HSMM_myo[HSMM_expressed_genes,]
my_genes <- row.names(subset(fData(HSMM_filtered),
          gene_short_name %in% c("CDK1", "MEF2C", "MYH3")))
cds_subset <- HSMM_filtered[my_genes,]
plot_genes_in_pseudotime(cds_subset, color_by = "Hours")

Alternative choices for ordering genes

Ordering based on genes that differ between clusters（上述提到的dpFeature算法,我用的是这种方法挑选的差异基因）

#我理解的意思就是通过算法使得细胞分成和实验差不多的结果（给出的例子就是同一时间收的细胞聚在一起，在这里我用的是seurat分群的结果作为依据），然后挑选出这些分群的基因，作为拟时间分析的基因。
#挑选出至少在5%的细胞中表达的基因
HSMM_myo=HSMM
HSMM_myo <- detectGenes(HSMM_myo, min_expr = 0.1)
fData(HSMM_myo)$use_for_ordering <-fData(HSMM_myo)$num_cells_expressed > 0.05 * ncol(HSMM_myo)

#然后就是进行PCA分析，根据下面的图选择合适的PC
plot_pc_variance_explained(HSMM_myo, return_all = F)

#然后重新降维，这里的num_dim我用的是seurat中选择的PC数
HSMM_myo <- reduceDimension(HSMM_myo,
                              max_components = 2,
                              norm_method = 'log',
                              num_dim = 3,
                              reduction_method = 'tSNE',
                              verbose = T)

#随后进行密度峰值聚类，鉴定cluster
#The densityPeak algorithm clusters cells based on each cell's local density (Ρ) and the nearest distance (Δ) of a cell to another cell with higher distance
#可以人为自己设定Ρ和Δ
#默认参数进行聚类
HSMM_myo <- clusterCells(HSMM_myo, verbose = F)
#查看聚类结果，主要是看一下用聚类分出来的结果和实验分群结果（seurat分群结果）是否一致，如果一致，那就可以进一步找这些cluster中的基因，作为分群的基因，不一致就要自己调参数
plot_cell_clusters(HSMM_myo, color_by = 'as.factor(Cluster)')
plot_cell_clusters(HSMM_myo, color_by = 'as.factor(Hours)') 

感觉用默认阈值聚类出来的效果不好，尝试自己设定阈值

#查看每一个细胞的Ρ和Δ，然后自己设定阈值，图片中的黑色点代表cluster的数量
plot_rho_delta(HSMM_myo, rho_threshold = 2, delta_threshold = 4)
#按照设定的阈值重新聚类
HSMM_myo <- clusterCells(HSMM_myo,
                 rho_threshold = 2,
                 delta_threshold = 10,
                 skip_rho_sigma = T,
                 verbose = F)
#查看聚类结果，两种结果图很一致，那么就找这些群的基因
plot_cell_clusters(HSMM_myo, color_by = 'as.factor(Cluster)')
plot_cell_clusters(HSMM_myo, color_by = 'as.factor(Hours)')

#这里的res.0.2是fData(HSMM_myo)的列名，可以根据自己的需要进行选择，我这里选择的是分辨率为0.2的seurat聚类结果。
#head(fData(HSMM_myo))

当我们觉得这个聚类make sense时，就可以进行下一步，找到区分这些cluster的基因

#找到差异基因
HSMM_expressed_genes <-  row.names(subset(fData(HSMM_myo),
                                          num_cells_expressed >= 10))
#这一步耗时很长
clustering_DEG_genes <-
    differentialGeneTest(HSMM_myo[HSMM_expressed_genes,],
          fullModelFormulaStr = '~Cluster',
          cores = 1)
#选取前1000个基因进行拟时间轨迹分析
#和之前的步骤一样，第一步确定合适的基因
HSMM_ordering_genes <-row.names(clustering_DEG_genes)[order(clustering_DEG_genes$qval)][1:1000]
#将基因添加到HSMM对象中
HSMM_myo <-setOrderingFilter(HSMM_myo,ordering_genes = HSMM_ordering_genes)
#第二步用DDRTree进行降维
HSMM_myo <-reduceDimension(HSMM_myo, method = 'DDRTree')
#第三步构建拟时间曲线
HSMM_myo <-orderCells(HSMM_myo)

HSMM_myo <-orderCells(HSMM_myo, root_state = GM_state(HSMM_myo))
#这一步出图
plot_cell_trajectory(HSMM_myo, color_by = "Hours")

小结：
构建拟时间轨迹一共分为三步，其中第一步选择较多，二三步基本一致：
第一步：确定合适的基因（这里仅针对前两种进行说明）；
· 简单地比较过程开始时和结束时收集的细胞，找出差异表达的基因；
· 基于不同的cluster找差异基因（官网推荐使用的）；
· 选择细胞间高度分散的基因(这里没有提到)；
第二步：用DDRTree算法降维
第三步：构建拟时间曲线