How to count and compare k-mer count vectors (and print the top ten highest contributions)?

Question

[nihaldm169]

Hello

For a project I got two FASTA samples of amino acid sequences. I should write a function that compares the two samples and then prints out into a file the k-mers with the top ten highest contributions (length of k = 12). How can I do this, maybe using k-mer count vectors? Or what method is useful for this?

Happy about any help, also links to webpages or tutorials!

Thanks!

[Joe]

How big are the sequences you need to assess? (I.e. does the speed of the solution matter?)

[nihaldm169]

The speed doesn't matter, no. And I don't know how big the sequences are... We got sent two files with the sequences yet in them. But I need to compare sequences of 12 letters length.

Answer 1

[Joe]

Here's a solution:

import re, sys
from collections import Counter
from Bio import SeqIO

k = 12
n = 10
kmers = re.compile("(?=(\w{%s}))" % k)

for rec in SeqIO.parse(sys.argv[1], 'fasta'):
    kseqs = kmers.findall(str(rec.seq))
    counts = Counter(kseqs)
    top_n = counts.most_common(n)
    print(top_n)

This prints the top n kmers of length k for each each sequence in the given fasta.

On a 5Mb bacterial genome this takes:

real    0m6.697s
user    0m6.060s
sys 0m0.624s

NB I haven't sanity checked the ouput, nor have I 'prettified' it, so use the above code advisedly and please do your own sanity checking. This is also not a super memory efficient script, so just beware.

I'm not sure what you can do to compare them personally. You will likely need to normalise the results for the length of the input sequence though, if they arent identical. You'll need to look for what tests you can do with categorical data though.