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

# load the stripe2 sequences from a FASTA file for motif enrichment
sequence = readDNAStringSet(system.file(package="PWMEnrich", 
  dir="extdata", file="stripe2.fa"))
sequence  

# perform motif enrichment!
res = motifEnrichment(sequence, PWMLogn.dm3.MotifDb.Dmel)

report = sequenceReport(res, 1)
report

# plot the top 30 most enriched motifs
plot(report[1:30], fontsize=7, id.fontsize=5)


###################################################
### code chunk number 2: stripe2visual
###################################################

# extract the 3 PWMs for the TFs we are interested in
ids = c("Dmelanogaster-FlyFactorSurvey-bcd_FlyReg_FBgn0000166", 
        "Dmelanogaster-FlyFactorSurvey-gt_FlyReg_FBgn0001150", 
        "Dmelanogaster-JASPAR_CORE-Kr-MA0452.1")
sel.pwms = PWMLogn.dm3.MotifDb.Dmel$pwms[ids]
names(sel.pwms) = c("bcd", "gt", "Kr")

# scan and get the raw scores
scores = motifScores(sequence, sel.pwms, raw.scores=TRUE)

# raw scores for the first (and only) input sequence
dim(scores[[1]])
head(scores[[1]])

# score starting at position 1 of forward strand
scores[[1]][1, "bcd"]
# score for the reverse complement of the motif, starting at the same position
scores[[1]][485, "bcd"]

# plot
plotMotifScores(scores, cols=c("green", "red", "blue"))


###################################################
### code chunk number 3: motifEcdf
###################################################
library(BSgenome.Dmelanogaster.UCSC.dm3)

# empirical distribution for the bcd motif
bcd.ecdf = motifEcdf(sel.pwms$bcd, Dmelanogaster, quick=TRUE)[[1]]

# find the score that is equivalent to the P-value of 1e-3
threshold.1e3 = log2(quantile(bcd.ecdf, 1 - 1e-3))
threshold.1e3 

# replot only the bcd motif hits with the P-value cutoff of 1e-3 (0.001)
plotMotifScores(scores, cols="green", sel.motifs="bcd", cutoff=threshold.1e3)

# P-value at each position
pvals = 1 - bcd.ecdf(scores[[1]][,"bcd"])

# position where the P-value is smaller that 1e-3
which(pvals < 1e-3)


###################################################
### code chunk number 4: tinman
###################################################
library(PWMEnrich.Dmelanogaster.background)

# load the pre-compiled lognormal background
data(PWMLogn.dm3.MotifDb.Dmel)

sequences = readDNAStringSet(system.file(package="PWMEnrich", 
  dir="extdata", file="tinman-early-top20.fa"))
 
res = motifEnrichment(sequences, PWMLogn.dm3.MotifDb.Dmel)

report = groupReport(res)
report

plot(report[1:10], fontsize=7, id.fontsize=5)


###################################################
### code chunk number 5: tinman-alt
###################################################
report.top = groupReport(res, by.top.motifs=TRUE)

report.top


###################################################
### code chunk number 6: tinman2
###################################################
res

# raw scores
res$sequence.nobg[1:5, 1:2]

# P-values
res$sequence.bg[1:5, 1:2]


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


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


###################################################
### code chunk number 9: bigmem
###################################################
useBigMemoryPWMEnrich(TRUE)


###################################################
### code chunk number 10: bigmemoff
###################################################
useBigMemoryPWMEnrich(FALSE)


###################################################
### code chunk number 11: readmotifs
###################################################
library(PWMEnrich.Dmelanogaster.background)

motifs.denovo = readMotifs(system.file(package="PWMEnrich", 
  dir="extdata", file="example.transfac"), remove.acc=TRUE)
  
motifs.denovo  

# convert count matrices into PWMs
genomic.acgt = getBackgroundFrequencies("dm3")
pwms.denovo = toPWM(motifs.denovo, prior=genomic.acgt)

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

# use new motifs for motif enrichment
res.denovo = motifEnrichment(sequences[1:5], bg.denovo)
groupReport(res.denovo)


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


###################################################
### code chunk number 13: custombg
###################################################
library(PWMEnrich.Dmelanogaster.background)
data(dm3.upstream2000)

# make a lognormal background for the two motifs using only first 20 promoters
bg.seq = dm3.upstream2000[1:20]

# the sequences are split into 100bp chunks and fitted
bg.custom = makeBackground(pwms.denovo, bg.seq=bg.seq, type="logn", bg.len=100,
	bg.source="20 promoters split into 100bp chunks")

bg.custom


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