軌跡分析系列：

• 單細(xì)胞之軌跡分析-1：RNA velocity

擬時(shí)間序列分析（Pseudotime分析）的字面意思是通過構(gòu)建細(xì)胞間的變化軌跡來重塑細(xì)胞隨著時(shí)間的變化過程该贾。從具體的分類分析和復(fù)雜程度來說，可以分為細(xì)胞軌跡分析和細(xì)胞譜系分析。
其原理此前已經(jīng)介紹過，參考單細(xì)胞測(cè)序的軌跡推斷俺驶。
細(xì)胞軌跡分析指的是簡單模型的細(xì)胞變化軌跡分析勃蜘，通常指的是細(xì)胞沿著某個(gè)過程有特定化的變化終點(diǎn)，軌跡具有簡單樹狀結(jié)構(gòu)糊识，一端是“根”玷氏，另一端是“葉”堵未；細(xì)胞譜系分析通常指的是某類祖源細(xì)胞，在特定條件下盏触，有多個(gè)發(fā)育軌跡和命運(yùn)渗蟹，變化過程類似復(fù)雜樹狀分支變化過程。因此赞辩，簡單細(xì)胞軌跡分析和細(xì)胞譜系分析原理上類似雌芽，復(fù)雜程度有所區(qū)別，當(dāng)然辨嗽，基于此的分析手法和方式也會(huì)有所不同世落。

偽時(shí)間是衡量單個(gè)細(xì)胞在細(xì)胞分化等過程中取得了多大進(jìn)展的指標(biāo)。在許多生物學(xué)過程中糟需，細(xì)胞并不是完全同步的屉佳。在細(xì)胞分化等過程的單細(xì)胞表達(dá)研究中，捕獲的細(xì)胞在分化方面可能分布廣泛洲押。也就是說忘古，在同一時(shí)間捕獲的細(xì)胞群中，有些細(xì)胞可能已經(jīng)很長時(shí)間了诅诱，而有些細(xì)胞甚至還沒有開始這個(gè)過程。當(dāng)想要了解在細(xì)胞從一種狀態(tài)轉(zhuǎn)換到另一種狀態(tài)時(shí)所發(fā)生的調(diào)節(jié)更改的順序時(shí)送朱，這種異步性會(huì)產(chǎn)生主要問題娘荡。跟蹤同時(shí)捕獲的細(xì)胞間的表達(dá)可以產(chǎn)生對(duì)基因動(dòng)力學(xué)一個(gè)大致的認(rèn)識(shí)，該基因表達(dá)的明顯變異性將非常高驶沼。Monocle根據(jù)每個(gè)cell在學(xué)習(xí)軌跡上的進(jìn)展對(duì)其進(jìn)行排序炮沐，從而緩解了由于異步而產(chǎn)生的問題。Monocle不是跟蹤表達(dá)式隨時(shí)間變化的函數(shù)回怜，而是跟蹤沿軌跡變化的函數(shù)大年，我們稱之為偽時(shí)間换薄。偽時(shí)間是一個(gè)抽象的分化單位：它只是一個(gè)cell到軌跡起點(diǎn)的距離，沿著最短路徑測(cè)量翔试。軌跡的總長度是由細(xì)胞從起始狀態(tài)移動(dòng)到結(jié)束狀態(tài)所經(jīng)歷的總轉(zhuǎn)錄變化量來定義的轻要。

生信技能樹此前推送過一篇文章：擬時(shí)序分析就是差異分析的細(xì)節(jié)剖析，指出擬時(shí)序分析就是差異分析的細(xì)節(jié)剖析垦缅，所有的大樣本量差異分析都可以轉(zhuǎn)為擬時(shí)序分析冲泥，可以加深我們對(duì)擬時(shí)序分析的理解。

Monocle2是做單細(xì)胞擬時(shí)分析最有名的R包壁涎。
相較還在持續(xù)開發(fā)中的Monocle3來說凡恍，Monocle2更穩(wěn)定且更傾向于半監(jiān)督的分析模式，更適合針對(duì)感興趣的細(xì)胞亞群做個(gè)性化分析怔球。

monocle2 文章鏈接：https://www.nature.com/articles/nmeth.4402

文章中的核心理論為：每個(gè)細(xì)胞都可以表示為高維空間中的一個(gè)點(diǎn)嚼酝，在高維空間中，每個(gè)維對(duì)應(yīng)著一個(gè)有序基因的表達(dá)水平竟坛。高維數(shù)據(jù)首先通過幾種降維方法闽巩，如PCA(默認(rèn))、擴(kuò)散映射等流码，投射到低維空間又官。Monocle 2然后在自動(dòng)選擇的一組數(shù)據(jù)質(zhì)心上構(gòu)造一棵生成樹（DDRTree算法）。然后漫试，該算法將細(xì)胞移動(dòng)到它們最近的樹的頂點(diǎn)六敬，更新頂點(diǎn)的位置以適應(yīng)細(xì)胞，學(xué)習(xí)新的生成樹驾荣，并迭代地繼續(xù)這個(gè)過程外构，直到樹和細(xì)胞的位置已經(jīng)收斂。在這個(gè)過程中播掷，Monocle 2保持了高維空間和低維空間之間的可逆映射审编，從而既學(xué)習(xí)了軌跡，又降低了數(shù)據(jù)的維數(shù)歧匈。一旦Monocle 2學(xué)會(huì)了樹垒酬，用戶就會(huì)選擇一個(gè)tip作為根。計(jì)算每個(gè)單元的偽時(shí)間作為其沿樹到根的測(cè)地線距離件炉，并根據(jù)主圖自動(dòng)分配其分枝勘究。因?yàn)閙onocle2學(xué)習(xí)樹結(jié)構(gòu)，與其他方法相比斟冕，分支結(jié)構(gòu)自動(dòng)出現(xiàn)口糕。當(dāng)它更新細(xì)胞位置并細(xì)化樹時(shí)，monocle2簡化了軌跡的結(jié)構(gòu)磕蛇，修剪了小的分支景描，這樣最終的軌跡只保留了描述細(xì)胞狀態(tài)顯著差異的分支十办。

