Example 1

#########################################################################
#### Practica de LCG para la materia de BioInformatica aplicada      ####
#### Contiene: instrucciones para la instalacion completa  de        ####
#### Bioconductor y un ejemplo de analisis de un data set de 8       ####
#### Affymetrixarrays. Seimplementan y comparan MAS 5.0,             ####
#### GC-RMA (PM-MM model) and GC-RMA (PM only model)                 ####
####                                Rosa Maria Gutierrez 15/02/07    ####
#########################################################################

## Datos bajados de http://www.affymetrix.com/analysis/download_center2.affx
## HUMAN GENOME U95 DATA SET

## instalar Bioconductor
source("http://www.bioconductor.org/biocLite.R")
biocLite()
getBioC("pkg1")  ## Para instalar paquetes sueltos
biocLite(c("pkg1", "pkg2"))## Para instalar paquetes sueltos


## Cargar librerias
library("affy") ## cargar libreria affy package
library("matchprobes")  ## la libreria gcrma necesita la libreria matchprobes 
library("gcrma")  ## cargar libreria  gcrma package

dataDir<-"Affymetrixarrays/data"
workingDir<-"Affymetrixarrays"
fileNames<-dir(dataDir)
fileNames<-paste(dataDir, "/", fileNames, sep="")
LatinData<-ReadAffy(filenames=fileNames)  
## 8 arreglos para el analisis: "1521m99hpp_av06.CEL" "1521n99hpp_av06.CEL"
"1521o99hpp_av06.CEL" "1521p99hpp_av06.CEL" 
## "1521q99hpp_av06.CEL" "1521r99hpp_av06.CEL" "1521s99hpp_av06.CEL"
"1521t99hpp_av06.CEL". Los 4 primeros son un juego de replicas
# y los otros cuatro son otro juego de replicas


## Calcular medidas de expresion por MAS 5.0
LatinData.mas5<-expresso(LatinData, pmcorrect.method="mas", bgcorrect.method="mas", 
summary.method="mas")
LatinData.mas5.norm<-affy.scalevalue.exprSet(LatinData.mas5)  
write.table(exprs(LatinData.mas5.norm), paste(workingDir,
"\\LatinData_mas5_norm.txt", sep=""), sep="\t")


## Calcular medidas de expresion por gcRMA (PM-only)
LatinData.gcrma<-gcrma(LatinData, type="affinities")
write.table(exprs(LatinData.gcrma), paste(workingDir, "\\LatinData_gcrma_noMM.txt",
sep=""), sep="\t")

## Calcular medidas de expresion por gcRMA (contain MM)
LatinData.gcrma.MM<-gcrma(LatinData, type="fullmodel")
write.table(exprs(LatinData.gcrma.MM), paste(workingDir, "\\LatinData_gcrma_MM.txt",
sep=""), sep="\t")

#################
## Definicion de funciones en R

evaluate<-function(exprs.all, names)
{
        par(mfrow=c(3,2))
        for(i in 1:length(exprs.all))
        {
                a<-matrix(NA, nrow(exprs.all[[i]]), 2)
                a[,1]<-apply(exprs.all[[i]], 1, mean, na.rm=TRUE)
                a[,2]<-apply(exprs.all[[i]], 1, sd, na.rm=TRUE)
                plot(c(0, 12), c(0, 3), type="n", xlab="mean log intensities", ylab="sd log
intensities")
                points(a[,1], a[,2], pch=".")
                title(names[i])
        }
}


evaluate1<-function(exprs.all, names)
{
        for(i in 1:length(exprs.all))
        {
                a<-cor(exprs.all[[i]], use="pairwise.complete.obs")
                print(paste(names[i], ":",  mean(a[a<1])))        
        }
}

## Obtener los valores de expresion de las primeras 4 replicas
LatinData.mas5.norm.matrix<-log2(exprs(LatinData.mas5.norm)[,1:4])
LatinData.gcrma.matrix<-exprs(LatinData.gcrma)[,1:4]
LatinData.gcrma.MM.matrix<-exprs(LatinData.gcrma.MM)[,1:4]

## Graficar evaluaciones. 
jpeg(filename=paste(workingDir, "\\eval_replicate_1_new.jpg", sep=""), width=1200,
height=1024)
exprs.all<-list(LatinData.mas5.norm.matrix, LatinData.gcrma.matrix,
LatinData.gcrma.MM.matrix)
evaluate(exprs.all, c("MAS5","GCRMA (no MM)", "GCRMA (MM)"))
dev.off()
evaluate1(exprs.all, c("MAS5", "GCRMA(no MM)", "GCRMA(MM)"))
## print out of evaluate1
##[1] "MAS5 : 0.893039182394166"
##[1] "GCRMA(no MM) : 0.9992816227403"
##[1] "GCRMA(MM) : 0.998767869268427"


## Obtener los valores de expresion de las primeras siguentes 4 replicas
LatinData.mas5.norm.matrix<-log2(exprs(LatinData.mas5.norm)[,5:8])
LatinData.gcrma.matrix<-exprs(LatinData.gcrma)[,5:8]
LatinData.gcrma.MM.matrix<-exprs(LatinData.gcrma.MM)[,5:8]

## Plot evaluation figure again. 
jpeg(filename=paste(workingDir, "\\eval_replicate_2_new.jpg", sep=""), width=1200,
height=1024)
exprs.all<-list(LatinData.mas5.norm.matrix, LatinData.gcrma.matrix,
LatinData.gcrma.MM.matrix)
evaluate(exprs.all, c("MAS5","GCRMA (no MM)", "GCRMA (MM)"))
dev.off()
evaluate1(exprs.all, c("MAS5","GCRMA(no MM)", "GCRMA(MM)"))

## print out of evaluate1
##[1] "MAS5 : 0.900234102618892"
##[1] "GCRMA(no MM) : 0.999446021907957"
##[1] "GCRMA(MM) : 0.998953186433171"