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

###################################################
### code chunk number 1: Init
###################################################
options(width=70)
options(useFancyQuotes=FALSE)
set.seed(0)
library(kebabs)
kebabsVersion <- packageDescription("kebabs")$Version
kebabsDateRaw <- packageDescription("kebabs")$Date
kebabsDateYear <- as.numeric(substr(kebabsDateRaw, 1, 4))
kebabsDateMonth <- as.numeric(substr(kebabsDateRaw, 6, 7))
kebabsDateDay <- as.numeric(substr(kebabsDateRaw, 9, 10))
kebabsDate <- paste(month.name[kebabsDateMonth], " ",
                     kebabsDateDay, ", ",
                     kebabsDateYear, sep="")


###################################################
### code chunk number 2: KeBABS.Rnw:136-138 (eval = FALSE)
###################################################
## source("http://bioconductor.org/biocLite.R")
## biocLite("kebabs")


###################################################
### code chunk number 3: KeBABS.Rnw:165-166
###################################################
library(kebabs)


###################################################
### code chunk number 4: KeBABS.Rnw:173-174 (eval = FALSE)
###################################################
## vignette("kebabs")


###################################################
### code chunk number 5: KeBABS.Rnw:179-180 (eval = FALSE)
###################################################
## help(kebabs)


###################################################
### code chunk number 6: KeBABS.Rnw:187-188
###################################################
data(TFBS)


###################################################
### code chunk number 7: KeBABS.Rnw:193-195
###################################################
enhancerFB
length(enhancerFB)


###################################################
### code chunk number 8: KeBABS.Rnw:200-202 (eval = FALSE)
###################################################
## hist(width(enhancerFB), breaks=30, xlab="Sequence Length",
##      main="Distribution of Sequence Lengths")


###################################################
### code chunk number 9: KeBABS.Rnw:205-209
###################################################
pdf("001.pdf")
hist(width(enhancerFB), breaks=30, main="Distribution of Sequence Lenghts",
        xlab="Sequence Length")
dev.off()


###################################################
### code chunk number 10: KeBABS.Rnw:220-221
###################################################
showAnnotatedSeq(enhancerFB, sel=3)


###################################################
### code chunk number 11: KeBABS.Rnw:226-227
###################################################
head(yFB)


###################################################
### code chunk number 12: KeBABS.Rnw:232-233
###################################################
table(yFB)


###################################################
### code chunk number 13: KeBABS.Rnw:240-244
###################################################
numSamples <- length(enhancerFB)
trainingFraction <- 0.7
train <- sample(1:numSamples, trainingFraction * numSamples)
test <- c(1:numSamples)[-train]


###################################################
### code chunk number 14: KeBABS.Rnw:251-252
###################################################
specK2 <- spectrumKernel(k=2)


###################################################
### code chunk number 15: KeBABS.Rnw:259-261
###################################################
model <- kbsvm(x=enhancerFB[train], y=yFB[train], kernel=specK2,
               pkg="e1071", svm="C-svc", cost=15)


###################################################
### code chunk number 16: KeBABS.Rnw:268-271
###################################################
pred <- predict(model, enhancerFB[test])
head(pred)
head(yFB[test])


###################################################
### code chunk number 17: KeBABS.Rnw:276-277
###################################################
evaluatePrediction(pred, yFB[test], allLabels=unique(yFB))


###################################################
### code chunk number 18: KeBABS.Rnw:284-290
###################################################
gappyK1M3  <- gappyPairKernel(k=1, m=3)
model <- kbsvm(x=enhancerFB[train], y=yFB[train],
               kernel=gappyK1M3, pkg="e1071", svm="C-svc",
               cost=15)
pred <- predict(model, enhancerFB[test])
evaluatePrediction(pred, yFB[test], allLabels=unique(yFB))