monocle2官網(wǎng)：http://cole-trapnell-lab.github.io/monocle-release/docs/

摘要

單細(xì)胞水平的研究使人們可以描述復(fù)雜生理過程和高度異質(zhì)性細(xì)胞群體的轉(zhuǎn)錄調(diào)控。這些研究有助于發(fā)現(xiàn)識(shí)別特定細(xì)胞亞型的基因超棺、標(biāo)記生物過程中間狀態(tài)的基因向族，以及在兩種不同的細(xì)胞命運(yùn)之間過渡態(tài)的基因。在許多單細(xì)胞研究中说搅，單個(gè)細(xì)胞以不同步的方式執(zhí)行基因表達(dá)過程炸枣。實(shí)際上，每個(gè)細(xì)胞都是正在研究的轉(zhuǎn)錄過程的一個(gè)瞬間弄唧。Monocle包是分析單細(xì)胞測(cè)序的工具适肠。
Monocle引入了在偽時(shí)間(擬時(shí)間)內(nèi)對(duì)單個(gè)細(xì)胞排序的策略，利用單個(gè)細(xì)胞的非同步進(jìn)程候引，將它們置于與細(xì)胞分化等生物學(xué)過程相對(duì)應(yīng)的軌跡上侯养。Monocle利用先進(jìn)的機(jī)器學(xué)習(xí)技術(shù)(反向圖嵌入)從單細(xì)胞數(shù)據(jù)中學(xué)習(xí)顯式的主圖（展現(xiàn)細(xì)胞轉(zhuǎn)錄特征相似性關(guān)系的圖，Monocle2使用DDTree降維圖澄干，Monocle3使用UMAP降維圖）來對(duì)細(xì)胞進(jìn)行排序逛揩，Monocle的機(jī)器學(xué)習(xí)算法可以依據(jù)上述降維圖形，學(xué)習(xí)描述細(xì)胞如何從一種狀態(tài)過渡到另一種狀態(tài)的軌跡麸俘。Monocle假設(shè)軌跡是樹狀結(jié)構(gòu)辩稽，一端是“根”，另一端是“葉”从媚。一個(gè)細(xì)胞在生物過程的開始逞泄，從根開始沿著主干進(jìn)行，直到它到達(dá)第一個(gè)分支拜效。然后喷众，該細(xì)胞必須選擇一條路徑，并沿著樹移動(dòng)越來越遠(yuǎn)紧憾，直到它到達(dá)一片葉子到千。一個(gè)細(xì)胞的假時(shí)間值是它返回根所需的距離。降維方面monocle與seurat的過程大同小異赴穗，首先進(jìn)行數(shù)據(jù)標(biāo)準(zhǔn)化憔四，其次選擇部分基因代表細(xì)胞轉(zhuǎn)錄特征 ，最后選用適當(dāng)?shù)乃惴ń稻S般眉。這可以強(qiáng)大而準(zhǔn)確地解決復(fù)雜的生物過程加矛。
Monocle也可以進(jìn)行聚類(即使用t-SNE和密度峰值聚類)和差異基因表達(dá)測(cè)試，使人們能夠識(shí)別在不同狀態(tài)下差異表達(dá)的基因煤篙，沿著生物過程以及不同的細(xì)胞命運(yùn)時(shí)基因表達(dá)的變化。Monocle是專為單細(xì)胞RNA-Seq研究設(shè)計(jì)的毁腿，但也可以用于其他分析辑奈。

Introduction

Monocle 2包括新的和改進(jìn)的算法用于細(xì)胞分類和計(jì)數(shù)苛茂，執(zhí)行細(xì)胞亞群之間的差異表達(dá)分析，以及細(xì)胞軌跡重建鸠窗。

Monocle主要可以進(jìn)行以下三種分析：

• 細(xì)胞的聚類妓羊、分類和計(jì)數(shù)。 Single-cell RNA-Seq experiments allow you to discover new (and possibly rare) subtypes of cells. Monocle helps you identify them.
• 重建單細(xì)胞軌跡。 In development, disease, and throughout life, cells transition from one state to another. Monocle helps you discover these transitions.
• 差異表達(dá)分析。 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使用一種簡單的村象、無偏的和高度可擴(kuò)展的統(tǒng)計(jì)程序來選擇具有軌跡進(jìn)展特征的基因。然后笋鄙，它采用了一類流形學(xué)習(xí)算法，旨在在高維單細(xì)胞RNA-seq數(shù)據(jù)中嵌入一個(gè)主圖。以前的方法是通過啟發(fā)式分析細(xì)胞之間的成對(duì)距離來推斷分支結(jié)構(gòu)淹父，而Monocle 2可以使用這張圖來直接識(shí)別發(fā)育的命運(yùn)決定。我們已經(jīng)通過廣泛的基準(zhǔn)測(cè)試證明怎虫，Monocle 2優(yōu)于其他工具暑认，如Wishbone，而不需要用戶指定軌跡的結(jié)構(gòu)大审。

Monocle的使用

分析流程示意圖：

1. 安裝

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

或者

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

2. 創(chuàng)建CellDataSet

