背景介紹
偽時間是衡量單個細胞在細胞分化等過程中取得了多大進展的指標。在許多生物學過程中，細胞并不是完全同步的。在細胞分化等過程的單細胞表達研究中，捕獲的細胞在分化方面可能分布廣泛。也就是說，在同一時間捕獲的細胞群中，有些細胞可能已經(jīng)很長時間了，而有些細胞甚至還沒有開始這個過程。當您想要了解在細胞從一種狀態(tài)轉(zhuǎn)換到另一種狀態(tài)時所發(fā)生的調(diào)節(jié)更改的順序時，這種異步性會產(chǎn)生主要問題。跟蹤同時捕獲的細胞間的表達可以產(chǎn)生對基因動力學一個大致的認識，該基因表達的明顯變異性將非常高。Monocle根據(jù)每個cell在學習軌跡上的進展對其進行排序，從而緩解了由于異步而產(chǎn)生的問題。Monocle不是跟蹤表達式隨時間變化的函數(shù)，而是跟蹤沿軌跡變化的函數(shù)，我們稱之為偽時間。偽時間是一個抽象的分化單位：它只是一個cell到軌跡起點的距離，沿著最短路徑測量。軌跡的總長度是由細胞從起始狀態(tài)移動到結(jié)束狀態(tài)所經(jīng)歷的總轉(zhuǎn)錄變化量來定義的。
還是老規(guī)矩，跟著官網(wǎng)學習Monocle2的使用方法。官網(wǎng)如下：
http://cole-trapnell-lab.github.io/monocle-release/docs/
一、Abstract
單細胞基因表達研究使人們能夠在復(fù)雜的生物過程和高度異質(zhì)性的細胞群中描述轉(zhuǎn)錄調(diào)控。這些研究有助于發(fā)現(xiàn)識別特定亞型細胞的基因，或標記生物過程中的中間狀態(tài)，以及在兩種不同的細胞命運之間過渡態(tài)的基因。在許多單細胞研究中，單個細胞是以不同步的方式執(zhí)行基因表達過程。實際上，每個細胞都是正在研究的轉(zhuǎn)錄過程中的一個瞬間。Monocle包是分析單細胞測序的工具。Monocle引入了在偽時間(擬時間)內(nèi)對單個細胞排序的策略，利用單個細胞的非同步進程，將它們置于與細胞分化等生物學過程相對應(yīng)的軌跡上。****Monocle利用先進的機器學習技術(shù)(反向圖嵌入)從單細胞數(shù)據(jù)中學習顯式的主圖來對細胞進行排序，這可以強大而準確地解決復(fù)雜的生物過程。Monocle也可以進行聚類(即使用t-SNE和密度峰值聚類)。Monocle也可以進行差異基因表達測試，使人們能夠識別在不同狀態(tài)下差異表達的基因，沿著生物過程以及不同的細胞命運時基因表達的變化。Monocle是專為單細胞RNA-Seq研究設(shè)計的，但也可以用于其他分析。
二、Introduction
Monocle 2包括新的和改進的算法，用于對細胞進行分類和計數(shù)，執(zhí)行細胞亞群之間的差異表達分析，以及重建細胞軌跡。Monocle 2也被重新設(shè)計，可以很好地完成包含數(shù)萬個或更多細胞的大型單細胞RNA-Seq數(shù)據(jù)分析。
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使用一種簡單的、無偏的和高度可擴展的統(tǒng)計程序來選擇具有軌跡進展特征的基因。然后，它采用了一類流形學習算法，旨在在高維單細胞RNA-seq數(shù)據(jù)中嵌入一個主圖。以前的方法是通過啟發(fā)式分析細胞之間的成對距離來推斷分支結(jié)構(gòu)，而Monocle 2可以使用這張圖來直接識別發(fā)育的命運決定。我們已經(jīng)通過廣泛的基準測試證明，Monocle 2優(yōu)于其他工具，如Wishbone，而不需要用戶指定軌跡的結(jié)構(gòu)。
三、Monocle2算法
Monocle 2的算法是再2017年發(fā)表在nature methods 上，文章鏈接如下
https://www.nature.com/articles/nmeth.4402
文章中的核心理論為：每個細胞都可以表示為高維空間中的一個點，在高維空間中，每個維對應(yīng)著一個有序基因的表達水平。高維數(shù)據(jù)首先通過幾種降維方法，如PCA(默認)、擴散映射等，投射到低維空間。Monocle 2然后在自動選擇的一組數(shù)據(jù)質(zhì)心上構(gòu)造一棵生成樹。然后，該算法將細胞移動到它們最近的樹的頂點，更新頂點的位置以適應(yīng)細胞，學習新的生成樹，并迭代地繼續(xù)這個過程，直到樹和細胞的位置已經(jīng)收斂。在這個過程中，Monocle 2保持了高維空間和低維空間之間的可逆映射，從而既學習了軌跡，又降低了數(shù)據(jù)的維數(shù)。****一旦Monocle 2學會了樹，用戶就會選擇一個tip作為根。計算每個單元的偽時間作為其沿樹到根的測地線距離，并根據(jù)主圖自動分配其分枝。因為monocle2學習樹結(jié)構(gòu)，與其他方法相比，分支結(jié)構(gòu)自動出現(xiàn)。當它更新細胞位置并細化樹時，monocle2簡化了軌跡的結(jié)構(gòu)，修剪了小的分支，這樣最終的軌跡只保留了描述細胞狀態(tài)顯著差異的分支。
四、Installing Monocle2
Monocle需要在R語言環(huán)境下運行。官網(wǎng)給出的都是R-3.4比較老的版本，現(xiàn)在都是R-3.6或者R-4.0以上版本。我們可以直接在Bioconductor中安裝Monocle 2。
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("monocle")
查看系統(tǒng)中安裝的這個包的版本文檔
browseVignettes("monocle")
當然也可以在GitHub上安裝最新版的Monocle2。
install.packages("devtools")
devtools::install_github("cole-trapnell-lab/monocle-release@develop")
有時我們添加的特性需要你安裝某些額外的軟件包。當您嘗試上面的命令時，可能會看到錯誤。您可以通過輸入(例如)來在錯誤消息中安裝包。
BiocManager::install(c("DDRTree", "pheatmap"))
最后檢查一下是否安裝成功
四、創(chuàng)建CellDataSet
方法一: 將Seurat object中數(shù)據(jù)提取來創(chuàng)建
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值通常是對數(shù)正態(tài)分布的，而UMIs或讀計數(shù)使用負二項更好地建模。要處理計數(shù)數(shù)據(jù)，請將負二項分布指定為newCellDataSet的expressionFamily參數(shù)：
稀疏矩陣用negbinomial.size()，
FPKM值用tobit()，
logFPKM值用gaussianff()
方法二: 直接讀取表達矩陣來創(chuàng)建
library(data.table)
##讀取數(shù)據(jù)
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)
##創(chuàng)建
pd <- new("AnnotatedDataFrame", data = pd)
fd <- new("AnnotatedDataFrame", data = fd)
HSMM <- newCellDataSet(as.matrix(data),
phenoData = pd, featureData = fd,
expressionFamily = tobit())
###如果數(shù)據(jù)量大，建議轉(zhuǎn)化為稀疏矩陣
HSMM <- newCellDataSet(as(as.matrix(data), "sparseMatrix"),
phenoData = pd,
featureData = fd,
expressionFamily = tobit())
方法三: 將Seurat對象直接轉(zhuǎn)化為CellDataSet對象
importCDS(scRNA.Osteoclastic)
如果我們想要從Seurat對象或SCESet中導(dǎo)入所有的插槽，我們可以設(shè)置參數(shù)'import_all'為TRUE。#(默認為FALSE或只保留最小數(shù)據(jù)集)
五、估計size factor和離散度
size facotr幫助我們標準化細胞之間的mRNA的差異。
離散度值可以幫助我們進行后續(xù)的差異分析。
monocle_cds <- estimateSizeFactors(monocle_cds)
monocle_cds <- estimateDispersions(monocle_cds)
六、過濾低質(zhì)量的細胞
大多數(shù)單細胞工作流程至少會包含一些由死細胞或空孔組成的庫。同樣重要的是要刪除doublets:由兩個或多個細胞意外生成的庫。這些細胞可以破壞下游步驟，如偽時間排序或聚類。要知道一個特定的基因有多少個表達，或者一個給定的細胞有多少個基因表達，通常是很方便的。Monocle提供了一個簡單的函數(shù)來計算這些統(tǒng)計數(shù)據(jù)。
monocle_cds <- detectGenes(monocle_cds, min_expr = 0.1)
print(head(fData(monocle_cds)))
七、細胞分類
Monocle官網(wǎng)教程提供了4個分類方法：
**Classifying cells by type ** Clustering cells without marker genes Clustering cells using marker genes Imputing cell type
1、Classifying cells by type
Monocle提供了一個簡單的系統(tǒng)，根據(jù)您選擇的標記基因的表達來標記細胞。您只需提供一組函數(shù)，Monocle就可以使用這些函數(shù)對每個單元格進行注釋。例如，您可以為幾種單元格類型中的每種類型提供一個函數(shù)。這些函數(shù)接受每個單元格的表達式數(shù)據(jù)作為輸入，并返回TRUE以告訴Monocle某個單元格滿足函數(shù)定義的條件。所以你可以有一個功能，對表達成肌細胞特異性基因的細胞來說是真的，另一個功能對成纖維細胞特異性基因來說是真的，等等。
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提供了一種算法，可以用來計算“未知”單元格的類型。該算法在函數(shù)clusterCells中實現(xiàn)，根據(jù)全局表達式輪廓將單元格分組在一起。這樣的話，如果你的細胞表達了許多特定于成肌細胞的基因，但恰好缺少MYF5，我們?nèi)匀豢梢宰R別它為成肌細胞。clusterCells可以在無監(jiān)督的方式下使用，也可以在半監(jiān)督的“”模式下使用，這允許使用一些專家知識來輔助算法。我們先來看看無監(jiān)督模式。
第一步是決定哪些基因用于細胞聚類。我們可以使用所有的基因，但是我們會包含很多基因，這些基因的表達水平不足以提供有意義的信號。包括它們只會給系統(tǒng)增加噪音。我們可以根據(jù)平均表達水平篩選基因，我們還可以選擇在細胞間異?？勺兊幕颉＿@些基因往往對細胞狀態(tài)有很高的信息量。
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函數(shù)標記了在隨后對clusterCells的調(diào)用中用于聚類的基因，盡管我們可以根據(jù)需要提供其他基因列表。plot_ordering_genes函數(shù)顯示了基因表達的變異性（離散性）如何依賴于細胞間的平均表達。紅線表示基于這種關(guān)系的單片體對色散的預(yù)期。我們標記用于聚類的基因顯示為黑點，而其他基因顯示為灰點。
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
首先，我們將選擇一組不同的基因用于細胞聚類。在我們挑選高表達和高度可變的基因之前?，F(xiàn)在，我們將挑選與標記共變異的基因。從某種意義上說，我們將建立一個大的基因列表來作為標記，這樣即使一個細胞沒有MYF5，它也可以被其他基因識別為成肌細胞。
marker_diff <- markerDiffTable(HSMM[expressed_genes,],
cth,
residualModelFormulaStr = "~Media + num_genes_expressed",
cores = 1)
函數(shù)markerDiffTable接受CellDataSet和CellTypeHierarchy，并根據(jù)提供的函數(shù)將所有單元格分類為類型。然后在識別不同基因類型之間差異表達的基因之前，它去除了所有“未知”和“不明確”的功能。同樣，您可以提供要從該測試中排除的影響的剩余模型。然后函數(shù)返回一個測試結(jié)果的數(shù)據(jù)幀，您可以使用它來選擇要用于聚類的基因。通常情況下，最好挑選最適合每種細胞類型的前10或20個基因。這確保了聚類基因不受一種細胞類型標記的支配。如果可能的話，通常需要為每種類型創(chuàng)建一個平衡的標記面板。Monocle提供了一個方便的功能，可以根據(jù)基因在每種類型中的表達受限程度來對基因進行排序。
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
注意，我們已經(jīng)減少了成肌細胞群中“污染”成纖維細胞的數(shù)量，反之亦然。但是“未知”細胞呢？如果為clusterCells提供CellTypeHierarcy，Monocle將使用它對整個簇進行分類，而不僅僅是單個細胞?；旧?，clusterCells的工作原理與以前完全相同，只是在構(gòu)建集群之后，它會統(tǒng)計每個集群中每種細胞類型的頻率。當一個簇由超過一定百分比（在本例中為10%）的特定類型組成時，簇中的所有單元格都被設(shè)置為該類型。如果一個集群由多個單元類型組成，則整個集群都被標記為“不明確”。如果沒有高于閾值的單元類型，則集群被標記為“未知”。因此，Monocle可以幫助您在缺少標記數(shù)據(jù)的情況下估算每個細胞的類型。
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())
八、構(gòu)建軌跡
在發(fā)育過程中，為了對刺激做出反應(yīng)，在整個生命周期中，細胞從一種功能“狀態(tài)”過渡到另一種功能“狀態(tài)”。處于不同狀態(tài)的細胞表達不同的基因，產(chǎn)生蛋白質(zhì)和代謝物的動態(tài)重復(fù)，從而完成它們的工作。當細胞在不同的狀態(tài)間移動時，會經(jīng)歷一個轉(zhuǎn)錄重組的過程，一些基因被沉默，而另一些則被激活。這些瞬態(tài)通常很難描述，因為在更穩(wěn)定的端點狀態(tài)之間凈化細胞可能很困難或不可能。單細胞RNA-Seq可以讓你看到這些狀態(tài)而不需要純化。然而，要做到這一點，我們必須確定每個細胞的可能狀態(tài)范圍。
Monocle介紹了利用RNA-Seq進行單細胞軌跡分析的策略。Monocle不是通過實驗將細胞凈化成離散狀態(tài)，而是使用一種算法來學習每個細胞作為動態(tài)生物過程的一部分必須經(jīng)歷的基因表達變化序列。一旦掌握了基因表達變化的整體“軌跡”，Monocle就可以將每個細胞放置在軌跡中的適當位置。然后，你可以使用Monocle的差異分析工具來尋找在軌跡過程中受調(diào)控的基因，如“尋找隨假時間變化的基因”一節(jié)所述。如果這個過程有多個結(jié)果，Monocle將重建一個“分支”軌跡。這些分支與細胞的“決定”相對應(yīng)，Monocle提供了強有力的工具來識別受它們影響并參與制造它們的基因。在“分析單細胞軌跡中的分支”部分中，可以看到如何分析分支。Monocle依靠一種稱為反向圖嵌入的機器學習技術(shù)來構(gòu)造單細胞軌跡.
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中，低水平表達的基因通常非常嘈雜，但有些基因包含有關(guān)細胞狀態(tài)的重要信息。單核細胞通過檢測這些基因在細胞群體中的表達模式來對細胞進行排序。Monocle尋找以“有趣的”（即不只是嘈雜的）方式變化的基因，并用它們來構(gòu)造數(shù)據(jù)。Monocle為您提供了多種工具來選擇基因，這些基因?qū)a(chǎn)生一個健壯、準確和具有生物學意義的軌跡。你可以使用這些工具來執(zhí)行一個完全“無監(jiān)督”的分析，在這個分析中，Monocle不知道你認為哪個基因是重要的?；蛘撸憧梢岳胢arker gene知識來定義生物學進展，從而形成Monocle的軌跡。我們認為這種模式是“半監(jiān)督”的分析。
Monocle官網(wǎng)教程提供了4個選擇方法：
選擇發(fā)育差異表達基因 選擇clusters差異表達基因 選擇離散程度高的基因 自定義發(fā)育marker基因
前三種都是無監(jiān)督分析方法，細胞發(fā)育軌跡生成完全不受人工干預(yù)；最后一種是半監(jiān)督分析方法，可以使用先驗知識輔助分析。
##使用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)
理想情況下，我們希望盡可能少地使用正在研究的系統(tǒng)生物學的先驗知識。我們希望從數(shù)據(jù)中發(fā)現(xiàn)重要的排序基因，而不是依賴于文獻和教科書，因為這可能會在排序中引入偏見。我們將從一種更簡單的方法開始，但是我們通常推薦一種更復(fù)雜的方法，稱為“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")
如果你的樹上有一堆狀態(tài)，那就很難弄清楚每個狀態(tài)都落在樹上的什么地方了。有時，可以方便地“分面”軌跡圖，以便更容易地查看每個狀態(tài)的位置.
plot_cell_trajectory(HSMM, color_by = "State") +
facet_wrap(~State, nrow = 1)
如果你沒有時間序列，你可能需要利用你對系統(tǒng)的生物學知識，根據(jù)特定標記基因的表達位置來設(shè)置根。例如，在這個實驗中，一個高度增殖的祖細胞群正在產(chǎn)生兩種類型的有絲分裂后細胞。因此，根部應(yīng)該有表達高水平增殖標記的細胞。我們可以使用抖動圖來判斷哪種狀態(tài)對應(yīng)于快速擴散：
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
九、差異分析
官方給出的差異分析有三大方法，我們重點關(guān)注第三個：****根據(jù)偽時間功能尋找差異基因
1、Basic Differential Analysis
2、Finding Genes that Distinguish Cell Type or State
3、Finding Genes that Change as a Function of Pseudotime
Monocle的主要工作是通過生物過程（如細胞分化）將細胞按順序排列，而不知道要提前查看哪些基因。一旦這樣做了，你就可以分析細胞，找到隨著細胞進展而變化的基因。例如，你可以發(fā)現(xiàn)隨著細胞“成熟”而顯著上調(diào)的基因。讓我們來看看一組對肌肉生成很重要的基因：
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函數(shù)指出Monocle應(yīng)該通過表達式值擬合自然樣條曲線，以幫助它將表達式的變化描述為進程的函數(shù)。讓我們加入基因注釋，這樣就很容易看出哪些基因是重要的。
十、根據(jù)偽時間表達pattern聚類基因
** 研究時間序列基因表達研究時出現(xiàn)的一個常見問題是：“哪些基因遵循相似的動力學趨勢”？Monocle可以通過對具有相似趨勢的基因進行分組來幫助你回答這個問題，所以你可以分析這些基因組，看看它們有什么共同點。Monocle提供了一種方便的方法來可視化所有假時間依賴性基因。函數(shù)plot_pseudotime_heatmap采用CellDataSet對象（通常只包含重要基因的子集），并生成平滑的表達曲線，非常類似于plot_genes_in_pseudotime。然后，它對這些基因進行聚類，并使用pheatmap軟件包進行繪圖。這讓你可以想象出基因模塊在不同的時間內(nèi)共同變化。**
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)
十一、單細胞軌跡的“分支”分析
通常，單細胞軌跡包括分支。這些分支的產(chǎn)生是因為細胞執(zhí)行不同的基因表達程序。在發(fā)育過程中，當細胞做出命運選擇時，分支出現(xiàn)在軌跡中：一個發(fā)育譜系沿著一條路徑前進，而另一個譜系產(chǎn)生第二條路徑。Monocle包含分析這些分支事件的廣泛功能。
以Steve Quake's的實驗室進行的實驗為例，Barbara Treutlein團隊從發(fā)育中的小鼠肺中捕獲了細胞。他們在發(fā)育早期捕獲細胞，后來當肺中同時含有兩種主要類型的上皮細胞（AT1和AT2）時，細胞就要決定變成AT1或AT2。Monocle可以將這個過程重構(gòu)為一個分支軌跡，允許您對決策點進行非常詳細的分析。下圖顯示了使用部分數(shù)據(jù)重建的Monocle軌跡。有一個分支，標記為“1”。當細胞從樹的左上角經(jīng)過樹枝的早期發(fā)育階段時，哪些基因會發(fā)生變化？哪些基因在分枝間有差異表達？為了回答這個問題，Monocle為您提供了一個特殊的統(tǒng)計測試：分支表達式分析建模，或BEAM。BEAM(Branched expression analysis modeling)是一種統(tǒng)計方法，用于尋找以依賴于分支的方式調(diào)控的基因。
plot_cell_trajectory(HSMM, color_by = "State")
BEAM進行統(tǒng)計分析
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)
轉(zhuǎn)載來自：
Monocle2 | 單細胞測序的擬時序分析（細胞軌跡分析）
相關(guān)資料：
官網(wǎng)
monocle2 擬時間分支點分析結(jié)果解讀
單細胞分析實錄(15): 基于monocle2的擬時序分析
Monocle2包學習筆記（四）：Differential Expression Analysis
|
