This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Biopython - Calculation Used In Motif Building

I was trying to build a motif using Biopython 1.58 and I tried to cross-check the amino acid frequency from one column out of curiosity. Why do the values not correspond to the 'exact' value? e.g. given for 1st column there are amino acids: N,S,P,S,V,P,S,N,I,I,I,H,V,V,L. The output of the motif object is: {... 'N': 0.1279761904761905, 'P': 0.1279761904761905, 'S': 0.19047619047619052, ...}. Shoudn't 'N' be 2/15 = 0.06667, 'P' be 2/15 = 0.1333, and 'S' be 3/15 = 0.2? May be I misunderstood, or is there some other transformation done to the amino acid frequency used to build the sequence motif? Grateful if someone could clarify.

Thanks

biopython motif
• 5,040 views • link
updated by Neilfws Neilfws • written by Olivier Olivier

According to your example, N should also be 2/15 and indeed has the same value as P. We also need more details such as: which Biopython module was used?

• 0 views • link
written by Neilfws Neilfws

You're right. My mistake. Thanks. I used the Bio.Motif module.

• 0 views • link
written by Olivier Olivier

Further to Neil's helpful comment, could you actually show us the code you are using with Biopython to get these numbers?

• 0 views • link
written by Peter Peter

Here's the script. Still learning..

from Bio.Motif import Motif
from Bio.Alphabet import IUPAC, Gapped
from Bio.Align.Applications import MuscleCommandline
from Bio import SeqIO, AlignIO, Entrez
import os

#preparing script for communicating with NCBI's Entrez 
os.environ['http_proxy'] = 'http://130.14.29.110:80'
Entrez.email = 'myemailaddress'

#searching for wrky sequence
search_handle = Entrez.esearch(term='wrky AND (tomato[orgn] OR potato[orgn] OR pepper[orgn] OR arabidopsis[orgn] OR)',
                        db='protein')

id_list = Entrez.read(search_handle)['IdList'][:15]

#Retrieving sequences
gb_handle = Entrez.efetch(db='protein', id=id_list, rettype='gb')

fasta_records = SeqIO.convert(gb_handle, 'gb', '15_wrky.fas', 'fasta')

#Aligning sequences
muscle_cline = MuscleCommandline(input='/home/olivier/15_wrky.fas',
                             out='/home/olivier/15_wrky.aln',
                             clwstrict=True)
stdout, stderr = muscle_cline()

#Loading the msa and building a sequence motif
alphabet = Gapped(IUPAC.protein)
alignment_file = AlignIO.read('/home/olivier/15_wrky.aln','clustal',
                              alphabet=alphabet)
m = Motif(alphabet=alphabet)

for sequence in alignment_file:
    m.add_instance(sequence.seq[367:430])

print(m.pwm())
• 0 views • link
updated by Neilfws Neilfws • written by Olivier Olivier

I edited this by indenting lines of code with 4 spaces (to display code properly).

• 0 views • link
written by Neilfws Neilfws

1 answer

is there some other transformation done to the amino acid frequency used to build the sequence motif

I believe that the answer is "yes" and a clue as to the transformation is found in the source code for the pwm method:

if laplace=True (default), pseudocounts equal to self.background multiplied by self.beta are added to all positions.

• 1 views • link
written by Neilfws Neilfws

That text is in the docstring, so doing help(m.pwm) or help(m) would have told you this (where m is a Motif object as in the example).

• 0 views • link
written by Peter Peter

Thanks Peter. I had not thought of it.

• 0 views • link
written by Olivier Olivier

If you find docstrings missing important information please let us know or suggest improvements. Thanks!

• 0 views • link
written by Peter Peter

i.e. If you want raw counts, use:

m.pwm(laplace=False)
• 0 views • link
updated by Neilfws Neilfws • written by Peter Peter

Thanks a lot Peter for clarifying. I'm still going through the Biopython tutorial as a preparation for my research project. Hope that my questions are not too basic, as I'm from a molecular biology background.

• 0 views • link
written by Olivier Olivier

Log in to answer this question.

More posts like this