Biopython使用4:PDB结构

Biopython 第十一章:走向3D:PDB模块

PDBParser和MMCIFParser

PDBParser 解析器

使用p=PDBParser(参数)来构建一个PDB解释器. 接收以下参数 (常用QUIET=True ) :

• PERMISSIVE - Evaluated as a Boolean. If false, exceptions in constructing the SMCRA data structure are fatal. If true (DEFAULT),the exceptions are caught, but some residues or atoms will be missing.
• get_header - unused argument kept for historical compatibilty.
• structure_builder - an optional user implemented StructureBuilder class.
• QUIET - Evaluated as a Boolean. If true, warnings issued in constructing the SMCRA data will be suppressed. If false (DEFAULT), they will be shown. These warnings might be indicative of problems in the PDB file!

>>> from Bio.PDB.PDBParser import PDBParser
>>> p = PDBParser()
>>> s = p.get_structure('1d3z','1D3Z.pdb')

解释器对象有以下方法和属性:

• .get_structure(id,file) : 读取结构, 返回Structure对象. 可以用文件名或句柄
• .get_header()/.get_trailer() : 获取PDB文件的文件头和文件尾. 文件头是一个字典. 文件尾是一个列表. 文件头和文件尾要在结构解释后才有信息.
• .PERMISSIVE, .QUIET, .structure_builder 同PDBParser的参数.
• .header与.trailer : 文件头和文件尾信息. 同相应get方法.
• .line_counter: 行数.

读取后获得的Structure对象有多项属性. 其中header属性同PDBParser获得. id属性也是解析时指定.

Bio.PDB.parse_pdb_header(file) 也可以用来读取文件头. 包括head, deposition_date, release_date, structure_method, resolution, structure_reference, journal_reference, author and compound.

读取mmCIF文件

类似PDBParser, 参数支持structure_builder与QUIET. 解释器支持get_structure方法, 有line_counter属性. 要获取类似文件头信息, 要用MMCIF2Dict 方法. 获取的也是一个字典.

>>> from Bio.PDB.MMCIFParser import MMCIFParser
>>> parser = MMCIFParser()
>>> structure = parser.get_structure('1fat', '1fat.cif')
### 要解释mmCIF的额外信息, 需要MMCIF2Dict
>>> from Bio.PDB.MMCIF2Dict import MMCIF2Dict
>>> mmcif_dict = MMCIF2Dict('1fat.cif')

PDBIO类对象保存文件

该类主要用于输出PDB文件. from Bio.PDB import PDBIO来加载后, 直接构建一个io对象. 用set_structure来指定结构后, 再用save方法来保存.

如果想保存部分结构而非全部结构, 可以用Select类来实现. Select支持四种方法来进行选择: accept_model(model), accept_chain(chain), accept_residue(residue), accept_atom(atom). 默认情况下, 每种方法的返回值都为1(即被保存), 可以继承该类自行构建子类来构建自己的选择器, 从而实现强大的功能. 当经常需要选择性保存时, 这种做法比较方便.

>>> from Bio.PDB import PDBIO
### 简单的读写
>>> io = PDBIO()
>>> io.set_structure(s)
>>> io.save('out.pdb')

#### 构建一个选择器来选择GLY. 
>>> class GlySelect(Select):
...     def accept_residue(self, residue):
...         if residue.get_name()=='GLY':
...             return True
...         else:
...             return False
...
>>> io.save('gly_only.pdb', GlySelect())

Structure类对象

Biopython的结构采用SMCRA体系架构（Structure/Model/Chain/Residue/Atom，结构/模型/链/残基/原子）.

一般地，一个实体子类（即原子，残基，链，模型）能通过id作为键来从父类（分别为残基，链，模型，结构）中提取。可以使用len(object)来获取实体的子类的数量.

可以使用child_list = parent_entity.get_list()来获取子对象的列表 (顺序有讲究), 也可以用get_parent()来获取父类. 对于残基和原子, 其id是一个元组, 比较奇怪, 建议获取全列表后再用索引获取.

在SMCRA的所有层次水平，你还可以提取一个完整id 。完整id是包含所有从顶层对象（结构）到当前对象的id的一个元组。一个残基对象的完整id可以这么得到：

>>> full_id = residue.get_full_id()
>>> print full_id
("1abc", 0, "A", ("", 10, "A"))
>>> atom0=residue.get_list()[0]
>>> print atom0.get_full_id()
("1abc", 0, "A", ("", 10, "A"), ('N', ' '))

一些通用的方法:

