[BioC] bug info

Sat Jun 11 01:34:20 CEST 2011

"doesn't work" is not a useful description of the situation.  give the
exact error message and
the exact code used.  what you have will not work at all as you do not
quote what seems to be a filename.
also give the result of sessionInfo() as noted in the posting guide

On Thu, Jun 9, 2011 at 12:09 PM, Mathieu Olivier <omathieu at lallemand.com> wrote:
> Hi,
>
> We are having problems writing a csv file of our microarray results. I used to do always the same tutorial and no problem. But now I changed computer and downloaded R 2.13.0 and my last step doesn't work anymore.
>
> Here is the command I copied in R : write.csv(topTable(ebfit,coef="TvsC",adjust="BH",number=nrow(ebfit)),EDL933againstneglimmaebresults.txt)
>
> Here is my tutorial and I used it often no problem until now. So can you guys help me ?
>
> Thank you very much and have a nice day.
>
> Olivier
> ####################################################################
>
>
> #Analysis of microarray results from tutorial
>
>
> ####################################################################
>
>
>
> #1
>
> ####################################################################
>
> # Inform R where to take and read the files
>
> # Using File, directory and choose folder
>
> ####################################################################
>
> Scroll down menu by command line
>
> File
>
> Change working directory
>
>
> #2
>
> ####################################################################
>
> #   First download the limma package from Bioconductor
>
> #  Do this only once
>
> ####################################################################
>
> source("http://bioconductor.org/biocLite.R")
>
> biocLite("limma")
>
>
> #3
>
> ####################################################################
>
> # Load limma library in R
>
> ####################################################################
>
>
> library(limma)
>
>
> #4 OPTIONAL
>
> ####################################################################
>
> # For more informations about limma package (optional)
>
> ####################################################################
>
>
> limmaUsersGuide()
>
>
> #5
>
> ####################################################################
>
> # Read the file with experiment design
>
> # The file must be previously created in Excel and saved in tab delimited format
>
> # Columns name are very important as well as capital letters
>
> ####################################################################
>
>
> targets <- readTargets("infoEDLvsnegarrayinv.txt")
>
>
> targets
>
>
> #6
>
> ####################################################################
>
> # Read the files with "read,maimages" function in limma
>
> # To specify where the target vector is in the FileName column summarized
>
> #  by targets$FileName and by specifying which column to read in the different files.
>
> #  Explanations of the command line:
>
> ####################################################################
>
>
> RG <- read.maimages(targets,source="quantarray",wt.fun=wtIgnore.Filter)
>
>
> attributes(RG)  # i.e. what is inside this object?
>
>
> #7
>
> ####################################################################
>
> # Read the newly created file to check if everything is all right
>
> ####################################################################
>
>
> show(RG)
>
>
> #8
>
> ####################################################################
>
> # Look at the 5 first lines of RG file
>
> ####################################################################
>
>
> summary(RG$R)
>
>
> RG[1:4,]
>
>
> #9     OPTIONAL
>
> ####################################################################
>
> # Read Gal file gene list and verify if info is there with: names(RG$genes)
>
> #  If answer is NULL, do the following command
>
> ####################################################################
>
>
> names(RG$genes)
>
>
> #10
>
> ####################################################################
>
> # Read gene list from Gal file
>
> ####################################################################
>
>
> RG$genes<-readGAL()
>
> names(RG$printer)
>
> RG$printer<-getLayout(RG$genes)
>
> names(RG$printer)
>
>
> #11
>
> ####################################################################
>
> # Verify gene frequency on the array
>
> ####################################################################
>
>
> freqTab = table(RG$genes$Name)
>
> head(freqTab)
>
>
> #12
>
> ####################################################################
>
> #  Verify signal intensity in channel R
>
> ####################################################################
>
>
> par(mfrow=c(3,2))
>
> for(i in 1:6){imageplot(log2(RG$R[,i]),RG$printer,low="white",high="red")}
>
>
> #13
>
> ####################################################################
>
> #  Verify signal intensity in channel G
>
> ####################################################################
>
>
> par(mfrow=c(3,2))
>
> for(i in 1:6){imageplot(log2(RG$G[,i]),RG$printer,low="white",high="green")}
>
>
> #14
>
> ####################################################################
>
> #  Verify background intensity in channel R
>
> ####################################################################
>
>
> par(mfrow=c(3,2))
>
>
> for(i in 1:6){imageplot(log2(RG$Rb[,i]),RG$printer,low="white",high="red")}
>
>
> #15
>
> ####################################################################
>
> #  Verify background intensity in channel G
>
> ####################################################################
>
>
> par(mfrow=c(3,2))
>
>
> for(i in 1:6){imageplot(log2(RG$Gb[,i]),RG$printer,low="white",high="green")}
>
>
> #16
>
> ####################################################################
>
> # Do a graph to verify correlation between background and spot signal in channel R
>
> ####################################################################
>
>
> par(mfrow=c(2,3))
>
> for(i in 1:6){ plot(log2(RG$Rb[,i]),log2(RG$R[,i]))}
>
>
> #17
>
> ####################################################################
>
> # Do a graph to verify correlation between background and spot signal in channel G
>
> ####################################################################
>
>
> par(mfrow=c(2,3))
>
> for(i in 1:6){ plot(log2(RG$Gb[,i]),log2(RG$G[,i]))}
>
>
> #18
>
> ####################################################################
>
> # It's important to identify blanks and empty spots
>
> #  SPECIFY A DIFFERENT SYMBOL FOR EACH IN M vs. A PLOTS
>
> ####################################################################
>
>
> spottypes <- readSpotTypes("spottype.txt")
>
>
> spottypes
>
>
> RG$genes$Status <- controlStatus(spottypes,RG)
>
>
> #19
>
> ####################################################################
>
> # Fluorescence boxplot in R channel
>
> ####################################################################
>
>
> par(mfrow=c(1,1))
>
> boxplot(RG$R~col(RG$R),names=colnames(RG$R),xlab="Array",
>
> ylab="fluorescence",main="Slide red fluorescence")
>
>
> #20
>
> ####################################################################
>
> # Fluorescence boxplot in G channel
>
> ####################################################################
>
>
> par(mfrow=c(1,1))
>
> boxplot(RG$G~col(RG$G),names=colnames(RG$G),xlab="Array",
>
> ylab="fluorescence",main="Slide green fluorescence")
>
>
>
> #21
>
> ####################################################################
>
> # Two channels intensity graph
>
> ####################################################################
>
>
> plotDensities(RG)
>
>
> #22
>
> ####################################################################
>
> # M vs A graph before normalization
>
> ####################################################################
>
>
> MA = normalizeWithinArrays(RG,method="none",bc.method="none")
>
> par(mfrow=c(2,3))
>
>
> for(i in 1:6) {
>
>  plotMA(MA[,i],legend=F)
>
>  abline(0,0,col="blue")
>
> }
>
>
> #23 OPTIONAL, slow down computer
>
> ####################################################################
>
> #  Following plots are optional
>
> #  Same M vs A plots as before but perhaps somewhat more aesthetic
>
> #  using smoothScatter from the geneplotter package
>
> ####################################################################
>
>
> source("http://bioconductor.org/biocLite.R")
>
> biocLite("geneplotter")
>
>
> library(geneplotter)
>
>
> par(mfrow=c(2,4))
>
>
> for(i in 1:8) {
>
>    smoothScatter(MA$A[,i],MA$M[,i],xlab = "A",ylab="M",main = rownames(MA$targets)[i])
>
>    abline(h=0,col="red")
>
> }
>
>
> #24       OPTIONAL
>
> ####################################################################
>
> #  MA plots with global loess normalization and simple background subtraction
>
> ####################################################################
>
>
> MA = normalizeWithinArrays(RG,method="loess",bc.method="subtract")
>
> par(mfrow=c(2,4))
>
>
> for(i in 1:5) {
>
>  plotMA(MA[,i],legend=F)
>
>  abline(0,0,col="blue")
>
> }
>
>
> #26
>
> ####################################################################
>
> # Global Lowess normalization without background subtraction
>
> ####################################################################
>
>
> par(mfrow=c(2,3))
>
> MA = normalizeWithinArrays(RG,method="loess",bc.method="none")
>
> for(i in 1:6) {
>
>  plotMA(MA[,i],legend=F)
>
>  abline(0,0,col="blue")
>
> }
>
>
> #27
>
> ####################################################################
>
> # Two channels intensity graph with normalization
>
> ####################################################################
>
>
> par(mfrow=c(1,1))
>
> plotDensities(MA)
>
>
> #28
>
> ####################################################################
>
> # Between array normalization
>
> ####################################################################
>
>
> MA = normalizeBetweenArrays(MA,method="Aquantile")
>
>
> #29
>
> ####################################################################
>
> # Verify the impact of normalization between arrays with another two
>
> # channel intensity graph
>
> ####################################################################
>
>
> plotDensities(MA)
>
>
> #30
>
> ####################################################################
>
> #  Boxplots of the M-values after normalization
>
> ####################################################################
>
>
> par(mfrow=c(1,1))
>
> boxplot(MA$M~col(MA$M),names=colnames(MA$M),xlab="Array",
>
> ylab="M-values",main="Slide specific M-values after aquantile normalization")
>
>
> #31
>
> ####################################################################
>
> # To perform statistical test, blank and empty must be removed
>
> ####################################################################
>
>
> MAfilter=MA[MA$genes$Name!="BLANK"&MA$genes$Name!="EMPTY",]
>
> MAfilter = MAfilter[order(MAfilter$genes$Name),]
>
> show(MAfilter)
>
>
> #32
>
> ####################################################################
>
> # Next the M and A values are calculated for each gene duplicate on the chip
>
> ####################################################################
>
>
> MAave=avereps(MAfilter,ID=MAfilter$genes$Name)
>
> show(MAave)
>
>
> #33
>
> ####################################################################
>
> # Matrix design to analyse the results specifying all the (log) ratios against control
>
> ####################################################################
>
>
> design <-modelMatrix(targets,ref="control")
>
> design  # incomplete design matrix sans dye effet
>
> designfull=cbind(1,design)  # spécification de l'intercept pour le modèle avec l'effet dye.
>
> colnames(designfull) = c("dye","TvsC")  # donne le nom des effet.
>
> designfull
>
>
> #34
>
> ####################################################################
>
> # Calulate model and generate a statistical test regular t-test for each gene
>
> ####################################################################
>
>
> fit <- lmFit(MAave,designfull)
>
> show(fit)
>
>
> #  create the regular linear model t-tests and corresponding P-values
>
> ordinary.t <- fit$coef / fit$stdev.unscaled / fit$sigma
>
> rawP = 2*pt(abs(ordinary.t),df=fit$df.residual,lower=F)
>
>
> #35
>
> ####################################################################
>
> #  Side-by-side histograms for the P-values of Dye and Treatment effects
>
> ####################################################################
>
>
> par(mfrow=c(1,2))
>
> hist(rawP[,1],xlab=" t-test P-values",main="Dye")
>
> hist(rawP[,2],xlab=" t-test P-values",main="Trt")
>
>
> #36
>
> ####################################################################
>
> #  Provide a shrinkage or empirical Bayes analysis for the data using LIMMA
>
> ####################################################################
>
>
> par(mfrow=c(1,1))
>
> ebfit=eBayes(fit)
>
> show(ebfit)
>
>
> #37
>
> ####################################################################
>
> #  Plot the gene-specific variances versus the shrinkage specific variance
>
> ####################################################################
>
>
> plot(ebfit$s2.post,fit$sigma^2,main="Gene-specific versus shrinkage variances")
>
> abline(0,1,col="red")  # add a reference line with 0 intercept and slope 1
>
> abline(h=mean(fit$sigma^2),col="yellow")  # add a horizontal reference line at the average variance
>
>
> #38
>
> ####################################################################
>
> #  Side-by-side histograms for the Shrinkage-Based P-values of Dye and Treatment effects
>
> ####################################################################
>
>
> par(mfrow=c(1,2))
>
> hist(ebfit$p.value[,1],xlab="shrinkage P-values",main="Dye")
>
> hist(ebfit$p.value[,2],xlab="shrinkage P-values",main="Trt")
>
>
> #39
>
> ####################################################################
>
> #  Provide side-by-side "volcano" plots for regular versus shrinkage
>
> #  based tests
>
> ####################################################################
>
>
> par(mfrow=c(1,2))
>
> plot(fit$coef[,2],-log10(rawP[,2]),ylab="-log10(p-value)"
>
> ,xlab="Trt effect on log2 scale",main = "Regular linear model")
>
> plot(ebfit$coef[,2],-log10(ebfit$p.value[,2]),ylab="-log10(p-value)"
>
> ,xlab="Trt effect on log2 scale",main ="Shrinkage analysis")
>
>
> #40
>
> ####################################################################
>
> #  Write out all of the results to a csv file.
>
> ####################################################################
>
>
> write.csv(topTable(ebfit,coef="TvsC",adjust="BH",number=nrow(ebfit)),EDL933againstneglimmaebresults.txt)
>
>
>
> *************************************
> Olivier Mathieu
> Assistant de recherche en biologie moléculaire
> Research assistant in molecular biology
> Institut Rosell-Lallemand
> Montréal, Québec
> P S.V.P., veuillez considérer l'Environnement avant d'imprimer ce courriel
>
>
>        [[alternative HTML version deleted]]
>
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
>