inst/doc/epigenomix.R

In epigenomix: Epigenetic and gene transcription data normalization and integration with mixture models

### R code from vignette source '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: RNAseqImport
###################################################
data(fpkm)
head(fpkm[c(-2,-8), ])


###################################################
### code chunk number 4: RNAseqToESet
###################################################
mat <- log2(as.matrix(fpkm[, c("CEBPA_WT", "CEBPA_KO")]))
rownames(mat) <- fpkm$tss_id
eSet.seq <- ExpressionSet(mat)
pData(eSet.seq)$CEBPA <- factor(c("wt", "ko"))
fData(eSet.seq)$chr <- fpkm$chr
fData(eSet.seq)$tss <- fpkm$tss


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


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


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


###################################################
### code chunk number 8: BiomaRt (eval = FALSE)
###################################################
## library("biomaRt")
## transcripts <- unique(unlist(probeToTrans))
## mart <- useMart("ENSEMBL_MART_ENSEMBL",dataset="mmusculus_gene_ensembl", host="www.ensembl.org")
## transToTSS <- getBM(attributes=c("ensembl_transcript_id",
##     "chromosome_name", "transcript_start",
##     "transcript_end", "strand"),
##     filters="ensembl_transcript_id",
##     values=transcripts, mart=mart)
## indNeg <- transToTSS$strand == -1
## transToTSS$transcript_start[indNeg] <- transToTSS$transcript_end[indNeg]
## transToTSS$transcript_end <- NULL


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


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


###################################################
### code chunk number 11: makeChIPseqSet_rna
###################################################
promoters.seq <- GRanges(seqnames=fData(eSet.seq)$chr,
                        ranges=IRanges(start=fData(eSet.seq)$tss, width=1),
                        probe=featureNames(eSet.seq))
promoters.seq <- resize(promoters.seq, width=3000, fix="center")
promoters.seq <- split(promoters.seq, elementMetadata(promoters.seq)$probe)


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


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


###################################################
### code chunk number 14: 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 15: integrateData
###################################################
eSet$CEBPA
colnames(chipSet)
chipSet$CEBPA <- factor(c("wt", "ko"))
colData(chipSet)

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


###################################################
### code chunk number 16: 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 17: mlMixModel2
###################################################
mlmm


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


###################################################
### code chunk number 19: 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), sdAlpha=1, itb=2000, nmc=8000, thin=5)


###################################################
### code chunk number 20: bayesMixModel2
###################################################
bayesmm


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


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


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