原理簡(jiǎn)介

單細(xì)胞轉(zhuǎn)錄組是通過(guò)油包水保留一個(gè)細(xì)胞的整個(gè)轉(zhuǎn)錄組信息丹锹，且細(xì)胞分化是一個(gè)連續(xù)的過(guò)程，所以單細(xì)胞轉(zhuǎn)錄組測(cè)序能夠捕捉到不同分化狀態(tài)的細(xì)胞。而分化狀態(tài)的確定是根據(jù)不同細(xì)胞之間基因表達(dá)的連續(xù)變化精钮。所以可以通過(guò)軌跡分析分析出細(xì)胞的分化軌跡或者細(xì)胞亞型的演化過(guò)程。
Monocle分析能夠通過(guò)機(jī)器學(xué)習(xí)構(gòu)建單細(xì)胞軌跡剃斧，并通過(guò)差異分析過(guò)去軌跡過(guò)程中受調(diào)控的基因轨香，但是缺少對(duì)關(guān)鍵驅(qū)動(dòng)基因的下游分析。
GeneSwitches通過(guò)對(duì)擬時(shí)序結(jié)果進(jìn)一步分析幼东，挖掘?qū)?xì)胞分化有重要影響的開(kāi)關(guān)基因臂容。

代碼詳解

GeneSwitches分析首先對(duì)分化軌跡中的基因進(jìn)行二值化分析科雳，將基因的表達(dá)值化為1（打開(kāi)）或0（關(guān)閉）。為了實(shí)現(xiàn)這個(gè)目標(biāo)脓杉，程序?qū)蓚€(gè)符合高斯分布的模型擬合到每個(gè)基因的表達(dá)中糟秘，計(jì)算二值化的特定閾值

sce_p1 <- binarize_exp(sce_p1, ncores = 3)

或者可以通過(guò)所有基因的表達(dá)直方圖，尋找零分布和非0分布的中斷球散，以確定全局閾值尿赚。

### 檢查二值化的閾值
h <- hist(assays(sce_p1)$expdata, breaks = 200, plot = FALSE)
{plot(h, freq = FALSE, xlim = c(0,2), ylim = c(0,1), main = "Histogram of gene expression",
xlab = "Gene expression", col = "darkgoldenrod2", border = "grey")
abline(v=0.2, col="blue")}

image.png

選擇0.2作為閾值

bn_cutoff <- 0.2
sce_p1 <- binarize_exp(sce_p1, fix_cutoff = TRUE, binarize_cutoff = bn_cutoff)

然后對(duì)這些潛在的開(kāi)關(guān)基因進(jìn)行邏輯回歸分析和McFadden’s Pseudo R^2擬時(shí)間相關(guān)性分析。通過(guò)邏輯回歸分析推算出每個(gè)開(kāi)關(guān)基因的開(kāi)關(guān)時(shí)間（Switching Time Point）

sce_p1 <- find_switch_logistic_fastglm(sce_p1, downsample = TRUE, show_warning = FALSE)

下面我們對(duì)排序后的開(kāi)關(guān)基因進(jìn)行可視化蕉堰，首先會(huì)過(guò)濾在細(xì)胞中表達(dá)小于10%凌净，F(xiàn)DR>0.05,和McFadden's Pseudo R^2（<0.03）(擬合不佳)的基因。然后選擇前15的擬合最佳的基因進(jìn)行可視化屋讶。也可以保留包含表面蛋白和轉(zhuǎn)錄因子相關(guān)的基因冰寻。

## 選擇前15個(gè)最佳擬合的基因
sg_allgenes <- filter_switchgenes(sce_p1, allgenes = TRUE, topnum = 15)
## 選擇前20個(gè)最佳擬合的蛋白或TF相關(guān)的基因
sg_gtypes <- filter_switchgenes(sce_p1, allgenes = FALSE, topnum = 20,
                                genelists = gs_genelists, genetype = c("Surface proteins", "TFs"))
##合并去重
sg_vis <- rbind(sg_gtypes, sg_allgenes[setdiff(rownames(sg_allgenes), rownames(sg_gtypes)),])

沿著偽時(shí)間線繪制所選基因。開(kāi)啟的基因繪制在線條上方皿渗，而關(guān)閉的基因繪制在線條下方斩芭。

