擬時（pseudotime）分析，又稱細(xì)胞軌跡（cell trajectory）分析慕淡，通過擬時分析可以推斷出發(fā)育過程細(xì)胞的分化軌跡或細(xì)胞亞型的演化過程，在發(fā)育相關(guān)研究中使用頻率較高蔚润。主要基于關(guān)鍵基因的表達(dá)模式训貌，在擬時間中對單個細(xì)胞進(jìn)行排序灌危，模擬出時間發(fā)育過程的動態(tài)變化康二。

既然講到所謂的擬時分析不過是一種排序而已，那么我們就要知道排序的基本要素：

• 對什么排序
• 如何判斷先后順序
• 如何尋找分支點(diǎn)（如果分支的話）

早在20世紀(jì)30年代勇蝙，數(shù)量生態(tài)學(xué)領(lǐng)域前蘇聯(lián)學(xué)者Ranensky就提出了排序的概念沫勿，并發(fā)展了一個簡單的排序方法（見Sobolev和Utekhin 1973），但只限于在前蘇聯(lián)傳播（Greig-Smith 1980）味混，Ramensky當(dāng)時應(yīng)用一個或兩個環(huán)境因子梯度去排列植物群落产雹，他用的名詞是德文“ordnung”。經(jīng)過多年的發(fā)展翁锡，排序技術(shù)依然是一種在低維空間排布高維數(shù)據(jù)的降維技術(shù)蔓挖。當(dāng)我們講到排序的時候，離不開降維盗誊；講到降維的時候时甚，離不開特征提取（或者選擇）哈踱。

那么對應(yīng)到擬時分析的描述中：1)關(guān)鍵基因就是特征選擇的結(jié)果荒适，2)擬時間就是排序空間，3)排序就是細(xì)胞的演化軌跡开镣。所有的擬時分析都離不開這三點(diǎn)刀诬。具體地：

Monocle引入了在偽時間(擬時間)內(nèi)對單個細(xì)胞排序的策略，利用單個細(xì)胞的非同步進(jìn)程邪财，將它們置于與細(xì)胞分化等生物學(xué)過程相對應(yīng)的軌跡上陕壹。Monocle利用先進(jìn)的機(jī)器學(xué)習(xí)技術(shù)(反向圖嵌入)從單細(xì)胞數(shù)據(jù)中學(xué)習(xí)顯式的主圖來對細(xì)胞進(jìn)行排序，該方法能夠可靠树埠、準(zhǔn)確地解決復(fù)雜的生物學(xué)過程糠馆。

Monocle2 可以做什么？

• Clustering, classifying, and counting cells
• Constructing single-cell trajectories.
• Differential expression analysis

當(dāng)然我們關(guān)心的是第二個功能了怎憋，但是不防也看一下它的其他功能又碌。算法原理可以參看文章Reversed graph embedding resolves complex single-cell developmental trajectories

Monocle 2通過反向圖嵌入學(xué)習(xí)單細(xì)胞軌跡的方法。每個細(xì)胞都可以表示為高維空間中的一個點(diǎn)绊袋，其中每個維對應(yīng)于有序基因的表達(dá)水平毕匀。高維數(shù)據(jù)首先通過PCA(默認(rèn)方法)、擴(kuò)散圖(diffusion maps)等幾種降維方法中的任何一種投影到較低的維空間癌别。然后Monocle 2在自動選擇的數(shù)據(jù)中心集上構(gòu)造一個生成樹皂岔。然后，該算法細(xì)胞移動到樹中最近的頂點(diǎn)展姐，不斷更新頂點(diǎn)的位置以“適合”細(xì)胞躁垛，學(xué)習(xí)新的生成樹，并迭代地繼續(xù)這個過程诞仓，直到樹和細(xì)胞的位置收斂為止缤苫。在整個過程中，Monocle 2保持了高維空間和低維空間之間的可逆映射墅拭，從而既學(xué)習(xí)了軌跡活玲，又降低了數(shù)據(jù)的維數(shù)。一旦Monocle 2習(xí)得這棵樹谍婉，并選擇一個tip作為“根”舒憾。每個細(xì)胞的偽時間計(jì)算為其沿樹到根的測地距離（Geodesic Distance，在圖論中穗熬，Geodesic Distance 就是圖中兩節(jié)點(diǎn)的最短路徑的距離）镀迂，并根據(jù)主圖自動分配其分支。

主要過程如下唤蔗，最難理解的應(yīng)該是反轉(zhuǎn)圖嵌入探遵，對算法感興趣的同學(xué)可以研究一下流行學(xué)習(xí)的知識窟赏。

• dpFeature selection， identifying genes that define biological process 確定和生物過程有關(guān)的基因
  selecting the genes differentially expressed between clusters of cells identified with tSNE dimension reduction followed by density peak clustering選擇差異表達(dá)基因箱季；
• 反向圖嵌入（GRE）：默認(rèn)DDRTree涯穷。它是一種流形學(xué)習(xí)算法，旨在將主要圖形嵌入高維單細(xì)胞RNA-seq數(shù)據(jù)中藏雏，找到高維基因表達(dá)空間與低維空間映射拷况，同時在這個縮小的空間中學(xué)習(xí)圖的結(jié)構(gòu)（K個質(zhì)心）來解決這個問題；
• 偽時間分配和分支計(jì)算：利用DDRTree構(gòu)建MST掘殴，將細(xì)胞投影到MST上赚瘦，MST遞歸計(jì)算偽時間；

算法層面適可而止奏寨。

seurat2monocle

library(monocle)
#<-importCDS(pbmc,import_all = TRUE)
#Load Seurat object
pbmc <- readRDS('D:\\Users\\Administrator\\Desktop\\RStudio\\single_cell\\filtered_gene_bc_matrices\\hg19pbmc_tutorial.rds')

