### R code from vignette source 'vignettes/Rolexa/inst/doc/Rolexa-vignette.Rnw' ################################################### ### code chunk number 1: Rolexa-vignette.Rnw:50-53 ################################################### library(Rolexa) rolenv = SetModel(idsep="_") GetModel(rolenv) ################################################### ### code chunk number 2: build ################################################### path = SolexaPath(system.file("extdata", package="ShortRead")) ################################################### ### code chunk number 3: load ################################################### (int = readIntensities(path,pattern="s_1_0001",withVariability=FALSE)) (seq = CombineReads(run=rolenv,path=path,pattern="s_1_0001_seq*")) (seq_fastq = readFastq(path)) ################################################### ### code chunk number 4: cross-talk ################################################### (theta=OptimizeAngle(int=int))[1:10,] int=DeCorrelateChannels(int=int,theta=theta) ################################################### ### code chunk number 5: dephase ################################################### (rate=OptimizeRate(int=int)) int=DeCorrelateCycles(int=int,rate=rate) ################################################### ### code chunk number 6: normalize ################################################### int2=TileNormalize(run=rolenv,int=int) ################################################### ### code chunk number 7: evalscore ################################################### (res=SeqScore(run=rolenv,int=int,seqInit=seq,cycles=1:36))$sread ################################################### ### code chunk number 8: filter ################################################### 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) ################################################### ### code chunk number 9: save (eval = FALSE) ################################################### ## SaveResults(run=rolenv,results=res2,outpath="./") ################################################### ### code chunk number 10: fork (eval = FALSE) ################################################### ## library(fork) ## ForkBatch(run=rolenv,path=path,outpath="./",prefix="rs_",nthreads=2,nfiles=5,lane=1,tiles=1,idsep="_") ################################################### ### code chunk number 11: onebatch (eval = FALSE) ################################################### ## OneBatch(run=rolenv,path=path,lane=1,tiles=tiles[n:m],outpath="./",prefix="rs_") ################################################### ### code chunk number 12: combined ################################################### 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)) ################################################### ### code chunk number 13: histo ################################################### ChannelHistogram(int=int,cycles=1,par=list(mfrow=c(2,2),mar=c(4, 4, 2, 1)+.1)) ################################################### ### code chunk number 14: plotcycles ################################################### par(mfrow=c(2,2),mar=c(4, 4, 2, 1)+.1) PlotCycles(run=rolenv,int=int,seq=seq,cycles=c(1,20)) ################################################### ### code chunk number 17: sessionInfo ################################################### toLatex(sessionInfo())