未测试 批量绘制配对样本差异基因的箱线图

参考生信星球

今天讲批量绘制GEO配对样本差异基因的箱线图。

1. 获取数据

这个例子是GSE5109，由6个配对样本组成，输入数据是表达矩阵、分组信息、配对信息以及limma差异分析得到的结果表格。

load("geo_pair_boxplot.Rdata")
expm[1:4,1:4]
#>        GSM101299 GSM101300 GSM101301 GSM101302
#> STK26   5.754235  6.659144  2.233099  6.353956
#> VMO1    9.222983  7.538749  5.291330  6.947587
#> CNKSR2  7.233939  7.306732  4.274562  7.785351
#> IL17RE  7.604526  7.677775  4.684485  7.409410
group_list
#> [1] post post pre  post pre  pre 
#> Levels: pre < post
pairinfo
#> [1] 1 2 1 3 2 3
#> Levels: 1 2 3
head(deg)
#>       logFC  AveExpr         t      P.Value adj.P.Val         B
#> 1  2.662621 4.373020  12.93505 0.0001471005 0.9278859 -4.320673
#> 2  2.467153 6.158566  10.93078 0.0002945367 0.9278859 -4.324235
#> 3  2.237728 5.859695  10.80083 0.0003093558 0.9278859 -4.324533
#> 4  2.039293 6.121906  10.60373 0.0003336256 0.9278859 -4.325005
#> 5  2.055491 5.700167  10.39701 0.0003616455 0.9278859 -4.325528
#> 6 -1.955322 6.733026 -10.06996 0.0004121430 0.9278859 -4.326419
#>       probe_id symbol change        v ENTREZID
#> 1  224407_s_at  STK26     up 3.832386    51765
#> 2  235751_s_at   VMO1     up 3.530861   284013
#> 3   1554607_at CNKSR2     up 3.509542    22866
#> 4    229401_at IL17RE     up 3.476741   132014
#> 5   1558323_at TMEM72     up 3.441717   643236
#> 6 1556411_s_at     KY   down 3.384952   339855

2.作图

我们以logFC绝对值最大的几个基因作箱线图，比较每个人在手术前后的基因表达量。做几张图就改n等于几，这里以10张图为例。

通过搜索和探索解决了以下几个问题：

1.实现顺序调整：上调和下调基因分别在一起，不管几个，都是先画下调，再画上调。

2.ggplot2批量绘图的传参问题，aes_string代替aes，就可以使用列名取子集得到的字符串了。

3.patchwork多子图拼图，可以将子图组成一个列表，作为wrap_plots的参数，这一点完美代替了gridarrange。

n = 10
library(ggplot2)
x = deg$symbol[order(head(deg$logFC,n))]
dat <- data.frame(pairinfo=pairinfo,group=group_list,t(expm[x,]))
pl = list()
for(i in 1:n){
  pl[[i]] = ggplot(dat, aes_string("group",colnames(dat)[i+2],fill="group")) +
    geom_boxplot() +
    geom_point(size=2, alpha=0.5) +
    geom_line(aes(group=pairinfo), color="black", linetype="11") +
    xlab("") +
    ylab(paste("Expression of",colnames(dat)[i+2]))+
    theme_classic()+
    theme(legend.position = "none")+
    scale_fill_manual(values = c("#00AFBB", "#E7B800"))
}
#拼图
library(patchwork)
pb_top10= wrap_plots(pl,nrow = 2)+plot_annotation(tag_levels = 'a')
pb_top10