#Extract data, phenotype data, and feature data from the SeuratObject
data <- as(as.matrix(pbmc@assays$RNA@counts), 'sparseMatrix')
pd <- new('AnnotatedDataFrame', data = pbmc@meta.data)
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
fd <- new('AnnotatedDataFrame', data = fData)

#Construct monocle cds
monocle_cds <- newCellDataSet(data,
                              phenoData = pd,
                              featureData = fd,
                              lowerDetectionLimit = 0.5,
                              expressionFamily = negbinomial.size())

如果你只有count矩陣可以這么讀入：

#do not run
HSMM_expr_matrix <- read.table("fpkm_matrix.txt")
HSMM_sample_sheet <- read.delim("cell_sample_sheet.txt")
HSMM_gene_annotation <- read.delim("gene_annotations.txt")

pd <- new("AnnotatedDataFrame", data = HSMM_sample_sheet)
fd <- new("AnnotatedDataFrame", data = HSMM_gene_annotation)
HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix),
    phenoData = pd, featureData = fd)

這里重點(diǎn)是理解CellDataSet 的數(shù)據(jù)結(jié)構(gòu)起意。

• exprs, a numeric matrix of expression values, where rows are genes, and columns are cells
• phenoData, an AnnotatedDataFrame object, where rows are cells, and columns are cell attributes** (such as cell type, culture condition, day captured, etc.)
• featureData, an AnnotatedDataFrame object, where rows are features (e.g. genes), and columns are gene attributes, such as biotype, gc content, etc.

官網(wǎng)給出許多其他建議以及讀入數(shù)據(jù)的注意事項(xiàng)，需要注意的是病瞳，如果是UMI數(shù)據(jù)杜恰，不應(yīng)該在創(chuàng)建CellDataSet之前自己對其進(jìn)行標(biāo)準(zhǔn)化。也不應(yīng)該試圖將UMI計(jì)數(shù)轉(zhuǎn)換為相對豐度(通過將其轉(zhuǎn)換為FPKM/TPM數(shù)據(jù))仍源。Monocle將在內(nèi)部執(zhí)行所有需要的標(biāo)準(zhǔn)化步驟心褐。自己將其正常化可能會破壞Monocle的一些關(guān)鍵步驟笼踩。

數(shù)據(jù)過濾（QC）

> HSMM<-monocle_cds
## 歸一化 
> HSMM <- estimateSizeFactors(HSMM)
> HSMM <- estimateDispersions(HSMM)
#Filtering low-quality cells
> HSMM <- detectGenes(HSMM, min_expr = 3 )
> print(head(fData(HSMM)))
              gene_short_name num_cells_expressed
AL627309.1         AL627309.1                   9
AP006222.2         AP006222.2                   3
RP11-206L10.2   RP11-206L10.2                   5
RP11-206L10.9   RP11-206L10.9                   3
LINC00115           LINC00115                  18
NOC2L                   NOC2L                 254
> expressed_genes <- row.names(subset(fData(HSMM),
+                                     num_cells_expressed >= 10))
> print(head(pData(HSMM)))
               orig.ident nCount_RNA nFeature_RNA percent.mt RNA_snn_res.0.5 seurat_clusters Size_Factor num_genes_expressed
AAACATACAACCAC     pbmc3k       2419          779  3.0177759               1               1          NA                 779
AAACATTGAGCTAC     pbmc3k       4903         1352  3.7935958               3               3          NA                1352
AAACATTGATCAGC     pbmc3k       3147         1129  0.8897363               1               1          NA                1129
AAACCGTGCTTCCG     pbmc3k       2639          960  1.7430845               2               2          NA                 960
AAACCGTGTATGCG     pbmc3k        980          521  1.2244898               6               6          NA                 521
AAACGCACTGGTAC     pbmc3k       2163          781  1.6643551               1               1          NA                 781

細(xì)胞分類（Classifying）

Monocle提供了一個簡單的系統(tǒng)逗爹，可以根據(jù)您選擇的marker基因的表達(dá)來標(biāo)記細(xì)胞。您只需提供Monocle可以用來注釋每個單元格的一組函數(shù)嚎于。

CDC20 <- row.names(subset(fData(HSMM), gene_short_name == "CDC20"))
GABPB2 <- row.names(subset(fData(HSMM),
                             gene_short_name == "GABPB2"))

cth <- newCellTypeHierarchy()
cth <- addCellType(cth, "CDC20", classify_func =
                     function(x) { x[CDC20,] >= 1 })
cth <- addCellType(cth, "GABPB2", classify_func = function(x){ x[GABPB2,] >= 1 })

HSMM <- classifyCells(HSMM, cth, 0.1)
table(pData(HSMM)$CellType)

 CDC20  GABPB2 Unknown 
      4      69    2565 

我只是隨意拿幾個基因來試一下

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())

聚類

第一步是決定用哪些基因來進(jìn)行聚類（特征選擇）掘而。我們可以使用所有的基因，但是我們將包括很多沒有表達(dá)到足夠高水平來提供有意義的信號的基因于购，它們只會給系統(tǒng)增加噪音袍睡。我們可以根據(jù)平均表達(dá)水平篩選基因，我們還可以選擇細(xì)胞間異常變異的基因肋僧。這些基因往往對細(xì)胞狀態(tài)有很高的信息含量斑胜。

#Clustering cells without marker genes 

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)

紅線表示單片基于這種關(guān)系對色散的期望。我們標(biāo)記用于聚類的基因用黑點(diǎn)表示嫌吠，其他的用灰點(diǎn)表示止潘。

# Plots the percentage of variance explained by the each component based on PCA from the normalized expression data using the same procedure used in reduceDimension function.
# HSMM@auxClusteringData[["tSNE"]]$variance_explained <- NULL
plot_pc_variance_explained(HSMM, return_all = F) # norm_method='log'

> HSMM <- reduceDimension(HSMM, max_components = 2, num_dim = 10,
+                         reduction_method = 'tSNE', verbose = T)
Remove noise by PCA ...
Reduce dimension by tSNE ...
> HSMM <- clusterCells(HSMM, num_clusters = 2)
Distance cutoff calculated to 2.589424 
> plot_cell_clusters(HSMM, 1, 2, color = "CellType",
+                    markers = c("CDC20", "GABPB2"))