2.1 將Seurat object中數(shù)據(jù)提取來創(chuàng)建??
數(shù)據(jù)準(zhǔn)備：pbmc3k數(shù)據(jù)下載（軌跡分析的前提是待分析的細(xì)胞有緊密的發(fā)育關(guān)系蘸际，PBMC細(xì)胞不是很好的的示例數(shù)據(jù)，在此僅作為演示徒扶。）
由于該數(shù)據(jù)集沒有對(duì)細(xì)胞類型進(jìn)行注釋粮彤，因此我們參考seurat標(biāo)準(zhǔn)流程對(duì)這個(gè)數(shù)據(jù)集進(jìn)行注釋。

library(dplyr)
library(Seurat)
library(patchwork) #用來做拼圖的包酷愧，后面的p1|p1|p3在一張圖上展示三個(gè)圖就是這個(gè)包的功勞

# Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
# Initialize the Seurat object with the raw (non-normalized data).
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
pbmc <- RunUMAP(pbmc, dims = 1:10)
new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono", 
    "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
pbmc[['cell_type']] <- pbmc@active.ident #將注釋結(jié)果添加到metadata
levels(pbmc)
saveRDS(pbmc, file = "pbmc.rds")

運(yùn)行完上述代碼驾诈，得到注釋好的數(shù)據(jù)集。

從Seurat對(duì)象中提取構(gòu)建CDS對(duì)象所需要的3個(gè)輸入文件：表達(dá)矩陣信息溶浴、基因信息和表型信息
??monocle輸入的是count矩陣（不建議使用data矩陣）

library(monocle)
pbmc <- readRDS("pbmc.rds") #導(dǎo)入注釋好的seurat對(duì)象（已注釋）

##提取表型信息--細(xì)胞信息(建議載入細(xì)胞的聚類或者細(xì)胞類型鑒定信息乍迄、實(shí)驗(yàn)條件等信息)
expr_matrix <- as(as.matrix(pbmc@assays$RNA@counts), 'sparseMatrix')
##提取表型信息到p_data(phenotype_data)里面 
p_data <- pbmc@meta.data 
p_data$celltype <- pbmc@active.ident  ##整合每個(gè)細(xì)胞的細(xì)胞鑒定信息到p_data里面。如果已經(jīng)添加則不必重復(fù)添加
##提取基因信息 如生物類型士败、gc含量等
f_data <- data.frame(gene_short_name = row.names(pbmc),row.names = row.names(pbmc))
##expr_matrix的行數(shù)與f_data的行數(shù)相同(gene number), expr_matrix的列數(shù)與p_data的行數(shù)相同(cell number)

#構(gòu)建CDS對(duì)象
pd <- new('AnnotatedDataFrame', data = p_data) 
fd <- new('AnnotatedDataFrame', data = f_data)
#將p_data和f_data從data.frame轉(zhuǎn)換AnnotatedDataFrame對(duì)象闯两。
cds <- newCellDataSet(expr_matrix,
                      phenoData = pd,
                      featureData = fd,
                      lowerDetectionLimit = 0.5,
                      expressionFamily = negbinomial.size())

FPKM/TPM值通常是對(duì)數(shù)正態(tài)分布的，而UMIs或讀計(jì)數(shù)使用負(fù)二項(xiàng)更好地建模谅将。要處理計(jì)數(shù)數(shù)據(jù)漾狼，需要將負(fù)二項(xiàng)分布指定為newCellDataSet的expressionFamily參數(shù)：

negbinomial.size()和negbinomial()：輸入的表達(dá)矩陣為UMI,一般適用于10x的數(shù)據(jù)；negbinomial()的結(jié)果更準(zhǔn)確饥臂，但是計(jì)算更耗時(shí)逊躁；一般建議采用negbinomial.size()。
稀疏矩陣用negbinomial.size()
tobit():適用于輸入的表達(dá)矩陣為FPKM或者TPM, 構(gòu)建monocle2的class時(shí)會(huì)自動(dòng)進(jìn)行l(wèi)og化計(jì)算
gaussianff():輸入為log化后的FPKM或者TPM隅熙。(目前在單細(xì)胞數(shù)據(jù)中稽煤，F(xiàn)PKM已不多用核芽，smart-seq2平臺(tái)數(shù)據(jù)一般采用TPM)

2.2 直接讀取表達(dá)矩陣來創(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())

2.3 將Seurat對(duì)象直接轉(zhuǎn)化為CellDataSet對(duì)象

importCDS(pbmc)

如果我們想要從Seurat對(duì)象或SCESet中導(dǎo)入所有的插槽酵熙，我們可以設(shè)置參數(shù)'import_all'為TRUE轧简。#(默認(rèn)為FALSE或只保留最小數(shù)據(jù)集)

3. 估計(jì)size factor和離散度

size facotr幫助我們標(biāo)準(zhǔn)化細(xì)胞之間的mRNA的差異。
離散度值可以幫助我們進(jìn)行后續(xù)的差異分析匾二。
（類似于seurat的數(shù)據(jù)歸一化處理）

cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)

與seurat把標(biāo)準(zhǔn)化后的表達(dá)矩陣保存在對(duì)象中不同哮独，monocle只保存一些中間結(jié)果在對(duì)象中，需要用時(shí)再用這些中間結(jié)果轉(zhuǎn)化察藐。經(jīng)過上面三個(gè)函數(shù)的計(jì)算皮璧，mycds對(duì)象中多了SizeFactors、Dipersions转培、num_cells_expressed和num_genes_expressed等信息恶导。

4. 過濾低質(zhì)量的細(xì)胞

