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

###################################################
### code chunk number 1: ImportLibrary
###################################################
library(DiffLogo)


###################################################
### code chunk number 2: ImportMotifsFromMotifDb
###################################################
library(MotifDb)

## import motifs
hitIndeces <- grep ('CTCF', 
                    values (MotifDb)$geneSymbol, 
                    ignore.case=TRUE)
list <- as.list(MotifDb[hitIndeces])

## get two motifs
pwm1 <- list$"Hsapiens-jolma2013-CTCF"
## trim and reverse complement
pwm2 <- list$"Hsapiens-JASPAR_CORE-CTCF-MA0139.1"[4:1, 18:2]


###################################################
### code chunk number 3: ImportMotifsFromDiffLogoPackage
###################################################
## import nine DNA motifs for transcription factor CTCF from matrix
motif_folder <- "extdata/pwm"
motif_names_dna = c(
  "GM12878", "H1-hESC", "HeLa-S3", "HepG2",   "HUVEC", 
  "K562",    "MCF7",    "NHEK",    "ProgFib")
motifs_dna = list()
for (name in motif_names_dna) {
  fileName <- paste(motif_folder,"/",name,".txt",sep="")
  file <- system.file(fileName, package = "DiffLogo")
  motifs_dna[[name]] <- as.matrix(read.delim(file, FALSE))
}
## import DNA motifs for three transcription factors from table
motif_folder <- "extdata/alignments"
motif_names_dna2 <- c("Mad", "Max", "Myc")
motifs_dna2 <- list()
for (name in motif_names_dna2) {
  fileName <- paste(motif_folder,"/",name,".txt",sep="")
  file <- system.file(fileName, package = "DiffLogo")
  fileContent <- readLines(file)
  fileContent <- unlist(lapply(
    X = fileContent, 
    FUN = function(x){ strsplit(x = x, split = "\t")[[1]][[1]] }))
  motifs_dna2[[name]] <- getPwmFromAlignment(fileContent, DNA, 1)
}
## import three ASN motifs for one protein domain from fasta files
motif_folder = "extdata/alignments"
motif_names_asn = c("F-box_fungi.seq", "F-box_metazoa.seq", 
                    "F-box_viridiplantae.seq")
motifs_asn = list()
for (name in motif_names_asn) {
  fileName = paste(motif_folder,"/",name,".fa",sep="")
  file = system.file(fileName, package = "DiffLogo")
  fileContent <- readLines(file)
  fileContent <- fileContent[seq(from = 2, by = 2, 
                   length.out = floor(length(fileContent)/2))]
  motifs_asn[[name]] <- getPwmFromAlignment(fileContent, ASN, 1)
}


###################################################
### code chunk number 4: PlotSequenceLogo
###################################################
## plot classic sequence logo
library(seqLogo)
seqLogo::seqLogo(pwm = pwm1)


###################################################
### code chunk number 5: PlotCustomSequenceLogo
###################################################
## plot custom sequence logo
par(mfrow=c(2,1), pin=c(3, 1), mar = c(2, 4, 1, 1))
DiffLogo::seqLogo(pwm = pwm1)
DiffLogo::seqLogo(pwm = pwm2, stackHeight = sumProbabilities)
par(mfrow=c(1,1), pin=c(1, 1), mar=c(5.1, 4.1, 4.1, 2.1))


###################################################
### code chunk number 6: PlotDiffLogo
###################################################
## plot DiffLogo
diffLogoFromPwm(pwm1 = pwm1, pwm2 = pwm2)

## diffLogoFromPwm is a convenience function for
diffLogoObj = createDiffLogoObject(pwm1 = pwm1, pwm2 = pwm2)
diffLogo(diffLogoObj)


###################################################
### code chunk number 7: PlotDiffLogoTable
###################################################
## plot table of difference logos for CTFC motifs (DNA)
diffLogoTable(PWMs = motifs_dna, )
## plot table of difference logos for E-Box motifs (DNA)
diffLogoTable(PWMs = motifs_dna2)
## plot table of difference logos for F-Box motifs (ASN)
diffLogoTable(PWMs = motifs_asn, alphabet = ASN)


###################################################
### code chunk number 8: ExportToFigure
###################################################
## parameters
widthToHeightRatio = 16/10;
size = length(motifs_dna) * 2
resolution <- 300
width <- size * widthToHeightRatio
height <- size

## export single DiffLogo as pdf document
fileName <- "Comparison_of_two_motifs.pdf"
pdf(file = fileName, width = width, height = height)
diffLogoFromPwm(pwm1 = pwm1, pwm2 = pwm2)
dev.off()

## export DiffLogo table as png image
fileName <- "Comparison_of_multiple_motifs.png"
png(
  filename = fileName, res = resolution,
  width = width * resolution, height = height * resolution)
diffLogoTable(PWMs = motifs_dna)
dev.off()