###################################################
### code chunk number 19: KeBABS.Rnw:295-301
###################################################
gappyK1M3  <- gappyPairKernel(k=1, m=3)
model <- kbsvm(x=enhancerFB[train], y=yFB[train],
               kernel=gappyK1M3, pkg="LiblineaR", svm="C-svc",
               cost=15)
pred <- predict(model, enhancerFB[test])
evaluatePrediction(pred, yFB[test], allLabels=unique(yFB))


###################################################
### code chunk number 20: KeBABS.Rnw:308-314
###################################################
gappyK1M3  <- gappyPairKernel(k=1, m=3)
model <- kbsvm(x=enhancerFB[train], y=yFB[train],
               kernel=gappyK1M3, pkg="LiblineaR", svm="l1rl2l-svc",
               cost=5)
pred <- predict(model, enhancerFB[test])
evaluatePrediction(pred, yFB[test], allLabels=unique(yFB))


###################################################
### code chunk number 21: KeBABS.Rnw:319-320 (eval = FALSE)
###################################################
## ?kbsvm


###################################################
### code chunk number 22: KeBABS.Rnw:418-419
###################################################
specK2 <- spectrumKernel(k=2)


###################################################
### code chunk number 23: KeBABS.Rnw:424-425
###################################################
specK2 <- spectrumKernel(k=2, normalized=FALSE)


###################################################
### code chunk number 24: KeBABS.Rnw:451-452
###################################################
mismK3M1 <- mismatchKernel(k=3, m=1)


###################################################
### code chunk number 25: KeBABS.Rnw:468-469
###################################################
gappyK1M2 <- gappyPairKernel(k=1, m=2)


###################################################
### code chunk number 26: KeBABS.Rnw:491-494
###################################################
motifCollection1 <- c("A[CG]T","C.G","C..G.T","G[^A][AT]",
                       "GT.A[CA].[CT]G")
motif1 <- motifKernel(motifCollection1)


###################################################
### code chunk number 27: KeBABS.Rnw:532-534
###################################################
gappyK1M2ps <- gappyPairKernel(k=1, m=2, distWeight=1,
                                 normalized=FALSE)


###################################################
### code chunk number 28: KeBABS.Rnw:541-544
###################################################
seq1 <- AAStringSet(c("GACGAGGACCGA","AGTAGCGAGGT","ACGAGGTCTTT",
                      "GGACCGAGTCGAGG"))
positionMetadata(seq1) <- c(3, 5, 2, 10)


###################################################
### code chunk number 29: KeBABS.Rnw:551-555
###################################################
wdK3 <- spectrumKernel(k=3, distWeight=1, mixCoef=c(0.5,0.33,0.17),
                         normalized=FALSE)
km <- getKernelMatrix(wdK3, seq1)
km


###################################################
### code chunk number 30: KeBABS.Rnw:560-563
###################################################
positionMetadata(seq1) <- NULL
km <- getKernelMatrix(wdK3, seq1)
km


###################################################
### code chunk number 31: KeBABS.Rnw:606-615 (eval = FALSE)
###################################################
## curve(linWeight(x, sigma=5), from=-25, to=25, xlab="p - q", ylab="weight",
##       main="Predefined Distance Weighting Functions", col="green")
## curve(expWeight(x, sigma=5), from=-25, to=25, col="blue", add=TRUE)
## curve(gaussWeight(x, sigma=5), from=-25, to=25, col="red", add=TRUE)
## curve(swdWeight(x), from=-25, to=25, col="orange", add=TRUE)
## legend('topright', inset=0.03, title="Weighting Functions", c("linWeight",
##        "expWeight", "gaussWeight", "swdWeight"), fill=c("green", "blue",
##        "red", "orange"))
## text(17, 0.70, "sigma = 5")


###################################################
### code chunk number 32: KeBABS.Rnw:618-629
###################################################
pdf("002.pdf")
curve(linWeight(x, sigma=5), from=-25, to=25, xlab="p - q", ylab="weight",
      main="Predefined Distance Weighting Functions", col="green")