大多數(shù)單細(xì)胞工作流程至少會(huì)包含一些由死細(xì)胞或空孔組成的庫。同樣重要的是要?jiǎng)h除doublets:由兩個(gè)或多個(gè)細(xì)胞意外生成的庫浸须。這些細(xì)胞可以破壞下游步驟惨寿，如偽時(shí)間排序或聚類。要知道一個(gè)特定的基因有多少個(gè)表達(dá)删窒，或者一個(gè)給定的細(xì)胞有多少個(gè)基因表達(dá)裂垦，通常是很方便的。Monocle提供了一個(gè)簡單的函數(shù)來計(jì)算這些統(tǒng)計(jì)數(shù)據(jù)肌索。
因?yàn)镾eurat已經(jīng)完成細(xì)胞過濾蕉拢，此步可省略
但由于Seurat是通過基因表達(dá)量來對(duì)細(xì)胞進(jìn)行過濾。因此诚亚，我們可以通過下面的代碼用表達(dá)某基因的細(xì)胞的數(shù)目對(duì)基因進(jìn)行過濾晕换，從而得到后續(xù)操作需要的基因。

cds <- detectGenes(cds, min_expr = 0.1) #這一操作會(huì)在fData(cds)中添加一列num_cells_expressed
print(head(fData(cds)))#此時(shí)有13714個(gè)基因
expressed_genes <- row.names(subset(fData(cds),
    num_cells_expressed >= 10)) #過濾掉在小于10個(gè)細(xì)胞中表達(dá)的基因站宗，還剩11095個(gè)基因闸准。

5. 細(xì)胞分類

Monocle官網(wǎng)教程提供了4個(gè)分類方法：
Classifying cells by type
Clustering cells without marker genes
Clustering cells using marker genes
Imputing cell type
建議先將細(xì)胞注釋好再進(jìn)行Monocle分析，不建議使用monocle做細(xì)胞分類梢灭。

6. 軌跡定義基因選擇及可視化和構(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: 選擇定義過程的基因
推斷單細(xì)胞軌跡是一個(gè)機(jī)器學(xué)習(xí)問題夷家。在單細(xì)胞RNA-Seq中，低水平表達(dá)的基因通常非常嘈雜敏释，但有些基因包含有關(guān)細(xì)胞狀態(tài)的重要信息库快。因此，軌跡推斷的第一步就是選擇Monocle將用作機(jī)器學(xué)習(xí)方法輸入的基因钥顽。這叫做特征選擇义屏，它對(duì)軌跡的形狀有很大的影響。
Monocle主要基于關(guān)鍵基因的表達(dá)模式，通過學(xué)習(xí)每個(gè)細(xì)胞必須經(jīng)歷的基因表達(dá)變化的序列湿蛔，根據(jù)擬時(shí)間值中對(duì)單個(gè)細(xì)胞進(jìn)行排序膀曾，模擬出時(shí)間發(fā)育過程的動(dòng)態(tài)變化。而這個(gè)排序技術(shù)表現(xiàn)是一種在低維空間排布高維數(shù)據(jù)的降維技術(shù)阳啥。
Monocle提供了多種工具來選擇基因，這些基因?qū)a(chǎn)生一個(gè)健壯财喳、準(zhǔn)確和具有生物學(xué)意義的軌跡察迟。你可以使用這些工具來執(zhí)行一個(gè)完全“無監(jiān)督”的分析，在這個(gè)分析中耳高，Monocle不知道你認(rèn)為哪個(gè)基因是重要的扎瓶。或者泌枪，你可以利用marker gene知識(shí)來定義生物學(xué)進(jìn)展概荷，從而形成Monocle的軌跡。我們認(rèn)為這種模式是“半監(jiān)督”的分析碌燕。

• Monocle官網(wǎng)教程提供了4個(gè)選擇方法：
  選擇發(fā)育差異表達(dá)基因
  選擇clusters差異表達(dá)基因
  選擇離散程度高的基因
  自定義發(fā)育marker基因

前三種都是無監(jiān)督分析方法误证，細(xì)胞發(fā)育軌跡生成完全不受人工干預(yù)；最后一種是半監(jiān)督分析方法修壕，可以使用先驗(yàn)知識(shí)輔助分析愈捅。

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

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

#這一步輸入的expressed_genes來自于步驟4芦昔。
#????后續(xù)分析使用的是該方法
#也可輸入seurat篩選出的高變基因：expressed_genes <- VariableFeatures(pbmc) 
diff <- differentialGeneTest(cds[expressed_genes,],fullModelFormulaStr="~cell_type",cores=1) 
#~后面是表示對(duì)誰做差異分析的變量，理論上可以為p_data的任意列名
head(diff)

##差異表達(dá)基因作為軌跡構(gòu)建的基因,差異基因的選擇標(biāo)準(zhǔn)是qval<0.01,decreasing=F表示按數(shù)值增加排序
deg <- subset(diff, qval < 0.01) #選出2829個(gè)基因
deg <- deg[order(deg$qval,decreasing=F),]
head(deg)

##差異基因的結(jié)果文件保存
write.table(deg,file="train.monocle.DEG.xls",col.names=T,row.names=F,sep="\t",quote=F)