• get_id() : 可以获取该实例的id. 和entity.id属性一致.
• get_full_id() : 可以获取完整id, 是一个元组. 解析后, 对象full_id属性就会有值.
• get_level() : 可以获取该对象的层级, 分别是SMCRA其中一种.
• get_parent() : 返回上级父类. 对于Structure层面返回空.
• get_list() : 获取子类的列表的copy. 获取的是直接子类. Atom层面没有该方法.
• get_iterator() : 获取子类的生成器. 在循环中比get_list()快.
• get_models/chains/residues/atoms() : 获取某层级所有子类的生成器.
• has_id(id) : 检查是否含有指定id的子类.
• transform(rot,tran) : 对实体原子坐标进行坐标转换.
• xtra : 该对象储存的额外信息. 字典.

以下属性方法基本都有, 但一般不用. 除非很清楚父子类间操作.

• child_dict与child_list : 下级子类储存位置. 用于列出子类. 一个字典, 一个列表.
• add(entity) : 前者添加子类对象到实体, 会自动控制父子间关系.
• copy() : 复制该对象(会去除父类), 会将子类也复制添加进去.
• insert(pos,entity) : 插入子类对象.
• detach_parent() : 移除父类. 父类设置为None.
• detach_child(id) : 通过 id 来移除一个子类.
• set_parent(entity) : 设置父类对象.

## transform(rot, tran) 用法
## rot是右乘旋转矩阵, 3x3 array
## tran是平移向量, 大小3的向量
>>> rotation = rotmat(pi, Vector(1, 0, 0))
>>> translation = array((0, 0, 1), 'f')
>>> entity.transform(rotation, translation)

Structure/Model/Chain

这三层次都比较高. Structure的ID靠指定字符串, Model的ID源自于文件位置, 从0开始的整数. 一般只有一个模型, 有些结构含有多个模型(例如NMR). Chain的ID一般是链的标识符, 单字符(一般是字母), 如'A'. 每个模型的链具有唯一ID. 有多个模型时会出现多个相同ID的链.

• Structure.header : Structure层面储存文件头信息, 字典.
• Model.serial_num : Model层级特有, 标识模型的序号, 从1开始计算. 这个值一般是id+1.
• Chain.get_unpacked_list() : Chain层级有, 返回undisordered残基的列表, 而不包含disordered.

Residue

残基的ID是一个三元组: (het, resseq, icode):

• het 是异质域, 例如'W'是水, 'H_残基名'是非标准氨基酸/核酸. 空值代表标准氨基酸和核酸.
• resseq是残基序列编号. 即在链中的编号. 整形.
• icode是插入码, 当残基含有disordered 信息(即其余同位置的可选残基)时, 会记录这种特点的插入码. 例如'A'.
• 如果第一和第三项为空或空白, 此时可以直接用残基序列编号访问残基, 如chain[60]. 如果含有水, 则不行了(第一项有het). 因为有时水的编号可能和某氨基酸一致…

残基存在disordered状态. 例如存在Thr 80 A, Ser 80 B, Asn 81, 这个80残基可能有突变, 有另外一种可能情况. 这个会在icode中反映出来.

残基类对象有这些特殊方法 :

• is_disordered()或disordered, 如果含有disordered原子, 就会返回True.
• flag_disordered() : 设置disordered flag. 一般不要动.
• get_resname()或resname, 获取残基名, 一般三个大写字母.
• get_segid()或segid : 返回SEGID, 例如"CHN1"
• get_unpacked_list() : 返回undisordered原子的列表.
• sort() : 将原子排序. N, CA,C, O总是最前, 后面根据字母数字顺序.

Bio.PDB.is_aa(residue, standard=False) 方法可以检查氨基酸. residue可以使残基对象, 也可以是3字母字符串. standard设置True, 只检查20个标准氨基酸. 例如FME默认True, 设置后就False.

Atom

原子属性有很多, 除了之前很多SMCRA对象的共有方法外, 还有很多方法获取的是对应的属性.

原子的ID就是他的名称name. 要注意唯一性. 一般就是PDB文件中原子名称去除空格来创建. 而实际在PDB原子名称是可以有空格的.例如'CA '代表的是钙而非Cα (' CA '). 此时为了防止歧义, 会保留空格 (所以要慎防这种奇葩情况的出现, 尽管Ca不常见).

