###################################################
### chunk number 1: 
###################################################
#line 51 "vignettes/Rolexa/inst/doc/Rolexa-vignette.Rnw"
library(Rolexa)
rolenv = SetModel(idsep="_")
GetModel(rolenv)


###################################################
### chunk number 2: build
###################################################
#line 86 "vignettes/Rolexa/inst/doc/Rolexa-vignette.Rnw"
path = SolexaPath(system.file("extdata", package="ShortRead"))


###################################################
### chunk number 3: load
###################################################
#line 90 "vignettes/Rolexa/inst/doc/Rolexa-vignette.Rnw"
(int = readIntensities(path,pattern="s_1_0001",withVariability=FALSE))
(seq = CombineReads(run=rolenv,path=path,pattern="s_1_0001_seq*"))
(seq_fastq = readFastq(path))


###################################################
### chunk number 4: cross-talk
###################################################
#line 102 "vignettes/Rolexa/inst/doc/Rolexa-vignette.Rnw"
(theta=OptimizeAngle(int=int))[1:10,]
int=DeCorrelateChannels(int=int,theta=theta)


###################################################
### chunk number 5: dephase
###################################################
#line 107 "vignettes/Rolexa/inst/doc/Rolexa-vignette.Rnw"
(rate=OptimizeRate(int=int))
int=DeCorrelateCycles(int=int,rate=rate)


###################################################
### chunk number 6: normalize
###################################################
#line 112 "vignettes/Rolexa/inst/doc/Rolexa-vignette.Rnw"
int2=TileNormalize(run=rolenv,int=int)


###################################################
### chunk number 7: evalscore
###################################################
#line 122 "vignettes/Rolexa/inst/doc/Rolexa-vignette.Rnw"
(res=SeqScore(run=rolenv,int=int,seqInit=seq,cycles=1:36))$sread


###################################################
### chunk number 8: filter
###################################################
#line 133 "vignettes/Rolexa/inst/doc/Rolexa-vignette.Rnw"
rolenv@MinimumTagLength = as.integer(1)
(res2 = FilterResults(run=rolenv,results=res))$sread
str = as.matrix(res$sread[241:253])
nt = DNA_ALPHABET
post.entropy = matrix(0,nrow=nrow(str),ncol=36)
post.entropy[which(str %in% nt[5:10])] = 1
post.entropy[which(str %in% nt[11:14])] = log2(3)
post.entropy[which(str == 'N')] = 2
matplot(1:36,y=apply(post.entropy,1,cumsum),t='l',lty=1,col=rainbow(6),ylim=c(0,30),xlim=c(1,36),xlab="cycles",ylab="cumulative entropy",main="Tag length cutoff")
lines(1:36,rolenv@IThresholds,t='l',lty=2,lwd=2,col="tomato")
abline(v=nchar(res2$sread[241:253]),col=rainbow(6),lty=2)
legend(x=0,y=30,res2$sread[241:253],col=rainbow(6),lty=1,bg="white",cex=.8)


###################################################
### chunk number 9: save eval=FALSE
###################################################
## #line 149 "vignettes/Rolexa/inst/doc/Rolexa-vignette.Rnw"
## SaveResults(run=rolenv,results=res2,outpath="./")


###################################################
### chunk number 10: fork eval=FALSE
###################################################
## #line 158 "vignettes/Rolexa/inst/doc/Rolexa-vignette.Rnw"
## library(fork)
## ForkBatch(run=rolenv,path=path,outpath="./",prefix="rs_",nthreads=2,nfiles=5,lane=1,tiles=1,idsep="_")


###################################################
### chunk number 11: onebatch eval=FALSE
###################################################
## #line 163 "vignettes/Rolexa/inst/doc/Rolexa-vignette.Rnw"
## OneBatch(run=rolenv,path=path,lane=1,tiles=tiles[n:m],outpath="./",prefix="rs_")


###################################################
### chunk number 12: combined
###################################################
#line 178 "vignettes/Rolexa/inst/doc/Rolexa-vignette.Rnw"
CombinedPlot(run=rolenv,int=int,seq=seq,scores=as(quality(seq_fastq),"matrix"),colonies=sample(1:nrow(int),4),par=list(mfrow=c(2,2),cex=.6,mar=c(4, 4, 2, 1)+.1))


###################################################
### chunk number 13: histo
###################################################
#line 185 "vignettes/Rolexa/inst/doc/Rolexa-vignette.Rnw"
ChannelHistogram(int=int,cycles=1,par=list(mfrow=c(2,2),mar=c(4, 4, 2, 1)+.1))


###################################################
### chunk number 14: plotcycles
###################################################
#line 188 "vignettes/Rolexa/inst/doc/Rolexa-vignette.Rnw"
par(mfrow=c(2,2),mar=c(4, 4, 2, 1)+.1)
PlotCycles(run=rolenv,int=int,seq=seq,cycles=c(1,20))


###################################################
### chunk number 17: sessionInfo
###################################################
#line 215 "vignettes/Rolexa/inst/doc/Rolexa-vignette.Rnw"
toLatex(sessionInfo())