inst/doc/mosaics-example.R

In mosaics: MOSAiCS (MOdel-based one and two Sample Analysis and Inference for ChIP-Seq)

### R code from vignette source 'mosaics-example.Rnw'

###################################################
### code chunk number 1: preliminaries
###################################################
options(prompt = "R> ")


###################################################
### code chunk number 2: mosaics-prelim
###################################################
library("mosaics")


###################################################
### code chunk number 3: mosaicsExample-prelim
###################################################
library("mosaicsExample")


###################################################
### code chunk number 4: constructBins
###################################################
constructBins( infile=system.file( file.path("extdata","wgEncodeSydhTfbsGm12878Stat1StdAlnRep1_chr22_sorted.bam"), package="mosaicsExample"), 
    fileFormat="bam", outfileLoc="./", 
    byChr=FALSE, useChrfile=FALSE, chrfile=NULL, excludeChr=NULL, 
    PET=FALSE, fragLen=200, binSize=200, capping=0 )
constructBins( infile=system.file( file.path("extdata","wgEncodeSydhTfbsGm12878InputStdAlnRep1_chr22_sorted.bam"), package="mosaicsExample"), 
    fileFormat="bam", outfileLoc="./", 
    byChr=FALSE, useChrfile=FALSE, chrfile=NULL, excludeChr=NULL, 
    PET=FALSE, fragLen=200, binSize=200, capping=0 )


###################################################
### code chunk number 5: io-readbin
###################################################
fileName <- c(
    "wgEncodeSydhTfbsGm12878Stat1StdAlnRep1_chr22_sorted.bam_fragL200_bin200.txt",
    "wgEncodeSydhTfbsGm12878InputStdAlnRep1_chr22_sorted.bam_fragL200_bin200.txt")

binTFBS <- readBins( type=c("chip","input"), fileName=fileName )


###################################################
### code chunk number 6: io-bindata-show
###################################################
binTFBS


###################################################
### code chunk number 7: io-bindata-print
###################################################
print(binTFBS)[ 90600:90614, ]


###################################################
### code chunk number 8: io-bindata-plot (eval = FALSE)
###################################################
## plot( binTFBS )
## plot( binTFBS, plotType="input" )


###################################################
### code chunk number 9: fig-bindata-plot-hist
###################################################
plot(binTFBS)


###################################################
### code chunk number 10: fig-bindata-plot-input
###################################################
plot( binTFBS, plotType="input" )


###################################################
### code chunk number 11: io-mosaicsfit
###################################################
fitTFBS <- mosaicsFit( binTFBS, analysisType="IO", bgEst="rMOM" )


###################################################
### code chunk number 12: io-mosaicsfit-show
###################################################
fitTFBS


###################################################
### code chunk number 13: io-mosaicsfit-plot (eval = FALSE)
###################################################
## plot(fitTFBS)


###################################################
### code chunk number 14: fig-mosaicsfit-plot
###################################################
plot(fitTFBS)


###################################################
### code chunk number 15: ts-mosaicspeak
###################################################
peakTFBS <- mosaicsPeak( fitTFBS, signalModel="2S", FDR=0.05, 
maxgap=200, minsize=50, thres=10 )


###################################################
### code chunk number 16: io-mosaicspeak-show
###################################################
peakTFBS


###################################################
### code chunk number 17: io-mosaicspeak-print
###################################################
head(print(peakTFBS))


###################################################
### code chunk number 18: io-mosaicspeak-export (eval = FALSE)
###################################################
## export( peakTFBS, type="txt", filename="peakTFBS.txt" )
## export( peakTFBS, type="bed", filename="peakTFBS.bed" )
## export( peakTFBS, type="gff", filename="peakTFBS.gff" )
## export( peakTFBS, type="narrowPeak", filename="peakTFBS.narrowPeak" )
## export( peakTFBS, type="broadPeak", filename="peakTFBS.broadPeak" )


###################################################
### code chunk number 19: generateWig
###################################################
generateWig( infile=system.file( file.path("extdata","wgEncodeSydhTfbsGm12878Stat1StdAlnRep1_chr22_sorted.bam"), package="mosaicsExample"), 
    fileFormat="bam", outfileLoc="./", 
    byChr=FALSE, useChrfile=FALSE, chrfile=NULL, excludeChr=NULL, 
    PET=FALSE, fragLen=200, span=200, capping=0, normConst=1 )


