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

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex(use.unsrturl=FALSE)


###################################################
### code chunk number 2: rtl-PFMatrix
###################################################
library(TFBSTools)
pfm = PFMatrix(ID="MA0004.1", name="Arnt", matrixClass="Zipper-Type", strand="+",
        bg=c(A=0.25, C=0.25, G=0.25, T=0.25),
        tags=list(family="Helix-Loop-Helix", species="10090", 
                  tax_group="vertebrates",medline="7592839", type="SELEX", 
                  ACC="P53762", pazar_tf_id="TF0000003",
                  TFBSshape_ID="11", TFencyclopedia_ID="580"),
        profileMatrix=matrix(c(4L,  19L, 0L,  0L,  0L,  0L,
                               16L, 0L,  20L, 0L,  0L,  0L,
                               0L,  1L,  0L,  20L, 0L,  20L,
                               0L,  0L,  0L,  0L,  20L, 0L),
                             byrow=TRUE, nrow=4, 
                             dimnames=list(c("A", "C", "G", "T")))
        )

## coerced to matrix
as.matrix(pfm)


## access the slots of pfm
ID(pfm)
name(pfm)
Matrix(pfm)


## conversion to pwm, icm
pwm = toPWM(pfm)
pwm 

icm = toICM(pfm)
icm

## get the reverse complment matrix with all the same information except the strand.
reverseComplement(pwm)


###################################################
### code chunk number 3: rtl-PFMatrixList
###################################################
pfm2 = pfm
pfmList = PFMatrixList(pfm1=pfm, pfm2=pfm2, use.names=TRUE)
pfmList
names(pfmList)


###################################################
### code chunk number 4: rtl-searchDB
###################################################
library(JASPAR2014)
opts = list()
opts[["species"]] = 9606
opts[["name"]] = "RUNX1"
#opts[["class"]] = "Ig-fold"
opts[["type"]] = "SELEX"
opts[["all_versions"]] = TRUE
PFMatrixList = getMatrixSet(JASPAR2014, opts)
PFMatrixList

opts2 = list()
opts2[["type"]] = "SELEX"
PFMatrixList2 = getMatrixSet(JASPAR2014, opts2)
PFMatrixList2


###################################################
### code chunk number 5: rtl-operateDb (eval = FALSE)
###################################################
## db = "myMatrixDb.sqlite"
## initializeJASPARDB(db)
## storeMatrix(db, pfm)
## deleteMatrixHavingID(db, "MA0003")


###################################################
### code chunk number 6: rtl-matrixMethods
###################################################
pwm = toPWM(pfm, pseudocounts=0.8)
pwm


###################################################
### code chunk number 7: rtl-matrixMethods
###################################################
icm = toICM(pfm, pseudocounts=0.8, schneider=TRUE)
icm


###################################################
### code chunk number 8: seqLogo1
###################################################
seqLogo(icm)


###################################################
### code chunk number 9: rtl-PFMSimi
###################################################
## one to one comparison
data(MA0003.2)
data(MA0004.1)
pfmSubject = MA0003.2
pfmQuery =  MA0004.1
PFMSimilarity(pfmSubject, pfmQuery)

## one to several comparsion
PFMSimilarity(pfmList, pfmQuery)

## align IUPAC string
IUPACString = "ACGTMRWSYKVHDBN"
PFMSimilarity(pfmList, IUPACString)


###################################################
### code chunk number 10: permuteMatrix
###################################################
#library(JASPAR2014)
#pfmSubject = getMatrixByID(JASPAR2014, ID="MA0043")
#pfmQuery = getMatrixByID(JASPAR2014, ID="MA0048")
#opts = list()
#opts[["class"]] = "Ig-fold"
#pfmList = getMatrixSet(JASPAR2014, opts)
permuteMatrix(pfmQuery)
permuteMatrix(pfmList, type="intra")
permuteMatrix(pfmList, type="inter")


