#! /usr/local/bin/ost
# -*- coding: utf-8 -*-
"""Translate models for Sergey from PDB + extra data into ModelCIF."""
# EXAMPLES for running:
"""
ost translate2modelcif.py ./raw_data ./raw_data/ptm_plddt.all.txt \
./web_dloads/pivot ./modelcif --prefix=F000347 \
--pdb-web-path=./web_dloads/pdb \
--refseq-path=./web_dloads/consensus_all.fasta
"""
# NOTE: add "--compress" for final runs
from timeit import default_timer as timer
import argparse
import filecmp
import gzip
import os
import pandas as pd
import shutil
import sys
import zipfile
import ihm
import ihm.citations
import modelcif
import modelcif.associated
import modelcif.dumper
import modelcif.model
import modelcif.protocol
import modelcif.reference
from ost import io
def _abort_msg(msg, exit_code=1):
"""Write error message and exit with exit_code."""
print(f"{msg}\nAborting.", file=sys.stderr)
sys.exit(exit_code)
def _warn_msg(msg):
"""Write a warning message to stdout."""
print(f"WARNING: {msg}")
def _check_file(file_path):
"""Make sure a file exists and is actually a file."""
if not os.path.exists(file_path):
_abort_msg(f"File not found: '{file_path}'.")
if not os.path.isfile(file_path):
_abort_msg(f"File path does not point to file: '{file_path}'.")
def _check_folder(dir_path):
"""Make sure a file exists and is actually a file."""
if not os.path.exists(dir_path):
_abort_msg(f"Path not found: '{dir_path}'.")
if not os.path.isdir(dir_path):
_abort_msg(f"Path does not point to a directory: '{dir_path}'.")
def _check_opts_folder(dir_path):
"""Remove trailing '/' (return fixed one) and check if path valid."""
if dir_path.endswith("/"):
dir_path = dir_path[:-1]
_check_folder(dir_path)
return dir_path
def _parse_args():
"""Parse command line arguments."""
parser = argparse.ArgumentParser(
formatter_class=argparse.RawDescriptionHelpFormatter,
description=__doc__,
)
parser.add_argument(
"model_base_dir",
type=str,
metavar="<MODEL BASE DIR>",
help="Directory with pub_data* directories with model PDBs.",
)
parser.add_argument(
"metadata_file",
type=str,
metavar="<METADATA FILE>",
help="Path to table with metadata.",
)
parser.add_argument(
"msa_data_dir",
type=str,
metavar="<MSA DIR>",
help="Directory with F*.fasta files for custom MSAs.",
)
parser.add_argument(
"out_dir",
type=str,
metavar="<OUTPUT DIR>",
help="Path to directory to store results.",
)
parser.add_argument(
"--prefix",
type=str,
metavar="<PREFIX>",
default="",
help="Only process families starting with given prefix. By default " \
+ "all families are processed.",
)
parser.add_argument(
"--compress",
default=False,
action="store_true",
help="Compress ModelCIF file with gzip "
"(note that acc. data is zipped either way).",
)
parser.add_argument(
"--pdb-web-path",
type=str,
metavar="<PDB WEB PATH>",
default=None,
help="Optional path to directory with F*.pdb files as available on " \
+ "NMPFamsDB web site. Used to check if top ranked model the same.",
)
parser.add_argument(
"--refseq-path",
type=str,
metavar="<PDB REFSEQ PATH>",
default=None,
help="Optional path to fasta file with all ref. seq. as available on " \
+ "NMPFamsDB web site. Used to check if it matches with MSA.",
)
opts = parser.parse_args()
# check input
opts.model_base_dir = _check_opts_folder(opts.model_base_dir)
_check_file(opts.metadata_file)
opts.msa_data_dir = _check_opts_folder(opts.msa_data_dir)
# check out_dir
if opts.out_dir.endswith("/"):
opts.out_dir = opts.out_dir[:-1]
if not os.path.exists(opts.out_dir):
os.makedirs(opts.out_dir)
if not os.path.isdir(opts.out_dir):
_abort_msg(f"Path '{opts.out_dir}' does not point to a directory.")
# check optional paths
if opts.pdb_web_path is not None:
opts.pdb_web_path = _check_opts_folder(opts.pdb_web_path)
if opts.refseq_path is not None:
_check_file(opts.refseq_path)
return opts
# pylint: disable=too-few-public-methods
class _GlobalPTM(modelcif.qa_metric.Global, modelcif.qa_metric.PTM):
"""Predicted accuracy according to the TM-score score in [0,1]"""
name = "pTM"
software = None
class _GlobalPLDDT(modelcif.qa_metric.Global, modelcif.qa_metric.PLDDT):
"""Predicted accuracy according to the CA-only lDDT in [0,100]"""
name = "pLDDT"
software = None
class _LocalPLDDT(modelcif.qa_metric.Local, modelcif.qa_metric.PLDDT):
"""Predicted accuracy according to the CA-only lDDT in [0,100]"""
name = "pLDDT"
software = None
class _NmpfamsdbTrgRef(modelcif.reference.TargetReference):
"""PDB as target reference."""
name = "Other"
other_details = "NMPFamsDB"
# pylint: enable=too-few-public-methods
def _get_res_num(r, use_auth=False):
"""Get res. num. from auth. IDs if reading from mmCIF files."""
if use_auth:
return int(r.GetStringProp("pdb_auth_resnum"))
else:
return r.number.num
def _get_ch_name(ch, use_auth=False):
"""Get chain name from auth. IDs if reading from mmCIF files."""
if use_auth:
return ch.GetStringProp("pdb_auth_chain_name")
else:
return ch.name
class _OST2ModelCIF(modelcif.model.AbInitioModel):
"""Map OST entity elements to ihm.model"""
def __init__(self, *args, **kwargs):
"""Initialise a model"""
for i in ["ost_entity", "asym", "scores_json"]:
if i not in kwargs:
raise TypeError(f"Required keyword argument '{i}' not found.")
self.ost_entity = kwargs.pop("ost_entity")
self.asym = kwargs.pop("asym")
self.scores_json = kwargs.pop("scores_json")
# use auth IDs for res. nums and chain names
self.use_auth = False
# fetch plddts per residue
self.plddts = []
for res in self.ost_entity.residues:
b_factors = [a.b_factor for a in res.atoms]
assert len(set(b_factors)) == 1 # must all be equal!
self.plddts.append(b_factors[0])
super().__init__(*args, **kwargs)
def get_atoms(self):
# ToDo [internal]: Take B-factor out since its not a B-factor?
# NOTE: this assumes that _get_res_num maps residue to pos. in seqres
# within asym
for atm in self.ost_entity.atoms:
yield modelcif.model.Atom(
asym_unit=self.asym[_get_ch_name(atm.chain, self.use_auth)],
seq_id=_get_res_num(atm.residue, self.use_auth),
atom_id=atm.name,
type_symbol=atm.element,
x=atm.pos[0],
y=atm.pos[1],
z=atm.pos[2],
het=atm.is_hetatom,
biso=atm.b_factor,
occupancy=atm.occupancy,
)
def add_scores(self):
"""Add QA metrics from AF2 scores."""
# global scores
self.qa_metrics.extend(
(
_GlobalPLDDT(self.scores_json["plddt_global"]),
_GlobalPTM(self.scores_json["ptm"]),
)
)
# local scores
i = 0
for chn_i in self.ost_entity.chains:
ch_name = _get_ch_name(chn_i, self.use_auth)
for res_i in chn_i.residues:
# local pLDDT
res_num = _get_res_num(res_i, self.use_auth)
self.qa_metrics.append(
_LocalPLDDT(
self.asym[ch_name].residue(res_num),
self.plddts[i],
)
)
i += 1
def _get_audit_authors():
"""Return the list of authors that produced this model."""
return (
"Pavlopoulos, Georgios A.",
"Baltoumas, Fotis A.",
"Liu, Sirui",
"Selvitopi, Oguz",
"Camargo, Antonio Pedro",
"Nayfach, Stephen",
"Azad, Ariful",
"Roux, Simon",
"Call, Lee",
"Ivanova, Natalia N.",
"Chen, I-Min",
"Paez-Espino, David",
"Karatzas, Evangelos",
"Novel Metagenome Protein Families Consortium",
"Iliopoulos, Ioannis",
"Konstantinidi, Konstantinos",
"Tiedje, James M.",
"Pett-Ridge, Jennifer",
"Baker, David",
"Visel, Axel",
"Ouzounis, Christos A.",
"Ovchinnikov, Sergey",
"Buluc, Aydin",
"Kyrpides, Nikos C.",
)
def _get_metadata(metadata_file):
"""Read csv file with metedata and prepare for next steps."""
metadata = pd.read_csv(metadata_file, sep=' ',
names=["ID", "mdl", "pTM", "pLDDT"])
return metadata
def _get_pdb_files(model_base_dir):
"""Collect PDB files from pub_data_* folders.
Returns dict with key = family name and value = list of paths to PDB files.
"""
pdb_files_split = dict() # to return
pdb_files_raw = set() # to check for duplicates
pub_paths = [f for f in os.listdir(model_base_dir) \
if f.startswith("pub_data_")]
# NOTE: we sort pub_paths to ensure that pub_data_02 is before _03
for pub_path in sorted(pub_paths):
sub_path = os.path.join(model_base_dir, pub_path)
pdb_files_new = [f for f in os.listdir(sub_path) \
if not f.startswith(".")]
for f in pdb_files_new:
f_path = os.path.join(sub_path, f)
f_name = f.split('_')[0]
if f_name in pdb_files_split:
pdb_files_split[f_name].append(f_path)
else:
pdb_files_split[f_name] = [f_path]
# check global list
pdb_files_new_set = set(pdb_files_new)
new_duplicates = pdb_files_raw.intersection(pdb_files_new_set)
if new_duplicates:
_warn_msg(f"{len(new_duplicates)} duplicated files found in " \
f"{sub_path}.")
pdb_files_raw = pdb_files_raw.union(pdb_files_new_set)
return pdb_files_split
def _get_config():
"""Define AF setup."""
msa_description = "MSA created by calculating the central or \"pivot\" " \
"sequence of each seed MSA, and refining each " \
"alignment using that sequence as the guide."
mdl_description = "Model generated using AlphaFold (v2.0.0 with models " \
"finetuned to return ptm weights) producing 5 models, " \
"without model relaxation, without templates, ranked " \
"by pLDDT, starting from a custom MSA."
af_config = {}
return {
"af_config": af_config,
"af_version": "2.0.0",
"mdl_description": mdl_description,
"msa_description": msa_description,
"use_templates": False,
"use_small_bfd": False,
"use_multimer": False,
}
def _get_protocol_steps_and_software(config_data):
"""Create the list of protocol steps with software and parameters used."""
protocol = []
# MSA step
step = {
"method_type": "coevolution MSA",
"name": None,
"details": config_data["msa_description"],
}
step["input"] = "target_sequences"
step["output"] = "MSA"
step["software"] = []
step["software_parameters"] = {}
protocol.append(step)
# modelling step
step = {
"method_type": "modeling",
"name": None,
"details": config_data["mdl_description"],
}
# get input data
# Must refer to data already in the JSON, so we try keywords
step["input"] = "target_sequences_and_MSA"
# get output data
# Must refer to existing data, so we try keywords
step["output"] = "model"
# get software
if config_data["use_multimer"]:
step["software"] = [{
"name": "AlphaFold-Multimer",
"classification": "model building",
"description": "Structure prediction",
"citation": ihm.Citation(
pmid=None,
title="Protein complex prediction with "
+ "AlphaFold-Multimer.",
journal="bioRxiv",
volume=None,
page_range=None,
year=2021,
authors=[
"Evans, R.",
"O'Neill, M.",
"Pritzel, A.",
"Antropova, N.",
"Senior, A.",
"Green, T.",
"Zidek, A.",
"Bates, R.",
"Blackwell, S.",
"Yim, J.",
"Ronneberger, O.",
"Bodenstein, S.",
"Zielinski, M.",
"Bridgland, A.",
"Potapenko, A.",
"Cowie, A.",
"Tunyasuvunakool, K.",
"Jain, R.",
"Clancy, E.",
"Kohli, P.",
"Jumper, J.",
"Hassabis, D.",
],
doi="10.1101/2021.10.04.463034",
),
"location": "https://github.com/deepmind/alphafold",
"type": "package",
"version": config_data["af_version"],
}]
else:
step["software"] = [{
"name": "AlphaFold",
"classification": "model building",
"description": "Structure prediction",
"citation": ihm.citations.alphafold2,
"location": "https://github.com/deepmind/alphafold",
"type": "package",
"version": config_data["af_version"],
}]
step["software_parameters"] = config_data["af_config"]
protocol.append(step)
return protocol
def _get_title(fam_name):
"""Get a title for this modelling experiment."""
return f"AlphaFold model for NMPFamsDB Family {fam_name}"
def _get_model_details(fam_name):
"""Get the model description."""
db_url = f"https://bib.fleming.gr/NMPFamsDB/family?id={fam_name}"
# TODO: check if ok to use HTML for the URL
db_url = f"<a href=\"{db_url}\" target=\"_blank\">{db_url}</a>"
return f"Model generated using AlphaFold (v2.0.0) for the " \
f"\"Representative Sequence\" of NMPFamsDB Metagenome / " \
f"Metatranscriptome Family {fam_name}.\n\nSee {db_url} for " \
f"additional details."
def _get_model_group_name():
"""Get a name for a model group."""
return None
def _get_sequence(chn, use_auth=False):
"""Get the sequence out of an OST chain incl. '-' for gaps in resnums."""
# initialise (add gaps if first is not at num. 1)
lst_rn = _get_res_num(chn.residues[0], use_auth)
idx = 1
sqe = "-" * (lst_rn - 1) \
+ chn.residues[0].one_letter_code
for res in chn.residues[idx:]:
lst_rn += 1
while lst_rn != _get_res_num(res, use_auth):
sqe += "-"
lst_rn += 1
sqe += res.one_letter_code
return sqe
def _get_entities(pdb_file, ref_seq, fam_name):
"""Gather data for the mmCIF (target) entities."""
_check_file(pdb_file)
ost_ent = io.LoadPDB(pdb_file)
if ost_ent.chain_count != 1:
raise RuntimeError(
f"Unexpected oligomer in {pdb_file}"
)
chn = ost_ent.chains[0]
sqe_gaps = _get_sequence(chn)
# NOTE: can have gaps to accommodate "X" in ref_seq
exp_seq = sqe_gaps.replace('-', 'X')
len_diff = len(ref_seq.string) - len(exp_seq)
if len_diff > 0:
exp_seq += 'X' * len_diff
if exp_seq != ref_seq.string:
raise RuntimeError(f"Sequence in {pdb_file} does not match ref_seq")
# TODO: waiting for input on whether they change handling of "X" in seq
# -> HERE: assuming that "X" were in modelled sequence and PDB can have gaps
cif_ent = {
"seqres": ref_seq.string, # HACK for testing: .replace('X', 'A')
"pdb_sequence": sqe_gaps,
"pdb_chain_id": [_get_ch_name(chn, False)],
"fam_name": fam_name,
"description": f"Representative Sequence of NMPFamsDB Family {fam_name}"
}
return [cif_ent], ost_ent
def _get_modelcif_entities(target_ents, asym_units, system):
"""Create ModelCIF entities and asymmetric units."""
for cif_ent in target_ents:
mdlcif_ent = modelcif.Entity(
# NOTE: sequence here defines residues in model!
cif_ent["seqres"],
description=cif_ent["description"],
source=None,
references=[
_NmpfamsdbTrgRef(
cif_ent["fam_name"],
cif_ent["fam_name"],
align_begin=1,
align_end=len(cif_ent["seqres"])
)
],
)
# NOTE: this assigns (potentially new) alphabetic chain names
for pdb_chain_id in cif_ent["pdb_chain_id"]:
asym_units[pdb_chain_id] = modelcif.AsymUnit(
mdlcif_ent, strand_id=pdb_chain_id,
)
system.target_entities.append(mdlcif_ent)
def _get_assoc_aln_file(fle_path):
"""Generate a modelcif.associated.File object pointing to FASTA formatted
file containing MSA.
"""
cfile = modelcif.associated.File(
fle_path,
details="Custom MSA for modelling",
data=modelcif.data.Data("Custom MSA for modelling")
)
cfile.file_format = "fasta"
cfile.file_content = "multiple sequence alignments"
return cfile
def _get_associated_files(mdl_name, arc_files):
"""Create entry for associated files."""
# package all into zip file
return modelcif.associated.Repository(
"",
[modelcif.associated.ZipFile(f"{mdl_name}.zip", files=arc_files)],
)
# NOTE: by convention MA expects zip file with same name as model-cif
def _assemble_modelcif_software(soft_dict):
"""Create a modelcif.Software instance from dictionary."""
return modelcif.Software(
soft_dict["name"],
soft_dict["classification"],
soft_dict["description"],
soft_dict["location"],
soft_dict["type"],
soft_dict["version"],
citation=soft_dict["citation"],
)
def _get_modelcif_protocol_software(js_step):
"""Assemble software entries for a ModelCIF protocol step."""
if js_step["software"]:
if len(js_step["software"]) == 1:
sftwre = _assemble_modelcif_software(js_step["software"][0])
else:
sftwre = []
for sft in js_step["software"]:
sftwre.append(_assemble_modelcif_software(sft))
sftwre = modelcif.SoftwareGroup(elements=sftwre)
if js_step["software_parameters"]:
params = []
for key, val in js_step["software_parameters"].items():
params.append(modelcif.SoftwareParameter(key, val))
if isinstance(sftwre, modelcif.SoftwareGroup):
sftwre.parameters = params
else:
sftwre = modelcif.SoftwareGroup(
elements=(sftwre,), parameters=params
)
return sftwre
return None
def _get_modelcif_protocol_data(data_label, target_entities, aln_data, model):
"""Assemble data for a ModelCIF protocol step."""
if data_label == "target_sequences":
data = modelcif.data.DataGroup(target_entities)
elif data_label == "MSA":
data = aln_data
elif data_label == "target_sequences_and_MSA":
data = modelcif.data.DataGroup(target_entities)
data.append(aln_data)
elif data_label == "model":
data = model
else:
raise RuntimeError(f"Unknown protocol data: '{data_label}'")
return data
def _get_modelcif_protocol(protocol_steps, target_entities, aln_data, model):
"""Create the protocol for the ModelCIF file."""
protocol = modelcif.protocol.Protocol()
for js_step in protocol_steps:
sftwre = _get_modelcif_protocol_software(js_step)
input_data = _get_modelcif_protocol_data(
js_step["input"], target_entities, aln_data, model
)
output_data = _get_modelcif_protocol_data(
js_step["output"], target_entities, aln_data, model
)
protocol.steps.append(
modelcif.protocol.Step(
input_data=input_data,
output_data=output_data,
name=js_step["name"],
details=js_step["details"],
software=sftwre,
)
)
protocol.steps[-1].method_type = js_step["method_type"]
return protocol
def _compress_cif_file(cif_file):
"""Compress cif file and delete original."""
with open(cif_file, 'rb') as f_in:
with gzip.open(cif_file + '.gz', 'wb') as f_out:
shutil.copyfileobj(f_in, f_out)
os.remove(cif_file)
def _package_associated_files(repo):
"""Compress associated files into single zip file and delete original."""
# zip settings tested for good speed vs compression
for archive in repo.files:
with zipfile.ZipFile(archive.path, "w", zipfile.ZIP_BZIP2) as cif_zip:
for zfile in archive.files:
cif_zip.write(zfile.path, arcname=zfile.path)
os.remove(zfile.path)
def _store_as_modelcif(data_json, ost_ent, out_dir, mdl_name, compress):
"""Mix all the data into a ModelCIF file."""
print(" generating ModelCIF objects...", end="")
pstart = timer()
# create system to gather all the data
system = modelcif.System(
title=data_json["title"],
id=data_json["mdl_id"].upper(),
model_details=data_json["model_details"],
)
# create an asymmetric unit and an entity per target sequence
asym_units = {}
_get_modelcif_entities(
data_json["target_entities"], asym_units, system
)
# audit_authors
system.authors.extend(data_json["audit_authors"])
# set up the model to produce coordinates
model = _OST2ModelCIF(
assembly=modelcif.Assembly(asym_units.values()),
asym=asym_units,
ost_entity=ost_ent,
scores_json=data_json,
name=data_json["mdl_name"],
)
print(f" ({timer()-pstart:.2f}s)")
print(" processing QA scores...", end="", flush=True)
pstart = timer()
model.add_scores()
print(f" ({timer()-pstart:.2f}s)")
model_group = modelcif.model.ModelGroup(
[model], name=data_json["model_group_name"]
)
system.model_groups.append(model_group)
# handle additional files
aln_file = _get_assoc_aln_file(data_json["aln_file_name"])
system.repositories.append(
_get_associated_files(mdl_name, [aln_file])
)
# get data and steps
protocol = _get_modelcif_protocol(
data_json["protocol"], system.target_entities, aln_file.data, model,
)
system.protocols.append(protocol)
# write modelcif System to file (NOTE: no PAE here!)
print(" write to disk...", end="", flush=True)
pstart = timer()
# copy aln file to compress them
shutil.copyfile(
data_json["aln_file_path"],
os.path.join(out_dir, data_json["aln_file_name"])
)
# NOTE: we change path and back while being exception-safe to handle zipfile
oldpwd = os.getcwd()
os.chdir(out_dir)
mdl_fle = f"{mdl_name}.cif"
try:
with open(mdl_fle, "w", encoding="ascii") as mmcif_fh:
modelcif.dumper.write(mmcif_fh, [system])
_package_associated_files(system.repositories[0])
if compress:
_compress_cif_file(mdl_fle)
mdl_fle += ".gz"
finally:
os.chdir(oldpwd)
print(f" ({timer()-pstart:.2f}s)")
def _translate2modelcif(f_name, opts, metadata_fam, pdb_files, ref_seq_check):
"""Convert a model with its accompanying data to ModelCIF."""
# TODO: unclear what to do with such cases; skipped for now
if len(metadata_fam) != 5:
_warn_msg(f"Unexpected number of {len(metadata_fam)} models in " \
f"metadata for family {f_name}. Skipping...")
return
#
mdl_id = f_name
# skip if done already
if opts.compress:
cifext = "cif.gz"
else:
cifext = "cif"
mdl_path = os.path.join(opts.out_dir, f"{mdl_id}.{cifext}")
if os.path.exists(mdl_path):
print(f" {mdl_id} already done...")
return
# go for it...
print(f" translating {mdl_id}...")
pdb_start = timer()
# get aln_data and ref. seq. for this entry
aln_file = f"{f_name}.fasta"
aln_path = os.path.join(opts.msa_data_dir, aln_file)
# TODO: if we need to handle files without ALN, this needs fixing
# -> e.g. 11 extra models in pub_data_* cannot be handled right now
if not os.path.exists(aln_path):
_warn_msg(f"Cannot deal with missing MSA for {f_name} (yet). " \
f"Skipping...")
return
#
aln = io.LoadAlignment(aln_path) # note: this checks that it's an actual MSA
ref_seq = aln.sequences[0]
if ref_seq_check is not None and ref_seq_check.string != ref_seq.string:
raise RuntimeError(f"Sequence mismatch for {f_name}")
# TODO: allow "X" (or whatever can be used to label unknown AA) if needed
if "X" in ref_seq.string:
_warn_msg(f"Cannot deal with 'X' in ref_seq for {f_name} (yet). " \
f"Skipping...")
return
#
# gather data into JSON-like structure
print(" preparing data...", end="")
pstart = timer()
config_data = _get_config()
mdlcf_json = {}
mdlcf_json["audit_authors"] = _get_audit_authors()
mdlcf_json["protocol"] = _get_protocol_steps_and_software(config_data)
mdlcf_json["config_data"] = config_data
mdlcf_json["mdl_id"] = mdl_id
mdlcf_json["aln_file_name"] = aln_file
mdlcf_json["aln_file_path"] = aln_path
# find model to process
# TODO: here just top pLDDT model processed; extend for more if needed...
top_metadata = metadata_fam.loc[metadata_fam.pLDDT.idxmax()]
pdb_list_sel = [f for f in pdb_files if top_metadata.mdl in f]
if len(pdb_list_sel) != 1:
# this should only happen if duplicated file in pub_data_*
# TODO: for now no warning shown and we just pick first hit
# -> first hit is from lowest "pub_data_*"
# -> unclear if we should worry about it...
pass
mdlcf_json["mdl_name"] = f"Top ranked model ({top_metadata.mdl})"
# get scores for this entry
mdlcf_json["plddt_global"] = top_metadata.pLDDT
mdlcf_json["ptm"] = top_metadata.pTM
# process coordinates
pdb_file = pdb_list_sel[0]
target_entities, ost_ent = _get_entities(pdb_file, ref_seq, f_name)
mdlcf_json["target_entities"] = target_entities
# sanity check (only for top ranked model!)
if opts.pdb_web_path is not None:
pdb_file_web = os.path.join(opts.pdb_web_path, f"{f_name}.pdb")
if not filecmp.cmp(pdb_file, pdb_file_web):
raise RuntimeError(f"PDB file mismatch web vs top-ranked for " \
f"{f_name}")
# fill annotations
mdlcf_json["title"] = _get_title(f_name)
mdlcf_json["model_details"] = _get_model_details(f_name)
mdlcf_json["model_group_name"] = _get_model_group_name()
print(f" ({timer()-pstart:.2f}s)")
# save ModelCIF
_store_as_modelcif(mdlcf_json, ost_ent, opts.out_dir, mdl_id, opts.compress)
# check if result can be read and has expected seq.
ent, ss = io.LoadMMCIF(mdl_path, seqres=True)
exp_seqs = [trg_ent["pdb_sequence"] \
for trg_ent in mdlcf_json["target_entities"]]
assert ent.chain_count == len(exp_seqs), f"Bad chain count {mdl_id}"
# here we expect auth = label IDs
ent_seq = "".join([_get_sequence(chn, False) for chn in ent.chains])
ent_seq_a = "".join([_get_sequence(chn, True) for chn in ent.chains])
assert ent_seq == ent_seq_a
assert ent_seq == "".join(exp_seqs), f"Bad seq. {mdl_id}"
ent_seqres = "".join([ss.FindSequence(chn.name).string \
for chn in ent.chains])
exp_seqres = "".join([trg_ent["seqres"] \
for trg_ent in mdlcf_json["target_entities"]])
assert ent_seqres == exp_seqres, f"Bad seqres {mdl_id}"
print(f" ... done with {mdl_id} ({timer()-pdb_start:.2f}s).")
def _main():
"""Run as script."""
opts = _parse_args()
# parse/fetch global data
metadata_full = _get_metadata(opts.metadata_file)
pdb_files_split = _get_pdb_files(opts.model_base_dir)
if opts.refseq_path is not None:
refseqs = io.LoadSequenceList(opts.refseq_path)
else:
refseqs = None
# get on with models
print(f"Working on {opts.prefix}*...")
# iterate over models
for f_name in sorted(pdb_files_split):
if f_name.startswith(opts.prefix):
_translate2modelcif(
f_name, opts,
metadata_full[metadata_full.ID == f_name],
pdb_files_split[f_name],
refseqs.FindSequence(f_name) if refseqs is not None else None,
)
print(f"... done with {opts.prefix}*.")
if __name__ == "__main__":
_main()