Agilent microarray data analysis end in no significant genes ?

Question

[DareDevil]

Hi I am analyzing microarray under agilent platform. I used following codes

library(limma)
workDir<- "home/dare_devil/arrays/data1"

annotation <- read.table(file.path(workDir, "annotation.txt"), 
                        header = TRUE, sep = '\t', quote = "",
                        stringsAsFactors = FALSE)
if (length(unique(annotation$Name)) != nrow(annotation)) {stop("The names in the annotation file must be unique!")}

rawData <- read.maimages(unique(file.path(workDir, annotation$File)),
                        source = "agilent", green.only = TRUE,
                        names = annotation$Name, annotation = c("ProbeName"))

# correct, normalize, and extract
bgData <- backgroundCorrect(rawData)
normData <- normalizeBetweenArrays(bgData, method = "quantile")
normEset <- normData$E; rownames(normEset) <- normData$genes$ProbeName

#create design
conditions<- paste(annotation$Group)
conditions <- factor(conditions, levels=unique(conditions))
design <- model.matrix(~0+ conditions)
colnames(design) <- levels(conditions)

#fit the design
fit <- lmFit(normEset, design)
#Make contrast
contrast.matrix <- makeContrasts(Grp1 - Grp2, levels=design)

#Fit contrast
fit.cont<- contrasts.fit(fit, contrast.matrix)
#ebayes Fit
fit.cont<- eBayes(fit.cont)

#check the results
results<-decideTests(fit.cont,adjust.method="fdr",p=0.05)
summary(results)

I got following results

      Grp1 - Grp2
Down             0
NotSig       62948
Up               0

But after running the code I did not get any significant genes. Am I doing anything wrong?

Answer 1

[Kevin Blighe]

I do not see anything majorly wrong.

Is this a single colour (single channel) array? - if 'yes', you could apply some of the extra filtering steps, here: A: How to process (seems) Agilent microarrry data?

What is Grp1 and Grp2 and what are the numbers of samples comprising each?

If you run the following, what are the top genes?

topTable(fit.cont, adjust.method = 'fdr', p = 1)

Kevin

[DareDevil]

Hi Kevin,

#fit the design
fit <- lmFit(normEset, design)

After above step I can perform directly ebayes, right?

fita <- eBayes(fit)

I get:

results<-decideTests(fita,adjust.method="fdr",p=0.05)
summary(results)
        Grp1  Grp2
Down      15     2
NotSig  1715  1694
Up     61221 61276

This is mentioned in the manual limma page 41-43

[Kevin Blighe]

Yes, eBayes() can be used in different ways

[DareDevil]

Can you look at this link. it gives different results

[DareDevil]

Grp1 is treated group and Grp2 is control group and each group has 2 replicates

topTable(fit.cont, adjust.method = 'fdr', p = 1)

    ID     logFC  AveExpr         t      P.Value adj.P.Val        B
49310 A_24_P304071 -5.390904 6.496221 -24.71072 1.554252e-05 0.1294594 2.859211
10623  A_23_P51487 -5.526571 7.021997 -23.98208 1.751334e-05 0.1294594 2.815474
22235 A_32_P452655 -5.437434 9.659100 -22.41388 2.293240e-05 0.1294594 2.709814
13117 A_32_P473302  4.523709 7.848596  20.38747 3.344458e-05 0.1294594 2.545139
25404  A_23_P10025  4.451570 7.280976  20.36668 3.358065e-05 0.1294594 2.543257
4959  A_23_P105862  4.432852 7.064774  20.23012 3.449170e-05 0.1294594 2.530789
6282  A_23_P204286 -4.564843 8.994801 -19.97466 3.628034e-05 0.1294594 2.506958
38709  A_23_P10025  4.529188 7.183408  19.87968 3.697471e-05 0.1294594 2.497926
39696  A_23_P10025  4.558354 7.186030  19.22008 4.228408e-05 0.1294594 2.432489
39786 A_23_P145874 -4.733944 5.160861 -19.16862 4.273711e-05 0.1294594 2.427177

[Kevin Blighe]

I see. The problem is the sample size. With 2 versus 2, you are limiting yourself to what you can actually discover. There was a recent related thread, here: A: RNA-seq power estimation using ssizeRNA program

If this work is just for training purposes, then I would not worry about it too much.

If you want, you could try some very strict / rigorous filtering steps, like I do in the other thread - this will lessen the 'harshness' of the FDR and may give you something in return.