Monocle2 | 單細胞測序的擬時序分析（細胞軌跡分析）

Amazing427 2022-04-24

展開全文

背景介紹

偽時間是衡量單個細胞在細胞分化等過程中取得了多大進展的指標。在許多生物學過程中，細胞并不是完全同步的。在細胞分化等過程的單細胞表達研究中，捕獲的細胞在分化方面可能分布廣泛。也就是說，在同一時間捕獲的細胞群中，有些細胞可能已經很長時間了，而有些細胞甚至還沒有開始這個過程。當您想要了解在細胞從一種狀態轉換到另一種狀態時所發生的調節更改的順序時，這種異步性會產生主要問題。跟蹤同時捕獲的細胞間的表達可以產生對基因動力學一個大致的認識，該基因表達的明顯變異性將非常高。Monocle根據每個cell在學習軌跡上的進展對其進行排序，從而緩解了由于異步而產生的問題。Monocle不是跟蹤表達式隨時間變化的函數，而是跟蹤沿軌跡變化的函數，我們稱之為偽時間。偽時間是一個抽象的分化單位：它只是一個cell到軌跡起點的距離，沿著最短路徑測量。軌跡的總長度是由細胞從起始狀態移動到結束狀態所經歷的總轉錄變化量來定義的。

還是老規矩，跟著官網學習Monocle2的使用方法。官網如下：

http://cole-trapnell-lab.github.io/monocle-release/docs/

一、Abstract

單細胞基因表達研究使人們能夠在復雜的生物過程和高度異質性的細胞群中描述轉錄調控。這些研究有助于發現識別特定亞型細胞的基因，或標記生物過程中的中間狀態，以及在兩種不同的細胞命運之間過渡態的基因。在許多單細胞研究中，單個細胞是以不同步的方式執行基因表達過程。實際上，每個細胞都是正在研究的轉錄過程中的一個瞬間。Monocle包是分析單細胞測序的工具。Monocle引入了在偽時間(擬時間)內對單個細胞排序的策略，利用單個細胞的非同步進程，將它們置于與細胞分化等生物學過程相對應的軌跡上。****Monocle利用先進的機器學習技術(反向圖嵌入)從單細胞數據中學習顯式的主圖來對細胞進行排序，這可以強大而準確地解決復雜的生物過程。Monocle也可以進行聚類(即使用t-SNE和密度峰值聚類)。Monocle也可以進行差異基因表達測試，使人們能夠識別在不同狀態下差異表達的基因，沿著生物過程以及不同的細胞命運時基因表達的變化。Monocle是專為單細胞RNA-Seq研究設計的，但也可以用于其他分析。

二、Introduction

Monocle 2包括新的和改進的算法，用于對細胞進行分類和計數，執行細胞亞群之間的差異表達分析，以及重建細胞軌跡。Monocle 2也被重新設計，可以很好地完成包含數萬個或更多細胞的大型單細胞RNA-Seq數據分析。

Monocle perform three main types of analysis:

Clustering, classifying, and counting cells. Single-cell RNA-Seq experiments allow you to discover new (and possibly rare) subtypes of cells. Monocle helps you identify them.
Constructing single-cell trajectories. In development, disease, and throughout life, cells transition from one state to another. Monocle helps you discover these transitions.
Differential expression analysis. Characterizing new cell types and states begins with comparing them to other, better understood cells. Monocle includes a sophisticated but easy to use system for differential expression.

首先，Monocle 2使用一種簡單的、無偏的和高度可擴展的統計程序來選擇具有軌跡進展特征的基因。然后，它采用了一類流形學習算法，旨在在高維單細胞RNA-seq數據中嵌入一個主圖。以前的方法是通過啟發式分析細胞之間的成對距離來推斷分支結構，而Monocle 2可以使用這張圖來直接識別發育的命運決定。我們已經通過廣泛的基準測試證明，Monocle 2優于其他工具，如Wishbone，而不需要用戶指定軌跡的結構。

三、Monocle2算法

Monocle 2的算法是再2017年發表在nature methods 上，文章鏈接如下

