Hi all,
I have just analyzed some Affymetrix Mouse Gene 2.0 ST arrays using limma and I found no significantly differentially expressed genes (according to adjusted p value being <0.05).  I tried a few different filtering methods, including no filter, and I tried to analyze males and females separately (making separate esets) as well as together. I adjusted for batch dates as well.  I have a sample size of 36 (12 per group).  I did the quality control analysis on Affymetrix Expression Console (not in R) and everything passed according to Affy's quick reference guide and then I proceeded to read in the CEL files and analyze using the code below.  Is there any step that I have left out, or anything else I can do to try to find significant differences?  Any quality control steps in R that I should do besides what I've done already, or should I try a different multiple comparison method besides BH?  My advisor refuses to believe that there could be no significant differences...
Thank you for your help.
Julia

#Load libraries
library(limma)
library(oligo)
#Load CEL files from working directory
celFiles <- list.celfiles()
affyRaw <- read.celfiles(celFiles)
#Normalize/background correct
librarypd.mogene.2.0.st)
eset <- rma(affyRaw)
#Save expression set to an output file
write.exprs(eset,file="maleLog2transformedNormalizeddata02202014.txt")
#Adding gene annotation to eset
library(mogene20sttranscriptcluster.db)
my_frame <- data.frame(exprs(eset))
mogene20sttranscriptcluster()
Annot <-
data.frame(ACCNUM=sapply(contents(mogene20sttranscriptclusterACCNUM), paste, collapse=", "),
SYMBOL=sapply(contents(mogene20sttranscriptclusterSYMBOL), paste, collapse=", "),
DESC=sapply(contents(mogene20sttranscriptclusterGENENAME), paste, collapse=", "))
all <- merge(Annot, my_frame, by.x=0, by.y=0, all=T)
write.table(all,file="maleannotateddata02202014.txt",sep="\t")
#Filter background (very low expression) probes
library(genefilter)
ind <- genefilter(eset, filterfun(kOverA(12, 6)))
eset.filt <- eset[ind,]
dim(eset.filt)
#Or filter by using antigenomic probesets as a measure of
background(filtered more genes out)
librarypd.mogene.2.0.st)
con <- dbpd.mogene.2.0.st)
dbGetQuery(con, "select * from type_dict;")
#type of probes in my array:
table(dbGetQuery(con, "select type from featureSet;")[,1])
antigm <- dbGetQuery(con, "select meta_fsetid from core_mps inner
join featureSet on core_mps.fsetid=featureSet.fsetid where
featureSet.type='4';")
bkg <- apply(exprs(eset)[as.character(antigm[,1]),], 2, quantile, probs=0.95)
minval <- max(bkg)
ind <- genefilter(eset, filterfun(kOverA(12, minval)))
eset.filtB<-eset[ind,]
dim(eset.filtB)
#Filter out control probes and background probes
library(affycoretools)
eset.filtC <- getMainProbes(eset)
library(genefilter)
ind <- genefilter(eset.filtC, filterfun(kOverA(12, 6)))
eset.dblfiltmale <- eset.filtC[ind,]
dim(eset)
dim(eset.filtC)
dim(eset.dblfilt)
#Read targets file into R
targets <- readTargets("Targets.txt", row.names="Name")
#Make design matrix with only diet as covariate
f <- factor(targets$PaternalDiet, levels=c("c","d","s"))
design <- model.matrix(~0+f)
colnames(design)<-c("c","d","s")
fit <- lmFit(eset.filtB, design)
#Make contrast matrix and models
contrast.matrix <- makeContrasts(d-c, s-d, s-c, levels=design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
#Show most significant differentially expressed genes- do for each coefficient separately
topTable(fit2, coef=1, adjust="BH")
#Show top F stats(sig. differences in any contrast)
topTableF(fit2, number=30)
# Make design matrix for 2 factors (sex and diet), adjusted for batch(date)
DS <- paste(targets$PaternalDiet, targets$Sex, sep=".")
DS<-factor(DS,
levels=c("c.female","d.female","s.female","c.male","d.male","s.male"))
design <- model.matrix(~0+DS+Date, targets)
colnames(design) <- gsub("targets", "", colnames(design))
colnames(design)[7:9] <- paste0("Date", 1:3)
fit <- lmFit(eset.dblfilt, design)
cont.matrix <- makeContrasts(
DvsCinMale=DSd.male-DSc.male,
SvsCinMale=DSs.male-DSc.male,
SvsDinMale=DSs.male-DSd.male,
DvsCinFemale=DSd.female-DSc.female,
SvsCinFemale=DSs.female-DSc.female,
SvsDinFemale=DSs.female-DSd.female,
levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)
topTable(fit2, coef=1, adjust="BH")
topTableF(fit2, number=30)

Dear Julia,

Have you made some plots of your data?  For example, plotMDS().  One should always do this, and it's not the same as an Affy Expression Console quality analysis.

Do the treated and untreated groups look different in an MDS plot?  If they do, there you will certainly get DE.  If they don't, then the plot will show you which samples are not doing what you think they should, and then you can follow up why this might be so.

Best wishes
Gordon

PS. There's no rule that says you must use FDR < 0.05.  Larger cutoffs can be used.

Dear Julia,

I see that you have read the help page for plotMDS(), which is good, and that you have tried to copy the usage line given in the help page into your R session.

This isn't at all necessary.  In practice, you just type

   plotMDS(eset.dblfiltmale)

perhaps adding a few arguments if you want to customize the plot.

I can see that the R style of documentation can be a bit counter-intuitive, but the purpose of the usage line on the help page is to describe the defaults for all the arguments.  The reason that the defaults are set is so that you don't have to set them yourself.  To get an idea of how the function is used in practice it is best to look at the example code given at the bottom of the help page or else at the case studies in the limma User's Guide.

Best wishes
Gordon