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]]

