Question: Anova P-Values
2
gravatar for Diana
6.0 years ago by
Diana770
Germany
Diana770 wrote:

Hi everyone,

I'm trying to do an ANOVA on normalized microarray data. My input file looks like this:

Probe_ID    control.CEL    sample1.CEL    sample2.CEL    sample3.CEL
AB112960    2.273088453    2.475903455    2.094961932    2.331040798
AF004856    2.466043833    2.809889315    2.257952048    2.492190838

I have replicates for each (control, sample1, sample2, sample3). I have run ANOVA in r using this script:

aovdata <- read.csv("microarray.csv", header = TRUE)
aovdata1 = aovdata[,2:13]
row.names(aovdata1)<-aovdata$Probe_ID
tissue<-gl(4,3,labels=c("control","aoi","ap","pp"))
aof <- function(x) 
  {m<-data.frame(tissue, x); 
  anova(aov(x ~ tissue))}
anovaresults <- apply(aovdata1, 1,aof)
str(anovaresults)

After performing ANOVA, I'm having trouble extracting the Probe Ids along with their p-value. How can I just extract the probe Ids and the corresponding p-value from ANOVA test? When I view summary(anovaresults), it gives me some different kind of information than str(anovaresults) ?

Thanks a lot!!!

R • 4.8k views
ADD COMMENTlink modified 6.0 years ago by zx87546.8k • written 6.0 years ago by Diana770

Is there a reason you're not using limma? That would be the more standard way to analyze microarray data and results in more convenient data.frames. Aside from that, anova returns an object, which is being coerced into a vector by as.vector. The apply command then returns a matrix of those vectors. You might just run the anova on a single probe and compare those results before and after sending them through as.vector(). That should give you an idea of how everything is formatted. Personally, I would store the anova results in something and then return some sort of meaningfully formatted output from it. Then, you would have better control over things.

ADD REPLYlink written 6.0 years ago by Devon Ryan88k

limma is definitively done for that. It also adds some very useful functionalities like FDR adjustment which might be relevant to use as you run multiple tests. http://www.statsci.org/smyth/pubs/limma-biocbook-reprint.pdf

ADD REPLYlink written 6.0 years ago by Manu Prestat3.9k
2
gravatar for zx8754
6.0 years ago by
zx87546.8k
London
zx87546.8k wrote:

If you are after Pvalues only then try something like this:

aof <- function(x) {
m<-data.frame(tissue, x); 
anova(aov(x ~ tissue)))[,"Pr(>F)"]
}

http://stackoverflow.com/questions/3366506/extract-p-value-from-aov

ADD COMMENTlink modified 6.0 years ago • written 6.0 years ago by zx87546.8k

Thanks a lot!!! but I'm getting this error: Error in summary(anova(aov(x ~ tissue)))[[1]][["Pr(>F)"]] : subscript out of bounds

ADD REPLYlink written 6.0 years ago by Diana770

I updated the answer, can you try again?

ADD REPLYlink written 6.0 years ago by zx87546.8k

Thanks. i think there is some problem with the code I have written because the function you've given me works but the output is like this: 0.184097265 0.908233888 0.359941167 0.975710791 0.774837269 NA NA NA NA NA with numbers in 1 row and NAs in 2nd row and its not loading completely. I have 38535 genes and 12 columns only. I don't understand why its doing that. Its not giving me the output as gene IDs in 1st column and p-values in 2nd column

ADD REPLYlink written 6.0 years ago by Diana770
1

The NAs are from the "Residuals" row (try running the anova on a single row and you'll see this). Just use [1, "Pr(>F)"] instead, keeping in mind that this could fail if you start making more complicated designs. The function, as written, will never give you gene names as output, since the "anova" function doesn't care about the row names and, thus, won't return them. Perhaps returning 'c(rownames(x), anova(lm(x~tissue))[1,"Pr(>F)"])' would work, though I can't say for sure.

ADD REPLYlink written 6.0 years ago by Devon Ryan88k

Thanks !!! it works!

ADD REPLYlink written 6.0 years ago by Diana770
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1183 users visited in the last hour