Monocle2包學(xué)習(xí)筆記（三）：Constructing Single Cell Trajectories

image

Monocle2使用反向圖嵌入（Reversed Graph Embedding）的機(jī)器學(xué)習(xí)算法，來(lái)對(duì)單細(xì)胞RNA-seq數(shù)據(jù)進(jìn)行單細(xì)胞發(fā)育軌跡推斷分析。Monocle2不是通過(guò)實(shí)驗(yàn)方法將細(xì)胞純化成不同的離散狀態(tài)，而是使用一種算法來(lái)學(xué)習(xí)每個(gè)細(xì)胞在動(dòng)態(tài)生物學(xué)過(guò)程中經(jīng)歷的基因表達(dá)變化的順序。一旦它了解了基因表達(dá)變化的整體“軌跡”，Monocle2就可以將每個(gè)細(xì)胞分配到發(fā)育軌跡中的適當(dāng)位置中。并且，我們還可以使用Monocle2的差異分析工具來(lái)查找在發(fā)育軌跡過(guò)程中受到調(diào)控的基因。如果這個(gè)發(fā)育過(guò)程中存在多個(gè)分支結(jié)果，Monocle2將重建一個(gè)“分支”軌跡。這些不同的分支與細(xì)胞的“決策”相對(duì)應(yīng)，Monocle2提供了強(qiáng)大的分析工具來(lái)識(shí)別受它們影響的基因，并參與這些基因的形成。

在整個(gè)生命的發(fā)育過(guò)程中，細(xì)胞會(huì)對(duì)不同的刺激做出相應(yīng)的響應(yīng)，從一種功能“狀態(tài)”轉(zhuǎn)變到另一種功能“狀態(tài)”。不同狀態(tài)下的細(xì)胞會(huì)表達(dá)不同的基因，產(chǎn)生不同的蛋白質(zhì)及其代謝產(chǎn)物，從而實(shí)現(xiàn)不同功能狀態(tài)下的工作。當(dāng)細(xì)胞在不同功能狀態(tài)之間轉(zhuǎn)變時(shí)，它們會(huì)經(jīng)歷一個(gè)轉(zhuǎn)錄重新配置的過(guò)程，其中一些基因的表達(dá)被沉默，而另一些基因則被激活。這些細(xì)胞的瞬時(shí)狀態(tài)通常很難被表征，因?yàn)樵诟€(wěn)定的端點(diǎn)狀態(tài)之間純化細(xì)胞可能是非常困難的或不可能的。單細(xì)胞RNA-Seq技術(shù)可以使我們?cè)诓恍枰兓?xì)胞的情況下，通過(guò)擬時(shí)序分析推斷出這些瞬時(shí)狀態(tài)下的細(xì)胞狀態(tài)。然而，要做到這一點(diǎn)，我們必須確定每個(gè)細(xì)胞在可能的狀態(tài)范圍內(nèi)的位置。

細(xì)胞排序的工作流程

Trajectory step 1: choose genes that define a cell's progress（選擇用于定義細(xì)胞進(jìn)程的基因）

首先，我們必須選擇用哪些基因去用于后續(xù)的細(xì)胞排序過(guò)程，這個(gè)過(guò)程稱為特征選擇（feature selection），篩選出的這些基因?qū)罄m(xù)細(xì)胞軌跡的形狀有著最重要的影響。因此，我們應(yīng)該盡可能的選擇那些最能反映細(xì)胞狀態(tài)的基因。
這里，我們需要選擇用哪些基因來(lái)定義細(xì)胞肌發(fā)生的進(jìn)程。我們最終希望能獲得一組基因集，這些基因的表達(dá)隨著研究的發(fā)育過(guò)程的進(jìn)展而增加（或減少）。篩選用于細(xì)胞排序的基因的一種有效方法是，將發(fā)育過(guò)程開(kāi)始時(shí)收集的細(xì)胞與結(jié)束時(shí)收集的細(xì)胞簡(jiǎn)單地進(jìn)行比較，并找到差異表達(dá)的基因。

