add a column on the result table with the log likelyhood of the genotypes

R/genotype.R

@@ -213,11 +213,12 @@ genotyping = function(sequencing,sequencing.error,phetero) {
   g.good = rownames(pnucs)[best.genotype]
   g.alt  = rownames(pnucs)[second.genotype]
 
-  score = apply(probas, MARGIN = 1, max) - apply(probas, MARGIN = 1, second)
+  lod = apply(probas, MARGIN = 1, max)
+  score = lod - apply(probas, MARGIN = 1, second)
   
   true = apply(sequencing[,1:2],MARGIN = 1,function(x) paste0(sort(x),collapse = ""))
    
-  sequencing = cbind(sequencing,probas,g.good,g.alt,score = exp(score),ok=true==g.good)
+  sequencing = cbind(sequencing,probas,g.good,g.alt,lod,score = exp(score),ok=true==g.good)
   
   return(sequencing)
 }

simulations.Rmd

+---
+title: "R Notebook"
+output: html_notebook
+---
+
+```{r}
+library(genotype)
+```
+
+
+```{r}
+n.sim = 10000
+
+homogenotypes = rep('A')
+```
+
+# Impacte de la profondeur de squenage
+
+```{r}
+p.hetero = 1/1000
+error.rate = 1/100
+merror = error_jc(rate = error.rate)
+couvertures =1:50
+```
+
+```{r}
+success.homo = numeric(length = length(couvertures))
+i = 0
+for (c in couvertures) {
+  s = simulate_sites(rep('A',n.sim),
+                     rep('A',n.sim),
+                     coverage = c,
+                     error = merror)
+  g = genotyping(s,sequencing.error = error.rate,phetero = p.hetero)
+  i = i+1
+  success.homo[i] = sum(g$ok) 
+}
+```