motif searching in a multi-fasta file

Question

I am having trouble trying to search for motifs in a multi fasta file. I have used two techniques one gives me the name of the sequence where the motif is found but doesn't give me the the motif and its position. The other method does not return anything. Please see the code below, any assistance will be highly appreciated.

This is the protein sequence: (Please note I have only put in the first 3 sequences, the actual file I am working with has 100 sequences)

>P1
MSDNGPQNQRNAPRITFGGPSDSTGSNQNGERSGARSKQRRPQGLPNNTASWFTALTQHG
KEDLKFPRGQGVPINTNSSPDDQIGYYRRATRRIRGGDGKMKDLSPRWYFYYLGTGPEAG
LPYGANKDGIIWVATEGALNTPKDHIGTRNPANNAAIVLQLPQGTTLPKGFYAEGSRGGS
QASSRSSSRSRNSSRNSTPGSSRGISPARMAGNGGDAALALLLLDRLNQLESKMSGKGQQ
QQGQXXTXKSAAEASKKPRQKRTATKAYNVTQAFGRRGPEQTQGNFGDQELIRQGTDYKH
WPQIAQFAPSASAFFGMSRIGMEVTPSGTWLTYTGAIKLDDKDPNFKDQVILLNKHIDAY
KTFPPTEPKKDKKKKADETQALPQRQKKQQTVTLLPAADLDDFSKQLQQSMSSADSTQA*

>P2
MSDNGPQNQRNAPRITFGGPSDSTGSNQNGERSGARSKQRRPQGLPNNTASWFTALTQHG
KEDLKFPRGQGVPINTNSSPDDQIGYYRRATRRIRGGDGKMKDLSPRWYFYYLGTGPEAG
LPYGANKDGIIWVATEGALNTPKDHIGTRNPANNAAIVLQLPQGTTLPKGFYAEGSRGGS
QASSRSSSRSRNSSRNSTPGSSRGISPARMAGNGGDAALALLLLDRLNQLESKMSGKGQQ
QQGQTVTKKSAAEASKKPRQKRTATKAYNVTQAFGRRGPEQTQGNFGDQELIRQGTDYKH
WPQIAQFAPSASAFFGMSRIGMEVTPSGTWLTYTGAIKLDDKDPNFKDQVILLNKHIDAY
KTFPPTEPKKDKKKKADETQALPQRQKKQQTVTLLPAADLDDFSKQLQQSMSSADSTQA*

>P3
MSDNGPQNQRNAPRITFGGPSDSTGSNQNGERSGARSKQRRPQGLPNNTASWFTALTQHG
KEDLKFPRGQGVPINTNSSPDDQIGYYRRATRRIRGGDGKMKDLSPRWYFYYLGTGPEAG
LPYGANKDGIIWVATEGALNTPKDHIGTRNPANNAAIVLQLPQGTTLPKGFYAEGSRGGS
QASSRSSSRSRNSSRNSTPGSSRGISPARMAGNGGDAALALLLLDRLNQLESKMSGKGQQ
QQGQTVTKKSAAEASKKPRQKRTATKAYNVTQAFGRRGPEQTQGNFGDQELIRQGTDYKH
WPQIAQFAPSASAFFGMSRIGMEVTPSGTWLTYTGAIKLDDKDPNFKDQVILLNKHIDAY
KTFPPTEPKKDKKKKADETQALPQRQKKQQTVTLLPAADLDDFSKQLQQSMSSADSTQA*

The first one is:

infile=open("sequence.fasta",'r')
out=open("Result.csv",'w')
pattern=re.compile(r"[YFWLIMVA]..[LMALVN].[AGSTCD]..[LAIVNFYMW]")
for line in infile:
  line = line.strip("\n")
  if line.startswith('>'):
     name=line
else:
   s=re.finditer(pattern,line)
   print('%s:%s' %(name,s))
   out.write('%s:\t%s\n' %(name,s))

