6
$\begingroup$

I'm looking for some software that can "efficiently" (time and memory) compute the exact percentage of non-unique k-mers in a genome for given k. I don't need the k-mers or the abundances itself, I just need the percentage.

Alternatively, the result could be:

  1. the number of different k-mers observed in the genome
  2. the number of unique k-mers, i.e. k-mers that occur only once in the genome, or the number of non-unique k-mers, i.e., k-mers that occur more than once in the genome.

which can be easily used to compute the requested percentage.

Requested features:

  • don't count k-mers with ambiguous nucleotides
  • count a k-mer and its reverse complement together
  • k upto 50 or 100
  • for large genomes (upto some GB)
  • efficiency

For the study in mind, I like to vary k. Hence, I will run the software several times with different k.

Any help is welcome.

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ Welcome to the site. What is your current solution? Which non-efficient programs have you found (To avoid proposing them)? Can't the other answers in the k-mer help you to make that program? $\endgroup$
    – llrs
    Commented Mar 14, 2018 at 10:14

1 Answer 1

Reset to default
5
$\begingroup$

You can do all of that with khmer. For example, abundance-dist-single.py produces a file with columns: k-mer abundance, k-mer count, cumulative count, and fraction of total distinct k-mers. So for question 1 you would sum column 2. For question 2 you would just get thek-mer countassociated with ak-mer abundance` of 1.

That package also provides a python API (see the read the docs link above) if you need to customize things more.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks a lot for that hint. I made a first test of khmer using A.thaliana: python abundance-dist-single.py -k <k> -b <genome> <out>. It looks great. It's fast and easy to use. However, I get an error for k > 32. Is there a way to go further than 32? $\endgroup$
    – Jens
    Commented Mar 20, 2018 at 21:49
  • $\begingroup$ I'm not sure it's possible with the command line script, but it appears to be possible in the API. $\endgroup$
    – Devon Ryan
    Commented Mar 21, 2018 at 7:43
  • $\begingroup$ Thanks for that hint, I'll check that. However, I also found out for some examples (and reading the documentation/papers afterwards ;) that it is not exact. Hence, I currently also use the parameters -N 20 -x 1E7 to avoid too many mistakes. Just in case anyone like to do a similar case study. In addition, I have some problems with some species where khmer states ERROR: abundance distribution is uniformly zero; nothing to report. Please verify that the input files are valid. The files seem to be okay and have been used in other analysis. Any ideas? $\endgroup$
    – Jens
    Commented Mar 21, 2018 at 8:45
  • $\begingroup$ No, sorry, I've never seen that error. $\endgroup$
    – Devon Ryan
    Commented Mar 21, 2018 at 12:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.