--- 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) } ``` ```{r} plot(couvertures,success.homo,type = "b",cex=0.5,log="x") ``` ```{r} s = simulate_sites(rep('A',n.sim), rep('A',n.sim), coverage = 10, error = merror) g = genotyping(s,sequencing.error = error.rate,phetero = p.hetero) hist(log10(g$score)) ``` ```{r} s = simulate_sites(rep('A',n.sim), rep('A',n.sim), coverage = 30, error = merror) g = genotyping(s,sequencing.error = error.rate,phetero = p.hetero) hist(log10(g$score)) ``` ```{r} success.hetero = numeric(length = length(couvertures)) i = 0 for (c in couvertures) { s = simulate_sites(rep('A',n.sim), rep('T',n.sim), coverage = c, error = merror) g = genotyping(s,sequencing.error = error.rate,phetero = p.hetero) i = i+1 success.hetero[i] = sum(g$ok) } ``` ```{r} plot(couvertures,success.hetero,type = "b",cex=0.5,col="blue",log="x") points(couvertures,success.homo,type = "b",cex=0.5,col="red") ``` ```{r} s = simulate_sites(rep('A',n.sim), rep('T',n.sim), coverage = 10, error = merror) g = genotyping(s,sequencing.error = error.rate,phetero = p.hetero) hist(log10(g$score)) ``` ```{r} s = simulate_sites(rep('A',n.sim), rep('T',n.sim), coverage = 30, error = merror) g = genotyping(s,sequencing.error = error.rate,phetero = p.hetero) hist(log10(g$score)) ``` # Impact du taux d'erreur ```{r} p.hetero = 1/1000 error.rate = 10^seq(from = -4,to = -1, length.out = 50) couverture = 10 ``` ```{r} success.homo = numeric(length = length(error.rate)) i = 0 for (e in error.rate) { merror = error_jc(rate = e) s = simulate_sites(rep('A',n.sim), rep('A',n.sim), coverage = couverture, error = merror) g = genotyping(s,sequencing.error = e,phetero = p.hetero) i = i+1 success.homo[i] = sum(g$ok) } ``` ```{r} plot(error.rate,success.homo,type = "b",cex=0.5,log="x") ``` ```{r} success.hetero = numeric(length = length(error.rate)) i = 0 for (e in error.rate) { merror = error_jc(rate = e) s = simulate_sites(rep('A',n.sim), rep('T',n.sim), coverage = couverture, error = merror) g = genotyping(s,sequencing.error = e,phetero = p.hetero) i = i+1 success.hetero[i] = sum(g$ok) } ``` ```{r} plot(error.rate,success.hetero,type = "b",cex=0.5,col="blue",log="xy",ylim=c(2000,10000)) points(error.rate,success.homo,type = "b",cex=0.5,col="red") ```