|
|
| (10 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 |
| | + | }} |
| | + | |
| | + | == Introduction == |
| | [[File:AAindexExample.png|200px|thumb|right|Hydrophobicity coloring with KYTJ820101]] | | [[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 14: |
Line 18: |
| | * 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 24: |
Line 28: |
| | 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 == |
| | | | |
| | <source lang="python"> | | <source lang="python"> |
| | import aaindex | | import aaindex |
| − | aaindex.init(path='.')
| |
| − |
| |
| | aaindex2b KYTJ820101 | | aaindex2b KYTJ820101 |
| | spectrum b, yellow_white_blue | | spectrum b, yellow_white_blue |
| Line 43: |
Line 46: |
| | </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'):
| + | * Protscale from [http://www.rubor.de/pymol_extensions_de.html rTools] does a similar job in coloring by amino acid properties |
| − | '''
| |
| − | 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
| |
| − | aaindex2b KYTJ820101
| |
| − | spectrumany b, white yellow forest
| |
| − | show surface
| |
| − | '''
| |
| − | from pymol import cmd, stored
| |
| − | entry = get(key)
| |
| − | median = entry.median()
| |
| − | | |
| − | if int(quiet) != 0:
| |
| − | print entry.desc.strip()
| |
| − | | |
| − | def lookup(resn):
| |
| − | 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
| |
| − | 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>
| |
| | | | |
| − | [[Category:Plugins]]
| |
| | [[Category:Script_Library]] | | [[Category:Script_Library]] |
| | [[Category:Biochemical_Scripts]] | | [[Category:Biochemical_Scripts]] |
| | [[Category:Structural_Biology_Scripts]] | | [[Category:Structural_Biology_Scripts]] |
| | + | [[Category:Pymol-script-repo]] |
