Agilent microarray data analysis end in no significant genes ?
1
0
Entering edit mode
10 months ago
dare_devil ★ 1.5k

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?

microarray r limma • 384 views
ADD COMMENT
1
Entering edit mode
10 months ago

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

ADD COMMENT
0
Entering edit mode

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
ADD REPLY
1
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

Yes, eBayes() can be used in different ways

ADD REPLY
0
Entering edit mode

Can you look at this link. it gives different results

ADD REPLY

Login before adding your answer.

Traffic: 3163 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6