## 軌跡構(gòu)建基因可視化
ordergene <- rownames(deg) 
cds <- setOrderingFilter(cds, ordergene)  
#這一步是很重要的胖腾，在我們得到想要的基因列表后烟零，我們需要使用setOrderingFilter將它嵌入cds對(duì)象，后續(xù)的一系列操作都要依賴于這個(gè)list咸作。
#setOrderingFilter之后锨阿，這些基因被儲(chǔ)存在cds@featureData@data[["use_for_ordering"]]，可以通過table(cds@featureData@data[["use_for_ordering"]])查看
pdf("train.ordergenes.pdf")
plot_ordering_genes(cds)
dev.off()
#出的圖黑色的點(diǎn)表示用來構(gòu)建軌跡的差異基因记罚，灰色表示背景基因墅诡。紅色的線是根據(jù)第2步計(jì)算的基因表達(dá)大小和離散度分布的趨勢(shì)(可以看到，找到的基因?qū)儆陔x散度比較高的基因)

??選擇的用于排序的基因數(shù)目一般在2000左右比較合適

gene數(shù)太多的話也可以選擇top基因

ordergene <- row.names(deg)[order(deg$qval)][1:400]

這樣用的文獻(xiàn)：腫瘤浸潤髓系細(xì)胞的pan-cancer單細(xì)胞轉(zhuǎn)錄圖譜

Step 2: 降維
一旦細(xì)胞有序排列，我們就可以在降維空間中可視化軌跡末早。所以首先選擇用于細(xì)胞排序的基因烟馅，然后使用反向圖嵌入(DDRTree)算法對(duì)數(shù)據(jù)進(jìn)行降維。

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

Step 3: 擬時(shí)間軸軌跡構(gòu)建和在擬時(shí)間內(nèi)排列細(xì)胞
將表達(dá)數(shù)據(jù)投射到更低的維度空間然磷，通過機(jī)器學(xué)習(xí)描述細(xì)胞如何從一種狀態(tài)過渡到另一種狀態(tài)的軌跡郑趁。假設(shè)軌跡具有樹狀結(jié)構(gòu)，一端是“根”姿搜，另一端是“葉”寡润。盡可能地將最佳樹與數(shù)據(jù)匹配起來。這項(xiàng)任務(wù)被稱為“歧管學(xué)習(xí)”舅柜，在生物過程的開始階段梭纹，細(xì)胞從根部開始，沿著主干前進(jìn)致份，直到到達(dá)第一個(gè)分支（如果有的話）变抽。然后，細(xì)胞必須選擇一條路徑氮块，沿著樹走得越來越遠(yuǎn)绍载，直到到達(dá)一片葉子。一個(gè)細(xì)胞的偽時(shí)間值是它回到根的距離雇锡。

根據(jù)order gene的表達(dá)趨勢(shì)逛钻，將細(xì)胞排序并完成軌跡構(gòu)建

cds <- orderCells(cds)
#??使用root_state參數(shù)可以設(shè)置擬時(shí)間軸的根，如下面的擬時(shí)間著色圖中可以看出锰提，左邊是根曙痘。根據(jù)state圖可以看出，根是State1立肘，若要想把另一端設(shè)為根边坤，可以按如下操作
#cds <- orderCells(cds, root_state = 5) #把State5設(shè)成擬時(shí)間軸的起始點(diǎn)

可視化：根據(jù)cds@phenoData@data中的表型信息(metadata)對(duì)細(xì)胞上色

1. 以pseudotime值上色 (Pseudotime是monocle2基于細(xì)胞基因表達(dá)信息計(jì)算的概率，表示時(shí)間的先后谅年。)

pdf("train.monocle.pseudotime.pdf",width = 7,height = 7)
plot_cell_trajectory(cds,color_by="Pseudotime", size=1,show_backbone=TRUE) 
dev.off()

圖中123是分叉點(diǎn)茧痒，不代表特別意義，也不代表時(shí)間先后融蹂。這張圖按照sudotime旺订，時(shí)間先后是從左往右，左邊是起點(diǎn)（root)超燃。起點(diǎn)的設(shè)置可以使用orderCells的root_state參數(shù)区拳，將右側(cè)設(shè)置為起點(diǎn)。

2. 以細(xì)胞類型上色

pdf("train.monocle.celltype.pdf",width = 7,height = 7)
plot_cell_trajectory(cds,color_by="cell_type", size=1,show_backbone=TRUE)
dev.off()

3. 以細(xì)胞狀態(tài)上色

pdf("train.monocle.state.pdf",width = 7,height = 7)
plot_cell_trajectory(cds, color_by = "State",size=1,show_backbone=TRUE)
dev.off()

state的多少是monocle算出來的意乓，不能調(diào)整樱调，與輸入的用于軌跡學(xué)習(xí)的基因有關(guān)。分叉和頂點(diǎn)之間或者頂點(diǎn)和頂點(diǎn)之間為一個(gè)state，與發(fā)育軌跡時(shí)間先后沒有關(guān)系笆凌，與細(xì)胞類型也不完全相關(guān)圣猎。

4. 按照seurat分群排序細(xì)胞

pdf("seurat.clusters.pdf",width = 7,height = 7)
plot_cell_trajectory(cds, color_by = "seurat_clusters")
dev.off()

5. 以細(xì)胞狀態(tài)上色（拆分）“分面”軌跡圖，以便更容易地查看每個(gè)狀態(tài)的位置乞而。

pdf("train.monocle.state.faceted.pdf",width = 10,height = 7)
plot_cell_trajectory(cds, color_by = "State") + facet_wrap("~State", nrow = 1)
dev.off()

??這一部分的配色可以直接對(duì)接ggsci送悔，也可以使用scale_color_manual()自己設(shè)置顏色

library(ggsci)
p1=plot_cell_trajectory(cds, color_by = "cell_type")  + scale_color_npg() 
p2=plot_cell_trajectory(cds, color_by = "State")  + scale_color_nejm()
colour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB","#98FB98","#F08080")
p3=plot_cell_trajectory(cds, color_by = "State")  + scale_color_manual(values = colour)
p1|p2|p3

