###################################################
### chunk number 1: 
###################################################
#line 52 "vignettes/seqbias/inst/doc/overview.Rnw"
library(seqbias)
library(Rsamtools)
ref_fn <- system.file( "extra/example.fa", package = "seqbias" )
ref_f <- FaFile( ref_fn )
open.FaFile( ref_f )


###################################################
### chunk number 2: 
###################################################
#line 61 "vignettes/seqbias/inst/doc/overview.Rnw"
reads_fn <- system.file( "extra/example.bam", package = "seqbias" )


###################################################
### chunk number 3: 
###################################################
#line 94 "vignettes/seqbias/inst/doc/overview.Rnw"
ref_seqs <- scanFaIndex( ref_f )


###################################################
### chunk number 4: 
###################################################
#line 103 "vignettes/seqbias/inst/doc/overview.Rnw"
I <- random.intervals( ref_seqs, n = 5, m = 100000 )


###################################################
### chunk number 5: 
###################################################
#line 112 "vignettes/seqbias/inst/doc/overview.Rnw"
seqs <- scanFa( ref_f, I )


###################################################
### chunk number 6: 
###################################################
#line 119 "vignettes/seqbias/inst/doc/overview.Rnw"
neg_idx <- as.logical( I@strand == '-' )
seqs[neg_idx] <- reverseComplement( seqs[neg_idx] )


###################################################
### chunk number 7: 
###################################################
#line 129 "vignettes/seqbias/inst/doc/overview.Rnw"
counts <- count.reads( reads_fn, I )


###################################################
### chunk number 8: 
###################################################
#line 142 "vignettes/seqbias/inst/doc/overview.Rnw"
freqs <- kmer.freq( seqs, counts )


###################################################
### chunk number 9: fig1
###################################################
#line 148 "vignettes/seqbias/inst/doc/overview.Rnw"
if( require(ggplot2) ) {
  P <- qplot( x = pos,
              y = freq,
              ylim = c(0.15,0.4),
              color = seq,
              data  = freqs,
              geom  = "line" )
  P <- P + facet_grid( seq ~ . )
  print(P)
} else {
  par( mar = c(5,1,1,1), mfrow = c(4,1) )
  with( subset( freqs, seq == "a" ),
        plot( freq ~ pos, ylim = c(0.15,0.4), sub = "a", type = 'l' ) )
  with( subset( freqs, seq == "c" ),
        plot( freq ~ pos, ylim = c(0.15,0.4), sub = "c", type = 'l' ) )
  with( subset( freqs, seq == "g" ),
        plot( freq ~ pos, ylim = c(0.15,0.4), sub = "g", type = 'l' ) )
  with( subset( freqs, seq == "t" ),
        plot( freq ~ pos, ylim = c(0.15,0.4), sub = "t", type = 'l' ) )
} 


###################################################
### chunk number 10: 
###################################################
#line 195 "vignettes/seqbias/inst/doc/overview.Rnw"
sb <- seqbias.fit( ref_fn, reads_fn, L = 5, R = 20 )


###################################################
### chunk number 11: 
###################################################
#line 210 "vignettes/seqbias/inst/doc/overview.Rnw"
bias <- seqbias.predict( sb, I )


###################################################
### chunk number 12: 
###################################################
#line 218 "vignettes/seqbias/inst/doc/overview.Rnw"
counts.adj <- mapply( FUN = `/`, counts, bias, SIMPLIFY = F )


###################################################
### chunk number 13: 
###################################################
#line 223 "vignettes/seqbias/inst/doc/overview.Rnw"
freqs.adj <- kmer.freq( seqs, counts.adj )


###################################################
### chunk number 14: fig2
###################################################
#line 228 "vignettes/seqbias/inst/doc/overview.Rnw"
if( require(ggplot2) ) {
  P <- qplot( x = pos,
              y = freq,
              ylim = c(0.15,0.4),
              color = seq,
              data  = freqs.adj,
              geom  = "line" )
  P <- P + facet_grid( seq ~ . )
  print(P)
} else {
  par( mar = c(5,1,1,1), mfrow = c(4,1) )
  with( subset( freqs.adj, seq == "a" ),
        plot( freq ~ pos, ylim = c(0.15,0.4), sub = "a", type = 'l' ) )
  with( subset( freqs.adj, seq == "c" ),
        plot( freq ~ pos, ylim = c(0.15,0.4), sub = "c", type = 'l' ) )
  with( subset( freqs.adj, seq == "g" ),
        plot( freq ~ pos, ylim = c(0.15,0.4), sub = "g", type = 'l' ) )
  with( subset( freqs.adj, seq == "t" ),
        plot( freq ~ pos, ylim = c(0.15,0.4), sub = "t", type = 'l' ) )
}


###################################################
### chunk number 15: 
###################################################
#line 266 "vignettes/seqbias/inst/doc/overview.Rnw"
seqbias.save( sb, "my_seqbias_model.yml" )

# load the model sometime later
sb <- seqbias.load( ref_fn, "my_seqbias_model.yml" )


###################################################
### chunk number 16: 
###################################################
#line 279 "vignettes/seqbias/inst/doc/overview.Rnw"
sessionInfo()