### 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:159-160
###################################################
library(kebabs)


###################################################
### code chunk number 4: kebabs.Rnw:167-168 (eval = FALSE)
###################################################
## vignette("kebabs")


###################################################
### code chunk number 5: kebabs.Rnw:173-174 (eval = FALSE)
###################################################
## help(kebabs)


###################################################
### code chunk number 6: kebabs.Rnw:181-182
###################################################
data(TFBS)


###################################################
### code chunk number 7: kebabs.Rnw:187-189
###################################################
enhancerFB
length(enhancerFB)


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


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


###################################################
### code chunk number 10: kebabs.Rnw:214-215
###################################################
showAnnotatedSeq(enhancerFB, sel=3)


###################################################
### code chunk number 11: kebabs.Rnw:220-221
###################################################
head(yFB)


###################################################
### code chunk number 12: kebabs.Rnw:226-227
###################################################
table(yFB)


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


###################################################
### code chunk number 14: kebabs.Rnw:245-246
###################################################
specK2 <- spectrumKernel(k=2)


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


###################################################
### code chunk number 16: kebabs.Rnw:262-265
###################################################
pred <- predict(model, enhancerFB[test])
head(pred)
head(yFB[test])


###################################################
### code chunk number 17: kebabs.Rnw:270-271
###################################################
evaluatePrediction(pred, yFB[test], allLabels=unique(yFB))


###################################################
### code chunk number 18: kebabs.Rnw:278-284
###################################################
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:289-295
###################################################
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:302-308
###################################################
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:313-314 (eval = FALSE)
###################################################
## ?kbsvm


###################################################
### code chunk number 22: kebabs.Rnw:412-413
###################################################
specK2 <- spectrumKernel(k=2)


###################################################
### code chunk number 23: kebabs.Rnw:418-419
###################################################
specK2 <- spectrumKernel(k=2, normalized=FALSE)


###################################################
### code chunk number 24: kebabs.Rnw:445-446
###################################################
mismK3M1 <- mismatchKernel(k=3, m=1)


###################################################
### code chunk number 25: kebabs.Rnw:462-463
###################################################
gappyK1M2 <- gappyPairKernel(k=1, m=2)


###################################################
### code chunk number 26: kebabs.Rnw:485-488
###################################################
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:526-528
###################################################
gappyK1M2ps <- gappyPairKernel(k=1, m=2, distWeight=1,
                                 normalized=FALSE)


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


###################################################
### code chunk number 29: kebabs.Rnw:545-549
###################################################
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:554-557
###################################################
positionMetadata(seq1) <- NULL
km <- getKernelMatrix(wdK3, seq1)
km