curve(expWeight(x, sigma=5), from=-25, to=25, col="blue", add=TRUE)
curve(gaussWeight(x, sigma=5), from=-25, to=25, col="red", add=TRUE)
curve(swdWeight(x), from=-25, to=25, col="orange", add=TRUE)
legend('topright', inset=0.03, title="Weighting Functions", c("linWeight",
       "expWeight", "gaussWeight", "swdWeight"), fill=c("green", "blue",
       "red", "orange"))
text(17, 0.70, "sigma = 5")
dev.off()


###################################################
### code chunk number 33: KeBABS.Rnw:646-655
###################################################
swdWeight <- function(d)
{
    if (missing(d))
        return(function(d) swdWeight(d))

    1 / (2 * (abs(d) + 1))
}
swdK3 <- spectrumKernel(k=3, distWeight=swdWeight(),
                        mixCoef=c(0.5,0.33,0.17))


###################################################
### code chunk number 34: KeBABS.Rnw:660-665
###################################################
data(TFBS)
names(enhancerFB) <- paste("Sample", 1:length(enhancerFB), sep="_")
enhancerFB
kmSWD <- getKernelMatrix(swdK3, x=enhancerFB, selx=1:5)
kmSWD[1:5,1:5]


###################################################
### code chunk number 35: KeBABS.Rnw:679-690
###################################################
udWeight <- function(d, base=2)
{
    if (missing(d))
        return(function(d) udWeight(d, base=base))

    return(base^(-d))
}

specudK3 <- spectrumKernel(k=3, distWeight=udWeight(base=4),
                                mixCoef=c(0,0.3,0.7))
kmud <- getKernelMatrix(specudK3, x=enhancerFB, selx=1:5)


###################################################
### code chunk number 36: KeBABS.Rnw:723-787 (eval = FALSE)
###################################################
## getGenesWithExonIntronAnnotation <- function(geneList, genomelib,
##                                              txlib)
## {
##     library(BSgenome)
##     library(genomelib, character.only=TRUE)
##     library(txlib, character.only=TRUE)
##     genome <- getBSgenome(genomelib)
##     txdb <- eval(parse(text=txlib))
##     exonsByGene <- exonsBy(txdb, by ="gene")
## 
##     ## generate exon/intron annotation
##     annot <- rep("", length(geneList))
##     geneRanges <- GRanges()
##     exonsSelGenes <- exonsByGene[geneList]
## 
##     if (length(exonsSelGenes) != length(geneList))
##         stop("some genes are not found")
## 
##     for (i in 1:length(geneList))
##     {
##         exons <- unlist(exonsSelGenes[i])
##         exonRanges <- ranges(exons)
##         chr <-  as.character(seqnames(exons)[1])
##         strand <- as.character(strand(exons)[1])
##         numExons <- length(width(exonRanges))
## 
##         for (j in 1:numExons)
##         {
##             annot[i] <-
##                 paste(annot[i],
##                       paste(rep("e", width(exonRanges)[j]),
##                             collapse=""), sep="")
## 
##             if (j < numExons)
##             {
##                 annot[i] <-
##                     paste(annot[i],
##                           paste(rep("i", start(exonRanges)[j+1] -
##                                 end(exonRanges)[j] - 1),
##                           collapse=""), sep="")
##             }
##         }
## 
##         geneRanges <-
##             c(geneRanges,
##               GRanges(seqnames=Rle(chr),
##                       strand=Rle(strand(strand)),
##                       ranges=IRanges(start=start(exonRanges)[1],
##                                      end=end(exonRanges)[numExons])))
##     }
## 
##     ## retrieve gene sequences
##     seqs <- getSeq(genome, geneRanges)
##     names(seqs) <- geneList
##     ## assign annotation
##     annotationMetadata(seqs, annCharset="ei") <- annot
##     seqs
## }
## 
## ## get gene sequences for HBA1 and HBA2 with exon/intron annotation
## ## 3039 and 3040 are the geneID values for HBA1 and HBA2
## hba <- getGenesWithExonIntronAnnotation(c("3039", "3040"),
##            "BSgenome.Hsapiens.UCSC.hg19",
##            "TxDb.Hsapiens.UCSC.hg19.knownGene")


