### 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.MotifDb.Dmel)

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


###################################################
### code chunk number 2: simple-affinity
###################################################
res

# PWMs enriched in both sequences (Lognormal P-value)
head(motifRankingForGroup(res))

# PWMs enriched in the first sequence (Lognormal P-values)
head(motifRankingForSequence(res, 1))
# PWMs enriched in the second sequence (Lognormal P-values)
head(motifRankingForSequence(res, 2))


###################################################
### code chunk number 3: simple-affinity
###################################################
head(motifRankingForGroup(res, id=TRUE))

library(MotifDb)
db = values(MotifDb)

db[db$providerName == "dimm_da_SANGER_5_FBgn0000413", "dataSource"]
db[db$providerName == "dimm_da_SANGER_5_FBgn0023091", "dataSource"]


###################################################
### code chunk number 4: plots
###################################################
plotTopMotifsGroup(res, 6, rows=2, cols=3, xmargin.scale=0.7)


###################################################
### code chunk number 5: counts
###################################################
data(PWMPvalueCutoff1e4.dm3.MotifDb.Dmel)

res.count = motifEnrichment(sequences, PWMPvalueCutoff1e4.dm3.MotifDb.Dmel)

#  PWMs sorted by number of motifs hits (P-value < 0.0001) in both sequences
head(motifRankingForGroup(res.count, bg=FALSE))

#  PWMs sorted by the Z-score of motif enrichment
head(motifRankingForGroup(res.count))

# First sequence, PWMs sorted by number of motif hits (P-value < 0.0001)
head(motifRankingForSequence(res.count, 1, bg=FALSE))

# Second sequence
head(motifRankingForSequence(res.count, 2, bg=FALSE))


###################################################
### 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(motifRankingForGroup(res.denovo))


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