[BioC] no significant differences in gene expression-limma

[Sabet, Julia A]

Thank you Gordon.  I looked up the R documentation for plotMDS and did the following code, using my normalized eset, and came up with the following error.  Could you or someone else tell me what I did wrong and if I am missing some code?  I am new to R and had never heard of plotMDS before...I don't really understand what dim.plot is.
Thank you!

> mds<-plotMDS(eset.dblfiltmale, top=500, labels=colnames(eset.dblfiltmale), col=NULL, cex=1, dim.plot=c(1,2), ndim=max(dim.plot), gene.selection="pairwise", xlab=paste("Dimension",dim.plot[1]), ylab=paste("Dimension",dim.plot[2]))
Error in plotMDS.default(eset.dblfiltmale, top = 500, labels = colnames(eset.dblfiltmale),  : 
  object 'dim.plot' not found



-----Original Message-----
From: Gordon K Smyth [mailto:smyth at wehi.EDU.AU] 
Sent: Saturday, February 22, 2014 2:41 AM
To: Sabet, Julia A
Cc: Bioconductor mailing list
Subject: no significant differences in gene expression-limma

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:6}}