KmerStream
KmerStream copied to clipboard
pmelsted
Inconsistent predictions of G from ground truth set
Pall,
I have a bacterial genome 2866389 bp (2.8 Mbp) which is finished/closed and a set of Illumina PE reads to 90x coverage.
I ran the following for Q = 20, 10 and 0 and then put it through KmerStreamEstimate
KmerStream -q Q -k 31,41,51,61,71,81,91,101,111,121 -o k.tab -t 36 --tsv R1.fq.gz R2.fq.gz
I am using git HEAD with the recent commits.
I was hoping to see some consistency with the G estimates, but it seems to be changing wrt k?
Q=20
Q k F0 f1 F1 F0-f1 G e lambda
20 31 3299585 493126 99828961 2806459.0 2805991.89701 0.00494906985143 35.5770667429
20 41 3344589 557769 81730492 2786820.0 2786365.13834 0.00683562471453 29.3322977938
20 51 3333071 592388 66978637 2740683.0 2740263.7089 0.00885695513175 24.4424055913
20 61 3295172 616007 54567924 2679165.0 2678777.41407 0.0113030219247 20.3704584462
20 71 3235748 625295 43968174 2610453.0 2610103.65399 0.0142374984179 16.845374678
20 81 3133043 612159 34858661 2520884.0 2520623.97029 0.0175775131132 13.8293777299
20 91 3003070 596893 27015406 2406177.0 2406397.84817 0.0220979239132 11.2264919205
20 101 2853560 593094 20294110 2260466.0 2263858.37302 0.0290888014808 8.96439028247
20 111 2650170 593007 14587743 2057163.0 2075764.26773 0.0395590939439 7.02764915398
20 121 2404494 630370 9798492 1774124.0 1849163.75171 0.0579818605305 5.29887739307
Q=10
Q k F0 f1 F1 F0-f1 G e lambda
10 31 23238820 19188481 211932108 4050339.0 3419501.83239 0.0965584341374 61.97748046
10 41 24951970 21071837 193657981 3880133.0 3289632.00191 0.114964377917 58.8691929333
10 51 25741248 22060720 175436035 3680528.0 3140455.77674 0.131954691927 55.8632400746
10 61 25728340 22204961 157291715 3523379.0 3057818.75638 0.14713136109 51.4391883665
10 71 24981307 21547726 139242131 3433581.0 3060875.69114 0.160134075816 45.4909460725
10 81 23760034 20483993 121304895 3276041.0 2974236.75104 0.173863978665 40.7852182439
10 91 22183540 19087225 103504595 3096315.0 2854603.09556 0.189099634088 36.2588393325
10 101 20059862 17058780 85868079 3001082.0 2826521.24778 0.202742260699 30.3794210171
10 111 17461606 14690693 68423270 2770913.0 2645626.03974 0.218375670316 25.8627897413
10 121 14471688 11808783 51195091 2662905.0 2587409.56694 0.233618232371 19.7862339438
Q=0
Q k F0 f1 F1 F0-f1 G e lambda
0 31 23280809 19248612 211932108 4032197.0 3391482.35372 0.0969371952455 62.4895210696
0 41 25008847 21113700 193657981 3895147.0 3305218.88169 0.115174411602 58.5915753031
0 51 25662770 21889372 175436035 3773398.0 3263244.82718 0.130631908465 53.7612236565
0 61 25668828 22092345 157291715 3576483.0 3126446.66034 0.14621551104 50.3100586987
0 71 24997343 21488823 139242131 3508520.0 3148709.00282 0.159523840018 44.221975062
0 81 23847969 20582358 121304895 3265611.0 2959234.83715 0.17475086963 40.9919799122
0 91 22182579 19119747 103504595 3062832.0 2816923.81421 0.189495578172 36.7438389629
0 101 20085270 17120760 85868079 2964510.0 2786033.1511 0.20355585144 30.8209107153
0 111 17491351 14696703 68423270 2794648.0 2670449.51692 0.218431506673 25.6223791412
0 121 14489590 11867384 51195091 2622206.0 2544646.33846 0.234843802093 20.1187450792
Thanks for reporting this. The error rate/Genome size model isn't doing to well with high coverage. This is mostly due to the assumption that k-mers at distance 2 from the truth would not be observed twice by chance. Once you have higher coverage you would need to model that explicitly.
I'm working on extending the model so that it can account better for this. Additionally I'll have to replace the brentp method that is not working in #11
I sense a Catch-22 however - how do I know I have too much depth until I know the genome size?