多基因分组条形图

生信喵

背景

画出PMID：29720944文章同款的多基因分组条形图，以及多个单基因小图拼接大图。

使用场景

用boxplot展示基因表达量分布，同时计算p value。

输入数据

require(ggplot2)
require(tidyr)
require(cowplot)

1. 基因表达量列表

easy_input_expr.csv，多个sample中所有基因的表达量。

包含stage1, stage2, stage3，各10个样本，stage以样本名后面 _数字标识。

getwd()
d <- read.csv("easy_input_expr.csv", row.names=1)
dim(d) #12403 30
head(d)

d前几行如下：

2. 差异表达的基因名列表

单独一列基因名；

或者第一列是基因名，后面是其他信息。

up<-read.table("easy_input_up.txt",header = T, as.is = T, row.names = 1)
down<-read.table("easy_input_down.txt",header = T, as.is = T, row.names = 1)# 各13个基因

开始画图

1. 多组画在同一个绘图区域

分别绘制上下调的基因，然后用cowplot来拼。

require(tidyr)
require(ggplot2)

box <- lapply(list(up=row.names(up), down=row.names(down)), function(ii) {
    dd <- d[ii,]
    dd$gene = rownames(dd)
  
    ## 宽变长，变成tidy data frame，好用于ggplot2画图
    d2 <- gather(dd, sample, expr, 1:30)
  
    ## 通过sample名，衍生出stagge变量
    d2$stage <- paste('stage', sub(".*_(\\d)$", "\\1", d2$sample))

    ## 用one way anova计算 p value
    pvalues <- sapply(d2$gene, function(x) {
        res <- aov(expr ~ stage, data = subset(d2, gene == x))
        summary(res)[[1]]$'Pr(>F)'[1] #
    })
    pv <- data.frame(gene = d2$gene, pvalue = pvalues)

    ## 画boxplot，并标识p-value
    ggplot(d2, aes(gene, expr, fill=stage)) + 
      geom_boxplot() + 
      geom_text(aes(gene, y=max(d2$expr) * 1.1, 
                    label=paste("p = ",round(pvalue, 2))),
                data=pv, #
                inherit.aes=F) + #
      xlab(NULL)+ylab("Relative expression (log2)")
})

require(cowplot)
plot_grid(plotlist=box, labels=c("A", "B"), ncol=1)
ggsave(file="box1.pdf")

如果想用“*”或“**”来标记是否<0.05或<0.001，这样做：

把pvalue变成sigcode，然后再下面画图代码中 label=paste("p = ",round(pvalue, 2))标记p值的地方，改成 label = sigcode即可。

box <- lapply(list(up=row.names(up), down=row.names(down)), function(ii) {
    dd <- d[ii,]
    dd$gene = rownames(dd)
    d2 <- gather(dd, sample, expr, 1:30)
    d2$stage <- paste('stage', sub(".*_(\\d)$", "\\1", d2$sample))
    pvalues <- sapply(d2$gene, function(x) {
        res <- aov(expr ~ stage, data = subset(d2, gene == x))
        summary(res)[[1]]$'Pr(>F)'[1] #
    })
    pv <- data.frame(gene = d2$gene, pvalue = pvalues)
  
    ## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    pv$sigcode <- cut(pv$pvalue, c(0, 0.001, 0.01, 0.05, 0.1, 1), 
                    labels=c('***', '**', '*', '.', ' '))

    ggplot(d2, aes(gene, expr, fill=stage)) + 
      geom_boxplot() + 
      geom_text(aes(gene, y=max(d2$expr) * 1.1, 
                    label=pv$sigcode),
                data=pv, inherit.aes=F) + 
      xlab(NULL)+ylab("Relative expression (log2)")
})

require(cowplot)
plot_grid(plotlist=box, labels=c("A", "B"), ncol=1)
ggsave(file="box2.pdf")

2. 多个小图拼成一个大图

box2 <- lapply(c(row.names(up), row.names(down)), function(i) {
    dd <- d[i,]
    dd$gene = rownames(dd)
    d2 <- gather(dd, sample, expr, 1:30)
  
    ## 通过sample名，衍生出stagge变量
    d2$stage <- paste('stage', sub(".*_(\\d)$", "\\1", d2$sample))

    ## 用one way anova计算 p value
    res <- aov(expr ~ stage, data = d2)
    pv1 <- summary(res)[[1]]$'Pr(>F)'[1]
    pv1.lab <- paste("one-way ANOVA p =", round(pv1,3))

    ## t test for comparing stage 3 vs stage 1
    pv2 <- t.test(expr~stage, data=subset(d2, stage!='stage 2'))$p.value
    pv2.lab <- paste("stage III vs stage I p =", round(pv2, 3))

    lab <- paste(dd$gene, pv1.lab, pv2.lab, sep="\n")

    label <- ggdraw() + draw_label(lab)
    plot_grid(label, 
              ggplot(d2, aes(stage, expr, fill=stage)) +
                geom_boxplot() +
                xlab("") +
                ylab(paste("Relative ",dd$gene," expression (log2)")) +
                theme(axis.title.y = element_text(size = 10)), #调整y轴标题字号
              ncol=1, rel_heights=c(.3, 1))
})

plot_grid(plotlist=box2, ncol=6)
ggsave(file="box3.pdf")

生信喵

文中所用数据可以关注公众号“生信喵实验柴”

发送关键词“20250324”获取