Calculate the occurrence of motif in sequences (including overlaps) using Python

Question

I need to calculate the occurrences of a motif (including overlaps) in sequences (motif is passed in the first line of standard input and FASTA sequences in subsequent lines). The sequence name starts with >, and after whitespace is just a comment about the sequence that needs to be neglected. The input of program is like:

 AT
 >seq1 Comment......
 AGGTATA
 TGGCGCC
 >seq2 Comment.....
 GGCCGGCGC
 GCCGGCGGG

The output should be:

seq1: 2
seq2: 0

I decided to save the first line as a motif, strip the comment from sequence name, join lines of sequence in one line and save sequence names (keys) and sequences (values) in a dictionary. I also wrote a function for motif_count and want to call it on dictionary values and then save it in a separate dictionary for final output. However, it does not work the way I want it as the method motif_count does not produce any result when called on my d values.

Here is my code:

 #!/usr/bin/env python3

 import sys

 sequence = sys.stdin.readlines()
 d = {}
 temp_genename = None
 temp_sequence = None
 motif = sequence[0]

def motif_count(motif, temp_sequence):
    count = 0
    start = -1
    while True:
        start = temp_sequence.find(motif, start +1)
    if start <0:
        return count
    count += 1

if sequence[1][0] != '>':

    print("ERROR")

exit(1)

for line in sequence[1:]:
    if line[0] == '>':
       temp_genename = line.split(' ')[0].strip()
       tmp_genename = temp_genename.replace('>', '')
       temp_sequence = ""

else:

      temp_sequence += line.strip()

d[tmp_genename] = temp_sequence
for key in d:
    motif_count(motif, key)

Answer 1

I tried to fix the indentation of the code you posted, a typo in some variable names, plus a few other issues that caused it to not work as expected (see the comment in the code):

#!/usr/bin/env python3

import sys

sequence = sys.stdin.readlines()
d = {}
temp_genename = None
temp_sequence = None
# The lines returned by readlines() contain end-of-line characters,
# this is why strip() is needed
motif = sequence[0].strip()

def motif_count(motif, temp_sequence):
    count = 0
    start = -1
    while True:
        start = temp_sequence.find(motif, start +1)
        if start <0:
            return count
        count += 1

if sequence[1][0] != '>':

    print("ERROR")

    exit(1)

for line in sequence[1:]:
    if line[0] == '>':
        # Record the previous one in the dictionary
        if temp_genename is not None:
            d[temp_genename] = temp_sequence
        temp_genename = line.split(' ')[0].strip()
        temp_genename = temp_genename.replace('>', '')
        temp_sequence = ""

    else:

        temp_sequence += line.strip()

# This should be done each time a new sequence is found,
# not only at the end
d[temp_genename] = temp_sequence
for key in d:
    # We need to use the result
    # motif_count(motif, key)
    # We want to count the motif in the sequence, not in its name
    count = motif_count(motif, d[key])
    print(f"{key}: {count}")

Testing it:

$ cat tmp.fa 
AT
>seq1 Comment......
AGGTATA
TGGCGCC
>seq2 Comment.....
GGCCGGCGC
GCCGGCGGG
$ cat tmp.fa | ./count_motif.py
seq1: 2
seq2: 0

Here is another way to do, using valmap to make a dict of the motif counts from the dict of sequences, and with a few other programming style changes:

#!/usr/bin/env python3

import sys
from cytoolz import valmap

motif = sys.stdin.readline().strip()


# We can make this a function of only one parameter
# because motif is fixed
def motif_count(temp_sequence):
    count = 0
    start = -1
    while True:
        start = temp_sequence.find(motif, start +1)
        if start <0:
            return count
        count += 1


seq_dict = {}
temp_genename = None
temp_sequence = None

for line in sys.stdin:
    if line[0] == '>':
        # Record the previous one in the dictionary
        if temp_genename is not None:
            seq_dict[temp_genename] = temp_sequence
        temp_genename, *_ = line[1:].strip().split()
        temp_sequence = ""
    else:
        temp_sequence += line.strip()

msg = f"""File format might be wrong"
motif: {motif}
temp_genename: {temp_genename}
temp_genename: {temp_sequence}"""

assert temp_sequence is not None and len(temp_sequence), msg

# Record the last sequence
seq_dict[temp_genename] = temp_sequence

# Make a new dict where the values are the motif counts
# valmap needs the function motif_count to take only one parameter
counts_dict = valmap(motif_count, seq_dict)
# Print the results, iterating over dict items
for (name, count) in counts_dict.items():
    print(f"{name}: {count}")

Answer 2

This is a good task for Biopython:

from Bio import SeqIO

with open('input_file') as f:
    motif = next(f).rstrip().lower()
    for r in SeqIO.parse(f, 'fasta'):
        print(r.id, ':', r.seq.lower().count_overlap(motif))

Note: in future I would recommend don't add the motif to the start of your FASTA file, keep it separate so as not to alter the standard file format.

Answer 3

I like @Chris_Rands tidy use of OO, objects are really important per se.

There appears a paradox with your assignment that you might want to refer back to your supervisor/tutors.

If "AT" is the motif, this cannot overlap with another "AT" motif, therefore the request for "overlapping motifs" makes this part of the code superfluous.

It would be better expressed if the motif was "ATA" for example. Thus if the sequence was ATATATA the motif is present 3 times, but only twice if the motif was contiguous.

More generally the code essentially searches for AT microsatellite motifs and formally performed using e.g. tandem repeat finder.