1.寫(xiě)在前面的話：

近年來(lái)吮成，由于細(xì)胞的異質(zhì)性及發(fā)育分化等相關(guān)的問(wèn)題越來(lái)越被研究者們所關(guān)注，單細(xì)胞轉(zhuǎn)錄組分析為研究異質(zhì)細(xì)胞群的復(fù)雜生物學(xué)過(guò)程提供了方法和工具巢价。每一個(gè)細(xì)胞進(jìn)行轉(zhuǎn)錄組測(cè)序時(shí)就是細(xì)胞發(fā)育過(guò)程中的快照倒得，單細(xì)胞擬時(shí)間分析軟件Monocle2是基于R語(yǔ)言的安裝包，其功能基于單細(xì)胞轉(zhuǎn)錄組的表達(dá)矩陣患蹂，通過(guò)無(wú)監(jiān)督學(xué)習(xí)（Reversed Graph Embedding算法）的方式將細(xì)胞置于發(fā)育軌跡的不同分支上或颊，從而模擬細(xì)胞群體生物學(xué)過(guò)程。也就是我們經(jīng)常說(shuō)的擬時(shí)序（pseudotime）分析传于，又稱細(xì)胞軌跡（cell trajectory）分析囱挑。通過(guò)擬時(shí)分析可以推斷出發(fā)育過(guò)程細(xì)胞的分化軌跡或細(xì)胞亞型的演化過(guò)程，在發(fā)育相關(guān)研究中使用頻率較高沼溜。
模擬細(xì)胞的分化軌跡的軟件平挑，最常用的軟件為Monocle2。該算法不僅能模擬細(xì)胞的發(fā)育軌跡系草，同時(shí)也能對(duì)細(xì)胞進(jìn)行聚類（t-SNE）通熄。通過(guò)聚類獲得不同狀態(tài)下的差異基因，分析影響分支形成的關(guān)鍵基因及其功能找都，對(duì)研究相關(guān)生物學(xué)問(wèn)題有指導(dǎo)性的作用唇辨。
Monocle2主要基于關(guān)鍵基因的表達(dá)模式，通過(guò)學(xué)習(xí)每個(gè)細(xì)胞必須經(jīng)歷的基因表達(dá)變化的序列能耻，根據(jù)擬時(shí)間值中對(duì)單個(gè)細(xì)胞進(jìn)行排序赏枚，模擬出時(shí)間發(fā)育過(guò)程的動(dòng)態(tài)變化。而這個(gè)排序技術(shù)表現(xiàn)是一種在低維空間排布高維數(shù)據(jù)的降維技術(shù)嚎京。（具體描敘請(qǐng)參考：https://www.meiwen.com.cn/subject/zlgsvqtx.html）嗡贺，目前monocle已經(jīng)升級(jí)至monocle3，但在結(jié)果分析上monocle3的可讀性不如monocle2鞍帝，因此本研究中主推monocle2這個(gè)版本诫睬。
在做擬時(shí)序分析之前，先要明白幾個(gè)研究目的：

1. 對(duì)什么類型的細(xì)胞群進(jìn)行擬時(shí)許分析帕涌？這類細(xì)胞群在不同的分化軌跡或細(xì)胞亞型的區(qū)別是否明顯摄凡？
2. monocle做擬時(shí)序分析時(shí)，可與seurat聚類結(jié)果進(jìn)行對(duì)接蚓曼，因此可對(duì)seurat的分群結(jié)果進(jìn)行定性的解釋亲澡。

3. 可分析比較組間的差異或者多組樣本間擬時(shí)間分析差別。
   分析流程圖如下所示：

   
   
   image.png

2. 軟件安裝

先保證服務(wù)器上裝有R 纫版，我這里安裝的版本是version 3.5.1床绪。

biocLite()
biocLite("monocle")
install.packages("devtools")
devtools::install_github("cole-trapnell-lab/monocle-release@develop")
biocLite(c("DDRTree", "pheatmap"))

3. 輸入文件

對(duì)于輸入文件，有3種類型的格式：
一種是單細(xì)胞運(yùn)行獲得的3個(gè)結(jié)果文件。格式如下表所示：
準(zhǔn)備需要進(jìn)行擬時(shí)間分析數(shù)據(jù)的三個(gè)文件：
1)表達(dá)量文件：exprs：基因在所有細(xì)胞中的count值矩陣癞己。
格式示例：