# 使用differentialGeneTest函數(shù)進(jìn)行差異分析
diff_test_res <- differentialGeneTest(HSMM_myo[expressed_genes,],
              fullModelFormulaStr = "~Media")
# 提取顯著的差異表達(dá)基因
ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
# 細(xì)胞篩選過(guò)濾
HSMM_myo <- setOrderingFilter(HSMM_myo, ordering_genes)
# 篩選結(jié)果可視化
plot_ordering_genes(HSMM_myo)
image

上圖中，那些黑色的點(diǎn)對(duì)應(yīng)的基因，即為篩選出的用于后續(xù)對(duì)細(xì)胞進(jìn)行排序的基因。

Trajectory step 2: reduce data dimensionality（降低數(shù)據(jù)的維度）

接下來(lái)，我們使用reduceDimension函數(shù)對(duì)篩選出的基因進(jìn)行數(shù)據(jù)降維處理，

# 使用reduceDimension函數(shù)進(jìn)行數(shù)據(jù)降維，并指定method = 'DDRTree'參數(shù)使用DDRTree方法
HSMM_myo <- reduceDimension(HSMM_myo, max_components = 2,
    method = 'DDRTree')

Trajectory step 3: order cells along the trajectory（根據(jù)發(fā)育軌跡對(duì)細(xì)胞進(jìn)行排序）

最后，對(duì)于降維后的數(shù)據(jù)，我們將使用orderCells函數(shù)對(duì)細(xì)胞按發(fā)育軌跡進(jìn)行排序。

# 使用orderCells函數(shù)進(jìn)行細(xì)胞排序
HSMM_myo <- orderCells(HSMM_myo)
# 使用plot_cell_trajectory函數(shù)對(duì)排序的細(xì)胞進(jìn)行可視化
plot_cell_trajectory(HSMM_myo, color_by = "Hours")

image

從上圖可以看出，該細(xì)胞發(fā)育軌跡呈樹(shù)狀結(jié)構(gòu)。其中，在0時(shí)收集的細(xì)胞大多位于樹(shù)的頂端之一附近，而其他的則分布在兩個(gè)“分支”之間。Monocle2不知道先驗(yàn)樹(shù)的哪個(gè)軌跡作為“beginning”，因此我們經(jīng)常不得不使用root_state參數(shù)再次調(diào)用orderCells函數(shù)來(lái)指定起點(diǎn)。

plot_cell_trajectory(HSMM_myo, color_by = "State")
image

這里我們使用細(xì)胞所處的狀態(tài)對(duì)細(xì)胞進(jìn)行著色。

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)
  }
}
# 使用root_state參數(shù)指定起點(diǎn)
HSMM_myo <- orderCells(HSMM_myo, root_state = GM_state(HSMM_myo))
# 根據(jù)偽時(shí)間對(duì)細(xì)胞進(jìn)行著色
plot_cell_trajectory(HSMM_myo, color_by = "Pseudotime")
image
# 根據(jù)細(xì)胞的狀態(tài)進(jìn)行分面展示
plot_cell_trajectory(HSMM_myo, color_by = "State") +
    facet_wrap(~State, nrow = 1)
image

如果我們沒(méi)有時(shí)間序列數(shù)據(jù)，則可能需要根據(jù)已有的生物學(xué)知識(shí)，利用某些marker標(biāo)記基因的表達(dá)來(lái)設(shè)置root。例如，在該實(shí)驗(yàn)中，高度增殖的祖細(xì)胞群正在產(chǎn)生兩種類(lèi)型的有絲分裂后細(xì)胞。因此，根部應(yīng)具有表達(dá)高水平增殖marker基因的細(xì)胞。我們可以使用抖動(dòng)圖來(lái)可視化marker基因的表達(dá)，找出marker基因的哪個(gè)狀態(tài)對(duì)應(yīng)于快速增殖的細(xì)胞：

# 選擇marker基因
blast_genes <- row.names(subset(fData(HSMM_myo),
                         gene_short_name %in% c("CCNB2", "MYOD1", "MYOG")))
