## Quantum Dot
#runs and reps are number of restarts (without and with removing DNA)
#kb is the number of bases run for
# par is the decay parameter
readQDot <- function(runs=5,kb=1,par=0.6) sum(replicate(runs,min(rexp(1,par),kb)))
readQDotQual <- function(reps=3,runs=5,kb=1,par=0.6) min(replicate(reps,readQDot()))
hQDot <- hist(replicate(100000,readQDot()),breaks=100)
hQDotQual <- hist(replicate(100000,readQDotQual()),breaks=100)
plot(hQDot$mids,hQDot$counts/sum(hQDot$counts),type="h",xlab="Read length",ylab="Prob",main="QDot Read Distribution")
lines(hQDotQual$mids,hQDotQual$counts/sum(3*hQDot$counts),type="h",col="red")
legend(min(hQDotQual$mids),max(hQDot$counts)/sum(hQDot$counts),c("All reads", "3-pass reads"),col=c("black","red"),pch=19)
## general Exponential 
readsExp <- function(N=100000,par=0.6) rexp(N,par)
plot((0:100)/10,dexp((0:100)/10,0.6),type="h",xlab="Read length",ylim=c(0,dexp(0,0.6)),ylab="Prob",main="PacBio or Nanopore Read Distribution",col="black")
lines((0:100)/10,10*dexp((0:100)/10,0.01),col="red")
legend(6,0.6,c("Expected","Hope Against Hope"),col=c("black","red"),pch=19)
## PacBio Multipass
# kb is the size of the SMRTbell
readsPacBioQual <- function(kb=0.2,N=100000,par=0.6) { temp <- readsExp(N,par)/kb; sapply(1:5,function(x) mean(temp > x))}
plot(1:5,readsPacBioQual(0.3),type="l",ylim=c(0,1),xlab="Passes",ylab="Prob",main="PacBio Multipass Coverage")
lines(1:5,readsPacBioQual(0.7),type="l",col="red")
lines(1:5,readsPacBioQual(1),type="l",col="green")
lines(1:5,readsPacBioQual(2),type="l",col="blue")
legend(4,1,c("300bp","700bp","1kb","2kb"),col=c("black","red","green","blue"),pch=19)