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
library(pd.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)
library(pd.mogene.2.0.st)
con <- db(pd.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)

	[[alternative HTML version deleted]]