覺得這種軌跡不太直觀的也可以畫成樹形圖

p1 <- plot_cell_trajectory(cds, x = 1, y = 2, color_by = "celltype") + 
  theme(legend.position='none',panel.border = element_blank()) + #去掉第一個(gè)的legend
  scale_color_manual(values = colour) 
p2 <- plot_complex_cell_trajectory(cds, x = 1, y = 2,
                                   color_by = "celltype")+
  scale_color_manual(values = colour) +
  theme(legend.title = element_blank()) 
p1|p2

還可以畫沿時(shí)間軸的細(xì)胞密度圖

library(ggpubr)
df <- pData(cds) 
## pData(cds)取出的是cds對(duì)象中cds@phenoData@data的內(nèi)容
View(df)
ggplot(df, aes(Pseudotime, colour = cell_type, fill=cell_type)) +
  geom_density(bw=0.5,size=1,alpha = 0.5)+theme_classic2()

從上面的矩陣中也能提取不同state的細(xì)胞，分別畫圖

手動(dòng)設(shè)置顏色注意

ClusterName_color_panel <- c(
  "Naive CD4 T" = "#DC143C", "Memory CD4 T" = "#0000FF", "CD14+ Mono" = "#20B2AA",
  "B" = "#FFA500", "CD8 T" = "#9370DB", "FCGR3A+ Mono" = "#98FB98",
  "NK" = "#F08080", "DC" = "#0000FF", "Platelet" = "#20B2AA"
)
ggplot(df, aes(Pseudotime, colour = cell_type, fill=cell_type)) +
  geom_density(bw=0.5,size=1,alpha = 0.5)+theme_classic2()+ scale_fill_manual(name = "", values = ClusterName_color_panel)+scale_color_manual(name = "", values = ClusterName_color_panel)

提取感興趣的細(xì)胞（進(jìn)行后續(xù)分析）

#比如對(duì)State7的細(xì)胞感興趣
pdata <- Biobase::pData(cds)
s.cells <- subset(pdata, State=="7") %>% rownames()
save(s.cells, file = "Monocle_state7.rda")

保存結(jié)果

write.csv(pData(cds), "pseudotime.csv")
save(cds, file = "cds.rda")

7. 指定基因的可視化

直接查看一些基因隨細(xì)胞狀態(tài)等的表達(dá)變化

##選擇前4個(gè)top基因并將其對(duì)象取出
keygenes <- head(ordergene,4)
cds_subset <- cds[keygenes,]
##可視化：以state/celltype/pseudotime進(jìn)行
p1 <- plot_genes_in_pseudotime(cds_subset, color_by = "State")
p2 <- plot_genes_in_pseudotime(cds_subset, color_by = "cell_type")
p3 <- plot_genes_in_pseudotime(cds_subset, color_by = "Pseudotime")
plotc <- p1|p2|p3
ggsave("Genes_pseudotimeplot.pdf", plot = plotc, width = 16, height = 8)

#指定基因
s.genes <- c("SELL","CCR7","IL7R", "CD84","CCL5","S100A4")
p1 <- plot_genes_jitter(cds[s.genes,], grouping = "State", color_by = "State")
p2 <- plot_genes_violin(cds[s.genes,], grouping = "State", color_by = "State")
p3 <- plot_genes_in_pseudotime(cds[s.genes,], color_by = "State")
plotc <- p1|p2|p3
ggsave("Genes_Jitterplot.pdf", plot = plotc, width = 16, height = 8)

由于state1-7并不代表分化時(shí)間先后爪模，因此這張圖還要結(jié)合其他的圖來看

??擬時(shí)序展示單個(gè)基因表達(dá)量

colnames(pData(cds))
pData(cds)$CCL5 = log2( exprs(cds)['CCL5',]+1)
p1=plot_cell_trajectory(cds, color_by = "CCL5")  + scale_color_gsea()
pData(cds)$S100A4 = log2(exprs(cds)['S100A4',]+1)
p2=plot_cell_trajectory(cds, color_by = "S100A4")    + scale_color_gsea()
library(patchwork)
p1+p2

8. 尋找擬時(shí)相關(guān)的基因（擬時(shí)差異基因）

使用回歸算法
注意：不要使用多核運(yùn)算放祟，經(jīng)常會(huì)出現(xiàn)警告
Monocle的主要工作是通過生物過程（如細(xì)胞分化）將細(xì)胞按順序排列，而不知道要提前查看哪些基因呻右。一旦這樣做了，你就可以分析細(xì)胞鞋喇，找到隨著細(xì)胞進(jìn)展而變化的基因声滥。

官方給出的差異分析有三大方法：
1、Basic Differential Analysis
2侦香、Finding Genes that Distinguish Cell Type or State
3落塑、Finding Genes that Change as a Function of Pseudotime
我們重點(diǎn)關(guān)注第三個(gè)：根據(jù)偽時(shí)間功能尋找差異基因

sm.ns函數(shù)指出Monocle應(yīng)該通過表達(dá)式值擬合自然樣條曲線，以幫助它將表達(dá)式的變化描述為進(jìn)程的函數(shù)罐韩。

• 尋找擬時(shí)差異基因（qvalue體現(xiàn)基因與擬時(shí)的密切程度）繪制熱圖

#這里是把排序基因（ordergene）提取出來做回歸分析憾赁，來找它們是否跟擬時(shí)間有顯著的關(guān)系
#如果不設(shè)置，就會(huì)用所有基因來做它們與擬時(shí)間的相關(guān)性
Time_diff <- differentialGeneTest(cds[ordergene,], cores = 1, 
                                  fullModelFormulaStr = "~sm.ns(Pseudotime)")