https://www.nature.com/articles/nmeth.4402

文章中的核心理論為：每個細胞都可以表示為高維空間中的一個點，在高維空間中，每個維對應著一個有序基因的表達水平。高維數據首先通過幾種降維方法，如PCA(默認)、擴散映射等，投射到低維空間。Monocle 2然后在自動選擇的一組數據質心上構造一棵生成樹。然后，該算法將細胞移動到它們最近的樹的頂點，更新頂點的位置以適應細胞，學習新的生成樹，并迭代地繼續這個過程，直到樹和細胞的位置已經收斂。在這個過程中，Monocle 2保持了高維空間和低維空間之間的可逆映射，從而既學習了軌跡，又降低了數據的維數。****一旦Monocle 2學會了樹，用戶就會選擇一個tip作為根。計算每個單元的偽時間作為其沿樹到根的測地線距離，并根據主圖自動分配其分枝。因為monocle2學習樹結構，與其他方法相比，分支結構自動出現。當它更新細胞位置并細化樹時，monocle2簡化了軌跡的結構，修剪了小的分支，這樣最終的軌跡只保留了描述細胞狀態顯著差異的分支。

圖片

四、Installing Monocle2

Monocle需要在R語言環境下運行。官網給出的都是R-3.4比較老的版本，現在都是R-3.6或者R-4.0以上版本。我們可以直接在Bioconductor中安裝Monocle 2。

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("monocle")

查看系統中安裝的這個包的版本文檔

browseVignettes("monocle")

圖片

當然也可以在GitHub上安裝最新版的Monocle2。

install.packages("devtools")
devtools::install_github("cole-trapnell-lab/monocle-release@develop")

有時我們添加的特性需要你安裝某些額外的軟件包。當您嘗試上面的命令時，可能會看到錯誤。您可以通過輸入(例如)來在錯誤消息中安裝包。

BiocManager::install(c("DDRTree", "pheatmap"))

最后檢查一下是否安裝成功

library(monocle)

四、創建CellDataSet

方法一: 將Seurat object中數據提取來創建

data <- as(as.matrix(scRNA.Osteoclastic@assays$RNA@counts), 'sparseMatrix')
pd <- new('AnnotatedDataFrame', data = scRNA.Osteoclastic@meta.data)
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
fd <- new('AnnotatedDataFrame', data = fData)
monocle_cds <- newCellDataSet(data,
                              phenoData = pd,
                              featureData = fd,
                              lowerDetectionLimit = 0.5,
                              expressionFamily = negbinomial.size())

FPKM/TPM值通常是對數正態分布的，而UMIs或讀計數使用負二項更好地建模。要處理計數數據，請將負二項分布指定為newCellDataSet的expressionFamily參數：

稀疏矩陣用negbinomial.size()，

FPKM值用tobit()，

logFPKM值用gaussianff()

圖片

方法二: 直接讀取表達矩陣來創建

library(data.table)
##讀取數據
data <- fread("fpkm.txt",data.table = F,header = T)
pd <-  fread("metadata.txt",data.table = F,header = T)
fd <- fread("gene_annotations.txt",data.table = F,header = T)
##創建
pd <- new("AnnotatedDataFrame", data = pd)
fd <- new("AnnotatedDataFrame", data = fd)
HSMM <- newCellDataSet(as.matrix(data),
                       phenoData = pd, featureData = fd,
                       expressionFamily = tobit())
###如果數據量大，建議轉化為稀疏矩陣
HSMM <- newCellDataSet(as(as.matrix(data), "sparseMatrix"),
                       phenoData = pd, 
                       featureData = fd,
                       expressionFamily = tobit())

方法三: 將Seurat對象直接轉化為CellDataSet對象

importCDS(scRNA.Osteoclastic)

如果我們想要從Seurat對象或SCESet中導入所有的插槽，我們可以設置參數'import_all'為TRUE。#(默認為FALSE或只保留最小數據集)

五、估計size factor和離散度

size facotr幫助我們標準化細胞之間的mRNA的差異。

