## ----Init,echo=FALSE,message=FALSE,results='hide'------------------- options(width=70) options(useFancyQuotes=FALSE) knitr::opts_knit$set(width=70) 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="") ## ----eval=F--------------------------------------------------------- # if (!requireNamespace("BiocManager", quietly=TRUE)) # install.packages("BiocManager") # BiocManager::install("kebabs") ## ------------------------------------------------------------------- library(kebabs) ## ----eval=FALSE----------------------------------------------------- # vignette("kebabs") ## ----eval=FALSE----------------------------------------------------- # help(kebabs) ## ------------------------------------------------------------------- data(TFBS) ## ------------------------------------------------------------------- enhancerFB length(enhancerFB) ## ----eval=FALSE----------------------------------------------------- # hist(width(enhancerFB), breaks=30, xlab="Sequence Length", # main="Distribution of Sequence Lengths", col="lightblue") ## ----echo=FALSE,results='hide'-------------------------------------- pdf("001.pdf", width=8, height=5.5) hist(width(enhancerFB), breaks=30, xlab="Sequence Length", main="Distribution of Sequence Lengths", col="lightblue") dev.off() ## ------------------------------------------------------------------- showAnnotatedSeq(enhancerFB, sel=3) ## ------------------------------------------------------------------- head(yFB) ## ------------------------------------------------------------------- table(yFB) ## ------------------------------------------------------------------- numSamples <- length(enhancerFB) trainingFraction <- 0.7 train <- sample(1:numSamples, trainingFraction * numSamples) test <- c(1:numSamples)[-train] ## ------------------------------------------------------------------- specK2 <- spectrumKernel(k=2) ## ------------------------------------------------------------------- model <- kbsvm(x=enhancerFB[train], y=yFB[train], kernel=specK2, pkg="e1071", svm="C-svc", cost=15) ## ------------------------------------------------------------------- pred <- predict(model, enhancerFB[test]) head(pred) head(yFB[test]) ## ------------------------------------------------------------------- evaluatePrediction(pred, yFB[test], allLabels=unique(yFB)) ## ------------------------------------------------------------------- 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)) ## ------------------------------------------------------------------- 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)) ## ------------------------------------------------------------------- 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)) ## ----eval=FALSE----------------------------------------------------- # ?kbsvm ## ------------------------------------------------------------------- specK2 <- spectrumKernel(k=2) ## ------------------------------------------------------------------- specK2 <- spectrumKernel(k=2, normalized=FALSE) ## ------------------------------------------------------------------- mismK3M1 <- mismatchKernel(k=3, m=1) ## ------------------------------------------------------------------- gappyK1M2 <- gappyPairKernel(k=1, m=2) ## ------------------------------------------------------------------- motifCollection1 <- c("A[CG]T","C.G","C..G.T","G[^A][AT]", "GT.A[CA].[CT]G") motif1 <- motifKernel(motifCollection1) ## ------------------------------------------------------------------- gappyK1M2ps <- gappyPairKernel(k=1, m=2, distWeight=1, normalized=FALSE) ## ------------------------------------------------------------------- seq1 <- AAStringSet(c("GACGAGGACCGA","AGTAGCGAGGT","ACGAGGTCTTT", "GGACCGAGTCGAGG")) positionMetadata(seq1) <- c(3, 5, 2, 10) ## ------------------------------------------------------------------- wdK3 <- spectrumKernel(k=3, distWeight=1, mixCoef=c(0.5,0.33,0.17), normalized=FALSE) km <- getKernelMatrix(wdK3, seq1) km ## ------------------------------------------------------------------- positionMetadata(seq1) <- NULL km <- getKernelMatrix(wdK3, seq1) km ## ----eval=F--------------------------------------------------------- # curve(linWeight(x, sigma=5), from=-25, to=25, xlab="p - q", ylab="weight", # main="Predefined Distance Weighting Functions", col="green", lwd=2) # curve(expWeight(x, sigma=5), from=-25, to=25, col="blue", lwd=2, add=TRUE) # curve(gaussWeight(x, sigma=5), from=-25, to=25, col="red", lwd=2, add=TRUE) # curve(swdWeight(x), from=-25, to=25, col="orange", lwd=2, add=TRUE) # legend("topright", inset=0.03, title="Weighting Functions", c("linWeight", # "expWeight", "gaussWeight", "swdWeight"), col=c("green", "blue", # "red", "orange"), lwd=2) # text(19, 0.55, expression(paste(sigma, " = 5"))) ## ----echo=FALSE,results='hide'-------------------------------------- pdf("002.pdf", width=8, height=5.5) curve(linWeight(x, sigma=5), from=-25, to=25, xlab="p - q", ylab="weight", main="Predefined Distance Weighting Functions", col="green", lwd=2) curve(expWeight(x, sigma=5), from=-25, to=25, col="blue", lwd=2, add=TRUE) curve(gaussWeight(x, sigma=5), from=-25, to=25, col="red", lwd=2, n=201, add=TRUE) curve(swdWeight(x), from=-25, to=25, col="orange", lwd=2, add=TRUE) legend("topright", inset=0.03, title="Weighting Functions", c("linWeight", "expWeight", "gaussWeight", "swdWeight"), col=c("green", "blue", "red", "orange"), lwd=2) text(19, 0.55, expression(paste(sigma, " = 5"))) dev.off() ## ------------------------------------------------------------------- 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)) ## ------------------------------------------------------------------- data(TFBS) names(enhancerFB) <- paste("Sample", 1:length(enhancerFB), sep="_") enhancerFB kmSWD <- getKernelMatrix(swdK3, x=enhancerFB, selx=1:5) kmSWD[1:5, 1:5] ## ------------------------------------------------------------------- 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) ## ----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") ## ----eval=FALSE----------------------------------------------------- # annotationCharset(hba) # showAnnotatedSeq(hba, sel=1, start=1, end=400) ## ----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] ## ----eval=FALSE----------------------------------------------------- # km <- linearKernel(erK2) # km # kma <- linearKernel(erK2a) # kma ## ------------------------------------------------------------------- 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)) ## ------------------------------------------------------------------- ## 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)) ## ------------------------------------------------------------------- ## 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 ## ------------------------------------------------------------------- perfData <- performance(result) perfData which(perfData$BACC[1, ] == max(perfData$BACC[1, ])) which(perfData$BACC[2, ] == max(perfData$BACC[2, ])) ## ----eval=FALSE----------------------------------------------------- # plot(result, sel="BACC") ## ----echo=FALSE,results='hide'-------------------------------------- pdf("005.pdf", width=7, height=5) plot(result, sel="BACC") dev.off() ## ------------------------------------------------------------------- ## 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)) ## ------------------------------------------------------------------- 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) ## ------------------------------------------------------------------- kernelList1 <- list(spectrumKernel(k=3), mismatchKernel(k=3, m=1), gappyPairKernel(k=2, m=4)) ## ------------------------------------------------------------------- kernelList2 <- c(spectrumKernel(k=2:4), gappyPairKernel(k=1, m=2:5)) ## ------------------------------------------------------------------- specK2 <- spectrumKernel(k=2) km <- getKernelMatrix(specK2, x=enhancerFB) class(km) dim(km) km[1:3, 1:3] ## ----eval=FALSE----------------------------------------------------- # heatmap(km, symm=TRUE) ## ------------------------------------------------------------------- specK2 <- spectrumKernel(k=2) km <- specK2(x=enhancerFB) km[1:3, 1:3] ## ------------------------------------------------------------------- km <- getKernelMatrix(specK2, x=enhancerFB, selx=c(1, 4, 25, 137, 300)) km ## ------------------------------------------------------------------- seqs1 <- enhancerFB[1:200] seqs2 <- enhancerFB[201:500] km <- getKernelMatrix(specK2, x=seqs1, y=seqs2) dim(km) km[1:4, 1:5] ## ------------------------------------------------------------------- km <- getKernelMatrix(specK2, x=enhancerFB, selx=1:200, y=enhancerFB, sely=201:500) dim(km) ## ------------------------------------------------------------------- specK2 <- spectrumKernel(k=2, normalized=FALSE) erd <- getExRep(enhancerFB, selx=1:5, kernel=specK2, sparse=FALSE) erd ## ------------------------------------------------------------------- specK6 <- spectrumKernel(k=6, normalized=FALSE) erd <- getExRep(enhancerFB, selx=1:5, kernel=specK6, sparse=FALSE) dim(erd) erd[, 1:6] ## ------------------------------------------------------------------- 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] ## ------------------------------------------------------------------- library(apcluster) gappyK1M4 <- gappyPairKernel(k=1, m=4) km <- getKernelMatrix(gappyK1M4, enhancerFB) apres <- apcluster(s=km, p=0.8) length(apres) ## ----eval=FALSE----------------------------------------------------- # aggres <- aggExCluster(km, apres) # plot(aggres) ## ----echo=FALSE,results='hide'-------------------------------------- aggres <- aggExCluster(km, apres) pdf("004.pdf") plot(aggres) dev.off() ## ------------------------------------------------------------------- exrep <- getExRep(enhancerFB, gappyK1M4, sparse=FALSE) apres1 <- apcluster(s=linearKernel, x=exrep, p=0.1) length(apres1) ## ------------------------------------------------------------------- exrep <- getExRep(x=enhancerFB, selx=1:5, gappyK1M4, sparse=FALSE) dim(exrep) erquad <- getExRepQuadratic(exrep) dim(erquad) erquad[1:5, 1:5] ## ------------------------------------------------------------------- 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)) ## ------------------------------------------------------------------- ## compute symmetric kernel matrix for training samples kmtrain <- getKernelMatrix(gappyK1M4, x=enhancerFB, selx=train) model1 <- kbsvm(x=kmtrain, y=yFB[train], kernel=gappyK1M4, pkg="e1071", 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) pred1 <- predict(model1, kmtest) evaluatePrediction(pred1, yFB[test], allLabels=unique(yFB)) ## ------------------------------------------------------------------- preddec <- predict(model, exrep[test, ], predictionType="decision") evaluatePrediction(pred, yFB[test], allLabels=unique(yFB), decValues=preddec) ## ------------------------------------------------------------------- perf <- evaluatePrediction(pred, yFB[test], allLabels=unique(yFB), decValues=preddec, print=FALSE) perf ## ----eval=FALSE----------------------------------------------------- # rocdata <- computeROCandAUC(preddec, yFB[test], unique(yFB)) # plot(rocdata, main="Receiver Operating Characteristics", col="red", lwd=2) ## ----echo=FALSE,results='hide'-------------------------------------- rocdata <- computeROCandAUC(preddec, yFB[test], unique(yFB)) pdf("008.pdf") plot(rocdata, main="Receiver Operating Characteristics", col="red", lwd=2) dev.off() ## ------------------------------------------------------------------- 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) ## ------------------------------------------------------------------- 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) ## ------------------------------------------------------------------- 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) ## ------------------------------------------------------------------- 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) ## ------------------------------------------------------------------- performance(cvResult(model))$foldErrors ## ------------------------------------------------------------------- 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) ## ------------------------------------------------------------------- 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)] ## ------------------------------------------------------------------- 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)] ## ------------------------------------------------------------------- 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) ## ------------------------------------------------------------------- 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 ## ----eval=FALSE----------------------------------------------------- # plot(predProf, sel=1, ylim=c(-0.4, 0.2), heptads=TRUE, annotate=TRUE) ## ----echo=FALSE,results='hide'-------------------------------------- pdf("006.pdf", width=7, height=5.5) plot(predProf, sel=1, ylim=c(-0.4, 0.2), heptads=TRUE, annotate=TRUE) dev.off() ## ----eval=FALSE----------------------------------------------------- # plot(predProf, sel=c(1, 2), ylim=c(-0.4, 0.2), heptads=TRUE, annotate=TRUE) ## ----echo=FALSE,results='hide'-------------------------------------- pdf("007.pdf", width=7, height=5.5) plot(predProf, sel=c(1, 2), ylim=c(-0.4, 0.2), heptads=TRUE, annotate=TRUE) dev.off() ## ------------------------------------------------------------------- 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)) ## ------------------------------------------------------------------- featWeights <- featureWeights(model) length(featWeights) featWeights[[1]][1:5] featWeights[[2]][1:5] featWeights[[3]][1:5] ## ------------------------------------------------------------------- predProf <- getPredictionProfile(enhancerFB, gappyK1M2, featureWeights(model)[[2]], modelOffset(model)[2]) predProf ## ----eval=FALSE----------------------------------------------------- # toBibtex(citation("kebabs"))