###################################################
### code chunk number 20: mosaicsFitHMM-pre
###################################################
constructBins( infile=system.file( file.path("extdata","wgEncodeBroadHistoneGm12878H3k4me3StdAlnRep1_chr22_sorted.bam"), package="mosaicsExample"), 
    fileFormat="bam", outfileLoc="./", 
    byChr=FALSE, useChrfile=FALSE, chrfile=NULL, excludeChr=NULL, 
    PET=FALSE, fragLen=200, binSize=200, capping=0 )
constructBins( infile=system.file( file.path("extdata","wgEncodeBroadHistoneGm12878ControlStdAlnRep1_chr22_sorted.bam"), package="mosaicsExample"), 
    fileFormat="bam", outfileLoc="./", 
    byChr=FALSE, useChrfile=FALSE, chrfile=NULL, excludeChr=NULL, 
    PET=FALSE, fragLen=200, binSize=200, capping=0 )

fileName <- c(
    "wgEncodeBroadHistoneGm12878H3k4me3StdAlnRep1_chr22_sorted.bam_fragL200_bin200.txt",
    "wgEncodeBroadHistoneGm12878ControlStdAlnRep1_chr22_sorted.bam_fragL200_bin200.txt")

binHM <- readBins( type=c("chip","input"), fileName=fileName)
fitHM <- mosaicsFit( binHM, analysisType="IO", bgEst="rMOM" )


###################################################
### code chunk number 21: mosaicsFitHMM
###################################################
hmmHM <- mosaicsFitHMM( fitHM, signalModel = "2S", 
	init="mosaics", init.FDR = 0.05, parallel=FALSE, nCore=8 )


###################################################
### code chunk number 22: mosaicshmmfit-show
###################################################
hmmHM


###################################################
### code chunk number 23: mosaicshmmfit-plot (eval = FALSE)
###################################################
## plot(hmmHM)


###################################################
### code chunk number 24: fig-mosaicshmmfit-plot
###################################################
plot(hmmHM)


###################################################
### code chunk number 25: mosaicsPeakHMM
###################################################
peakHM <- mosaicsPeakHMM( hmmHM, FDR = 0.05, decoding="posterior",
  thres=10, parallel=FALSE, nCore=8 )


###################################################
### code chunk number 26: summit-extractreads
###################################################
peakHM <- extractReads( peakHM,
  chipFile=system.file( file.path("extdata","wgEncodeBroadHistoneGm12878H3k4me3StdAlnRep1_chr22_sorted.bam"), package="mosaicsExample"),
  chipFileFormat="bam", chipPET=FALSE, chipFragLen=200,
  controlFile=system.file( file.path("extdata","wgEncodeBroadHistoneGm12878ControlStdAlnRep1_chr22_sorted.bam"), package="mosaicsExample"), 
  controlFileFormat="bam", controlPET=FALSE, controlFragLen=200, parallel=FALSE, nCore=8 )
peakHM


###################################################
### code chunk number 27: summit-findsummit
###################################################
peakHM <- findSummit( peakHM, parallel=FALSE, nCore=8 )


###################################################
### code chunk number 28: summit-findsummit-print
###################################################
head(print(peakHM))


###################################################
### code chunk number 29: postprocess-adjustboundary
###################################################
peakHM <- adjustBoundary( peakHM, parallel=FALSE, nCore=8 )
peakHM
head(print(peakHM))


###################################################
### code chunk number 30: postprocess-filterpeak
###################################################
peakHM <- filterPeak( peakHM, parallel=FALSE, nCore=8 )
peakHM
head(print(peakHM))


###################################################
### code chunk number 31: summit-findsummit-plot (eval = FALSE)
###################################################
## plot( peakHM, filename="./peakplot.pdf" )


###################################################
### code chunk number 32: summit-findsummit-plot-peaknum-print
###################################################
print(peakHM)[8,]
print(peakTFBS)[18,]


###################################################
### code chunk number 33: summit-findsummit-plot-peaknum (eval = FALSE)
###################################################
## plot( peakHM, peakNum=8 )


###################################################
### code chunk number 34: fig-summit-findsummit-plot-peaknum-hm
###################################################
plot( peakHM, peakNum=8 )