###################################################
### code chunk number 31: kebabs.Rnw:600-609 (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:612-623
###################################################
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:640-649
###################################################
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:654-659
###################################################
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:673-684
###################################################
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:717-781 (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:786-788 (eval = FALSE)
###################################################
## annotationCharset(hba)
## showAnnotatedSeq(hba, sel=1, start=1, end=400)


###################################################
### code chunk number 38: kebabs.Rnw:793-799 (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:804-808 (eval = FALSE)
###################################################
## km <- linearKernel(erK2)
## km
## kma <- linearKernel(erK2a)
## kma


###################################################
### code chunk number 40: kebabs.Rnw:817-832
###################################################
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:837-847
###################################################
## 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:852-861
###################################################
## 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:866-870
###################################################
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:875-876 (eval = FALSE)
###################################################
## plot(result, sel="BACC")


###################################################
### code chunk number 45: kebabs.Rnw:879-882
###################################################
pdf("005.pdf")
plot(result, sel="BACC")
dev.off()


###################################################
### code chunk number 46: kebabs.Rnw:900-945
###################################################
## position-independent spectrum kernel normalized
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:954-962
###################################################
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:967-969
###################################################
kernelList1 <- list(spectrumKernel(k=3), mismatchKernel(k=3,m=1),
                    gappyPairKernel(k=2,m=4))


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


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


###################################################
### code chunk number 51: kebabs.Rnw:1004-1005 (eval = FALSE)
###################################################
## heatmap(km, symm=TRUE)


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


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


###################################################
### code chunk number 54: kebabs.Rnw:1038-1043
###################################################
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:1048-1051
###################################################
km <- getKernelMatrix(specK2, x=enhancerFB, selx=1:200,
                      y=enhancerFB, sely=201:500)
dim(km)


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


###################################################
### code chunk number 57: kebabs.Rnw:1076-1081
###################################################
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:1086-1094
###################################################
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:1108-1113
###################################################
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:1118-1120 (eval = FALSE)
###################################################
## aggres <- aggExCluster(km, apres)
## plot(aggres)


###################################################
### code chunk number 61: kebabs.Rnw:1123-1127
###################################################
aggres <- aggExCluster(km, apres)
pdf("004.pdf")
plot(aggres)
dev.off()


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


###################################################
### code chunk number 63: kebabs.Rnw:1153-1158
###################################################
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:1173-1183
###################################################
gappyK1M4 <- gappyPairKernel(k=1, m=4)
exrep <- getExRep(enhancerFB, gappyK1M4, sparse=FALSE)
numSamples <- length(enhancerFB)
trainingFraction <- 0.8
train <- sample(1:numSamples, trainingFraction * numSamples)
test <- c(1:numSamples)[-train]
model <- kbsvm(x=exrep[train, ], y=yFB[train], kernel=gappyK1M4,
               pkg="kernlab", svm="C-svc", cost=15)
pred <- predict(model, exrep[test, ])
evaluatePrediction(pred, yFB[test], allLabels=unique(yFB))


###################################################
### code chunk number 65: kebabs.Rnw:1188-1197
###################################################
## compute quadratic kernel matrix for training samples
kmtrain <- getKernelMatrix(gappyK1M4, x=enhancerFB, selx=train)
model1 <- kbsvm(x=kmtrain, y=yFB[train], kernel=gappyK1M4,
                pkg="kernlab", svm="C-svc", cost=15)
## compute rectangular kernel matrix of test samples against support vectors
kmtest <- getKernelMatrix(gappyK1M4, x=enhancerFB, y=enhancerFB, selx=test, 
                          sely=train[SVindex(model1)])
pred1 <- predict(model1, kmtest)
evaluatePrediction(pred1, yFB[test], allLabels=unique(yFB))


###################################################
### code chunk number 66: kebabs.Rnw:1208-1210
###################################################
perf <- evaluatePrediction(pred, yFB[test], allLabels=unique(yFB), print=FALSE)
perf


###################################################
### code chunk number 67: kebabs.Rnw:1257-1267
###################################################
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 68: kebabs.Rnw:1275-1283
###################################################
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 69: kebabs.Rnw:1292-1303
###################################################
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 70: kebabs.Rnw:1315-1324
###################################################
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 71: kebabs.Rnw:1335-1341
###################################################
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 72: kebabs.Rnw:1346-1358
###################################################
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 73: kebabs.Rnw:1363-1369
###################################################
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 74: kebabs.Rnw:1374-1385
###################################################
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 75: kebabs.Rnw:1451-1468
###################################################
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 76: kebabs.Rnw:1473-1474 (eval = FALSE)
###################################################
## plot(predProf, sel=1, ylim=c(-0.4, 0.2), heptads=TRUE, annotate=TRUE)


###################################################
### code chunk number 77: kebabs.Rnw:1477-1480
###################################################
pdf("006.pdf")
plot(predProf, sel=1, ylim=c(-0.4, 0.2), heptads=TRUE, annotate=TRUE)
dev.off()


###################################################
### code chunk number 78: kebabs.Rnw:1491-1492 (eval = FALSE)
###################################################
## plot(predProf, sel=c(1,2), ylim=c(-0.4, 0.2), heptads=TRUE, annotate=TRUE)


###################################################
### code chunk number 79: kebabs.Rnw:1495-1498
###################################################
pdf("007.pdf")
plot(predProf, sel=c(1,2), ylim=c(-0.4, 0.2), heptads=TRUE, annotate=TRUE)
dev.off()


###################################################
### code chunk number 80: kebabs.Rnw:1519-1526
###################################################
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 81: kebabs.Rnw:1531-1536
###################################################
featWeights <- featureWeights(model)
length(featWeights)
featWeights[[1]][1:5]
featWeights[[2]][1:5]
featWeights[[3]][1:5]


###################################################
### code chunk number 82: kebabs.Rnw:1541-1545
###################################################
predProf <- getPredictionProfile(enhancerFB, gappyK1M2,
                                 featureWeights(model)[[2]],
                                 modelOffset(model)[2])
predProf


###################################################
### code chunk number 83: kebabs.Rnw:1627-1628 (eval = FALSE)
###################################################
## toBibtex(citation("kebabs"))