###################################################
### code chunk number 37: KeBABS.Rnw:792-794 (eval = FALSE)
###################################################
## annotationCharset(hba)
## showAnnotatedSeq(hba, sel=1, start=1, end=400)


###################################################
### code chunk number 38: KeBABS.Rnw:799-805 (eval = FALSE)
###################################################
## specK2 <- spectrumKernel(k=2)
## specK2a <- spectrumKernel(k=2, annSpec=TRUE)
## erK2 <- getExRep(hba, specK2, sparse=FALSE)
## erK2[,1:6]
## erK2a <- getExRep(hba, specK2a, sparse=FALSE)
## erK2a[,1:6]


###################################################
### code chunk number 39: KeBABS.Rnw:810-814 (eval = FALSE)
###################################################
## km <- linearKernel(erK2)
## km
## kma <- linearKernel(erK2a)
## kma


###################################################
### code chunk number 40: KeBABS.Rnw:823-838
###################################################
data(CCoil)
ccseq
ccannot[1:3]
head(yCC)
yCC <- as.numeric(yCC)
## delete annotation metadata
annotationMetadata(ccseq) <- NULL
annotationMetadata(ccseq)
gappy <- gappyPairKernel(k=1, m=10)
train <- sample(1:length(ccseq), 0.8*length(ccseq))
test <- c(1:length(ccseq))[-train]
model <- kbsvm(ccseq[train], y=yCC[train], kernel=gappy,
               pkg="LiblineaR", svm="C-svc", cost=100)
pred <- predict(model, ccseq[test])
evaluatePrediction(pred, yCC[test], allLabels=unique(yCC))


###################################################
### code chunk number 41: KeBABS.Rnw:843-853
###################################################
## assign annotation metadata
annotationMetadata(ccseq, annCharset="abcdefg") <- ccannot
annotationMetadata(ccseq)[1:5]
annotationCharset(ccseq)
showAnnotatedSeq(ccseq, sel=2)
gappya <- gappyPairKernel(k=1, m=10, annSpec=TRUE)
model <- kbsvm(ccseq[train], y=yCC[train], kernel=gappya,
               pkg="LiblineaR", svm="C-svc", cost=100)
pred <- predict(model, ccseq[test])
evaluatePrediction(pred, yCC[test], allLabels=unique(yCC))


###################################################
### code chunk number 42: KeBABS.Rnw:858-867
###################################################
## grid search with two kernels and 6 hyperparameter values
## using the balanced accuracy as performance objective
model <- kbsvm(ccseq[train], y=yCC[train],
               kernel=c(gappy, gappya), pkg="LiblineaR", svm="C-svc",
               cost=c(1,10,50,100,200,500), explicit="yes", cross=5,
               perfParameters="ALL", perfObjective="BACC",
               showProgress=TRUE)
result <- modelSelResult(model)
result


###################################################
### code chunk number 43: KeBABS.Rnw:872-876
###################################################
perfData <- performance(result)
perfData
which(perfData$BACC[1,] == max(perfData$BACC[1,]))
which(perfData$BACC[2,] == max(perfData$BACC[2,]))


###################################################
### code chunk number 44: KeBABS.Rnw:881-882 (eval = FALSE)
###################################################
## plot(result, sel="BACC")


###################################################
### code chunk number 45: KeBABS.Rnw:885-888
###################################################
pdf("005.pdf")
plot(result, sel="BACC")
dev.off()


