PEC1.R

###################################################
### 0.Preparación de los datos
###################################################

###################################################
### 0.1.Asignación de directorios
###################################################
setwd("/<...>/PEC1")
workingDir <-getwd()
dataDir <-file.path(workingDir, "data")
resultsDir <- file.path(workingDir,"results")
setwd(workingDir)

###################################################
### 0.2.Asignación de opciones
###################################################
options(width=80)
options(digits=5)

###################################################
### 0.3.Datos
###################################################
print(dataDir)
print(resultsDir)
print(workingDir)

###################################################
### 0.4.Preparación de herramientas Biobase/affy y
### asignación de nombres de variables y ficheros
###################################################
require(Biobase)
require(affy)
sampleInfo <- read.AnnotatedDataFrame(file.path(dataDir,"samples2.txt"), header = TRUE, row.names = 1, sep="\t")
fileNames <- rownames(pData(sampleInfo))
rawData <- read.affybatch(filenames=file.path(dataDir,fileNames),phenoData=sampleInfo)

###################################################
### 1.Exploración y visualización
###################################################

###################################################
### 1.0. Datos 
###################################################

### Revisión de que los datos se copiaron
### correctamente
#class(rawData)
#str(rawData)
#rawData
#phenoData(rawData)
#pData(rawData)
#sampleNames(rawData)
#gn<-geneNames(rawData)
#pm(rawData, gn[24])

###################################################
### 1.1. Histograma
###################################################
info <- data.frame(grupo=c(rep(1,3), rep (2,3), rep(3,3), rep(4,3)))
sampleNames <- pData(rawData)$Sample
png("histogram.png")
hist(rawData, main="Signal distribution", col=info$grupo, lty=1:ncol(info))
legend (x="topright", legend=sampleNames, col=info$grupo, lty=1:ncol(info))
dev.off()

###################################################
### 1.2. Boxplot
###################################################
par(mfrow = c(1, 1))
png("boxplot.png")
boxplot(rawData, cex.axis=0.5, col=info$grupo+1, las=2, names=sampleNames, main="Signal distribution for all chips")
dev.off()

###################################################
### 1.3. Gráfico de degradación
###################################################
par(mfrow = c(1, 1))
deg<-AffyRNAdeg(rawData, log.it=T)
summaryAffyRNAdeg(deg) 
png("degradacion.png")
plotAffyRNAdeg(deg,cols=1:12)
legend (x="bottomright", legend=sampleNames, col=1:12, lty=1:nrow(info), cex=0.7)
dev.off()

###################################################
### 1.4. Dendograma
###################################################
clust.euclid.average <- hclust(dist(t(exprs(rawData))),method="average")
png("dendograma.png")
plclust(clust.euclid.average, labels=sampleNames, main="Hierarchical clustering of samples",  hang=-1)
dev.off()





###################################################
### 1.5. Image plot
###################################################
png("image.png")
par(mfrow = c(4, 3))
image(rawData)
dev.off()

###################################################
### 1.6 MAPlot 
###################################################
x11(width=10, height=8)
png("MAplot.png")
MAplot(rawData, pairs = TRUE, plot.method = "smoothScatter")
dev.off()

###################################################
### 2. Control de calidad
###################################################

###################################################
### 2.1. QC Report
###################################################
stopifnot(require(affyQCReport))
QCReport(rawData,file=file.path(resultsDir,"QCReport.pdf"))


###################################################
### 2.2. affyPLM
###################################################
stopifnot(require(affyPLM))
Pset<- fitPLM(rawData)

x11(width=10, height=8)
par(mfrow = c(1, 2))

png("RLE.png")
RLE(Pset, main = "Relative Log Expression", names=sampleNames, las=2, col=info$grupo+1, cex.axis=0.6,ylim=c(-5,5))
dev.off()
png("NSU.png")
NUSE(Pset, main = "Normalized Unscaled Standard Errors", las=2, names=sampleNames, las=2, col=info$grupo+1, cex.axis=0.6, ylim=c(0.5,1.5))
dev.off()






###################################################
### 3. Normalización
###################################################

###################################################
### 3.1. RMA
###################################################
eset_rma <- expresso(rawData,
bgcorrect.method="rma",
normalize.method="quantiles",
pmcorrect.method="pmonly",
summary.method="medianpolish")

###################################################
### 4. Filtraje no específico
###################################################

###################################################
### 4.1 Filtrado con nsFilter
###################################################
if (!require(mouse4302.db)){
source("http://bioconductor.org/biocLite.R")
biocLite("mouse4302.db")
}else{
require("mouse4302.db")
}
require(genefilter)

filteredData <- nsFilter(eset_rma, require.entrez=TRUE,
require.GOBP=FALSE, require.GOCC=FALSE,
require.GOMF=FALSE, require.CytoBand=FALSE,
remove.dupEntrez=TRUE, var.func=IQR,
var.cutoff=0.25, var.filter=TRUE,
filterByQuantile=TRUE, feature.exclude="^AFFX")

### observación de datos filtrados
#filteredData
#names(filteredData)
#class(filteredData$eset)
#print(filteredData$filter.log)

eset_filtered <-filteredData$eset








###################################################
### Salvar datos normalizados (rma) y filtrados
###################################################

save(eset_rma, eset_filtered, file=file.path(resultsDir, "datos.normalizados.Rda"))

###################################################
### 5. Selección de genes diferencialmente
### expresados a través de modelos lineales
###################################################

###################################################
### 5.0. Datos
###################################################
if (!exists("eset_filtered")) load(file.path(resultsDir, "datos.normalizados.Rda"))
#class(eset_filtered)