• .element, .mass : 返回元素以及原子量. 这个没有get方法.
• .get_vector() 返回原子坐标的vector.
• .get_coord() / .coord , 返回array([x,y,z],dtype=float32).
• .get_name()/.name, 返回原子名, 一般去掉前后空格, 如'N'.
• .get_fullname()/.fullname, 返回全名(四个字符,带空格). 如' N '.
• .get_serial_number()/.serial_number : 原子序号(从1开始).
• .get_bfactor / .bfactor, 返回bfactor温度因子数据
• .get_altloc/.altloc 返回原子记录的altloc. 是紊乱原子的标识. 一般是' '.
• .get_occupancy/.occupancy : 返回occupancy因子.
• .get_anisou/.anisou_array : 获取anisou数据, 很可能是空的.
• get_sigatm()是原子坐标的标准差, get_siguij()是温度因子的标准差.
• set_XXXX(值), 可以设置一些属性, 如coord, bfactor, anisou, altloc, occupancy, sigatm, siguij, serial_number, parent.
• is_disordered()或disordered, 如果含有disordered原子, 就会返回True.
• flag_disordered() : 设置disordered flag. 一般不要动.

注意: 原子序号项和电荷没有存储!!!

Bio.PDB.Vector 模块及相关运算

• Vector.angle(other) : 返回两个向量的夹角.
• Vector.get_array() : 返回该向量的numpy.array数组.
• Vector.left_multiply(matrix) : 返回矩阵左乘结果, Matrix x Vector.
• Vector.right_multiply(matrix) : 返回矩阵右乘结果, Vector x Matrix.
• Vector.norm() : 返回向量的模(长度).
• Vector.normsq() : 返回向量的模的平方.
• Vector.normalized() : 返回模1的单位向量.
• Vector.normalize() : 将向量变为单位向量(模1), 原位修改, 无返回值.

除了上述方法, Bio.PDB里面也有一些和向量操作相关的方法:

• Bio.PDB.calc_angle(v1,v2,v3) : 计算3个点的夹角, 返回浮点型.
• Bio.PDB.calc_dihedral(v1, v2, v3, v4): 返回4个点的二面角. [-pi,pi]区间.
• Bio.PDB.m2rotaxis(m): 返回旋转矩阵的角度以及轴向量的二元元组. m是rotaxis获得的旋转矩阵.
• Bio.PDB.refmat(p,q) : 镜像p在q的基础上. 返回3x3 array, 左乘后获得p.
• Bio.PDB.rotaxis(theta,vector) : theta是角度,pi计算.旋转p在q的基础上, 返回3x3 array, 左乘后获得p.
• Bio.PDB.rotaxis2m : 同rotaxis.
• Bio.PDB.rotmat(p, q) : 将向量p根据q来旋转矩阵(左乘), 返回3x3 array.
• Bio.PDB.vector_to_axis(line, point) : 点到线的投影的向量. 两个参数都是向量, 后者代表点.
• Bio.PDB.vectors (模块, 内含上述的方法以及Vectors类)

####  获得CB的估算位置
## 获得原子坐标向量
>>> n = residue['N'].get_vector()
>>> c = residue['C'].get_vector()
>>> ca = residue['CA'].get_vector()
## 计算C-CA, N-CA的向量. 
>>> n = n - ca
>>> c = c - ca
## 将N-CA向量绕C-CA轴旋转-120度, 大致可以到CB的位置
# 计算旋转矩阵用于沿CA-C旋转-120度
>>> rot = rotaxis(-pi * 120.0/180.0, c)
# 将旋转矩阵用于N-CA, 从而获得旋转后N-CA向量
>>> cb_at_origin = n.left_multiply(rot)
# 根据旋转后向量, 加上CA的位置, 可以计算出CB的位置.
>>> cb = cb_at_origin+ca

紊乱Disordered残基和原子

紊乱残基和原子分别使用DisorderedResidue和DisorderedAtom来处理. 在列出链内残基时, 相应的会显示紊乱残基, 会选择其中一个残基显示. 但是, 当有A,B多个可选残基时, 他不一定选的是A的残基, 默认会选择最后的一个残基作为缺省!

最好获得链的残基的方法是使用get_unpacked_list()来获取残基, 此时显示的会是A链的残基且是Residue而不是DisorderedResidue.

一个PDB例子是1EN2. 该例子残基10, 14, 16, 80, 81均是紊乱残基, 而1, 90-93是het残基, 94-169是水.

对于紊乱残基, 会有一些特殊的方法:

• is_disordeded() 会返回2, 表明该实体是一个实例的集合.
• disordered_get_list() 返回该紊乱残基可能的残基的列表(元素是Residue)
• disordered_get_id_list() 返回可能的残基名(类型)的列表.
• selected_child 会返回当前紊乱残基选择的残基(默认残基).
• disordered_get(id=None) 不指定时, 返回当前选择的残基. 如果指定了id(残基类型), 就会返回对应的残基.
• disordered_select(id) 可以选择当前紊乱残基的默认残基. 可以使用三字母残基种类来选择, 如disordered_select(‘GLY’)
• disordered_has_id(id) : 检查该残基是否含有某种残基, 使用的id就是残基种类(大写三字母).
• disordered_add(residue) 可以将残基加入到紊乱残基当中, 并且加入相应残基名作为id. 默认会选择该残基作为选择的残基.

