Extends Wang, et al, 2007 Fig 1 (pmid 17600086) to show detection threshold at higher coverage

wang_fig1_extended.R

## Extends Wang, et al, 2007 Fig 1 (pmid 17600086) to show detection threshold at higher coverage
##
                          
library(lattice)

## mu = error rate
mu_hp <- 0.0044
mu_nhp <- 0.0007

thresh <- function(N, mu, p=0.001){
  x <- 0:round(sqrt(N)) ## likely to fail for large values of N
  100*min(which(1-cumsum(dpois(x, N*mu)) < p))/N
}

pdf('wang_fig1_extended.pdf')
coverage <- seq(50, 2000, 50)
xyplot(threshold~coverage,
       data=do.call(rbind,
           lapply(c(mu_hp, mu_nhp),
                  function(mu){
                    data.frame(coverage=coverage,
                               threshold=sapply(coverage, thresh, mu),
                               mu=rep(gettextf('%f',mu),length(coverage)))
                  })),
       groups=mu,
       panel=function(x,y,...){
         panel.abline(h=seq(0.5,8,0.5),col='grey')
         panel.xyplot(x,y,...)
       },
       type='b',
       auto.key=TRUE,
       ylab='Detection threshold (%)',
       sub='reproduces Wang et al 2007 (pmid 17600086), Fig 1'       
       )

coverage <- seq(100, 20000, 100)
xyplot(threshold~coverage,
       data=do.call(rbind,
           lapply(c(mu_hp, mu_nhp, 0.0001),
                  function(mu){
                    data.frame(coverage=coverage,
                               threshold=sapply(coverage, thresh, mu),
                               mu=rep(gettextf('%f',mu),length(coverage)))
                  })),
       groups=mu,
       panel=function(x,y,...){
         panel.xyplot(x,y,...)
         panel.abline(v=2000)
         panel.abline(h=seq(0.25/2,1.5,0.25/2),col='grey')
       },
       type='l',
       ylab='Detection threshold (%)',       
       lwd=2,
       auto.key=TRUE,
       scales=list(
           x=list(log=FALSE),
           y=list(limits=c(0,1.5))
           ),
       sub='extends Wang et al 2007 (pmid 17600086), Fig 1'
       )

dev.off()