Run ANOVA on every column of a dataframe (cols representing mRNA/proteins)

Question

In a data frame of features (eg proteins or mRNA transcripts), I want to calculate a one-way ANOVA for every column (representing expression values). Each row represents a sample (id) that belongs to a group (5 groups in example data).

Searching the web and Stackoverflow only gave me non-functional and partial answers, the expected result would be a table with the ANOVA p-value for every feature (eg a df with first col = feature, second col = ANOVA p-value). Below the example data I post a non-working solution.

For pairwise differential expression analysis between the groups, I am using LIMMA which delivers multiple-comparison corrected p-values. If I understand correctly, ANOVA does not care for mul.comp.cor but only shows for which features there is a difference between any of the 5 groups (and a post-test, eg LIMMA, shows me in which of the groups there is the difference). It would be great if you could confirm that this way of assessing differential expression between my 5 groups is correct.

The solution was taken from here: One-way ANOVA for loop: how do I iterate through multiple columns of a dataframe?

example data

# gene df and list 
genes <- paste("gene",1:1000,sep="")
x <- list(
  A = sample(genes,300), 
  B = sample(genes,525), 
  C = sample(genes,440),
  D = sample(genes,350)
)

# expression dataframe
crete_exp_df <- function(gene_nr, sample_nr){
  df <- replicate(sample_nr, rnorm(gene_nr))
  rownames(df) <- paste("Gene", c(1:nrow(df)))
  colnames(df) <- paste("Sample", c(1:ncol(df)))
  return(df)
}

df1 <- crete_exp_df(50, 20)
df1 <- as.data.frame(df1)
df1$fid <- rownames(df1)

# creator for ANOVA
df4ANOVA <- df1 %>% pivot_longer(-fid) %>%
  pivot_wider(names_from="fid", values_from="value") %>%
  rename(id=name)

df4ANOVA$group <- c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5)

solution (not working)

library(tidyverse)
library(broom)

df4ANOVA$group <- as.factor(df4ANOVA$group )
df4ANOVA$id <- NULL

# ANOVA
baseformula <- " ~ group"

ANOVA_group <- for (i in 2:ncol(exp4anov)) {
 formula <- paste(colnames(exp4anov)[i], baseformula, sep="")

 p <- summary(aov(as.formula(formula), data=exp4anov))[[1]][["Pr(>F)"]][1]

 print(paste(formula, ": p=", p, sep=""))

}

Answer 1

A possibility is the ffmanova package:

library(ffmanova)
res <- ffmanova(as.matrix((df4ANOVA[1:50]))~ df4ANOVA$group)

You will find the classical one-way p-values as the second row of res$pRaw. That is, res$pRaw[2, ].

By adding, for example, nSim = 9999, to the ffmanova call, adjusted alternatives to pRaw are also available (res$pAdjusted and res$pAdjFDR).

Answer 2

You can use lapply to loop over your df4ANOVA data frame, using grep to only identify columns with Gene:

aov_list <- lapply(df4ANOVA[grep("Gene", colnames(df4ANOVA))], 
                   function(x) aov(df4ANOVA$group ~ x))

This results in a list of length 50, each containing the full results of aov for the column.

If you want a final data frame with p values and gene names, you could first tweak the lapply statement to just extract p values, then do.call to put them together:

aov_plist <- lapply(df4ANOVA[grep("Gene", colnames(df4ANOVA))], 
                   function(x) summary(aov(df4ANOVA$group ~ x))[[1]][5][1,])

finaldat <- data.frame(pval = do.call(rbind, aov_plist))
finaldat$gene <- rownames(finaldat)

Output:

#                pval    gene
# Gene 1  0.551729974  Gene 1
# Gene 2  0.725349369  Gene 2
# Gene 3  0.983277774  Gene 3
# Gene 4  0.811234760  Gene 4
# Gene 5  0.765013222  Gene 5
# Gene 6  0.144887277  Gene 6
# ...

A "cleaner" way to do this may be to use dplyr::bind_rows and tibble::rownames_to_column:

data.frame(t(bind_rows(aov_plist))) %>% 
  tibble::rownames_to_column("Gene")

Both get you functionally the same resulting data frame.