其實在前面也提到過Pseudotime分析忆肾，擬時間序列分析（Pseudotime分析）是通過構(gòu)建細胞間的變化軌跡來重塑細胞隨著時間的變化過程。從具體的分類分析和復(fù)雜程度來說我抠，可以分為細胞軌跡分析和細胞譜系分析泳叠。

其中，Monocle2是做單細胞擬時分析最有名的R包某弦。對比還在持續(xù)開發(fā)中的Monocle3來說，Monocle2更穩(wěn)定且更傾向于半監(jiān)督的分析模式而克，更適合針對感興趣的細胞亞群做個性化分析靶壮。

====Monocle2的使用=====

官方教程：http://cole-trapnell-lab.github.io/monocle-release/docs/

安裝

BiocManager::install("monocle")

2.?創(chuàng)建CellDataSet

創(chuàng)建所需要的CellDataSet一般有以下幾種方法。

2.1 將Seurat object中數(shù)據(jù)提取來創(chuàng)建

數(shù)據(jù)準備：暫時用pbmc3k數(shù)據(jù)（軌跡分析的前提是待分析的細胞有緊密的發(fā)育關(guān)系员萍，PBMC細胞不是很好的的示例數(shù)據(jù)腾降，在此僅作為演示。）

由于該數(shù)據(jù)集沒有對細胞類型進行注釋碎绎，因此我們參考seurat標準流程對這個數(shù)據(jù)集進行注釋螃壤。

library(dplyr)

library(Seurat)

library(patchwork)?

pbmc.data <- Read10X(data.dir = "filtered_gene_bc_matrices/hg19/")

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

levels(pbmc)

saveRDS(pbmc, file = "pbmc.rds")

運行完上述代碼，得到注釋好的數(shù)據(jù)集筋帖。

從Seurat對象中提取構(gòu)建CDS對象所需要的3個輸入文件：表達矩陣信息奸晴、基因信息和表型信息。

monocle輸入的是count矩陣（不建議使用data矩陣）

library(monocle)

pbmc <- readRDS("pbmc.rds")?

expr_matrix <- as(as.matrix(pbmc@assays$RNA@counts), 'sparseMatrix')

p_data <- pbmc@meta.data?

p_data$celltype <- pbmc@active.ident??

f_data <- data.frame(gene_short_name = row.names(pbmc),row.names = row.names(pbmc))

pd <- new('AnnotatedDataFrame', data = p_data)?

fd <- new('AnnotatedDataFrame', data = f_data)

cds <- newCellDataSet(expr_matrix,? phenoData = pd,featureData = fd,lowerDetectionLimit = 0.5,expressionFamily = negbinomial.size())

注：一般來說幕随，F(xiàn)PKM/TPM值通常是對數(shù)正態(tài)分布的蚁滋，而UMIs或讀計數(shù)使用負二項分布能更好地建模。要處理計數(shù)數(shù)據(jù)赘淮，需要將負二項分布指定為newCellDataSet的expressionFamily參數(shù)：

從參數(shù)的釋義來看辕录，negbinomial.size()和negbinomial()：輸入的表達矩陣為UMI,一般適用于10x的數(shù)據(jù)；negbinomial()的結(jié)果更準確梢卸，但是計算更耗時走诞；一般建議采用negbinomial.size()。稀疏矩陣用negbinomial.size()

tobit():適用于輸入的表達矩陣為FPKM或者TPM, 構(gòu)建monocle2的class時會自動進行l(wèi)og化計算蛤高。

gaussianff():輸入為log化后的FPKM或者TPM蚣旱。(目前在單細胞數(shù)據(jù)中碑幅，F(xiàn)PKM已不多用，smart-seq2平臺數(shù)據(jù)一般采用TPM)

