<p>For the first part is easy,</p>
<pre><code>from Bio import SeqIO
seqs = SeqIO.to_dict(SeqIO.parse(&quot;myfile.fa&quot;, &quot;fasta&quot;))
print (len(seqs))
</code></pre>
<p>The second part is solvable Python's amazing slicing, this is the concept ....</p>
<pre><code>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) 
</code></pre>
<p>To print it:</p>
<pre><code>key, value in mydict.items():
    print(&quot;No. of&quot; + key + '=' + len(mydict[key][value]))
</code></pre>
<p><strong>Note 1</strong>, the second answer is concept only.</p>
<p><strong>Note 2</strong> there can be an issue between the <code>list</code> location and the sequence location ... that is <strong>super</strong> important otherwise it will not work. I've added <code>-1</code> from both locations of the amino acid and the location on the <code>list</code>, <em>but this needs checking</em>. This is because <code>list</code>s start at <code>0</code> whilst amino acids start at <code>1</code></p>
<hr />
<p><strong>Update</strong> <em>Advice on Python</em> requested by OP</p>
<ol>
<li>Answer to Q1 has been debugged, tested and works.</li>
<li>Use Q2 as the basis ... just debug it (see point 4.), it will be a bit fiddly because moving from <code>lists</code> to aa strings.</li>
<li>Use an IDE, I use Spyder: don't get too caught up in &quot;the best IDE&quot; they just make life easier.</li>
<li>Book: Mark Lutz &quot;Learning Python 5th Edition&quot; by O'Reilly section on <code>lists</code> (slicing) and <code>dictionaries</code>  <strong>*</strong>.</li>
<li>Online resources <code>GeekforGeeks</code> on <a href="https://www.geeksforgeeks.org/python-lists/" rel="nofollow noreferrer"><code>lists</code></a> (slicing) and <a href="https://www.geeksforgeeks.org/python-dictionary/" rel="nofollow noreferrer"><code>dictionaries</code></a>.</li>
<li><a href="https://biopython.org/wiki/SeqIO" rel="nofollow noreferrer">BioPython docs on SeqIO</a> search for <code>record.id</code> and <code>record.seq</code> ... you need to understand the SeqIO object.</li>
</ol>
<p><strong>*</strong> This is physically a huge book so an electronic version is strongly preferred just for logistics. This goes beyond the full modern implementation of <code>OOP</code> in bioinformatics (which is not very modern in truth and can become unreadable).</p>

CC BY-SA 4.0
2022-11-19T17:42:08.523Z
2022-11-22T12:09:17.830Z
2022-11-22T12:09:17.830Z
M__
<p><em>This question has also been asked on <a href="https://www.biostars.org/p/9545742/" rel="nofollow noreferrer">Biostars</a></em></p>
<p>I want to find the frequency of each motif-matching peptide.
The sequences are:</p>
<pre><code>&gt; Sequence_1
MPPRRSIVEVKVLDVQKRRVPNKHYVYIIRVTWSSGATEAIYRRYSKFFDLQMQMLDKFP
MEGGQKDPKQRIIPFLPGKILFRRSHIRDVAVKRLIPIDEYCKALIQLPPYISQCDEVLQ
FFETRPEDLNPPKEEHIGKKKSGNDPTSVDPMVLEQYVVVADYQKQESSEISLSVGQVVD
IIEKNESGWWFVSTAEEQGWVPATCLEGQDGVQDEFSLQPEEEEKYTVIYPYTARDQDEM
NLERGAVVEVVQKNLEGWWKIRYQGKEGWAPASYLKKNSGEPLPPKLGPSSPAHSGALDL
DGVSRHQNAMGREKELLNNQRDGRFEGRLVPDGDVKQRSPKMRQRPPPRRDMTIPRGLNL

