###################################################
### chunk number 1: init
###################################################
options(width=50)
library(lattice)
library(xtable)


###################################################
### chunk number 2: 
###################################################
library(methylumi)
samps <- read.table(system.file("extdata/samples.txt",
                                package = "methylumi"),sep="\t",header=TRUE)
mldat <- methylumiR(system.file('extdata/exampledata.samples.txt',package='methylumi'),
                    qcfile=system.file('extdata/exampledata.controls.txt',package="methylumi"),
                    sampleDescriptions=samps)


###################################################
### chunk number 3: 
###################################################
mldat


###################################################
### chunk number 4: 
###################################################
getAssayDataNameSubstitutions()


###################################################
### chunk number 5: 
###################################################
md <- cmdscale(dist(t(exprs(mldat)[fData(mldat)$CHROMOSOME=='X',])),2)


###################################################
### chunk number 6: 
###################################################
plot(md,pch=c('F','M')[pData(mldat)$Gender],col=c('red','blue')[pData(mldat)$Gender])


###################################################
### chunk number 7: 
###################################################
avgPval <- colMeans(pvals(mldat))
par(las=2)
barplot(avgPval,ylab="Average P-Value")


###################################################
### chunk number 8: 
###################################################
controlTypes(mldat)


###################################################
### chunk number 9: 
###################################################
qcplot(mldat,"FIRST HYBRIDIZATION")


###################################################
### chunk number 10: 
###################################################
toKeep <- (avgPval<0.05)
pData(mldat)$Gender[9] <- "F"
mldat.norm <- normalizeMethyLumiSet(mldat[,toKeep])


###################################################
### chunk number 11: 
###################################################
library(limma)
dm <- model.matrix(~1+Gender,data=pData(mldat.norm))
colnames(dm)
fit1 <- lmFit(exprs(mldat.norm),dm)
fit2 <- eBayes(fit1)
tt <- topTable(fit2,coef=2,genelist=fData(mldat.norm)[,c('SYMBOL','CHROMOSOME')],number=1536)
x <- aggregate(tt$adj.P.Val,by=list(tt$CHROMOSOME),median)
colnames(x) <- c('Chromosome','Median adjusted P-value')


###################################################
### chunk number 12: 
###################################################
library(xtable)
xt <- xtable(x,label="tab:chromosomepvals",caption="The median adjusted p-value for each chromosome showing that the X-chromosome is highly significantly different between males and females")
digits(xt) <- 6
print(xt,include.rownames=FALSE,align="cr")


###################################################
### chunk number 13: 
###################################################
pdf('gender_probes_by_chrom.pdf')
print(xyplot(-log10(adj.P.Val)~CHROMOSOME,
       tt,ylab="-log10(Adjusted P-value)",
       main="P-values for probes\ndistinguising males from females"))
invisible(dev.off())


###################################################
### chunk number 14: 
###################################################
toLatex(sessionInfo())