###################################################
### code chunk number 46: KeBABS.Rnw:906-951
###################################################
## position-independent spectrum kernel
specK2 <- spectrumKernel(k=3)
#
## annotation specific spectrum normalized
specK2a <- spectrumKernel(k=3, annSpec=TRUE)
#
## spectrum kernel with presence normalized
specK2p <- spectrumKernel(k=3, presence=TRUE)
#
## mixed spectrum normalized
specK2m <- spectrumKernel(k=3, mixCoef=c(0.5, 0.33, 0.17))
#
## position-specific spectrum normalized
specK2ps <- spectrumKernel(k=3, distWeight=1)
#
## mixed position-specific spectrum kernel normalized
## also called weighted degree kernel normalized
specK2wd <- spectrumKernel(k=3, dist=1,
                             mixCoef=c(0.5, 0.33, 0.17))
#
## distance-weighted spectrum normalized
specK2lin <- spectrumKernel(k=3, distWeight=linWeight(sigma=10))
specK2exp <- spectrumKernel(k=3, distWeight=expWeight(sigma=10))
specK2gs <- spectrumKernel(k=3, distWeight=gaussWeight(sigma=10))
#
## shifted weighted degree with equal position weighting normalized
specK2swd <- spectrumKernel(k=3, distWeight=swdWeight(),
                              mixCoef=c(0.5,0.33,0.17))
#
## distance-weighted spectrum kernel with user defined distance
## weighting
udWeight <- function(d, base=2)
{
    if (!(is.numeric(base) && length(base==1)))
        stop("parameter 'base' must be a single numeric value\n")

    if (missing(d))
        return(function(d) udWeight(d, base=base))

    if (!is.numeric(d))
        stop("'d' must be a numeric vector\n")

    return(base^(-d))
}
specK2ud <- spectrumKernel(k=3, distWeight=udWeight(b=2))


###################################################
### code chunk number 47: KeBABS.Rnw:960-968
###################################################
specK25 <- spectrumKernel(k=2:5)
specK25
train <- 1:100
model <- kbsvm(x=enhancerFB[train], y=yFB[train],
               kernel=specK25, pkg="LiblineaR", svm="C-svc",
               cost=c(1,5,10,20,50,100), cross=5, explicit="yes",
               showProgress=TRUE)
modelSelResult(model)


###################################################
### code chunk number 48: KeBABS.Rnw:973-975
###################################################
kernelList1 <- list(spectrumKernel(k=3), mismatchKernel(k=3,m=1),
                    gappyPairKernel(k=2,m=4))


###################################################
### code chunk number 49: KeBABS.Rnw:980-981
###################################################
kernelList2 <- c(spectrumKernel(k=2:4), gappyPairKernel(k=1, m=2:5))


###################################################
### code chunk number 50: KeBABS.Rnw:1000-1005
###################################################
specK2 <- spectrumKernel(k=2)
km <- getKernelMatrix(specK2, x=enhancerFB)
class(km)
dim(km)
km[1:3, 1:3]


###################################################
### code chunk number 51: KeBABS.Rnw:1010-1011 (eval = FALSE)
###################################################
## heatmap(km, symm=TRUE)


###################################################
### code chunk number 52: KeBABS.Rnw:1029-1032
###################################################
specK2 <- spectrumKernel(k=2)
km <- specK2(x=enhancerFB)
km[1:3, 1:3]


###################################################
### code chunk number 53: KeBABS.Rnw:1037-1039
###################################################
km <- getKernelMatrix(specK2, x=enhancerFB, selx=c(1,4,25,137,300))
km


###################################################
### code chunk number 54: KeBABS.Rnw:1044-1049
###################################################
seqs1 <- enhancerFB[1:200]
seqs2 <- enhancerFB[201:500]
km <- getKernelMatrix(specK2, x=seqs1, y=seqs2)
dim(km)
km[1:4,1:5]


###################################################
### code chunk number 55: KeBABS.Rnw:1054-1057
###################################################
km <- getKernelMatrix(specK2, x=enhancerFB, selx=1:200,
                      y=enhancerFB, sely=201:500)