&gt;Sequence_2
MAEVRKFTKRLSKPGTAAELRQSVSEAVRGSVVLEKAKLVEPLDYENVITQRKTQIYSDP
LRDLLMFPMEDISISVIGRQRRTVQSTVPEDAEKRAQSLFVKECIKTYSTDWHVVNYKYE
DFSGDFRMLPCKSLRPEKIPNHVFEIDEDCEKDEDSSSLCSQKGGVIKQGWLHKANVNST
ITVTMKVFKRRYFYLTQLPDGSYILNSYKDEKNSKESKGCIYLDACIDVVQCPKMRRHAF
ELKMLDKYSHYLAAETEQEMEEWLIMLKKIIQINTDSLVQEKKDTVEAIQEEETSSQGKA
ENIMASLERSMHPELMKYGRETEQLNKLSRGDGRQNLFSFDSEVQRLDFSGIEPDVKPFE
EKCNKRFMVNCHDLTFNILGHIGDNAKGPPTNVEPFFINLALFDVKNNCKISADFHVDLN
PPSVREMLWGTSTQLSNDGNAKGFSPESLIHGIAESQLCYIKQGIFSVTNPHPEIFLVVR

&gt;Sequence_3
GDDSEWLKLPVDQKCEHKLWKARLSGYEEALKIFQKIKDEKSPEWSKYLGLIKKFVTDS
NAVVQLKGLEAALVYVENAHVAGKTTGEVVSGVVSKVFNQPKAKAKELGIEICLMYVEIE
KGESVQEELLKGLDNKNPKIIVACIETLRKALSEFGSKIISLKPIIKVLPKLFESRDKAV
RDEAKLFAIEIYRWNRDAVKHTLQNINSVQLKELEEEWVKLPTGAPKPSRFLRSQQELEA
KLEQQQSAGGDAEGGGDDGDEVPQVDAYELLDAVEILSKLPKDFYDKIEAKKWQERKEAL
EAVEVLVKNPKLEAGDYADLVKALKKVVGKDTNVMLVALAAKCLTGLAVGLRKKFGQYAG
HVVPTILEKFKEKKPQVVQALQEAIDAIFLTTTLQNISEDVLAVMDNKNPTIKQQTSLFI
ARSFRHCTSSTLPKSLLKPFCAALLKHINDSAPEVRDAAFEALGTALKVVGEKSVNPFLA

&gt; Sequence_4
MPPRRSIVEVKVLDVQKRRVPNKHYVYIIRVTWSSGATEAIYRRYSKFFDLQMQMLDKFP
MEGGQKDPKQRIIPFLPGKILFRRSHIRDVAVKRLIPIDEYCKALIQLPPYISQCDEVLQ
FFETRPEDLNPPKEEHIGKKKSGNDPTSVDPMVLEQYVVVADYQKQESSEISLSVGQVVD
IIEKNESGWWFVSTAEEQGWVPATCLEGQDGVQDEFSLQPEEEEKYTVIYPYTARDQDEM
NLERGAVVEVVQKNLEGWWKIRYQGKEGWAPASYLKKNSGEPLPPKLGPSSPAHSGALDL
DGVSRHQNAMGREKELLNNQRDGRFEGRLVPDGDVKQRSPKMRQRPPPRRDMTIPRGLNL 


&gt;Sequence_5
GDDSEWLKLPVDQKCEHKLWKARLSGYEEALKIFQKIKDEKSPEWSKYLGLIKKFVTDS
NAVVQLKGLEAALVYVENAHVAGKTTGEVVSGVVSKVFNQPKAKAKELGIEICLMYVEIE
KGESVQEELLKGLDNKNPKIIVACIETLRKALSEFGSKIISLKPIIKVLPKLFESRDKAV
RDEAKLFAIEIYRWNRDAVKHTLQNINSVQLKELEEEWVKLPTGAPKPSRFLRSQQELEA
KLEQQQSAGGDAEGGGDDGDEVPQVDAYELLDAVEILSKLPKDFYDKIEAKKWQERKEAL
EAVEVLVKNPKLEAGDYADLVKALKKVVGKDTNVMLVALAAKCLTGLAVGLRKKFGQYAG
HVVPTILEKFKEKKPQVVQALQEAIDAIFLTTTLQNISEDVLAVMDNKNPTIKQQTSLFI
ARSFRHCTSSTLPKSLLKPFCAALLKHINDSAPEVRDAAFEALGTALKVVGEKSVNPFLA
</code></pre>
<p>The fasta file actually has more sequences, I just chose the top 5, to put on here.</p>
<p>The code used is below:</p>
<pre><code>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&quot;(P[A-Z]{2}P..)&quot;,sequence)
     for m in matches:
        print(name, m.start(), m.end(), m.group())