from Bio.PDB import PDBParser
p=PDBParser()
dm=p.get_structure('1en2','pdb1en2.ent')
dc=dm[0]['A']
dr=dc.get_list()[9] # 或者 dc[10]
dr.disordered_get_list()
dr.disordered_select('SER')

遍历SMCRA

>>> p = PDBParser()
>>> structure = p.get_structure('X', 'pdb1fat.ent')
>>> for model in structure:
...     for chain in model:
...         for residue in chain:
...             for atom in residue:
...                 print atom
>>> atoms = chain.get_atoms()
>>> for atom in atoms:
...     print atom

另外, Bio.PDB.Selection里的一些方法可以帮助选择出相应的一些实例.

例如:

## 选择一个结构的所有残基, 原子的话相应 'A'
## A=atom, R=residue, C=chain, M=model, S=structure
>>> res_list = Selection.unfold_entities(structure, 'R')
## 这个方法也可以帮助选择一些原子对应的唯一的残基或链等. 
>>> chain_list = Selection.unfold_entities(atom_list, 'C')
## 有专用于选择相应实例的父类实例的
>>> residue_list = Selection.get_unique_parents(atom_list)
## 这个方法用来选择出unique的实例
>>> unique_atom=Selection.uniqueify(atom_list)

多肽对象 Polypeptide

这是一个多肽相关的类. 多肽类实际是列表list继承的类, 储存一系列的残基对象.

该类有两个Builder, 一个是CaPPBuilder, 一个是PPBuilder (也可以在Bio.PDB内获得), 前者根据Cα-Cα距离来建立多肽, 后者通过C-N距离来构建多肽.

# Using C-N
>>> ppb=PPBuilder()
>>> for pp in ppb.build_peptides(structure):
...     print pp.get_sequence()
...
# Using CA-CA
>>> ppb=CaPPBuilder()
>>> for pp in ppb.build_peptides(structure):
...     print pp.get_sequence()

RCSB 蛋白结构数据库

Biopython的PDB模块还能在PDB数据库进行下载. 通过PDBList对象来实现.

PDBList(self, server='ftp://ftp.wwpdb.org', pdb='/Users/hom', obsolete_pdb=None, verbose=True)

• server 参数是PDB服务器.
• pdb 参数指定了下载的位置. 默认是用户目录.

下载PDB文件: retrieve_pdb_file(pdb_code, obsolete=False, pdir=None, file_format=None, overwrite=False)

• pdb_code : 下载的pdb编号, 四个字符.
• file_format 指定格式, 包括mmCif, pdb, xml, mmtf, bundle (适合大结构). 新版本缺省下载PDBx/mmCif格式.
• pdir : 下载目录. 缺省会创建PDB型文件树.

可以在命令行运行python PDBList.py 1fat 来下载. python PDBList.py all /data/pdb 可以下载整个PDB数据库到/data/pdb. 使用-d选项可以将下载所有文件到同一目录而非对应子目录.

>>> from Bio.PDB import PDBList
>>> pdbl = PDBList()
>>> pdbl.retrieve_pdb_file('1FAT', file_format='pdb')

以下是PDBList对象其他一些方法:

• .pdb_server,.local_pdb, .obsolete_pdb 分别是设置的服务器地址, 本地保存地址, 本地保存obsolete地址. .flat_tree是False/True, 是保存文件用PDB树还是一个文件夹储存.
• pdbl.download_pdb_files(pdb_codes, obsolete=False, pdir=None, file_format=None, overwrite=False) 可以下载一系列PDB文件, pdb_codes是PDB号的列表.
• pdbl.update_pdb() 可以更新PDB数据库.
• pdbl.download_entire_pdb(listfile=None, file_format=None) 可以下载整个PDB数据库. listfile 可以指定所有PDB code写入的文件.
• pdbl.download_obsolete_entries参数和download_entire_pdb一样, 只下载obsolete的数据, 支持格式少些.
• pdbl.get_all_entries() : 可以获取所有PDB数据的编号.
• pdbl.get_all_obsolete() : 获取所有过时的pdb编号的列表.
• pdbl.get_seqres_file(savefile='pdb_seqres.txt') : 检索并保存所有PDB数据的序列.
• pdbl.get_recent_changes() : 返回列表包含三个列表,分别是[新加的,修改的,过时的]. 从服务器获取最近变化的列表.
• pdbl.get_status_list(url) : 通过url来获取每周更新状态的pdb列表. 对于查询过往每周的变化时可能很有用.