finding the frequency of a motif-matching peptide

Question

This question has also been asked on Biostars

I want to find the frequency of each motif-matching peptide. The sequences are:

> Sequence_1
MPPRRSIVEVKVLDVQKRRVPNKHYVYIIRVTWSSGATEAIYRRYSKFFDLQMQMLDKFP
MEGGQKDPKQRIIPFLPGKILFRRSHIRDVAVKRLIPIDEYCKALIQLPPYISQCDEVLQ
FFETRPEDLNPPKEEHIGKKKSGNDPTSVDPMVLEQYVVVADYQKQESSEISLSVGQVVD
IIEKNESGWWFVSTAEEQGWVPATCLEGQDGVQDEFSLQPEEEEKYTVIYPYTARDQDEM
NLERGAVVEVVQKNLEGWWKIRYQGKEGWAPASYLKKNSGEPLPPKLGPSSPAHSGALDL
DGVSRHQNAMGREKELLNNQRDGRFEGRLVPDGDVKQRSPKMRQRPPPRRDMTIPRGLNL

>Sequence_2
MAEVRKFTKRLSKPGTAAELRQSVSEAVRGSVVLEKAKLVEPLDYENVITQRKTQIYSDP
LRDLLMFPMEDISISVIGRQRRTVQSTVPEDAEKRAQSLFVKECIKTYSTDWHVVNYKYE
DFSGDFRMLPCKSLRPEKIPNHVFEIDEDCEKDEDSSSLCSQKGGVIKQGWLHKANVNST
ITVTMKVFKRRYFYLTQLPDGSYILNSYKDEKNSKESKGCIYLDACIDVVQCPKMRRHAF
ELKMLDKYSHYLAAETEQEMEEWLIMLKKIIQINTDSLVQEKKDTVEAIQEEETSSQGKA
ENIMASLERSMHPELMKYGRETEQLNKLSRGDGRQNLFSFDSEVQRLDFSGIEPDVKPFE
EKCNKRFMVNCHDLTFNILGHIGDNAKGPPTNVEPFFINLALFDVKNNCKISADFHVDLN
PPSVREMLWGTSTQLSNDGNAKGFSPESLIHGIAESQLCYIKQGIFSVTNPHPEIFLVVR

>Sequence_3
GDDSEWLKLPVDQKCEHKLWKARLSGYEEALKIFQKIKDEKSPEWSKYLGLIKKFVTDS
NAVVQLKGLEAALVYVENAHVAGKTTGEVVSGVVSKVFNQPKAKAKELGIEICLMYVEIE
KGESVQEELLKGLDNKNPKIIVACIETLRKALSEFGSKIISLKPIIKVLPKLFESRDKAV
RDEAKLFAIEIYRWNRDAVKHTLQNINSVQLKELEEEWVKLPTGAPKPSRFLRSQQELEA
KLEQQQSAGGDAEGGGDDGDEVPQVDAYELLDAVEILSKLPKDFYDKIEAKKWQERKEAL
EAVEVLVKNPKLEAGDYADLVKALKKVVGKDTNVMLVALAAKCLTGLAVGLRKKFGQYAG
HVVPTILEKFKEKKPQVVQALQEAIDAIFLTTTLQNISEDVLAVMDNKNPTIKQQTSLFI
ARSFRHCTSSTLPKSLLKPFCAALLKHINDSAPEVRDAAFEALGTALKVVGEKSVNPFLA

> Sequence_4
MPPRRSIVEVKVLDVQKRRVPNKHYVYIIRVTWSSGATEAIYRRYSKFFDLQMQMLDKFP
MEGGQKDPKQRIIPFLPGKILFRRSHIRDVAVKRLIPIDEYCKALIQLPPYISQCDEVLQ
FFETRPEDLNPPKEEHIGKKKSGNDPTSVDPMVLEQYVVVADYQKQESSEISLSVGQVVD
IIEKNESGWWFVSTAEEQGWVPATCLEGQDGVQDEFSLQPEEEEKYTVIYPYTARDQDEM
NLERGAVVEVVQKNLEGWWKIRYQGKEGWAPASYLKKNSGEPLPPKLGPSSPAHSGALDL
DGVSRHQNAMGREKELLNNQRDGRFEGRLVPDGDVKQRSPKMRQRPPPRRDMTIPRGLNL 


