cobindR.R
In cobindR: Finding Co-occuring motifs of transcription factor binding sites

### R code from vignette source 'cobindR.Rnw'

###################################################
### code chunk number 1: load.libraries
###################################################
library(cobindR)
library(Biostrings)


###################################################
### code chunk number 2: make.configuration.from.file
###################################################
#run cobindR
cfg <- cobindRConfiguration( fname = 
             system.file('extdata/config_default.yml',
             package='cobindR'))


###################################################
### code chunk number 3: make.configuration.from.script (eval = FALSE)
###################################################
## #run cobindR
## cfg <- cobindRConfiguration()


###################################################
### code chunk number 4: set.PFM.path
###################################################
pfm_path(cfg) <- 
  system.file('extdata/pfms',package='cobindR')


###################################################
### code chunk number 5: define.pairs
###################################################
pairs(cfg) <- c('ES_Sox2_1_c1058 ES_Oct4_1_c570')


###################################################
### code chunk number 6: set.PFM.path (eval = FALSE)
###################################################
## pfm_path(cfg) <- 'MotifDb'
## pairs(cfg) <- c('JASPAR_CORE:CREB1 JASPAR_CORE:KLF4',
## 		'JASPAR_CORE:CREB1 JASPAR_CORE:KLF4')


###################################################
### code chunk number 7: configure.sequence.input
###################################################
sequence_type(cfg) <- 'fasta'
sequence_source(cfg) <- system.file('extdata/sox_oct_example_vignette_seqs.fasta',
                                            package='cobindR')
sequence_origin(cfg) <- 'Mouse Embryonic Stem Cell Example ChIP-Seq Oct4 Peak Sequences'
species(cfg) <- 'Mus musculus'


###################################################
### code chunk number 8: configure.geneid.sequence.input (eval = FALSE)
###################################################
## tmp.geneid.file <- tempfile(pattern = "cobindR_sample_seq", 
##                      tmpdir = tempdir(), fileext = ".txt")
## write(c('#cobindR Example Mouse Genes','ENSMUSG00000038518',
## 		'ENSMUSG00000042596','ENSMUSG00000025927'),
## 		file = tmp.geneid.file)
## species(cfg) <- 'Mus musculus'
## sequence_type(cfg) <- 'geneid'
## sequence_source(cfg) <- tmp.geneid.file
## sequence_origin(cfg) <- 'ENSEMBL genes'
## upstream(cfg) <- downstream(cfg) <- 500


###################################################
### code chunk number 9: configure.chipseq.sequence.input (eval = FALSE)
###################################################
## sequence_type(cfg) <- 'chipseq'
## sequence_source(cfg) <- 
##   system.file('extdata/ucsc_example_ItemRGBDemo.bed',
## 				package='cobindR')
## sequence_origin(cfg) <- 'UCSC bedfile example'
## species(cfg) <- 'BSgenome.Mmusculus.UCSC.mm9'


###################################################
### code chunk number 10: define.background.model
###################################################
bg_sequence_type(cfg) <-  'markov.3.200'# or 'ushuffle.3.3000' or 'local.10.1000'


###################################################
### code chunk number 11: perform.search
###################################################
cobindR.bs <- cobindr( cfg, 
              name='cobind test using sampled sequences')
cobindR.bs <- search.pwm(cobindR.bs, min.score = '80')
cobindR.bs <- find.pairs(cobindR.bs, n.cpu = 3)


###################################################
### code chunk number 12: perform.rtfbssearch (eval = FALSE)
###################################################
## cobindR.bs = rtfbs(cobindR.bs)


###################################################
### code chunk number 13: perform.gadem (eval = FALSE)
###################################################
## cobindR.bs = search.gadem(cobindR.bs, deNovo=TRUE)


###################################################
### code chunk number 14: make.seq.background
###################################################
cobindR.bs <- generate.background(cobindR.bs)


###################################################
### code chunk number 15: perform.background.search
###################################################
cobindR.bs <- search.pwm(cobindR.bs, min.score='80', 
                        background_scan = TRUE)
cobindR.bs <- find.pairs(cobindR.bs, background_scan = TRUE, 
                        n.cpu = 3)


