Back

Pymol加载对称单元

有时候有一些蛋白是对称单元组成的,如铁蛋白,但是晶体结构pdb坐标只有一部分,另外一部分写在信息里,我们如果想看全部或者需要全部的pdb坐标,可以使用如下方法实现:

Assembly

该功能需要PyMOL 1.8以上 设置对称单元,如下:

1
2
3
4
5
6
7
set assembly, ""
fetch 3bw1, asu, async=0

set assembly, 1
fetch 3bw1, assembly1, async=0

set grid_mode

BiologicalUnit/Quat

该方法可以重构pdb的另外一边

python代码

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
'''
(c) 2010-2011 Thomas Holder, MPI for Developmental Biology

Module for reading REMARK records from PDB files and in particular
generate quaterny structure from REMARK 350.
'''

import sys, os
from pymol import cmd, stored

local_mirror_divided = '/mnt/bio/db/pdb.divided'

def pdbremarks(filename):
    '''
    Read REMARK lines from PDB file. Return dictionary with remarkNum as key
    and list of lines as value.
    '''
    remarks = dict()
    if not isinstance(filename, basestring):
        f = filename
    elif filename[-3:] == '.gz':
        import gzip
        f = gzip.open(filename)
    else:
        f = open(filename)
    for line in f:
        recname = line[0:6]
        if recname == 'REMARK':
            num = int(line[7:10])
            lstring = line[11:]
            remarks.setdefault(num, []).append(lstring)
    return remarks

def quat350(rem350):
    '''
    Get transformation matrices for biomolecule 1 from REMARK 350.
    '''
    biomt = dict()
    chains = tuple()
    seenbiomolecule = False
    for line in rem350:
        if line.startswith('BIOMOLECULE:'):
            if seenbiomolecule:
                break
            seenbiomolecule = True
        elif line.startswith('APPLY THE FOLLOWING TO CHAINS:'):
            chains = tuple(chain.strip() for chain in line[30:].split(','))
        elif line.startswith('                   AND CHAINS:'):
            chains += tuple(chain.strip() for chain in line[30:].split(','))
        elif line.startswith('  BIOMT'):
            row = int(line[7])
            num = int(line[8:12])
            vec = line[12:].split()
            vec = map(float, vec)
            biomt.setdefault(chains, dict()).setdefault(num, []).extend(vec)
    return biomt

def quat(name=None, filename=None, prefix=None, quiet=0):
    '''
DESCRIPTION

    Read REMARK 350 from `filename` and create biological unit
    (quaternary structure)

USAGE

    quat [name [, filename [, prefix]]]

ARGUMENTS

    name = string: name of object and basename of PDB file, if
    filename is not given {default: first loaded object}

    filename = string: file path {default: <name>.pdb}

    prefix = string: prefix for new objects {default: <name>}

EXAMPLE

    fetch 1rmv, type=pdb
    quat 1rmv
    '''
    quiet = int(quiet)
    if name is None:
        name = cmd.get_object_list()[0]
    if prefix is None:
        prefix = name
    if filename is None:
        candidates = [
            '%s.pdb' % (name),
            '%s/%s.pdb' % (cmd.get('fetch_path'), name),
            '%s/%s/pdb%s.ent.gz' % (local_mirror_divided, name[1:3], name),
        ]
        for filename in candidates:
            if os.path.exists(filename):
                break
        else:
            print 'please provide filename'
            return
        if not quiet:
            print 'loading from %s' % (filename)
    remarks = pdbremarks(filename)
    if 350 not in remarks:
        print 'There is no REMARK 350 in', filename
        return
    quat = quat350(remarks[350])
    for chains in quat:
        matrices = quat[chains]
        for num in matrices:
            mat = matrices[num][0:12]
            mat.extend([0,0,0,1])
            copy = '%s_%d' % (prefix, num)
            if not quiet:
                print 'creating %s' % (copy)
            cmd.create(copy, '/%s//%s' % (name, '+'.join(chains)))
            cmd.alter(copy, 'segi="%d"' % (num))
            cmd.transform_object(copy, mat)
    cmd.disable(name)
    cmd.group('%s_quat' % (prefix), '%s_*' % (prefix))

cmd.extend('quat', quat)

# vi:expandtab:smarttab

使用方法非常简单:

1
2
3
fetch 3bw1, type=pdb
quat 3bw1
as cartoon

以下是前后的区别: 前:

后:

参考资料:BiologicalUnit

Licensed under CC BY-NC-SA 4.0