+++ /dev/null
-# Copyright Yair Benita Y.Benita@pharm.uu.nl
-# Biopython (http://biopython.org) license applies
-
-import sys
-import ProtParamData, IsoelectricPoint
-from ProtParamData import kd # Added by Iddo to enable the gravy method
-from Bio.Seq import Seq
-from Bio.Alphabet import IUPAC
-from Bio.Data import IUPACData
-#from BioModule import
-
-class ProteinAnalysis:
- """
- This class contains methods for protein analysis. The class init method takes
- only one argument, the protein sequence as a string and build a sequence
- object using the Bio.Seq module. This is done just to make sure the sequence
- is a protein sequence and not anything else.
-
- methods:
-
- count_amino_acids:
-
- Simply counts the number times an amino acid is repeated in the protein
- sequence. Returns a dictionary {AminoAcid:Number} and also stores the
- dictionary in self.amino_acids_content.
-
- get_amino_acids_percent:
-
- The same as count_amino_acids only returns the Number in percentage of entire
- sequence. Returns a dictionary and stores the dictionary in
- self.amino_acids_content_percent.
-
- molecular_weight:
- Calculates the molecular weight of a protein.
-
- aromaticity:
-
- Calculates the aromaticity value of a protein according to Lobry, 1994. It is
- simply the relative frequency of Phe+Trp+Tyr.
-
-
- instability_index:
-
- Implementation of the method of Guruprasad et al. (Protein Engineering
- 4:155-161,1990). This method tests a protein for stability. Any value above 40
- means the protein is unstable (=has a short half life).
-
- flexibility:
- Implementation of the flexibility method of Vihinen et al. (Proteins. 1994 Jun;19(2):141-9).
-
- isoelectric_point:
- This method uses the module IsoelectricPoint to calculate the pI of a protein.
-
- secondary_structure_fraction:
- This methods returns a list of the fraction of amino acids which tend to be in Helix, Turn or Sheet.
- Amino acids in helix: V, I, Y, F, W, L.
- Amino acids in Turn: N, P, G, S.
- Amino acids in sheet: E, M, A, L.
- The list contains 3 values: [Helix, Turn, Sheet].
-
-
- protein_scale(Scale, WindwonSize, Edge):
-
- An amino acid scale is defined by a numerical value assigned to each type of
- amino acid. The most frequently used scales are the hydrophobicity or
- hydrophilicity scales and the secondary structure conformational parameters
- scales, but many other scales exist which are based on different chemical and
- physical properties of the amino acids. You can set several parameters that
- control the computation of a scale profile, such as the window size and the
- window edge relative weight value. WindowSize: The window size is the length
- of the interval to use for the profile computation. For a window size n, we
- use the i- ( n-1)/2 neighboring residues on each side of residue it compute
- the score for residue i. The score for residue is the sum of the scale values
- for these amino acids, optionally weighted according to their position in the
- window. Edge: The central amino acid of the window always has a weight of 1.
- By default, the amino acids at the remaining window positions have the same
- weight, but you can make the residue at the center of the window have a
- larger weight than the others by setting the edge value for the residues at
- the beginning and end of the interval to a value between 0 and 1. For
- instance, for Edge=0.4 and a window size of 5 the weights will be: 0.4, 0.7,
- 1.0, 0.7, 0.4. The method returns a list of values which can be plotted to
- view the change along a protein sequence. Many scales exist. Just add your
- favorites to the ProtParamData modules.
- """
- def __init__(self, ProtSequence):
- if ProtSequence.islower():
- self.sequence = Seq(ProtSequence.upper(), IUPAC.protein)
- else:
- self.sequence = Seq(ProtSequence, IUPAC.protein)
- self.amino_acids_content = None
- self.amino_acids_percent = None
- self.length = len(self.sequence)
-
- def count_amino_acids(self):
- ProtDic = dict([ (k, 0) for k in IUPACData.protein_letters])
- for i in ProtDic.keys():
- ProtDic[i]=self.sequence.count(i)
- self.amino_acids_content = ProtDic
- return ProtDic
-
- """Calculate the amino acid content in percents.
- input is the dictionary from CountAA.
- output is a dictionary with AA as keys."""
- def get_amino_acids_percent(self):
- if not self.amino_acids_content:
- self.count_amino_acids()
-
- PercentAA = {}
- for i in self.amino_acids_content.keys():
- if self.amino_acids_content[i] > 0:
- PercentAA[i]=self.amino_acids_content[i]/float(self.length)
- else:
- PercentAA[i] = 0
- self.amino_acids_percent = PercentAA
- return PercentAA
-
- # Calculate MW from Protein sequence
- # Calculate MW from Protein sequence
- def molecular_weight (self):
- # make local dictionary for speed
- MwDict = {}
- # remove a molecule of water from the amino acid weight.
- for i in IUPACData.protein_weights.keys():
- MwDict[i] = IUPACData.protein_weights[i] - 18.02
- MW = 18.02 # add just one water molecule for the whole sequence.
- for i in self.sequence:
- MW += MwDict[i]
- return MW
-
- # calculate the aromaticity according to Lobry, 1994.
- # Arom=sum of relative frequency of Phe+Trp+Tyr
- def aromaticity(self):
- if not self.amino_acids_percent:
- self.get_amino_acids_percent()
-
- Arom= self.amino_acids_percent['Y']+self.amino_acids_percent['W']+self.amino_acids_percent['F']
- return Arom
-
- # a function to calculate the instability index according to:
- # Guruprasad K., Reddy B.V.B., Pandit M.W. Protein Engineering 4:155-161(1990).
- def instability_index(self):
- #make the dictionary local for speed.
- DIWV=ProtParamData.DIWV.copy()
- score=0.0
- for i in range(self.length - 1):
- DiPeptide=DIWV[self.sequence[i]][self.sequence[i+1]]
- score += DiPeptide
- return (10.0/self.length) * score
-
- # Calculate the flexibility according to Vihinen, 1994.
- # No argument to change window size because parameters are specific for a window=9.
- # the parameters used are optimized for determining the flexibility.
- def flexibility(self):
- Flex = ProtParamData.Flex.copy()
- Window=9
- Weights=[0.25,0.4375,0.625,0.8125,1]
- List=[]
- for i in range(self.length - Window):
- SubSeq=self.sequence[i:i+Window]
- score = 0.0
- for j in range(Window/2):
- score += (Flex[SubSeq[j]]+Flex[SubSeq[Window-j-1]]) * Weights[j]
- score += Flex[SubSeq[Window/2+1]]
- List.append(score/5.25)
- return List
-
- # calculate the gravy according to kyte and doolittle.
- def gravy(self):
- ProtGravy=0.0
- for i in self.sequence:
- ProtGravy += kd[i]
-
- return ProtGravy/self.length
-
- # this method is used to make a list of relative weight of the
- # window edges compared to the window center. The weights are linear.
- # it actually generates half a list. For a window of size 9 and edge 0.4
- # you get a list of [0.4, 0.55, 0.7, 0.85].
- def _weight_list(self, window, edge):
- unit = ((1.0-edge)/(window-1))*2
- list = [0.0]*(window/2)
- for i in range(window/2):
- list[i] = edge + unit * i
- return list
-
- # this method allows you to compute and represent the profile produced
- # by any amino acid scale on a selected protein.
- # Similar to expasy's ProtScale: http://www.expasy.org/cgi-bin/protscale.pl
- # The weight list returns only one tail. If the list should be [0.4,0.7,1.0,0.7,0.4]
- # what you actually get from _weights_list is [0.4,0.7]. The correct calculation is done
- # in the loop.
- def protein_scale(self, ParamDict, Window, Edge=1.0):
- # generate the weights
- weight = self._weight_list(Window,Edge)
- list = []
- # the score in each Window is divided by the sum of weights
- sum_of_weights = 0.0
- for i in weight: sum_of_weights += i
- # since the weight list is one sided:
- sum_of_weights = sum_of_weights*2+1
-
- for i in range(self.length-Window+1):
- subsequence = self.sequence[i:i+Window]
- score = 0.0
- for j in range(Window/2):
- # walk from the outside of the Window towards the middle.
- # Iddo: try/except clauses added to avoid raising an exception on a non-standad amino acid
- try:
- score += weight[j] * ParamDict[subsequence[j]] + weight[j] * ParamDict[subsequence[Window-j-1]]
- except KeyError:
- sys.stderr.write('warning: %s or %s is not a standard amino acid.\n' %
- (subsequence[j],subsequence[Window-j-1]))
-
- # Now add the middle value, which always has a weight of 1.
- if subsequence[Window/2] in ParamDict:
- score += ParamDict[subsequence[Window/2]]
- else:
- sys.stderr.write('warning: %s is not a standard amino acid.\n' % (subsequence[Window/2]))
-
- list.append(score/sum_of_weights)
- return list
-
- # calculate the isoelectric point.
- def isoelectric_point(self):
- if not self.amino_acids_content:
- self.count_amino_acids()
- X = IsoelectricPoint.IsoelectricPoint(self.sequence, self.amino_acids_content)
- return X.pi()
-
- # calculate fraction of helix, turn and sheet
- def secondary_structure_fraction (self):
- if not self.amino_acids_percent:
- self.get_amino_acids_percent()
- Helix = self.amino_acids_percent['V'] + self.amino_acids_percent['I'] + self.amino_acids_percent['Y'] + self.amino_acids_percent['F'] + self.amino_acids_percent['W'] + self.amino_acids_percent['L']
- Turn = self.amino_acids_percent['N'] + self.amino_acids_percent['P'] + self.amino_acids_percent['G'] + self.amino_acids_percent['S']
- Sheet = self.amino_acids_percent['E'] + self.amino_acids_percent['M'] + self.amino_acids_percent['A'] + self.amino_acids_percent['L']
- return Helix, Turn, Sheet
-
-#---------------------------------------------------------#
-"""
-X = ProteinAnalysis("MAEGEITTFTALTEKFNLPPGNYKKPKLLYCSNGGHFLRILPDGTVDGTRDRSDQHIQLQLSAESVGEVYIKSTETGQYLAMDTSGLLYGSQTPSEECLFLERLEENHYNTYTSKKHAEKNWFVGLKKNGSCKRGPRTHYGQKAILFLPLPV")
-print X.count_amino_acids()
-print X.get_amino_acids_percent()
-print X.molecular_weight()
-print X.aromaticity()
-print X.instability_index()
-print X.flexibility()
-print X.pi()
-print X.secondary_structure_fraction()
-print X.protein_scale(ProtParamData.kd, 9, 0.4)
-"""
git://source.jalview.org / jabaws.git / blobdiff