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="Tvs C",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)),E DL933againstneglimmaebresults.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]]

"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="T vsC",adjust="BH",number=nrow(ebfit)),EDL933againstneglimmaebresults.tx t) > > 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 >