This one above returns, the output below. It gives me the name of the sequence which is P100, but it does not give me (1) The motif found (2) the position of the motif. And the Results.csv file that it generates is blank.

#output
>P100:<callable_iterator object at 0x000001D0611B7AF0

I then tried a different technique:

infile=open("sequence.fasta",'r')
open=open("result.csv",'w')
pattern=re.compile(r"(P[A-Z]{2}P")
for line in infile:
  line=line.strip("\n")
  if line.startswith('>'):
     name=line
else:
   s = re.finditer(pattern,line)
   for match_obj in s:
       print(match_obj)
       print(match_obj.group())
  out.write('%s:\t%s\t%s\n' %(name, match_obj,match_obj.group()))

The one above does not return anything at all. the Results.csv file is also blank. I'm still new in this python, and some of the techniques I used here I found on bioinformatics stackexchange.

The output I desire, is where I can get (1) Name of sequence, (2) Motif found and (3) position of motif (where it starts and ends) as shown in the example below in the context of every sequence in the multi fasta file

>P100
motif:
position:

Any format is fine as long as I can get something similar to the format above. Could I also be assistet with the indention as that cofuses me.

Answer 1

First of all, you can't do this this way. Fasta sequences can, and indeed yours do, span multiple lines so in order to find a pattern, you need to join those lines into one to make sure you don't miss any that are found across more than one line.

Next, as a general rule, it is a bad idea to hardcode output file names. That means the user cannot control the output, and that any existing files will be overwritten, so just print to standard output and redirect to a file:

foo.py > out_file

You also don't want to hardcode the name of the input file since that just makes the script harder to use. Instead, read the input file name as an argument. With all this in mind, here is a working version of your script:

#!/usr/bin/env python3
import re
import sys

def print_matches(seq_name, matches):
     if matches:
          for match in matches:
               print(name,match.start(), match.end(),match.group(),sep='\t' )
      

pattern = re.compile(r"[YFWLIMVA]..[LMALVN].[AGSTCD]..[LAIVNFYMW]")
seq = ""

try:
     with open(sys.argv[1],'r') as infile:
          for line in infile:
               line = line.strip("\n")
               if line.startswith('>'):
                    name=line.lstrip('>')
                    if seq != "":
                         matches = re.finditer(pattern, seq)
                         print_matches(name, matches)
                         seq = ""
               else:
                    seq += line
          else:
               matches = re.finditer(pattern, seq)
               print_matches(name, matches)
                    
except IOError as e:
     print("Failed to open file %s with error %s" % (sys.argv[1], e))
     sys.exit(1)

Running that on your example gives:

$ foo.py sequence.fasta 
P2  44  53  LPNNTASWF
P2  122 131 YGANKDGII
P2  209 218 MAGNGGDAA
P2  219 228 ALLLLDRLN
P2  304 313 AQFAPSASA
P2  391 400 VTLLPAADL
P2  410 419 MSSADSTQA
P3  44  53  LPNNTASWF
P3  122 131 YGANKDGII
P3  209 218 MAGNGGDAA
P3  219 228 ALLLLDRLN
P3  304 313 AQFAPSASA
P3  391 400 VTLLPAADL
P3  410 419 MSSADSTQA
P3  44  53  LPNNTASWF
P3  122 131 YGANKDGII
P3  209 218 MAGNGGDAA
P3  219 228 ALLLLDRLN
P3  304 313 AQFAPSASA
P3  391 400 VTLLPAADL
P3  410 419 MSSADSTQA

The output is 4 tab-separated columns per match with the sequence name, the match's start position, its end position and the actual sequence matched.

Answer 2

Read the sequences with Biopython, reverse translate the sequences and the pattern peptides, construct pattern regexes, search (find_matches()) for the patterns for each sequence in a loop using DNA Chisel.

Disclaimer: I'm the maintainer of DNA Chisel