我發(fā)現(xiàn)指定cluster為5的時候，只能畫出來4個,為什么辫诅？

HSMM <- clusterCells(HSMM, num_clusters = 5)
plot_cell_clusters(HSMM, 1, 2)

Monocle允許我們減去“無趣的”變異源的影響凭戴，以減少它們對集群的影響。您可以使用到clusterCells和其他幾個Monocle函數(shù)的residualmodelformula astr參數(shù)來實(shí)現(xiàn)這一點(diǎn)炕矮。此參數(shù)接受一個R模型公式字符串么夫，該字符串指定要在cluster之前減去者冤。

HSMM <- reduceDimension(HSMM, max_components = 2, num_dim = 2,
                        reduction_method = 'tSNE',
                        residualModelFormulaStr = "~Size_Factor + num_genes_expressed",
                        verbose = T)
HSMM <- clusterCells(HSMM, num_clusters = 5)
plot_cell_clusters(HSMM, 1, 2, color = "Cluster")  # 圖不放了

plot_cell_clusters(HSMM, 1, 2, color = "Cluster") +
  facet_wrap(~CellType)

monocle 支持半監(jiān)督聚類Clustering cells using marker genes來定義細(xì)胞以及外部導(dǎo)入細(xì)胞類型（Imputing cell type）（cell與type對應(yīng)文件），這里均不再展示档痪。

構(gòu)建軌跡

在發(fā)育過程中譬嚣，細(xì)胞會對刺激做出反應(yīng)，并在整個生命過程中钞它，從一種功能“狀態(tài)”轉(zhuǎn)變?yōu)榱硪环N功能“狀態(tài)”。不同狀態(tài)的細(xì)胞表達(dá)不同的基因殊鞭，產(chǎn)生蛋白質(zhì)和代謝物的動態(tài)重復(fù)序列遭垛，從而完成它們的工作。當(dāng)細(xì)胞在不同狀態(tài)間移動時操灿，會經(jīng)歷一個轉(zhuǎn)錄重組的過程锯仪，其中一些基因被沉默，另一些基因被激活趾盐。這些瞬時狀態(tài)通常很難描述庶喜，因?yàn)樵诟€(wěn)定的端點(diǎn)狀態(tài)之間純化細(xì)胞可能是困難的或不可能的。

單細(xì)胞RNA-Seq可以使您在不需要純化的情況下看到這些狀態(tài)救鲤。然而久窟，要做到這一點(diǎn)，我們必須確定每個細(xì)胞的可能狀態(tài)區(qū)間本缠。

Monocle介紹了利用RNA-Seq進(jìn)行單細(xì)胞軌跡分析的策略斥扛。Monocle不是通過實(shí)驗(yàn)將細(xì)胞純化成離散狀態(tài)，而是使用一種算法來學(xué)習(xí)每個細(xì)胞必須經(jīng)歷的基因表達(dá)變化序列丹锹，作為動態(tài)生物學(xué)過程的一部分稀颁。一旦它了解了基因表達(dá)變化的整體“軌跡”，Monocle就可以將每個細(xì)胞置于軌跡中的適當(dāng)位置楣黍。

然后匾灶，您可以使用Monocle的微分分析工具包來查找在軌跡過程中受到調(diào)控的基因，如查找作為偽時間函數(shù)變化的基因租漂。如果這個過程有多個結(jié)果阶女，Monocle將重建一個“分支”軌跡。這些分支與細(xì)胞的“決策”相對應(yīng)哩治，Monocle提供了強(qiáng)大的工具來識別受它們影響的基因张肾，并參與這些基因的形成以及如何進(jìn)行分析分支。Monocle依靠一種叫做反向圖嵌入的機(jī)器學(xué)習(xí)技術(shù)來構(gòu)建單細(xì)胞軌跡锚扎。

在開始之前您也需要問吞瞪，什么是擬時（偽時間）分析。

偽時間是衡量單個細(xì)胞在細(xì)胞分化等過程中取得了多大進(jìn)展的指標(biāo)驾孔。在許多生物學(xué)過程中芍秆，細(xì)胞并不是完全同步的惯疙。在細(xì)胞分化等過程的單細(xì)胞表達(dá)研究中，捕獲的細(xì)胞在分化方面可能分布廣泛妖啥。也就是說霉颠，在同一時間捕獲的細(xì)胞群中，有些細(xì)胞可能已經(jīng)很長時間了荆虱，而有些細(xì)胞甚至還沒有開始這個過程蒿偎。

當(dāng)您想要了解在細(xì)胞從一種狀態(tài)轉(zhuǎn)換到另一種狀態(tài)時所發(fā)生的調(diào)節(jié)更改的順序時，這種異步性會產(chǎn)生主要問題怀读。跟蹤同時捕獲的細(xì)胞間的表達(dá)可以產(chǎn)生對基因動力學(xué)一個大致的認(rèn)識诉位，該基因表達(dá)的明顯變異性將非常高。Monocle根據(jù)每個cell在學(xué)習(xí)軌跡上的進(jìn)展對其進(jìn)行排序菜枷，從而緩解了由于異步而產(chǎn)生的問題苍糠。

Monocle不是跟蹤表達(dá)式隨時間變化的函數(shù)，而是跟蹤沿軌跡變化的函數(shù)啤誊，我們稱之為偽時間岳瞭。偽時間是一個抽象的分化單位：它只是一個cell到軌跡起點(diǎn)的距離，沿著最短路徑測量蚊锹。軌跡的總長度是由細(xì)胞從起始狀態(tài)移動到結(jié)束狀態(tài)所經(jīng)歷的總轉(zhuǎn)錄變化量來定義的瞳筏。

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

基因的選擇