</code></pre>
<p>This is giving me my desired output, which is below:</p>
<pre><code>Sequence1 73 79 PFLPGK
Sequence1 281 287 PLPPKL
Sequence1 288 294 PSSPAH
Sequence4 73 79 PFLPGK
Sequence4 281 287 PLPPKL
Sequence4 288 294 PSSPAH
</code></pre>
<p>Now I want to:</p>
<ol>
<li><p>Calculate the number of sequences in that file. I know in linux environment I would have to use <code>grep -c &quot;&gt;&quot; Sequences.fasta</code> to get the total number of sequences in the fasta file, but how do I do it in python?.</p>
</li>
<li><p>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 &amp; 4, so the <strong>frequency</strong> is 2, meaning that two sequences have PFLPGK. How do I find that frequency for each of the peptides using python.</p>
</li>
</ol>
<p>Many thanks for your assistance.</p>

CC BY-SA 4.0
2022-11-20T00:10:50.207Z
Hi @thole the historic protocol on the site was to ask an OP what code has already been written to address in this instance Q 1 and 2. I personally think there should be a grace period. However, please avoid asking the site to write your code base from scratch for future reference.
M__
CC BY-SA 4.0
2022-11-21T10:15:13.343Z
Yes, pelase don't treat this site as a free script writing service. We are happy to _help_ you with your work, but we won't just write it for you like this.
terdon
CC BY-SA 4.0
2022-11-22T00:42:58.613Z
Just to add, this just for Q2. Q1 is a fair question by itself any time and I like the solution.
M__
CC BY-SA 4.0
2022-11-22T09:27:57.903Z
@M__ thank you, I was able to finally figure out how to work on Q1 without using Biopython. As for Q2 I am still working on it. Do you by any chance have any suggestions on books or free sites that I can use to continue learn python? Once again thank you for your suggestions
thole
CC BY-SA 4.0
2022-11-22T12:00:10.633Z
@thole post updated as requested, the information sort is below the `------`. There was a bug in my implementation of `Q1` it now works.
M__
<phylogenetics><proteins><python><motifs><homework>
finding the frequency of a motif-matching peptide
# finding the frequency of a motif-matching peptide

(Asked by: username_0 on 2022-11-19)

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:

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?.

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.

**Tags:** phylogenetics, proteins, python, motifs, homework
> **Comments:**
> **username_1**: Hi @thole the historic protocol on the site was to ask an OP what code has already been written to address in this instance Q 1 and 2. I personally think there should be a grace period. However, please avoid asking the site to write your code base from scratch for future reference.  
> **username_2**: Yes, pelase don't treat this site as a free script writing service. We are happy to _help_ you with your work, but we won't just write it for you like this.  
> **username_1**: Just to add, this just for Q2. Q1 is a fair question by itself any time and I like the solution.  
> **username_0**: @M__ thank you, I was able to finally figure out how to work on Q1 without using Biopython. As for Q2 I am still working on it. Do you by any chance have any suggestions on books or free sites that I can use to continue learn python? Once again thank you for your suggestions  
> **username_1**: @thole post updated as requested, the information sort is below the `------`. There was a bug in my implementation of `Q1` it now works.  

---

## 1 answers

### Answer by username_1 on 2022-11-19. Score: 1 | Accepted answer

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

Answer to Q1 has been debugged, tested and works.
Use Q2 as the basis ... just debug it (see point 4.), it will be a bit fiddly because moving from lists to aa strings.
Use an IDE, I use Spyder: don't get too caught up in "the best IDE" they just make life easier.
Book: Mark Lutz "Learning Python 5th Edition" by O'Reilly section on lists (slicing) and dictionaries *.
Online resources GeekforGeeks on lists (slicing) and dictionaries.
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).


---

