### R code from vignette source '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:141-143 (eval = FALSE) ################################################### ## source("http://bioconductor.org/biocLite.R") ## biocLite("kebabs") ################################################### ### code chunk number 3: kebabs.Rnw:164-165 ################################################### library(kebabs) ################################################### ### code chunk number 4: kebabs.Rnw:172-173 (eval = FALSE) ################################################### ## vignette("kebabs") ################################################### ### code chunk number 5: kebabs.Rnw:178-179 (eval = FALSE) ################################################### ## help(kebabs) ################################################### ### code chunk number 6: kebabs.Rnw:186-187 ################################################### data(TFBS) ################################################### ### code chunk number 7: kebabs.Rnw:192-194 ################################################### enhancerFB length(enhancerFB) ################################################### ### code chunk number 8: kebabs.Rnw:199-201 (eval = FALSE) ################################################### ## hist(width(enhancerFB), breaks=30, xlab="Sequence Length", ## main="Distribution of Sequence Lengths") ################################################### ### code chunk number 9: kebabs.Rnw:204-208 ################################################### pdf("001.pdf") hist(width(enhancerFB), breaks=30, main="Distribution of Sequence Lenghts", xlab="Sequence Length") dev.off() ################################################### ### code chunk number 10: kebabs.Rnw:219-220 ################################################### showAnnotatedSeq(enhancerFB, sel=3) ################################################### ### code chunk number 11: kebabs.Rnw:225-226 ################################################### head(yFB) ################################################### ### code chunk number 12: kebabs.Rnw:231-232 ################################################### table(yFB) ################################################### ### code chunk number 13: kebabs.Rnw:239-243 ################################################### numSamples <- length(enhancerFB) trainingFraction <- 0.7 train <- sample(1:numSamples, trainingFraction * numSamples) test <- c(1:numSamples)[-train] ################################################### ### code chunk number 14: kebabs.Rnw:250-251 ################################################### specK2 <- spectrumKernel(k=2) ################################################### ### code chunk number 15: kebabs.Rnw:258-260 ################################################### model <- kbsvm(x=enhancerFB[train], y=yFB[train], kernel=specK2, pkg="e1071", svm="C-svc", cost=15) ################################################### ### code chunk number 16: kebabs.Rnw:267-270 ################################################### pred <- predict(model, enhancerFB[test]) head(pred) head(yFB[test]) ################################################### ### code chunk number 17: kebabs.Rnw:275-276 ################################################### evaluatePrediction(pred, yFB[test], allLabels=unique(yFB)) ################################################### ### code chunk number 18: kebabs.Rnw:283-289 ################################################### 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:294-300 ################################################### 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:307-313 ################################################### 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:318-319 (eval = FALSE) ################################################### ## ?kbsvm ################################################### ### code chunk number 22: kebabs.Rnw:417-418 ################################################### specK2 <- spectrumKernel(k=2) ################################################### ### code chunk number 23: kebabs.Rnw:423-424 ################################################### specK2 <- spectrumKernel(k=2, normalized=FALSE) ################################################### ### code chunk number 24: kebabs.Rnw:450-451 ################################################### mismK3M1 <- mismatchKernel(k=3, m=1) ################################################### ### code chunk number 25: kebabs.Rnw:467-468 ################################################### gappyK1M2 <- gappyPairKernel(k=1, m=2) ################################################### ### code chunk number 26: kebabs.Rnw:490-493 ################################################### 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:531-533 ################################################### gappyK1M2ps <- gappyPairKernel(k=1, m=2, distWeight=1, normalized=FALSE) ################################################### ### code chunk number 28: kebabs.Rnw:540-543 ################################################### seq1 <- AAStringSet(c("GACGAGGACCGA","AGTAGCGAGGT","ACGAGGTCTTT", "GGACCGAGTCGAGG")) positionMetadata(seq1) <- c(3, 5, 2, 10) ################################################### ### code chunk number 29: kebabs.Rnw:550-554 ################################################### 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:559-562 ################################################### positionMetadata(seq1) <- NULL km <- getKernelMatrix(wdK3, seq1) km ################################################### ### code chunk number 31: kebabs.Rnw:605-614 (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:617-628 ################################################### 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:645-654 ################################################### 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:659-664 ################################################### 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:678-689 ################################################### 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:722-786 (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:791-793 (eval = FALSE) ################################################### ## annotationCharset(hba) ## showAnnotatedSeq(hba, sel=1, start=1, end=400) ################################################### ### code chunk number 38: kebabs.Rnw:798-804 (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:809-813 (eval = FALSE) ################################################### ## km <- linearKernel(erK2) ## km ## kma <- linearKernel(erK2a) ## kma ################################################### ### code chunk number 40: kebabs.Rnw:822-837 ################################################### 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:842-852 ################################################### ## 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:857-866 ################################################### ## 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:871-875 ################################################### 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:880-881 (eval = FALSE) ################################################### ## plot(result, sel="BACC") ################################################### ### code chunk number 45: kebabs.Rnw:884-887 ################################################### pdf("005.pdf") plot(result, sel="BACC") dev.off() ################################################### ### code chunk number 46: kebabs.Rnw:905-950 ################################################### ## 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:959-967 ################################################### 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:972-974 ################################################### kernelList1 <- list(spectrumKernel(k=3), mismatchKernel(k=3,m=1), gappyPairKernel(k=2,m=4)) ################################################### ### code chunk number 49: kebabs.Rnw:979-980 ################################################### kernelList2 <- c(spectrumKernel(k=2:4), gappyPairKernel(k=1, m=2:5)) ################################################### ### code chunk number 50: kebabs.Rnw:999-1004 ################################################### specK2 <- spectrumKernel(k=2) km <- getKernelMatrix(specK2, x=enhancerFB) class(km) dim(km) km[1:3, 1:3] ################################################### ### code chunk number 51: kebabs.Rnw:1009-1010 (eval = FALSE) ################################################### ## heatmap(km, symm=TRUE) ################################################### ### code chunk number 52: kebabs.Rnw:1028-1031 ################################################### specK2 <- spectrumKernel(k=2) km <- specK2(x=enhancerFB) km[1:3, 1:3] ################################################### ### code chunk number 53: kebabs.Rnw:1036-1038 ################################################### km <- getKernelMatrix(specK2, x=enhancerFB, selx=c(1,4,25,137,300)) km ################################################### ### code chunk number 54: kebabs.Rnw:1043-1048 ################################################### 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:1053-1056 ################################################### km <- getKernelMatrix(specK2, x=enhancerFB, selx=1:200, y=enhancerFB, sely=201:500) dim(km) ################################################### ### code chunk number 56: kebabs.Rnw:1073-1076 ################################################### specK2 <- spectrumKernel(k=2, normalized=FALSE) erd <- getExRep(enhancerFB, selx=1:5, kernel=specK2, sparse=FALSE) erd ################################################### ### code chunk number 57: kebabs.Rnw:1081-1086 ################################################### 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:1091-1099 ################################################### 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:1113-1118 ################################################### 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:1123-1125 (eval = FALSE) ################################################### ## aggres <- aggExCluster(km, apres) ## plot(aggres) ################################################### ### code chunk number 61: kebabs.Rnw:1128-1132 ################################################### aggres <- aggExCluster(km, apres) pdf("004.pdf") plot(aggres) dev.off() ################################################### ### code chunk number 62: kebabs.Rnw:1143-1146 ################################################### exrep <- getExRep(enhancerFB, gappyK1M4, sparse=FALSE) apres1 <- apcluster(s=linearKernel, x=exrep, p=0.1) length(apres1) ################################################### ### code chunk number 63: kebabs.Rnw:1158-1163 ################################################### 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:1178-1188 ################################################### 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:1193-1202 ################################################### ## 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)) ################################################### ### code chunk number 66: kebabs.Rnw:1213-1215 ################################################### preddec <- predict(model, exrep[test, ], predictionType="decision") evaluatePrediction(pred, yFB[test], allLabels=unique(yFB), decValues=preddec) ################################################### ### code chunk number 67: kebabs.Rnw:1220-1223 ################################################### perf <- evaluatePrediction(pred, yFB[test], allLabels=unique(yFB), decValues=preddec, print=FALSE) perf ################################################### ### code chunk number 68: kebabs.Rnw:1228-1230 (eval = FALSE) ################################################### ## rocdata <- computeROCandAUC(preddec, yFB[test], unique(yFB)) ## plot(rocdata, main="Receiver Operating Characteristics") ################################################### ### code chunk number 69: kebabs.Rnw:1233-1237 ################################################### rocdata <- computeROCandAUC(preddec, yFB[test], unique(yFB)) pdf("008.pdf") plot(rocdata, main="Receiver Operating Characteristics") dev.off() ################################################### ### code chunk number 70: kebabs.Rnw:1291-1301 ################################################### 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 71: kebabs.Rnw:1309-1317 ################################################### 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 72: kebabs.Rnw:1326-1337 ################################################### 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 73: kebabs.Rnw:1349-1358 ################################################### 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 74: kebabs.Rnw:1364-1365 ################################################### performance(cvResult(model))$foldErrors ################################################### ### code chunk number 75: kebabs.Rnw:1375-1381 ################################################### 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 76: kebabs.Rnw:1386-1398 ################################################### 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 77: kebabs.Rnw:1403-1409 ################################################### 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 78: kebabs.Rnw:1414-1425 ################################################### 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 79: kebabs.Rnw:1491-1508 ################################################### 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 80: kebabs.Rnw:1513-1514 (eval = FALSE) ################################################### ## plot(predProf, sel=1, ylim=c(-0.4, 0.2), heptads=TRUE, annotate=TRUE) ################################################### ### code chunk number 81: kebabs.Rnw:1517-1520 ################################################### pdf("006.pdf") plot(predProf, sel=1, ylim=c(-0.4, 0.2), heptads=TRUE, annotate=TRUE) dev.off() ################################################### ### code chunk number 82: kebabs.Rnw:1531-1532 (eval = FALSE) ################################################### ## plot(predProf, sel=c(1,2), ylim=c(-0.4, 0.2), heptads=TRUE, annotate=TRUE) ################################################### ### code chunk number 83: kebabs.Rnw:1535-1538 ################################################### pdf("007.pdf") plot(predProf, sel=c(1,2), ylim=c(-0.4, 0.2), heptads=TRUE, annotate=TRUE) dev.off() ################################################### ### code chunk number 84: kebabs.Rnw:1559-1566 ################################################### 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 85: kebabs.Rnw:1571-1576 ################################################### featWeights <- featureWeights(model) length(featWeights) featWeights[[1]][1:5] featWeights[[2]][1:5] featWeights[[3]][1:5] ################################################### ### code chunk number 86: kebabs.Rnw:1581-1585 ################################################### predProf <- getPredictionProfile(enhancerFB, gappyK1M2, featureWeights(model)[[2]], modelOffset(model)[2]) predProf ################################################### ### code chunk number 87: kebabs.Rnw:1825-1826 (eval = FALSE) ################################################### ## toBibtex(citation("kebabs"))