2.2 直接讀取表達矩陣來創(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對象直接轉(zhuǎn)化為CellDataSet對象

importCDS(pbmc)

3.?估計size factor和離散度

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

cds <- estimateSizeFactors(cds)

cds <- estimateDispersions(cds)

注：與seurat把標準化后的表達矩陣保存在對象中不同异吻，monocle只保存一些中間結(jié)果在對象中裹赴，需要用時再用這些中間結(jié)果轉(zhuǎn)化。經(jīng)過上面三個函數(shù)的計算诀浪，mycds對象中多了SizeFactors棋返、Dipersions、num_cells_expressed和num_genes_expressed等信息雷猪。

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

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

因為Seurat已經(jīng)完成細胞過濾嚷辅，此步可省略。

但由于Seurat是通過基因表達量來對細胞進行過濾距误。因此簸搞，我們可以通過下面的代碼用表達某基因的細胞的數(shù)目對基因進行過濾，從而得到后續(xù)操作需要的基因准潭。

cds <- detectGenes(cds, min_expr = 0.1)? //將會在fData(cds)中添加一列num_cells_expressed

print(head(fData(cds)))

expressed_genes <- row.names(subset(fData(cds),num_cells_expressed >= 10))? //過濾掉在小于10個細胞中表達的基因

5. 細胞分類

Monocle官網(wǎng)教程提供了4個分類方法：

Classifying cells by type

Clustering cells without marker genes

Clustering cells using marker genes

Imputing cell type

但是：很多博主建議先將細胞注釋好再進行Monocle分析趁俊，不建議使用monocle做細胞分類。

6. 軌跡定義基因選擇及可視化和構(gòu)建軌跡

一般來說刑然，主要流程如下：

???Step 1: choosing genes that define progress

???Step 2: reducing the dimensionality of the data

???Step 3: ordering the cells in pseudotime

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

推斷單細胞軌跡是一個機器學(xué)習(xí)問題寺擂。在單細胞RNA-Seq中，低水平表達的基因通常非常嘈雜泼掠，但有些基因包含有關(guān)細胞狀態(tài)的重要信息怔软。因此，軌跡推斷的第一步就是選擇Monocle將用作機器學(xué)習(xí)方法輸入的基因择镇。這叫做特征選擇挡逼，它對軌跡的形狀有很大的影響。

Monocle主要基于關(guān)鍵基因的表達模式腻豌，通過學(xué)習(xí)每個細胞必須經(jīng)歷的基因表達變化的序列家坎，根據(jù)擬時間值中對單個細胞進行排序嘱能，模擬出時間發(fā)育過程的動態(tài)變化。而這個排序技術(shù)表現(xiàn)是一種在低維空間排布高維數(shù)據(jù)的降維技術(shù)虱疏。

Monocle提供了多種工具來選擇基因惹骂，這些基因?qū)a(chǎn)生一個健壯、準確和具有生物學(xué)意義的軌跡做瞪。你可以使用這些工具來執(zhí)行一個完全“無監(jiān)督”的分析析苫，在這個分析中，Monocle不知道你認為哪個基因是重要的穿扳●媒模或者，你可以利用marker gene知識來定義生物學(xué)進展矛物，從而形成Monocle的軌跡茫死。我們認為這種模式是“半監(jiān)督”的分析。

Monocle官網(wǎng)教程提供了4個選擇方法：

??? 選擇發(fā)育差異表達基因

??? 選擇clusters差異表達基因

??? 選擇離散程度高的基因

??? 自定義發(fā)育marker基因

前三種都是無監(jiān)督分析方法履羞，細胞發(fā)育軌跡生成完全不受人工干預(yù)峦萎；最后一種是半監(jiān)督分析方法，可以使用先驗知識輔助分析忆首。

//使用seurat選擇的高變基因

express_genes <- VariableFeatures(pbmc)

cds<- setOrderingFilter(cds, express_genes)

plot_ordering_genes(cds)

//使用clusters差異表達基因

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)

但是爱榔，理想情況下，我們希望盡可能少地使用生物學(xué)的先驗知識糙及。我們希望從數(shù)據(jù)中發(fā)現(xiàn)重要的排序基因详幽，而不是依賴于文獻，因為這可能會在排序中引入偏見浸锨。我們將從一種更簡單的方法開始唇聘，但是我們通常推薦一種更復(fù)雜的方法，稱為“dpFeature”柱搜。

//這一步輸入的expressed_genes來自于步驟4迟郎。

//后續(xù)分析使用的是該方法

//也可輸入seurat篩選出的高變基因：expressed_genes<- VariableFeatures(pbmc)

diff <-differentialGeneTest(cds[expressed_genes,],fullModelFormulaStr="~cell_type",cores=1)

//注：~后面是表示對誰做差異分析的變量，理論上可以為p_data的任意列名

head(diff)

