Hello everyone, I am analysing a study in GEO database with a number:GSE41364, I'm downloading raw data from GEO site. Unfortunately I have stumbled upon a problem. I don't know if I'm doing something wrong or things are just ok, but when it comes to summary results after fitting there seems to be a problem. Intercept shows that there are 0 genes that are down regulated and NotSig. But then EXP2 and EXP3 shows that there are a few thousands gene that are NoSig, but 0 genes that are down regulated and up regulated. But on the other hand the data form 'treatment5 μM of 5-Aza' and 'treatment10 μM of 5-Aza' in result after fitting seems to be just fine.

Thanks for help :)

library(limma)
library(affy)
library(hgu133plus2.db)
library(annotate)

#The data files are contained in the extdata directory of the estrogen package:
datadir = dir()

targets = readTargets(file ='Targets.txt', path = NULL, sep='\t', row.name='FileName')
targets

ab <- ReadAffy(filenames=targets$FileName, celfile.path=NULL)
eset <- rma(ab)
dim(eset)
plotMDS(eset)

ID <- featureNames(eset)
Symbol <- getSYMBOL(ID,"hgu133plus2.db")
fData(eset) <- data.frame(Symbol=Symbol)

#We remove all probes that do not have an Entrez Gene ID and Symbol:
HasSymbol <- !is.na(fData(eset)$Symbol)
eset <- eset[HasSymbol,]
dim(eset)
plotMDS(eset)

#Differential expression
Exp <- factor(targets$Experiment)
treatment <- factor(targets$treatment, levels=c("0 μM of 5-Aza","5 μM of 5-Aza", "10 μM of 5-Aza"))
design <- model.matrix(~Exp+treatment)

#fitting 
fit <- lmFit(eset, design)
fit <- eBayes(fit, trend=TRUE, robust=TRUE)
results <- decideTests(fit)
summary(results)

#Output of my code:

datadir = dir()
> targets = readTargets(file ='Targets.txt', path = NULL, sep='\t', row.name='FileName')
> targets
                                             FileName      treatment Experiment
GSM1015541_HT29_P12_D1.CEL GSM1015541_HT29_P12_D1.CEL  0 μM of 5-Aza          1
GSM1015542_HT29_P12_D2.CEL GSM1015542_HT29_P12_D2.CEL  0 μM of 5-Aza          2
GSM1015543_HT29_P12_D3.CEL GSM1015543_HT29_P12_D3.CEL  0 μM of 5-Aza          3
GSM1015544_HT29_P12_A1.CEL GSM1015544_HT29_P12_A1.CEL  5 μM of 5-Aza          1
GSM1015545_HT29_P12_A2.CEL GSM1015545_HT29_P12_A2.CEL  5 μM of 5-Aza          2
GSM1015546_HT29_P12_A3.CEL GSM1015546_HT29_P12_A3.CEL  5 μM of 5-Aza          3
GSM1015547_HT29_P14_B2.CEL GSM1015547_HT29_P14_B2.CEL 10 μM of 5-Aza          1
GSM1015548_HT29_P14_B4.CEL GSM1015548_HT29_P14_B4.CEL 10 μM of 5-Aza          2
GSM1015549_HT29_P14_B5.CEL GSM1015549_HT29_P14_B5.CEL 10 μM of 5-Aza          3

> ab <- ReadAffy(filenames=targets$FileName, celfile.path=NULL)
> eset <- rma(ab)

Background correcting
Normalizing
Calculating Expression


> dim(eset)
Features  Samples 
   54675        9 
> plotMDS(eset)

> ID <- featureNames(eset)
> Symbol <- getSYMBOL(ID,"hgu133plus2.db")
> fData(eset) <- data.frame(Symbol=Symbol)
> HasSymbol <- !is.na(fData(eset)$Symbol)
> eset <- eset[HasSymbol,]
> dim(eset)
Features  Samples 
   43138        9 
> plotMDS(eset)


> Exp <- factor(targets$Experiment)
> treatment <- factor(targets$treatment, levels=c("0 μM of 5-Aza","5 μM of 5-Aza", "10 μM of 5-Aza"))
> design <- model.matrix(~Exp+treatment)
> fit <- lmFit(eset, design)
> fit <- eBayes(fit, trend=TRUE, robust=TRUE)
> results <- decideTests(fit)
> summary(results)
       (Intercept)  Exp2  Exp3 treatment5 μM of 5-Aza treatment10 μM of 5-Aza
Down             0     0     0                   4123                    5264
NotSig           0 43138 43138                  34424                   33781
Up           43138     0     0                   4591                    4093

There isn't any problem here. The results look quite normal. The Intercept doesn't compare any two groups so you don't need to test it equal to zero. You just ignore the Intercept results.

It isn't clear that there were actually three distinct experiments in this study. You could probably just use design <- model.matrix(~treatment).