離散度值可以幫助我們進行后續的差異分析。


monocle_cds <- estimateSizeFactors(monocle_cds)
monocle_cds <- estimateDispersions(monocle_cds)

六、過濾低質量的細胞

大多數單細胞工作流程至少會包含一些由死細胞或空孔組成的庫。同樣重要的是要刪除doublets:由兩個或多個細胞意外生成的庫。這些細胞可以破壞下游步驟，如偽時間排序或聚類。要知道一個特定的基因有多少個表達，或者一個給定的細胞有多少個基因表達，通常是很方便的。Monocle提供了一個簡單的函數來計算這些統計數據。

monocle_cds <- detectGenes(monocle_cds, min_expr = 0.1)
print(head(fData(monocle_cds)))

圖片

七、細胞分類

Monocle官網教程提供了4個分類方法：

**Classifying cells by type **
Clustering cells without marker genes
Clustering cells using marker genes
Imputing cell type

1、Classifying cells by type

Monocle提供了一個簡單的系統，根據您選擇的標記基因的表達來標記細胞。您只需提供一組函數，Monocle就可以使用這些函數對每個單元格進行注釋。例如，您可以為幾種單元格類型中的每種類型提供一個函數。這些函數接受每個單元格的表達式數據作為輸入，并返回TRUE以告訴Monocle某個單元格滿足函數定義的條件。所以你可以有一個功能，對表達成肌細胞特異性基因的細胞來說是真的，另一個功能對成纖維細胞特異性基因來說是真的，等等。


MYF5_id <- row.names(subset(fData(HSMM), gene_short_name == "MYF5"))
ANPEP_id <- row.names(subset(fData(HSMM),
    gene_short_name == "ANPEP"))

cth <- newCellTypeHierarchy()
cth <- addCellType(cth, "Myoblast", classify_func =
    function(x) { x[MYF5_id,] >= 1 })
cth <- addCellType(cth, "Fibroblast", classify_func = function(x)
{ x[MYF5_id,] < 1 & x[ANPEP_id,] > 1 })
HSMM <- classifyCells(HSMM, cth, 0.1)
table(pData(HSMM)$CellType)
###畫圖
pie <- ggplot(pData(HSMM),
aes(x = factor(1), fill = factor(CellType))) + geom_bar(width = 1)
pie + coord_polar(theta = "y") +
theme(axis.title.x = element_blank(), axis.title.y = element_blank())

2、Clustering cells without marker genes

Monocle提供了一種算法，可以用來計算“未知”單元格的類型。該算法在函數clusterCells中實現，根據全局表達式輪廓將單元格分組在一起。這樣的話，如果你的細胞表達了許多特定于成肌細胞的基因，但恰好缺少MYF5，我們仍然可以識別它為成肌細胞。clusterCells可以在無監督的方式下使用，也可以在半監督的“”模式下使用，這允許使用一些專家知識來輔助算法。我們先來看看無監督模式。

第一步是決定哪些基因用于細胞聚類。我們可以使用所有的基因，但是我們會包含很多基因，這些基因的表達水平不足以提供有意義的信號。包括它們只會給系統增加噪音。我們可以根據平均表達水平篩選基因，我們還可以選擇在細胞間異?？勺兊幕?。這些基因往往對細胞狀態有很高的信息量。

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

圖片

setOrderingFilter函數標記了在隨后對clusterCells的調用中用于聚類的基因，盡管我們可以根據需要提供其他基因列表。plot_ordering_genes函數顯示了基因表達的變異性（離散性）如何依賴于細胞間的平均表達。紅線表示基于這種關系的單片體對色散的預期。我們標記用于聚類的基因顯示為黑點，而其他基因顯示為灰點。

plot_pc_variance_explained(HSMM, return_all = F)

圖片

Monocle使用t-SNE來給細胞聚類可視化

HSMM <- reduceDimension(HSMM, max_components = 2, num_dim = 6,
                reduction_method = 'tSNE', verbose = T)
HSMM <- clusterCells(HSMM, num_clusters = 2)
plot_cell_clusters(HSMM, 1, 2, color = "CellType",
    markers = c("MYF5", "ANPEP"))

