差異基因表達(dá)分析是RNA-Seq實(shí)驗(yàn)中的常見(jiàn)任務(wù)，Monocle2可以幫助我們找到在不同組細(xì)胞之間差異表達(dá)的基因，并評(píng)估這些基因差異表達(dá)變化的統(tǒng)計(jì)學(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ù)對(duì)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)
# 對(duì)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
接下來(lái)，我們根據(jù)幾個(gè)關(guān)鍵的marker基因?qū)⒊杉〖?xì)胞(myoblasts)與成纖維細(xì)胞(fibroblasts)區(qū)分開(kāi)來(lái)，使用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ù)對(duì)不同細(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還可以基于擬時(shí)序分析的結(jié)果查看隨著細(xì)胞進(jìn)展而變化的基因。同樣，我們需要指定用于差異分析的模型。這個(gè)模型比我們之前用來(lái)查看CellType之間差異的模型要復(fù)雜一些。Monocle2為每個(gè)細(xì)胞分配一個(gè)“偽時(shí)間”值，該值記錄了其在實(shí)驗(yàn)過(guò)程中的進(jìn)度。Monocle2使用VGAM軟件包將基因的表達(dá)水平建模為偽時(shí)間的平滑非線性函數(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ù)對(duì)基因的表達(dá)水平進(jìn)行可視化
plot_genes_in_pseudotime(cds_subset,
color_by = "Hours")
Clustering Genes by Pseudotemporal Expression Pattern
Monocle2可以對(duì)具有相似表達(dá)趨勢(shì)的基因進(jìn)行分組，并對(duì)其進(jìn)行可視化展示。Monocle2使用plot_pseudotime_heatmap函數(shù)基于CellDataSet對(duì)象（通常僅包含重要基因的子集）生成平滑的表達(dá)曲線，類似于plot_genes_in_pseudotime函數(shù)。然后，將這些基因進(jìn)行聚類，并使用pheatmap包對(duì)其進(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á)分析，這可以幫助我們減去某些因素以查看其他因素的影響。為此，我們必須同時(shí)指定完整模型和簡(jiǎn)化模型。這里，我們指定完整模型捕獲CellType和Hours的效果，而簡(jiǎn)化模型僅捕獲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ā)育軌跡會(huì)存在不同的分支，出現(xiàn)分支是因?yàn)榧?xì)胞執(zhí)行其他不同的功能。當(dāng)細(xì)胞做出命運(yùn)選擇時(shí)，其發(fā)育軌跡中就會(huì)出現(xiàn)分支：其中，一個(gè)發(fā)育世系沿著一條路徑前進(jìn)，而另一個(gè)世系沿著第二條路徑前進(jìn)。Monocle2提供了一種特殊的統(tǒng)計(jì)檢驗(yàn)方法：branched expression analysis modeling（BEAM），可以對(duì)不同的分支事件進(jìn)行分析。
lung <- load_lung()
plot_cell_trajectory(lung, color_by = "Time")
BEAM方法接受一個(gè)已通過(guò)orderCells函數(shù)排序的CellDataSet以及發(fā)育軌跡中分支點(diǎn)的名稱作為輸入，返回每個(gè)基因的顯著性得分表。得分顯著的基因在其表達(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á)水平的變化。此熱圖同時(shí)顯示了兩個(gè)譜系中基因表達(dá)水平的變化，它還要求選擇要檢查的分支點(diǎn)。其中，列是偽時(shí)間中的點(diǎn)，行是基因，偽時(shí)間的起始點(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)
參考來(lái)源：http://cole-trapnell-lab.github.io/monocle-release/docs/