# 使用plot_genes_jitter函數(shù)可視化marker基因的表達(dá)
plot_genes_jitter(HSMM_myo[blast_genes,],
    grouping = "State",
    min_expr = 0.1)
image

為了確認(rèn)細(xì)胞排序的準(zhǔn)確性，我們可以選擇成肌過(guò)程的幾個(gè)marker基因，繪制它們的表達(dá)在偽時(shí)間過(guò)程中的變化。

HSMM_expressed_genes <-  row.names(subset(fData(HSMM_myo),
                                   num_cells_expressed >= 10))
HSMM_filtered <- HSMM_myo[HSMM_expressed_genes,]
# 篩選marker基因
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函數(shù)可視化marker基因的表達(dá)
plot_genes_in_pseudotime(cds_subset, color_by = "Hours")
image

Alternative choices for ordering genes

Ordering based on genes that differ between clusters

我們建議用戶使用通過(guò)無(wú)監(jiān)督程序“dpFeature”選擇的基因用于后續(xù)的細(xì)胞排序。
要使用dpFeature，我們首先篩選出至少在5％的所有細(xì)胞中都表達(dá)的基因。

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)

接下來(lái)，我們執(zhí)行PCA分析，以識(shí)別每個(gè)PC所能解釋的方差，可以通過(guò)繪制碎石圖來(lái)查看選擇需要的PC。

plot_pc_variance_explained(HSMM_myo, return_all = F)
image

然后，我們將那些主要的PC使用t-SNE降維的方法進(jìn)行降維處理，并將其投影到二維層面進(jìn)行可視化展示。

# 使用reduceDimension函數(shù)進(jìn)行數(shù)據(jù)降維
HSMM_myo <- reduceDimension(HSMM_myo,
                              max_components = 2,
                              norm_method = 'log',
                              num_dim = 3,
                              reduction_method = 'tSNE',
                              verbose = T)
# 使用clusterCells函數(shù)對(duì)細(xì)胞進(jìn)行聚類(lèi)分群
HSMM_myo <- clusterCells(HSMM_myo, verbose = F)
# 使用plot_cell_clusters函數(shù)對(duì)聚類(lèi)后的結(jié)果進(jìn)行可視化
plot_cell_clusters(HSMM_myo, color_by = 'as.factor(Cluster)')
plot_cell_clusters(HSMM_myo, color_by = 'as.factor(Hours)')
image

After we confirm the clustering makes sense, we can then perform differential gene expression test as a way to extract the genes that distinguish them.

# 使用differentialGeneTest函數(shù)進(jìn)行差異表達(dá)分析
clustering_DEG_genes <- differentialGeneTest(HSMM_myo[HSMM_expressed_genes,],
          fullModelFormulaStr = '~Cluster',
          cores = 1)

We will then select the top 1000 significant genes as the ordering genes.

# 選擇top 1000的基因
HSMM_ordering_genes <- row.names(clustering_DEG_genes)[order(clustering_DEG_genes$qval)][1:1000]

HSMM_myo <- setOrderingFilter(HSMM_myo, ordering_genes = HSMM_ordering_genes)
# 數(shù)據(jù)降維
HSMM_myo <- reduceDimension(HSMM_myo, method = 'DDRTree')
# 細(xì)胞排序
HSMM_myo <- orderCells(HSMM_myo)
# 設(shè)置root狀態(tài)
HSMM_myo <- orderCells(HSMM_myo, root_state = GM_state(HSMM_myo))
# 排序結(jié)果可視化
plot_cell_trajectory(HSMM_myo, color_by = "Hours")
image

Selecting genes with high dispersion across cells

變化很大的基因通常對(duì)于識(shí)別細(xì)胞亞群或沿著軌跡排序細(xì)胞很有幫助。在單細(xì)胞RNA-Seq中，基因的變異通常取決于其均值，因此我們必須謹(jǐn)慎選擇如何根據(jù)它們的變化來(lái)選擇基因。

disp_table <- dispersionTable(HSMM_myo)
ordering_genes <- subset(disp_table,
                  mean_expression >= 0.5 &
                  dispersion_empirical >= 1 * dispersion_fit)$gene_id

Ordering cells using known marker genes