圖片

指定10個cluster

HSMM <- clusterCells(HSMM, num_clusters = 10)
plot_cell_clusters(HSMM)

圖片

3、Clustering cells using marker genes
首先，我們將選擇一組不同的基因用于細胞聚類。在我們挑選高表達和高度可變的基因之前。現在，我們將挑選與標記共變異的基因。從某種意義上說，我們將建立一個大的基因列表來作為標記，這樣即使一個細胞沒有MYF5，它也可以被其他基因識別為成肌細胞。

marker_diff <- markerDiffTable(HSMM[expressed_genes,],
            cth,
            residualModelFormulaStr = "~Media + num_genes_expressed",
            cores = 1)

函數markerDiffTable接受CellDataSet和CellTypeHierarchy，并根據提供的函數將所有單元格分類為類型。然后在識別不同基因類型之間差異表達的基因之前，它去除了所有“未知”和“不明確”的功能。同樣，您可以提供要從該測試中排除的影響的剩余模型。然后函數返回一個測試結果的數據幀，您可以使用它來選擇要用于聚類的基因。通常情況下，最好挑選最適合每種細胞類型的前10或20個基因。這確保了聚類基因不受一種細胞類型標記的支配。如果可能的話，通常需要為每種類型創建一個平衡的標記面板。Monocle提供了一個方便的功能，可以根據基因在每種類型中的表達受限程度來對基因進行排序。

candidate_clustering_genes <-
    row.names(subset(marker_diff, qval < 0.01))
marker_spec <-
  calculateMarkerSpecificity(HSMM[candidate_clustering_genes,], cth)
head(selectTopMarkers(marker_spec, 3))

semisup_clustering_genes <- unique(selectTopMarkers(marker_spec, 500)$gene_id)
HSMM <- setOrderingFilter(HSMM, semisup_clustering_genes)
plot_ordering_genes(HSMM)

plot_pc_variance_explained(HSMM, return_all = F)

HSMM <- reduceDimension(HSMM, max_components = 2, num_dim = 3,
  norm_method = 'log',
  reduction_method = 'tSNE',
  residualModelFormulaStr = "~Media + num_genes_expressed",
  verbose = T)
HSMM <- clusterCells(HSMM, num_clusters = 2)
plot_cell_clusters(HSMM, 1, 2, color = "CellType")

4、Imputing cell type
注意，我們已經減少了成肌細胞群中“污染”成纖維細胞的數量，反之亦然。但是“未知”細胞呢？如果為clusterCells提供CellTypeHierarcy，Monocle將使用它對整個簇進行分類，而不僅僅是單個細胞。基本上，clusterCells的工作原理與以前完全相同，只是在構建集群之后，它會統計每個集群中每種細胞類型的頻率。當一個簇由超過一定百分比（在本例中為10%）的特定類型組成時，簇中的所有單元格都被設置為該類型。如果一個集群由多個單元類型組成，則整個集群都被標記為“不明確”。如果沒有高于閾值的單元類型，則集群被標記為“未知”。因此，Monocle可以幫助您在缺少標記數據的情況下估算每個細胞的類型。

HSMM <- clusterCells(HSMM,
              num_clusters = 2,
              frequency_thresh = 0.1,
              cell_type_hierarchy = cth)
plot_cell_clusters(HSMM, 1, 2, color = "CellType",
    markers = c("MYF5", "ANPEP"))

pie <- ggplot(pData(HSMM), aes(x = factor(1), fill =
    factor(CellType))) +
geom_bar(width = 1)
pie + coord_polar(theta = "y") +
theme(axis.title.x = element_blank(), axis.title.y = element_blank())

八、構建軌跡

