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) 
}
plot(couvertures,success.homo,type = "b",cex=0.5,log="x")

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))

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))

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) 
}
plot(couvertures,success.hetero,type = "b",cex=0.5,col="blue",log="x")
points(couvertures,success.homo,type = "b",cex=0.5,col="red")

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))

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

p.hetero = 1/1000
error.rate = 10^seq(from = -4,to = -1, length.out = 50)
couverture = 10
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) 
}
plot(error.rate,success.homo,type = "b",cex=0.5,log="x")

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) 
}
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")
LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3J9CmxpYnJhcnkoZ2Vub3R5cGUpCmBgYAoKCmBgYHtyfQpuLnNpbSA9IDEwMDAwCgpob21vZ2Vub3R5cGVzID0gcmVwKCdBJykKYGBgCgojIEltcGFjdGUgZGUgbGEgcHJvZm9uZGV1ciBkZSBz6XF1ZW7nYWdlCgpgYGB7cn0KcC5oZXRlcm8gPSAxLzEwMDAKZXJyb3IucmF0ZSA9IDEvMTAwCm1lcnJvciA9IGVycm9yX2pjKHJhdGUgPSBlcnJvci5yYXRlKQpjb3V2ZXJ0dXJlcyA9MTo1MApgYGAKCmBgYHtyfQpzdWNjZXNzLmhvbW8gPSBudW1lcmljKGxlbmd0aCA9IGxlbmd0aChjb3V2ZXJ0dXJlcykpCmkgPSAwCmZvciAoYyBpbiBjb3V2ZXJ0dXJlcykgewogIHMgPSBzaW11bGF0ZV9zaXRlcyhyZXAoJ0EnLG4uc2ltKSwKICAgICAgICAgICAgICAgICAgICAgcmVwKCdBJyxuLnNpbSksCiAgICAgICAgICAgICAgICAgICAgIGNvdmVyYWdlID0gYywKICAgICAgICAgICAgICAgICAgICAgZXJyb3IgPSBtZXJyb3IpCiAgZyA9IGdlbm90eXBpbmcocyxzZXF1ZW5jaW5nLmVycm9yID0gZXJyb3IucmF0ZSxwaGV0ZXJvID0gcC5oZXRlcm8pCiAgaSA9IGkrMQogIHN1Y2Nlc3MuaG9tb1tpXSA9IHN1bShnJG9rKSAKfQpgYGAKCmBgYHtyfQpwbG90KGNvdXZlcnR1cmVzLHN1Y2Nlc3MuaG9tbyx0eXBlID0gImIiLGNleD0wLjUsbG9nPSJ4IikKYGBgCmBgYHtyfQpzID0gc2ltdWxhdGVfc2l0ZXMocmVwKCdBJyxuLnNpbSksCiAgICAgICAgICAgICAgICAgICAgIHJlcCgnQScsbi5zaW0pLAogICAgICAgICAgICAgICAgICAgICBjb3ZlcmFnZSA9IDEwLAogICAgICAgICAgICAgICAgICAgICBlcnJvciA9IG1lcnJvcikKZyA9IGdlbm90eXBpbmcocyxzZXF1ZW5jaW5nLmVycm9yID0gZXJyb3IucmF0ZSxwaGV0ZXJvID0gcC5oZXRlcm8pCgpoaXN0KGxvZzEwKGckc2NvcmUpKQpgYGAKCmBgYHtyfQpzID0gc2ltdWxhdGVfc2l0ZXMocmVwKCdBJyxuLnNpbSksCiAgICAgICAgICAgICAgICAgICAgIHJlcCgnQScsbi5zaW0pLAogICAgICAgICAgICAgICAgICAgICBjb3ZlcmFnZSA9IDMwLAogICAgICAgICAgICAgICAgICAgICBlcnJvciA9IG1lcnJvcikKZyA9IGdlbm90eXBpbmcocyxzZXF1ZW5jaW5nLmVycm9yID0gZXJyb3IucmF0ZSxwaGV0ZXJvID0gcC5oZXRlcm8pCgpoaXN0KGxvZzEwKGckc2NvcmUpKQpgYGAKCmBgYHtyfQpzdWNjZXNzLmhldGVybyA9IG51bWVyaWMobGVuZ3RoID0gbGVuZ3RoKGNvdXZlcnR1cmVzKSkKaSA9IDAKZm9yIChjIGluIGNvdXZlcnR1cmVzKSB7CiAgcyA9IHNpbXVsYXRlX3NpdGVzKHJlcCgnQScsbi5zaW0pLAogICAgICAgICAgICAgICAgICAgICByZXAoJ1QnLG4uc2ltKSwKICAgICAgICAgICAgICAgICAgICAgY292ZXJhZ2UgPSBjLAogICAgICAgICAgICAgICAgICAgICBlcnJvciA9IG1lcnJvcikKICBnID0gZ2Vub3R5cGluZyhzLHNlcXVlbmNpbmcuZXJyb3IgPSBlcnJvci5yYXRlLHBoZXRlcm8gPSBwLmhldGVybykKICBpID0gaSsxCiAgc3VjY2Vzcy5oZXRlcm9baV0gPSBzdW0oZyRvaykgCn0KYGBgCgoKYGBge3J9CnBsb3QoY291dmVydHVyZXMsc3VjY2Vzcy5oZXRlcm8sdHlwZSA9ICJiIixjZXg9MC41LGNvbD0iYmx1ZSIsbG9nPSJ4IikKcG9pbnRzKGNvdXZlcnR1cmVzLHN1Y2Nlc3MuaG9tbyx0eXBlID0gImIiLGNleD0wLjUsY29sPSJyZWQiKQoKYGBgCgoKYGBge3J9CnMgPSBzaW11bGF0ZV9zaXRlcyhyZXAoJ0EnLG4uc2ltKSwKICAgICAgICAgICAgICAgICAgICAgcmVwKCdUJyxuLnNpbSksCiAgICAgICAgICAgICAgICAgICAgIGNvdmVyYWdlID0gMTAsCiAgICAgICAgICAgICAgICAgICAgIGVycm9yID0gbWVycm9yKQpnID0gZ2Vub3R5cGluZyhzLHNlcXVlbmNpbmcuZXJyb3IgPSBlcnJvci5yYXRlLHBoZXRlcm8gPSBwLmhldGVybykKCmhpc3QobG9nMTAoZyRzY29yZSkpCmBgYAoKCmBgYHtyfQpzID0gc2ltdWxhdGVfc2l0ZXMocmVwKCdBJyxuLnNpbSksCiAgICAgICAgICAgICAgICAgICAgIHJlcCgnVCcsbi5zaW0pLAogICAgICAgICAgICAgICAgICAgICBjb3ZlcmFnZSA9IDMwLAogICAgICAgICAgICAgICAgICAgICBlcnJvciA9IG1lcnJvcikKZyA9IGdlbm90eXBpbmcocyxzZXF1ZW5jaW5nLmVycm9yID0gZXJyb3IucmF0ZSxwaGV0ZXJvID0gcC5oZXRlcm8pCgpoaXN0KGxvZzEwKGckc2NvcmUpKQpgYGAKCgojIEltcGFjdCBkdSB0YXV4IGQnZXJyZXVyCgpgYGB7cn0KcC5oZXRlcm8gPSAxLzEwMDAKZXJyb3IucmF0ZSA9IDEwXnNlcShmcm9tID0gLTQsdG8gPSAtMSwgbGVuZ3RoLm91dCA9IDUwKQpjb3V2ZXJ0dXJlID0gMTAKYGBgCgoKYGBge3J9CnN1Y2Nlc3MuaG9tbyA9IG51bWVyaWMobGVuZ3RoID0gbGVuZ3RoKGVycm9yLnJhdGUpKQppID0gMApmb3IgKGUgaW4gZXJyb3IucmF0ZSkgewogIG1lcnJvciA9IGVycm9yX2pjKHJhdGUgPSBlKQogIHMgPSBzaW11bGF0ZV9zaXRlcyhyZXAoJ0EnLG4uc2ltKSwKICAgICAgICAgICAgICAgICAgICAgcmVwKCdBJyxuLnNpbSksCiAgICAgICAgICAgICAgICAgICAgIGNvdmVyYWdlID0gY291dmVydHVyZSwKICAgICAgICAgICAgICAgICAgICAgZXJyb3IgPSBtZXJyb3IpCiAgZyA9IGdlbm90eXBpbmcocyxzZXF1ZW5jaW5nLmVycm9yID0gZSxwaGV0ZXJvID0gcC5oZXRlcm8pCiAgaSA9IGkrMQogIHN1Y2Nlc3MuaG9tb1tpXSA9IHN1bShnJG9rKSAKfQoKYGBgCgpgYGB7cn0KcGxvdChlcnJvci5yYXRlLHN1Y2Nlc3MuaG9tbyx0eXBlID0gImIiLGNleD0wLjUsbG9nPSJ4IikKYGBgCgoKYGBge3J9CnN1Y2Nlc3MuaGV0ZXJvID0gbnVtZXJpYyhsZW5ndGggPSBsZW5ndGgoZXJyb3IucmF0ZSkpCmkgPSAwCmZvciAoZSBpbiBlcnJvci5yYXRlKSB7CiAgbWVycm9yID0gZXJyb3JfamMocmF0ZSA9IGUpCiAgcyA9IHNpbXVsYXRlX3NpdGVzKHJlcCgnQScsbi5zaW0pLAogICAgICAgICAgICAgICAgICAgICByZXAoJ1QnLG4uc2ltKSwKICAgICAgICAgICAgICAgICAgICAgY292ZXJhZ2UgPSBjb3V2ZXJ0dXJlLAogICAgICAgICAgICAgICAgICAgICBlcnJvciA9IG1lcnJvcikKICBnID0gZ2Vub3R5cGluZyhzLHNlcXVlbmNpbmcuZXJyb3IgPSBlLHBoZXRlcm8gPSBwLmhldGVybykKICBpID0gaSsxCiAgc3VjY2Vzcy5oZXRlcm9baV0gPSBzdW0oZyRvaykgCn0KCmBgYAoKYGBge3J9CnBsb3QoZXJyb3IucmF0ZSxzdWNjZXNzLmhldGVybyx0eXBlID0gImIiLGNleD0wLjUsY29sPSJibHVlIixsb2c9Inh5Iix5bGltPWMoMjAwMCwxMDAwMCkpCnBvaW50cyhlcnJvci5yYXRlLHN1Y2Nlc3MuaG9tbyx0eXBlID0gImIiLGNleD0wLjUsY29sPSJyZWQiKQpgYGAK