>Sequence_5
GDDSEWLKLPVDQKCEHKLWKARLSGYEEALKIFQKIKDEKSPEWSKYLGLIKKFVTDS
NAVVQLKGLEAALVYVENAHVAGKTTGEVVSGVVSKVFNQPKAKAKELGIEICLMYVEIE
KGESVQEELLKGLDNKNPKIIVACIETLRKALSEFGSKIISLKPIIKVLPKLFESRDKAV
RDEAKLFAIEIYRWNRDAVKHTLQNINSVQLKELEEEWVKLPTGAPKPSRFLRSQQELEA
KLEQQQSAGGDAEGGGDDGDEVPQVDAYELLDAVEILSKLPKDFYDKIEAKKWQERKEAL
EAVEVLVKNPKLEAGDYADLVKALKKVVGKDTNVMLVALAAKCLTGLAVGLRKKFGQYAG
HVVPTILEKFKEKKPQVVQALQEAIDAIFLTTTLQNISEDVLAVMDNKNPTIKQQTSLFI
ARSFRHCTSSTLPKSLLKPFCAALLKHINDSAPEVRDAAFEALGTALKVVGEKSVNPFLA

The fasta file actually has more sequences, I just chose the top 5, to put on here.

The code used is below:

import Bio
import regex

from Bio import SeqIO
input_file = 'Sequences.fasta'
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
for fasta in fasta_sequences:
     name, sequence = fasta.id, str(fasta.seq)
     matches=regex.finditer(r"(P[A-Z]{2}P..)",sequence)
     for m in matches:
        print(name, m.start(), m.end(), m.group())

This is giving me my desired output, which is below:

Sequence1 73 79 PFLPGK
Sequence1 281 287 PLPPKL
Sequence1 288 294 PSSPAH
Sequence4 73 79 PFLPGK
Sequence4 281 287 PLPPKL
Sequence4 288 294 PSSPAH

Now I want to:

1. Calculate the number of sequences in that file. I know in linux environment I would have to use grep -c ">" Sequences.fasta to get the total number of sequences in the fasta file, but how do I do it in python?.

2. Once I know how many sequences are in file in total, I want to count how many of those sequences have the motif-matching peptides above. For example how many sequences have PFLPGK at positions 73 79? In the example above PFLPGK its in sequence 1 & 4, so the frequency is 2, meaning that two sequences have PFLPGK. How do I find that frequency for each of the peptides using python.

Many thanks for your assistance.

Answer 1

For the first part is easy,

from Bio import SeqIO
seqs = SeqIO.to_dict(SeqIO.parse("myfile.fa", "fasta"))
print (len(seqs))

The second part is solvable Python's amazing slicing, this is the concept ....

from Bio import SeqIO
loci = {'PFLPGK':[1,10], 'PFLPGK':[50:60] }   
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
mydict = {}
for key, values in loci.items():
    myrange = '-'.join(values)
    mydict = {key:myrange} # you need to initiate the dictionary without defaultdict
    for record in fasta_sequences:
        if record.seq[values[0]-1:values[1]-1] == key:
             mydict[key][myrange].append(record.id) 

To print it:

key, value in mydict.items():
    print("No. of" + key + '=' + len(mydict[key][value]))

Note 1, the second answer is concept only.

Note 2 there can be an issue between the list location and the sequence location ... that is super important otherwise it will not work. I've added -1 from both locations of the amino acid and the location on the list, but this needs checking. This is because lists start at 0 whilst amino acids start at 1

---

Update Advice on Python requested by OP

1. Answer to Q1 has been debugged, tested and works.
2. Use Q2 as the basis ... just debug it (see point 4.), it will be a bit fiddly because moving from lists to aa strings.
3. Use an IDE, I use Spyder: don't get too caught up in "the best IDE" they just make life easier.
4. Book: Mark Lutz "Learning Python 5th Edition" by O'Reilly section on lists (slicing) and dictionaries *.
5. Online resources GeekforGeeks on lists (slicing) and dictionaries.
6. BioPython docs on SeqIO search for record.id and record.seq ... you need to understand the SeqIO object.

* This is physically a huge book so an electronic version is strongly preferred just for logistics. This goes beyond the full modern implementation of OOP in bioinformatics (which is not very modern in truth and can become unreadable).