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列表. 对于查询过往每周的变化时可能很有用.
