This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How To Compare All The K-Mers Of A Given Length With All The K-Length Sub Strings Of A Dna Sequence?

if you have a set of k-mers of a given length, how you can compare each k-mer with each k-length sub string of a DNA sequence?

For example,

k-mers, (for k=4)

AAAA AAAT AAAG AAAC . . . . TTTT

sequence =ATGCCCATCAAAGGCTCATTGCGACC

3 answers

In R, after installing the Biostrings Bioconductor package:

library(Biostrings)
s = DNAString('ATGCCCATCAAAGGCTCATTGCGACC')
kmercounts = oligonucleotideFrequency(s,4)
head(kmercounts)
kmercounts[kmercounts>0]

The last line above returns:

AAAG AAGG AGGC ATCA ATGC ATTG CAAA CATC CATT CCAT CCCA CGAC CTCA GACC GCCC GCGA 
   1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1 
GCTC GGCT TCAA TCAT TGCC TGCG TTGC 
   1    1    1    1    1    1    1

If you want to know the kmer count for a specific kmer, you can do this:

x['AAAG']

which returns:

AAAG 
   1

This is great. Thanks

Jellyfish is a great program for this purpose.

Jellyfish has pretty much streamlined it, as far as I've heard

If it's a short sequence you can do it in python like this:

seq = 'AGATAGATAGACACAGAAATGGGACCACAC'
kmers = {}
k = 4
for i in range(len(seq) - k + 1):
   kmer = seq[i:i+k]
   if kmers.has_key(kmer):
      kmers[kmer] += 1
   else:
      kmers[kmer] = 1

for kmer, count in kmers.items():
   print kmer + "\t" + str(count)

If it's longer, like a whole genome, I would use jellyfish like Alastair suggested.

if you want a sorted list of kmers you can append this to the above dode:

import operator
sortedKmer = kmers.items()
sortedKmer.sort(key = operator.itemgetter(1), reverse = True)
for item in sortedKmer:
   print item[0] + "\t" + str(item[1])

thanks.Can you please mention how to do it in C?

You can use a collections.Counter to more efficiently count the kmers after they are generated by your substring iterator.

Log in to answer this question.