### 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))