###################################################
### 5.1. Generación de factores para el modelo
### diseño de la matriz
###################################################
targets <- pData(eset_filtered)
stopifnot(require(limma))
TS <- paste(targets$Factor1, targets$Factor2, sep=".")
TS <- as.factor(TS)
print(TS)

design <- model.matrix(~0+TS)
colnames(design) <- levels(TS)
print(design)

###################################################
### 5.2. Definición de los contrastes a evaluar
###################################################
cont.matrix <- makeContrasts (Age=(Aged.LPS-Aged.Medium)-(Young.LPS-Young.Medium),
                              AgedbyTreat=Aged.LPS-Aged.Medium,
                              levels = design)
print(cont.matrix)

###################################################
### 5.3. Cálculo del modelo lineal
### y extracción de contrastes deseados
###################################################
fit<-lmFit(eset_filtered, design)
fit.main<-contrasts.fit(fit, cont.matrix)
fit.eBayes<-eBayes(fit.main)


###################################################
### 5.3.1 Volcano plot
###################################################
coefnum = 1
png("volcano.png")
opt <- par(cex.lab = 0.7)
volcanoplot(fit.eBayes, coef=coefnum, highlight=10, names=fit.eBayes$ID, main=paste("Differentially expressed genes", colnames(cont.matrix)[coefnum], sep="\n"))
abline(v=c(-1.5,1.5))
abline(h=-log10(0.05))
par(opt)
dev.off()

###################################################
### 5.3 Top tables
###################################################
###################################################
### 5.3.1 Edad (Aged vs Young)
###################################################
topAge <- topTable (fit.eBayes, number=nrow(fit.eBayes), coef="Age", adjust="fdr", p.value=0.05,lfc=1.5)

### Ordenamos según p-value
topAge <-topAge[order(topAge$P.Value,decreasing=TRUE),]
summary(topAge)
mouse4302()

### Escribimos resultados en tabla

selectedIDs<-rownames(topAge)

# convertimos a identificadores Entrez y a Gene Symbols
Symbol.topAge <- unlist(mget(topAge$ID, mouse4302SYMBOL))
EntrezID.topAge <- unlist(mget(topAge$ID,mouse4302ENTREZID))

# Haremos la columna de Entrez sea hiperenlazable
paraEnlace <- list (misgenes=EntrezID.topAge)

# Preparamos el data.frame con el que  crear ́ el archivo de resultados
otherNames = data.frame(selectedIDs, Symbol.topAge, topAge)

names(otherNames) = c("Affy ID", "Gene Symbol", names(topAge))

# Creación de html con los resultados
htmlpage(paraEnlace,
         filename =file.path(resultsDir, "SelectedGenesAge.html") ,
         title = "DEG de la comparacion entre ratones jovenes y viejos ",
         othernames = otherNames,
         table.head = c("Entrez IDs", names(otherNames)),
         table.center = TRUE,
         repository=list("en"))

###################################################
### 5.3.2 Treatment (LPS vs Medium) para Aged
###################################################
topAgedbyTreat  <- topTable (fit.eBayes, number=nrow(fit.eBayes) , coef="AgedbyTreat", adjust="fdr",p.value=0.05,lfc=1.5)

### Escribimos resultados en tabla
selectedIDs<-rownames(topAgedbyTreat)

# convertimos a identificadores Entrez y a Gene Symbols
Symbol.topAgedbyTreat <- unlist(mget(topAgedbyTreat$ID, mouse4302SYMBOL))
Symbol.topAgedbyTreat
EntrezID.topAgedbyTreat <- unlist(mget(topAgedbyTreat$ID,mouse4302ENTREZID))
EntrezID.topAgedbyTreat

# Haremos la columna de Entrez sea hiperenlazable
paraEnlace <- list (misgenes=EntrezID.topAgedbyTreat)
paraEnlace

# Preparamos el data.frame con el que  crear ́ el archivo de resultados
otherNames = data.frame(selectedIDs, Symbol.topAgedbyTreat, topAgedbyTreat)

names(otherNames) = c("Affy ID", "Gene Symbol", names(topAgedbyTreat))

# Creación de html con los resultados
htmlpage(paraEnlace,
         filename =file.path(resultsDir, "SelectedGenesAgedbyTreat.html") ,
         title = "DEG del efecto del LPS en ratones viejos ",
         othernames = otherNames,
         table.head = c("Entrez IDs", names(otherNames)),
         table.center = TRUE,
         repository=list("en"))

###################################################
### 5.4 Decide tests
###################################################
stopifnot(require(annotate))
Symbol.fit.eBayes <- unlist(mget(fit.eBayes$ID, mouse4302SYMBOL))
res<-decideTests(fit.eBayes, method="separate", adjust.method="fdr", p.value=0.05,lfc=1.5)

sum.res.rows<-apply(abs(res),1,sum)
res.selected<-res[sum.res.rows!=0,]
print(summary(res))



###################################################
### 5.3.2 Venn diagram
###################################################
png("vennDiagram.png")
vennDiagram (res.selected[,1:2], main="Genes in common #1", cex=0.9)
dev.off()

###################################################
### 5.3.3 Heatmap plot
###################################################
probeNames<-rownames(res)
probeNames.selected<-probeNames[sum.res.rows!=0]
exprs2cluster <-exprs(res)[probeNames.selected,]

color.map <- function(horas) { if (horas< 20) "yellow" else "red" }
grupColors <- unlist(lapply(pData(res)$Sample, color.map))

x11(width=10, height=10)
png("heatmap.png")
heatmap(exprs2cluster,
        col=rainbow(100),
        ColSideColors=grupColors,
        cexCol=0.9)
dev.off()