0.起因

最近在群里看到一个图，棒棒糖图的点用饼图代替了。

然后又看到一个文章里，用饼图代替了富集分析的点图，妙啊~

这个地球上有什么漂亮的图能躲过我的小胖手呢？画他！

1.搜索

试了几个不同的关键词，最后这样搜到了。

答案指向一个包，名叫scatterpie。一看作者，啊，这不是Y叔嘛~

2. 跑示例代码

2.1 示例数据

set.seed(123) long <- rnorm(50, sd=100) lat <- rnorm(50, sd=50) d <- data.frame(long=long, lat=lat) d <- with(d, d[abs(long) < 150 & abs(lat) < 70,]) n <- nrow(d) d$region <- factor(1:n) d$A <- abs(rnorm(n, sd=1)) d$B <- abs(rnorm(n, sd=2)) d$C <- abs(rnorm(n, sd=3)) d$D <- abs(rnorm(n, sd=4)) d[1, 4:7] <- d[1, 4:7] * 3 head(d) ## long lat region A B C D ## 1 -56.047565 12.665926 1 2.13121969 8.663359 3.928711 8.676792 ## 2 -23.017749 -1.427338 2 0.25688371 1.403569 1.375096 4.945092 ## 4 7.050839 68.430114 3 0.24669188 0.524395 3.189978 5.138863 ## 5 12.928774 -11.288549 4 0.34754260 3.144288 3.789556 2.295894 ## 8 -126.506123 29.230687 5 0.95161857 3.029335 1.048951 2.471943 ## 9 -68.685285 6.192712 6 0.04502772 3.203072 2.596539 4.439393

2.2 示例图片

library(ggplot2) library(scatterpie) ggplot() + geom_scatterpie(aes(x=long, y=lat), data=d, cols=LETTERS[1:4]) + coord_equal()

3.仿造示例数据并画图

3.1找个数据做差异分析

library(tinyarray) library(dplyr) library(org.Hs.eg.db) gse = "GSE42872" geo = geo_download(gse) Group = rep(c("control","treat"),each = 3) Group = factor(Group) find_anno(geo$gpl) ## [1] "`library(hugene10sttranscriptcluster.db);ids <- toTable(hugene10sttranscriptclusterSYMBOL)` and `ids <- AnnoProbe::idmap('GPL6244')` are both avaliable" ids <- AnnoProbe::idmap(geo$gpl,destdir = tempdir()) deg = get_deg(geo$exp,Group,ids,logFC_cutoff = 2) head(deg) ## logFC AveExpr t P.Value adj.P.Val B probe_id ## 1 5.780170 7.370282 82.94833 3.495205e-12 1.163798e-07 16.32898 8133876 ## 2 -4.212683 9.106625 -68.40113 1.437468e-11 2.393169e-07 15.71739 7965335 ## 3 5.633027 8.763220 57.61985 5.053466e-11 4.431880e-07 15.04752 7972259 ## 4 -3.801663 9.726468 -57.21112 5.324059e-11 4.431880e-07 15.01709 7972217 ## 5 3.263063 10.171635 50.51733 1.324638e-10 8.821294e-07 14.45166 8129573 ## 6 -3.843247 9.667077 -45.87910 2.681063e-10 1.487856e-06 13.97123 8015806 ## symbol change ENTREZID ## 1 CD36 up 948 ## 2 DUSP6 down 1848 ## 3 DCT up 1638 ## 4 SPRY2 down 10253 ## 5 MOXD1 up 26002 ## 6 ETV4 down 2118 table(deg$change) ## ## down stable up ## 48 18092 82

3.2 用差异基因做富集分析

gene_up = deg$symbol[deg$change=="up"] gene_down = deg$symbol[deg$change=="down"] gene_diff = c(gene_up,gene_down) length(gene_diff) ## [1] 130 head(gene_diff) ## [1] "CD36" "DCT" "MOXD1" "NUPR1" "ST6GALNAC2" ## [6] "TXNIP" g = quick_enrich(gene_diff,kkgo_file = "kk_gop.Rdata") dat = g$go@result %>% arrange(pvalue) %>% filter(onTOLOGY=="BP") %>% head(10) %>% mutate(Description = factor(Description,levels = Description)) p1 = ggplot(dat)+ geom_point(aes(x = -log10(pvalue),y = Description))+ theme_bw() p1

3.3 计算每条term的上下调基因数量

library(stringr) p = str_split(dat$geneID,"/") names(p) = dat$ID get_num = function(x){ c(sum(x %in% gene_up),sum(x %in% gene_down)) } re = sapply(p, get_num) %>% t() %>% as.data.frame() rownames(re) = names(p) colnames(re) = c("up","down")

3.4 仿制示例数据并画图

试了一下，如果纵坐标是分类变量是会报错的。。。

所以曲线救国，先生成一组连续型数据，然后改刻度值。

而需要保持饼图的形状是圆的，又需要固定坐标系,横纵坐标刻度是一样的。

这时候难度就来了，需要生成的纵坐标得是和横坐标范围一致的等差数列，用seq函数。

dp = cbind(lp = -log10(dat$pvalue), Description = dat$Description, re) dp$lab = seq(min(dp$lp), max(dp$lp), length.out = nrow(dp)) head(dp) ## lp Description up ## GO:0071900 7.264160 regulation of protein serine/threonine kinase activity 5 ## GO:0050673 6.929666 epithelial cell proliferation 8 ## GO:0043405 6.351294 regulation of MAP kinase activity 5 ## GO:0009314 5.748852 response to radiation 6 ## GO:0010165 5.691160 response to X-ray 0 ## GO:0071214 5.635287 cellular response to abiotic stimulus 6 ## down lab ## GO:0071900 10 5.525307 ## GO:0050673 8 5.718513 ## GO:0043405 5 5.911719 ## GO:0009314 8 6.104925 ## GO:0010165 5 6.298131 ## GO:0071214 6 6.491337 p2 = ggplot() + geom_scatterpie(aes(x=lp, y=lab), data=dp, cols=colnames(re), pie_scale = 1.5) + theme_bw()+ scale_fill_manual(values = c("#f87669", "#2874C5"))+ scale_y_continuous(labels = dp$Description,breaks = dp$lab)+ labs(x = "-log10(pvalue)",y = "Description")+ coord_equal() p2

3.5 要棒棒糖还不简单嘛

p2 + geom_segment(data = dp, aes(x=as.integer(min(dp$lp)), xend=lp, y=lab, yend=lab))