無(wú)監(jiān)督聚類(lèi)對(duì)細(xì)胞進(jìn)行排序是可取的，因?yàn)樗苊饬嗽诜治鲋幸肫?。但是，無(wú)監(jiān)督的機(jī)器學(xué)習(xí)方法也會(huì)存在一定的問(wèn)題。例如，當(dāng)我們使用無(wú)監(jiān)督學(xué)習(xí)時(shí)，處于不同細(xì)胞周期的細(xì)胞會(huì)對(duì)發(fā)育軌跡的形狀產(chǎn)生重大的影響。但是，如果我們希望專注于生物學(xué)過(guò)程中與細(xì)胞周期無(wú)關(guān)的影響，該怎么辦？Monocle2的“半監(jiān)督”聚類(lèi)排序方法可以幫助我們實(shí)現(xiàn)相應(yīng)的處理。

以半監(jiān)督聚類(lèi)方式對(duì)細(xì)胞進(jìn)行排序非常簡(jiǎn)單。首先使用CellTypeHierchy系統(tǒng)定義標(biāo)記細(xì)胞進(jìn)展的基因，這與我們將其用于細(xì)胞類(lèi)型分類(lèi)的方式非常相似。然后，使用它來(lái)選擇與這些標(biāo)記共同變化的有序基因。最后，就像我們?cè)跓o(wú)監(jiān)督的排序中一樣，根據(jù)這些基因?qū)?xì)胞進(jìn)行排序。因此，無(wú)監(jiān)督排序和半監(jiān)督排序之間的唯一區(qū)別是我們使用哪些基因進(jìn)行排序。

# 選擇marker基因
CCNB2_id <- row.names(subset(fData(HSMM_myo), gene_short_name == "CCNB2"))
MYH3_id <- row.names(subset(fData(HSMM_myo), gene_short_name == "MYH3"))
# 初始化CellTypeHierarchy對(duì)象
cth <- newCellTypeHierarchy()
# 根據(jù)marker基因的表達(dá)標(biāo)記不同類(lèi)型的細(xì)胞
cth <- addCellType(cth,
                   "Cycling myoblast",
                   classify_func = function(x) { x[CCNB2_id,] >= 1 })

cth <- addCellType(cth,
                   "Myotube",
                   classify_func = function(x) { x[MYH3_id,] >= 1 })

cth <- addCellType(cth,
                   "Reserve cell",
                   classify_func = function(x) { x[MYH3_id,] == 0 & x[CCNB2_id,] == 0 })
# 使用classifyCells函數(shù)對(duì)細(xì)胞進(jìn)行分類(lèi)
HSMM_myo <- classifyCells(HSMM_myo, cth)

Now we select the set of genes that co-vary (in either direction) with these two "bellweather" genes:

marker_diff <- markerDiffTable(HSMM_myo[HSMM_expressed_genes,],
                       cth,
                       cores = 1)
#semisup_clustering_genes <- row.names(subset(marker_diff, qval < 0.05))
# 選擇top 1000的基因
semisup_clustering_genes <-
    row.names(marker_diff)[order(marker_diff$qval)][1:1000]

Using the top 1000 genes for ordering produces a trajectory that's highly similar to the one we obtained with unsupervised methods, but it's a little "cleaner".

HSMM_myo <- setOrderingFilter(HSMM_myo, semisup_clustering_genes)
#plot_ordering_genes(HSMM_myo)
# 數(shù)據(jù)降維
HSMM_myo <- reduceDimension(HSMM_myo, max_components = 2,
                            method = 'DDRTree', norm_method = 'log')
# 對(duì)細(xì)胞進(jìn)行排序
HSMM_myo <- orderCells(HSMM_myo)
HSMM_myo <- orderCells(HSMM_myo, root_state = GM_state(HSMM_myo))
# 細(xì)胞排序結(jié)果可視化
plot_cell_trajectory(HSMM_myo, color_by = "CellType") +
                     theme(legend.position = "right")
image

參考來(lái)源：http://cole-trapnell-lab.github.io/monocle-release/docs/

image
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書(shū)系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容