[BioC] no significant differences in gene expression-limma

[Gordon K Smyth]

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.

> Date: Fri, 21 Feb 2014 01:32:06 +0000
> From: "Sabet, Julia A" <julia.sabet at tufts.edu>
> To: "bioconductor at r-project.org" <bioconductor at r-project.org>
> Subject: [BioC] no significant differences in gene expression-limma
>
> 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)

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}