Monocle可以進(jìn)行下面三種類型的分析，可以與seurat無縫連接：

1 對細(xì)胞進(jìn)行聚類、分組和計算；
2.進(jìn)行擬時間分析；
3.差異表達(dá)分析；

2.Constructing Single Cell Trajectories

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

理論上希望找到一組能夠隨著研究過程的進(jìn)展而增加(或減少)表達(dá)的基因。同時盡可能少地使用正在研究的系統(tǒng)生物學(xué)的先驗知識。我們希望從數(shù)據(jù)中發(fā)現(xiàn)重要的排序基因，而不是依賴于文獻(xiàn)和教科書，因為這可能會在排序中引入偏見。將從一種簡單的方式開始，但是推薦使用另外一種更有效的方法“dpFeature”(后面會提到)。
首先，找到一組基因的有效方法之一是簡單地比較過程開始時和結(jié)束時收集的細(xì)胞（前提是你有這一組信息），并找到差異表達(dá)的基因。

#這一種情況適用于有一個明確的時間節(jié)點(diǎn)（這里就是換培養(yǎng)基）
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")

根據(jù)上面的擬時間軌跡，可以看出0時收集的細(xì)胞聚集在樹的一端，其他時期分布在數(shù)的兩條分支上。Monocle是不知道樹的哪邊是起始端，因此往往需要再次調(diào)用orderCells函數(shù)，使用root_state參數(shù)去指定起始端。

#下面就是進(jìn)行的再次調(diào)用orderCells函數(shù)的步驟
#首先按照樹的分段進(jìn)行顏色區(qū)分（state參數(shù)）
plot_cell_trajectory(HSMM_myo, color_by = "State")

可以看到，一共有5種狀態(tài)，就是有五個顏色

#下面這個函數(shù)用來識別哪一段樹枝包含時間0的細(xì)胞最多?這里不是太理解
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")

#如果你細(xì)胞數(shù)太多，可以分開看每個細(xì)胞亞群所處的位置
plot_cell_trajectory(HSMM_myo, color_by = "State") +
    facet_wrap(~State, nrow = 1)

如果沒有時間序列，就需要根據(jù)生物學(xué)知識找到root。就這個實(shí)驗而言，一個高度增殖的祖細(xì)胞群產(chǎn)生兩種有絲分裂的細(xì)胞。因此，root就應(yīng)該為高表達(dá)增殖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)

#為了確認(rèn)順序是正確的，我們可以選擇幾個肌源性進(jìn)展過程中的標(biāo)記
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算法,我用的是這種方法挑選的差異基因）

#我理解的意思就是通過算法使得細(xì)胞分成和實(shí)驗差不多的結(jié)果（給出的例子就是同一時間收的細(xì)胞聚在一起，在這里我用的是seurat分群的結(jié)果作為依據(jù)），然后挑選出這些分群的基因，作為擬時間分析的基因。
#挑選出至少在5%的細(xì)胞中表達(dá)的基因
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)

#然后就是進(jìn)行PCA分析，根據(jù)下面的圖選擇合適的PC
plot_pc_variance_explained(HSMM_myo, return_all = F)

#然后重新降維，這里的num_dim我用的是seurat中選擇的PC數(shù)
HSMM_myo <- reduceDimension(HSMM_myo,
                              max_components = 2,
                              norm_method = 'log',
                              num_dim = 3,
                              reduction_method = 'tSNE',
                              verbose = T)

#隨后進(jìn)行密度峰值聚類，鑒定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
#可以人為自己設(shè)定Ρ和Δ
#默認(rèn)參數(shù)進(jìn)行聚類
HSMM_myo <- clusterCells(HSMM_myo, verbose = F)
#查看聚類結(jié)果，主要是看一下用聚類分出來的結(jié)果和實(shí)驗分群結(jié)果（seurat分群結(jié)果）是否一致，如果一致，那就可以進(jìn)一步找這些cluster中的基因，作為分群的基因，不一致就要自己調(diào)參數(shù)
plot_cell_clusters(HSMM_myo, color_by = 'as.factor(Cluster)')
plot_cell_clusters(HSMM_myo, color_by = 'as.factor(Hours)') 

感覺用默認(rèn)閾值聚類出來的效果不好，嘗試自己設(shè)定閾值

#查看每一個細(xì)胞的Ρ和Δ，然后自己設(shè)定閾值，圖片中的黑色點(diǎn)代表cluster的數(shù)量
plot_rho_delta(HSMM_myo, rho_threshold = 2, delta_threshold = 4)
#按照設(shè)定的閾值重新聚類
HSMM_myo <- clusterCells(HSMM_myo,
                 rho_threshold = 2,
                 delta_threshold = 10,
                 skip_rho_sigma = T,
                 verbose = F)
#查看聚類結(jié)果，兩種結(jié)果圖很一致，那么就找這些群的基因
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)的列名，可以根據(jù)自己的需要進(jìn)行選擇，我這里選擇的是分辨率為0.2的seurat聚類結(jié)果。
#head(fData(HSMM_myo))

當(dāng)我們覺得這個聚類make sense時，就可以進(jìn)行下一步，找到區(qū)分這些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個基因進(jìn)行擬時間軌跡分析
#和之前的步驟一樣，第一步確定合適的基因
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進(jìn)行降維
HSMM_myo <-reduceDimension(HSMM_myo, method = 'DDRTree')
#第三步構(gòu)建擬時間曲線
HSMM_myo <-orderCells(HSMM_myo)

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

小結(jié)：
構(gòu)建擬時間軌跡一共分為三步，其中第一步選擇較多，二三步基本一致：
第一步：確定合適的基因（這里僅針對前兩種進(jìn)行說明）；
· 簡單地比較過程開始時和結(jié)束時收集的細(xì)胞，找出差異表達(dá)的基因；
· 基于不同的cluster找差異基因（官網(wǎng)推薦使用的）；
· 選擇細(xì)胞間高度分散的基因(這里沒有提到)；
第二步：用DDRTree算法降維
第三步：構(gòu)建擬時間曲線