dim(km)


###################################################
### code chunk number 56: KeBABS.Rnw:1074-1077
###################################################
specK2 <- spectrumKernel(k=2, normalized=FALSE)
erd <- getExRep(enhancerFB, selx=1:5, kernel=specK2, sparse=FALSE)
erd


###################################################
### code chunk number 57: KeBABS.Rnw:1082-1087
###################################################
specK6 <- spectrumKernel(k=6, normalized=FALSE)
erd <- getExRep(enhancerFB, selx=1:5, kernel=specK6,
                sparse=FALSE)
dim(erd)
erd[,1:6]


###################################################
### code chunk number 58: KeBABS.Rnw:1092-1100
###################################################
specK6 <- spectrumKernel(k=6, normalized=FALSE)
erd <- getExRep(enhancerFB, kernel=specK6, sparse=FALSE)
dim(erd)
object.size(erd)
ers <- getExRep(enhancerFB, kernel=specK6, sparse=TRUE)
dim(ers)
object.size(ers)
ers[1:5, 1:6]


###################################################
### code chunk number 59: KeBABS.Rnw:1114-1119
###################################################
library(apcluster)
gappyK1M4 <- gappyPairKernel(k=1, m=4)
km <- getKernelMatrix(gappyK1M4, enhancerFB)
apres <- apcluster(s=km, p=0.8)
length(apres)


###################################################
### code chunk number 60: KeBABS.Rnw:1124-1126 (eval = FALSE)
###################################################
## aggres <- aggExCluster(km, apres)
## plot(aggres)


###################################################
### code chunk number 61: KeBABS.Rnw:1129-1133
###################################################
aggres <- aggExCluster(km, apres)
pdf("004.pdf")
plot(aggres)
dev.off()


###################################################
### code chunk number 62: KeBABS.Rnw:1144-1147
###################################################
exrep <- getExRep(enhancerFB, gappyK1M4, sparse=FALSE)
apres1 <- apcluster(s=linearKernel, x=exrep, p=0.1)
length(apres1)


###################################################
### code chunk number 63: KeBABS.Rnw:1159-1164
###################################################
exrep <- getExRep(x=enhancerFB, selx=1:5, gappyK1M4, sparse=FALSE)
dim(exrep)
erquad <- getExRepQuadratic(exrep)
dim(erquad)
erquad[1:5,1:5]


###################################################
### code chunk number 64: KeBABS.Rnw:1179-1189
###################################################
gappyK1M4 <- gappyPairKernel(k=1, m=4)
exrep <- getExRep(enhancerFB, gappyK1M4, sparse=TRUE)
numSamples <- length(enhancerFB)
trainingFraction <- 0.7
train <- sample(1:numSamples, trainingFraction * numSamples)
test <- c(1:numSamples)[-train]
model <- kbsvm(x=exrep[train, ], y=yFB[train], kernel=gappyK1M4,
               pkg="LiblineaR", svm="C-svc", cost=15)
pred <- predict(model, exrep[test, ])
evaluatePrediction(pred, yFB[test], allLabels=unique(yFB))


###################################################
### code chunk number 65: KeBABS.Rnw:1200-1202
###################################################
preddec <- predict(model, exrep[test, ], predictionType="decision")
evaluatePrediction(pred, yFB[test], allLabels=unique(yFB), decValues=preddec)


###################################################
### code chunk number 66: KeBABS.Rnw:1207-1209 (eval = FALSE)
###################################################
## rocdata <- computeROCandAUC(preddec, yFB[test], unique(yFB))
## plot(rocdata, main="Receiver Operating Characteristics")


###################################################
### code chunk number 67: KeBABS.Rnw:1212-1216
###################################################
rocdata <- computeROCandAUC(preddec, yFB[test], unique(yFB))
pdf("008.pdf")
plot(rocdata, main="Receiver Operating Characteristics")
dev.off()