###################################################
### code chunk number 16: plot_gc_test (eval = FALSE)
###################################################
## tmp <- testCpG(cobindR.bs, do.plot=T)


###################################################
### code chunk number 17: plot_gc_profile (eval = FALSE)
###################################################
## plot.gc(cobindR.bs, wind.size=200, frac = 2) 


###################################################
### code chunk number 18: plot_predicted_logos (eval = FALSE)
###################################################
## pred.motifs <- predicted2pwm(cobindR.bs, as.pfm=TRUE)
## # normalized column sums as required by seqLogo
## pred.norm.motifs <- lapply(pred.motifs, function(x) x / colSums(x)[1]) 
## 
## # load sequence logo plot function
## plot.tfbslogo(x=cobindR.bs,c('ES_Sox2_1_c1058','ES_Oct4_1_c570'))


###################################################
### code chunk number 19: plot_positions_profile (eval = FALSE)
###################################################
## plot.positionprofile(cobindR.bs)


###################################################
### code chunk number 20: plot_positions_simple (eval = FALSE)
###################################################
## plot.positions.simple(cobindR.bs)


###################################################
### code chunk number 21: plot_heatmap (eval = FALSE)
###################################################
## plot.tfbs.heatmap(cobindR.bs, include.empty.seqs=FALSE)


###################################################
### code chunk number 22: plot_venn (eval = FALSE)
###################################################
## plot.tfbs.venndiagram(cobindR.bs, pwms = c('ES_Sox2_1_c1058', 'ES_Oct4_1_c570'), include.empty.seqs=FALSE)


###################################################
### code chunk number 23: plot_pairdistance (eval = FALSE)
###################################################
## plot.pairdistance(cobindR.bs, pwm1='ES_Sox2_1_c1058',
##                   pwm2='ES_Oct4_1_c570')


###################################################
### code chunk number 24: plot_pairdistribution (eval = FALSE)
###################################################
## plot.pairdistribution(cobindR.bs, pwm1='ES_Sox2_1_c1058',
##                       pwm2='ES_Oct4_1_c570')


###################################################
### code chunk number 25: plot_detrending (eval = FALSE)
###################################################
## plot.detrending(cobindR.bs, pwm1='ES_Sox2_1_c1058',
## 				pwm2='ES_Oct4_1_c570', bin_length=10, abs.distance=FALSE, overlap=4)


###################################################
### code chunk number 26: write.detrending
###################################################
tmp.sig.pairs = get.significant.pairs(x = cobindR.bs, pwm1='ES_Sox2_1_c1058',pwm2='ES_Oct4_1_c570', bin_length=10, abs.distance=FALSE,overlap=4)

tmp.resultbs.file <- tempfile(pattern = "cobindR_detrening_result_bindingsites", tmpdir = tempdir(), fileext = ".tsv")
write.table(tmp.sig.pairs[[1]], file=tmp.resultbs.file, sep="\t", quote=F)
system(paste('head',tmp.resultbs.file))

tmp.resultcp.file = gsub("bindingsites","candidates_pairs", tmp.resultbs.file)
write.table(tmp.sig.pairs[[2]], file=tmp.resultcp.file, sep="\t", quote=F)
system(paste('head',tmp.resultcp.file))


###################################################
### code chunk number 27: write.bindingsites
###################################################
tmp.result.bs.file <- tempfile(pattern = "cobindR_bindingsite_pred", 
                     tmpdir = tempdir(), fileext = ".txt")
write.bindingsites(cobindR.bs, file=tmp.result.bs.file, background=FALSE)
system(paste('head',tmp.result.bs.file))


###################################################
### code chunk number 28: write.bindingsites
###################################################
tmp.inseq.file <- tempfile(pattern = "cobindR_input_sequences", 
                     tmpdir = tempdir(), fileext = ".fasta")
# slotname = 'bg_sequences' to obtain the background sequences
write.sequences(cobindR.bs, file=tmp.inseq.file,
				slotname= "sequences")
#system(paste('head',tmp.inseq.file,'n=10'))


###################################################
### code chunk number 29: cleanup
###################################################
try(unlink(tmp.result.file))
try(unlink(tmp.result.bs.file))
try(unlink(tmp.inseq.file))