## ----installme, eval=FALSE-----------------------------------------------
## source("http://bioconductor.org/biocLite.R")
## biocLite("polyester")

## ----builtinex, warning=FALSE, message=FALSE-----------------------------
library(polyester)
library(Biostrings)

fasta_file = system.file('extdata', 'chr22.fa', package='polyester')
fasta = readDNAStringSet(fasta_file)
small_fasta = fasta[1:20]
writeXStringSet(small_fasta, 'chr22_small.fa')
fold_changes = c(4, 4, 1/4, 1/4, rep(1, 16))
outdir = 'simulated_reads'

# ~20x coverage ----> reads per transcript = length/readlength * 20
# "width" is operating on a DNAStringSet (from Biostrings)
readspertx = round(20 * width(small_fasta) / 100)
simulate_experiment('chr22_small.fa', reads_per_transcript=readspertx, 
    num_reps=10, fold_changes=fold_changes, outdir=outdir) 

## ----countmat------------------------------------------------------------
# set up matrix:
num_timepoints = 12
countmat = matrix(readspertx, nrow=length(small_fasta), ncol=num_timepoints)

# add spikes in expression at certain timepoints to certain transcripts:
up_early = c(1,2) 
up_late = c(3,4)
countmat[up_early, 2] = 3*countmat[up_early, 2]
countmat[up_early, 3] = round(1.5*countmat[up_early, 3])
countmat[up_late, 10] = 6*countmat[up_late, 10]
countmat[up_late, 11] = round(1.2*countmat[up_late, 11])

# simulate reads:
simulate_experiment_countmat('chr22_small.fa', readmat=countmat, 
    outdir='timecourse_reads') 

## ----installbg, eval=FALSE-----------------------------------------------
## source("http://bioconductor.org/biocLite.R")
## biocLite("ballgown")

## ----loadbg, warning=FALSE, message=FALSE--------------------------------
library(ballgown)

## ----datasim-------------------------------------------------------------
data(bg)
countmat = fpkm_to_counts(bg, threshold=0.01, mean_rps=400000)
params = get_params(countmat)
Ntranscripts = 50
Nsamples = 10
custom_readmat = create_read_numbers(mu=params$mu, fit=params$fit,
    p0=params$p0, m=Ntranscripts, n=Nsamples, seed=103)

## ----info, results='markup'----------------------------------------------
sessionInfo()