## Test script for course Introduction to R and Bioconductor for gene
## expression analysis
##
## Run after RBioCOct2010_install.R
##
## Alexander.Ploner@ki.se 2010-10-12

WORKDIR = tempdir()
setwd(WORKDIR)

datadir = system.file("extdata", package="estrogen")
setwd(datadir)

library(affy)   # Affymetrix pre-processing
library(limma)  # two-color pre-processing; differential expression

phenoData <- read.AnnotatedDataFrame("estrogen.txt", sep="")
eset <- justRMA(filenames=rownames(pData(phenoData)), phenoData=phenoData)

design <- model.matrix(~ estrogen, pData(eset)) # describe model to be fit
fit <- lmFit(eset, design)  # fit each probeset to model
efit <- eBayes(fit)        # empirical Bayes adjustment
topTable(efit, coef=2)      # table of differentially expressed probesets

allest = topTable(efit, coef=2, num=nrow(eset))
plot(logFC~AveExpr, data=allest, pch=".")
topest = allest[1:451, ]
points(logFC~AveExpr, data=topest)
plot(-log10(adj.P.Val)~logFC, data=allest, pch=".")
points(-log10(adj.P.Val)~logFC, data=topest)

library(annaffy)
anncols = aaf.handler()[c(1:3, 8:9, 11:13)]
anntable = aafTableAnn(topest$ID, "hgu95av2.db", anncols)
stattab  = aafTable(logFC=topest$logFC, tstat=topest$t, adjPval=topest$adj.P.Val)
tab = merge(anntable, stattab)
saveHTML(tab, "test.html", title = "Estrogen Coef=2")


estrogen = ReadAffy(filenames=rownames(pData(phenoData)), phenoData=phenoData)
estrogen

RNAdeg = AffyRNAdeg(estrogen)
summaryAffyRNAdeg(RNAdeg)

library(simpleaffy)
QCest = qc(estrogen)
avbg(QCest) 
SFS = sfs(QCest)
SFS
names(SFS) = rownames(pData(phenoData))
round(outer(SFS, SFS, "/"),2)
percent.present(QCest)
ratios(QCest)

library(affyPLM)
PLMest = fitPLM(estrogen)
image(PLMest, which=1, type="resid")
RLE(PLMest)
NUSE(PLMest)

library(arrayQualityMetrics)
rpt <- arrayQualityMetrics(estrogen, outdir=file.path(WORKDIR, "temp"))
rpt2 = arrayQualityMetrics(eset, outdir=file.path(WORKDIR, "temp2"))

library(genefilter)
filt_estrogen <- nsFilter(eset, require.entrez=TRUE, require.GOBP=TRUE, remove.dupEntrez=TRUE, feature.exclude="^AFFX", var.cutoff=0.5)
fit_f <- lmFit(filt_estrogen$eset, design)  # fit each probeset to model
efit_f <- eBayes(fit_f)        # empirical Bayes adjustment
topTable(efit_f, coef=2)      # table of differentially expressed probesets
tt = topTable(efit_f, coef=2, num=nrow(topest))
table(tt$ID %in% topest$ID)