inst/doc/findsegments.R

In tilingArray: Transcript mapping with high-density oligonucleotide tiling arrays

### R code from vignette source 'findsegments.Rnw'

###################################################
### code chunk number 1: defGenData
###################################################
genData = function(lenx, nrSeg, nrRep=1, stddev=0.1) {
  x  = matrix(as.numeric(NA), nrow=lenx, ncol=nrRep)
  cp = sort(sample(1:floor(lenx/15), nrSeg-1) * 15)
  cpb = c(1, cp, lenx+1)
  s  = 0
  for (j in 2:length(cpb)) {
    sel = cpb[j-1]:(cpb[j]-1)
    s = (.5+runif(1))*sign(rnorm(1))+s
    x[sel, ] = rnorm(length(sel)*nrRep, mean=s, sd=stddev)
  }
  return(list(x=x, cp=cp))
}


###################################################
### code chunk number 2: plotData
###################################################
set.seed(4711)
lenx = 1000
nrSeg = 10
gd = genData(lenx, nrSeg)
plot(gd$x, pch=".")
abline(v=gd$cp, col="blue")


###################################################
### code chunk number 3: loadlib
###################################################
library("tilingArray")


###################################################
### code chunk number 4: findSegments
###################################################
maxseg = 12
maxk = 500
seg  = segment(gd$x, maxk=maxk, maxseg=maxseg)
seg
seg@breakpoints[nrSeg+(-1:1)]
gd$cp


###################################################
### code chunk number 5: checkClaim
###################################################
stopifnot(all(gd$cp==seg@breakpoints[[nrSeg]]))


###################################################
### code chunk number 6: testCP
###################################################
for(i in 1:20){
  gd  = genData(lenx, nrSeg)
  seg = segment(gd$x, maxk=maxk, maxseg=maxseg)
  stopifnot(seg@breakpoints[[nrSeg]][, "estimate"] == gd$cp)
}


###################################################
### code chunk number 7: modelSelect1
###################################################
nrSeg = 22
gd = genData(lenx, nrSeg, nrRep=2, stddev=1/3)
s = segment(gd$x, maxk=lenx, maxseg=as.integer(nrSeg * 2.5))


###################################################
### code chunk number 8: plotSimSeg
###################################################
par(mai=c(1,1,0.1,0.01))
plot(row(gd$x), gd$x, pch=".")


###################################################
### code chunk number 9: penLLsim
###################################################
par(mai=c(1,1,0.1,0.01))
plotPenLL(s)


###################################################
### code chunk number 10: printMaxima
###################################################
which.max(logLik(s, penalty="AIC"))
which.max(logLik(s, penalty="BIC"))