### R code from vignette source 'vignettes/epigenomix/inst/doc/epigenomix.Rnw'

###################################################
### code chunk number 1: epigenomix.Rnw:35-37
###################################################
options(width=60)
options(continue=" ")


###################################################
### code chunk number 2: GEdata
###################################################
library(epigenomix)
data(eSet)
pData(eSet)


###################################################
### code chunk number 3: ChIPdata
###################################################
data(mappedReads)
names(mappedReads)


###################################################
### code chunk number 4: ProbeToTranscript
###################################################
probeToTrans <- fData(eSet)$transcript
probeToTrans <- strsplit(probeToTrans, ",")
names(probeToTrans) <- featureNames(eSet)


###################################################
### code chunk number 5: TranscriptToTSS
###################################################
data(transToTSS)
head(transToTSS)


###################################################
### code chunk number 6: BiomaRt (eval = FALSE)
###################################################
## library("biomaRt")
## transcripts <- unique(unlist(transToTSS))
## mart <- useMart("ensembl", dataset="mmusculus_gene_ensembl")
## transToTSS <- getBM(attributes=c("ensembl_transcript_id",
##     "chromosome_name", "transcript_start",
##     "transcript_end", "strand"),
##     filters="ensembl_transcript_id",
##     values=transcripts, mart=mart)


###################################################
### code chunk number 7: dataMatching
###################################################
promoters <- matchProbeToPromoter(probeToTrans,
    transToTSS, promWidth=6000, mode="union")
promoters[["10345616"]]


###################################################
### code chunk number 8: makeChIPseqSet
###################################################
chipSetRaw <- summarizeReads(mappedReads, promoters, summarize="add")
chipSetRaw
head(chipVals(chipSetRaw))


###################################################
### code chunk number 9: normalizeData
###################################################
chipSet <- normalizeChIP(chipSetRaw, method="quantile")


###################################################
### code chunk number 10: normalizationPlot
###################################################
par(mfrow=c(1,2))
plot(chipVals(chipSetRaw)[,1], chipVals(chipSetRaw)[,2],
     xlim=c(1,600), ylim=c(1,600), main="Raw")
plot(chipVals(chipSet)[,1], chipVals(chipSet)[,2],
     xlim=c(1,600), ylim=c(1,600), main="Quantile")


###################################################
### code chunk number 11: integrateData
###################################################
eSet$CEBPA
sampleNames(chipSet)
chipSet$CEBPA <- factor(c("wt", "ko"))
pData(chipSet)

intData <- integrateData(eSet, chipSet, 
    factor="CEBPA", reference="wt")
head(intData)


###################################################
### code chunk number 12: mlMixModel1
###################################################
mlmm = mlMixModel(intData[,"z"],
    normNull=c(2, 3), expNeg=1, expPos=4, 
    sdNormNullInit=c(0.5, 1), rateExpNegInit=0.5, rateExpPosInit=0.5,
    pi=rep(1/4, 4))


###################################################
### code chunk number 13: mlMixModel2
###################################################
mlmm


###################################################
### code chunk number 14: mlMixModelPlot
###################################################
par(mfrow=c(1,2))
plotComponents(mlmm, xlim=c(-2, 2), ylim=c(0, 3))
plotClassification(mlmm)


###################################################
### code chunk number 15: bayesMixModel1
###################################################
set.seed(1515)
bayesmm = bayesMixModel(intData[,"z"],
    normNull=c(2, 3), expNeg=1, expPos=4, 
    sdNormNullInit=c(0.5, 1), rateExpNegInit=0.5, rateExpPosInit=0.5,
    shapeNorm0=c(10, 10), scaleNorm0=c(10, 10), shapeExpNeg0=0.01,
    scaleExpNeg0=0.01, shapeExpPos0=0.01, scaleExpPos0=0.01,
    pi=rep(1/4, 4), itb=2000, nmc=8000, thin=5)


###################################################
### code chunk number 16: bayesMixModel2
###################################################
bayesmm


###################################################
### code chunk number 17: bayesMixModelPlot
###################################################
par(mfrow=c(1,2))
plotComponents(bayesmm, xlim=c(-2, 2), ylim=c(0, 3))
plotClassification(bayesmm, method="mode")


###################################################
### code chunk number 18: compMixModels
###################################################
table(classification(mlmm, method="maxDens"),
      classification(bayesmm, method="mode"))


###################################################
### code chunk number 19: annotation (eval = FALSE)
###################################################
## posProbes <- rownames(intData)[classification(bayesmm, method="mode") == 4]
## library("mogene10sttranscriptcluster.db")
## unlist(mget(posProbes, mogene10sttranscriptclusterSYMBOL))