首先，我們必須決定我們將使用哪些基因來定義細(xì)胞在肌生成過程中的進(jìn)展牡昆。我們最終想要的是一組基因乏矾，能夠隨著我們研究過程的進(jìn)展而增加(或減少)表達(dá)。

理想情況下迁杨，我們希望盡可能少地使用正在研究的系統(tǒng)生物學(xué)的先驗(yàn)知識钻心。我們希望從數(shù)據(jù)中發(fā)現(xiàn)重要的排序基因，而不是依賴于文獻(xiàn)和教科書铅协，因?yàn)檫@可能會在排序中引入偏見捷沸。我們將從一種更簡單的方法開始，但是我們通常推薦一種更復(fù)雜的方法狐史，稱為“dpFeature”痒给。

diff_test_res <- differentialGeneTest(HSMM[expressed_genes,],
                                      fullModelFormulaStr = "~percent.mt")
ordering_genes <- row.names (subset(diff_test_res, qval < 0.1)) ## 不要也寫0.1 ，而是要寫0.01骏全。

HSMM <- setOrderingFilter(HSMM, ordering_genes)
plot_ordering_genes(HSMM)

根據(jù)時間點(diǎn)的差異分析來選擇基因通常是非常有效的苍柏，但是如果我們沒有時間序列數(shù)據(jù)呢?如果細(xì)胞在我們的生物過程中是異步分化的(通常是這樣)，Monocle通辰保可以從同時捕獲的單個種群重建它們的軌跡试吁。下面是兩種選擇基因的方法，它們完全不需要知道實(shí)驗(yàn)的設(shè)計(jì)。

一旦我們有了一個用于排序的基因id列表熄捍，我們就需要在HSMM對象中設(shè)置它們烛恤，因?yàn)榻酉聛淼膸讉€函數(shù)將依賴于它們。

降維

#Trajectory step 2: reduce data dimensionality  
HSMM <- reduceDimension(HSMM, max_components = 2,
                            method = 'DDRTree')

排序

#Trajectory step 3: order cells along the trajectory  
HSMM <- orderCells(HSMM)
plot_cell_trajectory(HSMM, color_by = "seurat_clusters")

有了樹結(jié)構(gòu)分類顏色是可以侄自己指定的余耽。軌跡呈樹形結(jié)構(gòu)缚柏，Monocle事先不知道應(yīng)該調(diào)用樹的哪個軌跡中的哪個來調(diào)用“開始”，所以我們經(jīng)常不得不再次使用root_state參數(shù)調(diào)用orderCells來指定開始碟贾。首先币喧，我們繪制軌跡，這次按“狀態(tài)”給細(xì)胞上色袱耽。

plot_cell_trajectory(HSMM, color_by = "State")

plot_cell_trajectory(HSMM, color_by = "Pseudotime")

“狀態(tài)”只是自定義的某術(shù)語杀餐，表示一段分支。下面的函數(shù)便于識別包含seurat_clusters編號為零的大多數(shù)細(xì)胞的狀態(tài)扛邑。然后我們可以把它傳遞給orderCells。這個函數(shù)只是為了演示如何自定義root_state 铐然，更簡單的方法是?orderCells 看如何設(shè)置蔬崩。

GM_state <- function(cds){
  if (length(unique(pData(cds)$State)) > 1){
    T0_counts <- table(pData(cds)$State, pData(cds)$seurat_clusters)[,"0"]
    return(as.numeric(names(T0_counts)[which
                                       (T0_counts == max(T0_counts))]))
  } else {
    return (1)
  }
}
HSMM_1 <- orderCells(HSMM, root_state = GM_state(HSMM))
plot_cell_trajectory(HSMM_1, color_by = "Pseudotime")

可以看出不同的排序條件擬時方向是不一樣的。同時如果想看每個狀態(tài)下細(xì)胞的分化可以指定分面搀暑，當(dāng)然你也可以指定其他的分面方式沥阳。

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

如果沒有timeseries，可能需要根據(jù)某些標(biāo)記基因的表達(dá)位置來設(shè)置根自点，使用您對系統(tǒng)的生物學(xué)知識桐罕。例如，在這個實(shí)驗(yàn)中桂敛，高度增殖的祖細(xì)胞群產(chǎn)生兩種有絲分裂后的細(xì)胞功炮。所以根應(yīng)該有表達(dá)高水平增殖標(biāo)記的細(xì)胞。我們可以使用抖動圖找出對應(yīng)于快速增殖的狀態(tài)

blast_genes <- row.names(subset(fData(HSMM),
                                gene_short_name %in% c("CCNB2", "NOC2L", "CDC20")))
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", "ICAM2", "ICAM2")))
cds_subset <- HSMM_filtered[my_genes,]
plot_genes_in_pseudotime(cds_subset, color_by = "seurat_clusters")

差異分析

差異基因表達(dá)分析是RNA-Seq實(shí)驗(yàn)中的一項(xiàng)常見任務(wù)薪伏。Monocle可以幫助你找到不同細(xì)胞群間差異表達(dá)的基因，并評估這些變化的統(tǒng)計(jì)顯著性粗仓。這些比較要求您有一種方法將細(xì)胞收集到兩個或更多組中嫁怀。這些組由每個CellDataSet的表現(xiàn)型數(shù)據(jù)表中的列定義。Monocle將評估不同細(xì)胞群中每個基因表達(dá)水平的顯著性借浊。

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 = "~percent.mt")
> head(diff_test_res)
      status           family       pval      qval gene_short_name num_cells_expressed use_for_ordering
MEF2D     OK negbinomial.size 0.10544931 0.4590520           MEF2D                   1            FALSE
CCNB1     OK negbinomial.size 0.45277970 0.7043240           CCNB1                   1            FALSE
MEF2C     OK negbinomial.size 0.28733226 0.5028315           MEF2C                   2            FALSE
TPM2      OK negbinomial.size 0.94348541 0.9696949            TPM2                   0            FALSE
CDK1      OK negbinomial.size 0.12768261 0.4590520            CDK1                   0            FALSE
CCND1     OK negbinomial.size 0.04021331 0.4590520           CCND1                   0            FALSE


