### R code from vignette source 'vignettes/PWMEnrich/inst/doc/PWMEnrich.Rnw'

###################################################
### code chunk number 1: simple
###################################################
library(PWMEnrich)
library(PWMEnrich.Dmelanogaster.background)

# load the pre-compiled lognormal background
data(PWMLogn.dm3.jaspar.insects)

# scan two sequences from a FASTA file for motif enrichment
sequences = read.DNAStringSet(system.file(package="PWMEnrich", 
  dir="extdata", file="example.fa"))
sequences  
res = motifEnrichment(sequences, PWMLogn.dm3.jaspar.insects)


###################################################
### code chunk number 2: simple-affinity
###################################################
# PWMs enriched in first sequence (Lognormal P-values)
head(sort(res$sequence.bg[1,]))
# PWMs enriched in second sequence (Lognormal P-values)
head(sort(res$sequence.bg[2,]))

# most enriched in the group of two sequences (Lognormal P-value)
head(sort(res$group.bg))


###################################################
### code chunk number 3: plots
###################################################
sorted.motif.names = names(sort(res$group.bg))
pwms = PWMLogn.dm3.jaspar.insects$pwms
plotMultipleMotifs(pwms[sorted.motif.names], sorted.motif.names, 3, 2, 
	xmargin.scale=0.8, ymargin.scale=0.8, xfontsize=11, yfontsize=11,
	titlefontsize=14)


###################################################
### code chunk number 4: counts
###################################################
data(PWMPvalueCutoff1e3.dm3.jaspar.insects)

res.count = motifEnrichment(sequences, PWMPvalueCutoff1e3.dm3.jaspar.insects)

# First sequence, PWMs sorted by number of motif hits with P-value < 0.001
head(sort(res.count$sequence.nobg[1,], decreasing=TRUE))

# Second sequence
head(sort(res.count$sequence.nobg[2,], decreasing=TRUE))

# the whole group, PWMs sorted by Z-score for the number of motif hits
head(sort(res.count$group.bg, decreasing=TRUE), 8)


###################################################
### code chunk number 5: diff
###################################################
res.diff = motifDiffEnrichment(sequences[1], sequences[2], 
    PWMLogn.dm3.jaspar.insects)

# motifs differentially enriched in the first sequence (large positive values)
head(sort(res.diff$group.bg, decreasing=TRUE))

# motifs differentially enriched in the second sequence (large negative values)
head(sort(res.diff$group.bg))


###################################################
### code chunk number 6: parallel
###################################################
registerCoresPWMEnrich(4)


###################################################
### code chunk number 7: parallel-stop
###################################################
registerCoresPWMEnrich(NULL)


###################################################
### code chunk number 8: read
###################################################
motifs.denovo = readMotifs(system.file(package="PWMEnrich", 
  dir="extdata", file="example.transfac"), remove.acc=TRUE)

motifs.denovo

bg.denovo = makeBackground(motifs.denovo, organism="dm3", type="logn", quick=TRUE)

res.denovo = motifEnrichment(sequences, bg.denovo)
head(sort(res.denovo$group.bg))


###################################################
### code chunk number 9: bg-investigate
###################################################
bg.denovo
bg.denovo$bg.mean


###################################################
### code chunk number 10: pfmtopwm
###################################################
library("BSgenome.Dmelanogaster.UCSC.dm3")
# make a lognormal background for the two motifs using only first 20 promoters
bg.seq = Dmelanogaster$upstream2000[1:20]
# the sequences are split into 100bp chunks and fitted
bg.custom = makePWMLognBackground(bg.seq, motifs.denovo, bg.len=100, 
    bg.source="20 promoters split into 100bp chunks")
bg.custom


###################################################
### code chunk number 11: sessionInfo
###################################################
toLatex(sessionInfo())