注意：該文件的表頭必須以“%%**\n%”的形式出現(xiàn)膀斋，否則就會(huì)出現(xiàn)報(bào)錯(cuò)。
2)表型文件 phenoData:
格式示例痹雅，該文件即為每個(gè)細(xì)胞的barcode信息仰担。：

3. featureData：
   該文件為所有基因名信息，有兩種格式绩社，一種是Ensemble Id和Symbol Id的組合忱屑，一種均是Symbol Id,但要保證第二列一定是Symbol Id蛋哭。
   格式示例：

   
   
   
   

   第二種輸入數(shù)據(jù)的格式是直接轉(zhuǎn)換seurat對(duì)象為CellDataSet格式送滞。
   第三種輸入數(shù)據(jù)的格式是只有基因和細(xì)胞構(gòu)成的表達(dá)量矩陣异雁。

4. 特征向量轉(zhuǎn)化為CellDataSet對(duì)象

方法一：3個(gè)特征文件轉(zhuǎn)換成CellDataSet格式阱扬。命令行如下所示：

gbm <- load_cellranger_matrix(my_dir) ###my_dir中包括有三個(gè)文件另玖。
HSMM <- newCellDataSet(exprs(gbm),phenoData = new("AnnotatedDataFrame", data = pData(gbm)),featureData = new("AnnotatedDataFrame", data = my_feat),lowerDetectionLimit = 0.5,expressionFamily = negbinomial.size())

方法二：直接通過(guò)seurat生成的rds文件的方式進(jìn)行轉(zhuǎn)換袁串。

由于seurat軟件包的升級(jí)刨沦，大量的函數(shù)命名發(fā)生了變化鸠真，對(duì)2和3這兩種類型的版本的操作也是不一樣悯仙。

spleen <-readRDS("seurat.analysis.rds") 

Seurat2.4的版本：

monocle_cds <-importCDS(spleen,import_all=TRUE) 

Seurat3.0及更高版本：

data <- as(as.matrix(spleen@assays$RNA@counts), 'sparseMatrix')

此處需要提供絕對(duì)的count值，即UMI矩陣吠卷。

pd <- new('AnnotatedDataFrame', data =data) 
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data) 
fd <- new('AnnotatedDataFrame', data = fData)    
HSMM <- newCellDataSet(data, 
?            phenoData = pd, 
?            featureData = fd, 
?            lowerDetectionLimit = 0.5, 
?            expressionFamily = negbinomial.size())

方法三：僅僅只有表達(dá)量數(shù)據(jù)的時(shí)候如何進(jìn)行處理锡垄？

如果你沒(méi)有3個(gè)標(biāo)準(zhǔn)格式的文件，也沒(méi)有生成rds文件祭隔，僅僅只有表達(dá)量矩陣货岭，具體原理和有rds文件一致，可按照如下方法進(jìn)行處理疾渴。

HSMM_expr_matrix <- read.csv("fpkm_matrix.csv",header=T)
pd <- new('AnnotatedDataFrame', data =data.frame(colnames(HSMM_expr_matrix)))
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
fd <- new('AnnotatedDataFrame', data = fData)    
HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix) , phenoData = pd, featureData = fd, lowerDetectionLimit = 0.5, expressionFamily = negbinomial.size())

5. 對(duì)monocle對(duì)象進(jìn)行歸一化千贯。

如果你采用的是第二種輸入數(shù)據(jù)進(jìn)行分析，即使你原先的Seurat對(duì)象已經(jīng)歸一化過(guò)了搞坝，官方推薦在monocle中重新進(jìn)行歸一化處理搔谴。

HSMM <- estimateSizeFactors(HSMM)   
HSMM <- estimateDispersions(HSMM)
HSMM <- detectGenes(HSMM, min_expr = 0.1) ###過(guò)濾掉低質(zhì)量的基因。

6.數(shù)據(jù)降維桩撮，特征基因的選擇敦第。

在進(jìn)行降維聚類之前，我們先獲得高表達(dá)的基因集作為每個(gè)聚類的用來(lái)order的Feature基因店量。當(dāng)然我們可以使用所有的基因芜果，但一些表達(dá)量特別低的基因提供的聚類信號(hào)往往會(huì)制造分析噪音，F(xiàn)eature基因的選擇性很多融师，一種是我們可以根據(jù)基因的平均表達(dá)水平來(lái)進(jìn)行篩選右钾，另外我們也可以選擇細(xì)胞間異常變異的基因。這些基因往往能較好地反映不同細(xì)胞的狀態(tài)。以平均表達(dá)量高于0.1來(lái)進(jìn)行篩選舀射。

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