Time_diff <- Time_diff[,c(5,2,3,4,1,6,7)] #把gene放前面散吵，也可以不改
write.csv(Time_diff, "Time_diff_all.csv", row.names = F)
Time_genes <- Time_diff %>% pull(gene_short_name) %>% as.character()
p=plot_pseudotime_heatmap(cds[Time_genes,], num_clusters=4, show_rownames=T, return_heatmap=T)
ggsave("Time_heatmapAll.pdf", p, width = 5, height = 10)

?? ????前面通過設(shè)置num_clusters將熱圖聚成了四個(gè)cluster龙考，如過想要把每個(gè)cluster的基因單獨(dú)提出來做分析

p$tree_row
# Call:
# hclust(d = d, method = method)
# Cluster method   : ward.D2 
# Number of objects: 2829 
clusters <- cutree(p$tree_row, k = 4)
clustering <- data.frame(clusters)
clustering[,1] <- as.character(clustering[,1])
colnames(clustering) <- "Gene_Clusters"
table(clustering)
# 1    2    3    4 
# 570 1031  506  722 
write.csv(clustering, "Time_clustering_all.csv", row.names = F)

這樣就把每個(gè)基因?qū)儆谀膫€(gè)cluster提取出來了，后續(xù)可以做每個(gè)cluster的富集分析矾睦。參考：擬時(shí)序分析的熱圖提取基因問題

• 擬時(shí)差異基因熱圖繪制（提取了前100個(gè)）

Time_genes <- top_n(Time_diff, n = 100, desc(qval)) %>% pull(gene_short_name) %>% as.character()
p = plot_pseudotime_heatmap(cds[Time_genes,], num_clusters=4, show_rownames=T, return_heatmap=T)
ggsave("Time_heatmapTop100.pdf", p, width = 5, height = 10)

這個(gè)圖的橫軸是擬時(shí)間晦款，cluster1的基因是在擬時(shí)排序起點(diǎn)高表達(dá)的基因，cluster2的基因則是在擬時(shí)排序的重點(diǎn)高表達(dá)的枚冗。cluster數(shù)的多少是由plot_pseudotime_heatmap函數(shù)中的num_clusters參數(shù)定義的

注：plot_pseudotime_heatmap函數(shù)可以來可視化所有monocle的假時(shí)間依賴性基因缓溅。plot_pseudotime_heatmap采用CellDataSet對(duì)象（通常只包含重要基因的子集），并生成平滑的表達(dá)曲線赁温，非常類似于plot_genes_in_pseudotime坛怪。然后它對(duì)這些基因進(jìn)行聚類，并使用pheatmap軟件包進(jìn)行繪圖股囊。繪出的熱圖可以讓我們觀測(cè)到假時(shí)間依賴性基因中的不同基因模塊在不同的時(shí)間內(nèi)共同變化袜匿，能比較好的回答時(shí)間序列基因表達(dá)中“哪些基因遵循相似的動(dòng)力學(xué)趨勢(shì)”這一常見問題。

• 顯著差異基因按熱圖結(jié)果排序并保存

hp.genes <- p$tree_row$labels[p$tree_row$order]
Time_diff_sig <- Time_diff[hp.genes, c("gene_short_name", "pval", "qval")]
write.csv(Time_diff_sig, "Time_diff_sig.csv", row.names = F)

另外：也可手動(dòng)選擇基因來繪制熱圖毁涉，查看其表達(dá)模式

marker_genes <- row.names(subset(fData(cds),
                                 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(cds[marker_genes,],
                                      fullModelFormulaStr = "~sm.ns(Pseudotime)")
sig_gene_names <- row.names(subset(diff_test_res, qval < 0.1))
plot_pseudotime_heatmap(cds[sig_gene_names,],
                        num_clusters = 6,
                        cores = 1,
                        show_rownames = T)

9. 單細(xì)胞軌跡的“分支”分析

??上一步尋找擬時(shí)相關(guān)的基因是全局的沉帮，找擬時(shí)起點(diǎn)和終點(diǎn)相關(guān)的基因，這一步則是尋找和分叉點(diǎn)相關(guān)的基因。

單細(xì)胞軌跡常常包括分支穆壕。這些分支的產(chǎn)生是因?yàn)榧?xì)胞執(zhí)行不同的基因表達(dá)程序待牵。在發(fā)育過程中，當(dāng)細(xì)胞做出命運(yùn)選擇時(shí)喇勋，分支出現(xiàn)在軌跡中：一個(gè)發(fā)育譜系沿著一條路徑前進(jìn)缨该，而另一個(gè)譜系產(chǎn)生第二條路徑。Monocle包含分析這些分支事件的廣泛功能川背。Monocle提供了一個(gè)特殊的統(tǒng)計(jì)測(cè)試：分支表達(dá)式分析建模贰拿，或BEAM。
BEAM(Branched expression analysis modeling)是一種統(tǒng)計(jì)方法熄云，用于尋找以依賴于分支的方式調(diào)控的基因膨更。

plot_cell_trajectory(cds, color_by = "State")

結(jié)合前面color_by="Pseudotime"的圖，這個(gè)圖默認(rèn)的是最左邊是root(pseudotime 0)缴允，所以這個(gè)圖是status1是最早的荚守，在分叉2處分為status1和status7，status2又隨后在分叉3處分為status3和status4练般，status4隨后又在分叉1處分出status6和status5矗漾。

BEAM進(jìn)行統(tǒng)計(jì)分析

