Question: Tukeys'S Hsd After Anova Using R
2
deschisandy60 wrote:
``````The dataset has 1000 genes and contains 24 samples with two mouse strains tested (129 and B6) and six brain regions. There are two replicates for each region.
``````

The ANOVA was performed as follows:

``````      sdata<-read.table("http://www.chibi.ubc.ca/wp-content/uploads/2013/02/
strain <- gl(2,12,24, label=c("129","bl6"))
region <- gl(6,2,24, label=c("ag", "cb", "cx", "ec", "hp", "mb"))
# define ANOVA function ##  ANOVA
aof <- function(x) {
m<-data.frame(strain,region, x);
anova(aov(x ~ strain + region + strain*region, m))
}
# apply analysis to the data and get the pvalues.
anovaresults <- apply(sdata, 1, aof)

pvalues<-data.frame( lapply(anovaresults, function(x) { x["Pr(>F)"][1:3,] })    )

# Get the genes with good region effect pvalues.
reg.hi.p <-t(data.frame(pvalues[2, pvalues[2,] < 0.0001 & pvalues[3,] > 0.1]))
reg.hi.pdata <- sdata[ row.names(reg.hi.p), ]
``````

A significant p-value resulting from a 1-way ANOVA test would indicate that a gene is differentially expressed in at least one of the groups analyzed. Now that there are more than two groups being analyzed, however, the 1-way ANOVA does not specifically indicate which pair of groups exhibits statistical differences. I know that Post Hoc tests can be applied in this specific situation to determine which specific pair/pairs are differentially expressed in each of the regions ( irrespective of the strains). I would like to know how to apply the Tukey's HSD using R in this case to find out which of these genes ( the ones with good region effect pvalues) are expressed in which region ( for instance in which brain region like "ag","cb","cx and so on).

I tried doing it this way but generates the following error:

``````                  mcps <- TukeyHSD(anovaresults, ordered=TRUE)
Error in UseMethod("TukeyHSD") :  no applicable method for 'TukeyHSD' applied to an object of class "list"
``````
R differential-expression • 4.2k views
modified 6.7 years ago • written 6.7 years ago by deschisandy60
2
David W4.7k wrote:

see `?TukeyHSD`. The function can take an anova fit (as returned by `aov`) but not a list or a an ANOVA table which is what you have in your list. You'll need to edit the function you're using to run your ANOVAs so that they return the `aov` object.

Note though, a p-value < 0.05 doesn't guarantee you'll have significant pairwise-comparisons. Also, these aren't "one way" ANOVAs -- you have two predictor-variables.

I noticed that it is a list object but am unsure how to create a anova fit now.