1
$\begingroup$

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)
$\endgroup$
4
  • 3
    $\begingroup$ Hello Dari, this looks a bit like a homework. It's ok to ask, but it's nice to disclaim it. $\endgroup$ May 8, 2019 at 19:06
  • $\begingroup$ Maybe it's just an error when copying your actual code in your question, but in motif_count, the while loop does not contain the test that can trigger the return, nor the incrementation of count. exit(1) should also probably be re-indented. I tried to fix this. Hopefully this corresponds to your actual code now. $\endgroup$
    – bli
    May 9, 2019 at 13:10
  • $\begingroup$ Your code does nothing with the values returned by motif_count. To create a new dictionary with values modified by a function, the valmap function from the cytoolz module could be useful. $\endgroup$
    – bli
    May 9, 2019 at 13:16
  • $\begingroup$ @bli please don't correct code in questions! We can't know if an error is because the OP didn't copy/paste correctly or if that's the actual code being used, so when you fix errors in the question, it could end up causing everyone to waste their time trying to find the error in code that is now correct. $\endgroup$
    – terdon
    May 9, 2019 at 13:19

3 Answers 3

1
$\begingroup$

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}")
$\endgroup$
1
  • 1
    $\begingroup$ The lack of "exit (1)" indentation was an important bug $\endgroup$
    – M__
    May 9, 2019 at 21:18
3
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ The OP seems to want to read the input from stdin. Using a file is a far superior approach, of course, but maybe stdin is a requirement. $\endgroup$
    – terdon
    May 8, 2019 at 20:22
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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