###################################################
### code chunk number 35: mosaicsRunAll (eval = FALSE)
###################################################
##      mosaicsRunAll( 
##          chipFile="/scratch/eland/STAT1_ChIP_eland_results.txt", 
##          chipFileFormat="eland_result", 
##          controlFile="/scratch/eland/STAT1_control_eland_results.txt", 
##          controlFileFormat="eland_result", 
##          binfileDir="/scratch/bin/", 
##          peakFile=c("/scratch/peak/STAT1_peak_list.bed", 
##                     "/scratch/peak/STAT1_peak_list.gff"),
##          peakFileFormat=c("bed", "gff"),
##          reportSummary=TRUE, 
##          summaryFile="/scratch/reports/mosaics_summary.txt", 
##          reportExploratory=TRUE, 
##          exploratoryFile="/scratch/reports/mosaics_exploratory.pdf", 
##          reportGOF=TRUE, 
##          gofFile="/scratch/reports/mosaics_GOF.pdf",
##          byChr=FALSE, useChrfile=FALSE, chrfile=NULL, excludeChr="chrM",
##          PET=FALSE, FDR=0.05, fragLen=200, binSize=fragLen, capping=0,
##          bgEst="rMOM", signalModel="BIC", parallel=TRUE, nCore=8 )


###################################################
### code chunk number 36: ts-readbin
###################################################
exampleBinData <- readBins( type=c("chip","input","M","GC","N"),
    fileName=c( system.file( file.path("extdata","chip_chr21.txt"), package="mosaicsExample"),
    system.file( file.path("extdata","input_chr21.txt"), package="mosaicsExample"),
    system.file( file.path("extdata","M_chr21.txt"), package="mosaicsExample"),
    system.file( file.path("extdata","GC_chr21.txt"), package="mosaicsExample"),
    system.file( file.path("extdata","N_chr21.txt"), package="mosaicsExample") ) )


###################################################
### code chunk number 37: ts-bindata-plot
###################################################
plot( exampleBinData, plotType="M" )
plot( exampleBinData, plotType="GC" )
plot( exampleBinData, plotType="M|input" )
plot( exampleBinData, plotType="GC|input" )  


###################################################
### code chunk number 38: fig-bindata-plot-M
###################################################
plot( exampleBinData, plotType="M" )


###################################################
### code chunk number 39: fig-bindata-plot-GC
###################################################
plot( exampleBinData, plotType="GC" )


###################################################
### code chunk number 40: fig-bindata-plot-M-input
###################################################
plot( exampleBinData, plotType="M|input" )


###################################################
### code chunk number 41: fig-bindata-plot-GC-input
###################################################
plot( exampleBinData, plotType="GC|input" )


###################################################
### code chunk number 42: ts-mosaicsfit (eval = FALSE)
###################################################
## exampleFit <- mosaicsFit( exampleBinData, analysisType="TS", bgEst="automatic" )


###################################################
### code chunk number 43: os-mosaicspeak (eval = FALSE)
###################################################
## OneSamplePeak <- mosaicsPeak( OneSampleFit, signalModel="2S", FDR=0.05,
## maxgap=200, minsize=50, thres=10 )


###################################################
### code chunk number 44: os-readbin (eval = FALSE)
###################################################
## exampleBinData <- readBins( type=c("chip","M","GC","N"),
##     fileName=c( system.file( file.path("extdata","chip_chr21.txt"), package="mosaicsExample"),
##     system.file( file.path("extdata","M_chr21.txt"), package="mosaicsExample"),
##     system.file( file.path("extdata","GC_chr21.txt"), package="mosaicsExample"),
##     system.file( file.path("extdata","N_chr21.txt"), package="mosaicsExample") ) )


###################################################
### code chunk number 45: os-mosaicsfit (eval = FALSE)
###################################################
## exampleFit <- mosaicsFit( exampleBinData, analysisType="OS", bgEst="automatic" )


###################################################
### code chunk number 46: os-mosaicspeak (eval = FALSE)
###################################################
## examplePeak <- mosaicsPeak( exampleFit, signalModel="2S", FDR=0.05,
## maxgap=200, minsize=50, thres=10 )


###################################################
### code chunk number 47: tuning-case1 (eval = FALSE)
###################################################
## exampleFit <- mosaicsFit( exampleBinData, analysisType="IO", bgEst="rMOM" )


###################################################
### code chunk number 48: tuning-case2 (eval = FALSE)
###################################################
## exampleFit <- mosaicsFit( exampleBinData, analysisType="IO", 
## 	bgEst="automatic", truncProb=0.80 )