###################################################
### chunk number 1: 
###################################################
options(SweaveHooks=list(fig=function()
          par(mar=c(3.1,2.1,4.1,2.1))))


###################################################
### chunk number 2: phyloHmm1
###################################################
require("rphast")
seqnames <- c("hg18", "panTro2", "ponAbe2", "rheMac2", "equCab2", 
              "canFam2", "dasNov2", "mm9", "rn4", "galGal3", "xenTro2") 
exampleArchive <- system.file("extdata", "examples.zip", package="rphast")
unzip(exampleArchive, "placentalMammals.mod")
neutralMod <- read.tm("placentalMammals.mod")
neutralMod$tree <- prune.tree(neutralMod$tree, seqs=seqnames, all.but=TRUE)


###################################################
### chunk number 3: phyloHmm2
###################################################
unzip(system.file("extdata", "NRSF.zip", package="rphast"))
mafFiles <- list.files("NRSF", pattern="*.maf", full.names=TRUE)
nrsfNames <- basename(mafFiles)
nrsfNames <- substr(nrsfNames, 1, nchar(nrsfNames)-4) # remove .maf suffix
nrsfSites <- read.feat("NRSF/NRSF.gff")
motifLen <- 21   
msaList <- list()
for (i in 1:length(mafFiles)) {
  smallMsa <- read.msa(mafFiles[i])  
  smallMsa <- strip.gaps.msa(smallMsa)
  smallMsa$offset <- 0 
  feat <- nrsfSites[which(nrsfSites$feature == nrsfNames[i]),]
  if (feat$strand == "-") 
    smallMsa <- reverse.complement(smallMsa)
  if (ncol.msa(smallMsa) != motifLen)
    cat("wrong number of columns in ", nrsfNames[i], "\n")
  msaList[[nrsfNames[i]]] <- smallMsa
}
aggMsa <- concat.msa(msaList)


###################################################
### chunk number 4: phyloHmm3
###################################################
feats <- feat(seqname="hg18", src=as.character(sapply(nrsfNames, rep, motifLen)),
              feature=rep(sprintf("site.%i", 1:motifLen), length(mafFiles)),
              start=1:(motifLen*length(mafFiles)),
              end=2:(motifLen*length(mafFiles)+1))
mods <- phyloFit(aggMsa, init.mod=neutralMod, no.rates=TRUE, features=feats,
                 scale.only=TRUE, ninf.sites=10)


###################################################
### chunk number 5: lrByMotifPosition
###################################################
nullLike <- numeric()
lr  <- numeric()

for (i in 1:21) {
  nullLike[i] <- likelihood.msa(aggMsa[,(0:(length(mafFiles)-1))*21+i], neutralMod)
  lr[i] <- mods[[i]]$likelihood - nullLike[i]
}
barplot(lr, names.arg=1:21, ylab="likelihood ratio")


###################################################
### chunk number 6: seqLogoPlot
###################################################
require("seqLogo")
m <- read.table("NRSF/NRSF.mtf")
pwm <- makePWM(t(m))
seqLogo(pwm, xfontsize=10)


###################################################
### chunk number 7: seqLogoPlotBackgd
###################################################
m <- matrix(0, nrow=21, ncol=4)
for (i in 1:21)
  m[i,] <- mods[[i]]$backgd
pwm <- makePWM(t(m))
seqLogo(pwm, xfontsize=10)


###################################################
### chunk number 8: phyloHmm7
###################################################
mods[["neutral"]] <- neutralMod

get.trans.mat <- function(lambda, state.names, motifLen) {
  trans.mat <- matrix(0, nrow=length(state.names), ncol=length(state.names),
                      dimnames=list(state.names, state.names))
  trans.mat["neutral", "site.1"] <- lambda
  lastSite <- sprintf("site.%i", motifLen)
  trans.mat["neutral", "neutral"] <- 1 - lambda
  for (i in 1:(motifLen-1)) 
    trans.mat[sprintf("site.%i", i), sprintf("site.%i", i+1)] <- 1
  trans.mat[lastSite, "neutral"] <- 1.0
  trans.mat
}

lambda <- 0.0001
trans.mat <- get.trans.mat(lambda, names(mods), motifLen)
eq.freq <- ifelse(names(mods) == "neutral", 0.999, (1-0.999)/motifLen)

nrsf.hmm <- hmm(trans.mat, eq.freq, begin.freq=eq.freq)


###################################################
### chunk number 9: phyloHmm8
###################################################
simLength <- 100000
simData <- simulate.msa(mods, simLength, hmm=nrsf.hmm, get.features=TRUE)


###################################################
### chunk number 10: phyloHmm9
###################################################
hmmScores <- score.hmm(msa=simData$msa, mod=mods, hmm=nrsf.hmm, viterbi=TRUE,
                       states=sprintf("site.%i", 1:motifLen))


###################################################
### chunk number 11: phyloHmm10
###################################################
predicted <- hmmScores$in.states
correct <- simData$feats[substr(simData$feats$feature, 1, 4)=="site",]

numSitePredicted <- coverage.feat(predicted)
numSiteCorrect <- coverage.feat(correct)
numSiteOverlap <- coverage.feat(predicted, correct)
cat(numSitePredicted, numSiteCorrect, numSiteOverlap, "\n")


###################################################
### chunk number 12: vignette3BrowserView
###################################################

wholeRegion <- feat("hg18", src=".", feature="all",start=1, end=simLength)
sensitivity <- coverage.feat(correct, predicted)/coverage.feat(correct)
specificity <- coverage.feat(correct, predicted, not=c(TRUE, TRUE), 
                            region.bounds=wholeRegion) /
  coverage.feat(correct, not=TRUE, region.bounds=wholeRegion)
ppv <- coverage.feat(correct, predicted, region.bounds=wholeRegion) /
  coverage.feat(predicted, region.bounds=wholeRegion)
cat(specificity, sensitivity, ppv, "\n")

tracks <- list(feat.track(correct, "actual binding sites"),
               feat.track(predicted, "predicted binding sites", col="red"),
               wig.track(hmmScores$post.prob.wig$coord,
                         hmmScores$post.prob.wig$post.prob,
                         "hmm Posterior probilities",
                         smooth=FALSE, col="red", ylim=c(0,1)))
plot.track(tracks)


###################################################
### chunk number 13: phyloHmm12
###################################################
unlink("NRSF", recursive=TRUE)