繪制用于order基因和非order的平均表達(dá)量的分布點(diǎn)圖灭将。這里不做展示。

7.聚類

1.首先根據(jù)TSNE進(jìn)行降維處理

HSMM_myo <- reduceDimension(HSMM_myo, max_components= max_components, num_dim = 3, norm_method = 'log', reduction_method = 'tSNE',verbose = T)

對(duì)細(xì)胞進(jìn)行聚類后控，在Seurat中采用的是分辨率來(lái)確定cluster的數(shù)目庙曙。而monocle中可以直接指定聚類數(shù)目。主要指出的這里所聚類獲得的cluster數(shù)目比我們賦值的要少一個(gè)浩淘。即當(dāng)num_clusters=3時(shí)捌朴，你只獲得了2個(gè)cluster。具體解釋有些難懂张抄，在這里不做過(guò)多的解釋砂蔽。

2.獲得對(duì)細(xì)胞排序起到關(guān)鍵作用的基因集。

擬時(shí)間分析不僅是要對(duì)分析的細(xì)胞群進(jìn)行排序署惯，還要獲得覺(jué)得細(xì)胞排序的關(guān)鍵基因集左驾。這種基因集有兩種情況，在有先驗(yàn)知識(shí)的情況下极谊，我們可以從系統(tǒng)生物學(xué)的角度獲得一系列與細(xì)胞發(fā)育相關(guān)的基因集诡右，從而對(duì)細(xì)胞進(jìn)行排序，這種方式是最為保險(xiǎn)的轻猖，但局限性是對(duì)未知的發(fā)育情況不能進(jìn)行分析帆吻。另外一種情況就是使用無(wú)監(jiān)督聚類方式計(jì)算關(guān)鍵基因集。接下來(lái)我們采用differentialGeneTest方式獲得clustering_DEG_genes（與聚類相關(guān)的差異基因集）

clustering_DEG_genes <-differentialGeneTest(HSMM_myo[HSMM_expressed_genes,],fullModelFormulaStr = '~Cluster', cores = 4)

上一步過(guò)程是對(duì)所有的細(xì)胞進(jìn)行無(wú)監(jiān)督訓(xùn)練咙边，運(yùn)行時(shí)間與細(xì)胞數(shù)和基因數(shù)相關(guān)猜煮，一般會(huì)花很長(zhǎng)的時(shí)間“苄恚可以根據(jù)cores的數(shù)目進(jìn)行并行王带。
differentialGeneTest這個(gè)函數(shù)測(cè)試每個(gè)基因的差異表達(dá)，取決于偽時(shí)間或根據(jù)指定的其他協(xié)變量市殷。 “ differentialGeneTest”是Monocle的主要差異分析常規(guī)愕撰， 它接受一個(gè)CellDataSet和兩個(gè)模型公式作為輸入，指定由實(shí)現(xiàn)的廣義譜系模型“ VGAM”包被丧。 也就是說(shuō)我們可以根據(jù)指定’~cluster’或者擬時(shí)間值來(lái)獲得差異基因盟戏。差異基因的結(jié)果如下表所示：

在這個(gè)表格中，我們先看一下表頭對(duì)應(yīng)的關(guān)鍵列甥桂。第一列沒(méi)有列名柿究，為基因名稱。pval,qval 為差異基因的顯著檢驗(yàn)打分黄选。num_cells_expressed為表達(dá)這個(gè)基因的細(xì)胞數(shù)蝇摸。use_for_ordering該基因是否是作為排序的差異基因婶肩。

3. 細(xì)胞排序，一般情況下使用qval排在前1000個(gè)基因進(jìn)行排序貌夕。

HSMM_ordering_genes <-row.names(clustering_DEG_genes)[order(clustering_DEG_genes$qval)]
HSMM_myo <-setOrderingFilter(HSMM_myo,ordering_genes = HSMM_ordering_genes)
HSMM_myo <-reduceDimension(HSMM_myo, method = 'DDRTree')
HSMM_myo <-orderCells(HSMM_myo)

根據(jù)排序好細(xì)胞進(jìn)行結(jié)果可視化律歼。
命令行如下所示：