###################################################
### code chunk number 68: KeBABS.Rnw:1270-1280
###################################################
data(CCoil)
ccseq
head(yCC)
head(ccgroups)
gappyK1M6 <- gappyPairKernel(k=1, m=6)
model <- kbsvm(x=ccseq, y=as.numeric(yCC), kernel=gappyK1M6,
               pkg="LiblineaR", svm="C-svc", cost=30, cross=3,
               noCross=2, groupBy=ccgroups, perfObjective="BACC",
               perfParameters=c("ACC", "BACC"))
cvResult(model)


###################################################
### code chunk number 69: KeBABS.Rnw:1288-1296
###################################################
specK24 <- spectrumKernel(k=2:4)
gappyK1M24 <- gappyPairKernel(k=1, m=2:4)
gridKernels <- c(specK24, gappyK1M24)
cost <- c(1,10, 100, 1000, 10000)
model <- kbsvm(x=enhancerFB, y=yFB, kernel=gridKernels,
               pkg="LiblineaR", svm="C-svc", cost=cost, cross=3,
               explicit="yes", showProgress=TRUE)
modelSelResult(model)


###################################################
### code chunk number 70: KeBABS.Rnw:1305-1316
###################################################
specK34 <- spectrumKernel(k=3:4)
gappyK1M34 <- gappyPairKernel(k=1, m=3:4)
gridKernels <- c(specK34, gappyK1M34)
pkgs <- c("e1071", "LiblineaR", "LiblineaR")
svms <- c("C-svc","C-svc","l1rl2l-svc")
cost <- c(50, 50, 12)
model <- kbsvm(x=enhancerFB, y=yFB, kernel=gridKernels,
               pkg=pkgs, svm=svms, cost=cost, cross=10,
               explicit="yes", showProgress=TRUE,
               showCVTimes=TRUE)
modelSelResult(model)


###################################################
### code chunk number 71: KeBABS.Rnw:1328-1337
###################################################
specK34 <- spectrumKernel(k=3:4)
gappyK1M34 <- gappyPairKernel(k=1, m=3:4)
gridKernels <- c(specK34, gappyK1M34)
cost <- c(10, 50, 100)
model <- kbsvm(x=enhancerFB, y=yFB, kernel=gridKernels,
               pkg="LiblineaR", svm="C-svc", cost=cost, cross=10,
               explicit="yes", nestedCross=4)
modelSelResult(model)
cvResult(model)


###################################################
### code chunk number 72: KeBABS.Rnw:1348-1354
###################################################
head(yReg)
gappyK1M2 <- gappyPairKernel(k=1, m=2)
model <- kbsvm(x=enhancerFB, y=yReg, kernel=gappyK1M2,
               pkg="e1071", svm="nu-svr", nu=c(0.5,0.6,0.7,0.8),
               cross=10, showProgress=TRUE)
modelSelResult(model)


###################################################
### code chunk number 73: KeBABS.Rnw:1359-1371
###################################################
numSamples <- length(enhancerFB)
trainingFraction <- 0.7
train <- sample(1:numSamples, trainingFraction * numSamples)
test <- c(1:numSamples)[-train]
model <- kbsvm(x=enhancerFB[train], y=yReg[train],
               kernel=gappyK1M2, pkg="e1071", svm="nu-svr",
               nu=0.7)
pred <- predict(model, enhancerFB[test])
mse <- sum((yReg[test] - pred)^2)/length(test)
mse
featWeights <- featureWeights(model)
colnames(featWeights)[which(featWeights > 0.4)]


###################################################
### code chunk number 74: KeBABS.Rnw:1376-1382
###################################################
model <- kbsvm(x=enhancerFB[train], y=yReg[train],
               kernel=spectrumKernel(k=2), pkg="e1071", svm="nu-svr",
               nu=0.7)
pred <- predict(model, enhancerFB[test])
featWeights <- featureWeights(model)
colnames(featWeights)[which(featWeights > 0.4)]


