diff --git a/R/genotype.R b/R/genotype.R index f940670..e0a7f2a 100644 --- a/R/genotype.R +++ b/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) } diff --git a/simulations.Rmd b/simulations.Rmd new file mode 100644 index 0000000..b70cb44 --- /dev/null +++ b/simulations.Rmd @@ -0,0 +1,38 @@ +--- +title: "R Notebook" +output: html_notebook +--- + +```{r} +library(genotype) +``` + + +```{r} +n.sim = 10000 + +homogenotypes = rep('A') +``` + +# Impacte de la profondeur de séquençage + +```{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) +} +``` diff --git a/simulations.nb.html b/simulations.nb.html new file mode 100644 index 0000000..86b8c00 --- /dev/null +++ b/simulations.nb.html @@ -0,0 +1,329 @@ + + + + + + + + + + + + + +R Notebook + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + + + + + + + + + + + + + + + + +
library(genotype)
+ + + + + + +
n.sim = 10000
+
+homogenotypes = rep('A')
+ + + +
+

Impacte de la profondeur de séquençage

+ + + +
p.hetero = 1/1000
+error.rate = 1/100
+merror = error_jc(rate = error.rate)
+couvertures =1:50
+ + + + + + +
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) 
+}
+ + +
+ +
LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3J9CmxpYnJhcnkoZ2Vub3R5cGUpCmBgYAoKCmBgYHtyfQpuLnNpbSA9IDEwMDAwCgpob21vZ2Vub3R5cGVzID0gcmVwKCdBJykKYGBgCgojIEltcGFjdGUgZGUgbGEgcHJvZm9uZGV1ciBkZSBz6XF1ZW7nYWdlCgpgYGB7cn0KcC5oZXRlcm8gPSAxLzEwMDAKZXJyb3IucmF0ZSA9IDEvMTAwCm1lcnJvciA9IGVycm9yX2pjKHJhdGUgPSBlcnJvci5yYXRlKQpjb3V2ZXJ0dXJlcyA9MTo1MApgYGAKCmBgYHtyfQpzdWNjZXNzLmhvbW8gPSBudW1lcmljKGxlbmd0aCA9IGxlbmd0aChjb3V2ZXJ0dXJlcykpCmkgPSAwCmZvciAoYyBpbiBjb3V2ZXJ0dXJlcykgewogIHMgPSBzaW11bGF0ZV9zaXRlcyhyZXAoJ0EnLG4uc2ltKSwKICAgICAgICAgICAgICAgICAgICAgcmVwKCdBJyxuLnNpbSksCiAgICAgICAgICAgICAgICAgICAgIGNvdmVyYWdlID0gYywKICAgICAgICAgICAgICAgICAgICAgZXJyb3IgPSBtZXJyb3IpCiAgZyA9IGdlbm90eXBpbmcocyxzZXF1ZW5jaW5nLmVycm9yID0gZXJyb3IucmF0ZSxwaGV0ZXJvID0gcC5oZXRlcm8pCiAgaSA9IGkrMQogIHN1Y2Nlc3MuaG9tb1tpXSA9IHN1bShnJG9rKSAKfQpgYGAK
+ + + +
+ + + + + + + +