文章轉(zhuǎn)自：
作者：周運(yùn)來(lái)就是我
鏈接：http://www.reibang.com/p/66c387e1de3d
來(lái)源：簡(jiǎn)書
著作權(quán)歸作者所有。商業(yè)轉(zhuǎn)載請(qǐng)聯(lián)系作者獲得授權(quán)，非商業(yè)轉(zhuǎn)載請(qǐng)注明出處。
Seurat2monocle
library(monocle)
#<-importCDS(pbmc,import_all = TRUE)
#Load Seurat object
pbmc <- readRDS('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)
對(duì)于輸入文件终蒂，有3種類型的格式：
一種是單細(xì)胞運(yùn)行獲得的3個(gè)結(jié)果文件仆潮。
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.
- 準(zhǔn)備需要進(jìn)行擬時(shí)間分析數(shù)據(jù)的三個(gè)文件：
表達(dá)量文件：exprs：基因在所有細(xì)胞中的count值矩陣
2表型文件 phenoData:
該文件即為每個(gè)細(xì)胞的barcode信息。
3)featureData：
該文件為所有基因名信息肯污，有兩種格式，一種是Ensemble Id和Symbol Id的組合吨枉，一種均是Symbol Id,但要保證第二列一定是Symbol Id蹦渣。
官網(wǎng)給出許多其他建議以及讀入數(shù)據(jù)的注意事項(xiàng)，需要注意的是貌亭，如果是UMI數(shù)據(jù)柬唯，不應(yīng)該在創(chuàng)建CellDataSet之前自己對(duì)其進(jìn)行標(biāo)準(zhǔn)化。也不應(yīng)該試圖將UMI計(jì)數(shù)轉(zhuǎn)換為相對(duì)豐度(通過(guò)將其轉(zhuǎn)換為FPKM/TPM數(shù)據(jù))圃庭。Monocle將在內(nèi)部執(zhí)行所有需要的標(biāo)準(zhǔn)化步驟锄奢。自己將其正常化可能會(huì)破壞Monocle的一些關(guān)鍵步驟剧腻。
數(shù)據(jù)質(zhì)控
> 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提供了一個(gè)簡(jiǎn)單的系統(tǒng)拘央，可以根據(jù)您選擇的marker基因的表達(dá)來(lái)標(biāo)記細(xì)胞。您只需提供Monocle可以用來(lái)注釋每個(gè)單元格的一組函數(shù)书在。
聚類
第一步是決定用哪些基因來(lái)進(jìn)行聚類（特征選擇）灰伟。我們可以使用所有的基因，但是我們將包括很多沒(méi)有表達(dá)到足夠高水平來(lái)提供有意義的信號(hào)的基因儒旬，它們只會(huì)給系統(tǒng)增加噪音栏账。我們可以根據(jù)平均表達(dá)水平篩選基因，我們還可以選擇細(xì)胞間異常變異的基因义矛。這些基因往往對(duì)細(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)系對(duì)色散的期望。我們標(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的時(shí)候，只能畫出來(lái)4個(gè),為什么制轰？
HSMM <- clusterCells(HSMM, num_clusters = 5)
plot_cell_clusters(HSMM, 1, 2)
Monocle允許我們減去“無(wú)趣的”變異源的影響前计，以減少它們對(duì)集群的影響。您可以使用到clusterCells和其他幾個(gè)Monocle函數(shù)的residualmodelformula astr參數(shù)來(lái)實(shí)現(xiàn)這一點(diǎn)垃杖。此參數(shù)接受一個(gè)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來(lái)定義細(xì)胞以及外部導(dǎo)入細(xì)胞類型（Imputing cell type）（cell與type對(duì)應(yīng)文件）调俘，這里均不再展示伶棒。
構(gòu)建軌跡
在發(fā)育過(guò)程中旺垒，細(xì)胞會(huì)對(duì)刺激做出反應(yīng)，并在整個(gè)生命過(guò)程中肤无，從一種功能“狀態(tài)”轉(zhuǎn)變?yōu)榱硪环N功能“狀態(tài)”先蒋。不同狀態(tài)的細(xì)胞表達(dá)不同的基因，產(chǎn)生蛋白質(zhì)和代謝物的動(dòng)態(tài)重復(fù)序列宛渐，從而完成它們的工作竞漾。當(dāng)細(xì)胞在不同狀態(tài)間移動(dòng)時(shí)，會(huì)經(jīng)歷一個(gè)轉(zhuǎn)錄重組的過(guò)程窥翩，其中一些基因被沉默业岁，另一些基因被激活。這些瞬時(shí)狀態(tài)通常很難描述寇蚊，因?yàn)樵诟€(wěn)定的端點(diǎn)狀態(tài)之間純化細(xì)胞可能是困難的或不可能的笔时。
單細(xì)胞RNA-Seq可以使您在不需要純化的情況下看到這些狀態(tài)。然而仗岸，要做到這一點(diǎn)糊闽，我們必須確定每個(gè)細(xì)胞的可能狀態(tài)區(qū)間。
Monocle介紹了利用RNA-Seq進(jìn)行單細(xì)胞軌跡分析的策略爹梁。Monocle不是通過(guò)實(shí)驗(yàn)將細(xì)胞純化成離散狀態(tài)，而是使用一種算法來(lái)學(xué)習(xí)每個(gè)細(xì)胞必須經(jīng)歷的基因表達(dá)變化序列提澎，作為動(dòng)態(tài)生物學(xué)過(guò)程的一部分姚垃。一旦它了解了基因表達(dá)變化的整體“軌跡”，Monocle就可以將每個(gè)細(xì)胞置于軌跡中的適當(dāng)位置盼忌。
然后积糯，您可以使用Monocle的微分分析工具包來(lái)查找在軌跡過(guò)程中受到調(diào)控的基因放刨，如查找作為偽時(shí)間函數(shù)變化的基因明肮。如果這個(gè)過(guò)程有多個(gè)結(jié)果，Monocle將重建一個(gè)“分支”軌跡鹏秋。這些分支與細(xì)胞的“決策”相對(duì)應(yīng)跨嘉，Monocle提供了強(qiáng)大的工具來(lái)識(shí)別受它們影響的基因川慌，并參與這些基因的形成以及如何進(jìn)行分析分支。Monocle依靠一種叫做反向圖嵌入的機(jī)器學(xué)習(xí)技術(shù)來(lái)構(gòu)建單細(xì)胞軌跡祠乃。
在開始之前您也需要問(wèn)梦重，什么是擬時(shí)（偽時(shí)間）分析。
偽時(shí)間是衡量單個(gè)細(xì)胞在細(xì)胞分化等過(guò)程中取得了多大進(jìn)展的指標(biāo)亮瓷。在許多生物學(xué)過(guò)程中琴拧，細(xì)胞并不是完全同步的。在細(xì)胞分化等過(guò)程的單細(xì)胞表達(dá)研究中嘱支，捕獲的細(xì)胞在分化方面可能分布廣泛蚓胸。也就是說(shuō)挣饥，在同一時(shí)間捕獲的細(xì)胞群中，有些細(xì)胞可能已經(jīng)很長(zhǎng)時(shí)間了沛膳，而有些細(xì)胞甚至還沒(méi)有開始這個(gè)過(guò)程扔枫。
當(dāng)您想要了解在細(xì)胞從一種狀態(tài)轉(zhuǎn)換到另一種狀態(tài)時(shí)所發(fā)生的調(diào)節(jié)更改的順序時(shí)，這種異步性會(huì)產(chǎn)生主要問(wè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)生的問(wèn)題搓侄。
Monocle不是跟蹤表達(dá)式隨時(shí)間變化的函數(shù)，而是跟蹤沿軌跡變化的函數(shù)话速，我們稱之為偽時(shí)間讶踪。偽時(shí)間是一個(gè)抽象的分化單位：它只是一個(gè)cell到軌跡起點(diǎn)的距離，沿著最短路徑測(cè)量泊交。軌跡的總長(zhǎng)度是由細(xì)胞從起始狀態(tài)移動(dòng)到結(jié)束狀態(tài)所經(jīng)歷的總轉(zhuǎn)錄變化量來(lái)定義的乳讥。
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
基因的選擇
首先，我們必須決定我們將使用哪些基因來(lái)定義細(xì)胞在肌生成過(guò)程中的進(jìn)展廓俭。我們最終想要的是一組基因云石，能夠隨著我們研究過(guò)程的進(jìn)展而增加(或減少)表達(dá)。
理想情況下研乒，我們希望盡可能少地使用正在研究的系統(tǒng)生物學(xué)的先驗(yàn)知識(shí)汹忠。我們希望從數(shù)據(jù)中發(fā)現(xiàn)重要的排序基因，而不是依賴于文獻(xiàn)和教科書雹熬，因?yàn)檫@可能會(huì)在排序中引入偏見宽菜。我們將從一種更簡(jiǎ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ù)時(shí)間點(diǎn)的差異分析來(lái)選擇基因通常是非常有效的阵幸，但是如果我們沒(méi)有時(shí)間序列數(shù)據(jù)呢?如果細(xì)胞在我們的生物過(guò)程中是異步分化的(通常是這樣)，Monocle通逞渴溃可以從同時(shí)捕獲的單個(gè)種群重建它們的軌跡侨嘀。下面是兩種選擇基因的方法，它們完全不需要知道實(shí)驗(yàn)的設(shè)計(jì)捂襟。
一旦我們有了一個(gè)用于排序的基因id列表咬腕，我們就需要在HSMM對(duì)象中設(shè)置它們，因?yàn)榻酉聛?lái)的幾個(gè)函數(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)用樹的哪個(gè)軌跡中的哪個(gè)來(lái)調(diào)用“開始”，所以我們經(jīng)常不得不再次使用root_state參數(shù)調(diào)用orderCells來(lái)指定開始举反。首先懊直，我們繪制軌跡，這次按“狀態(tài)”給細(xì)胞上色火鼻。
···
plot_cell_trajectory(HSMM, color_by = "State")
···
plot_cell_trajectory(HSMM, color_by = "Pseudotime")
“狀態(tài)”只是單片樹的術(shù)語(yǔ)室囊，表示這段樹。下面的函數(shù)便于識(shí)別包含時(shí)間為零的大多數(shù)細(xì)胞的狀態(tài)魁索。然后我們可以把它傳遞給orderCells
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")
可以看出不同的排序條件擬時(shí)方向是不一樣的融撞。同時(shí)如果想看每個(gè)狀態(tài)下細(xì)胞的分化可以指定分面，當(dāng)然你也可以指定其他的分面方式粗蔚。
plot_cell_trajectory(HSMM_myo, color_by = "State") +
facet_wrap(~State, nrow = 1)
如果沒(méi)有timeseries尝偎，可能需要根據(jù)某些標(biāo)記基因的表達(dá)位置來(lái)設(shè)置根，使用您對(duì)系統(tǒng)的生物學(xué)知識(shí)鹏控。例如致扯，在這個(gè)實(shí)驗(yàn)中，高度增殖的祖細(xì)胞群產(chǎn)生兩種有絲分裂后的細(xì)胞当辐。所以根應(yīng)該有表達(dá)高水平增殖標(biāo)記的細(xì)胞抖僵。我們可以使用抖動(dòng)圖找出對(duì)應(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)
單個(gè)基因的時(shí)間變化（我可以隨意選擇基因而你不可以，你要選有意義的生活）
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á)的基因裆针，并評(píng)估這些變化的統(tǒng)計(jì)顯著性。這些比較要求您有一種方法將細(xì)胞收集到兩個(gè)或更多組中寺晌。這些組由每個(gè)CellDataSet的表現(xiàn)型數(shù)據(jù)表中的列定義。Monocle將評(píng)估不同細(xì)胞群中每個(gè)基因表達(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)
選擇對(duì)區(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中的差異分析過(guò)程非常靈活:您在測(cè)試中使用的模型公式可以包含pData表中作為列存在的任何項(xiàng)呻征，包括Monocle在其他分析步驟中添加的列。例如罢浇，如果您使用集群細(xì)胞陆赋，您可以使用集群作為模型公式來(lái)測(cè)試cluster之間不同的基因。
Finding Genes that Change as a Function of Pseudotime
Monocle的主要工作是將細(xì)胞按照生物過(guò)程(如細(xì)胞分化)的順序排列嚷闭，而不事先知道要查看哪些基因攒岛。一旦這樣做了，你就可以分析細(xì)胞胞锰，找到隨著細(xì)胞進(jìn)步而變化的基因灾锯。例如，你可以發(fā)現(xiàn)當(dāng)細(xì)胞“成熟”時(shí)嗅榕，基因顯著上調(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
在研究時(shí)間序列基因表達(dá)研究時(shí)吵聪，一個(gè)常見的問(wèn)題是:“哪些基因遵循相似的動(dòng)力學(xué)趨勢(shì)?”Monocle可以通過(guò)將具有相似趨勢(shì)的基因分組來(lái)幫助你回答這個(gè)問(wèn)題，這樣你就可以分析這些基因組兼雄，看看它們有什么共同之處吟逝。Monocle提供了一種方便的方法來(lái)可視化所有偽時(shí)間相關(guān)的基因。函數(shù)plot_pseudotime_heatmap接受一個(gè)CellDataSet對(duì)象(通常只包含重要基因的子集)赦肋，并生成與plot_genes_in_pseudotime類似的平滑表達(dá)曲線.然后块攒，它將這些基因聚類并使用pheatmap包繪制它們。這允許您可視化基因模塊佃乘，這些模塊在偽時(shí)間內(nèi)共同變化囱井。
注意下面的num_clusters 指的是基因可以聚成幾個(gè)類，而不是細(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可以在多個(gè)因素存在的情況下進(jìn)行差異分析琅绅，這可以幫助你減去一些因素來(lái)看到其他因素的影響。在下面的簡(jiǎn)單例子中鹅巍，Monocle測(cè)試了三個(gè)基因在成肌細(xì)胞和成纖維細(xì)胞之間的差異表達(dá)千扶，同時(shí)減去percent.mt的影響。為此骆捧，我們必須同時(shí)指定完整模型和簡(jiǎn)化模型澎羞。完整的模型同時(shí)捕捉細(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ā)育過(guò)程中枫攀，當(dāng)細(xì)胞做出命運(yùn)的選擇時(shí)括饶，分支就會(huì)出現(xiàn)在軌跡中:一個(gè)發(fā)育譜系沿著一條路徑前進(jìn)，而另一個(gè)譜系產(chǎn)生第二條路徑来涨。Monocle包含用于分析這些分支事件的廣泛功能图焰。
芭芭拉·特雷特琳和她的同事們?cè)谑返俜颉た说膶?shí)驗(yàn)室里進(jìn)行了一項(xiàng)實(shí)驗(yàn)，他們從正在發(fā)育的老鼠肺中獲取細(xì)胞蹦掐。他們?cè)诩?xì)胞發(fā)育的早期捕獲細(xì)胞技羔，之后當(dāng)肺包含兩種主要類型的上皮細(xì)胞(AT1和AT2)，以及即將決定成為AT1或AT2的細(xì)胞時(shí)卧抗。Monocle可以將這個(gè)過(guò)程重構(gòu)為一個(gè)分支軌跡藤滥，讓你可以非常詳細(xì)地分析決策點(diǎn)。下圖顯示了使用Monocle的一些數(shù)據(jù)重建的一般方法社裆。有一個(gè)單獨(dú)的分支拙绊，標(biāo)記為“1”。當(dāng)細(xì)胞從樹的左上方通過(guò)樹枝的早期發(fā)育階段通過(guò)時(shí)，哪些基因發(fā)生了變化?哪些基因在分支間有差異表達(dá)?為了回答這個(gè)問(wèn)題时呀，Monocle為您提供了一個(gè)特殊的統(tǒng)計(jì)測(cè)試:分支表達(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")]
使用一種特殊類型的熱圖谨娜，您可以可視化所有明顯依賴于分支的基因的變化航攒。這張熱圖同時(shí)顯示了兩種血統(tǒng)的變化。它還要求您選擇要檢查的分支點(diǎn)趴梢。列是偽時(shí)間中的點(diǎn)漠畜，行是基因，偽時(shí)間的開始在熱圖的中間坞靶。當(dāng)你從熱圖的中間讀到右邊的時(shí)候憔狞，你正在跟隨一個(gè)偽時(shí)間譜系。當(dāng)你讀到左邊時(shí)彰阴，另一個(gè)瘾敢。這些基因是分層聚類的，因此您可以可視化具有類似的依賴于序列的基因模塊.
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)
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)
