BiSeq manual giving incompatible objects

Question

[dec986]

Hello,

I'm attempting to use the BiSeq package, and following the examples in their manual.

However, I am running into trouble with the compareTwoSamples function, which is supposed to take a BSrel object, but following the manual, this ends up being an S4 object.

Does anyone have any script that they have run BiSeq successfully with? I am completely stuck

library(methods)#base library, ensures correct execution from command line

b_files <- list.files(path = '.', pattern = "BL.*RRBS.*background.1.*.class.1.*")
f_files <- list.files(path = '.', pattern = "FL.*RRBS.*background.1.*.class.1.*")

all_files <- c(b_files, f_files)
rm(b_files, f_files)
library(BiSeq)#will need for DataFrame
groups <- c( rep('B',10), rep('F', 10) )
colData <- DataFrame(group = c( rep('B',10), rep('F', 10) ), row.names = all_files)
rrbs <- readBismark(all_files, colData)
rrbs.clust.unlim <- clusterSites(object = rrbs, groups = factor(colData(rrbs)$group), perc.samples = 1, min.sites = 20, max.dist = 100)
#colData(rrbs)$group must be run as a factor, which the manual doesn't say http://seqanswers.com/forums/archive/index.php/t-42504.html
clusterSitesToGR(rrbs.clust.unlim)
ind.cov <- totalReads(rrbs.clust.unlim) > 0
quant <- quantile(totalReads(rrbs.clust.unlim)[ind.cov], 0.9)
rrbs.clust.lim <- limitCov(rrbs.clust.unlim, maxCov = quant)
predictedMeth <- predictMeth(object = rrbs.clust.lim)#fairly slow (~2 min on RRBS data)
#b <- predictedMeth[, colData(predictedMeth)$group == 'B']
#f <- predictedMeth[, colData(predictedMeth)$group == 'F']
mc.cores <- 4
save.image()
print("about to do betaRegression for betaResults")
betaResults <- betaRegression(formula = ~group, link = "probit", object = predictedMeth, type = "BR")
data(betaResults)
save.image()
print("about to make betaResultsNull")
#betaResultsNull <- betaRegression(formula = ~group.null, link = "probit", object = predictedMethNull, type="BR")
#data(betaResultsNull)
save.image()
print('about to do makeVariogram')
vario <- makeVariogram(betaResults)
save.image()
data(vario)
print('about to do vario.sm')
save.image()
vario.sm <- smoothVariogram(vario, sill = 0.9)
print('about to make vario.aux')
save.image()
vario.aux <- makeVariogram(betaResults, make.variogram=FALSE)
print('about to make vario.sm$pValsList')
save.image()
vario.sm$pValsList <- vario.aux$pValsList
print('about to make locCor')
save.image()
locCor <- estLocCorvario.sm)
print('about to make clusters.rej')
save.image()
clusters.rej <- testClusters(locCor, FDR.cluster = 0.1)
print('about to make clusters.trimmed')
save.image()
clusters.trimmed <- trimClusters(clusters.rej, FDR.loc = 0.05)
print('about to make DMRs')
DMRs <- findDMRs(clusters.trimmed, max.dist = 100, diff.dir = TRUE)
save.image()