在發育過程中，為了對刺激做出反應，在整個生命周期中，細胞從一種功能“狀態”過渡到另一種功能“狀態”。處于不同狀態的細胞表達不同的基因，產生蛋白質和代謝物的動態重復，從而完成它們的工作。當細胞在不同的狀態間移動時，會經歷一個轉錄重組的過程，一些基因被沉默，而另一些則被激活。這些瞬態通常很難描述，因為在更穩定的端點狀態之間凈化細胞可能很困難或不可能。單細胞RNA-Seq可以讓你看到這些狀態而不需要純化。然而，要做到這一點，我們必須確定每個細胞的可能狀態范圍。
Monocle介紹了利用RNA-Seq進行單細胞軌跡分析的策略。Monocle不是通過實驗將細胞凈化成離散狀態，而是使用一種算法來學習每個細胞作為動態生物過程的一部分必須經歷的基因表達變化序列。一旦掌握了基因表達變化的整體“軌跡”，Monocle就可以將每個細胞放置在軌跡中的適當位置。然后，你可以使用Monocle的差異分析工具來尋找在軌跡過程中受調控的基因，如“尋找隨假時間變化的基因”一節所述。如果這個過程有多個結果，Monocle將重建一個“分支”軌跡。這些分支與細胞的“決定”相對應，Monocle提供了強有力的工具來識別受它們影響并參與制造它們的基因。在“分析單細胞軌跡中的分支”部分中，可以看到如何分析分支。Monocle依靠一種稱為反向圖嵌入的機器學習技術來構造單細胞軌跡.

The ordering workflow

Step 1: choosing genes that define progress
Step 2: reducing the dimensionality of the data
Step 3: ordering the cells in pseudotime

Step 1: 選擇定義過程的基因

推斷單細胞軌跡是一個機器學習問題。第一步是選擇Monocle將用作機器學習方法輸入的基因。這叫做特征選擇，它對軌跡的形狀有很大的影響。在單細胞RNA-Seq中，低水平表達的基因通常非常嘈雜，但有些基因包含有關細胞狀態的重要信息。單核細胞通過檢測這些基因在細胞群體中的表達模式來對細胞進行排序。Monocle尋找以“有趣的”（即不只是嘈雜的）方式變化的基因，并用它們來構造數據。Monocle為您提供了多種工具來選擇基因，這些基因將產生一個健壯、準確和具有生物學意義的軌跡。你可以使用這些工具來執行一個完全“無監督”的分析，在這個分析中，Monocle不知道你認為哪個基因是重要的?；蛘?，你可以利用marker gene知識來定義生物學進展，從而形成Monocle的軌跡。我們認為這種模式是“半監督”的分析。

Monocle官網教程提供了4個選擇方法：

選擇發育差異表達基因
選擇clusters差異表達基因
選擇離散程度高的基因
自定義發育marker基因

前三種都是無監督分析方法，細胞發育軌跡生成完全不受人工干預；最后一種是半監督分析方法，可以使用先驗知識輔助分析。

##使用clusters差異表達基因
deg.cluster <- FindAllMarkers(scRNA.Osteoclastic)
diff.genes <- subset(deg.cluster,p_val_adj<0.05)$gene
HSMM <- setOrderingFilter(HSMM, diff.genes)
plot_ordering_genes(HSMM)
##使用seurat選擇的高變基因
var.seurat <- VariableFeatures(scRNA.Osteoclastic)
HSMM <- setOrderingFilter(HSMM, var.genes)
plot_ordering_genes(HSMM)
##使用monocle選擇的高變基因
disp_table <- dispersionTable(HSMM)
disp.genes <- subset(disp_table, mean_expression >= 0.1 & dispersion_empirical >= 1 * dispersion_fit)$gene_id
HSMM <- setOrderingFilter(HSMM, disp.genes)
plot_ordering_genes(HSMM)

圖片

理想情況下，我們希望盡可能少地使用正在研究的系統生物學的先驗知識。我們希望從數據中發現重要的排序基因，而不是依賴于文獻和教科書，因為這可能會在排序中引入偏見。我們將從一種更簡單的方法開始，但是我們通常推薦一種更復雜的方法，稱為“dpFeature”。


diff_test_res <- differentialGeneTest(HSMM_myo[expressed_genes,],
              fullModelFormulaStr = "~Media")
ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
HSMM_myo <- setOrderingFilter(HSMM_myo, ordering_genes)
plot_ordering_genes(HSMM_myo)

