Biopython---part 1

PDBParser和MMCIFParser

PDBParser 解析器
读取pdb文件，使用p=PDBParser(参数)来构建一个PDB解释器.

from Bio.PDB.PDBParser import PDBParser
p = PDBParser()
structure = p.get_structure(“mmdb_6LU7”,“mmdb_6LU7.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属性也是解析时指定.

MMCIFParser解析器
读取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’)

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

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

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

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

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 : 该对象储存的额外信息. 字典.
以下属性方法基本都有, 但一般不用. 除非很清楚父子类间操作.

• Structure/Model/Chain
  这三层次都比较高. Structure的ID靠指定字符串, Model的ID源自于文件位置, 从0开始的整数. 一般只有一个模型, 有些结构含有多个模型. 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类)

获得原子坐标向量

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

遍历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

获得原子坐标

res_1_CA = structure[0]['A'][56]['CA']
print(res_1_CA.coord)
res_2_CA = structure[0]['A'][327]['CA']
print(res_2_CA.coord)

计算距离

cutoff = 10
binding_residue = []
for res in structure[0].get_residues():
     # skip the LDP residue
    if res == LDP:
        continue
    # skip non-amino acid residues
    elif res.id[0].startswith("H"):
        continue
    else:
        try:
            acarbon = res['CA']
        except:
            continue
        distances = []
        for atom in LDP:
            diff = acarbon.coord - atom.coord
            distances.append(np.sqrt(np.sum(diff*diff)))
        if min(distances) < cutoff:
            binding_residue.append(res)

参考：
https://blog.csdn.net/platinhom/article/details/87964427
https://youtu.be/mL8NPpRxgJA
https://www.youtube.com/watch?v=dxVKG2gNSos