Counting the number of k-mers like monomers, dimers to hexamers from the fasta file

Question

Anyone here who uses Python for solving bioinformatics problems. This is the code I have written for counting the number of k-mers like monomers, dimers to hexamers from the fasta file. You just have to give the ncbi accession number for the fasta sequence and then it counts the number of k-mers. If you've time, please check the code as I think it's a bit long one and I have used try/except for solving IndexError. Your suggestions would be valuable. Thanks.

    from Bio import Entrez
    Entrez.email = '[email protected]'
    monomers = list('ATGC')
    dimers = []
    for i in monomers:
        for j in monomers:
            dimers.append(i+j)
    
    trimers = []
    for i in monomers:
        for j in monomers:
            for k in monomers:
                trimers.append(i+j+k)
     
    tetramers = []
    for i in monomers:
        for j in monomers:
            for k in monomers:
                for l in monomers:
                    tetramers.append(i+j+k+l)
     
    pentamers = []
    for i in monomers:
        for j in monomers:
            for k in monomers:
                for l in monomers:
                    for m in monomers:
                        pentamers.append(i+j+k+l+m)
     
    hexamers = []
    for i in monomers:
        for j in monomers:
            for k in monomers:
                for l in monomers:
                    for m in monomers:
                        for n in monomers:
                            hexamers.append(i+j+k+l+m+n)
     
    file = input('Enter the ncbi accession number: ')
    handle = Entrez.efetch(db = 'nucleotide', id = file,rettype="fasta", retmode="text")
    record = handle.read()
    fasta_string = ''.join(record.split('\n')[1:])
    k = int(input('Enter the value of k: '))
    print('The sequence is',fasta_string)
    fasta_list = []
    if k == 1:
        a = True
        while a:
            try:
                for i in range(0,len(fasta_string),1):
                    fasta_list.append(fasta_string[i])
            except:
                break
            a = False
        for i in monomers:
                    print('count of' ,i, 'is' , fasta_list.count(i))
    elif k == 2:
        a = True
        while a:
            try:
                for i in range(0,len(fasta_string),2):
                    fasta_list.append(fasta_string[i]+fasta_string[i+1])
            except:
                break
            a = False
        for i in dimers:
                    print('count of' ,i, 'is' , fasta_list.count(i))
        
    elif k == 3:
        a = True
        while a:
            try:
                for i in range(0,len(fasta_string),3):
                    fasta_list.append(fasta_string[i]+fasta_string[i+1]+fasta_string[i+2])
            except:
                break
            a = False
        for i in trimers:
                    print('count of' ,i, 'is' , fasta_list.count(i))
     
    elif k == 4:
        a = True
        while a:
            try:
                for i in range(0,len(fasta_string),4):
                    fasta_list.append(fasta_string[i]+fasta_string[i+1]+fasta_string[i+2]+fasta_string[i+3])
            except:
                break
            a = False
        for i in tetramers:
                    print('count of' ,i, 'is' , fasta_list.count(i))
     
    elif k == 5:
        a = True
        while a:
            try:
                for i in range(0,len(fasta_string),5):
                    fasta_list.append(fasta_string[i]+fasta_string[i+1]+fasta_string[i+2]+fasta_string[i+3]+fasta_string[i+4])
            except:
                break
            a = False
        for i in pentamers:
                    print('count of' ,i, 'is' , fasta_list.count(i))
    elif k == 6:
        a = True
        while a:
            try:
                for i in range(0,len(fasta_string),6):
                    fasta_list.append(fasta_string[i]+fasta_string[i+1]+fasta_string[i+2]+fasta_string[i+3]+fasta_string[i+4]+fasta_string[i+5])
            except:
                break
            a = False
        for i in hexamers:
                    print('count of' ,i, 'is' , fasta_list.count(i))

Counting the number of k-mers like monomers, dimers to hexamers from the fasta file

Answer 1

The code can be simplified quite a bit.

Using itertools.product, the code like this:

trimers = []
    for i in monomers:
        for j in monomers:
            for k in monomers:
                trimers.append(i+j+k)

can be reduced to:

k_mers = list(''.join(t) for t in itertools.product('ACGT', repeat=k))

A common Python idiom for grouping a sequence is

zip(*[iter(sequence)]*k)

it generates k-tuples from the sequence. Which can be counted using a collections.Counter. So this code:

a = True
while a:
    try:
        for i in range(0,len(fasta_string),3):
            fasta_list.append(fasta_string[i]+fasta_string[i+1]+fasta_string[i+2])

    except:
        break
    a = False
for i in trimers:
            print('count of' ,i, 'is' , fasta_list.count(i))

can be simplified to:

counts = Counter(''.join(t) for t in zip(*[iter(fasta_string)]*k))

The code asks for k, so it doesn't make sense to generate all the other k-mers.

The final code could look like:

from collections import Counter
from itertools import product

file = input('Enter the ncbi accession number: ')
k = int(input('Enter the value of k: '))

handle = Entrez.efetch(db = 'nucleotide', id = file,rettype="fasta", retmode="text")
record = handle.read()
fasta_string = ''.join(record.split('\n')[1:])