###################################################
### code chunk number 11: rtl-PWMSimilarity
###################################################
  data(MA0003.2)
  data(MA0004.1)
  pwm1 = toPWM(MA0003.2, type="prob")
  pwm2 = toPWM(MA0004.1, type="prob")
  PWMSimilarity(pwm1, pwm2, method="Euclidean")
  PWMSimilarity(pwm1, pwm2, method="Pearson")
  PWMSimilarity(pwm1, pwm2, method="KL")


###################################################
### code chunk number 12: rtl-search
###################################################
library(Biostrings)
data(MA0003.2)
data(MA0004.1)
pwmList = PWMatrixList(MA0003.2=toPWM(MA0003.2), MA0004.1=toPWM(MA0004.1), 
                       use.names=TRUE)
subject = DNAString("GAATTCTCTCTTGTTGTAGTCTCTTGACAAAATG")
siteset = searchSeq(pwm, subject, seqname="seq1", min.score="60%", strand="*")

sitesetList = searchSeq(pwmList, subject, seqname="seq1",
                        min.score="60%", strand="*")


## generate gff style output
head(writeGFF3(siteset))
head(writeGFF3(sitesetList))
head(writeGFF2(siteset))

## get the relative score
relScore(siteset)
relScore(sitesetList)


###################################################
### code chunk number 13: rtl-searchAln
###################################################
library(Biostrings)
data(MA0003.2)
pwm = toPWM(MA0003.2)
aln1 = DNAString("ACTTCACCAGCTCCCTGGCGGTAAGTTGATC---AAAGG---AAACGCAAAGTTTTCAAG")
aln2 = DNAString("GTTTCACTACTTCCTTTCGGGTAAGTAAATATATAAATATATAAAAATATAATTTTCATC")
sitePairSet = searchAln(pwm, aln1, aln2, seqname1="seq1", seqname2="seq2",
                        min.score="50%", cutoff=0.5, 
                        strand="*", type="any")
## generate gff style output
head(writeGFF3(sitePairSet))
head(writeGFF2(sitePairSet))

## search the Axt alignment
library(CNEr)
axtFilesHg19DanRer7 <- file.path(system.file("extdata", package="CNEr"),
                                 "hg19.danRer7.net.axt")
axtHg19DanRer7 <- readAxt(axtFilesHg19DanRer7)
sitePairSet <-  searchAln(pwm, axtHg19DanRer7, min.score="80%",
                          windowSize=51L, cutoff=0.7, strand="*",
                          type="any", conservation=NULL, mc.cores=2)
GRangesTFBS <- toGRangesList(sitePairSet, axtHg19DanRer7)
GRangesTFBS$targetTFBS
GRangesTFBS$queryTFBS


###################################################
### code chunk number 14: searchBSgenome (eval = FALSE)
###################################################
## library(rtracklayer)
## library(JASPAR2014)
## library(BSgenome.Hsapiens.UCSC.hg19)
## library(BSgenome.Mmusculus.UCSC.mm10)
## pfm = getMatrixByID(JASPAR2014, ID="MA0004.1")
## pwm=toPWM(pfm)
## chain = import.chain("Downloads/hg19ToMm10.over.chain")
## sitePairSet = searchPairBSgenome(pwm, BSgenome.Hsapiens.UCSC.hg19,
##                                  BSgenome.Mmusculus.UCSC.mm10,
##                                  chr1="chr1", chr2="chr1",
##                                  min.score="90%", strand="+", chain=chain)


###################################################
### code chunk number 15: MEME-wrapper (eval = FALSE)
###################################################
## motifSet = runMEME(file.path(system.file("extdata", 
##                                          package="TFBSTools"), "crp0.s"),
##                    binary="meme",
##                    arguments=list("-nmotifs"=3)
##                   )
## sitesSeq(motifSet, type="all")
## sitesSeq(motifSet, type="none")
## consensusMatrix(motifSet)


###################################################
### code chunk number 16: session-info
###################################################
  sessionInfo()