ggplot的boxplot添加顯著性 | Add P-values and Significance Levels to ggplots | 方差分析

參考：Add P-values and Significance Levels to ggplots

 

ggpubr的包比較侷限，能用的test也比較侷限，可是作起來快速簡單。

當狀況特殊時ggpubr就不能用了，能夠本身作了顯著性test以後再顯示在圖上。

# show lable in facet grid plot
dat_text <- data.frame()
for (i in names(paired_list)) {
    # Compute t-test
    res <- t.test(value ~ group, data = paired_list[[i]], paired = TRUE)
    dat_text <- rbind(dat_text, data.frame(variable=i, pvalue=res$p.value))
}

dat_text$label <- paste("P", round(dat_text$pvalue, 3), sep="=")

dat_text[dat_text$pvalue<0.05 & dat_text$pvalue>0.01,]$label <- paste("*", 
                                  dat_text[dat_text$pvalue<0.05 & dat_text$pvalue>0.01,]$label, sep=" ")

dat_text[dat_text$pvalue<0.01,]$label <- paste("**", "P<0.01", sep=" ")

    library(ggplot2) 
    options(repr.plot.width=8, repr.plot.height=12) # 8x8
    g2 <- ggplot(data=genes_expr_melt, aes(x=pseudotime, y=value, fill=group, color=group)) +
        geom_point(size=0.01, alpha=0.5, aes(color=group, fill=group)) +
        labs(x = "Pseudotime", y = "Relative expression", title = "Neuronal lineage") +
        geom_smooth(method = 'loess',se=F,size=0.15,span = 0.7) + # ,alpha=0.05, weight=0.1,
        facet_wrap( ~ variable, ncol=3, labeller = label_context, scales = "free_y") + # 
        geom_text(size    = 5, data    = dat_text, mapping = aes(x = Inf, y = Inf, label = label),
            hjust   = 1.05,vjust   = 1.5, color=ifelse(dat_text$pvalue < 0.05,'red','black')) +
        # themes
        theme(strip.background = element_blank(),
             panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
             panel.spacing=unit(.4, "lines"),panel.border = element_rect(color = "black", fill = NA, size = 0.5))+
        theme(axis.text.x  = element_text(face="plain", angle=0, size = 10, color = "black", vjust=0.5),
            axis.text.y  = element_text(face="plain", size = 10, color = "black"),
            axis.title =element_text(size = 15)) +
        theme(strip.background = element_rect(fill = "gray90", color = NA))+
          # theme(legend.position = "none") + # must remove legend
          theme(strip.placement = "outside", strip.text.x = element_text(face="plain", size = 13),
                  strip.text.y = element_text(face="plain", size = 11)) +
        theme(strip.text.x = element_text(margin = margin(1,0,1,0, "mm"))) +
        scale_color_manual(values=c("deepskyblue","red","gray50")) +
        scale_fill_manual(values=c("deepskyblue","red","gray50"))
    plot(g2)
# }

 

多組比較，挑選感興趣的顯示顯著性。

data("ToothGrowth")
head(ToothGrowth)　　

library(ggpubr)
my_comparisons <- list( c("0.5", "1"), c("1", "2"), c("0.5", "2") )
options(repr.plot.width=4, repr.plot.height=4)
ggplot(ToothGrowth, aes(x=as.character(dose), y=len, fill=dose)) +
  geom_boxplot(outlier.size=NA, size=0.01, outlier.shape = NA) +
  geom_jitter(width = 0.3, size=0.01) +# , aes(color=supp) +
  stat_compare_means(comparisons = my_comparisons)+ # Add pairwise comparisons p-value
  stat_compare_means(label.y = 50, label.x = 1.5)     # Add global p-value

還能夠設定一個ref group來顯示顯著性差別，只須要改一下設定。

  stat_compare_means(method = "anova", label.y = 1.3, label.x = 3)+ # Add pairwise comparisons p-value
  # # Add global p-value
  stat_compare_means(label = "p.signif", method = "t.test", ref.group = "hNP-D20", label.y = 1.1) + 

　　

 

生物學的強烈推薦看看Y叔的公衆號裏的統計相關的文章，很是的基礎和實用。

統計

• Five things biologists should know about statistics
• 什麼是T檢驗
• 富集基因之註釋缺失
• 落入窠臼
• 你昨天才作的分析，多是幾年前的結果！
• 掐架的額外收穫
• boxplot
• 如何告別單身
• 主成分分析
• 一文解決RT-PCR的統計分析

 

代碼例子：

options(repr.plot.width=7, repr.plot.height=6)
# facet boxplot
bp <- ggplot(expr_data2, aes(x=group, y=expression, fill=NA)) + 
  geom_boxplot(outlier.size=NA, size=0.01, outlier.shape = NA) + 
  geom_jitter(width = 0.3, size=0.01, aes(color=cluster)) +
  # + geom_boxplot( + 
  facet_grid( cluster ~ gene, switch="y") + # , scales = "free"
  theme_bw() + 
  stat_compare_means(aes(group = group, label = ..p.signif..), label.x = 1.3,label.y = 1.3, 
                     method = "wilcox.test", hide.ns = T) + # label = "p.format",
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(x = "", y = "", title = "") +
  theme(panel.spacing=unit(.3, "lines"),panel.border = element_rect(color = "black", fill = NA, size = 0.2)) +
  theme(axis.ticks.x = element_blank(), axis.ticks = element_line(size = 0.1), 
        axis.text.x  = element_text(face="plain", angle=90, size = 8, color = "black", vjust=0.5),
        axis.text.y  = element_text(face="plain", size = 4, color = "black"),
        axis.title =element_text(size = 12)) +
  theme(strip.background = element_rect(fill = "gray97", color = NA))+
  theme(legend.position = "none") +
  theme(strip.placement = "outside", strip.text.x = element_text(face="italic", size = 11),
          strip.text.y = element_text(face="plain", size = 11)) +
  scale_y_continuous(position="right", limits = c(-0.5,1.5)) +
  scale_fill_manual(values=brewer.pal(8,"Set2")[c(2,3,7,1,5,6)]) +
  scale_color_manual(values=brewer.pal(8,"Set2")[c(2,3,7,1,5,6)])
bp