kmer counts with cython implementation

Question

I have this function implemented in Cython:

def count_kmers_cython(str string, list alphabet, int kmin, int kmax):
"""
Count occurrence of kmers in a given string.
"""
counter = {}
cdef int i
cdef int j
cdef int N = len(string)
limits = range(kmin, kmax + 1)
for i in range(0, N - kmax + 1):
    for j in limits:
        kmer = string[i:i+j]
        counter[kmer] = counter.get(kmer, 0) + 1
return counter

Can I do better with cython? Or Can I have any away to improve it? I am new to cython, that is my first attempt.

I will use this to count kmers in DNA with alphabet restrict to 'ACGT'. The length of the general input string is the average bacterial genomes (130 kb to over 14 Mb, where each 1 kb = 1000 bp). The size of the kmers will be 3 < kmer < 16.

I wish to know if I could go further and maybe use cython in this function to:

def compute_kmer_stats(kmer_list, counts, len_genome, max_e):
"""
This function computes the z_score to find under/over represented kmers 
according to a cut off e-value.
Inputs:
    kmer_list - a list of kmers
    counts - a dictionary-type with k-mers as keys and counts as values.
    len_genome - the total length of the sequence(s).
    max_e - cut off e-values to report under/over represented kmers.
Outputs:
    results - a list of lists as [k-mer, observed count, expected count, z-score, e-value]
"""
print(colored('Starting to compute the kmer statistics...
',
              'red',
              attrs=['bold']))
results = []
# number of tests, used to convert p-value to e-value.
n = len(list(kmer_list))
for kmer in kmer_list:
    k = len(kmer)
    prefix, sufix, center = counts[kmer[:-1]], counts[kmer[1:]], counts[kmer[1:-1]]
    # avoid zero division error
    if center == 0:
        expected = 0
    else:
        expected = (prefix * sufix) // center
        observed = counts[kmer]
        sigma = math.sqrt(expected * (1 - expected / (len_genome - k + 1)))
        # avoid zero division error
        if sigma == 0.0:
            z_score = 0.0
        else:
            z_score = ((observed - expected) / sigma)
            # pvalue for all kmers/palindromes under represented
            p_value_under = (math.erfc(-z_score / math.sqrt(2)) / 2)
            # pvalue for all kmers/palindromes over represented
            p_value_over = (math.erfc(z_score / math.sqrt(2)) / 2)
            # evalue for all kmers/palindromes under represented
            e_value_under = (n * p_value_under)
            # evalue for all kmers/palindromes over represented
            e_value_over = (n * p_value_over)
            if e_value_under <= max_e:
                results.append([kmer, observed, expected, z_score, p_value_under, e_value_under])
            elif e_value_over <= max_e:
                results.append([kmer, observed, expected, z_score, p_value_over, e_value_over])
return results

OBS - Thank you CodeSurgeon by the help. I know there are other tools to count kmer efficiently but I am learning Python so I am trying to write my own functions and code.

kmer counts with cython implementation

0 Answers