This is a read-only mirror of pymolwiki.org

Difference between revisions of "AAindex"

From PyMOL Wiki
Jump to navigation Jump to search
(see also)
m (1 revision)
 
(8 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]]
  
Line 20: 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>
  
Line 35: Line 44:
 
spectrum b, yellow_white_blue
 
spectrum b, yellow_white_blue
 
show surface
 
show surface
</source>
 
 
== The Script ==
 
<source lang="python">
 
'''
 
(c) 2010-2011 Thomas Holder, MPI for Developmental Biology
 
 
Python parser for AAindex: Amino Acid Index Database
 
http://www.genome.jp/aaindex/
 
 
PyMOL commands:
 
 
    aaindex2b
 
    pmf
 
'''
 
 
import sys, os
 
 
_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 __getitem__(self, aai):
 
        return self.get(aai)
 
    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 __getitem__(self, aaij):
 
        return self.get(aaij[0], aaij[1])
 
    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[key]
 
 
def _float_or_None(x):
 
    if x == 'NA' or x == '-':
 
        return None
 
    return float(x)
 
 
def init(path=None, index='13'):
 
    '''
 
    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.
 
    '''
 
    index = str(index)
 
    if path is None:
 
        for path in [os.path.split(__file__)[0], '.', cmd.get('fetch_path')]:
 
            if os.path.exists(os.path.join(path, 'aaindex' + index[0])):
 
                break
 
        print >> sys.stderr, 'path =', path
 
    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.
 
    '''
 
    if not os.path.exists(filename):
 
        import urllib
 
        url = 'ftp://ftp.genome.jp/pub/db/community/aaindex/' + os.path.split(filename)[1]
 
        print 'Downloading "%s"' % (url)
 
        filename = urllib.urlretrieve(url, filename)[0]
 
        print 'Saved to "%s"' % (filename)
 
    f = open(filename)
 
 
    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
 
    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
 
 
    # 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])
 
 
    idmap = dict()
 
    cmd.iterate_state(state, '(%s) or (%s)' % (selection1, selection2),
 
            'idmap[model,index] = [(resn,name),(x,y,z)]', space={'idmap': idmap})
 
    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 = idmap[id1]
 
        a2 = 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>
 
</source>
  
Line 462: Line 50:
 
* Protscale from [http://www.rubor.de/pymol_extensions_de.html rTools] does a similar job in coloring by amino acid properties
 
* Protscale from [http://www.rubor.de/pymol_extensions_de.html rTools] does a similar job in coloring by amino acid properties
  
[[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]]

Latest revision as of 03:58, 4 March 2019

Type Python Script
Download aaindex.py
Author(s) Thomas Holder
License BSD
This code has been put under version control in the project Pymol-script-repo

Introduction

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/

This script is a python parser for the AAindex flat files which will be downloaded from ftp://ftp.genome.jp/pub/db/community/aaindex/

The script provides two PyMOL commands (but can also be used without PyMOL).

  • aaindex2b: Loads numerical indices from aaindex1 as b-factors into your structure
  • pmf: Potential of Mean Force (aaindex3)

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.

import aaindex
aaindex.init(path='.')

aaindex.grep('volume')
x = aaindex.get('KRIW790103')
print(x)
print(x.get('A'))

aaindex.init(index='2')
aaindex.grep('blosum')
x = aaindex.get('HENS920102')
print(x.get('A', 'K'))

PyMOL Example

import aaindex
aaindex2b KYTJ820101
spectrum b, yellow_white_blue
show surface

See Also

  • Protscale from rTools does a similar job in coloring by amino acid properties