MYOG_ID1 <- HSMM[row.names(subset(fData(HSMM),
                                      gene_short_name %in% c("CCND1", "CCNB2"))),]
plot_genes_jitter(MYOG_ID1, grouping = "seurat_clusters", ncol= 2)

選擇對區(qū)分細(xì)胞類型有意義的基因（Finding Genes that Distinguish Cell Type or State ）

#Finding Genes that Distinguish Cell Type or State 
to_be_tested <- row.names(subset(fData(HSMM),
                                 gene_short_name %in% c("CDC20", "NCAM1", "ANPEP")))
cds_subset <- HSMM[to_be_tested,]


diff_test_res <- differentialGeneTest(cds_subset,
                                      fullModelFormulaStr = "~CellType")
diff_test_res[,c("gene_short_name", "pval", "qval")]

plot_genes_jitter(cds_subset,
                  grouping = "CellType",
                  color_by = "CellType",
                  nrow= 1,
                  ncol = NULL,
                  plot_trend = TRUE)

full_model_fits <-
  fitModel(cds_subset,  modelFormulaStr = "~CellType")
reduced_model_fits <- fitModel(cds_subset, modelFormulaStr = "~1")
diff_test_res <- compareModels(full_model_fits, reduced_model_fits)
diff_test_res

      status           family        pval       qval
CDC20     OK negbinomial.size 0.007822883 0.02346865
NCAM1     OK negbinomial.size 0.791131484 0.88906005
ANPEP     OK negbinomial.size 0.889060052 0.88906005

Monocle中的差異分析過程非常靈活:您在測試中使用的模型公式可以包含pData表中作為列存在的任何項(xiàng)塘淑，包括Monocle在其他分析步驟中添加的列。例如蚂斤，如果您使用集群細(xì)胞存捺，您可以使用集群作為模型公式來測試cluster之間不同的基因。

Finding Genes that Change as a Function of Pseudotime

Monocle的主要工作是將細(xì)胞按照生物過程(如細(xì)胞分化)的順序排列曙蒸，而不事先知道要查看哪些基因召噩。一旦這樣做了母赵，你就可以分析細(xì)胞，找到隨著細(xì)胞進(jìn)步而變化的基因具滴。例如凹嘲，你可以發(fā)現(xiàn)當(dāng)細(xì)胞“成熟”時，基因顯著上調(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")]

      gene_short_name         pval         qval
MEF2C           MEF2C 2.996841e-02 3.995789e-02
CCNB2           CCNB2 5.328721e-03 1.065744e-02
MYH3             MYH3 6.371156e-01 6.371156e-01
TNNT1           TNNT1 1.634704e-12 6.538818e-12

plot_genes_in_pseudotime(cds_subset, color_by = "seurat_clusters")

Clustering Genes by Pseudotemporal Expression Pattern

在研究時間序列基因表達(dá)研究時周蹭，一個常見的問題是:“哪些基因遵循相似的動力學(xué)趨勢?”Monocle可以通過將具有相似趨勢的基因分組來幫助你回答這個問題，這樣你就可以分析這些基因組疲恢，看看它們有什么共同之處凶朗。Monocle提供了一種方便的方法來可視化所有偽時間相關(guān)的基因。函數(shù)plot_pseudotime_heatmap接受一個CellDataSet對象(通常只包含重要基因的子集)显拳，并生成與plot_genes_in_pseudotime類似的平滑表達(dá)曲線.然后棚愤，它將這些基因聚類并使用pheatmap包繪制它們。這允許您可視化基因模塊杂数，這些模塊在偽時間內(nèi)共同變化宛畦。

注意下面的num_clusters 指的是基因可以聚成幾個類，而不是細(xì)胞揍移。

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)

Multi-Factorial Differential Expression Analysis

Monocle可以在多個因素存在的情況下進(jìn)行差異分析次和，這可以幫助你減去一些因素來看到其他因素的影響。在下面的簡單例子中那伐，Monocle測試了三個基因在成肌細(xì)胞和成纖維細(xì)胞之間的差異表達(dá)踏施，同時減去percent.mt的影響。為此罕邀，我們必須同時指定完整模型和簡化模型畅形。完整的模型同時捕捉細(xì)胞類型和percent.mt的影響。

to_be_tested <-
  row.names(subset(fData(HSMM),
                   gene_short_name %in% c("TPM1", "MYH3", "CCNB2", "GAPDH")))

cds_subset <- HSMM[to_be_tested,]

diff_test_res <- differentialGeneTest(cds_subset,
                                      fullModelFormulaStr = "~CellType + percent.mt",
                                      reducedModelFormulaStr = "~percent.mt")
diff_test_res[,c("gene_short_name", "pval", "qval")]

      gene_short_name       pval      qval
GAPDH           GAPDH 0.07990737 0.1598147
CCNB2           CCNB2 0.04148172 0.1598147
TPM1             TPM1 0.90861287 0.9086129
MYH3             MYH3 0.77085745 0.9086129

plot_genes_jitter(cds_subset,
                  grouping = "seurat_clusters", color_by = "CellType", plot_trend = TRUE) +
  facet_wrap( ~ feature_label, scales= "free_y")

Analyzing Branches in Single-Cell Trajectories

通常诉探，單細(xì)胞的軌跡包括分支束亏。分支的產(chǎn)生是因?yàn)榧?xì)胞執(zhí)行替代基因表達(dá)程序。在發(fā)育過程中阵具，當(dāng)細(xì)胞做出命運(yùn)的選擇時碍遍，分支就會出現(xiàn)在軌跡中:一個發(fā)育譜系沿著一條路徑前進(jìn)，而另一個譜系產(chǎn)生第二條路徑阳液。Monocle包含用于分析這些分支事件的廣泛功能怕敬。

