Блог

У меня есть цикл, который создает список ggplots для любого заданного количества анализируемых веществ (в этом примере 3). Выполняется ANOVA, генерируется и отображается p-значение для каждого ggplot.

Тем не менее, я также хочу пометить эти графики скобками с p-значением Tukey HSD. Если возможно, я хочу визуализировать только скобки, которые находятся ниже скорректированного p < 0,05.

Я написал два скрипта, которые выполняют эти действия отдельно; выводом для моего цикла является список списков (элементов ggplot), а выводом из моего скрипта Tukey HSD является фрейм данных.

Скобки, на которые я ссылаюсь, будут выглядеть примерно так, как в этом примере ggplot:

Ниже приведен мой код. Во-первых, данные:

 set.seed(10)
Label<-as.data.frame(c("Baseline","Baseline","Baseline","Baseline","A","B","B","C","D"))
Label<-do.call("rbind", replicate(10, Label, simplify = FALSE))
colnames(Label)<-"Label"
Analyte1<- rnorm(90, mean = 1, sd = 1)
Analyte2<- rnorm(90, mean = 1, sd = 1)
Analyte3<- rnorm(90, mean = .2, sd = 1)
df<-cbind(Label,Analyte1,Analyte2,Analyte3)
 

Ниже приведен мой цикл ANOVA:

 library(dplyr)
library(ggplot2)
library(ggpubr)
library(rstatix)
library(rlist)

# Select numeric columns to obtain df length.
sample_list <-
  colnames(select_if(df, is.numeric))
sample_list

# Create a list where the plots will be saved.
ANOVA.plots <- list()

# The for loop.
for (i in 1:length(sample_list)) {
  ANOVA.ggplot <-
    ggplot(
      df,
      aes_string(
        x = "Label",
        y = sample_list[i],
        fill = "Label",
        title = sample_list[i],
        outlier.shape = NA
      )
    )  
    border(color = "black", size = 2.5)  
    geom_jitter(
      aes(fill = Label),
      shape = 21,
      size = 4,
      color = "black",
      stroke = 1,
      position = position_jitter()
    )  
    geom_boxplot(alpha = 0.7,
                 linetype = 1,
                 size = 1.1)  
    theme(legend.position = "none")  
    scale_fill_brewer(palette = "Set2")  
    stat_boxplot(geom = "errorbar",
                 size = 1.1,
                 width = 0.4)  
    theme (axis.title.y = element_blank())  
    theme (axis.title.x = element_blank())  
    stat_compare_means(
      inherit.aes = TRUE,
      data = df,
      method =  "anova",
      paired = FALSE ,
      method.args = list(var.equal = TRUE),
      fontface = "bold.italic",
      size = 5,
      vjust = 0,
      hjust = 0
    )
  ANOVA.plots[[i]] <-  ANOVA.ggplot
}

# Print each ggplot element from the list of lists.
ANOVA.plots
 

В настоящее время этот фрагмент генерирует графики, похожие на этот из Analyte 2.

Наконец, этот фрагмент предназначен для получения фрейма данных Tukey HSD.

 TUKEY.length<-length(sample_list)

TUKEY.list <-
  vector(mode = "list", length = TUKEY.length) # Empty list for looping, where Tukey HSD results are stored.

for (i in 1:TUKEY.length   1) {
  TUKEY.list[[i]] <-
    aov(df[[i]] ~ df[[1]]) %>% tukey_hsd() # Obtain Tukey HSD p-values from comparing groups.
}

TUKEY.list <-
  TUKEY.list[lengths(TUKEY.list) > 0L] # Remove empty lists, artifacts from piping.

p.value.threshold <- 0.05

TUKEY.df<-list.rbind(TUKEY.list)
Analytes <-
  colnames(df[, sapply(df, class) %in% c('integer', 'numeric')]) 
# Stores the analyte names for Tukey. Only pulls numeric colnames.

