## ---- eval=TRUE, warnings = FALSE, echo=TRUE,message=FALSE-------------------- library(ChIPanalyser) #Load data data(ChIPanalyserData) # Loading DNASequenceSet from BSgenome object # We recommend using the latest version of the genome # Please ensure that all your data is aligned to the same version of the genome library(BSgenome.Dmelanogaster.UCSC.dm3) DNASequenceSet <-getSeq(BSgenome.Dmelanogaster.UCSC.dm3) #Loading Position Frequency Matrix PFM <- file.path(system.file("extdata",package="ChIPanalyser"),"BCDSlx.pfm") #Checking if correctly loaded ls() ## ----eval=TRUE, warnings = FALSE---------------------------------------------- eveLocusChip<-processingChIP(profile=eveLocusChip, loci=eveLocus, cores=1) eveLocusChip ## ---- eval =TRUE-------------------------------------------------------------- # PFMs are automatically converted to PWM when build genomicProfiles GP<-genomicProfiles(PFM=PFM,PFMFormat="raw") GP ## ---- eval=FALSE-------------------------------------------------------------- # GP<-genomicProfiles(PWM=PositionWeightMatrix) ## ----eval=TRUE,warnings=FALSE------------------------------------------------- ## surpress dependency warnings optimal<-suppressWarnings(computeOptimal(genomicProfiles=GP, DNASequenceSet=DNASequenceSet, ChIPScore=eveLocusChip, chromatinState=Access)) ## ---- eval=TRUE, echo=TRUE---------------------------------------------------- ## Lambda Values seq(0.25,5,by=0.25) ## Bound Molecule Values c(1, 10, 20, 50, 100, 200, 500,1000,2000, 5000,10000,20000,50000, 100000, 200000, 500000, 1000000) ## ---- eval =T----------------------------------------------------------------- optimalParam<-optimal$Optimal optimalParam$OptimalParameters ## ---- eval=TRUE, warnings = FALSE, fig.width=10, fig.height=8----------------- # Plotting Optimal heat maps par(oma=c(0,0,3,0)) layout(matrix(1:8,ncol=4, byrow=T),width=c(6,1.5,6,1.5),height=c(1,1)) plotOptimalHeatMaps(optimalParam,layout=FALSE) ## ----eval=TRUE---------------------------------------------------------------- optimalParam<-searchSites(optimal,lambdaPWM=1.25,BoundMolecules=10000) ## ---- eval=TRUE,fig.width=12, fig.height=4.5---------------------------------- plotOccupancyProfile(predictedProfile=optimalParam$ChIPProfiles, ChIPScore=eveLocusChip, chromatinState=Access, occupancy=optimalParam$Occupancy, goodnessOfFit=optimalParam$goodnessOfFit) ## ---- eval =TRUE-------------------------------------------------------------- ## Suggested Parameters to start with. parameterOptions() ## Changing parameters PO <- parameterOptions(noiseFilter="sigmoid",chipSd=150,chipMean=150,lociWidth=30000) ## ----eval=FALSE--------------------------------------------------------------- # ## Top 50 loci based on ChIP score # processingChIP(profile="/path/to/ChIP", # loci=NULL, # reduce=50, # parameterOptions=PO) # # ## Top 50 loci ALSO containing peaks # processingChIP(profile="/path/to/ChIP", # loci=NULL, # reduce=50, # peaks="/path/to/peaks", # parameterOptions=PO) # # ## Top 50 loci containing BOTH peaks and Accessible DNA # processingChIP(profile="/path/to/ChIP", # loci=NULL, # reduce=50, # peaks="/path/to/peaks", # chromatinState="/path/to/chromatinState" # parameterOptions=PO) # ## ---- eval=TRUE--------------------------------------------------------------- str(genomicProfiles()) GP <- genomicProfiles(PFM=PFM, PFMFormat="raw", BPFrequency=DNASequenceSet) GP ## ---- eval=FALSE-------------------------------------------------------------- # ## Parsing pre computed parameters (processingChIP function) # GP<-genomicProfiles(PFM=PFM, PFMFormat="raw", BPFrequency=DNASequenceSet, # ChIPScore=ChIPScore) # # ## Parsing pre assigned function (parameterOptions) # parameterOptions<-parameterOptions(lambdaPWM=c(1,2,3), # boundMolecules=c(5,50,500)) # GP<-genomicProfiles(PFM=PFM, PFMFormat="raw", BPFrequency=DNASequenceSet, # parameterOptions=parameterOptions) # # ## Direct parameter assignement # # GP<-genomicProfiles(PFM=PFM, PFMFormat="raw", BPFrequency=DNASequenceSet, # lambdaPWM=c(1,2,3),boundMolecules=c(4,500,8000)) ## ---- eval=FALSE-------------------------------------------------------------- # ## Setting custom parameters # OP<-parameterOptions(lambdaPWM=seq(1,10,by=0.5), # boundMolecules=seq(1,100000, length.out=20)) # # ## Computing ONLY Optimal Parameters and MSE as goodness Of Fit metric # optimal<-computeOptimal(genomicProfiles=GP, # DNASequenceSet=DNASequenceSet, # ChIPScore=eveLocusChip, # chromatinState=Access, # parameterOptions=OP, # optimalMethod="MSE", # returnAll=FALSE) # # ### Computing ONLY Optimal Parameters and using Rank slection method # optimal<-computeOptimal(genomicProfiles=GP, # DNASequenceSet=DNASequenceSet, # ChIPScore=eveLocusChip, # chromatinState=Access, # parameterOptions=OP, # optimalMethod="all", # rank=TRUE) # ## ---- eval=FALSE-------------------------------------------------------------- # ## Extracted Optimal Parameters # optimalParam<-optimal$Optimal # # ## Plotting heat maps # plotOptimalHeatMaps(optimalParam,overlay=TRUE) ## ----eval=TRUE,warnings=FALSE------------------------------------------------- ## Creating genomic Profiles object with PFMs and associated parameters GP <- genomicProfiles(PFM=PFM,PFMFormat="raw",BPFrequency=DNASequenceSet, lambdaPWM=1, boundMolecules=58794) ## Computing Genome Wide Score required GW <- computeGenomeWideScores(genomicProfiles=GP, DNASequenceSet=DNASequenceSet, chromatinState=Access) GW ## Computing PWM score above threshold pwm <- computePWMScore(genomicProfiles=GW, DNASequenceSet=DNASequenceSet, loci=eveLocusChip, chromatinState=Access) pwm ## Computing Occupancy of sites above threshold occup <- computeOccupancy(genomicProfiles=pwm) occup ## Compute ChIP seq like profiles chip <- computeChIPProfile(genomicProfiles=occup, loci=eveLocusChip) chip ## Compute goodness Of Fit of model accu <- profileAccuracyEstimate(genomicProfiles=chip, ChIPScore=eveLocusChip) accu ## ---- eval=TRUE,fig.width=12, fig.height=4.5---------------------------------- plotOccupancyProfile(predictedProfile=chip, ChIPScore=eveLocusChip, chromatinState=Access, occupancy=occup, goodnessOfFit=accu, geneRef=geneRef) ## ---- eval=TRUE,echo=TRUE----------------------------------------------------- parameterOptions() ## ---- eval=T, echo=T---------------------------------------------------------- str(genomicProfiles()) ## ---- eval=F, echo=T---------------------------------------------------------- # ## Accessors and Setters for parameterOptions and genomicProfiles # avrageExpPWMScore(obj) # backgroundSignal(obj) # backgroundSignal(obj)<-value # boundMolecules(obj) # boundMolecules(obj)<-value # BPFrequency(obj) # BPFrequency(obj)<-value # chipMean(obj) # chipMean(obj)<-value # chipSd(obj) # chipSd(obj)<-value # chipSmooth(obj) # chipSmooth(obj)<-value # DNASequenceLength(obj) # drop(obj) # lambdaPWM(obj) # lambdaPWM(obj)<-value # lociWidth(obj) # lociWidth(obj)<-value # maxPWMScore(obj) # maxSignal(obj) # maxSignal(obj)<-value # minPWMScore(obj) # naturalLog(obj) # naturalLog(obj)<-value # noiseFilter(obj) # noiseFilter(obj)<-value # noOfSites(obj) # noOfSites(obj)<-value # PFMFormat(obj) # PFMFormat(obj)<-value # ploidy(obj) # ploidy(obj)<-value # PositionFrequencyMatrix(obj) # PositionFrequencyMatrix(obj)<-value # PositionWeightMatrix(obj) # PositionWeightMatrix<-value # profiles(obj) # PWMpseudocount(obj) # PWMpseudocount(obj)<-value # PWMThreshold(obj) # PWMThreshold(obj)<-value # removeBackground(obj) # removeBackground(obj)<-value # # stepSize(obj) # stepSize(obj)<-value # strandRule(obj) # strandRule(obj)<-value # whichstrand(obj) # whichstrand(obj)<-value # # ## ChIPScore slots accessors # loci(obj) # scores(obj) ## ----eval=T------------------------------------------------------------------- sessionInfo()