###################################################
### code chunk number 75: KeBABS.Rnw:1387-1398
###################################################
model <- kbsvm(x=enhancerFB[train], y=yReg[train],
               kernel=spectrumKernel(k=2), pkg="e1071", svm="nu-svr",
               nu=c(0.5,0.55,0.6), cross=10, nestedCross=5)
modelSelResult(model)
cvResult(model)
model <- kbsvm(x=enhancerFB[train], y=yReg[train],
               kernel=gappyPairKernel(k=1,m=2), pkg="e1071",
               svm="nu-svr", nu=c(0.6,0.65,0.7), cross=10,
               nestedCross=5)
modelSelResult(model)
cvResult(model)


###################################################
### code chunk number 76: KeBABS.Rnw:1464-1481
###################################################
data(CCoil)
gappya <- gappyPairKernel(k=1,m=11, annSpec=TRUE)
model <- kbsvm(x=ccseq, y=as.numeric(yCC), kernel=gappya,
               pkg="e1071", svm="C-svc", cost=15)
featureWeights(model)[,1:5]
GCN4 <- AAStringSet(c("MKQLEDKVEELLSKNYHLENEVARLKKLV",
                      "MKQLEDKVEELLSKYYHTENEVARLKKLV"))
names(GCN4) <- c("GCN4wt", "GCN_N16Y,L19T")
annCharset <- annotationCharset(ccseq)
annot <- c("abcdefgabcdefgabcdefgabcdefga",
           "abcdefgabcdefgabcdefgabcdefga")

annotationMetadata(GCN4, annCharset=annCharset) <- annot
predProf <- getPredictionProfile(GCN4, gappya,
                                 featureWeights(model),
                                 modelOffset(model))
predProf


###################################################
### code chunk number 77: KeBABS.Rnw:1486-1487 (eval = FALSE)
###################################################
## plot(predProf, sel=1, ylim=c(-0.4, 0.2), heptads=TRUE, annotate=TRUE)


###################################################
### code chunk number 78: KeBABS.Rnw:1490-1493
###################################################
pdf("006.pdf")
plot(predProf, sel=1, ylim=c(-0.4, 0.2), heptads=TRUE, annotate=TRUE)
dev.off()


###################################################
### code chunk number 79: KeBABS.Rnw:1504-1505 (eval = FALSE)
###################################################
## plot(predProf, sel=c(1,2), ylim=c(-0.4, 0.2), heptads=TRUE, annotate=TRUE)


###################################################
### code chunk number 80: KeBABS.Rnw:1508-1511
###################################################
pdf("007.pdf")
plot(predProf, sel=c(1,2), ylim=c(-0.4, 0.2), heptads=TRUE, annotate=TRUE)
dev.off()


###################################################
### code chunk number 81: KeBABS.Rnw:1532-1539
###################################################
table(yMC)
gappyK1M2 <- gappyPairKernel(k=1, m=2)
model <- kbsvm(x=enhancerFB[train], y=yMC[train],
               kernel=gappyK1M2, pkg="LiblineaR",
               svm="C-svc", cost=300)
pred <- predict(model, enhancerFB[test])
evaluatePrediction(pred, yMC[test], allLabels=unique(yMC))


###################################################
### code chunk number 82: KeBABS.Rnw:1544-1549
###################################################
featWeights <- featureWeights(model)
length(featWeights)
featWeights[[1]][1:5]
featWeights[[2]][1:5]
featWeights[[3]][1:5]


###################################################
### code chunk number 83: KeBABS.Rnw:1554-1558
###################################################
predProf <- getPredictionProfile(enhancerFB, gappyK1M2,
                                 featureWeights(model)[[2]],
                                 modelOffset(model)[2])
predProf


###################################################
### code chunk number 84: KeBABS.Rnw:1594-1595 (eval = FALSE)
###################################################
## toBibtex(citation("kebabs"))