BEAM_res <- BEAM(cds[ordergene,], branch_point = 1, cores = 2) 
#這里用的是ordergene，也就是第六步dpFeature找出來的基因薄料。如果前面用的是seurat的marker基因敞贡，記得改成express_genes
#BEAM_res <- BEAM(cds, branch_point = 1, cores = 2) #對(duì)2829個(gè)基因進(jìn)行排序，運(yùn)行慢
BEAM_res <- BEAM_res[order(BEAM_res$qval),]
BEAM_res <- BEAM_res[,c("gene_short_name", "pval", "qval")]
head(BEAM_res)
#           gene_short_name         pval         qval
# CD79A               CD79A 2.757676e-73 7.782161e-70
# TCL1A               TCL1A 1.574889e-65 2.222168e-62
# IGLL5               IGLL5 2.356778e-64 2.216942e-61
# S100A9             S100A9 1.504319e-58 1.061297e-55
# S100A8             S100A8 6.028175e-57 3.402302e-54
# LINC00926       LINC00926 3.180527e-55 1.495908e-52
write.csv(BEAM_res, "BEAM_res.csv", row.names = F)
plot_genes_branched_heatmap(cds[row.names(subset(BEAM_res,
                                                  qval < 1e-4)),],
                            branch_point = 1, #繪制的是哪個(gè)分支
                            num_clusters = 4, #分成幾個(gè)cluster摄职，根據(jù)需要調(diào)整
                            cores = 1,
                            use_gene_short_name = T,
                            show_rownames = T)#有632個(gè)gene誊役，太多了

該熱圖顯示的是同一時(shí)間點(diǎn)兩個(gè)譜系的變化，熱圖的列是偽時(shí)間的點(diǎn)琳钉，行是基因势木。這張圖最上面的條條，灰色的代表分叉前歌懒，左邊紅色代表左邊這個(gè)cell fate啦桌，右邊藍(lán)色代表右邊這個(gè)cell fate，從熱圖中間往右讀及皂，是偽時(shí)間的一個(gè)譜系甫男，往左是另一個(gè)譜系⊙樯眨基因是被按照等級(jí)聚類的板驳，需要結(jié)合生物學(xué)知識(shí)來進(jìn)行解讀。

????因?yàn)榍懊嬖O(shè)置的branch_point = 1碍拆，根據(jù)按照State繪制的trajectory圖若治，在1這個(gè)分枝上分出的是5和6這兩個(gè)status慨蓝，所以上圖比較的是state5和state6這兩個(gè)branch，而且上圖中的pre-branch使用的是4, 2, 7, 1這四個(gè)status（BEAM tries to traverse backward from the cell on the branch point all the way back to the root cell (the cell with pseudotime 0) and use all those cells as the the pre-branch）端幼。而且礼烈，上圖左邊的Cell fate 1指的是state5，Cell fate 2指的是state6（Cell fate 1 corresponds to the state with small id while cell fate 2 corresponds to sate with bigger id ）婆跑。

#選前100個(gè)基因可視化
BEAM_genes <- top_n(BEAM_res, n = 100, desc(qval)) %>% pull(gene_short_name) %>% as.character()
p <- plot_genes_branched_heatmap(cds[BEAM_genes,],  branch_point = 1, 
                                 num_clusters = 3, show_rownames = T, return_heatmap = T)
ggsave("BEAM_heatmap.pdf", p$ph_res, width = 6.5, height = 10)

這個(gè)圖的解讀：
· How to map cell fate to branches?
· monocle2 擬時(shí)間分支點(diǎn)分析結(jié)果解讀

#顯著差異基因(top100)按熱圖結(jié)果排序并保存
##如果要所有的差異基因此熬，就把前面所632個(gè)基因的熱圖存為p
hp.genes <- p$ph_res$tree_row$labels[p$ph_res$tree_row$order]
BEAM_sig <- BEAM_res[hp.genes, c("gene_short_name", "pval", "qval")]
write.csv(BEAM_sig, "BEAM_sig.csv", row.names = F)

選擇上面熱圖中或顯著差異基因中感興趣的基因進(jìn)行可視化

head(BEAM_sig)
#        gene_short_name         pval         qval
# STX11            STX11 2.928177e-20 1.180474e-18
# CEBPD            CEBPD 2.137741e-22 1.058369e-20
# TYROBP          TYROBP 3.157539e-54 1.272939e-51
# FCER1G          FCER1G 1.094272e-45 2.573364e-43
# SRGN              SRGN 7.564014e-24 4.356255e-22
# CTSD              CTSD 2.666673e-17 7.601364e-16
genes <- row.names(subset(fData(cds),
                          gene_short_name %in% c( "SITX11", "CEBPD", "TYROBP")))

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

可以看到，在pre-branch細(xì)胞向 5/6這兩個(gè)分叉分化時(shí)滑进，這兩個(gè)基因都是在state5中高表達(dá)的犀忱，也就是state5的分支相關(guān)基因。

??注意：不同分組間的細(xì)胞盡量不要放在一起做軌跡分析扶关，同一組的生物學(xué)重復(fù)可以一起分析阴汇。明顯沒有生物學(xué)相關(guān)性的細(xì)胞也不要放在一起做軌跡分析。

10. 文獻(xiàn)里的monocle

有時(shí)候細(xì)胞群不會(huì)完美的按照分叉排列节槐。如下鲫寄，只要一類細(xì)胞占某一個(gè)分叉細(xì)胞的大部分也可以。

Single-cell reconstruction of the adult human heart during heart failure and recovery reveals the cellular landscape underlying cardiac function, Extended Data Fig. 8a-e

分析過程中使用的主要函數(shù)匯總