Question: How to count and compare k-mer count vectors (and print the top ten highest contributions)?
0
11 months ago by
nihaldm1690 wrote:

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?

Thanks!

modified 11 months ago by Joe17k • written 11 months ago by nihaldm1690

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

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.

0
11 months ago by
Joe17k
United Kingdom
Joe17k wrote:

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.