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

  • $\begingroup$ Hi @Lipsa, please provide a test example and your output together with any possible error messages. This appears certainly solvable. Is this method in Biopython? $\endgroup$ – Michael Nov 21 at 9:52
  • $\begingroup$ I've had a guess at what I think the question should read $\endgroup$ – Michael Nov 21 at 10:12

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.

| improve this answer | |

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
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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