Step 2: 降維

下一步，我們將把空間縮小到一個二維空間，當Monocle對細胞進行排序時，我們將能夠很容易地可視化和解釋這些空間。

HSMM <- reduceDimension(HSMM, max_components = 2,
    method = 'DDRTree')

Step 3: 按照軌跡排序細胞

HSMM <- orderCells(HSMM)
plot_cell_trajectory(HSMM, color_by = "seurat_clusters")

圖片

plot_cell_trajectory(HSMM, color_by = "State")

圖片

plot_cell_trajectory(HSMM, color_by = "Pseudotime")

圖片

如果你的樹上有一堆狀態，那就很難弄清楚每個狀態都落在樹上的什么地方了。有時，可以方便地“分面”軌跡圖，以便更容易地查看每個狀態的位置.

plot_cell_trajectory(HSMM, color_by = "State") +
  facet_wrap(~State, nrow = 1)

圖片

如果你沒有時間序列，你可能需要利用你對系統的生物學知識，根據特定標記基因的表達位置來設置根。例如，在這個實驗中，一個高度增殖的祖細胞群正在產生兩種類型的有絲分裂后細胞。因此，根部應該有表達高水平增殖標記的細胞。我們可以使用抖動圖來判斷哪種狀態對應于快速擴散：

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

圖片

HSMM_expressed_genes <-  row.names(subset(fData(HSMM),
                                          num_cells_expressed >= 10))
HSMM_filtered <- HSMM[HSMM_expressed_genes,]
my_genes <- row.names(subset(fData(HSMM_filtered),
                             gene_short_name %in% c("YWHAB", "GAPDH", "TNNC1")))
cds_subset <- HSMM_filtered[my_genes,]
plot_genes_in_pseudotime(cds_subset, color_by = "seurat_clusters")

圖片

plot_genes_in_pseudotime(cds_subset, color_by =  "State")

圖片


genes <- c("TNNT2", "TNNC1", "CDK1")
p1 <- plot_genes_jitter(HSMM[genes,], grouping = "State", color_by = "State")
p2 <- plot_genes_violin(HSMM[genes,], grouping = "State", color_by = "State")
p3 <- plot_genes_in_pseudotime(HSMM[genes,], color_by = "State")
plotc <- p1|p2|p3

圖片

九、差異分析

官方給出的差異分析有三大方法，我們重點關注第三個：****根據偽時間功能尋找差異基因

1、Basic Differential Analysis
2、Finding Genes that Distinguish Cell Type or State
3、Finding Genes that Change as a Function of Pseudotime

Monocle的主要工作是通過生物過程（如細胞分化）將細胞按順序排列，而不知道要提前查看哪些基因。一旦這樣做了，你就可以分析細胞，找到隨著細胞進展而變化的基因。例如，你可以發現隨著細胞“成熟”而顯著上調的基因。讓我們來看看一組對肌肉生成很重要的基因：

to_be_tested <- row.names(subset(fData(HSMM),
                                 gene_short_name %in% c("MYH3", "MEF2C", "CCNB2", "TNNT1")))
cds_subset <- HSMM[to_be_tested,]

diff_test_res <- differentialGeneTest(cds_subset,
                                      fullModelFormulaStr = "~sm.ns(Pseudotime)")

diff_test_res[,c("gene_short_name", "pval", "qval")]
plot_genes_in_pseudotime(cds_subset, color_by ="State")

這個sm.ns函數指出Monocle應該通過表達式值擬合自然樣條曲線，以幫助它將表達式的變化描述為進程的函數。讓我們加入基因注釋，這樣就很容易看出哪些基因是重要的。

圖片

十、根據偽時間表達pattern聚類基因