//差異表達基因作為軌跡構(gòu)建的基因,差異基因的選擇標準是qval<0.01,decreasing=F表示按數(shù)值增加排序

deg <- subset(diff, qval < 0.01)? ? ? ? ? ? #選出2829個基因

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對象，后續(xù)的一系列操作都要依賴于這個list健爬。

#setOrderingFilter之后控乾，這些基因被儲存在cds@featureData@data[["use_for_ordering"]]，可以通過table(cds@featureData@data[["use_for_ordering"]])查看

plot_ordering_genes(cds)

注：圖中黑色的點表示用來構(gòu)建軌跡的差異基因浑劳，灰色表示背景基因阱持。紅色的線是根據(jù)第2步計算的基因表達大小和離散度分布的趨勢(可以看到，找到的基因?qū)儆陔x散度比較高的基因)

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

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

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

Step 2: 降維

一旦細胞有序排列魔熏，我們就可以在降維空間中可視化軌跡衷咽。所以首先選擇用于細胞排序的基因鸽扁，然后使用反向圖嵌入(DDRTree)算法對數(shù)據(jù)進行降維。

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

Step 3: 擬時間軸軌跡構(gòu)建和在擬時間內(nèi)排列細胞

將表達數(shù)據(jù)投射到更低的維度空間镶骗，通過機器學(xué)習(xí)描述細胞如何從一種狀態(tài)過渡到另一種狀態(tài)的軌跡桶现。假設(shè)軌跡具有樹狀結(jié)構(gòu)，一端是“根”鼎姊，另一端是“葉”骡和。盡可能地將最佳樹與數(shù)據(jù)匹配起來。這項任務(wù)被稱為“歧管學(xué)習(xí)”相寇，在生物過程的開始階段慰于，細胞從根部開始，沿著主干前進唤衫，直到到達第一個分支（如果有的話）婆赠。然后，細胞必須選擇一條路徑佳励，沿著樹走得越來越遠休里，直到到達一片葉子。一個細胞的偽時間值是它回到根的距離赃承。

根據(jù)order gene的表達趨勢妙黍，將細胞排序并完成軌跡構(gòu)建。

cds <- orderCells(cds)

//使用root_state參數(shù)可以設(shè)置擬時間軸的根瞧剖，如下面的擬時間著色圖中可以看出拭嫁，左邊是根。根據(jù)state圖可以看出筒繁，根是State1噩凹，若要想把另一端設(shè)為根，可以按如下操作毡咏。

#cds <- orderCells(cds, root_state = 5)#把State5設(shè)成擬時間軸的起始點

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

1. 以pseudotime值上色 (Pseudotime是monocle2基于細胞基因表達信息計算的概率，表示時間的先后逮刨。)

plot_cell_trajectory(cds,color_by="Pseudotime",size=1,show_backbone=TRUE)

2. 以細胞類型上色

plot_cell_trajectory(cds,color_by="cell_type",size=1,show_backbone=TRUE)

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

plot_cell_trajectory(cds, color_by = "State", size=1,show_backbone=TRUE)

4.?按照seurat分群排序細胞

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

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

同樣滴修己，我們也可以使用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

我們同樣也可以用樹形圖來展示：

plot_complex_cell_trajectory(cds, x = 1, y = 2,color_by = "celltype")

還可以畫沿時間軸的細胞密度圖

library(ggpubr)

df <- pData(cds)?

ggplot(df, aes(Pseudotime, colour = cell_type, fill=cell_type)) +

? geom_density(bw=0.5,size=1,alpha = 0.5)+theme_classic2()

手動進行顏色設(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ù)分析）

比如對State7的細胞感興趣：

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. 指定基因的可視化

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

//選擇前4個top基因并將其對象取出

keygenes <- head(ordergene,4)

cds_subset <- cds[keygenes,]

//可視化：以state/celltype/pseudotime進行

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

p1|p2|p3

//也可以自己指定基因

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

p1|p2|p3

//擬時序展示單個基因表達量

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. 尋找擬時相關(guān)的基因（擬時差異基因）

使用回歸算法 （注意：不要使用多核運算睬愤，經(jīng)常會出現(xiàn)警告）

Monocle的主要工作是通過生物過程（如細胞分化）將細胞按順序排列片仿，而不知道要提前查看哪些基因。一旦這樣做了尤辱，你就可以分析細胞砂豌，找到隨著細胞進展而變化的基因厢岂。

