博客

Boxplot绘制步骤

1. 导入数据–> 2. 数据预处理 –> 3. 统计检验+pvalue校正 –> 4. pvalue加到boxplot中

1. 导入数据

   library(openxlsx)
   setwd("/Users/mac/Desktop/Imperial\ College\ of\ London/1--data")

2. 数据预处理

##去除data中Group.name.3中的空值
data = data[-which(is.na(data$Group.name.3)),]
​
##将行名中的"-"," "变成”_“如果是"-"," "会影响后面做图 
colnames(data) <- gsub("-","_",colnames(data))
colnames(data) <- gsub(" ","_",colnames(data)

3. 统计检验+pvalue校正

##定义p_value校正的function 选择方法 默认为wilcox.test统计检验 与 FDR校正
get_adj_p <- function(data, .col, .grp = "Sample", comparisons = NULL,
                      method = "wilcox.test", p.adjust.method = "fdr", p.digits = 3L, ...)
 {
  # Compute p-values
  comparison.formula <- paste0(.col, "~", .grp) %>%
    as.formula()
  pvalues <- ggpubr::compare_means(
    formula = comparison.formula, data = data,
    method = method,
    p.adjust.method = p.adjust.method,
    ...
  )
  
  # If a comparison list is provided, extract the comparisons of interest for plotting
  if (!is.null(comparisons)) {
    pvalues <- purrr::map_df(comparisons, ~ pvalues %>% dplyr::filter(group1 == .x[1] & group2 == .x[2]))
  }
  
  # P-value y coordinates
  y.max <- data %>%
    dplyr::pull(.col) %>%
    max(na.rm = TRUE)
  p.value.y.coord <- rep(y.max, nrow(pvalues))
  
  step.increase <- (1:nrow(pvalues)) * (y.max / 10)
  p.value.y.coord <- p.value.y.coord + step.increase
  
  pvalues <- pvalues %>%
    dplyr::mutate(
      y.position = p.value.y.coord,
      p.adj = format.pval(.data$p.adj, digits = p.digits)
    )
  
  pvalues
}



​####定义比较的组 按自己需求
my_comparisons <- list( c("Unhealthy_Control", "Healthy_Control"), c("Unhealthy_Control", "Healthy_Athletes"), c("Healthy_Control","Healthy_Athletes"))
####获得p_adj 统计方法可以指定

p_adj <- get_adj_p(data, .col = "XNLD_Tyr", .grp = "Group.name.3",method = "wilcox.test", p.adjust.method = "fdr")
p_adj
# A tibble: 3 × 9#  .y.      group1            group2               p p.adj p.format p.signif method   y.position
#  <chr>    <chr>             <chr>            <dbl> <chr> <chr>    <chr>    <chr>         <dbl>
#1 XNLD_Tyr Unhealthy_Control Healthy_Athletes 0.203 0.30  0.20     ns       Wilcoxon     84751.
#2 XNLD_Tyr Unhealthy_Control Healthy_Control  0.110 0.30  0.11     ns       Wilcoxon     92455.
#3 XNLD_Tyr Healthy_Athletes  Healthy_Control  0.566 0.57  0.57     ns       Wilcoxon    100160.                        

4. pvalue加到boxplot中

##作图根据Group.name.3做图,p.adj会添加到图上
##如果想用*代替value 可以更改代码中stat_pvalue_manual(p_adj, label = "p.adj")
##-->stat_pvalue_manual(p_adj, label = "p.signif")
plot <- ggplot(data,aes_string(x= "Group.name.3" ,color="Group.name.3",y="XNLD_Tyr"))+ geom_boxplot() + 
    geom_jitter(width = 0.2,shape = 20, size =2.5)+
    scale_color_manual(values=c("lightgreen","lightblue","red"))+
    stat_boxplot(geom = "errorbar",size =0.5, width= 0.5)+ theme_classic()+ stat_pvalue_manual(p_adj, label = "p.adj")

左图p.adj 右图p.signif

PS.批量做图

for (i in col.name[8:length(col.name)]){
##nn为了将图保存在不同的变量下
  nn <- paste0('p_', i)
  p_adj <- get_adj_p(data, .col = i, .grp = "Group.name.3")
##aes_string才可以接受"XNLD_Tyr"，aes()直接就是XNLD_Tyr不需要加引号
##但我们是循环所以i是带引号的
  plot <- ggplot(data,aes_string(x= "Group.name.3" ,color="Group.name.3",y=i))+ geom_boxplot() + geom_jitter(width = 0.2,shape = 20, size =2.5)+
    scale_color_manual(values=c("lightgreen","lightblue" ,"red"))+
    stat_boxplot(geom = "errorbar",size =0.5, width= 0.5)+ theme_classic()+ stat_pvalue_manual(p_adj, label = "p.signif")
  #+ stat_compare_means(comparisons = my_comparisons,method = "wilcox.test",p.adjust.method ="BH") 
  ##as.ggplot是为了后面可以plot_grid
  plot <- as.ggplot(plot)
  assign(nn, plot)
  result <- paste0(result,",`",nn,"`")
}
result
##作图
plot_grid(`p_D_Ala`,`p_D_Ser`,`p_D_His`,`p_D_Cit`,`p_D_Pro`,`p_D2_ABA`,`p_D_Val`
          ,nrow =3, axis = "l", align = "v",labels="AUTO")

如果不足，欢迎批评指正