How to count the kmer occurrence in FASTA file considering overlapping and reverse complement?

Question

I am using count_overlap() for counting the kmers from Biopython .

Does it take the reverse complement of kmer into account?

I need to count reverse complement too

Answer 1

There's a lot detail missing here, but at a guess could you try,

from Bio import SeqIO
from Bio.Seq import Seq
for record in SeqIO.parse("ls_orchid.gbk", "genbank"):
    reversec = record.seq.reverse_complement())
    revcount = reversec.count_overlap('AT')
    print (revcount + count_overlap('AT'))

Essentially count_overlap only counts the forward strand if this is the Biopython method. Not debugged the above.

Answer 2

Here's a C++ kmer counter I wrote for my own purposes. It can count canonical and reverse complement kmers from FASTA files up to length 31. Run kmer-counter --help or read the README for count options.

For longer kmers, you'd want to look at an approximation tool like Jellyfish.

Answer 3

I would NOT recommend using biopython for a task like this. First, it's slow, but also it's really not necessary. There are several other existing tools optimised for exactly this purpose.

My favourite is KMC. It's very fast, and by default, it reports canonical kmers (Parameter -b would turn this off). So, assuming you have KMC installed, this is how you get all the kmers from a range of occurances (1 - 10000).

# just making a directory for temprary files
mkdir tmp

# building a kmer database with k = 21, 
# 16 cores and 
# up to 64GB of memory. 
# Counting occurances between 1 to 10000, 
# `-fa` specifies that the input is fasta, you can also input other formats
# the database name will be kmer_counts
# & saving temporary data to tmp directory; 
kmc -k21 -t16 -m64 -ci1 -cs10000 -fa my_seq.fasta kmer_counts tmp # run kmc

# extract kmers from the database in a text file
kmc_dump -ci1 -cx10000 kmer_counts kmer_k21.dump