官方給出的差異分析有三大方法：

1、Basic DifferentialAnalysis

2阳距、Finding Genes thatDistinguish Cell Type or State

3塔粒、Finding Genes thatChange as a Function of Pseudotime

我們重點關(guān)注第三個：根據(jù)偽時間功能尋找差異基因

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

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

//這里是把排序基因（ordergene）提取出來做回歸分析，來找它們是否跟擬時間有顯著的關(guān)系,如果不設(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()

plot_pseudotime_heatmap(cds[Time_genes,],num_clusters=4, show_rownames=T, return_heatmap=T)

注：圖中橫軸是擬時間，cluster1的基因是在擬時排序起點高表達的基因馍管，cluster2的基因則是在擬時排序的重點高表達的郭赐。cluster數(shù)的多少是由plot_pseudotime_heatmap函數(shù)中的num_clusters參數(shù)定義的。

plot_pseudotime_heatmap函數(shù)可以來可視化所有monocle的假時間依賴性基因咽斧。plot_pseudotime_heatmap采用CellDataSet對象（通常只包含重要基因的子集）堪置，并生成平滑的表達曲線，非常類似于plot_genes_in_pseudotime张惹。然后它對這些基因進行聚類舀锨，并使用pheatmap軟件包進行繪圖。繪出的熱圖可以讓我們觀測到假時間依賴性基因中的不同基因模塊在不同的時間內(nèi)共同變化宛逗，能比較好的回答時間序列基因表達中“哪些基因遵循相似的動力學(xué)趨勢”這一常見問題坎匿。

顯著差異基因按熱圖結(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)

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. 單細胞軌跡的“分支”分析

上一步尋找擬時相關(guān)的基因是全局的，找擬時起點和終點相關(guān)的基因雷激，這一步則是尋找和分叉點相關(guān)的基因替蔬。

單細胞軌跡常常包括分支。這些分支的產(chǎn)生是因為細胞執(zhí)行不同的基因表達程序屎暇。在發(fā)育過程中承桥，當(dāng)細胞做出命運選擇時，分支出現(xiàn)在軌跡中：一個發(fā)育譜系沿著一條路徑前進根悼，而另一個譜系產(chǎn)生第二條路徑凶异。Monocle包含分析這些分支事件的廣泛功能。Monocle提供了一個特殊的統(tǒng)計測試：分支表達式分析建模挤巡，或BEAM剩彬。BEAM(Branched expression analysis modeling)是一種統(tǒng)計方法，用于尋找以依賴于分支的方式調(diào)控的基因矿卑。

plot_cell_trajectory(cds, color_by = "State")

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

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) #對所有基因進行排序轻黑，運行慢

BEAM_res <-BEAM_res[order(BEAM_res$qval),]

BEAM_res <-BEAM_res[,c("gene_short_name", "pval", "qval")]

head(BEAM_res)

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, num_clusters = 4, cores = 1,use_gene_short_name= T,show_rownames = T)

該熱圖顯示的是同一時間點兩個譜系的變化糊肤，熱圖的列是偽時間的點，行是基因苔悦。這張圖最上面的條條轩褐，灰色的代表分叉前，左邊紅色代表左邊這個cell fate玖详，右邊藍色代表右邊這個cell fate把介，從熱圖中間往右讀，是偽時間的一個譜系蟋座，往左是另一個譜系拗踢。基因是被按照等級聚類的向臀，需要結(jié)合生物學(xué)知識來進行解讀巢墅。但是因為我們在參數(shù)設(shè)置的時候，選擇branch_point = 1是1這個分叉點的兩個cell的狀態(tài)券膀。

#也可以選前100個基因可視化

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)

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

//如果要所有的差異基因君纫，就把前面所有基因的熱圖存為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)

//選擇上面熱圖中或顯著差異基因中感興趣的基因進行可視化

head(BEAM_sig)

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)

從圖中可以看出，這2個基因都是在左邊的cell fate中高表達的芹彬。

10. 文獻里的monocle

有時候細胞群不會完美的按照分叉排列蓄髓。如下，只要一類細胞占某一個分叉細胞的大部分也可以舒帮。

本文使用 文章同步助手 同步