芭芭拉·特雷特琳和她的同事們在史蒂夫·奎克的實(shí)驗(yàn)室里進(jìn)行了一項(xiàng)實(shí)驗(yàn)，他們從正在發(fā)育的老鼠肺中獲取細(xì)胞帘皿。他們在細(xì)胞發(fā)育的早期捕獲細(xì)胞东跪，之后當(dāng)肺包含兩種主要類型的上皮細(xì)胞(AT1和AT2)，以及即將決定成為AT1或AT2的細(xì)胞時。Monocle可以將這個過程重構(gòu)為一個分支軌跡虽填，讓你可以非常詳細(xì)地分析決策點(diǎn)丁恭。下圖顯示了使用Monocle的一些數(shù)據(jù)重建的一般方法。有一個單獨(dú)的分支斋日，標(biāo)記為“1”牲览。當(dāng)細(xì)胞從樹的左上方通過樹枝的早期發(fā)育階段通過時，哪些基因發(fā)生了變化?哪些基因在分支間有差異表達(dá)?為了回答這個問題恶守，Monocle為您提供了一個特殊的統(tǒng)計(jì)測試:分支表達(dá)式分析建模第献，或BEAM。

lung <- load_lung()
plot_cell_trajectory(lung, color_by = "Time")

# 圖如下

BEAM_res <- BEAM(lung, branch_point = 1, cores = 1)
BEAM_res <- BEAM_res[order(BEAM_res$qval),]
BEAM_res <- BEAM_res[,c("gene_short_name", "pval", "qval")]

使用一種特殊類型的熱圖兔港，您可以可視化所有明顯依賴于分支的基因的變化庸毫。這張熱圖同時顯示了兩種血統(tǒng)的變化。它還要求您選擇要檢查的分支點(diǎn)衫樊。列是偽時間中的點(diǎn)飒赃，行是基因，偽時間的開始在熱圖的中間科侈。當(dāng)你從熱圖的中間讀到右邊的時候资厉，你正在跟隨一個偽時間譜系前计。當(dāng)你讀到左邊時得封，另一個汽久。這些基因是分層聚類的羡洛，因此您可以可視化具有類似的依賴于序列的基因模塊.

plot_genes_branched_heatmap(lung[row.names(subset(BEAM_res,
                                                  qval < 1e-4)),],
                            branch_point = 1,
                            num_clusters = 4,
                            cores = 1,
                            use_gene_short_name = T,
                            show_rownames = T)

那么問題來了挂脑，這個圖應(yīng)該怎么看呢？

github上面是這樣回答的：

https://github.com/cole-trapnell-lab/monocle-release/issues/219
monocle2 擬時間分支點(diǎn)分析結(jié)果解讀

lung_genes <- row.names(subset(fData(lung),
                               gene_short_name %in% c("Ccnd2", "Sftpb", "Pdpn")))
plot_genes_branched_pseudotime(lung[lung_genes,],
                               branch_point = 1,
                               color_by = "Time",
                               ncol = 1)

