#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (C) 2014 The ProteinDF development team.
# see also AUTHORS and README if provided.
#
# This file is a part of the ProteinDF software package.
#
# The ProteinDF is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# The ProteinDF is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with ProteinDF. If not, see <http://www.gnu.org/licenses/>.
from __future__ import annotations
from typing import Optional
from .atomgroup import AtomGroup
from .atom import Atom
from .position import Position
from .functions import split_residue_key
import sys
import os
import optparse
import re
import copy
import logging
logger = logging.getLogger(__name__)
[ドキュメント]
class Pdb(object):
""" """
AmberToolsVer = 22
def __init__(self, file_path=None, mode=None):
"""
create empty PDB object
mode: None or amber
"""
self._data = {}
self._ssbonds = []
if file_path:
self.load(file_path)
self._mode = mode
if isinstance(self._mode, str):
self._mode = self._mode.upper()
# modpdb table
self._modpdb_amber_atm_tbl = [
{
"name": "NA",
"symbol": "Na",
"rename": "Na+",
},
{
"name": "CL",
"symbol": "Cl",
"rename": "Cl-",
},
]
self._modpdb_formal_atm_tbl = [
{
"name": "NA",
"symbol": "Na",
"rename": "NA",
},
{
"name": "CL",
"symbol": "Cl",
"rename": "CL",
},
]
self._modpdb_amber_res_tbl = {"NA": "Na+", "CL": "Cl-"}
self._modpdb_formal_res_tbl = {"NA": "NA ", "CL": "CL "}
self._modpdb_amber_resatom_table = {}
self._modpdb_formal_resatom_table = {}
if self.AmberToolsVer < 22:
# leap does not treat NH2 as H, so it ends up adding a duplicate NH2
self._modpdb_amber_resatom_table["NME"] = {
"HN2": "H",
"H1": "HH31",
"H2": "HH32",
"H3": "HH33",
}
# reduce adds NH2 if it is missing
self._modpdb_formal_resatom_table["NME"] = {
"H": "HN2",
"HH31": "H1",
"HH32": "H2",
"HH33": "H3",
}
else:
# leap does not treat NH2 as H, so it ends up adding a duplicate NH2
self._modpdb_amber_resatom_table["NME"] = {
"HN2": "H",
"CH3": "C",
}
# reduce adds NH2 if it is missing
self._modpdb_formal_resatom_table["NME"] = {
"H": "HN2",
"C": "CH3",
}
def __get_debug(self):
if not "_debug" in self.__dict__:
self._debug = False
return self._debug
def __set_debug(self, yn):
self._debug = yn
debug = property(__get_debug, __set_debug)
[ドキュメント]
def renumber(self):
for model_serial, model in self._data.items():
for index in range(len(model)):
model[index]["serial"] = index + 1
[ドキュメント]
def load(self, file_path):
if os.path.isfile(file_path) != True:
logger.critical("file not found: {}".format(file_path))
raise FileNotFoundError("file not found: {}".format(file_path))
model_serial = 1
chain_serial = 0
self._data.setdefault(model_serial, [])
with open(file_path, "r") as fin:
while True:
line = fin.readline()
if len(line) == 0:
break
line = line.rstrip("\n")
# if (len(line) != 80):
# continue
if self.debug == True:
logger.debug(line)
record_name = line[0:6]
if record_name == "SSBOND":
serNum = line[7:10]
chainID1 = str(line[15])
seqNum1 = int(line[17:21])
icode1 = str(line[21])
chainID2 = str(line[29])
seqNum2 = int(line[31:35])
icode2 = str(line[35])
ssbond = (
{"chain_id": chainID1, "seq_num": seqNum1, "i_code": icode1},
{"chain_id": chainID2, "seq_num": seqNum2, "i_code": icode2},
)
self._ssbonds.append(ssbond)
elif (record_name == "ATOM ") or (record_name == "HETATM"):
if len(line) < 80:
line = line + (" " * (80 - len(line)))
serial = int(line[6:11])
name4 = line[12:16]
name = name4.strip()
alt_loc = line[16]
res_name = line[17:20]
chain_id = line[21]
if (chain_id == " ") and (res_name != "WAT"):
chain_id = chr(ord("A") + (chain_serial % 26))
res_seq = line[22:26]
i_code = line[26]
coord_x = line[30:38]
coord_y = line[38:46]
coord_z = line[46:54]
occupancy = line[54:60].strip()
temp_factor = line[60:66].strip()
element = line[76:78].strip()
charge = line[78:80].strip()
item = {}
item["serial"] = serial
item["record_name"] = record_name
item["name"] = name
item["alt_loc"] = alt_loc
# res_name
if res_name in ["HID", "HIE", "HIP"]:
# rename AMBER residue name dialect
res_name = "HIS"
item["res_name"] = res_name
item["chain_id"] = chain_id
item["res_seq"] = int(res_seq)
item["i_code"] = i_code
item["coord"] = [float(coord_x), float(coord_y), float(coord_z)]
if len(occupancy) != 0:
item["occupancy"] = float(occupancy)
else:
item["occupancy"] = 1.0
if len(temp_factor) != 0:
item["temp_factor"] = float(temp_factor)
else:
item["temp_factor"] = 0.0
if len(element) != 0:
# TODO: create an atom conversion table
element = element.strip()
if element == "D":
element = "H"
item["element"] = element
else:
# see https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/tutorials/pdbintro.html
# TODO: change to use a lookup table for the conversion
name4 = name4.upper()
name2 = name4[0:2]
name2s = name2.strip()
if (len(name4.strip()) == 4) and (name2[0] == "H"):
element = "H"
elif name2[0].isnumeric() == True:
element = name2[1]
elif len(name2s) == 2:
element = name2[0] + name2[1].lower()
else:
element = name2[1]
# if name2 == "CL":
# element = "Cl"
# elif name2 == "NA":
# element = "Na"
# elif name2 == "MG":
# element = "Mg"
# elif name2 == "FE":
# element = "Fe"
# else:
# element = name[0]
item["element"] = element
if len(charge) != 0:
charge_last_char = charge[-1]
if (charge_last_char == "+") or (charge_last_char == "-"):
charge = charge_last_char + charge[0:-1]
item["charge"] = charge
else:
item["charge"] = " "
self._data[model_serial].append(item)
continue
elif record_name == "MODEL ":
serial = int(line[10:14])
model_serial = serial
self._data.setdefault(model_serial, [])
chain_serial = 0
continue
elif record_name == "TER ":
if len(line) < 27:
line = line + (" " * (27 - len(line)))
serial = int(line[6:11]) if line[6:11].isdigit() else 0
resname = line[17:20]
chain_id = line[21]
if (chain_id == " ") and (res_name != "WAT"):
chain_id = chr(ord("A") + (chain_serial % 26))
chain_serial += 1
res_seq = line[22:26]
i_code = line[26]
item = {}
item["serial"] = serial
item["record_name"] = record_name
item["res_name"] = res_name
item["chain_id"] = chain_id
item["res_seq"] = res_seq
item["i_code"] = i_code
self._data[model_serial].append(item)
continue
[ドキュメント]
def get_atomgroup(self, select_model: Optional[int] = None, select_altloc: str = "A") -> AtomGroup:
"""
return AtomGroup object
"""
root = AtomGroup()
for model_serial, model_items in self._data.items():
if (select_model == None) or (int(select_model) == int(model_serial)):
model_name = "model_%d" % (model_serial)
model = AtomGroup()
model.name = model_name
for index in range(len(model_items)):
item = model_items[index]
record_name = item["record_name"]
serial = item["serial"]
if (record_name == "ATOM ") or (record_name == "HETATM"):
name = item["name"]
alt_loc = item["alt_loc"]
res_name = item["res_name"]
chain_id = item["chain_id"]
res_seq = int(item["res_seq"])
i_code = item["i_code"]
coord = item["coord"]
occupancy = item.get("occupancy", 1.0)
temp_factor = item.get("temp_factor", 0.0)
element = item.get("element", "X")
charge = item.get("charge", 0.0)
if charge == " ":
charge = 0.0
if chain_id == " ":
chain_id = "_"
if model.has_group(chain_id) == False:
chain = AtomGroup()
chain.name = chain_id
model.set_group(chain_id, chain)
res_key = "%d%s" % (res_seq, i_code.strip())
if model[chain_id].has_group(res_key) == False:
residue = AtomGroup()
residue.name = res_name
model[chain_id].set_group(res_key, residue)
# create atom object -------------------------------
# print(name, res_name, coord)
atom = Atom()
atom.symbol = element
atom.xyz = Position(coord)
atom.name = name
atom.charge = charge
atom_key = "%d_%s" % (serial, name)
# set the atom object ------------------------------
if (alt_loc == " ") or (alt_loc == select_altloc):
model[chain_id][res_key].set_atom(atom_key, atom)
else:
logger.debug(
'skip alt_loc="{alt_loc}" atom: {atom_str}'.format(alt_loc=alt_loc, atom_str=str(atom))
)
for ssbond in self._ssbonds:
chain_id1 = ssbond[0]["chain_id"]
seq_num1 = ssbond[0]["seq_num"]
icode1 = ssbond[0].get("i_code", " ").strip()
res_key1 = "%d%s" % (seq_num1, icode1)
chain_id2 = ssbond[1]["chain_id"]
seq_num2 = ssbond[1]["seq_num"]
icode2 = ssbond[1].get("i_code", " ").strip()
res_key2 = "%d%s" % (seq_num2, icode2)
if (
model.has_group(chain_id1)
and model[chain_id1].has_group(res_key1)
and model[chain_id1][res_key1].has_atomname("SG")
and model.has_group(chain_id2)
and model[chain_id2].has_group(res_key2)
and model[chain_id2][res_key2].has_atomname("SG")
):
SG1 = model[chain_id1][res_key1]["SG"]
SG2 = model[chain_id2][res_key2]["SG"]
model.add_bond(SG1, SG2, 1)
root.set_group(model_name, model)
return root
[ドキュメント]
def set_by_atomgroup(self, atomgroup: AtomGroup, is_charge2tempfactor: bool = False) -> None:
if not isinstance(atomgroup, AtomGroup):
raise TypeError("Expected AtomGroup, got {}".format(type(atomgroup).__name__))
atomgroup = self.get_modpdb_atomgroup(atomgroup)
re_model_serial = re.compile(r"^model_(\d+)")
re_atom_serial = re.compile(r"^(\d+)")
self._data = {}
item = {}
model_serial = 1
for model_key, model in atomgroup.groups():
match_obj = re_model_serial.match(model_key)
if match_obj != None:
model_serial = int(match_obj.group(1))
self._data.setdefault(model_serial, [])
serial = 1
for chain_id, chain in model.groups():
if chain_id != "_":
item["chain_id"] = chain_id
else:
item["chain_id"] = " "
for res_key, residue in chain.groups():
try:
res_seq, i_code_part = split_residue_key(res_key)
i_code = i_code_part if i_code_part != "" else " "
except ValueError:
res_seq = 0
i_code = " "
item["res_seq"] = res_seq
# resname
item["res_name"] = residue.name
has_OXT = False
for key, atom in residue.atoms():
item["record_name"] = "ATOM "
item["serial"] = serial
serial += 1
item["alt_loc"] = " "
item["i_code"] = i_code
# name
name = atom.name
if name.strip() == "OXT":
has_OXT = True
item["name"] = name
item["coord"] = [atom.xyz.x, atom.xyz.y, atom.xyz.z]
item["element"] = atom.symbol
item["charge"] = atom.charge
# assign charge to temperature factor (B-factor)
if is_charge2tempfactor:
item["temp_factor"] = atom.charge
self._data[model_serial].append(copy.copy(item))
if has_OXT == True:
item["record_name"] = "TER "
item["serial"] = serial
serial += 1
self._data[model_serial].append(copy.copy(item))
# TER
if len(self._data[model_serial]) > 0 and (self._data[model_serial][-1]["record_name"] != "TER "):
item["record_name"] = "TER "
item["serial"] = serial
serial += 1
self._data[model_serial].append(copy.copy(item))
self._sort_by_serial()
def _sort_by_serial(self):
for model_serial, model in self._data.items():
model.sort(key=lambda x: int(x["serial"]))
def __str__(self):
occupancy = 1.0
output = ""
for model_serial, model in self._data.items():
output += "MODEL %4d\n" % (model_serial)
for index in range(len(model)):
item = model[index]
record_name = item["record_name"]
serial = int(item["serial"])
if (record_name == "ATOM ") or (record_name == "HETATM"):
name = item["name"]
element = item["element"]
# see. https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/tutorials/pdbintro.html
if len(name) < 4:
if len(element) == 1:
name = f" {name:<3s}"
else:
name = f"{name:<4s}"
alt_loc = item["alt_loc"]
res_name = item["res_name"]
chain_id = item["chain_id"]
res_seq = int(item["res_seq"])
i_code = item["i_code"]
coord = item["coord"]
occupancy = item.setdefault("occupancy", 1.0)
temp_factor = item.setdefault("temp_factor", 1.0)
element = item.setdefault("element", " ")
charge = int(item.setdefault("charge", 0))
if charge == 0:
charge_str = " "
else:
charge_str = f"{charge:+1d}"
element_str = element.upper()
x, y, z = coord[0], coord[1], coord[2]
output += f"ATOM {serial:>5d} {name:>4s}{alt_loc:1s}{res_name:>3s} {chain_id:1s}{res_seq:>4d}{i_code:1s} {x:>8.3f}{y:>8.3f}{z:>8.3f}{occupancy:6.2f}{temp_factor:6.2f} {element_str:>2s}{charge_str:2s}\n"
elif record_name == "TER ":
output += f"TER {serial:>5d} {res_name:3s} {chain_id:1s}{res_seq:>4d}{i_code:1s}\n"
return output
[ドキュメント]
def get_modpdb_atomgroup(self, ag_protein):
"""Change the name of a residue or atom.
Called by 'set_by_atomgroup()'
Args:
ag_protein (AtomGroup): AtomGroup object of protein
Returns:
AtomGroup: renamed object
"""
if not isinstance(ag_protein, AtomGroup):
raise TypeError("Expected AtomGroup, got {}".format(type(ag_protein).__name__))
mode = self._mode
retval = AtomGroup(ag_protein)
for model_key, model in retval.groups():
for chain_key, chain in model.groups():
for res_key, res in chain.groups():
res = self._modpdb_res(res, mode)
for atom_key, atom in res.atoms():
atom = self._modpdb_resatom(res, atom, mode)
atom = self._modpdb_atom(atom, mode)
return retval
def _modpdb_atom(self, atom, mode=None):
new_name = atom.name
atomname = atom.name.strip().upper()
symbol = atom.symbol
if mode == "AMBER":
for item in self._modpdb_amber_atm_tbl:
if (item["name"] == atomname) and (item["symbol"] == symbol):
new_name = item["rename"]
else:
# "FORMAL"
for item in self._modpdb_formal_atm_tbl:
if (item["name"] == atomname) and (item["symbol"] == symbol):
new_name = item["rename"]
atom.name = new_name
return atom
def _modpdb_res(self, res, mode=None):
# rename HIS name in the AMBER mode
if mode == "AMBER":
res = self._rename_to_amber_dialect(res)
# rename res.name by using residue name table
resname = res.name.upper()
resname = resname.strip()
if mode == "AMBER":
if resname in self._modpdb_amber_res_tbl:
res.name = self._modpdb_amber_res_tbl[resname]
else:
if resname in self._modpdb_formal_res_tbl:
res.name = self._modpdb_formal_res_tbl[resname]
return res
def _modpdb_resatom(self, res, atom, mode=None):
resname = res.name.upper()
resname = resname.strip().lstrip()
atom_name = atom.name.strip().lstrip().upper()
if mode == "AMBER":
if resname in self._modpdb_amber_resatom_table:
if atom_name in self._modpdb_amber_resatom_table[resname]:
atom.name = self._modpdb_amber_resatom_table[resname][atom_name]
logger.debug(":{}@{} -> :{}@{}".format(resname, atom_name, resname, atom.name))
else:
if resname in self._modpdb_formal_resatom_table:
if atom_name in self._modpdb_formal_resatom_table[resname]:
atom.name = self._modpdb_formal_resatom_table[resname][atom_name]
return atom
def _rename_to_amber_dialect(self, res):
"""translate HIS to HID, HIE or HIP"""
if not isinstance(res, AtomGroup):
raise TypeError("Expected AtomGroup, got {}".format(type(res).__name__))
if res.name == "HIS":
# check kinds of "HIS"
has_delta_H = False
has_epsilon_H = False
if res.has_atomname("HD1") and res.has_atomname("HD2"):
has_delta_H = True
if res.has_atomname("HE1") and res.has_atomname("HE2"):
has_epsilon_H = True
if has_delta_H and has_epsilon_H:
res.name = "HIP"
elif has_delta_H:
res.name = "HID"
elif has_epsilon_H:
res.name = "HIE"
logger.debug("found HIS: rename HIS to {}".format(res.name))
return res
[ドキュメント]
def main():
# initialize
# parse args
parser = optparse.OptionParser(usage="%prog [options] PDB_FILE", version="%prog 1.0")
parser.add_option("-o", "--output", dest="output_path", help="PDB output file", metavar="FILE")
parser.add_option("-v", "--verbose", dest="verbose", action="store_false", default=False, help="print message")
(opts, args) = parser.parse_args()
if len(args) == 0:
parser.print_help()
sys.exit(1)
# setting
file_path = args[0]
verbose = opts.verbose
#
pdb_obj = Pdb(file_path)
print(pdb_obj)
# end
if __name__ == "__main__":
main()