make_error_models.R

In alyssafrazee/polyester: Simulate RNA-seq reads

## create error models to be released wtih polyester
## from the GemSIM/GemErr estimation model

modfolder = '../polyester_paper/error_models'
platforms = c('ill100v4_mate1', 'ill100v4_mate2', 
    'ill100v4_single', 'ill100v5_mate1', 'ill100v5_mate2', 
    'ill100v5_single', 'r454ti_single')
for(platform in platforms){
    i = which(platforms == platform)
    model = read.table(paste(modfolder, platform, sep='/'), header=TRUE)
    names(model)[1] = 'refbase'
    model$refbase = substr(model$refbase, 4, 4)
    eval(parse(text=paste0('model', i, ' <- model')))
    save(list=paste0('model',i), file=paste0('data/', platform, '.rda'), 
        compress='xz')
}

## figures for paper:
## put one example in true manuscript,
## include the others in supplementary data.
library(RSkittleBrewer)
library(usefulstuff)
colrs = RSkittleBrewer('tropical')
getColor = function(nt){
    switch(nt, A='black', C=colrs[1], G=colrs[2], T='deeppink', N=colrs[4])
}
nts = c('A','C','G','T','N')

plot_nt = function(model, nt){
    d = subset(model, refbase==nt)
    wrongnts = nts[-which(nts==nt)]
    errCols = paste0('read', wrongnts)
    mnum = as.matrix(model[,2:6])
    ymax = max(mnum[mnum<0.8])
    plot(1:100, as.matrix(d[errCols[1]])[-1], col=makeTransparent(getColor(wrongnts[1])),
        pch=19, cex=0.5, xlab='Read Position', ylab='Error Probability', 
        type='o', ylim=c(0,ymax))
    for(i in 2:4){
        points(1:100, as.matrix(d[errCols[i]])[-1], pch=19, cex=0.5, type='o',
            col=makeTransparent(getColor(wrongnts[i])))
    }
    title(nt)
    legend('topleft', wrongnts, pch=19, col=sapply(wrongnts, getColor), cex=0.7)
}

plot_model = function(model, file){
    pdf(file)
    par(mfrow=c(2,2))
    for(nt in c('A','C','G','T')){
        plot_nt(model, nt)
    }
    dev.off()
}

plot_model(model1, 'illumina4_mate1.pdf')
plot_model(model2, 'illumina4_mate2.pdf')
plot_model(model3, 'illumina4_single.pdf')
plot_model(model4, 'illumina5_mate1.pdf')
plot_model(model5, 'illumina5_mate2.pdf')
plot_model(model6, 'illumina5_single.pdf')