0
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – M__
    Nov 20, 2022 at 0:10
  • $\begingroup$ 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. $\endgroup$
    – terdon
    Nov 21, 2022 at 10:15
  • $\begingroup$ Just to add, this just for Q2. Q1 is a fair question by itself any time and I like the solution. $\endgroup$
    – M__
    Nov 22, 2022 at 0:42
  • $\begingroup$ @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 $\endgroup$
    – thole
    Nov 22, 2022 at 9:27
  • $\begingroup$ @thole post updated as requested, the information sort is below the ------. There was a bug in my implementation of Q1 it now works. $\endgroup$
    – M__
    Nov 22, 2022 at 12:00

1 Answer 1

0
$\begingroup$

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

$\endgroup$

Your Answer

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

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