** 研究時間序列基因表達研究時出現的一個常見問題是：“哪些基因遵循相似的動力學趨勢”？Monocle可以通過對具有相似趨勢的基因進行分組來幫助你回答這個問題，所以你可以分析這些基因組，看看它們有什么共同點。Monocle提供了一種方便的方法來可視化所有假時間依賴性基因。函數plot_pseudotime_heatmap采用CellDataSet對象（通常只包含重要基因的子集），并生成平滑的表達曲線，非常類似于plot_genes_in_pseudotime。然后，它對這些基因進行聚類，并使用pheatmap軟件包進行繪圖。這讓你可以想象出基因模塊在不同的時間內共同變化。**

marker_genes <- row.names(subset(fData(HSMM),
                                 gene_short_name %in% c("MEF2C", "MEF2D", "MYF5",
                                                        "ANPEP", "PDGFRA","MYOG",
                                                        "TPM1",  "TPM2",  "MYH2",
                                                        "MYH3",  "NCAM1", "TNNT1",
                                                        "TNNT2", "TNNC1", "CDK1",
                                                        "CDK2",  "CCNB1", "CCNB2",
                                                        "CCND1", "CCNA1", "ID1")))

diff_test_res <- differentialGeneTest(HSMM[marker_genes,],
                                      fullModelFormulaStr = "~sm.ns(Pseudotime)")
sig_gene_names <- row.names(subset(diff_test_res, qval < 0.1))
plot_pseudotime_heatmap(HSMM[sig_gene_names,],
                        num_clusters = 6,
                        cores = 1,
                        show_rownames = T)

圖片

十一、單細胞軌跡的“分支”分析

通常，單細胞軌跡包括分支。這些分支的產生是因為細胞執行不同的基因表達程序。在發育過程中，當細胞做出命運選擇時，分支出現在軌跡中：一個發育譜系沿著一條路徑前進，而另一個譜系產生第二條路徑。Monocle包含分析這些分支事件的廣泛功能。

以Steve Quake's的實驗室進行的實驗為例，Barbara Treutlein團隊從發育中的小鼠肺中捕獲了細胞。他們在發育早期捕獲細胞，后來當肺中同時含有兩種主要類型的上皮細胞（AT1和AT2）時，細胞就要決定變成AT1或AT2。Monocle可以將這個過程重構為一個分支軌跡，允許您對決策點進行非常詳細的分析。下圖顯示了使用部分數據重建的Monocle軌跡。有一個分支，標記為“1”。當細胞從樹的左上角經過樹枝的早期發育階段時，哪些基因會發生變化？哪些基因在分枝間有差異表達？為了回答這個問題，Monocle為您提供了一個特殊的統計測試：分支表達式分析建模，或BEAM。BEAM(Branched expression analysis modeling)是一種統計方法，用于尋找以依賴于分支的方式調控的基因。

plot_cell_trajectory(HSMM, color_by = "State")

圖片

BEAM進行統計分析

BEAM_res <- BEAM(HSMM, branch_point = 1, cores = 8)
BEAM_res <- BEAM_res[order(BEAM_res$qval),]
BEAM_res <- BEAM_res[,c("gene_short_name", "pval", "qval")]
plot_genes_branched_heatmap(HSMM[row.names(subset(BEAM_res,
                                                  qval < 1e-4)),],
                            branch_point = 1,
                            num_clusters = 4,
                            cores = 1,
                            use_gene_short_name = T,
                            show_rownames = T)

該熱圖顯示的是同一時間點兩個譜系的變化。熱圖的列是偽時間的點，行是基因。從熱圖中間往右讀，是偽時間的一個譜系；往左是另一個譜系。基因是被按照等級聚類的，所以你看到的基因表達模式和譜系的表達模式是非常相似的。

圖片

We can plot a couple of these genes, such as Pdpn and Sftpb (both known markers of fate in this system), using the plot_genes_branched_pseudotime function, which works a lot like the plot_genes_in_pseudotime function, except it shows two kinetic trends, one for each lineage, instead of one.

genes <- row.names(subset(fData(HSMM),
                          gene_short_name %in% c( "MEF2C", "CCNB2", "TNNT1")))

plot_genes_branched_pseudotime(HSMM[genes,],
                               branch_point = 1,
                               color_by = "State",
                               ncol = 1)