> str(HSMM)
Formal class 'CellDataSet' [package "monocle"] with 19 slots
  ..@ reducedDimS          : num [1:2, 1:2638] -2.613 -0.652 -1.784 0.745 -2.376 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2638] "AAACATACAACCAC" "AAACATTGAGCTAC" "AAACATTGATCAGC" "AAACCGTGCTTCCG" ...
  ..@ reducedDimW          : num [1:281, 1:2] -0.0384 0.0812 -0.0465 0.0681 0.0655 ...
  ..@ reducedDimA          : num [1:2, 1:2638] 13.57 29.4 -2.74 -25.16 -9.3 ...
  ..@ reducedDimK          : num [1:2, 1:126] -2.305 -0.692 -1.711 0.538 -1.935 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:126] "Y_1" "Y_2" "Y_3" "Y_4" ...
  ..@ minSpanningTree      :List of 10
  .. ..$ :List of 1
  .. .. ..$ Y_1: 'igraph.vs' Named int [1:2] 47 101
  .. .. .. ..- attr(*, "names")= chr [1:2] "Y_47" "Y_101"
  .. .. .. ..- attr(*, "env")=<weakref> 
  .. .. .. ..- attr(*, "graph")= chr "484d230c-94de-11e9-9f7f-07844984f7db"
  .. ..$ :List of 1
  .. .. ..$ Y_2: 'igraph.vs' Named int [1:2] 62 90
  .. .. .. ..- attr(*, "names")= chr [1:2] "Y_62" "Y_90"
  .. .. .. ..- attr(*, "env")=<weakref> 
  .. .. .. ..- attr(*, "graph")= chr "484d230c-94de-11e9-9f7f-07844984f7db"
  .. ..$ :List of 1
  .. .. ..$ Y_3: 'igraph.vs' Named int [1:2] 90 126
  .. .. .. ..- attr(*, "names")= chr [1:2] "Y_90" "Y_126"
  .. .. .. ..- attr(*, "env")=<weakref> 
  .. .. .. ..- attr(*, "graph")= chr "484d230c-94de-11e9-9f7f-07844984f7db"
  .. ..$ :List of 1
  .. .. ..$ Y_4: 'igraph.vs' Named int 106
  .. .. .. ..- attr(*, "names")= chr "Y_106"
  .. .. .. ..- attr(*, "env")=<weakref> 
  .. .. .. ..- attr(*, "graph")= chr "484d230c-94de-11e9-9f7f-07844984f7db"
  .. ..$ :List of 1
  .. .. ..$ Y_5: 'igraph.vs' Named int [1:2] 36 65
  .. .. .. ..- attr(*, "names")= chr [1:2] "Y_36" "Y_65"
  .. .. .. ..- attr(*, "env")=<weakref> 
  .. .. .. ..- attr(*, "graph")= chr "484d230c-94de-11e9-9f7f-07844984f7db"
  .. ..$ :List of 1
  .. .. ..$ Y_6: 'igraph.vs' Named int [1:2] 39 53
  .. .. .. ..- attr(*, "names")= chr [1:2] "Y_39" "Y_53"
  .. .. .. ..- attr(*, "env")=<weakref> 
  .. .. .. ..- attr(*, "graph")= chr "484d230c-94de-11e9-9f7f-07844984f7db"
  .. ..$ :List of 1
  .. .. ..$ Y_7: 'igraph.vs' Named int [1:2] 34 98
  .. .. .. ..- attr(*, "names")= chr [1:2] "Y_34" "Y_98"
  .. .. .. ..- attr(*, "env")=<weakref> 
  .. .. .. ..- attr(*, "graph")= chr "484d230c-94de-11e9-9f7f-07844984f7db"
  .. ..$ :List of 1
  .. .. ..$ Y_8: 'igraph.vs' Named int [1:2] 49 126
  .. .. .. ..- attr(*, "names")= chr [1:2] "Y_49" "Y_126"
  .. .. .. ..- attr(*, "env")=<weakref> 
  .. .. .. ..- attr(*, "graph")= chr "484d230c-94de-11e9-9f7f-07844984f7db"
  .. ..$ :List of 1
  .. .. ..$ Y_9: 'igraph.vs' Named int [1:2] 42 93
  .. .. .. ..- attr(*, "names")= chr [1:2] "Y_42" "Y_93"
  .. .. .. ..- attr(*, "env")=<weakref> 
  .. .. .. ..- attr(*, "graph")= chr "484d230c-94de-11e9-9f7f-07844984f7db"
  .. ..$ :List of 1
  .. .. ..$ Y_10: 'igraph.vs' Named int [1:2] 11 25
  .. .. .. ..- attr(*, "names")= chr [1:2] "Y_11" "Y_25"
  .. .. .. ..- attr(*, "env")=<weakref> 
  .. .. .. ..- attr(*, "graph")= chr "484d230c-94de-11e9-9f7f-07844984f7db"
  .. ..- attr(*, "class")= chr "igraph"
  ..@ cellPairwiseDistances: num [1:126, 1:126] 0 1.37 1.41 10.62 5.52 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:126] "Y_1" "Y_2" "Y_3" "Y_4" ...
  .. .. ..$ : chr [1:126] "Y_1" "Y_2" "Y_3" "Y_4" ...
  ..@ expressionFamily     :Formal class 'vglmff' [package "VGAM"] with 22 slots
  .. .. ..@ blurb             : chr [1:6] "Negative-binomial distribution with size known\n\n" "Links:    " "loglink(mu)" "\n" ...
  .. .. ..@ constraints       :  expression({  constraints <- cm.zero.VGAM(constraints, x = x, NULL, M = M, predictors.names = predictors.names, M1 = 2) })
  .. .. ..@ deviance          :function ()  
  .. .. ..@ fini              :  expression({ })
  .. .. ..@ first             :  expression({ })
  .. .. ..@ infos             :function (...)  
  .. .. ..@ initialize        :  expression({  M1 <- 1  temp5 <- w.y.check(w = w, y = y, Is.nonnegative.y = TRUE, Is.integer.y = TRUE, ncol.w.max | __truncated__
  .. .. ..@ last              :  expression({  misc$link <- rep_len("loglink", NOS)  names(misc$link) <- mynames1  misc$earg <- vector("list", M) | __truncated__
  .. .. ..@ linkfun           :function ()  
  .. .. ..@ linkinv           :function (eta, extra = NULL)  
  .. .. ..@ loglikelihood     :function (mu, y, w, residuals = FALSE, eta, extra = NULL, summation = TRUE)  
  .. .. ..@ middle            :  expression({ })
  .. .. ..@ middle2           :  expression({ })
  .. .. ..@ summary.dispersion: logi(0) 
  .. .. ..@ vfamily           : chr "negbinomial.size"
  .. .. ..@ validparams       :function (eta, y, extra = NULL)  
  .. .. ..@ validfitted       :function ()  
  .. .. ..@ simslot           :function (object, nsim)  
  .. .. ..@ hadof             :function ()  
  .. .. ..@ charfun           :function ()  
  .. .. ..@ deriv             :  expression({  eta <- cbind(eta)  M <- ncol(eta)  kmat <- matrix(Inf, n, M, byrow = TRUE)  newemu <- list(theta = | __truncated__
  .. .. ..@ weight            :  expression({  ned2l.dmunb2 <- 1/mu - 1/(mu + kmat)  wz <- ned2l.dmunb2 * dmu.deta^2  c(w) * wz })
  ..@ lowerDetectionLimit  : num 1
  ..@ dispFitInfo          :<environment: 0x000000023e5e8110> 
  ..@ dim_reduce_type      : chr "DDRTree"
  ..@ auxOrderingData      :<environment: 0x000000002c5323c8> 
  ..@ auxClusteringData    :<environment: 0x00000000240fc8a0> 
  ..@ experimentData       :Formal class 'MIAME' [package "Biobase"] with 13 slots
  .. .. ..@ name             : chr ""
  .. .. ..@ lab              : chr ""
  .. .. ..@ contact          : chr ""
  .. .. ..@ title            : chr ""
  .. .. ..@ abstract         : chr ""
  .. .. ..@ url              : chr ""
  .. .. ..@ pubMedIds        : chr ""
  .. .. ..@ samples          : list()
  .. .. ..@ hybridizations   : list()
  .. .. ..@ normControls     : list()
  .. .. ..@ preprocessing    : list()
  .. .. ..@ other            : list()
  .. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
  .. .. .. .. ..@ .Data:List of 2
  .. .. .. .. .. ..$ : int [1:3] 1 0 0
  .. .. .. .. .. ..$ : int [1:3] 1 1 0
  ..@ assayData            :<environment: 0x000000005575c7a0> 
  ..@ phenoData            :Formal class 'AnnotatedDataFrame' [package "Biobase"] with 4 slots
  .. .. ..@ varMetadata      :'data.frame': 17 obs. of  1 variable:
  .. .. .. ..$ labelDescription: chr [1:17] NA NA NA NA ...
  .. .. ..@ data             :'data.frame': 2638 obs. of  17 variables:
  .. .. .. ..$ orig.ident                     : Factor w/ 1 level "pbmc3k": 1 1 1 1 1 1 1 1 1 1 ...
  .. .. .. ..$ nCount_RNA                     : num [1:2638] 2419 4903 3147 2639 980 ...
  .. .. .. ..$ nFeature_RNA                   : int [1:2638] 779 1352 1129 960 521 781 782 790 532 550 ...
  .. .. .. ..$ percent.mt                     : num [1:2638] 3.02 3.79 0.89 1.74 1.22 ...
  .. .. .. ..$ RNA_snn_res.0.5                : Factor w/ 9 levels "0","1","2","3",..: 2 4 2 3 7 2 5 5 1 6 ...
  .. .. .. ..$ seurat_clusters                : Factor w/ 9 levels "0","1","2","3",..: 2 4 2 3 7 2 5 5 1 6 ...
  .. .. .. ..$ Size_Factor                    : num [1:2638] 1.108 2.245 1.441 1.208 0.449 ...
  .. .. .. ..$ num_genes_expressed            : int [1:2638] 126 181 149 153 34 103 104 112 64 51 ...
  .. .. .. ..$ Cluster                        : Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
  .. .. .. ..$ peaks                          : logi [1:2638] FALSE FALSE FALSE FALSE FALSE FALSE ...
  .. .. .. ..$ halo                           : logi [1:2638] FALSE FALSE FALSE FALSE FALSE FALSE ...
  .. .. .. ..$ delta                          : num [1:2638] 0.401 0.442 1.031 0.205 0.812 ...
  .. .. .. ..$ rho                            : num [1:2638] 36 33.3 44.8 47.7 31.7 ...
  .. .. .. ..$ nearest_higher_density_neighbor: logi [1:2638] NA NA NA NA NA NA ...
  .. .. .. ..$ Pseudotime                     : num [1:2638] 12.277 11.115 13.023 0.723 16.159 ...
  .. .. .. ..$ State                          : Factor w/ 5 levels "1","2","3","4",..: 5 2 5 1 5 4 5 5 5 1 ...
  .. .. .. ..$ CellType                       : Factor w/ 3 levels "CDC20","GABPB2",..: 3 3 3 3 3 3 3 3 3 3 ...
  .. .. ..@ dimLabels        : chr [1:2] "sampleNames" "sampleColumns"
  .. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
  .. .. .. .. ..@ .Data:List of 1
  .. .. .. .. .. ..$ : int [1:3] 1 1 0
  ..@ featureData          :Formal class 'AnnotatedDataFrame' [package "Biobase"] with 4 slots
  .. .. ..@ varMetadata      :'data.frame': 3 obs. of  1 variable:
  .. .. .. ..$ labelDescription: chr [1:3] NA NA NA
  .. .. ..@ data             :'data.frame': 13714 obs. of  3 variables:
  .. .. .. ..$ gene_short_name    : Factor w/ 13714 levels "7SK.2","A1BG",..: 527 715 9394 9395 5971 7345 5750 8158 9747 4917 ...
  .. .. .. ..$ num_cells_expressed: int [1:13714] 0 0 0 0 0 1 0 0 0 6 ...
  .. .. .. ..$ use_for_ordering   : logi [1:13714] FALSE FALSE FALSE FALSE FALSE FALSE ...
  .. .. ..@ dimLabels        : chr [1:2] "featureNames" "featureColumns"
  .. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
  .. .. .. .. ..@ .Data:List of 1
  .. .. .. .. .. ..$ : int [1:3] 1 1 0
  ..@ annotation           : chr(0) 
  ..@ protocolData         :Formal class 'AnnotatedDataFrame' [package "Biobase"] with 4 slots
  .. .. ..@ varMetadata      :'data.frame': 0 obs. of  1 variable:
  .. .. .. ..$ labelDescription: chr(0) 
  .. .. ..@ data             :'data.frame': 2638 obs. of  0 variables
  .. .. ..@ dimLabels        : chr [1:2] "sampleNames" "sampleColumns"
  .. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
  .. .. .. .. ..@ .Data:List of 1
  .. .. .. .. .. ..$ : int [1:3] 1 1 0
  ..@ .__classVersion__    :Formal class 'Versions' [package "Biobase"] with 1 slot
  .. .. ..@ .Data:List of 5
  .. .. .. ..$ : int [1:3] 3 5 2
  .. .. .. ..$ : int [1:3] 2 42 0
  .. .. .. ..$ : int [1:3] 1 3 0
  .. .. .. ..$ : int [1:3] 1 0 0
  .. .. .. ..$ : int [1:3] 1 2 0

STREAM:一款scRNA-seq擬時分析工具
Clustering by fast search and find of density peaks
發(fā)表在 Science 上的一種新聚類算法
重磅消息：諾禾致源10X Genomics單細(xì)胞轉(zhuǎn)錄組產(chǎn)品全新升級
基迪奧10x單細(xì)胞轉(zhuǎn)錄組產(chǎn)品新升級
An overview of algorithms for estimating pseudotime in single-cell RNA-seq data
A set of tools supporting the development, execution, and benchmarking of trajectory inference methods.
hemberg-lab.github.io
Diffusion pseudotime robustly reconstructs lineage branching
dynverse
Monocle
A new statistical method allows researchers to infer different developmental processes from RNA-Seq data
Pseudotime estimation: deconfounding single cell time series
Your top 3 single-cell RNA sequencing analysis tools
Single-cell trajectory analysis
STREAM
MIA: Luca Pinello, Huidong Chen, Single-cell trajectories from omics; Jonathan Hsu, CRISPR tiling
使用monocle包進(jìn)行pseudotime分析
關(guān)于細(xì)胞分化軌跡學(xué)習(xí)小筆記