[ -- Internship Proposals -- ]

Background

An online tool to calculate errors in Polymerase Chain Reaction (adapted from this work and using Opencpu).

Error rates are from Potapov V, Ong JL PLOS ONE 2017 .

The R code (speed improvements, comments) can be found on my GitHub.

Enzyme

Taq Q5 Phusion Deep Vent Pfu PrimeSTAR KOD Kapa HIFI DV Exo-

Length

100 200 400 800 1000 1500 2000 3000 5000

Cycles

15 20 22 25 30

# u<-1.5E-4 # l<- 2000 # numcycles<-25 mut<-10 N<-1000 lambda<-c(rep(0.872,20),rep(0.743,5),rep(0.146,5)) data<-data.frame(matrix(0,ncol = mut , nrow = (numcycles+1))) vec_p<-rep(0,mut) vec_p[1]<-N data[1,]<-vec_p shift <- function(x, n, invert=FALSE, default=0){ stopifnot(length(x)>=n) if(n==0){ return(x) } n <- ifelse(invert, n*(-1), n) return(c(rep(default, n), x[seq_len(length(x)-n)])) } pr_mutation<-sapply(1:mut, function(i){dmultinom(c((i-1),l-i+1),l,prob=c(u,1-u))}) n<-round(N*lambda[1]) vec<-round(pr_mutation*n) if (sum(vec) != n){ imax<-which.max(abs(round(pr_mutation*n)-pr_mutation*n)) if (sum(vec) > n){ vec[imax]<-vec[imax]-1 } else if (sum(vec) < n){ vec[imax]<-vec[imax]+1 } } vec<-vec+vec_p molecules<-vec data[2,]<-molecules for (k in 1:(numcycles-1)) { n<-round(sum(molecules)*lambda[k+1]) vec<-round(pr_mutation*n) if (sum(vec) != n){ imax<-which.max(abs(round(pr_mutation*n)-pr_mutation*n)) if (sum(vec) > n){ molecules[imax]<- molecules[imax]-1 } else if (sum(vec) < n){ molecules[imax]<- molecules[imax]+1 } } prob_multi<-sapply(1:mut,function(i){molecules[i]/sum(molecules)}) Error<-data.frame(matrix(0,ncol = mut , nrow = mut)) Error_shifted<-data.frame(matrix(0,ncol = mut , nrow = mut)) Error[1,]<-rmultinom(1,size=vec[1],prob=prob_multi) Error_shifted[1,]<-Error[1,] for (i in 2:mut) { Error[i,]<-rmultinom(1,size=vec[i],prob=prob_multi) Error_shifted[i,]<-shift(Error[i,],(i-1)) } SUM<-apply(Error_shifted,2,sum) molecules<-molecules+SUM data[(k+2),]<-molecules } res<-signif(100*molecules/sum(molecules),3) paste("[Code still under test] % molecules with 0 error:", res[1], ', 1 error: ', res[2], ', 2 errors: ', res[3])
Submit

.../...