两分组的交集差集可视化最常用的方法就是韦恩图,复杂点如果分组超过5组可以换成upset
图展示。虽然也是两分组,但如果有多个分类时韦恩图就不是好的选择了。此时,回到最简单的展示方式,用堆叠条形图确实是一个很好的选择,有种言简意赅的感觉,恰到好处地展示了数据。
比如,单细胞数据有多个细胞亚群,分别使用两种方法分析各亚群的差异基因,然后比较方法在各亚群间的结果情况。这个时候用堆叠条形图来展示就是不错的选择。
str(dge1)
'data.frame': 3837 obs. of 7 variables:
$ p_val : num 1.87e-139 1.60e-138 4.97e-138 4.61e-137 7.74e-125 ...
$ avg_log2FC: num -14.94 -17.67 -19.12 -2.69 1.79 ...
$ pct.1 : num 0.809 0.832 0.858 0.991 1 1 0.994 0.726 0.493 0.991 ...
$ pct.2 : num 0.951 0.966 0.976 1 0.999 1 1 0.926 0.839 0.998 ...
$ p_val_adj : num 1.14e-136 9.72e-136 3.02e-135 2.80e-134 4.71e-122 ...
$ cluster : Factor w/ 9 levels "Naive CD4 T",..: 1 1 1 1 1 1 1 1 1 1 ...
$ gene : chr "RFX5" "RFXANK" "RFXAP" "REL" ...
> str(dge2)
'data.frame': 4068 obs. of 7 variables:
$ p_val : num 1.21e-239 8.03e-239 8.03e-239 8.03e-239 8.03e-239 ...
$ avg_log2FC: num 0.949 0.817 0.817 0.817 0.817 ...
$ pct.1 : num 0.997 0.999 0.999 0.999 0.999 0.999 0.999 0.999 0.999 0.999 ...
$ pct.2 : num 0.481 0.639 0.639 0.639 0.639 0.639 0.639 0.639 0.639 0.639 ...
$ p_val_adj : num 6.69e-237 4.45e-236 4.45e-236 4.45e-236 4.45e-236 ...
$ cluster : Factor w/ 9 levels "Naive CD4 T",..: 1 1 1 1 1 1 1 1 1 1 ...
$ gene : chr "ATF5" "NR2C2" "NRL" "UBTF" ...
比如,上面展示了两种方法得到的9个细胞亚群差异基因,想要可视化展示结果差异,首先要根据差异结果整理需要的绘图数据。准备数据的方式有很多种,这里感受一下R语言一套组合拳的厉害:
shared <- stack(mapply(function(x,y){intersect(x$gene, y$gene)}, split(dge1, ~ cluster), split(dge2, ~ cluster)))
method1 <- stack(mapply(function(x,y){setdiff(x$gene, y$gene)}, split(dge1, ~ cluster), split(dge2, ~ cluster)))
method2 <- stack(mapply(function(x,y){setdiff(x$gene, y$gene)}, split(dge2, ~ cluster), split(dge1, ~ cluster)))
pdata <- rbind(method1, method2, shared)
colnames(pdata) <- c('gene', 'cluster')
pdata$type <- rep(c('method1', 'method2', 'shared'), time=c(nrow(method1), nrow(method2), nrow(shared)))
简短几行代码就搞定了两种方法各亚群差异基因的交集及差集,R语言里面的各种循环函数真的很香,简单一套代码组合拳下来数据就整理好了,谁看了不迷糊。
提醒一下这里使用的R-4.2
,如果使用低版本的R,split
函数调用的时候可能需要稍微改一改。下面就是借助ggplot2
画个堆叠条形图的事了:
library(ggplot2)
ggplot(pdata, aes(cluster, fill=type)) +
geom_bar(stat='count', width=0.6) +
theme_test() + xlab('') +
theme(text=element_text(size=15), axis.text.x=element_text(angle=30, hjust=1, vjust=1))