Monocle擬時間分析

Monocle可以進行下面三種類型的分析沪袭，可以與seurat無縫連接：

1 對細胞進行聚類叫倍、分組和計算枣抱；
2.進行擬時間分析；
3.差異表達分析蚊惯；

2.Constructing Single Cell Trajectories

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

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

#這一種情況適用于有一個明確的時間節(jié)點（這里就是換培養(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時收集的細胞聚集在樹的一端寓调，其他時期分布在數(shù)的兩條分支上。Monocle是不知道樹的哪邊是起始端锄码，因此往往需要再次調(diào)用orderCells函數(shù)夺英，使用root_state參數(shù)去指定起始端。

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

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

#下面這個函數(shù)用來識別哪一段樹枝包含時間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")

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

如果沒有時間序列重窟，就需要根據(jù)生物學(xué)知識找到root载萌。就這個實驗而言，一個高度增殖的祖細胞群產(chǎn)生兩種有絲分裂的細胞。因此扭仁，root就應(yīng)該為高表達增殖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算法,我用的是這種方法挑選的差異基因）

#我理解的意思就是通過算法使得細胞分成和實驗差不多的結(jié)果（給出的例子就是同一時間收的細胞聚在一起乖坠，在這里我用的是seurat分群的結(jié)果作為依據(jù)）搀突，然后挑選出這些分群的基因，作為擬時間分析的基因瓤帚。
#挑選出至少在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分析描姚，根據(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)

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

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

#查看每一個細胞的Ρ和Δ，然后自己設(shè)定閾值冰蘑，圖片中的黑色點代表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ù)自己的需要進行選擇武氓，我這里選擇的是分辨率為0.2的seurat聚類結(jié)果。
#head(fData(HSMM_myo))

當我們覺得這個聚類make sense時仇箱，就可以進行下一步县恕，找到區(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個基因進行擬時間軌跡分析
#和之前的步驟一樣，第一步確定合適的基因
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')
#第三步構(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)建擬時間軌跡一共分為三步剂桥，其中第一步選擇較多忠烛，二三步基本一致：
第一步：確定合適的基因（這里僅針對前兩種進行說明）；
· 簡單地比較過程開始時和結(jié)束時收集的細胞权逗，找出差異表達的基因美尸；
· 基于不同的cluster找差異基因（官網(wǎng)推薦使用的）；
· 選擇細胞間高度分散的基因(這里沒有提到)斟薇；
第二步：用DDRTree算法降維
第三步：構(gòu)建擬時間曲線