Simulating uncertainty about randomly sampling a fraction of a small number of mtDNA molecules.

simSamp.R

# Sampling noise

simulateSample = function(mtDNAPerWell=1000, mutantFraction=0.95, fracVolumeSampled=1.0/3.0){
 # Create a synthetic population of mtDNA molecules
 moleculeMutant = rbinom(n = round(mtDNAPerWell*fracVolumeSampled), size = 1, prob = mutantFraction)
 # Calculate mutant fraction in sampled molecules
 sum(moleculeMutant)/length(moleculeMutant)
}

simulateUncertainty = function(mtDNAPerWell=1000, mutantFraction=0.95, fracVolumeSampled=1.0/3.0, reps = 5000){
  experiment = replicate(reps,simulateSample(mtDNAPerWell=mtDNAPerWell, mutantFraction=mutantFraction, fracVolumeSampled=fracVolumeSampled))
  mlab = paste("mtDNAperWell:",mtDNAPerWell,"mutantFraction:",mutantFraction,"\nfracVolumeSampled:",formatC(fracVolumeSampled,3))
  hist(experiment,xlim=c(0,1),breaks=seq(0,1,0.01),xlab="Observed Fraction",main=mlab,cex.lab=2.55)
  abline(v=mutantFraction,col="red",lwd=3)
  legend("topleft",c("True Fraction"),col="red",lwd=2)
}

pdf("simulatingSampling.pdf",width=14,height=7)
 op=par(mfrow=c(1,2),mar = c(5,5,3,1))
 for(mtdna in c(10,20,50,100,200,500,1000)){
    simulateUncertainty(mtDNAPerWell=mtdna, mutantFraction=0.5, fracVolumeSampled=1.0/3.0, reps = 50000)
    simulateUncertainty(mtDNAPerWell=mtdna, mutantFraction=0.5, fracVolumeSampled=1.0, reps = 50000)
 }
 par(op)
dev.off()