|
|
| (12 intermediate revisions by 3 users not shown) |
| Line 1: |
Line 1: |
| | {{Infobox script-repo |
| | |type = script |
| | |filename = aaindex.py |
| | |author = [[User:Speleo3|Thomas Holder]] |
| | |license = BSD |
| | }} |
| | |
| | {{Infobox psico |
| | |module = psico.aaindex |
| | }} |
| | |
| | == Introduction == |
| | [[File:AAindexExample.png|200px|thumb|right|Hydrophobicity coloring with KYTJ820101]] |
| | |
| AAindex is a database of numerical indices representing various physicochemical and biochemical properties of amino acids and pairs of amino acids. See http://www.genome.jp/aaindex/ | | AAindex is a database of numerical indices representing various physicochemical and biochemical properties of amino acids and pairs of amino acids. See http://www.genome.jp/aaindex/ |
|
| |
|
| This script is a python parser for the AAindex flat files which you have to download: | | This script is a python parser for the AAindex flat files which will be downloaded from ftp://ftp.genome.jp/pub/db/community/aaindex/ |
| | |
| * ftp://ftp.genome.jp/pub/db/community/aaindex/aaindex1
| |
| * ftp://ftp.genome.jp/pub/db/community/aaindex/aaindex2
| |
| * ftp://ftp.genome.jp/pub/db/community/aaindex/aaindex3
| |
|
| |
|
| The script provides two PyMOL commands (but can also be used without PyMOL). | | The script provides two PyMOL commands (but can also be used without PyMOL). |
| Line 12: |
Line 22: |
| * pmf: Potential of Mean Force (aaindex3) | | * pmf: Potential of Mean Force (aaindex3) |
|
| |
|
| = Python Example = | | == Python Example == |
|
| |
|
| Consider the script is called aaindex.py, it is placed somewhere in your PYTHONPATH and the aaindex flatfiles are found in the current directory. | | Consider the script is called aaindex.py, it is placed somewhere in your PYTHONPATH and the aaindex flatfiles are found in the current directory. |
| Line 22: |
Line 32: |
| aaindex.grep('volume') | | aaindex.grep('volume') |
| x = aaindex.get('KRIW790103') | | x = aaindex.get('KRIW790103') |
| print x | | print(x) |
| print x.get('A') | | print(x.get('A')) |
|
| |
|
| | aaindex.init(index='2') |
| aaindex.grep('blosum') | | aaindex.grep('blosum') |
| x = aaindex.get('HENS920102') | | x = aaindex.get('HENS920102') |
| print x.get('A', 'K') | | print(x.get('A', 'K')) |
| </source> | | </source> |
|
| |
|
| = PyMOL Example = | | == PyMOL Example == |
| | |
| | Solvent-accessible surface coloring by [https://www.genome.jp/dbget-bin/www_bget?aaindex:KYTJ820101 Hydropathy index (Kyte-Doolittle, 1982)] |
|
| |
|
| <source lang="python"> | | <source lang="python"> |
| import aaindex
| | run aaindex.py |
| aaindex.init(path='.')
| |
| | |
| aaindex2b KYTJ820101 | | aaindex2b KYTJ820101 |
| spectrum b, yellow_white_blue | | spectrum b, white yellow forest |
| | set surface_solvent |
| show surface | | show surface |
| </source> | | </source> |
|
| |
|
| = The Script = | | == See Also == |
| <source lang="python">
| |
| '''
| |
| (c) 2010 Thomas Holder
| |
| | |
| Python parser for AAindex: Amino Acid Index Database
| |
| http://www.genome.jp/aaindex/
| |
| | |
| PyMOL commands:
| |
| | |
| aaindex2b
| |
| pmf
| |
| '''
| |
| | |
| _aaindex = dict()
| |
| _pymol_auto_arg_update = lambda: None
| |
| | |
| def search(pattern, searchtitle=True, casesensitive=False):
| |
| '''
| |
| Search for pattern in description and title (optional) of all records and
| |
| return matched records as list. By default search case insensitive.
| |
| '''
| |
| whatcase = lambda i: i
| |
| if not casesensitive:
| |
| pattern = pattern.lower()
| |
| whatcase = lambda i: i.lower()
| |
| matches = []
| |
| for record in _aaindex.itervalues():
| |
| if pattern in whatcase(record.desc) or searchtitle and pattern in whatcase(record.title):
| |
| matches.append(record)
| |
| return matches
| |
| | |
| def grep(pattern):
| |
| '''
| |
| Search for pattern in title and description of all records (case
| |
| insensitive) and print results on standard output.
| |
| '''
| |
| for record in search(pattern):
| |
| print record
| |
| | |
| class Record:
| |
| '''
| |
| Amino acid index (AAindex) Record
| |
| '''
| |
| aakeys = 'ARNDCQEGHILKMFPSTWYV'
| |
| def __init__(self):
| |
| self.key = None
| |
| self.desc = ''
| |
| self.ref = ''
| |
| self.authors = ''
| |
| self.title = ''
| |
| self.journal = ''
| |
| self.correlated = dict()
| |
| self.index = dict()
| |
| self.comment = ''
| |
| def extend(self, row):
| |
| i = len(self.index)
| |
| for x in row:
| |
| self.index[self.aakeys[i]] = x
| |
| i += 1
| |
| def get(self, aai, aaj=None, d=None):
| |
| assert aaj is None
| |
| return self.index.get(aai, d)
| |
| def median(self):
| |
| x = sorted(filter(None, self.index.values()))
| |
| half = len(x)/2
| |
| if len(x) % 2 == 1:
| |
| return x[half]
| |
| return (x[half-1] + x[half])/2.0
| |
| def __str__(self):
| |
| desc = self.desc.replace('\n', ' ').strip()
| |
| return '%s(%s: %s)' % (self.__class__.__name__, self.key, desc)
| |
| | |
| class MatrixRecord(Record):
| |
| '''
| |
| Matrix record for mutation matrices or pair-wise contact potentials
| |
| '''
| |
| def __init__(self):
| |
| Record.__init__(self)
| |
| self.index = []
| |
| self.rows = dict()
| |
| self.cols = dict()
| |
| def extend(self, row):
| |
| self.index.append(row)
| |
| def _get(self, aai, aaj):
| |
| i = self.rows[aai]
| |
| j = self.cols[aaj]
| |
| return self.index[i][j]
| |
| def get(self, aai, aaj, d=None):
| |
| try:
| |
| return self._get(aai, aaj)
| |
| except:
| |
| pass
| |
| try:
| |
| return self._get(aaj, aai)
| |
| except:
| |
| return d
| |
| def median(self):
| |
| x = []
| |
| for y in self.index:
| |
| x.extend(filter(None, y))
| |
| x.sort()
| |
| if len(x) % 2 == 1:
| |
| return x[len(x)/2]
| |
| return sum(x[len(x)/2-1:len(x)/2+1])/2.0
| |
| | |
| def get(key):
| |
| '''
| |
| Get record for key
| |
| '''
| |
| if len(_aaindex) == 0:
| |
| init()
| |
| return _aaindex.get(key)
| |
| | |
| def _float_or_None(x):
| |
| if x == 'NA' or x == '-':
| |
| return None
| |
| return float(x)
| |
| | |
| def init(path=None, index='123'):
| |
| '''
| |
| Read in the aaindex files. You need to run this (once) before you can
| |
| access any records. If the files are not within the current directory,
| |
| you need to specify the correct directory path. By default all three
| |
| aaindex files are read in.
| |
| '''
| |
| if path is None:
| |
| import sys
| |
| from os.path import sep
| |
| path = __file__.rsplit(sep, 1)[0]
| |
| print >> sys.stderr, 'path =', path
| |
| index = str(index)
| |
| if '1' in index:
| |
| _parse(path + '/aaindex1', Record)
| |
| if '2' in index:
| |
| _parse(path + '/aaindex2', MatrixRecord)
| |
| if '3' in index:
| |
| _parse(path + '/aaindex3', MatrixRecord)
| |
| _pymol_auto_arg_update()
| |
| | |
| def init_from_file(filename, type=Record):
| |
| _parse(filename, type)
| |
| | |
| def _parse(filename, rec, quiet=True):
| |
| '''
| |
| Parse aaindex input file. `rec` must be `Record` for aaindex1 and
| |
| `MarixRecord` for aaindex2 and aaindex3.
| |
| '''
| |
| try:
| |
| f = file(filename, 'rU')
| |
| except:
| |
| print ('Cannot open "%s", please download it from ' + \
| |
| 'ftp://ftp.genome.jp/pub/db/community/aaindex/') % filename
| |
| return
| |
| | |
| current = rec()
| |
| lastkey = None
| |
| | |
| for line in f:
| |
| key = line[0:2]
| |
| if key[0] == ' ':
| |
| key = lastkey
| |
| | |
| if key == '//':
| |
| _aaindex[current.key] = current
| |
| current = rec()
| |
| elif key == 'H ':
| |
| current.key = line[2:].strip()
| |
| elif key == 'R ':
| |
| current.ref += line[2:]
| |
| elif key == 'D ':
| |
| current.desc += line[2:]
| |
| elif key == 'A ':
| |
| current.authors += line[2:]
| |
| elif key == 'T ':
| |
| current.title += line[2:]
| |
| elif key == 'J ':
| |
| current.journal += line[2:]
| |
| elif key == '* ':
| |
| current.comment += line[2:]
| |
| elif key == 'C ':
| |
| a = line[2:].split()
| |
| for i in range(0, len(a), 2):
| |
| current.correlated[a[i]] = float(a[i+1])
| |
| elif key == 'I ':
| |
| a = line[1:].split()
| |
| if a[0] != 'A/L':
| |
| current.extend(map(_float_or_None, a))
| |
| elif list(Record.aakeys) != [i[0] for i in a] + [i[-1] for i in a]:
| |
| print 'Warning: wrong amino acid sequence for', current.key
| |
| else:
| |
| try:
| |
| assert list(Record.aakeys[:10]) == [i[0] for i in a]
| |
| assert list(Record.aakeys[10:]) == [i[2] for i in a]
| |
| except:
| |
| print 'Warning: wrong amino acid sequence for', current.key
| |
| elif key =='M ':
| |
| a = line[2:].split()
| |
| if a[0] == 'rows':
| |
| if a[4] == 'rows':
| |
| a.pop(4)
| |
| assert a[3] == 'cols' and len(a) == 6
| |
| i = 0
| |
| for aa in a[2]:
| |
| current.rows[aa] = i
| |
| i += 1
| |
| i = 0
| |
| for aa in a[5]:
| |
| current.cols[aa] = i
| |
| i += 1
| |
| else:
| |
| current.extend(map(_float_or_None, a))
| |
| elif not quiet:
| |
| print 'Warning: line starts with "%s"' % (key)
| |
| | |
| lastkey = key
| |
| | |
| ########## PYMOL ###########
| |
| | |
| # from Bio.SCOP.Raf import to_one_letter_code
| |
| # See also http://www.pymolwiki.org/index.php/Aa_codes
| |
| to_one_letter_code = {'PAQ': 'Y', 'AGM': 'R', 'ILE': 'I', 'PR3': 'C',
| |
| 'GLN': 'Q', 'DVA': 'V', 'CCS': 'C', 'ACL': 'R', 'GLX': 'Z', 'GLY': 'G',
| |
| 'GLZ': 'G', 'DTH': 'T', 'OAS': 'S', 'C6C': 'C', 'NEM': 'H', 'DLY': 'K',
| |
| 'MIS': 'S', 'SMC': 'C', 'GLU': 'E', 'NEP': 'H', 'BCS': 'C', 'ASQ': 'D',
| |
| 'ASP': 'D', 'SCY': 'C', 'SER': 'S', 'LYS': 'K', 'SAC': 'S', 'PRO': 'P',
| |
| 'ASX': 'B', 'DGN': 'Q', 'DGL': 'E', 'MHS': 'H', 'ASB': 'D', 'ASA': 'D',
| |
| 'NLE': 'L', 'DCY': 'C', 'ASK': 'D', 'GGL': 'E', 'STY': 'Y', 'SEL': 'S',
| |
| 'CGU': 'E', 'ASN': 'N', 'ASL': 'D', 'LTR': 'W', 'DAR': 'R', 'VAL': 'V',
| |
| 'CHG': 'A', 'TPO': 'T', 'CLE': 'L', 'GMA': 'E', 'HAC': 'A', 'AYA': 'A',
| |
| 'THR': 'T', 'TIH': 'A', 'SVA': 'S', 'MVA': 'V', 'SAR': 'G', 'LYZ': 'K',
| |
| 'BNN': 'A', '5HP': 'E', 'IIL': 'I', 'SHR': 'K', 'HAR': 'R', 'FME': 'M',
| |
| 'PYX': 'C', 'ALO': 'T', 'PHI': 'F', 'ALM': 'A', 'PHL': 'F', 'MEN': 'N',
| |
| 'TPQ': 'A', 'GSC': 'G', 'PHE': 'F', 'ALA': 'A', 'MAA': 'A', 'MET': 'M',
| |
| 'UNK': 'X', 'LEU': 'L', 'ALY': 'K', 'SET': 'S', 'GL3': 'G', 'TRG': 'K',
| |
| 'CXM': 'M', 'TYR': 'Y', 'SCS': 'C', 'DIL': 'I', 'TYQ': 'Y', '3AH': 'H',
| |
| 'DPR': 'P', 'PRR': 'A', 'CME': 'C', 'IYR': 'Y', 'CY1': 'C', 'TYY': 'Y',
| |
| 'HYP': 'P', 'DTY': 'Y', '2AS': 'D', 'DTR': 'W', 'FLA': 'A', 'DPN': 'F',
| |
| 'DIV': 'V', 'PCA': 'E', 'MSE': 'M', 'MSA': 'G', 'AIB': 'A', 'CYS': 'C',
| |
| 'NLP': 'L', 'CYQ': 'C', 'HIS': 'H', 'DLE': 'L', 'CEA': 'C', 'DAL': 'A',
| |
| 'LLP': 'K', 'DAH': 'F', 'HMR': 'R', 'TRO': 'W', 'HIC': 'H', 'CYG': 'C',
| |
| 'BMT': 'T', 'DAS': 'D', 'TYB': 'Y', 'BUC': 'C', 'PEC': 'C', 'BUG': 'L',
| |
| 'CYM': 'C', 'NLN': 'L', 'CY3': 'C', 'HIP': 'H', 'CSO': 'C', 'TPL': 'W',
| |
| 'LYM': 'K', 'DHI': 'H', 'MLE': 'L', 'CSD': 'A', 'HPQ': 'F', 'MPQ': 'G',
| |
| 'LLY': 'K', 'DHA': 'A', 'DSN': 'S', 'SOC': 'C', 'CSX': 'C', 'OMT': 'M',
| |
| 'DSP': 'D', 'PTR': 'Y', 'TRP': 'W', 'CSW': 'C', 'EFC': 'C', 'CSP': 'C',
| |
| 'CSS': 'C', 'SCH': 'C', 'OCS': 'C', 'NMC': 'G', 'SEP': 'S', 'BHD': 'D',
| |
| 'KCX': 'K', 'SHC': 'C', 'C5C': 'C', 'HTR': 'W', 'ARG': 'R', 'TYS': 'Y',
| |
| 'ARM': 'R', 'DNP': 'A'}
| |
| | |
| def aaindex2b(key='KYTJ820101', selection='(all)', quiet=0, var='b'):
| |
| '''
| |
| DESCRIPTION
| |
| | |
| "aaindex" looks up the Amino Acid Index from
| |
| http://www.genome.jp/aaindex/
| |
| for the given key and assignes b-factors to the given selection. Unknown
| |
| residues get the average index value assigned.
| |
| | |
| USAGE
| |
| | |
| aaindex2b [key [, selection]]
| |
| | |
| ARGUMENTS
| |
| | |
| key = string: Key of AAindex entry
| |
| | |
| selection = string: atoms to assign b-factors {default: (all)}
| |
| | |
| EXAMPLE
| |
|
| |
|
| # Hydropathy index by Kyte-Doolittle
| | * Protscale from [http://www.rubor.de/pymol_extensions_de.html rTools] does a similar job in coloring by amino acid properties |
| aaindex2b KYTJ820101
| |
| spectrum b, yellow_white_blue
| |
| show surface
| |
| '''
| |
| from pymol import cmd, stored
| |
| entry = get(key)
| |
| median = entry.median()
| |
|
| |
|
| if int(quiet) != 0:
| | [[Category:Script_Library]] |
| print entry.desc.strip()
| | [[Category:Biochemical_Scripts]] |
| | | [[Category:Structural_Biology_Scripts]] |
| def lookup(resn):
| | [[Category:Pymol-script-repo]] |
| one_letter = to_one_letter_code.get(resn, 'X')
| |
| value = entry.get(one_letter)
| |
| if value is None:
| |
| return median
| |
| return value
| |
| stored.aaindex = lookup
| |
| | |
| cmd.alter(selection, var + '=stored.aaindex(resn)')
| |
| | |
| def pmf(key, cutoff=7.0, selection1='(name CB)', selection2='', state=1, quiet=1):
| |
| '''
| |
| DESCRIPTION
| |
| | |
| Potential of Mean Force
| |
| | |
| ARGUMENTS
| |
| | |
| key = string: aaindex key
| |
| | |
| cutoff = float: distance cutoff {default: 7.0}
| |
| cutoff = (float, float): distance shell
| |
| | |
| selection1 = string: atom selection {default: (name CB)}
| |
| | |
| selection2 = string: atom selection {default: selection1}
| |
| | |
| NOTES
| |
| | |
| Does also support a list of keys and a list of cutoffs to deal with
| |
| multiple distance shells.
| |
| | |
| EXAMPLES
| |
| | |
| # load databases
| |
| aaindex.init('/path/to/aaindex/directory')
| |
| | |
| # distance dependent c-beta contact potentials
| |
| pmf SIMK990101, 5, /2x19//A//CB
| |
| pmf SIMK990102, [5, 7.5], /2x19//A//CB
| |
| pmf [SIMK990101, SIMK990102, SIMK990103], [0, 5, 7.5, 10], /2x19//A//CB
| |
| | |
| # interface potential
| |
| pseudo_sidechaincenters 2x19_scc, 2x19
| |
| pmf KESO980102, 7.0, /2x19_scc//A, /2x19_scc//B
| |
| distance /2x19_scc//A, /2x19_scc//B, cutoff=7.0
| |
| '''
| |
| from pymol import cmd, stored
| |
| from chempy import cpv
| |
| if cmd.is_string(key):
| |
| if key.lstrip().startswith('['):
| |
| key = cmd.safe_alpha_list_eval(key)
| |
| else:
| |
| key = [key]
| |
| if cmd.is_string(cutoff):
| |
| cutoff = eval(cutoff)
| |
| if not cmd.is_sequence(cutoff):
| |
| cutoff = [cutoff]
| |
| if len(cutoff) == len(key):
| |
| cutoff = [0.0] + list(cutoff)
| |
| if len(cutoff) != len(key) + 1:
| |
| print 'Error: Number of keys and number of cutoffs inconsistent'
| |
| return
| |
| state = int(state)
| |
| quiet = int(quiet)
| |
| if len(selection2) == 0:
| |
| selection2 = selection1
| |
| if not quiet and len(key) > 1:
| |
| print 'Distance shells:'
| |
| for i in range(len(key)):
| |
| print '%s %.1f-%.1f' % (key[i], cutoff[i], cutoff[i+1])
| |
| | |
| stored.idmap = dict()
| |
| cmd.iterate_state(state, '(%s) or (%s)' % (selection1, selection2),
| |
| 'stored.idmap[model,index] = [(resn,name),(x,y,z)]')
| |
| twoN = cmd.count_atoms(selection1) + cmd.count_atoms(selection2)
| |
| pairs = cmd.find_pairs(selection1, selection2, cutoff=max(cutoff),
| |
| state1=state, state2=state)
| |
| if len(pairs) == 0:
| |
| print 'Empty pair list'
| |
| return 0.0
| |
| | |
| matrix = map(get, key)
| |
| for i in matrix:
| |
| assert isinstance(i, MatrixRecord)
| |
| | |
| i_list = range(len(key))
| |
| u_sum = 0
| |
| count = 0
| |
| for id1, id2 in pairs:
| |
| a1 = stored.idmap[id1]
| |
| a2 = stored.idmap[id2]
| |
| r = cpv.distance(a1[1], a2[1])
| |
| for i in i_list:
| |
| if cutoff[i] <= r and r < cutoff[i+1]:
| |
| try:
| |
| aa1 = to_one_letter_code[a1[0][0]]
| |
| aa2 = to_one_letter_code[a2[0][0]]
| |
| u_sum += matrix[i].get(aa1, aa2)
| |
| count += 1
| |
| except:
| |
| print 'Failed for', a1[0], a2[0]
| |
| | |
| value = float(u_sum) / twoN
| |
| if not quiet:
| |
| print 'PMF: %.4f (%d contacts, %d residues)' % (value, count, twoN)
| |
| return value
| |
| | |
| try:
| |
| from pymol import cmd
| |
| cmd.extend('aaindex2b', aaindex2b)
| |
| cmd.extend('pmf', pmf)
| |
| def pymol_auto_arg_update():
| |
| aaindexkey_sc = cmd.Shortcut(_aaindex.keys())
| |
| cmd.auto_arg[0].update({
| |
| 'aaindex2b' : [ aaindexkey_sc , 'aaindexkey' , ', ' ],
| |
| 'pmf' : [ aaindexkey_sc , 'aaindexkey' , ', ' ],
| |
| })
| |
| cmd.auto_arg[1].update({
| |
| 'aaindex2b' : [ cmd.selection_sc , 'selection' , '' ],
| |
| })
| |
| cmd.auto_arg[2].update({
| |
| 'pmf' : [ cmd.selection_sc , 'selection' , '' ],
| |
| })
| |
| cmd.auto_arg[3].update({
| |
| 'pmf' : [ cmd.selection_sc , 'selection' , '' ],
| |
| })
| |
| _pymol_auto_arg_update = pymol_auto_arg_update
| |
| except:
| |
| pass
| |
| | |
| # vi: ts=4:sw=4:smarttab:expandtab
| |
| </source>
| |