print('The sequence is',fasta_string)

counts = Counter(''.join(t) for t in zip(*[iter(fasta_string)]*k))

for k_mer in (''.join(t) for t in itertools.product('ACGT', repeat=k)):
    print(f"count of {k_mer} is {counts[k_mer]}")

Answer 2

A Quick Preface

A monomer can mean different things in different contexts; it's just a way of referring to the most relevant "unit" element of the current context. This usually means amino acids when you're doing sequence alignments, but I suppose it could also mean nucleotides¹, although I've never seen that myself.

If you have a FASTA file with the base pairs, though, you're usually² trying to parse the order and type of the codons in the sequence. From there, you either compare the sequence against others to determine the impact of mutations. Some mutations might be harmless, since different base codons sometimes code for the same amino acid, but others can be extremely problematic, to say the least. For example, a deletion or an insertion can cause a frameshift, moving the entire sequence forward or backwards.

There are two reasons I bring this up. First, I think calling each base a k-mer will be confusing to researchers expecting a monomer to represent codons/amino acids. Second, since codons are three bases long, searching for all substrings of a length that isn't a multiple of three won't be very helpful.

The key point is that the counts of each substring of length three aren't really what matters. Each of these substrings is called a codon, which in turn represents an amino acid or stop command; it's the combination of multiple amino acids in a particular order that result in the production of a specific protein.

The last point I want to make about the code before providing feedback on the actual code itself is that if you limit your search to only substrings of length 3, you could implement the parsing mechanism as a deterministic finite automaton. There are only twenty-two amino acids and three stop codons, so writing a state transition table wouldn't take too long, and it would reduce the runtime complexity of the parsing the sequence to \$O\left(n\right)\$, since it would depend only on how long the single pass takes, which itself is a factor of only the length of the input sequence.

Anyways, on to the actual review.

---

Recommendations

The following recommendations are focused on your actual code, not the suggestions made above.

Defer the Preprocessing

Since you are searching for only one kind of \$k\$-mer, but you don't know the value of \$k\$ until the user chooses, I would suggest creating all of the possible permutations of length \$k\$ beforehand is a waste of effort.

Don't Print the Sequence

FASTA files can be ginormous (meaning several Gigabytes long), so printing the sequence is not very practical, since it would take both a ton of time and a ton of memory.

It also isn't super useful, since no one is going to be checking all several Gigabytes of the sequence to make sure it's the right one. Certainly not while it's scrolling by in the console at lightspeed. The user selected the sequence by its sequence number, so I would assume they know what sequence they wanted.

Use argparse Instead of Standard Input

Bioinformatics happens on the central supercomputer, not the researchers' computers a lot of the time. Since you can't interact with the program as its running (you usually submit a slurm request via a bash script with the execution parameters), you're better off relying on the argv contents so the script execution can be defined when the request is submitted.

You could also define the input using a redirection operator, but I like the argparse route better, although this really is just a personal preference. It seems less "clean" to me, but if it works, it works.

Use a Separate Resource Download Script

I'm not sure if the Entrez.Bio package includes a built-in caching mechanism, where it will know not to re-download a file you previously requested, but I also kind of feel like that's irrelevant.

Bioinformatics research depends on access to the supercomputer, and I've even seen grants come in the form of not money but the amount of computing hours that the grant money would have cost. In other words, I would not waste that precious time downloading a file, when you can do that for free and just include it in the slurm request.

Not to mention, (I'm breaking out in cold sweat even considering this possibility) can you imagine submitting a job request with an incorrect sequence ID? These jobs can take days to complete (and we're talking about programs written in C and/or Fortran, who knows about Python?), so accidentally submitting an incorrect job request... I'm not saying the PI would murder you for it, but if they did, a jury of bioinformatics researchers would probably not convict them for it.

Furthermore, there are a lot of things that have to be done before actual analysis can take place. You need to have done some analysis before hand to be able to conduct some kind of regression testing on the results you get back.

Separating the downloading and processing of a sequence file allows you to be able to analyze arbitrary files, even contrived ones you wrote yourself. This then allows you to perform basic unit testing on the script, to ensure you didn't accidentally count adenine twice and forget guanine or something.

Don't Build the k-mers in Memory

Since by the time you start parsing the input sequence you already know the value of \$k\$, there's no need to actually build a list of \$k\$-mers. What I would do is use an input buffer \$k\$-characters long and then print out the \$k\$-mer once the buffer is full³.

More specifically, I would open an output file and write out the specific \$k\$-mer found. Remember, it's the ordering, not necessarily the counts, of the coding sequences that matter.

---

1. To be clear, however, a nucleotide is not the same thing as the A/T/G/C bases we're parsing here. When these bases combine with a five-carbon sugar, they form a nucleoside, which is itself still only a subcomponent of a nucleotide.
2. I've never seen anything else, but I'm not a microbiologist. I was just an intern studying math, so feedback from actual experts is always welcome.
3. If you were parsing codons, as I suggest in the first section, I would use a lookup table here to output the resulting amino acid/codon character. Each one has a single-character representation, so you can output a result file 1/3 the length of the input, allowing for easier post-processing.