plot_timeline_ggplot(sg_vis, timedata = sce_p1$Pseudotime, txtsize = 3)

image.png

可以使用monocle的降維聚類來(lái)可視化基因表達(dá)和邏輯回歸擬合圖。

plot_gene_exp(sce_p1, gene = "VIM", reduction = "monocleRD", downsample = F)

image.png

GeneSwitches也可以用于基因集通路進(jìn)行排序羹奉。包括MSigDB標(biāo)志性（Liberzon等秒旋，2015年）提供的通路，KEGG的C2和C5基因集诀拭。首先對(duì)超幾何檢驗(yàn)應(yīng)用于提取在軌跡變化中顯著過(guò)表示的路徑迁筛。然后通過(guò)該路徑中基因的中位數(shù)開(kāi)關(guān)時(shí)間來(lái)確定該通路的開(kāi)關(guān)時(shí)間，繪制沿著pseudotime時(shí)間軸排列的山脊圖耕挨。

## 過(guò)濾r^2小于0.1的基因
sg_pw <- filter_switchgenes(sce_p1, allgenes = TRUE, r2cutoff = 0.1)
## 使用超幾何分布细卧，確定時(shí)間點(diǎn)
switch_pw <- find_switch_pathway(rowData(sce_p1), sig_FDR = 0.05,
                                 pathways = msigdb_h_c2_c5, sg_pw)
## 刪除冗余的路徑
switch_pw_reduce <- reduce_pathways(switch_pw, pathways = msigdb_h_c2_c5, 
                                    redundant_pw_rate = 0.8)
##
plot_pathway_density(switch_pw_reduce[1:10,], sg_pw, orderbytime = TRUE)

image.png

還可以根據(jù)自己的研究選擇特定的term進(jìn)行展示

sg_vis <- filter_switchgenes(sce_p1, topnum = 50, pathway_name = c("HALLMARK_MYOGENESIS",
                                                                "GO_CARDIAC_MUSCLE_TISSUE_DEVELOPMENT"))
plot_timeline_ggplot(sg_vis, timedata=sce_p1$Pseudotime, txtsize=3)

image.png

具有兩個(gè)分支的軌跡情況

和之前一樣的操作，確定每個(gè)基因的時(shí)間點(diǎn)筒占，過(guò)濾不符合要求的基因等

sce_p2 <- binarize_exp(sce_p2)
sce_p2 <- find_switch_logistic_fastglm(sce_p2, downsample = TRUE, show_warnings = FALSE)
sg_p1 <- filter_switchgenes(sce_p1, allgenes = TRUE, r2cutoff = 0.03)
sg_p2 <- filter_switchgenes(sce_p2, allgenes = TRUE, r2cutoff = 0.03)

然后贪庙，繪制兩條路徑之間的常見(jiàn)轉(zhuǎn)換基因以比較它們的順序。

sg_com <- common_genes(sg_p1, sg_p2, r2cutoff = 0.4,
                       path1name = "Definitive CM", path2name = "non-contractile")
common_genes_plot(sg_com, timedata = sce_p1$Pseudotime)

image.png

sg_disgs_scale <- distinct_genes(sg_p1, sg_p2, r2cutoff = 0.52, 
                                 path1name = "Definitive CM", path2name = "non-contractile",
                                 path1time = sce_p1$Pseudotime, path2time = sce_p2$Pseudotime, 
                                 scale_timeline = T, bin = 100)
# timedata need to be 1 to (number of bins)
plot_timeline_ggplot(sg_disgs_scale, timedata = 1:100, color_by = "Paths", 
                     iffulltml = FALSE, txtsize = 3)

image.png

這兩張針對(duì)不同轉(zhuǎn)換基因的圖僅顯示了一系列發(fā)生轉(zhuǎn)換事件的偽時(shí)間線翰苫。這個(gè)范圍實(shí)際上處于軌跡的末端止邮，而共同基因大多處于早期（共同基因圖）。
同樣奏窑，我們可以檢查兩條路徑的基因表達(dá)圖导披。

gn <- "DCN"
p <- plot_gene_exp(sce_p1, gene = gn, reduction = "monocleRD", 
                   downsample = FALSE, fitting = TRUE)
p <- plot_gene_exp(sce_p2, gene = gn, reduction = "monocleRD", 
                   downsample = FALSE, fitting = TRUE)

image.png

image.png