# The following was done to properly bind the analyte IDs to the Tukey data frame.
Analytes.df<-as.data.frame(Analytes)
Analytes.df$Value<-c(1:nrow(Analytes.df))
Analytes.df<-do.call("rbind", replicate(10, Analytes.df, simplify = FALSE))

Analytes.df <-
  Analytes.df[order(Analytes.df$Value), ] 
# Orders the dataframe so that the rownames can be assigned to the TUKEY HSD.
Analytes.df<-as.data.frame(Analytes.df[,-2])

toDrop <- c("^term", "^null") # Columns to drop, left over from Tukey loop.

TUKEY.df<-TUKEY.df[,!grepl(paste(toDrop, collapse = "|"),names(TUKEY.df))]
TUKEY.df<-cbind(Analytes.df, TUKEY.df)
TUKEY.df<-filter(TUKEY.df, p.adj <= p.value.threshold)
 

Желаемый результат по-прежнему будет отображать p-значение ANOVA и содержать скобки для тех групп, у которых скорректированное значение p ниже 0,05 после Tukey HSD.

В этом примере только анализируемый материал 2, группы A и B, имеет скорректированное значение p ниже 0,05 (т. Е. p = 0,000754), поэтому над этими двумя графиками должна появиться скобка. Однако код для отображения этих скобок Tukey-HSD должен быть способен работать с большим количеством переменных (т.Е. Более 20 анализируемых веществ) и фиксировать все попарные сравнения, которые являются значимыми для каждого графика, а не только для того, который я привел.

Вы могли бы использовать ggsignif для пользовательских аннотаций. В примере используется p < 0,2 в качестве ограничения для отображения нескольких полос ошибок:

 library(ggplot2)
library(ggsignif)
library(cowplot)
library(data.table)

set.seed(10)
df <- data.table(Label=rep(c("Baseline","Baseline","Baseline","Baseline","A","B","B","C","D"), 10),
                 `Analyte 1` = rnorm(90, mean = 1, sd = 1), 
                 `Analyte 2` = rnorm(90, mean = 1, sd = 1), 
                 `Analyte 3` = rnorm(90, mean = .2, sd = 1))
df[, Label := factor(Label, unique(Label))]
df <- melt(df, id.vars = "Label", variable.name = "Analyte")
dfl <- split(df, df$Analyte, drop=TRUE)

doPlots <- function(x, signif.cutoff=.05){
    set.seed(123)
    p1 <- ggplot(dfl[[x]], aes(x=Label, y=value, fill=Label))  
        geom_boxplot(alpha = 0.7, linetype = 1, size = 1.1)  
        theme(legend.position = "none")  
        scale_fill_brewer(palette = "Set2")  
        stat_boxplot(geom = "errorbar", size = 1.1, width = 0.4)  
        geom_jitter(aes(fill = Label), shape = 21, size = 4, stroke = 1)  
        ggtitle(x)
    a <- stats::TukeyHSD(stats::aov(value ~ Label, data = dfl[[x]]))[[1]]
    a <- stats::setNames(
        data.frame(
            do.call(rbind, strsplit(rownames(a), "-")),
            a[, "p adj"]
        ), c("Var1", "Var2", "p") )
    a <- subset(a[stats::complete.cases(a),], p < signif.cutoff)
    if (nrow(a) == 0) return(p1) else {
        a$p <- formatC(signif(a$p, digits = 3),
                          digits = 3, format = "g", flag = "#")
        keep.tests <- unname(t(apply(a[, -3], 1, sort)))
        keep.tests <- unname(split(keep.tests, seq(dim(keep.tests)[1])))
    ts <- list(a = a, keep.tests = keep.tests)
    p1   ggsignif::geom_signif(
        comparisons=ts$keep.tests,
        test="TukeyHSD",
        annotations=ts$a$p, 
        step_increase=0.1)
    }
}

plot_grid(plotlist=lapply(names(dfl), doPlots, signif.cutoff=.2), ncol=3)
 

^{Создано 2021-01-29 пакетом reprex (v1.0.0)}