plot_cell_trajectory(HSMM_myo, color_by = "State",cell_size = 0.75) ####以state進(jìn)行著色
plot_cell_trajectory(HSMM_myo, color_by = "State",cell_size = 0.75)+facet_wrap(~State, nrow = Cluster_num) ##繪制state的分面圖
plot_cell_trajectory(HSMM_myo, color_by = "Pseudotime",cell_size = 0.75) ##根據(jù)擬時(shí)間值著色
plot_cell_trajectory(HSMM_myo, color_by = "Cluster",cell_size = 0.75) ##根據(jù)細(xì)胞聚類的進(jìn)行著色
plot_cell_trajectory(HSMM_myo, color_by = "Cluster",cell_size = 0.75)+facet_wrap(~Cluster, nrow =Cluster_num) ###繪制clusster的分面圖
plot_cell_trajectory(HSMM_myo, color_by = "seurat_clusters",cell_size = 0.75) +scale_color_manual(values=col) ##如果有Seurat生的rds文件的話，按照seurat中分的群進(jìn)行著色啡专，如果不想用ggplot的默認(rèn)色险毁，可以提供顏色列表col list。
plot_cell_trajectory(HSMM_myo, color_by = "seurat_clusters", cell_size = 0.75) + facet_wrap (~seurat_clusters, nrow =Cluster_num) ###按照seurat中分的群繪制分面圖们童。
plot_cell_trajectory(HSMM_myo, color_by = "stim",cell_size = 0.75) ###按照樣本進(jìn)行著色
plot_cell_trajectory(HSMM_myo, color_by = "stim",cell_size = 0.75)+facet_wrap(~stim, ncol = 2)
##繪制樣本著色的分面圖畔况。

上述部分結(jié)果如下圖所示（不包括分面圖）

上圖表示在主成分中的細(xì)胞聚類分布的tsne圖。不同顏色代表細(xì)胞群的不同細(xì)胞命運(yùn)的分支慧库。

上圖表示依據(jù)seurat的cluster ID映射到擬時(shí)間的排序上齐板。
接下來(lái)可視化展示細(xì)胞state相關(guān)的差異基因的表達(dá)量分布情況吵瞻，可以根據(jù)這些基因來(lái)確定細(xì)胞的發(fā)育方向，下圖僅展示qvalue值排在前6個(gè)基因甘磨。橫縱左邊意義見(jiàn)坐標(biāo)軸的描述橡羞。

8 分支點(diǎn)的差異分析

在monocle中的分析中發(fā)現(xiàn)宣羊，細(xì)胞群能從一個(gè)軌跡分叉成不同的分支點(diǎn)，在科研中汰蜘，我們會(huì)比較關(guān)注發(fā)生分支的原因是什么仇冯，即分析分支點(diǎn)之間的差異。Monocle采用分支表達(dá)式分析建模族操，主要是BEAM函數(shù)苛坚，可以將分叉過(guò)程重構(gòu)為一個(gè)分支軌跡，從而分析不同細(xì)胞命運(yùn)下的差異。
命令行如下：

my_pseudotime_de <- differentialGeneTest(HSMM_myo, fullModelFormulaStr= " ~sm.ns(Pseudotime)" , cores = 4) ###這里是對(duì)擬時(shí)間值來(lái)計(jì)算差異基因泼舱。依然花的時(shí)間比較長(zhǎng)等缀。

舉個(gè)例子：我們分析branch_point = 1這個(gè)分支處的細(xì)胞命名分叉是如何進(jìn)行的。

BEAM_res <- BEAM(HSMM_myo, branch_point = 1, cores = 4) ###選擇branch_point =1娇昙。

即下圖中所示尺迂。下圖只有1個(gè)分支點(diǎn)，即分析state1,2,3 這三個(gè)state的差異冒掌。

對(duì)影響分析的基因根據(jù)qvalue進(jìn)行排序枪狂。

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

繪制與分支相關(guān)的基因熱圖。

my_branched_heatmap <- plot_genes_branched_heatmap (HSMM_myo[row.names (subset (BEAM_res, qval < 1e-4)), ], branch_point = 1, num_clusters = 4, cores = 4, use_gene_short_name = TRUE,show_rownames = FALSE,return_heatmap = TRUE)

關(guān)鍵參數(shù)為subset (BEAM_res, qval < 1e-4))宋渔，挑選基因進(jìn)行熱圖繪制州疾，也可以設(shè)置成其他的閾值。
branch_point=1皇拣，分支點(diǎn)選為1严蓖。num_clusters=4，將基因根據(jù)表達(dá)相似性分成4個(gè)模塊氧急。結(jié)果如下圖所示：