Feature requests: different mixture model; combined ska build & cov

[rderelle]

Hi John,

Thanks a lot for implementing the calculation of the optimal kmer count threshold (ska cov). I tested it using ska 0.3.1 and a paired-end sample (fastq files available here: https://drive.google.com/drive/folders/1HVO-6mOd7bh7CPOjXhA3lWAT0GA_8SC8?usp=sharing). My command line was:

./ska_0.3.1 cov reads_kraken/ERR8158003_1.fastq.gz reads_kraken/ERR8158003_2.fastq.gz -k 31 > cov.txt

Given the distribution of kmer counts, the cutoff should be around 4-6, but ska outputs a value of 17. I think this high treeshold is related to the "Error" message below obtained for low kmer sizes during the model fittng.

This is the beginning of the file 'cov.txt':

Count K_mers Mixture_density Component 1 25838523 3.2408096206707826e-1 Error 2 1993278 1.620404810335391e-1 Error 3 129884 5.4013493677846365e-2 Error 4 97141 1.3503373419461602e-2 Error 5 28053 2.7006746838923196e-3 Error 6 35628 4.5011244731538646e-4 Error 7 24702 6.430177818791246e-5 Error 8 35486 8.03772227348904e-6 Error 9 30571 8.93080252609934e-7 Error 10 39312 8.930802526126602e-8 Error 11 40110 8.118911389065957e-9 Error 12 44332 6.765759585478362e-10 Error 13 43658 5.20443537193471e-11 Error 14 49357 3.7176916177683685e-12 Error 15 47325 2.4891833425171915e-13 Error 16 51653 2.009020801325787e-14 Error 17 52514 1.921693705685228e-14 Coverage 18 53472 6.883935794424811e-14 Coverage 19 52873 2.4488922538072723e-13 Coverage 20 53588 8.282019181080119e-13 Coverage 21 54322 2.667583813189892e-12 Coverage 22 54876 8.201562459958256e-12 Coverage 23 55197 2.411959247393554e-11 Coverage 24 54703 6.797667675764525e-11 Coverage 25 54701 1.8391668285808207e-10 Coverage 26 54189 4.784636868075234e-10 Coverage 27 52711 1.1986335327845915e-9 Coverage 28 52430 2.895540187815414e-9 Coverage 29 51131 6.753560645868026e-9 Coverage 30 51195 1.522694413626287e-8 Coverage 31 49749 3.322402658582442e-8 Coverage 32 48342 7.02268990941827e-8 Coverage 33 48280 1.4394306882794722e-7 Coverage 34 47080 2.863604561107574e-7 Coverage 35 46010 5.534089852758182e-7 Coverage 36 45616 1.039788261971188e-6 Coverage 37 45274 1.9008348747272694e-6 Coverage 38 44701 3.383467426815247e-6 Coverage 39 44463 5.868114750185855e-6 Coverage 40 43372 9.922927345843403e-6 Coverage 41 43026 1.637031965870808e-5 Coverage

Also, would it be possible to perform the calculation of the estimated cutoff 'on the fly' while extracting the kmers from fastq files? As it is, we have to run ska twice (first to estimate the cutoff and then to build the kmer dictionary), which kind of defeat the purpose (twice the runtimes).

Thanks and best, Romain

[rderelle]

Here is the coverage histogram

[johnlees]

Thanks for the testing and files. The cov.txt file is actually fine, the coverage/error column is actually just to assist with the plotting.

It looks like the default mixture model just doesn't fit this data that well. I'll have a play around with some other types of mixture for the second component to see if we can offer an alternative in these cases.

I'll also think about adding an auto option rather than running ska cov separately, that might be easier after fixing #45. Of course you do end up counting k-mers twice this way if you run on every read set, but it's still only about ~60s per sample.