差異基因表達(dá)分析是RNA-Seq實(shí)驗(yàn)中的常見任務(wù)钱骂，Monocle2可以幫助我們找到在不同組細(xì)胞之間差異表達(dá)的基因，并評估這些基因差異表達(dá)變化的統(tǒng)計學(xué)意義掠归。
Basic Differential Analysis
這里吼砂，我們查看成肌細(xì)胞（myogenesis）數(shù)據(jù)在不同培養(yǎng)條件下的差異表達(dá)基因厉颤。
# 選擇myogenesis相關(guān)的基因集
marker_genes <- row.names(subset(fData(HSMM_myo),
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")))
# 使用differentialGeneTest函數(shù)對myogenesis基因集中的基因進(jìn)行差異表達(dá)分析
diff_test_res <- differentialGeneTest(HSMM_myo[marker_genes,],
fullModelFormulaStr = "~Media")
# Select genes that are significant at an FDR < 10%
# 篩選顯著的差異表達(dá)基因
sig_genes <- subset(diff_test_res, qval < 0.1)
# 對marker基因的表達(dá)進(jìn)行可視化
MYOG_ID1 <- HSMM_myo[row.names(subset(fData(HSMM_myo),
gene_short_name %in% c("MYOG", "CCNB2"))),]
plot_genes_jitter(MYOG_ID1,
grouping = "Media",
ncol= 2)
Finding Genes that Distinguish Cell Type or State
接下來绩鸣，我們根據(jù)幾個關(guān)鍵的marker基因?qū)⒊杉〖?xì)胞(myoblasts)與成纖維細(xì)胞(fibroblasts)區(qū)分開來怀大，使用Monocle2查看成纖維細(xì)胞和成肌細(xì)胞之間差異表達(dá)的基因。
# 選擇marker基因
to_be_tested <- row.names(subset(fData(HSMM),
gene_short_name %in% c("UBC", "NCAM1", "ANPEP")))
cds_subset <- HSMM[to_be_tested,]
# 使用differentialGeneTest函數(shù)對不同細(xì)胞類型進(jìn)行差異表達(dá)分析
diff_test_res <- differentialGeneTest(cds_subset,
fullModelFormulaStr = "~CellType")
diff_test_res[,c("gene_short_name", "pval", "qval")]
# 基因表達(dá)可視化
plot_genes_jitter(cds_subset,
grouping = "CellType",
color_by = "CellType",
nrow= 1,
ncol = NULL,
plot_trend = TRUE)
Finding Genes that Change as a Function of Pseudotime
Monocle2還可以基于擬時序分析的結(jié)果查看隨著細(xì)胞進(jìn)展而變化的基因呀闻。同樣化借，我們需要指定用于差異分析的模型。這個模型比我們之前用來查看CellType之間差異的模型要復(fù)雜一些总珠。Monocle2為每個細(xì)胞分配一個“偽時間”值屏鳍，該值記錄了其在實(shí)驗(yàn)過程中的進(jìn)度勘纯。Monocle2使用VGAM軟件包將基因的表達(dá)水平建模為偽時間的平滑非線性函數(shù)局服。
to_be_tested <- row.names(subset(fData(HSMM),
gene_short_name %in% c("MYH3", "MEF2C", "CCNB2", "TNNT1")))
cds_subset <- HSMM_myo[to_be_tested,]
# 使用differentialGeneTest函數(shù)進(jìn)行差異表達(dá)分析
diff_test_res <- differentialGeneTest(cds_subset,
fullModelFormulaStr = "~sm.ns(Pseudotime)")
# 使用plot_genes_in_pseudotime函數(shù)對基因的表達(dá)水平進(jìn)行可視化
plot_genes_in_pseudotime(cds_subset,
color_by = "Hours")
Clustering Genes by Pseudotemporal Expression Pattern
Monocle2可以對具有相似表達(dá)趨勢的基因進(jìn)行分組，并對其進(jìn)行可視化展示驳遵。Monocle2使用plot_pseudotime_heatmap函數(shù)基于CellDataSet對象（通常僅包含重要基因的子集）生成平滑的表達(dá)曲線淫奔，類似于plot_genes_in_pseudotime函數(shù)。然后堤结，將這些基因進(jìn)行聚類唆迁，并使用pheatmap包對其進(jìn)行熱圖繪制。
# 差異表達(dá)分析
diff_test_res <- differentialGeneTest(HSMM_myo[marker_genes,],
fullModelFormulaStr = "~sm.ns(Pseudotime)")
# 選擇顯著的差異基因
sig_gene_names <- row.names(subset(diff_test_res, qval < 0.1))
# 使用plot_pseudotime_heatmap函數(shù)進(jìn)行聚類分組可視化
plot_pseudotime_heatmap(HSMM_myo[sig_gene_names,],
num_clusters = 3,
cores = 1,
show_rownames = T)
Multi-Factorial Differential Expression Analysis
Monocle2還可以在存在多種因素的情況下執(zhí)行差異表達(dá)分析竞穷，這可以幫助我們減去某些因素以查看其他因素的影響唐责。為此，我們必須同時指定完整模型和簡化模型瘾带。這里鼠哥，我們指定完整模型捕獲CellType和Hours的效果，而簡化模型僅捕獲Hours的影響看政。
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 + Hours",
reducedModelFormulaStr = "~Hours")
plot_genes_jitter(cds_subset,
grouping = "Hours",
color_by = "CellType",
plot_trend = TRUE) +
facet_wrap( ~ feature_label, scales= "free_y")
Analyzing Branches in Single-Cell Trajectories
通常朴恳，細(xì)胞的發(fā)育軌跡會存在不同的分支，出現(xiàn)分支是因?yàn)榧?xì)胞執(zhí)行其他不同的功能允蚣。當(dāng)細(xì)胞做出命運(yùn)選擇時于颖，其發(fā)育軌跡中就會出現(xiàn)分支：其中，一個發(fā)育世系沿著一條路徑前進(jìn)嚷兔，而另一個世系沿著第二條路徑前進(jìn)森渐。Monocle2提供了一種特殊的統(tǒng)計檢驗(yàn)方法：branched expression analysis modeling（BEAM）做入，可以對不同的分支事件進(jìn)行分析。
lung <- load_lung()
plot_cell_trajectory(lung, color_by = "Time")
BEAM方法接受一個已通過orderCells函數(shù)排序的CellDataSet以及發(fā)育軌跡中分支點(diǎn)的名稱作為輸入同衣，返回每個基因的顯著性得分表母蛛。得分顯著的基因在其表達(dá)中被認(rèn)為是分支依賴性的。
# 使用BEAM函數(shù)進(jìn)行分支表達(dá)建模分析
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")]
我們可以使用plot_genes_branched_heatmap函數(shù)繪制特殊類型的熱圖乳怎，可視化所有明顯依賴于分支的基因的表達(dá)水平的變化彩郊。此熱圖同時顯示了兩個譜系中基因表達(dá)水平的變化，它還要求選擇要檢查的分支點(diǎn)蚪缀。其中秫逝，列是偽時間中的點(diǎn)，行是基因询枚，偽時間的起始點(diǎn)位于熱圖的中間违帆。這些基因按層次結(jié)構(gòu)聚類，可以可視化具有相似譜系依賴性表達(dá)模式的基因模塊金蜀。
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)
參考來源：http://cole-trapnell-lab.github.io/monocle-release/docs/
