diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 6476bce78f8fc3708d528125e9090ba470b0e0dc..68e4d2b6354787cb88fa10453e4ea535e004d3b0 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -47,7 +47,7 @@ test:ubuntu22.04:
- make -j 2
- echo " ... done building OST."
- echo " Downloading chemical compounds..."
- - wget ftp://ftp.wwpdb.org/pub/pdb/data/monomers/components.cif.gz
+ - wget https://files.wwpdb.org/pub/pdb/data/monomers/components.cif.gz
- echo " ... done downloading compounds"
- echo " Building a compound library..."
- stage/bin/chemdict_tool create components.cif.gz compounds.chemlib pdb -i
diff --git a/CHANGELOG.txt b/CHANGELOG.txt
index ef6d7297c4ad593b5988e81ecbd21a7e2c58ca3e..f5ed1821edbc59b0a5222f77bbee25603b04d223 100644
--- a/CHANGELOG.txt
+++ b/CHANGELOG.txt
@@ -1,3 +1,39 @@
+Changes in Release 2.9.0
+--------------------------------------------------------------------------------
+
+ * lDDT-PLI now penalizes added model contacts by default.
+ * Updated unassigned reasons (model and target) to better reflect changes
+ to ligand scoring made in OST 2.8.0.
+ * Added CSV output (--output-format csv) and by model (rather than target)
+ ligand output (--by-model-ligand-output) to compare-ligand-structures
+ action.
+ * Improved logging and output of compare-ligand-structures action.
+ * Residue SetChemType() is now exposed in Python.
+ * Allow reading of BIRD compounds in PRDCC format in the compound library.
+ Compounds in PRD formats cannot be read and are rejected cleanly.
+ * The chemdict_tool executable now exits with an error status if no compound
+ were imported.
+ * The SDF writer populates the program name and time line (2) as per SDF
+ specification.
+ * The mmCIF reader now refuses to read in files with multiple data blocks
+ (except in fault tolerant mode), and warns about mmCIF files containing
+ more than one model (atom_site.pdbx_PDB_model_num).
+ * GDT produces slightly higher scores. It's an optimization problem in the
+ end to find the highest possible fraction of CA atoms that superpose within
+ a specified threshold. OpenStructure GDTTS scores are now on average 0.21
+ points lower (range [0, 100]) than LGA when computed on all CASP15 TS models.
+ * Enable DockQ on protein-nucleotide and nucleotide-nucleotide interfaces.
+ * Enable per-atom lDDT scores.
+ * Optionally compute ICS/IPS scores on trimmed model. Reason for that is
+ scoring with incomplete target structures. The full ICS/IPS scores may
+ consider contacts between residues that are not resolved as false positives,
+ even though there is no experimental evidence. The approach with trimmed
+ models removes any model residue that cannot be mapped to the target
+ structure before scoring.
+ * The SDF reader is now aware of IOProfiles and can read files with invalid
+ bond types (order) from RDKit in fault tolerant mode.
+ * Several bug fixes and improvements.
+
Changes in Release 2.8.0
--------------------------------------------------------------------------------
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1ff94f97da8a9eced92f710db9139ef91150ea2e..188ade39134d121af61b5cc7c57b5390a2d6d2a2 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -14,7 +14,7 @@ cmake_policy(SET CMP0060 NEW)
project(OpenStructure CXX C)
set (CMAKE_EXPORT_COMPILE_COMMANDS 1)
set (OST_VERSION_MAJOR 2)
-set (OST_VERSION_MINOR 8)
+set (OST_VERSION_MINOR 9)
set (OST_VERSION_PATCH 0)
set (OST_VERSION_STRING ${OST_VERSION_MAJOR}.${OST_VERSION_MINOR}.${OST_VERSION_PATCH} )
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${PROJECT_SOURCE_DIR}/cmake_support)
diff --git a/actions/ost-compare-ligand-structures b/actions/ost-compare-ligand-structures
index b991037828f9b39a2e15924cd6cbefc2b5a669ca..3c58acd5951f393ec3142219407045d99b9bba4c 100644
--- a/actions/ost-compare-ligand-structures
+++ b/actions/ost-compare-ligand-structures
@@ -32,8 +32,10 @@ definition, and have properly named residues and atoms, in order for
ligand connectivity to be loaded correctly. Ligands loaded from SDF files
are exempt from this restriction, meaning any arbitrary ligand can be assessed.
-Output is written in JSON format (default: out.json). In case of no additional
-options, this is a dictionary with three keys:
+Output can be written in two format: JSON (default) or CSV, controlled by the
+--output-format/-of argument.
+
+Without additional options, the JSON ouput is a dictionary with four keys:
* "model_ligands": A list of ligands in the model. If ligands were provided
explicitly with --model-ligands, elements of the list will be the paths to
@@ -46,6 +48,7 @@ options, this is a dictionary with three keys:
* "status": SUCCESS if everything ran through. In case of failure, the only
content of the JSON output will be \"status\" set to FAILURE and an
additional key: "traceback".
+ * "ost_version": The OpenStructure version used for computation.
Each score is opt-in and the respective results are available in three keys:
@@ -70,9 +73,44 @@ This is a list of data items for each pair of model/reference ligands. The data
items follow the same structure as in "assigned_scores". If no score for a
specific pair of ligands could be computed, "score" and "coverage" are set to
null and a key "reason" is added giving an educated guess why this happened.
+
+CSV output is a table of comma-separated values, with one line for each
+reference ligand (or one model ligand if the --by-model-ligand-output flag was
+set).
+
+The following column is always available:
+
+ * reference_ligand/model_ligand: If reference ligands were provided explicitly
+ with --reference-ligands, elements of the list will be the paths to the
+ ligand SDF file(s). Otherwise, they will be the chain name, residue number
+ and insertion code of the ligand, separated by a dot. If the
+ --by-model-ligand-output flag was set, this will be model ligand instead,
+ following the same rules.
+
+If lDDT-PLI was enabled with --lddt-pli, the following columns are added:
+
+ * "lddt_pli", "lddt_pli_coverage" and "lddt_pli_(model|reference)_ligand"
+ are the lDDT-PLI score result, the corresponding coverage and assigned model
+ ligand (or reference ligand if the --by-model-ligand-output flag was set)
+ if an assignment was found, respectively, empty otherwise.
+ * "lddt_pli_unassigned" is empty if an assignment was found, otherwise it
+ lists the short reason this reference ligand was unassigned.
+
+If BiSyRMSD was enabled with --rmsd, the following columns are added:
+
+ * "rmsd", "rmsd_coverage". "lddt_lp" "bb_rmsd" and
+ "rmsd_(model|reference)_ligand" are the BiSyRMSD, the corresponding
+ coverage, lDDT-LP, backbone RMSD and assigned model ligand (or reference
+ ligand if the --by-model-ligand-output flag was set) if an assignment
+ was found, respectively, empty otherwise.
+ * "rmsd_unassigned" is empty if an assignment was found, otherwise it
+ lists the short reason this reference ligand was unassigned.
+
"""
import argparse
+import csv
+from io import StringIO
import json
import os
import sys
@@ -128,9 +166,9 @@ def _ParseArgs():
"--out",
"--output",
dest="output",
- default="out.json",
- help=("Output file name. The output will be saved as a JSON file. "
- "default: out.json"))
+ default=None,
+ help=("Output file name. "
+ "Default depends on format: out.json or out.csv"))
parser.add_argument(
"-mf",
@@ -154,6 +192,43 @@ def _ParseArgs():
help=("Format of reference file. pdb reads pdb but also pdb.gz, same "
"applies to cif/mmcif. Inferred from filepath if not given."))
+ parser.add_argument(
+ "-of",
+ "--out-format",
+ "--output-format",
+ dest="output_format",
+ choices=["json", "csv"],
+ default="json",
+ help=("Output format, JSON or CSV, in lowercase. "
+ "default: json"))
+
+ parser.add_argument(
+ "-csvm",
+ "--by-model-ligand",
+ "--by-model-ligand-output",
+ dest="output_by_model_ligand",
+ default=False,
+ action="store_true",
+ help=("For CSV output, this flag changes the output so that each line "
+ "reports one model ligand, instead of a reference ligand. "
+ "Has no effect with JSON output."))
+
+ parser.add_argument(
+ "--csv-extra-header",
+ dest="csv_extra_header",
+ default=None,
+ type=str,
+ help=("Extra header prefix for CSV output. This allows adding "
+ "additional annotations (such as target ID, group, etc) to the "
+ "output"))
+
+ parser.add_argument(
+ "--csv-extra-data",
+ dest="csv_extra_data",
+ default=None,
+ type=str,
+ help=("Additional data (columns) for CSV output."))
+
parser.add_argument(
"-mb",
"--model-biounit",
@@ -208,6 +283,7 @@ def _ParseArgs():
"--coverage-delta",
dest="coverage_delta",
default=0.2,
+ type=float,
help=("Coverage delta for partial ligand assignment."))
parser.add_argument(
@@ -215,8 +291,8 @@ def _ParseArgs():
'--verbosity',
dest="verbosity",
type=int,
- default=3,
- help="Set verbosity level. Defaults to 3 (INFO).")
+ default=2,
+ help="Set verbosity level. Defaults to 2 (Script).")
parser.add_argument(
"--full-results",
@@ -239,14 +315,22 @@ def _ParseArgs():
"--lddt-pli-radius",
dest="lddt_pli_radius",
default=6.0,
+ type=float,
help=("lDDT inclusion radius for lDDT-PLI."))
parser.add_argument(
- "--lddt-pli-amc",
- dest="lddt_pli_amc",
- default=False,
+ "--lddt-pli-add-mdl-contacts",
+ dest="lddt_pli_add_mdl_contacts",
+ default=True,
action="store_true",
- help=("Add model contacts (amc) when computing lDDT-PLI."))
+ help=("Add model contacts when computing lDDT-PLI."))
+
+ parser.add_argument(
+ "--no-lddt-pli-add-mdl-contacts",
+ dest="lddt_pli_add_mdl_contacts",
+ default=True,
+ action="store_false",
+ help=("DO NOT add model contacts when computing lDDT-PLI."))
# arguments relevant for rmsd
@@ -261,6 +345,7 @@ def _ParseArgs():
"--radius",
dest="radius",
default=4.0,
+ type=float,
help=("Inclusion radius to extract reference binding site that is used "
"for RMSD computation. Any residue with atoms within this "
"distance of the ligand will be included in the binding site."))
@@ -269,6 +354,7 @@ def _ParseArgs():
"--lddt-lp-radius",
dest="lddt_lp_radius",
default=15.0,
+ type=float,
help=("lDDT inclusion radius for lDDT-LP."))
parser.add_argument(
@@ -280,7 +366,20 @@ def _ParseArgs():
help=("Enumerate all potential binding sites in the model when "
"searching rigid superposition for RMSD computation"))
- return parser.parse_args()
+ parser.add_argument(
+ "-ms",
+ "--max--symmetries",
+ dest="max_symmetries",
+ default=1e4,
+ type=int,
+ help=("If more than that many isomorphisms exist for a target-ligand "
+ "pair, it will be ignored and reported as unassigned."))
+
+ args = parser.parse_args()
+ if args.output is None:
+ args.output = "out.%s" % args.output_format
+
+ return args
def _CheckCompoundLib():
@@ -320,8 +419,6 @@ def _LoadStructure(structure_path, sformat, fault_tolerant, bu_id):
The returned structure has structure_path attached as structure name
"""
- # increase loglevel, as we would pollute the info log with weird stuff
- ost.PushVerbosityLevel(ost.LogLevel.Error)
# Load the structure
if sformat == "mmcif":
if bu_id is not None:
@@ -346,7 +443,6 @@ def _LoadStructure(structure_path, sformat, fault_tolerant, bu_id):
raise Exception(f"No residues found in file: {structure_path}")
# restore old loglevel and return
- ost.PopVerbosityLevel()
entity.SetName(structure_path)
return entity
@@ -366,7 +462,11 @@ def _LoadLigand(file):
"""
Load a single ligand from file names. Return an entity.
"""
- return ost.io.LoadEntity(file, format="sdf")
+ ent = ost.io.LoadEntity(file, format="sdf")
+ ed = ent.EditXCS()
+ ed.RenameChain(ent.chains[0], file)
+ ed.UpdateICS()
+ return ent
def _CleanupStructure(entity):
@@ -425,7 +525,8 @@ def _SetupLDDTPLIScorer(model, model_ligands, reference, reference_ligands, args
substructure_match = args.substructure_match,
coverage_delta = args.coverage_delta,
lddt_pli_radius = args.lddt_pli_radius,
- add_mdl_contacts = args.lddt_pli_amc)
+ add_mdl_contacts = args.lddt_pli_add_mdl_contacts,
+ max_symmetries = args.max_symmetries)
def _SetupSCRMSDScorer(model, model_ligands, reference, reference_ligands, args):
return ligand_scoring_scrmsd.SCRMSDScorer(model, reference,
@@ -436,7 +537,9 @@ def _SetupSCRMSDScorer(model, model_ligands, reference, reference_ligands, args)
substructure_match = args.substructure_match,
coverage_delta = args.coverage_delta,
bs_radius = args.radius,
- lddt_lp_radius = args.lddt_lp_radius)
+ lddt_lp_radius = args.lddt_lp_radius,
+ full_bs_search = args.full_bs_search,
+ max_symmetries = args.max_symmetries)
def _Process(model, model_ligands, reference, reference_ligands, args):
@@ -527,6 +630,7 @@ def _Process(model, model_ligands, reference, reference_ligands, args):
##################
if args.lddt_pli:
+ LogScript("Computing lDDT-PLI scores")
out["lddt_pli"] = dict()
out["lddt_pli"]["assigned_scores"] = list()
for lig_pair in lddtpli_scorer.assignment:
@@ -563,8 +667,8 @@ def _Process(model, model_ligands, reference, reference_ligands, args):
for ref_lig_idx in range(shape[0]):
for mdl_lig_idx in range(shape[1]):
state = int(lddtpli_scorer.state_matrix[(ref_lig_idx, mdl_lig_idx)])
- target_key = out["reference_ligands"][lig_pair[0]]
- model_key = out["model_ligands"][lig_pair[1]]
+ target_key = out["reference_ligands"][ref_lig_idx]
+ model_key = out["model_ligands"][mdl_lig_idx]
if state == 0:
score = float(lddtpli_scorer.score_matrix[(ref_lig_idx, mdl_lig_idx)])
coverage = float(lddtpli_scorer.coverage_matrix[(ref_lig_idx, mdl_lig_idx)])
@@ -589,6 +693,7 @@ def _Process(model, model_ligands, reference, reference_ligands, args):
if args.rmsd:
+ LogScript("Computing BiSyRMSD scores")
out["rmsd"] = dict()
out["rmsd"]["assigned_scores"] = list()
for lig_pair in scrmsd_scorer.assignment:
@@ -611,7 +716,8 @@ def _Process(model, model_ligands, reference, reference_ligands, args):
aux_data["bs_ref_res_mapped"]],
"bs_mdl_res_mapped": [_QualifiedResidueNotation(r) for r in
aux_data["bs_mdl_res_mapped"]],
- "inconsistent_residues": [_QualifiedResidueNotation(r) for r in
+ "inconsistent_residues": [_QualifiedResidueNotation(r[0]) + \
+ "-" +_QualifiedResidueNotation(r[1]) for r in
aux_data["inconsistent_residues"]],
"transform": [transform_data[i:i + 4]
for i in range(0, len(transform_data), 4)]})
@@ -634,8 +740,8 @@ def _Process(model, model_ligands, reference, reference_ligands, args):
for ref_lig_idx in range(shape[0]):
for mdl_lig_idx in range(shape[1]):
state = int(scrmsd_scorer.state_matrix[(ref_lig_idx, mdl_lig_idx)])
- target_key = out["reference_ligands"][lig_pair[0]]
- model_key = out["model_ligands"][lig_pair[1]]
+ target_key = out["reference_ligands"][ref_lig_idx]
+ model_key = out["model_ligands"][mdl_lig_idx]
if state == 0:
score = float(scrmsd_scorer.score_matrix[(ref_lig_idx, mdl_lig_idx)])
coverage = float(scrmsd_scorer.coverage_matrix[(ref_lig_idx, mdl_lig_idx)])
@@ -654,7 +760,8 @@ def _Process(model, model_ligands, reference, reference_ligands, args):
aux_data["bs_ref_res_mapped"]],
"bs_mdl_res_mapped": [_QualifiedResidueNotation(r) for r in
aux_data["bs_mdl_res_mapped"]],
- "inconsistent_residues": [_QualifiedResidueNotation(r) for r in
+ "inconsistent_residues": [_QualifiedResidueNotation(r[0]) + \
+ "-" +_QualifiedResidueNotation(r[1]) for r in
aux_data["inconsistent_residues"]],
"transform": [transform_data[i:i + 4]
for i in range(0, len(transform_data), 4)]})
@@ -669,22 +776,97 @@ def _Process(model, model_ligands, reference, reference_ligands, args):
return out
+def _WriteCSV(out, args):
+ csv_dict = {}
+
+ if args.output_by_model_ligand:
+ ligand_by = "model_ligand"
+ ligand_other = "reference_ligand"
+ else:
+ ligand_by = "reference_ligand"
+ ligand_other = "model_ligand"
+
+ # Always fill-in basic reference ligand info
+ fieldnames = [ligand_by]
+ for ligand in out["%ss" % ligand_by]:
+ csv_dict[ligand] = {
+ ligand_by: ligand,
+ }
+
+ if args.lddt_pli:
+ fieldnames.extend(["lddt_pli", "lddt_pli_coverage",
+ "lddt_pli_%s" % ligand_other, "lddt_pli_unassigned"])
+ for score in out["lddt_pli"]["assigned_scores"]:
+ csv_dict[score[ligand_by]].update({
+ ligand_by: score[ligand_by],
+ "lddt_pli": score["score"],
+ "lddt_pli_coverage": score["coverage"],
+ "lddt_pli_%s" % ligand_other: score[ligand_other],
+ })
+ for ligand, reason in out["lddt_pli"][
+ "%s_unassigned_reason" % ligand_by].items():
+ csv_dict[ligand].update({
+ ligand_by: ligand,
+ "lddt_pli_unassigned": reason[0],
+ })
+
+ if args.rmsd:
+ fieldnames.extend(["rmsd", "lddt_lp", "bb_rmsd", "rmsd_coverage",
+ "rmsd_%s" % ligand_other, "rmsd_unassigned"])
+ for score in out["rmsd"]["assigned_scores"]:
+ csv_dict[score[ligand_by]].update({
+ ligand_by: score[ligand_by],
+ "rmsd": score["score"],
+ "lddt_lp": score["lddt_lp"],
+ "bb_rmsd": score["bb_rmsd"],
+ "rmsd_coverage": score["coverage"],
+ "rmsd_%s" % ligand_other: score[ligand_other],
+ })
+ for ligand, reason in out["rmsd"][
+ "%s_unassigned_reason" % ligand_by].items():
+ csv_dict[ligand].update({
+ ligand_by: ligand,
+ "rmsd_unassigned": reason[0],
+ })
+
+ if args.csv_extra_header or args.csv_extra_data:
+
+ extra_csv = StringIO(
+ args.csv_extra_header + os.linesep + args.csv_extra_data)
+ reader = csv.DictReader(extra_csv)
+ extra_data = next(iter(reader))
+ if None in extra_data:
+ raise ValueError("Not enough columns in --csv-extra-header")
+ fieldnames = reader.fieldnames + fieldnames
+ for ligand, row in csv_dict.items():
+ row.update(extra_data)
+
+ with open(args.output, 'w', newline='') as csvfile:
+ writer = csv.DictWriter(csvfile, fieldnames=fieldnames)
+ writer.writeheader()
+ for row in csv_dict.values():
+ writer.writerow(row)
def _Main():
args = _ParseArgs()
ost.PushVerbosityLevel(args.verbosity)
if args.verbosity < 4:
- sys.tracebacklimit = 100
+ # Hide tracebacks by default
+ # Run script with -v 4 (Verbose) or higher to display them
+ sys.tracebacklimit = 0
_CheckCompoundLib()
try:
# Load structures
+ LogScript("Loading data")
+ LogInfo("Loading reference structure")
reference_format = _GetStructureFormat(args.reference,
sformat=args.reference_format)
reference = _LoadStructure(args.reference,
sformat=reference_format,
bu_id=args.reference_biounit,
fault_tolerant = args.fault_tolerant)
+ LogInfo("Loading model structure")
model_format = _GetStructureFormat(args.model,
sformat=args.model_format)
model = _LoadStructure(args.model,
@@ -693,13 +875,19 @@ def _Main():
fault_tolerant = args.fault_tolerant)
# Load ligands
- model_ligands = _LoadLigands(args.model_ligands)
+ LogInfo("Loading reference ligands")
reference_ligands = _LoadLigands(args.reference_ligands)
+ LogInfo("Loading model ligands")
+ model_ligands = _LoadLigands(args.model_ligands)
# Cleanup
+ LogVerbose("Cleaning up reference structure")
cleaned_reference = _CleanupStructure(reference)
+ LogVerbose("Cleaning up model structure")
cleaned_model = _CleanupStructure(model)
+ LogVerbose("Cleaning up reference ligands")
cleaned_reference_ligands = _CleanupLigands(reference_ligands)
+ LogVerbose("Cleaning up model ligands")
cleaned_model_ligands = _CleanupLigands(model_ligands)
# Validate
@@ -712,18 +900,27 @@ def _Main():
cleaned_reference, cleaned_reference_ligands,
args)
+ out["ost_version"] = ost.__version__
out["status"] = "SUCCESS"
- with open(args.output, 'w') as fh:
- json.dump(out, fh, indent=4, sort_keys=False)
+ if args.output_format == "json":
+ with open(args.output, 'w') as fh:
+ json.dump(out, fh, indent=4, sort_keys=False)
+ else:
+ _WriteCSV(out, args)
+ LogScript("Saved results in %s" % args.output)
except Exception as exc:
- out = dict()
- out["status"] = "FAILURE"
- out["traceback"] = traceback.format_exc(limit=1000)
- etype, evalue, tb = sys.exc_info()
- out["exception"] = " ".join(traceback.format_exception_only(etype, evalue))
- with open(args.output, 'w') as fh:
- json.dump(out, fh, indent=4, sort_keys=False)
+ if args.output_format == "json":
+ out = dict()
+ out["status"] = "FAILURE"
+ out["traceback"] = traceback.format_exc(limit=1000)
+ etype, evalue, tb = sys.exc_info()
+ out["exception"] = " ".join(traceback.format_exception_only(etype, evalue))
+ with open(args.output, 'w') as fh:
+ json.dump(out, fh, indent=4, sort_keys=False)
+ LogWarning("Error information saved in %s" % args.output)
+ else:
+ LogScript("Error encountered, no output saved")
raise
diff --git a/actions/ost-compare-structures b/actions/ost-compare-structures
index 6b0269a4f1b66278811c4625d2ebe88e5a36e6a8..067d154df5ff73a1f9b592f6b0527c64d9508293 100644
--- a/actions/ost-compare-structures
+++ b/actions/ost-compare-structures
@@ -60,7 +60,9 @@ results:
* "min_pep_length"
* "min_nuc_length"
* "lddt_add_mdl_contacts"
+ * "lddt_inclusion_radius"
* "dockq_capri_peptide"
+ * "ost_version"
The pairwise sequence alignments are computed with Needleman-Wunsch using
BLOSUM62 (NUC44 for nucleotides). Many benchmarking scenarios preprocess the
@@ -246,6 +248,25 @@ def _ParseArgs():
"counterparts. Atoms specified in there follow the following "
"format: <chain_name>.<resnum>.<resnum_inscode>.<atom_name>"))
+ parser.add_argument(
+ "--aa-local-lddt",
+ dest="aa_local_lddt",
+ default=False,
+ action="store_true",
+ help=("Compute per-atom lDDT scores with default parameterization "
+ "and store as key \"aa_local_lddt\". Score for each atom is "
+ "accessible by key "
+ "<chain_name>.<resnum>.<resnum_inscode>.<aname>. "
+ "Alpha carbon from residue with number 42 in chain X can be "
+ "extracted with: data[\"aa_local_lddt\"][\"X.42..CA\"]. "
+ "If there is a residue insertion code, lets say A, the atom key "
+ "becomes \"X.42.A.CA\". "
+ "Stereochemical irregularities affecting lDDT are reported as "
+ "keys \"model_clashes\", \"model_bad_bonds\", "
+ "\"model_bad_angles\" and the respective reference "
+ "counterparts. Atoms specified in there follow the following "
+ "format: <chain_name>.<resnum>.<resnum_inscode>.<atom_name>"))
+
parser.add_argument(
"--bb-lddt",
dest="bb_lddt",
@@ -344,7 +365,8 @@ def _ParseArgs():
help=("Compute DockQ scores and its components. Relevant interfaces "
"with at least one contact (any atom within 5A) of the reference "
"structure are available as key \"dockq_reference_interfaces\". "
- "Only interfaces between peptide chains are considered here! "
+ "Protein-protein, protein-nucleotide and nucleotide-nucleotide "
+ "interfaces are considered. "
"Key \"dockq_interfaces\" is a subset of "
"\"dockq_reference_interfaces\" that contains interfaces that "
"can be mapped to the model. They are stored as lists in format "
@@ -382,7 +404,8 @@ def _ParseArgs():
"8A in combination with only considering CB atoms for "
"protein-peptide interactions. "
"Note that the resulting DockQ is not evaluated for these "
- "slightly updated fnat and irmsd (lrmsd stays the same)."
+ "slightly updated fnat and irmsd (lrmsd stays the same). "
+ "Raises an error if reference contains nucleotide chains. "
"This flag has no influence on patch_dockq scores."))
parser.add_argument(
@@ -411,6 +434,32 @@ def _ParseArgs():
"available as keys \"per_interface_ics_precision\", "
"\"per_interface_ics_recall\" and \"per_interface_ics\"."))
+ parser.add_argument(
+ "--ics-trimmed",
+ dest="ics_trimmed",
+ default=False,
+ action="store_true",
+ help=("Computes interface contact similarity (ICS) related scores but "
+ "on a trimmed model. That means that a mapping between model and "
+ "reference is performed and all model residues without reference "
+ "counterpart are removed. As a consequence, model contacts for "
+ "which we have no experimental evidence do not affect the score. "
+ "The effect of these added model contacts without mapping to "
+ "target would be decreased precision and thus lower ics. Recall is "
+ "not affected. Enabling this flag adds the following keys: "
+ "\"ics_trimmed\", \"ics_precision_trimmed\", "
+ "\"ics_recall_trimmed\", \"model_contacts_trimmed\". "
+ "The reference contacts and reference interfaces are the same "
+ "as for ics and available as keys: \"reference_contacts\", "
+ "\"contact_reference_interfaces\". "
+ "All these measures are also available on a per-interface basis "
+ "for each interface in the reference structure that are defined "
+ "as chain pairs with at least one contact (available as key "
+ " \"contact_reference_interfaces\"). The respective metrics are "
+ "available as keys \"per_interface_ics_precision_trimmed\", "
+ "\"per_interface_ics_recall_trimmed\" and "
+ "\"per_interface_ics_trimmed\"."))
+
parser.add_argument(
"--ips",
dest="ips",
@@ -438,6 +487,12 @@ def _ParseArgs():
"available as keys \"per_interface_ips_precision\", "
"\"per_interface_ips_recall\" and \"per_interface_ips\"."))
+ parser.add_argument(
+ "--ips-trimmed",
+ dest="ips_trimmed",
+ default=False,
+ action="store_true",
+ help=("The IPS equivalent of ICS on trimmed models."))
parser.add_argument(
"--rigid-scores",
@@ -451,7 +506,10 @@ def _ParseArgs():
"8.0] given these positions and transformation, \"oligo_gdtha\": "
"same with thresholds [0.5, 1.0, 2.0, 4.0], \"rmsd\": RMSD given "
"these positions and transformation, \"transform\": the used 4x4 "
- "transformation matrix that superposes model onto reference."))
+ "transformation matrix that superposes model onto reference, "
+ "\"rigid_chain_mapping\": equivalent of \"chain_mapping\" which "
+ "is used for rigid scores (optimized for RMSD instead of "
+ "QS-score/lDDT)."))
parser.add_argument(
"--patch-scores",
@@ -565,7 +623,7 @@ def _ParseArgs():
dest="verbosity",
type=int,
default=2,
- help="Set verbosity level. Defaults to 3 (Script).")
+ help="Set verbosity level. Defaults to 2 (Script).")
parser.add_argument(
"--lddt-add-mdl-contacts",
@@ -581,6 +639,14 @@ def _ParseArgs():
"not necessarily in the reference. No contact will "
"be added if the respective atom pair is not "
"resolved in the reference."))
+
+ parser.add_argument(
+ "--lddt-inclusion-radius",
+ dest="lddt_inclusion_radius",
+ type = float,
+ default=15.0,
+ help=("Passed to lDDT scorer. Affects all lDDT scores but not "
+ "chain mapping."))
return parser.parse_args()
@@ -706,6 +772,17 @@ def _LocalScoresToJSONDict(score_dict):
json_dict[f"{ch}.{num.num}.{ins_code}"] = _RoundOrNone(s)
return json_dict
+def _LocalAAScoresToJSONDict(score_dict):
+ """ Convert ResNums and atom names to str for JSON serialization
+ """
+ json_dict = dict()
+ for ch, ch_scores in score_dict.items():
+ for num, res_scores in ch_scores.items():
+ ins_code = num.ins_code.strip("\u0000")
+ for a, s in res_scores.items():
+ json_dict[f"{ch}.{num.num}.{ins_code}.{a}"] = _RoundOrNone(s)
+ return json_dict
+
def _InterfaceResiduesToJSONList(interface_dict):
""" Convert ResNums to str for JSON serialization.
@@ -780,7 +857,8 @@ def _Process(model, reference, args, model_format, reference_format):
min_pep_length = args.min_pep_length,
min_nuc_length = args.min_nuc_length,
lddt_add_mdl_contacts = args.lddt_add_mdl_contacts,
- dockq_capri_peptide = args.dockq_capri_peptide)
+ dockq_capri_peptide = args.dockq_capri_peptide,
+ lddt_inclusion_radius = args.lddt_inclusion_radius)
ir = _GetInconsistentResidues(scorer.aln)
if len(ir) > 0 and args.enforce_consistency:
@@ -810,7 +888,10 @@ def _Process(model, reference, args, model_format, reference_format):
if args.local_lddt:
out["local_lddt"] = _LocalScoresToJSONDict(scorer.local_lddt)
- if args.lddt or args.local_lddt:
+ if args.aa_local_lddt:
+ out["aa_local_lddt"] = _LocalAAScoresToJSONDict(scorer.aa_local_lddt)
+
+ if args.lddt or args.local_lddt or args.aa_local_lddt:
out["model_clashes"] = [x.ToJSON() for x in scorer.model_clashes]
out["model_bad_bonds"] = [x.ToJSON() for x in scorer.model_bad_bonds]
out["model_bad_angles"] = [x.ToJSON() for x in scorer.model_bad_angles]
@@ -845,8 +926,13 @@ def _Process(model, reference, args, model_format, reference_format):
[_RoundOrNone(x) for x in scorer.per_interface_qs_best]
if args.ics or args.ips:
- out["reference_contacts"] = scorer.native_contacts
out["model_contacts"] = scorer.model_contacts
+
+ if args.ics_trimmed or args.ips_trimmed:
+ out["model_contacts_trimmed"] = scorer.trimmed_model_contacts
+
+ if args.ics or args.ips or args.ics_trimmed or args.ips_trimmed:
+ out["reference_contacts"] = scorer.native_contacts
out["contact_reference_interfaces"] = scorer.contact_target_interfaces
if args.ics:
@@ -871,6 +957,28 @@ def _Process(model, reference, args, model_format, reference_format):
out["per_interface_ips"] = \
[_RoundOrNone(x) for x in scorer.per_interface_ips]
+ if args.ics_trimmed:
+ out["ics_trimmed"] = _RoundOrNone(scorer.ics_trimmed)
+ out["ics_precision_trimmed"] = _RoundOrNone(scorer.ics_precision_trimmed)
+ out["ics_recall_trimmed"] = _RoundOrNone(scorer.ics_recall_trimmed)
+ out["per_interface_ics_precision_trimmed"] = \
+ [_RoundOrNone(x) for x in scorer.per_interface_ics_precision_trimmed]
+ out["per_interface_ics_recall_trimmed"] = \
+ [_RoundOrNone(x) for x in scorer.per_interface_ics_recall_trimmed]
+ out["per_interface_ics_trimmed"] = \
+ [_RoundOrNone(x) for x in scorer.per_interface_ics_trimmed]
+
+ if args.ips_trimmed:
+ out["ips_trimmed"] = _RoundOrNone(scorer.ips_trimmed)
+ out["ips_precision_trimmed"] = _RoundOrNone(scorer.ips_precision_trimmed)
+ out["ips_recall_trimmed"] = _RoundOrNone(scorer.ips_recall_trimmed)
+ out["per_interface_ips_precision_trimmed"] = \
+ [_RoundOrNone(x) for x in scorer.per_interface_ips_precision_trimmed]
+ out["per_interface_ips_recall_trimmed"] = \
+ [_RoundOrNone(x) for x in scorer.per_interface_ips_recall_trimmed]
+ out["per_interface_ips_trimmed"] = \
+ [_RoundOrNone(x) for x in scorer.per_interface_ips_trimmed]
+
if args.dockq:
out["dockq_reference_interfaces"] = scorer.dockq_target_interfaces
out["dockq_interfaces"] = scorer.dockq_interfaces
@@ -890,8 +998,9 @@ def _Process(model, reference, args, model_format, reference_format):
out["oligo_gdtts"] = _RoundOrNone(scorer.gdtts)
out["oligo_gdtha"] = _RoundOrNone(scorer.gdtha)
out["rmsd"] = _RoundOrNone(scorer.rmsd)
- data = scorer.transform.data
+ data = scorer.rigid_transform.data
out["transform"] = [data[i:i + 4] for i in range(0, len(data), 4)]
+ out["rigid_chain_mapping"] = scorer.rigid_mapping.GetFlatMapping()
if args.patch_scores:
out["model_interface_residues"] = \
@@ -924,6 +1033,8 @@ def _Main():
args = _ParseArgs()
ost.PushVerbosityLevel(args.verbosity)
if args.verbosity < 4:
+ # Hide tracebacks by default
+ # Run script with -v 4 (Verbose) or higher to display them
sys.tracebacklimit = 0
_CheckCompoundLib()
try:
@@ -960,14 +1071,16 @@ def _Main():
out["min_pep_length"] = args.min_pep_length
out["min_nuc_length"] = args.min_nuc_length
out["lddt_add_mdl_contacts"] = args.lddt_add_mdl_contacts
+ out["lddt_inclusion_radius"] = args.lddt_inclusion_radius
out["dockq_capri_peptide"] = args.dockq_capri_peptide
+ out["ost_version"] = ost.__version__
out["status"] = "SUCCESS"
with open(args.output, 'w') as fh:
json.dump(out, fh, indent=4, sort_keys=False)
except Exception as exc:
out = dict()
out["status"] = "FAILURE"
- out["traceback"] = traceback.format_exc()
+ out["traceback"] = traceback.format_exc(limit=1000)
etype, evalue, tb = sys.exc_info()
out["exception"] = " ".join(traceback.format_exception_only(etype, evalue))
with open(args.output, 'w') as fh:
diff --git a/docker/Dockerfile b/docker/Dockerfile
index b7bbf68ab362aed955d0279f68455b31aadcc90a..ec5a420b47cdf07e2c4821220cf84d78133eb73c 100644
--- a/docker/Dockerfile
+++ b/docker/Dockerfile
@@ -2,12 +2,10 @@ FROM ubuntu:22.04
# ARGUMENTS
###########
-ARG OPENSTRUCTURE_VERSION="2.8.0"
+ARG OPENSTRUCTURE_VERSION="2.9.0"
ARG SRC_FOLDER="/usr/local/src"
ARG CPUS_FOR_MAKE=2
-ARG OPENMM_VERSION="7.7.0"
-ARG OPENMM_INCLUDE_PATH="/usr/local/openmm/include/"
-ARG OPENMM_LIB_PATH="/usr/local/openmm/lib/"
+
ARG DEBIAN_FRONTEND=noninteractive
# INSTALL SYSTEM DEPS
@@ -29,29 +27,15 @@ RUN apt-get update -y && apt-get install -y cmake \
python3-scipy \
python3-pandas \
python3-networkx \
- doxygen \
- swig \
clustalw \
- cython3 \
voronota \
+ libopenmm-dev \
+ libopenmm-plugins \
+ libparasail-dev \
locales && \
# CLEANUP
rm -rf /var/lib/apt/lists/*
-
-# INSTALL OPENMM
-################
-RUN cd ${SRC_FOLDER} && \
- wget -O openmm-${OPENMM_VERSION}.tar.gz -nc https://github.com/pandegroup/openmm/archive/${OPENMM_VERSION}.tar.gz && \
- mkdir ${SRC_FOLDER}/openmm-${OPENMM_VERSION} && \
- tar xf openmm-${OPENMM_VERSION}.tar.gz -C ${SRC_FOLDER}/openmm-${OPENMM_VERSION} --strip-components=1 && \
- mkdir -p ${SRC_FOLDER}/openmm-${OPENMM_VERSION}/build && \
- cd ${SRC_FOLDER}/openmm-${OPENMM_VERSION}/build && \
- cmake .. && make -j $CPUS_FOR_MAKE && make install && \
- cd ${SRC_FOLDER}/openmm-${OPENMM_VERSION}/build/python && \
- python3 setup.py build && python3 setup.py install && \
- rm ${SRC_FOLDER}/openmm-${OPENMM_VERSION}.tar.gz && \
- rm -rf ${SRC_FOLDER}/openmm-${OPENMM_VERSION}
-
+
# INSTALL OST
#############
RUN cd ${SRC_FOLDER} && \
@@ -61,20 +45,17 @@ RUN cd ${SRC_FOLDER} && \
tar xf openstructure-${OPENSTRUCTURE_VERSION}.tar.gz -C ${SRC_FOLDER}/openstructure-${OPENSTRUCTURE_VERSION} --strip-components=1 && \
mkdir -p ${SRC_FOLDER}/openstructure-${OPENSTRUCTURE_VERSION}/build && \
cd ${SRC_FOLDER}/openstructure-${OPENSTRUCTURE_VERSION}/build && \
- cmake .. -DOPTIMIZE=ON \
- -DENABLE_MM=ON \
+ cmake .. -DOPTIMIZE=1 \
+ -DENABLE_MM=1 \
+ -DOPEN_MM_PLUGIN_DIR=/lib/x86_64-linux-gnu/openmm/plugins \
+ -DENABLE_PARASAIL=1 \
-DCOMPILE_TMTOOLS=1 \
- -DOPEN_MM_LIBRARY=$OPENMM_LIB_PATH/libOpenMM.so \
- -DOPEN_MM_INCLUDE_DIR=$OPENMM_INCLUDE_PATH \
- -DOPEN_MM_PLUGIN_DIR=$OPENMM_LIB_PATH/plugins \
- -DENABLE_GFX=ON \
- -DENABLE_GUI=OFF \
- -DENABLE_INFO=OFF \
- -DCMAKE_C_FLAGS="-isystem /usr/include/boost/ -isystem ${OPENMM_INCLUDE_PATH}/include" \
- -DCMAKE_CXX_FLAGS="-isystem /usr/include/boost/ -isystem ${OPENMM_INCLUDE_PATH}/include" && \
+ -DENABLE_GFX=1 \
+ -DENABLE_GUI=0 \
+ -DENABLE_INFO=0 && \
make -j ${CPUS_FOR_MAKE} && \
- wget ftp://ftp.wwpdb.org/pub/pdb/data/monomers/components.cif.gz && \
- stage/bin/chemdict_tool create components.cif.gz compounds.chemlib pdb && stage/bin/chemdict_tool update ../modules/conop/data/charmm.cif compounds.chemlib charmm && \
+ wget https://files.wwpdb.org/pub/pdb/data/monomers/components.cif.gz && \
+ stage/bin/chemdict_tool create components.cif.gz compounds.chemlib pdb -i && stage/bin/chemdict_tool update ../modules/conop/data/charmm.cif compounds.chemlib charmm && \
cmake .. -DCOMPOUND_LIB=${SRC_FOLDER}/openstructure-${OPENSTRUCTURE_VERSION}/build/compounds.chemlib && \
make -j ${CPUS_FOR_MAKE} && make check && make install && \
rm ${SRC_FOLDER}/openstructure-${OPENSTRUCTURE_VERSION}.tar.gz && \
diff --git a/modules/base/src/test_utils/compare_files.cc b/modules/base/src/test_utils/compare_files.cc
index 265914ca1036d08e45c9b957656da7b184206e6c..098bf41b504c0b2bc1caa8f6eef504b839bd24d1 100644
--- a/modules/base/src/test_utils/compare_files.cc
+++ b/modules/base/src/test_utils/compare_files.cc
@@ -22,7 +22,7 @@
#include "compare_files.hh"
namespace ost {
-bool compare_files(const String& test, const String& gold_standard)
+bool compare_files(const String& test, const String& gold_standard, const std::unordered_set<int>& ignore_line)
{
std::ifstream test_stream(test.c_str());
if (!test_stream) {
@@ -36,6 +36,7 @@ bool compare_files(const String& test, const String& gold_standard)
return false;
}
String test_line, gold_line;
+ int line_num = 1;
while (true) {
bool test_read = static_cast<bool>(std::getline(test_stream, test_line));
bool gold_read = static_cast<bool>(std::getline(gold_stream, gold_line));
@@ -54,10 +55,15 @@ bool compare_files(const String& test, const String& gold_standard)
return false;
}
if (gold_line!=test_line) {
+ if (ignore_line.find(line_num) != ignore_line.end()) {
+ continue;
+ }
std::cerr << "line mismatch:" << std::endl << "test: " << test_line
- << std::endl << "gold: " << gold_line << std::endl;
+ << std::endl << "gold: " << gold_line << std::endl
+ << "line: " << line_num << std::endl;
return false;
}
+ ++line_num;
}
return true;
}
diff --git a/modules/base/src/test_utils/compare_files.hh b/modules/base/src/test_utils/compare_files.hh
index 57ab1f6c21e02d5c9b4c9df91663c84425c4f7e9..f8dad2d095fc6c3194ae12d9f74dd44c0d3da762 100644
--- a/modules/base/src/test_utils/compare_files.hh
+++ b/modules/base/src/test_utils/compare_files.hh
@@ -22,11 +22,13 @@
#include <ost/base.hh>
#include <ost/module_config.hh>
+#include <unordered_set>
namespace ost {
-bool DLLEXPORT_OST_BASE compare_files(const String& test,
- const String& gold_standard);
+bool DLLEXPORT_OST_BASE compare_files(const String& test,
+ const String& gold_standard,
+ const std::unordered_set<int>& ignore_line = {});
}
diff --git a/modules/conop/doc/connectivity.rst b/modules/conop/doc/connectivity.rst
index e379494601a31f27a15673f6b7992e140102c287..e18baf87427bab4bb68c5cddf520a80a70dd0424 100644
--- a/modules/conop/doc/connectivity.rst
+++ b/modules/conop/doc/connectivity.rst
@@ -40,14 +40,19 @@ Processors
The exact behaviour for a processor is implementation-specific. So far, two
classes implement the processor interface: A heuristic and a rule-based
processor. The processors mainly differ in the source of their connectivity
-information. The `HeuristicProcessor` uses a hard-coded heuristic connectivity
+information.
+
+The `HeuristicProcessor` uses a hard-coded heuristic connectivity
table for the 20 standard amino acids as well as nucleotides. For other
compounds such as ligands the `HeuristicProcessor` runs a distance-based
-connectivity algorithm that connects two atoms if they are closer than a certain
-threshold. The `RuleBasedProcessor` uses the
-:doc:`compound library <compoundlib>`, a connectivity library containing all
-molecular components present in the PDB files on PDB.org. The library can easily
-be extended with custom connectivity information, if required.
+connectivity algorithm that connects two atoms if they belong to the same or
+two consecutive residues, and are within a
+:func:`reasonable distance <ost.conop.IsBondFeasible>` of each other.
+
+The `RuleBasedProcessor` uses the :doc:`compound library <compoundlib>`,
+a connectivity library containing all molecular components present in the
+PDB files on PDB.org. The library can easily be extended with custom
+connectivity information, if required.
.. class:: Processor
diff --git a/modules/conop/doc/functions.rst b/modules/conop/doc/functions.rst
index 7b3618a9f3bc15ba5d5dcdc4318bc1d0f5096012..456ab200f733bfb565d0f167d46c86abbd76307d 100644
--- a/modules/conop/doc/functions.rst
+++ b/modules/conop/doc/functions.rst
@@ -33,9 +33,12 @@ Conop Functions
.. function:: IsBondFeasible(atom_a, atom_b)
- :return: True, if *atom_a* and *atom_b* are within a reasonable distance for a
- bond to be present. Depends on :attr:`~ost.mol.AtomHandle.radius` of
- atoms and heuristic formulas.
+ :return: True, if *atom_a* and *atom_b* are within a reasonable distance for
+ a bond to be present, namely if the distance between the two atoms
+ is between 0.0625 and 0.375 times the square of the sum of the
+ :attr:`van der Waals radii <ost.mol.AtomHandle.radius>` of the two
+ atom.
+
:rtype: :class:`bool`
:param atom_a: Atom to be checked.
:type atom_a: :class:`~ost.mol.AtomHandle`
diff --git a/modules/conop/src/chemdict_tool.cc b/modules/conop/src/chemdict_tool.cc
index 28e0a2d703d1918e873ac4ca9095c24d801ff132..08b11ae0d931f6c450b8647f8d224304ea04c89e 100644
--- a/modules/conop/src/chemdict_tool.cc
+++ b/modules/conop/src/chemdict_tool.cc
@@ -45,6 +45,7 @@ void PrintUsage()
std::cout << " -i - ignore compounds reserved by the PDB (01-99, DRG, INH, LIG)" << std::endl;
std::cout << " -o - ignore obsolete compounds" << std::endl;
std::cout << " -v - be more verbose" << std::endl;
+ std::cout << " -q - be more quiet (last of -v or -q is applied)" << std::endl;
}
int main(int argc, char const *argv[])
@@ -74,9 +75,15 @@ int main(int argc, char const *argv[])
ignore_obsolete=true;
} else if (param=="-v") {
Logger::Instance().PushVerbosityLevel(4);
- } else {
+ } else if (param=="-q") {
+ Logger::Instance().PushVerbosityLevel(0);
+ } else if (param=="-h") {
PrintUsage();
return 0;
+ } else {
+ std::cout << "Unrecognized argument '" << param << "'" << std::endl;
+ std::cout << "Try 'chemdict_tool -h' for more information" << std::endl;
+ return 1;
}
}
boost::iostreams::filtering_stream<boost::iostreams::input> filtered_istream;
@@ -89,34 +96,45 @@ int main(int argc, char const *argv[])
if (boost::iequals(".gz", boost::filesystem::extension(argv[2]))) {
filtered_istream.push(boost::iostreams::gzip_decompressor());
}
- filtered_istream.push(istream);
- io::ChemdictParser cdp(filtered_istream, dialect, ignore_reserved, ignore_obsolete);
- conop::CompoundLibPtr compound_lib;
- bool in_mem=false;
- if (!strncmp(argv[1], "create", 6)) {
- compound_lib=conop::CompoundLib::Create(":memory:");
- in_mem=true;
- } else if (!strncmp(argv[1], "update", 6)) {
- compound_lib=conop::CompoundLib::Load(argv[3]);
- } else {
- PrintUsage();
- return 0;
- }
- if (!compound_lib) {
- return 0;
- }
- assert(compound_lib);
- conop::CompoundLibPtr in_mem_lib=in_mem ? compound_lib :
- compound_lib->Copy(":memory:");
- compound_lib.reset();
- cdp.SetCompoundLib(in_mem_lib);
- cdp.Parse();
- in_mem_lib->SetChemLibInfo();
- conop::CompoundLibPtr copy = in_mem_lib->Copy(argv[3]);
- if (! copy) {
- std::cout << "Cannot save " << argv[3] << ": [Errno " << errno << "] "
- << strerror(errno) << std::endl;
+ filtered_istream.push(istream);
+ try {
+ io::ChemdictParser cdp(filtered_istream, dialect, ignore_reserved, ignore_obsolete);
+ conop::CompoundLibPtr compound_lib;
+ bool in_mem=false;
+ if (!strncmp(argv[1], "create", 6)) {
+ compound_lib=conop::CompoundLib::Create(":memory:");
+ in_mem=true;
+ } else if (!strncmp(argv[1], "update", 6)) {
+ compound_lib=conop::CompoundLib::Load(argv[3]);
+ } else {
+ PrintUsage();
+ return 0;
+ }
+ if (!compound_lib) {
+ return 0;
+ }
+ assert(compound_lib);
+ conop::CompoundLibPtr in_mem_lib=in_mem ? compound_lib :
+ compound_lib->Copy(":memory:");
+ compound_lib.reset();
+ cdp.SetCompoundLib(in_mem_lib);
+ cdp.Parse();
+ if (cdp.GetImportedCount() == 0)
+ {
+ std::cout << "No compound imported from " << argv[2] << std::endl;
return 1;
+ }
+ in_mem_lib->SetChemLibInfo();
+ conop::CompoundLibPtr copy = in_mem_lib->Copy(argv[3]);
+ if (! copy) {
+ std::cout << "Cannot save " << argv[3] << ": [Errno " << errno << "] "
+ << strerror(errno) << std::endl;
+ return 1;
+ }
+ }
+ catch (ost::Error const &err) {
+ std::cerr << "Error: " << err.what() << std::endl;
+ return 1;
}
return 0;
}
diff --git a/modules/conop/src/compound.hh b/modules/conop/src/compound.hh
index 1b0eb1cd5e54f80e7d0bcab44a7b88ca0c36c6d2..07916c231b868266dce8642d85d7e94172c36147 100644
--- a/modules/conop/src/compound.hh
+++ b/modules/conop/src/compound.hh
@@ -23,6 +23,8 @@
#include <map>
#include <boost/shared_ptr.hpp>
#include <ost/string_ref.hh>
+#include <ost/message.hh>
+#include <ost/log.hh>
#include <ost/conop/module_config.hh>
#include <ost/mol/chem_class.hh>
@@ -54,19 +56,27 @@ struct DLLEXPORT_OST_CONOP Date {
static Date FromString(const StringRef& str)
{
std::vector<StringRef> parts=str.split('-');
- assert(parts.size()==3);
+ if (parts.size() != 3) {
+ std::stringstream msg;
+ msg << "Invalid date string: '" << str << "'";
+ throw ost::Error(msg.str());
+ }
std::pair<bool, int> year=parts[0].to_int();
std::pair<bool, int> month=parts[1].to_int();
std::pair<bool, int> day=parts[2].to_int();
- assert(year.first); assert(month.first); assert(day.first);
+ if (! year.first || ! month.first || ! day.first) {
+ std::stringstream msg;
+ msg << "Invalid date string: '" << str << "'";
+ throw ost::Error(msg.str());
+ }
return Date(year.second, month.second, day.second);
}
-
+
String ToString() const;
int year;
int month;
- int day;
+ int day;
};
struct DLLEXPORT_OST_CONOP AtomSpec {
@@ -169,6 +179,12 @@ public:
const String& GetID() const {
return tlc_;
}
+
+ /// \brief set three-letter code that is unique for every compound
+ void SetID(const String& id) {
+ tlc_ = id;
+ }
+
Dialect GetDialect() const { return dialect_; }
String GetDialectAsString() const {
diff --git a/modules/conop/src/compound_lib.cc b/modules/conop/src/compound_lib.cc
index c84b0e9b7837036a66c737fc86c93d9d929bcf93..6782828f33bf2cca80cefed5c47586562778f9ea 100644
--- a/modules/conop/src/compound_lib.cc
+++ b/modules/conop/src/compound_lib.cc
@@ -264,13 +264,23 @@ void CompoundLib::AddCompound(const CompoundPtr& compound)
sqlite3_bind_text(stmt, 6, compound->GetFormula().c_str(),
compound->GetFormula().length(), NULL);
Date crea_date=compound->GetCreationDate();
+ if (crea_date == Date()) {
+ sqlite3_bind_null(stmt, 7);
+ }
+ else {
+ crea_date_str=crea_date.ToString();
+ sqlite3_bind_text(stmt, 7, crea_date_str.c_str(), crea_date_str.length(),
+ NULL);
+ }
Date modi_date=compound->GetModificationDate();
- crea_date_str=crea_date.ToString();
- modi_date_str=modi_date.ToString();
- sqlite3_bind_text(stmt, 7, crea_date_str.c_str(), crea_date_str.length(),
- NULL);
- sqlite3_bind_text(stmt, 8, modi_date_str.c_str(), modi_date_str.length(),
- NULL);
+ if (modi_date == Date()) {
+ sqlite3_bind_null(stmt, 8);
+ }
+ else {
+ modi_date_str=modi_date.ToString();
+ sqlite3_bind_text(stmt, 8, modi_date_str.c_str(), modi_date_str.length(),
+ NULL);
+ }
sqlite3_bind_text(stmt, 9, compound->GetName().c_str(),
compound->GetName().length(), NULL);
sqlite3_bind_text(stmt, 10, compound->GetInchi().c_str(),
diff --git a/modules/conop/tests/test_complib.py b/modules/conop/tests/test_complib.py
index f0138109fce400506c1f96150340b66ffc810773..f22332ae8a59be153204d966999f1976d04e98c3 100644
--- a/modules/conop/tests/test_complib.py
+++ b/modules/conop/tests/test_complib.py
@@ -11,10 +11,10 @@ def CreateComplib(compound_dict_path, chemlib_out_path, extra_args=None):
chemdict_tool_path = os.path.join(prefix_path, "bin", "chemdict_tool")
if not os.path.exists(chemdict_tool_path):
raise RuntimeError("Expect chemdict_tool:", chemdict_tool_path)
- cmd = [chemdict_tool_path, "create", compound_dict_path, chemlib_out_path]
+ cmd = [chemdict_tool_path, "create", compound_dict_path, chemlib_out_path, "-q"]
if extra_args:
cmd += extra_args
- subprocess.run(cmd, stdout=subprocess.PIPE)
+ subprocess.run(cmd, stdout=subprocess.PIPE, check=True)
class TestCompLib(unittest.TestCase):
@@ -118,6 +118,32 @@ class TestCompLib(unittest.TestCase):
# OX is obsolete
assert complib_no_obsolete.FindCompound("OX") is None
+ def test_bird(self):
+ """ the test file in test_compounds.cif contains two entries from
+ BIRD:
+
+ - PRD_900116 is not in a compound format and should not be added
+ to the compound lib.
+ - PRDCC_900116 is in a valid format and should be read in
+ """
+ complib = self.complib
+ # PRD_900116 not added
+ comp = complib.FindCompound("PRD_900116")
+ self.assertIsNone(comp)
+
+ # PRD_900116 is added
+ comp = complib.FindCompound("PRDCC_900116")
+ self.assertIsNotNone(comp)
+ self.assertEqual(
+ comp.smiles,
+ "C1[C@H]([C@@H]([C@H]([C@@H](O1)O[C@@H]2CO[C@H]([C@@H]([C@H]2O)O)O)O)O)O")
+ # The release flag is REF_ONLY and compound should not be marked as
+ # obsolete
+ self.assertFalse(comp.obsolete)
+ # Contains atoms and bonds
+ self.assertEqual(len(comp.atom_specs), 37)
+ self.assertEqual(len(comp.bond_specs), 38)
+
if __name__ == "__main__":
from ost import testutils
diff --git a/modules/conop/tests/testfiles/test_compounds.cif b/modules/conop/tests/testfiles/test_compounds.cif
index 340a904a1cabf1ff4e8af63557a9e8ec37ee29e1..bb9f6423c7648d13fc908f77cde9afd38d17968a 100644
--- a/modules/conop/tests/testfiles/test_compounds.cif
+++ b/modules/conop/tests/testfiles/test_compounds.cif
@@ -2158,3 +2158,240 @@ _chem_comp.pdbx_model_coordinates_db_code ?
_chem_comp.pdbx_subcomponent_list ?
_chem_comp.pdbx_processing_site RCSB
##
+
+# PRDCC_900116 is a BIRD entry in PRDCC format and should read in correctly in
+# the compound lib
+data_PRDCC_900116
+#
+_chem_comp.id PRDCC_900116
+_chem_comp.name 4-O-beta-D-xylopyranosyl-beta-D-xylopyranose
+_chem_comp.type saccharide
+_chem_comp.formula "C10 H18 O9"
+_chem_comp.pdbx_release_status REF_ONLY
+_chem_comp.formula_weight 282.245
+_chem_comp.pdbx_type ?
+_chem_comp.mon_nstd_parent_comp_id ?
+_chem_comp.pdbx_synonyms ?
+_chem_comp.pdbx_formal_charge ?
+_chem_comp.pdbx_initial_date ?
+_chem_comp.pdbx_modified_date ?
+_chem_comp.pdbx_ambiguous_flag ?
+_chem_comp.pdbx_replaced_by ?
+_chem_comp.pdbx_replaces ?
+_chem_comp.one_letter_code ?
+_chem_comp.three_letter_code ?
+_chem_comp.pdbx_model_coordinates_details ?
+_chem_comp.pdbx_model_coordinates_missing_flag ?
+_chem_comp.pdbx_ideal_coordinates_details ?
+_chem_comp.pdbx_ideal_coordinates_missing_flag ?
+_chem_comp.pdbx_model_coordinates_db_code ?
+_chem_comp.pdbx_subcomponent_list ?
+_chem_comp.pdbx_processing_site ?
+#
+loop_
+_chem_comp_atom.comp_id
+_chem_comp_atom.atom_id
+_chem_comp_atom.alt_atom_id
+_chem_comp_atom.type_symbol
+_chem_comp_atom.charge
+_chem_comp_atom.pdbx_align
+_chem_comp_atom.pdbx_aromatic_flag
+_chem_comp_atom.pdbx_leaving_atom_flag
+_chem_comp_atom.pdbx_stereo_config
+_chem_comp_atom.model_Cartn_x
+_chem_comp_atom.model_Cartn_y
+_chem_comp_atom.model_Cartn_z
+_chem_comp_atom.pdbx_model_Cartn_x_ideal
+_chem_comp_atom.pdbx_model_Cartn_y_ideal
+_chem_comp_atom.pdbx_model_Cartn_z_ideal
+_chem_comp_atom.pdbx_component_comp_id
+_chem_comp_atom.pdbx_residue_numbering
+_chem_comp_atom.pdbx_component_atom_id
+_chem_comp_atom.pdbx_polymer_type
+_chem_comp_atom.pdbx_ref_id
+_chem_comp_atom.pdbx_component_id
+_chem_comp_atom.pdbx_ordinal
+_chem_comp_atom.pdbx_backbone_atom_flag
+_chem_comp_atom.pdbx_n_terminal_atom_flag
+_chem_comp_atom.pdbx_c_terminal_atom_flag
+PRD_900116 O1A O9 O 0 1 N N N 62.369 63.326 22.057 -5.145 -0.343 1.296 XYP 1 O4A polymer 1 1 1 ? ? ?
+PRD_900116 C1A C9 C 0 1 N N R 62.245 64.245 21.048 -3.966 -0.610 0.535 XYP 1 C1B polymer 1 1 2 ? ? ?
+PRD_900116 C2A C10 C 0 1 N N R 63.637 64.692 20.579 -3.534 0.662 -0.201 XYP 1 C2B polymer 1 1 3 ? ? ?
+PRD_900116 C3A C7 C 0 1 N N R 63.577 65.572 19.393 -2.237 0.384 -0.966 XYP 1 C3B polymer 1 1 4 ? ? ?
+PRD_900116 C4A C6 C 0 1 N N R 62.774 65.007 18.298 -1.171 -0.108 0.019 XYP 1 C4B polymer 1 1 5 ? ? ?
+PRD_900116 C5A C8 C 0 1 N N N 61.391 64.538 18.774 -1.694 -1.350 0.746 XYP 1 C5B polymer 1 1 6 ? ? ?
+PRD_900116 O2A O8 O 0 1 N N N 64.259 65.423 21.647 -4.556 1.052 -1.120 XYP 1 O2B polymer 1 1 7 ? ? ?
+PRD_900116 O3A O6 O 0 1 N N N 64.944 65.808 18.894 -1.787 1.585 -1.596 XYP 1 O3B polymer 1 1 8 ? ? ?
+PRD_900116 O4A O5 O 0 1 N N N 62.601 66.032 17.300 0.022 -0.440 -0.695 XYP 1 O4B polymer 1 1 9 ? ? ?
+PRD_900116 O5A O7 O 0 1 N N N 61.492 63.643 19.942 -2.918 -1.031 1.411 XYP 1 O5B polymer 1 1 10 ? ? ?
+PRD_900116 C1B C5 C 0 1 N N S 62.805 65.537 16.116 1.220 -0.261 0.063 XYP 2 C1B polymer 1 1 11 ? ? ?
+PRD_900116 C2B C3 C 0 1 N N R 62.263 66.524 15.005 2.417 -0.761 -0.750 XYP 2 C2B polymer 1 1 12 ? ? ?
+PRD_900116 C3B C2 C 0 1 N N S 62.560 66.033 13.636 3.703 -0.511 0.043 XYP 2 C3B polymer 1 1 13 ? ? ?
+PRD_900116 C4B C1 C 0 1 N N R 63.947 65.755 13.420 3.806 0.981 0.372 XYP 2 C4B polymer 1 1 14 ? ? ?
+PRD_900116 C5B C4 C 0 1 N N N 64.542 64.882 14.493 2.556 1.414 1.142 XYP 2 C5B polymer 1 1 15 ? ? ?
+PRD_900116 O2B O3 O 0 1 N N N 60.863 66.689 15.161 2.274 -2.161 -1.000 XYP 2 O2B polymer 1 1 16 ? ? ?
+PRD_900116 O3B O2 O 0 1 N N N 62.135 67.086 12.661 4.832 -0.907 -0.738 XYP 2 O3B polymer 1 1 17 ? ? ?
+PRD_900116 O4B O1 O 0 1 N N N 64.114 65.051 12.120 4.965 1.215 1.174 XYP 2 O4B polymer 1 1 18 ? ? ?
+PRD_900116 O5B O4 O 0 1 N N N 64.275 65.381 15.859 1.394 1.126 0.361 XYP 2 O5B polymer 1 1 19 ? ? ?
+PRD_900116 HO1A H18 H 0 0 N N N 62.856 63.709 22.777 -5.476 -1.105 1.790 XYP 1 HO4A polymer 1 1 20 ? ? ?
+PRD_900116 H1A1 H15 H 0 0 N N N 61.704 65.134 21.405 -4.172 -1.397 -0.191 XYP 1 H1B polymer 1 1 21 ? ? ?
+PRD_900116 H2A H16 H 0 1 N N N 64.231 63.798 20.340 -3.368 1.461 0.521 XYP 1 H2B polymer 1 1 22 ? ? ?
+PRD_900116 H3A H11 H 0 1 N N N 63.139 66.536 19.692 -2.415 -0.381 -1.722 XYP 1 H3B polymer 1 1 23 ? ? ?
+PRD_900116 H4A H10 H 0 1 N N N 63.304 64.147 17.862 -0.956 0.675 0.745 XYP 1 H4B polymer 1 1 24 ? ? ?
+PRD_900116 H5A1 H14 H 0 0 N N N 60.895 63.999 17.953 -0.958 -1.679 1.479 XYP 1 H5B1 polymer 1 1 25 ? ? ?
+PRD_900116 H5A2 H13 H 0 0 N N N 60.791 65.418 19.050 -1.870 -2.147 0.023 XYP 1 H5B2 polymer 1 1 26 ? ? ?
+PRD_900116 HO2A H17 H 0 0 N N N 64.308 64.874 22.421 -5.408 1.244 -0.704 XYP 1 HO2B polymer 1 1 27 ? ? ?
+PRD_900116 HO3A H12 H 0 0 N N N 65.477 66.174 19.590 -2.414 1.953 -2.234 XYP 1 HO3B polymer 1 1 28 ? ? ?
+PRD_900116 H1B H9 H 0 1 N N N 62.310 64.564 15.980 1.151 -0.826 0.992 XYP 2 H1B polymer 1 1 29 ? ? ?
+PRD_900116 H2B H5 H 0 1 N N N 62.766 67.492 15.149 2.464 -0.225 -1.698 XYP 2 H2B polymer 1 1 30 ? ? ?
+PRD_900116 H3B H3 H 0 1 N N N 61.970 65.124 13.450 3.679 -1.087 0.969 XYP 2 H3B polymer 1 1 31 ? ? ?
+PRD_900116 H4B H2 H 0 1 N N N 64.510 66.699 13.386 3.882 1.553 -0.553 XYP 2 H4B polymer 1 1 32 ? ? ?
+PRD_900116 H5B1 H8 H 0 0 N N N 65.631 64.836 14.343 2.603 2.484 1.341 XYP 2 H5B1 polymer 1 1 33 ? ? ?
+PRD_900116 H5B2 H7 H 0 0 N N N 64.116 63.872 14.399 2.503 0.870 2.085 XYP 2 H5B2 polymer 1 1 34 ? ? ?
+PRD_900116 HO2B H6 H 0 0 N N N 60.676 67.003 16.038 1.479 -2.392 -1.499 XYP 2 HO2B polymer 1 1 35 ? ? ?
+PRD_900116 HO3B H4 H 0 0 N N N 61.215 67.284 12.789 5.682 -0.771 -0.299 XYP 2 HO3B polymer 1 1 36 ? ? ?
+PRD_900116 HO4B H1 H 0 0 N N N 65.035 64.866 11.976 5.093 2.141 1.420 XYP 2 HO4B polymer 1 1 37 ? ? ?
+#
+loop_
+_chem_comp_bond.comp_id
+_chem_comp_bond.atom_id_1
+_chem_comp_bond.atom_id_2
+_chem_comp_bond.value_order
+_chem_comp_bond.pdbx_aromatic_flag
+_chem_comp_bond.pdbx_stereo_config
+_chem_comp_bond.pdbx_ordinal
+PRD_900116 O4B C4B SING N N 1
+PRD_900116 O3B C3B SING N N 2
+PRD_900116 C4B C3B SING N N 3
+PRD_900116 C4B C5B SING N N 4
+PRD_900116 C3B C2B SING N N 5
+PRD_900116 C5B O5B SING N N 6
+PRD_900116 C2B O2B SING N N 7
+PRD_900116 C2B C1B SING N N 8
+PRD_900116 O5B C1B SING N N 9
+PRD_900116 C1B O4A SING N N 10
+PRD_900116 O4A C4A SING N N 11
+PRD_900116 C4A C5A SING N N 12
+PRD_900116 C4A C3A SING N N 13
+PRD_900116 C5A O5A SING N N 14
+PRD_900116 O3A C3A SING N N 15
+PRD_900116 C3A C2A SING N N 16
+PRD_900116 O5A C1A SING N N 17
+PRD_900116 C2A C1A SING N N 18
+PRD_900116 C2A O2A SING N N 19
+PRD_900116 C1A O1A SING N N 20
+PRD_900116 O4B HO4B SING N N 21
+PRD_900116 C4B H4B SING N N 22
+PRD_900116 C3B H3B SING N N 23
+PRD_900116 O3B HO3B SING N N 24
+PRD_900116 C2B H2B SING N N 25
+PRD_900116 O2B HO2B SING N N 26
+PRD_900116 C5B H5B2 SING N N 27
+PRD_900116 C5B H5B1 SING N N 28
+PRD_900116 C1B H1B SING N N 29
+PRD_900116 C4A H4A SING N N 30
+PRD_900116 C3A H3A SING N N 31
+PRD_900116 O3A HO3A SING N N 32
+PRD_900116 C5A H5A2 SING N N 33
+PRD_900116 C5A H5A1 SING N N 34
+PRD_900116 C1A H1A1 SING N N 35
+PRD_900116 C2A H2A SING N N 36
+PRD_900116 O2A HO2A SING N N 37
+PRD_900116 O1A HO1A SING N N 38
+#
+loop_
+_pdbx_chem_comp_descriptor.comp_id
+_pdbx_chem_comp_descriptor.type
+_pdbx_chem_comp_descriptor.program
+_pdbx_chem_comp_descriptor.program_version
+_pdbx_chem_comp_descriptor.descriptor
+PRD_900116 SMILES ACDLabs 12.01 "OC1C(C(O)C(OC1)OC2COC(O)C(C2O)O)O"
+PRD_900116 InChI InChI 1.03 "InChI=1S/C10H18O9/c11-3-1-18-10(8(15)5(3)12)19-4-2-17-9(16)7(14)6(4)13/h3-16H,1-2H2/t3-,4-,5+,6+,7-,8-,9-,10+/m1/s1"
+PRD_900116 InChIKey InChI 1.03 LGQKSQQRKHFMLI-SJYYZXOBSA-N
+PRD_900116 SMILES_CANONICAL CACTVS 3.385 "O[C@@H]1CO[C@@H](O[C@@H]2CO[C@@H](O)[C@H](O)[C@H]2O)[C@H](O)[C@H]1O"
+PRD_900116 SMILES CACTVS 3.385 "O[CH]1CO[CH](O[CH]2CO[CH](O)[CH](O)[CH]2O)[CH](O)[CH]1O"
+PRD_900116 SMILES_CANONICAL "OpenEye OEToolkits" 2.0.7 "C1[C@H]([C@@H]([C@H]([C@@H](O1)O[C@@H]2CO[C@H]([C@@H]([C@H]2O)O)O)O)O)O"
+PRD_900116 SMILES "OpenEye OEToolkits" 2.0.7 "C1C(C(C(C(O1)OC2COC(C(C2O)O)O)O)O)O"
+#
+loop_
+_pdbx_chem_comp_identifier.comp_id
+_pdbx_chem_comp_identifier.type
+_pdbx_chem_comp_identifier.program
+_pdbx_chem_comp_identifier.program_version
+_pdbx_chem_comp_identifier.identifier
+PRD_900116 "SYSTEMATIC NAME" ACDLabs 12.01 4-O-beta-D-xylopyranosyl-beta-D-xylopyranose
+PRD_900116 "SYSTEMATIC NAME" "OpenEye OEToolkits" 2.0.7 "(2~{S},3~{R},4~{S},5~{R})-2-[(3~{R},4~{R},5~{R},6~{R})-4,5,6-tris(oxidanyl)oxan-3-yl]oxyoxane-3,4,5-triol"
+#
+
+# PRD_900116 is a BIRD entry in PRD format. It doesn't have _chem_comp or anything useful
+# and should not end up in the compound lib
+data_PRD_900116
+#
+_pdbx_reference_molecule.prd_id PRD_900116
+_pdbx_reference_molecule.name 4beta-beta-xylobiose
+_pdbx_reference_molecule.represent_as branched
+_pdbx_reference_molecule.type Oligosaccharide
+_pdbx_reference_molecule.type_evidence_code ?
+_pdbx_reference_molecule.class metabolism
+_pdbx_reference_molecule.class_evidence_code ?
+_pdbx_reference_molecule.formula "C10 H18 O9"
+_pdbx_reference_molecule.formula_weight 282.245
+_pdbx_reference_molecule.chem_comp_id ?
+_pdbx_reference_molecule.release_status REL
+_pdbx_reference_molecule.replaces ?
+_pdbx_reference_molecule.replaced_by ?
+_pdbx_reference_molecule.compound_details ?
+_pdbx_reference_molecule.description oligosaccharide
+_pdbx_reference_molecule.representative_PDB_id_code 6RV8
+#
+_pdbx_reference_entity_list.prd_id PRD_900116
+_pdbx_reference_entity_list.ref_entity_id 1
+_pdbx_reference_entity_list.component_id 1
+_pdbx_reference_entity_list.type branched
+_pdbx_reference_entity_list.details ?
+#
+_pdbx_reference_entity_poly.prd_id PRD_900116
+_pdbx_reference_entity_poly.ref_entity_id 1
+_pdbx_reference_entity_poly.type oligosaccharide
+_pdbx_reference_entity_poly.db_code ?
+_pdbx_reference_entity_poly.db_name ?
+#
+_pdbx_reference_entity_sequence.prd_id PRD_900116
+_pdbx_reference_entity_sequence.ref_entity_id 1
+_pdbx_reference_entity_sequence.type saccharide
+_pdbx_reference_entity_sequence.NRP_flag N
+#
+loop_
+_pdbx_reference_entity_poly_seq.prd_id
+_pdbx_reference_entity_poly_seq.ref_entity_id
+_pdbx_reference_entity_poly_seq.num
+_pdbx_reference_entity_poly_seq.mon_id
+_pdbx_reference_entity_poly_seq.parent_mon_id
+_pdbx_reference_entity_poly_seq.hetero
+_pdbx_reference_entity_poly_seq.observed
+PRD_900116 1 1 XYP ? N Y
+PRD_900116 1 2 XYP ? N Y
+#
+_pdbx_reference_entity_poly_link.prd_id PRD_900116
+_pdbx_reference_entity_poly_link.ref_entity_id 1
+_pdbx_reference_entity_poly_link.component_id 1
+_pdbx_reference_entity_poly_link.link_id 1
+_pdbx_reference_entity_poly_link.atom_id_1 O4
+_pdbx_reference_entity_poly_link.comp_id_1 XYP
+_pdbx_reference_entity_poly_link.entity_seq_num_1 1
+_pdbx_reference_entity_poly_link.atom_id_2 C1
+_pdbx_reference_entity_poly_link.comp_id_2 XYP
+_pdbx_reference_entity_poly_link.entity_seq_num_2 2
+_pdbx_reference_entity_poly_link.value_order SING
+_pdbx_reference_entity_poly_link.details ?
+_pdbx_reference_entity_poly_link.insert_code_1 ?
+_pdbx_reference_entity_poly_link.insert_code_2 ?
+#
+loop_
+_pdbx_prd_audit.prd_id
+_pdbx_prd_audit.date
+_pdbx_prd_audit.processing_site
+_pdbx_prd_audit.action_type
+PRD_900116 2020-05-08 RCSB "Create molecule"
+PRD_900116 2020-07-29 RCSB "Initial release"
+#
\ No newline at end of file
diff --git a/modules/doc/actions.rst b/modules/doc/actions.rst
index 635f9524de9ad7cc098c1682748143f52c8ba51c..56afd04dcda6294dfeec0eefd5d64e1bced8fcd8 100644
--- a/modules/doc/actions.rst
+++ b/modules/doc/actions.rst
@@ -37,20 +37,22 @@ Details on the usage (output of ``ost compare-structures --help``):
[-mb MODEL_BIOUNIT] [-rb REFERENCE_BIOUNIT]
[-rna] [-ec] [-d] [-ds DUMP_SUFFIX] [-ft]
[-c CHAIN_MAPPING [CHAIN_MAPPING ...]] [--lddt]
- [--local-lddt] [--bb-lddt] [--bb-local-lddt]
- [--ilddt] [--cad-score] [--local-cad-score]
- [--cad-exec CAD_EXEC]
+ [--local-lddt] [--aa-local-lddt] [--bb-lddt]
+ [--bb-local-lddt] [--ilddt] [--cad-score]
+ [--local-cad-score] [--cad-exec CAD_EXEC]
[--usalign-exec USALIGN_EXEC]
[--override-usalign-mapping] [--qs-score]
[--dockq] [--dockq-capri-peptide] [--ics]
- [--ips] [--rigid-scores] [--patch-scores]
- [--tm-score] [--lddt-no-stereochecks]
+ [--ics-trimmed] [--ips] [--ips-trimmed]
+ [--rigid-scores] [--patch-scores] [--tm-score]
+ [--lddt-no-stereochecks]
[--n-max-naive N_MAX_NAIVE]
[--dump-aligned-residues] [--dump-pepnuc-alns]
[--dump-pepnuc-aligned-residues]
[--min-pep-length MIN_PEP_LENGTH]
[--min-nuc-length MIN_NUC_LENGTH] [-v VERBOSITY]
[--lddt-add-mdl-contacts]
+ [--lddt-inclusion-radius LDDT_INCLUSION_RADIUS]
Evaluate model against reference
@@ -113,7 +115,9 @@ Details on the usage (output of ``ost compare-structures --help``):
* "min_pep_length"
* "min_nuc_length"
* "lddt_add_mdl_contacts"
+ * "lddt_inclusion_radius"
* "dockq_capri_peptide"
+ * "ost_version"
The pairwise sequence alignments are computed with Needleman-Wunsch using
BLOSUM62 (NUC44 for nucleotides). Many benchmarking scenarios preprocess the
@@ -202,6 +206,19 @@ Details on the usage (output of ``ost compare-structures --help``):
counterparts. Atoms specified in there follow the
following format:
<chain_name>.<resnum>.<resnum_inscode>.<atom_name>
+ --aa-local-lddt Compute per-atom lDDT scores with default
+ parameterization and store as key "aa_local_lddt".
+ Score for each atom is accessible by key
+ <chain_name>.<resnum>.<resnum_inscode>.<aname>. Alpha
+ carbon from residue with number 42 in chain X can be
+ extracted with: data["aa_local_lddt"]["X.42..CA"]. If
+ there is a residue insertion code, lets say A, the
+ atom key becomes "X.42.A.CA". Stereochemical
+ irregularities affecting lDDT are reported as keys
+ "model_clashes", "model_bad_bonds", "model_bad_angles"
+ and the respective reference counterparts. Atoms
+ specified in there follow the following format:
+ <chain_name>.<resnum>.<resnum_inscode>.<atom_name>
--bb-lddt Compute global lDDT score with default
parameterization and store as key "bb_lddt". lDDT in
this case is only computed on backbone atoms: CA for
@@ -252,28 +269,29 @@ Details on the usage (output of ``ost compare-structures --help``):
--dockq Compute DockQ scores and its components. Relevant
interfaces with at least one contact (any atom within
5A) of the reference structure are available as key
- "dockq_reference_interfaces". Only interfaces between
- peptide chains are considered here! Key
- "dockq_interfaces" is a subset of
- "dockq_reference_interfaces" that contains interfaces
- that can be mapped to the model. They are stored as
- lists in format [ref_ch1, ref_ch2, mdl_ch1, mdl_ch2].
- The respective DockQ scores for "dockq_interfaces" are
- available as key "dockq". It's components are
- available as keys: "fnat" (fraction of reference
- contacts which are also there in model) "irmsd"
- (interface RMSD), "lrmsd" (ligand RMSD). The DockQ
- score is strictly designed to score each interface
- individually. We also provide two averaged versions to
- get one full model score: "dockq_ave", "dockq_wave".
- The first is simply the average of "dockq_scores", the
- latter is a weighted average with weights derived from
- number of contacts in the reference interfaces. These
- two scores only consider interfaces that are present
- in both, the model and the reference. "dockq_ave_full"
- and "dockq_wave_full" add zeros in the average
- computation for each interface that is only present in
- the reference but not in the model.
+ "dockq_reference_interfaces". Protein-protein,
+ protein-nucleotide and nucleotide-nucleotide
+ interfaces are considered. Key "dockq_interfaces" is a
+ subset of "dockq_reference_interfaces" that contains
+ interfaces that can be mapped to the model. They are
+ stored as lists in format [ref_ch1, ref_ch2, mdl_ch1,
+ mdl_ch2]. The respective DockQ scores for
+ "dockq_interfaces" are available as key "dockq". It's
+ components are available as keys: "fnat" (fraction of
+ reference contacts which are also there in model)
+ "irmsd" (interface RMSD), "lrmsd" (ligand RMSD). The
+ DockQ score is strictly designed to score each
+ interface individually. We also provide two averaged
+ versions to get one full model score: "dockq_ave",
+ "dockq_wave". The first is simply the average of
+ "dockq_scores", the latter is a weighted average with
+ weights derived from number of contacts in the
+ reference interfaces. These two scores only consider
+ interfaces that are present in both, the model and the
+ reference. "dockq_ave_full" and "dockq_wave_full" add
+ zeros in the average computation for each interface
+ that is only present in the reference but not in the
+ model.
--dockq-capri-peptide
Flag that changes two things in the way DockQ and its
underlying scores are computed which is proposed by
@@ -289,8 +307,9 @@ Details on the usage (output of ``ost compare-structures --help``):
only considering CB atoms for protein-peptide
interactions. Note that the resulting DockQ is not
evaluated for these slightly updated fnat and irmsd
- (lrmsd stays the same).This flag has no influence on
- patch_dockq scores.
+ (lrmsd stays the same). Raises an error if reference
+ contains nucleotide chains. This flag has no influence
+ on patch_dockq scores.
--ics Computes interface contact similarity (ICS) related
scores. A contact between two residues of different
chains is defined as having at least one heavy atom
@@ -314,6 +333,29 @@ Details on the usage (output of ``ost compare-structures --help``):
metrics are available as keys
"per_interface_ics_precision",
"per_interface_ics_recall" and "per_interface_ics".
+ --ics-trimmed Computes interface contact similarity (ICS) related
+ scores but on a trimmed model. That means that a
+ mapping between model and reference is performed and
+ all model residues without reference counterpart are
+ removed. As a consequence, model contacts for which we
+ have no experimental evidence do not affect the score.
+ The effect of these added model contacts without
+ mapping to target would be decreased precision and
+ thus lower ics. Recall is not affected. Enabling this
+ flag adds the following keys: "ics_trimmed",
+ "ics_precision_trimmed", "ics_recall_trimmed",
+ "model_contacts_trimmed". The reference contacts and
+ reference interfaces are the same as for ics and
+ available as keys: "reference_contacts",
+ "contact_reference_interfaces". All these measures are
+ also available on a per-interface basis for each
+ interface in the reference structure that are defined
+ as chain pairs with at least one contact (available as
+ key "contact_reference_interfaces"). The respective
+ metrics are available as keys
+ "per_interface_ics_precision_trimmed",
+ "per_interface_ics_recall_trimmed" and
+ "per_interface_ics_trimmed".
--ips Computes interface patch similarity (IPS) related
scores. They focus on interface residues. They are
defined as having at least one contact to a residue
@@ -337,6 +379,7 @@ Details on the usage (output of ``ost compare-structures --help``):
metrics are available as keys
"per_interface_ips_precision",
"per_interface_ips_recall" and "per_interface_ips".
+ --ips-trimmed The IPS equivalent of ICS on trimmed models.
--rigid-scores Computes rigid superposition based scores. They're
based on a Kabsch superposition of all mapped CA
positions (C3' for nucleotides). Makes the following
@@ -346,7 +389,9 @@ Details on the usage (output of ``ost compare-structures --help``):
thresholds [0.5, 1.0, 2.0, 4.0], "rmsd": RMSD given
these positions and transformation, "transform": the
used 4x4 transformation matrix that superposes model
- onto reference.
+ onto reference, "rigid_chain_mapping": equivalent of
+ "chain_mapping" which is used for rigid scores
+ (optimized for RMSD instead of QS-score/lDDT).
--patch-scores Local interface quality score used in CASP15. Scores
each model residue that is considered in the interface
(CB pos within 8A of any CB pos from another chain (CA
@@ -408,7 +453,7 @@ Details on the usage (output of ``ost compare-structures --help``):
sequences can be problematic as they may produce high
sequence identity alignments by pure chance.
-v VERBOSITY, --verbosity VERBOSITY
- Set verbosity level. Defaults to 3 (Script).
+ Set verbosity level. Defaults to 2 (Script).
--lddt-add-mdl-contacts
Only using contacts in lDDT thatare within a certain
distance threshold in the reference does not penalize
@@ -418,6 +463,9 @@ Details on the usage (output of ``ost compare-structures --help``):
necessarily in the reference. No contact will be added
if the respective atom pair is not resolved in the
reference.
+ --lddt-inclusion-radius LDDT_INCLUSION_RADIUS
+ Passed to lDDT scorer. Affects all lDDT scores but not
+ chain mapping.
.. _ost compare ligand structures:
@@ -438,14 +486,20 @@ Details on the usage (output of ``ost compare-ligand-structures --help``):
-r REFERENCE
[-rl [REFERENCE_LIGANDS ...]] [-o OUTPUT]
[-mf {pdb,cif,mmcif}]
- [-rf {pdb,cif,mmcif}] [-mb MODEL_BIOUNIT]
+ [-rf {pdb,cif,mmcif}] [-of {json,csv}]
+ [-csvm]
+ [--csv-extra-header CSV_EXTRA_HEADER]
+ [--csv-extra-data CSV_EXTRA_DATA]
+ [-mb MODEL_BIOUNIT]
[-rb REFERENCE_BIOUNIT] [-ft] [-rna]
[-sm] [-cd COVERAGE_DELTA] [-v VERBOSITY]
[--full-results] [--lddt-pli]
[--lddt-pli-radius LDDT_PLI_RADIUS]
- [--lddt-pli-amc] [--rmsd]
+ [--lddt-pli-add-mdl-contacts]
+ [--no-lddt-pli-add-mdl-contacts] [--rmsd]
[--radius RADIUS]
[--lddt-lp-radius LDDT_LP_RADIUS] [-fbs]
+ [-ms MAX_SYMMETRIES]
Evaluate model with non-polymer/small molecule ligands against reference.
@@ -480,8 +534,10 @@ Details on the usage (output of ``ost compare-ligand-structures --help``):
ligand connectivity to be loaded correctly. Ligands loaded from SDF files
are exempt from this restriction, meaning any arbitrary ligand can be assessed.
- Output is written in JSON format (default: out.json). In case of no additional
- options, this is a dictionary with three keys:
+ Output can be written in two format: JSON (default) or CSV, controlled by the
+ --output-format/-of argument.
+
+ Without additional options, the JSON ouput is a dictionary with four keys:
* "model_ligands": A list of ligands in the model. If ligands were provided
explicitly with --model-ligands, elements of the list will be the paths to
@@ -494,6 +550,7 @@ Details on the usage (output of ``ost compare-ligand-structures --help``):
* "status": SUCCESS if everything ran through. In case of failure, the only
content of the JSON output will be "status" set to FAILURE and an
additional key: "traceback".
+ * "ost_version": The OpenStructure version used for computation.
Each score is opt-in and the respective results are available in three keys:
@@ -519,6 +576,38 @@ Details on the usage (output of ``ost compare-ligand-structures --help``):
specific pair of ligands could be computed, "score" and "coverage" are set to
null and a key "reason" is added giving an educated guess why this happened.
+ CSV output is a table of comma-separated values, with one line for each
+ reference ligand (or one model ligand if the --by-model-ligand-output flag was
+ set).
+
+ The following column is always available:
+
+ * reference_ligand/model_ligand: If reference ligands were provided explicitly
+ with --reference-ligands, elements of the list will be the paths to the
+ ligand SDF file(s). Otherwise, they will be the chain name, residue number
+ and insertion code of the ligand, separated by a dot. If the
+ --by-model-ligand-output flag was set, this will be model ligand instead,
+ following the same rules.
+
+ If lDDT-PLI was enabled with --lddt-pli, the following columns are added:
+
+ * "lddt_pli", "lddt_pli_coverage" and "lddt_pli_(model|reference)_ligand"
+ are the lDDT-PLI score result, the corresponding coverage and assigned model
+ ligand (or reference ligand if the --by-model-ligand-output flag was set)
+ if an assignment was found, respectively, empty otherwise.
+ * "lddt_pli_unassigned" is empty if an assignment was found, otherwise it
+ lists the short reason this reference ligand was unassigned.
+
+ If BiSyRMSD was enabled with --rmsd, the following columns are added:
+
+ * "rmsd", "rmsd_coverage". "lddt_lp" "bb_rmsd" and
+ "rmsd_(model|reference)_ligand" are the BiSyRMSD, the corresponding
+ coverage, lDDT-LP, backbone RMSD and assigned model ligand (or reference
+ ligand if the --by-model-ligand-output flag was set) if an assignment
+ was found, respectively, empty otherwise.
+ * "rmsd_unassigned" is empty if an assignment was found, otherwise it
+ lists the short reason this reference ligand was unassigned.
+
options:
-h, --help show this help message and exit
-m MODEL, --mdl MODEL, --model MODEL
@@ -530,8 +619,8 @@ Details on the usage (output of ``ost compare-ligand-structures --help``):
-rl [REFERENCE_LIGANDS ...], --ref-ligands [REFERENCE_LIGANDS ...], --reference-ligands [REFERENCE_LIGANDS ...]
Path to reference ligand files.
-o OUTPUT, --out OUTPUT, --output OUTPUT
- Output file name. The output will be saved as a JSON
- file. default: out.json
+ Output file name. Default depends on format: out.json
+ or out.csv
-mf {pdb,cif,mmcif}, --mdl-format {pdb,cif,mmcif}, --model-format {pdb,cif,mmcif}
Format of model file. pdb reads pdb but also pdb.gz,
same applies to cif/mmcif. Inferred from filepath if
@@ -540,6 +629,19 @@ Details on the usage (output of ``ost compare-ligand-structures --help``):
Format of reference file. pdb reads pdb but also
pdb.gz, same applies to cif/mmcif. Inferred from
filepath if not given.
+ -of {json,csv}, --out-format {json,csv}, --output-format {json,csv}
+ Output format, JSON or CSV, in lowercase. default:
+ json
+ -csvm, --by-model-ligand, --by-model-ligand-output
+ For CSV output, this flag changes the output so that
+ each line reports one model ligand, instead of a
+ reference ligand. Has no effect with JSON output.
+ --csv-extra-header CSV_EXTRA_HEADER
+ Extra header prefix for CSV output. This allows adding
+ additional annotations (such as target ID, group, etc)
+ to the output
+ --csv-extra-data CSV_EXTRA_DATA
+ Additional data (columns) for CSV output.
-mb MODEL_BIOUNIT, --model-biounit MODEL_BIOUNIT
Only has an effect if model is in mmcif format. By
default, the asymmetric unit (AU) is used for scoring.
@@ -564,13 +666,16 @@ Details on the usage (output of ``ost compare-ligand-structures --help``):
-cd COVERAGE_DELTA, --coverage-delta COVERAGE_DELTA
Coverage delta for partial ligand assignment.
-v VERBOSITY, --verbosity VERBOSITY
- Set verbosity level. Defaults to 3 (INFO).
+ Set verbosity level. Defaults to 2 (Script).
--full-results Outputs scoring results for all model/reference ligand
pairs and store as key "full_results"
--lddt-pli Compute lDDT-PLI scores and store as key "lddt_pli".
--lddt-pli-radius LDDT_PLI_RADIUS
lDDT inclusion radius for lDDT-PLI.
- --lddt-pli-amc Add model contacts (amc) when computing lDDT-PLI.
+ --lddt-pli-add-mdl-contacts
+ Add model contacts when computing lDDT-PLI.
+ --no-lddt-pli-add-mdl-contacts
+ DO NOT add model contacts when computing lDDT-PLI.
--rmsd Compute RMSD scores and store as key "rmsd".
--radius RADIUS Inclusion radius to extract reference binding site
that is used for RMSD computation. Any residue with
@@ -582,3 +687,8 @@ Details on the usage (output of ``ost compare-ligand-structures --help``):
Enumerate all potential binding sites in the model
when searching rigid superposition for RMSD
computation
+ -ms MAX_SYMMETRIES, --max--symmetries MAX_SYMMETRIES
+ If more than that many isomorphisms exist for a
+ target-ligand pair, it will be ignored and reported as
+ unassigned.
+
diff --git a/modules/io/doc/io.rst b/modules/io/doc/io.rst
index 87bd4d3c808e7f6e6a2e990d1353c859ed3684a8..933f33cc3639d771808b3be71c64a26b5d23f7be 100644
--- a/modules/io/doc/io.rst
+++ b/modules/io/doc/io.rst
@@ -78,12 +78,15 @@ behaviour.
even if one is available in the IO Profile.
:param pdb_string: A PDB file as a string.
+ :type pdb_str: :class:`str`
:param profile: The IO Profile to read the entity with. For more information
on the IO Profiles available, see :doc:`profile`.
+ :type profile: :class:`ost.io.IOProfile`
:param process: If set to True, run the :class:`~ost.conop.Processor`
contained in the IO Profile.
+ :type process: :class:`bool`
:rtype: :class:`~ost.mol.EntityHandle`.
@@ -92,26 +95,42 @@ behaviour.
.. code-block:: python
- with open('protein.pdb') as pdb_fd:
- pdb_str = pdb.read()
+ with open('protein.pdb') as pdb_fh:
+ pdb_str = pdb_fh.read()
ent = io.PDBStrToEntity(pdb_str, ost.io.profiles['DEFAULT'], True)
-.. function:: LoadSDF(filename)
+.. function:: MMCifStrToEntity(mmcif_str, profile=IOProfile(), process=False)
- Load an SDF file and return an entity.
+ Load entity from a string in mmCIF format. By default the entity is loaded with
+ an empty IO Profile and is not processed with the :class:`~ost.conop.Processor`,
+ even if one is available in the IO Profile.
- :param filename: File to be loaded
- :type filename: :class:`str`
+ :param mmcif_str: mmCIF file as a string
+ :type mmcif_str: :class:`str`
+
+ :param profile: The IO Profile to read the entity with. For more information
+ on the IO Profiles available, see :doc:`profile`.
+ :type profile: :class:`ost.io.IOProfile`
+
+ :param process: If set to True, run the :class:`~ost.conop.Processor`
+ contained in the IO Profile.
+ :type process: :class:`bool`
+
+ :rtype: :class:`~ost.mol.EntityHandle`.
- :rtype: :class:`~ost.mol.EntityHandle`
+.. autofunction:: ost.io.LoadSDF
-.. function:: SDFStrToEntity(sdf_string)
+.. function:: SDFStrToEntity(sdf_string, profile=IOProfile())
Load entity from a string in SDF format.
:param pdb_string: A SDF file as a string.
+ :param profile: The IO Profile to read the entity with. For more information
+ on the IO Profiles available, see :doc:`profile`.
+ :type profile: :class:`ost.io.IOProfile`
+
:rtype: :class:`~ost.mol.EntityHandle`.
.. class:: ost.io.OMF
diff --git a/modules/io/doc/mmcif.rst b/modules/io/doc/mmcif.rst
index 1c7d17c98ed50a92261d15ba90175914fd39d1ed..23bd2a530692788147b3f248578710a3773d5e34 100644
--- a/modules/io/doc/mmcif.rst
+++ b/modules/io/doc/mmcif.rst
@@ -105,6 +105,9 @@ Notes:
It is a known limitation of the mmCIF format to allow ambiguous identifiers
for waters (and ligands to some extent) and so we have to require these
additional identifiers.
+* An mmCIF file can contain several models (``atom_site.pdbx_PDB_model_num``).
+ Only the first model occurring in the mmCIF file is read (regardless of the
+ actual model number). If extra models are ignored, a warning is logged.
Info Classes
diff --git a/modules/io/doc/structure_formats.rst b/modules/io/doc/structure_formats.rst
index 2fe0eaf676f9bc5f146225563a398011e0614a71..6903ed50f015523be82296cd0b61302bada00a86 100644
--- a/modules/io/doc/structure_formats.rst
+++ b/modules/io/doc/structure_formats.rst
@@ -51,12 +51,13 @@ radii.
SDF - Structured Data File
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-Chemical table (Ctab) file format (V2000; read-only V3000 experimental support),
-aka MDL Molfile.
-The SDF format does not support residues, chains or atom names natively.
+Chemical table (Ctab) file format, aka SDF, aka MDL Molfile, originally named
+after Molecular Design Limited, Inc. Follows the "BIOVIA Databases 2020" format
+specification (V2000; read-only V3000 experimental support).
The SDF importer supports loading gzipped files, which are auto-detected by the
.gz file extension.
+The SDF format does not support residues, chains or atom names natively.
The reader assigns 1-based atom indices as atom names.
SDF files containing several molecules are loaded into distinct chains,
named after the molecule name in the MOLfile header with a numerical prefix.
@@ -66,5 +67,5 @@ Chains are written as separate molecules. If a chain contains more than one
residue, they will be merged into a single molecule.
*Recognized File Extensions*
- .sdf, .sdf.gz
+ .sdf, .sdf.gz, .mol, .mol.gz
diff --git a/modules/io/pymod/CMakeLists.txt b/modules/io/pymod/CMakeLists.txt
index ead77d54b9d4cc0285845a166a52287774d73f67..02b392754d6b44ec9948b621ef86903ebf3ce060 100644
--- a/modules/io/pymod/CMakeLists.txt
+++ b/modules/io/pymod/CMakeLists.txt
@@ -4,6 +4,7 @@ set(OST_IO_PYMOD_SOURCES
export_mmcif_io.cc
export_omf_io.cc
export_map_io.cc
+ export_sdf_io.cc
)
set(OST_IO_PYMOD_MODULES
diff --git a/modules/io/pymod/__init__.py b/modules/io/pymod/__init__.py
index 66a518def1a71f81b638090da38f3a743584e2ae..dca962c0bcd74611721a5e439e28d6eb1511f2d0 100644
--- a/modules/io/pymod/__init__.py
+++ b/modules/io/pymod/__init__.py
@@ -348,9 +348,10 @@ def LoadCHARMMTraj(crd, dcd_file=None, profile='CHARMM',
def LoadMMCIF(filename, fault_tolerant=None, calpha_only=None,
profile='DEFAULT', remote=False, seqres=False, info=False):
"""
- Load a mmCIF file and return one or more entities. Several options allow to
- customize the exact behaviour of the mmCIF import. For more information on
- these options, see :doc:`profile`.
+ Load an mmCIF file and return the first model as an entity.
+
+ Several options allow to customize the exact behaviour of the mmCIF import.
+ For more information on these options, see :doc:`profile`.
Residues are flagged as ligand if they are not waters nor covered by an
``entity_poly`` record (ie. they are non-polymer entities in
@@ -555,3 +556,46 @@ def _PDBize(biounit, asu, seqres=None, min_polymer_size=None,
return pdb_bu
MMCifInfoBioUnit.PDBize = _PDBize
+
+
+def LoadSDF(filename, fault_tolerant=None, profile='DEFAULT'):
+ """
+ Load SDF file from disk and return an entity.
+
+ :param filename: File to be loaded
+ :type filename: :class:`str`
+
+ :param fault_tolerant: Enable/disable fault-tolerant import. If set, overrides
+ the value of :attr:`IOProfile.fault_tolerant`.
+ :type fault_tolerant: :class:`bool`
+
+ :param profile: Aggregation of flags and algorithms to control import and
+ processing of molecular structures. Can either be a
+ :class:`str` specifying one of the default profiles
+ ['DEFAULT', 'SLOPPY', 'CHARMM', 'STRICT'] or an actual object
+ of type :class:`ost.io.IOProfile`.
+ See :doc:`profile` for more info.
+ :type profile: :class:`str`/:class:`ost.io.IOProfile`
+
+ :raises: :exc:`~ost.io.IOException` if the import fails due to an erroneous or
+ inexistent file
+ """
+ def _override(val1, val2):
+ if val2 != None:
+ return val2
+ else:
+ return val1
+
+ if isinstance(profile, str):
+ prof = profiles.Get(profile)
+ elif isinstance(profile, IOProfile):
+ prof = profile.Copy()
+ else:
+ raise TypeError('profile must be of type string or IOProfile, ' + \
+ 'instead of %s' % type(profile))
+ prof.fault_tolerant = _override(prof.fault_tolerant, fault_tolerant)
+
+ reader = SDFReader(filename, prof)
+ ent = mol.CreateEntity()
+ reader.Import(ent)
+ return ent
diff --git a/modules/io/pymod/export_sdf_io.cc b/modules/io/pymod/export_sdf_io.cc
new file mode 100644
index 0000000000000000000000000000000000000000..08de790e2850229a568920c790767992491c075c
--- /dev/null
+++ b/modules/io/pymod/export_sdf_io.cc
@@ -0,0 +1,64 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2020 by the OpenStructure authors
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation; either version 3.0 of the License, or (at your option)
+// any later version.
+// This library 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 Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this library; if not, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+//------------------------------------------------------------------------------
+#include <boost/python.hpp>
+#include <boost/shared_ptr.hpp>
+#include <boost/python/suite/indexing/vector_indexing_suite.hpp>
+#include <boost/python/suite/indexing/map_indexing_suite.hpp>
+using namespace boost::python;
+
+#include <ost/export_helper/pair_to_tuple_conv.hh>
+#include <ost/io/mol/entity_io_sdf_handler.hh>
+#include <ost/io/mol/io_profile.hh>
+#include <ost/io/mol/sdf_reader.hh>
+#include <ost/io/mol/sdf_writer.hh>
+#include <ost/io/sdf_str.hh>
+using namespace ost;
+using namespace ost::io;
+using namespace ost::mol;
+
+String (*sdf_str_a)(const mol::EntityHandle&)=&EntityToSDFString;
+String (*sdf_str_b)(const mol::EntityView&)=&EntityToSDFString;
+
+void (*save_sdf_handle)(const mol::EntityHandle& entity, const String& filename)=&SaveSDF;
+void (*save_sdf_view)(const mol::EntityView& entity, const String& filename)=&SaveSDF;
+
+void (SDFWriter::*write_handle)(const mol::EntityHandle&)=&SDFWriter::Write;
+void (SDFWriter::*write_view)(const mol::EntityView&)=&SDFWriter::Write;
+
+void export_sdf_io()
+{
+ class_<SDFReader, boost::noncopyable>("SDFReader", init<String, const IOProfile&>())
+ .def("Import", &SDFReader::Import)
+ ;
+
+ class_<SDFWriter, boost::noncopyable>("SDFWriter", init<String>())
+ .def("Write", write_handle)
+ .def("Write", write_view)
+ ;
+
+ // def("LoadSDF", &LoadSDF);
+ def("SaveSDF", save_sdf_view);
+ def("SaveSDF", save_sdf_handle);
+
+ def("EntityToSDFStr", sdf_str_a);
+ def("EntityToSDFStr", sdf_str_b);
+
+ def("SDFStrToEntity", &SDFStringToEntity, (arg("SDF_string"),
+ arg("profile")=IOProfile()));
+}
diff --git a/modules/io/pymod/wrap_io.cc b/modules/io/pymod/wrap_io.cc
index bae73b1e0fd9168c7673fc39da52ce678d6878d9..b7cade00c33e39eefd2ecc982b7e2c4e2d5a2a6f 100644
--- a/modules/io/pymod/wrap_io.cc
+++ b/modules/io/pymod/wrap_io.cc
@@ -30,9 +30,7 @@ using namespace boost::python;
#include <ost/io/mol/entity_io_crd_handler.hh>
#include <ost/io/mol/entity_io_pqr_handler.hh>
#include <ost/io/mol/entity_io_mae_handler.hh>
-#include <ost/io/mol/entity_io_sdf_handler.hh>
#include <ost/io/mol/pdb_reader.hh>
-#include <ost/io/mol/sdf_str.hh>
#include <ost/io/mol/dcd_io.hh>
#include <ost/io/stereochemical_params_reader.hh>
using namespace ost;
@@ -66,18 +64,13 @@ BOOST_PYTHON_FUNCTION_OVERLOADS(save_entity_view_ov,
ost::mol::alg::StereoChemicalProps (*read_props_a)(String filename, bool check) = &ReadStereoChemicalPropsFile;
ost::mol::alg::StereoChemicalProps (*read_props_b)(bool check) = &ReadStereoChemicalPropsFile;
-String (*sdf_str_a)(const mol::EntityHandle&)=&EntityToSDFString;
-String (*sdf_str_b)(const mol::EntityView&)=&EntityToSDFString;
-
-void (*save_sdf_handle)(const mol::EntityHandle& entity, const String& filename)=&SaveSDF;
-void (*save_sdf_view)(const mol::EntityView& entity, const String& filename)=&SaveSDF;
-
}
void export_pdb_io();
void export_mmcif_io();
void export_omf_io();
void export_map_io();
+void export_sdf_io();
BOOST_PYTHON_MODULE(_ost_io)
{
class_<IOManager, boost::noncopyable>("IOManager", no_init)
@@ -122,14 +115,6 @@ BOOST_PYTHON_MODULE(_ost_io)
(arg("seq_list"), arg("filename"), arg("format")="auto"));
def("SaveSequence", &SaveSequence,
(arg("sequence"), arg("filename"), arg("format")="auto"));
- def("LoadSDF", &LoadSDF);
- def("SaveSDF", save_sdf_view);
- def("SaveSDF", save_sdf_handle);
-
- def("EntityToSDFStr", sdf_str_a);
- def("EntityToSDFStr", sdf_str_b);
-
- def("SDFStrToEntity", &SDFStringToEntity);
def("LoadCRD", &LoadCRD);
def("LoadCHARMMTraj_", &LoadCHARMMTraj, (arg("ent"), arg("trj_filename"),
@@ -149,6 +134,7 @@ BOOST_PYTHON_MODULE(_ost_io)
export_mmcif_io();
export_omf_io();
export_map_io();
+ export_sdf_io();
def("SaveCHARMMTraj", &SaveCHARMMTraj,
(arg("traj"), arg("pdb_filename"), arg("dcd_filename"), arg("stride")=1,
arg("profile")=IOProfile()));
diff --git a/modules/io/src/mol/chemdict_parser.cc b/modules/io/src/mol/chemdict_parser.cc
index 2213225d696f832b460799d00cda331f8bd24354..f1b2610ad549f78090b6ff2fb14ce47af92d4413 100644
--- a/modules/io/src/mol/chemdict_parser.cc
+++ b/modules/io/src/mol/chemdict_parser.cc
@@ -7,6 +7,8 @@ using namespace ost::conop;
bool ChemdictParser::OnBeginData(const StringRef& data_name)
{
+ // Create empty compound skeleton
+ valid_compound_ = false;
compound_.reset(new Compound(data_name.str()));
compound_->SetDialect(dialect_);
if (last_!=data_name[0]) {
@@ -14,6 +16,7 @@ bool ChemdictParser::OnBeginData(const StringRef& data_name)
std::cout << last_ << std::flush;
}
atom_map_.clear();
+ valid_atom_ = false;
return true;
}
@@ -94,6 +97,14 @@ void ChemdictParser::OnDataRow(const StarLoopDesc& header,
void ChemdictParser::OnDataItem(const StarDataItem& item)
{
if (item.GetCategory()==StringRef("chem_comp", 9)) {
+ if (item.GetName()==StringRef("id", 2)) {
+ if (compound_->GetID() != item.GetValue().str()) {
+ LOG_INFO("_chem_comp.id '" << item.GetValue() << "' doesn't match"
+ << "ID from data block '" << compound_->GetID() << "'");
+ compound_->SetID(item.GetValue().str());
+ }
+ valid_compound_ = true;
+ }
if (item.GetName()==StringRef("type", 4)) {
// convert type to uppercase
String type=item.GetValue().str();
@@ -144,21 +155,28 @@ void ChemdictParser::OnDataItem(const StarDataItem& item)
compound_->SetObsolete(true);
}
} else if (item.GetName()==StringRef("pdbx_replaced_by", 16)) {
- String replaced_by = item.GetValue().str();
- if (replaced_by != "?") {
- compound_->SetReplacedBy(replaced_by);
+ StringRef replaced_by = item.GetValue();
+ if (! IsUndefined(replaced_by)) {
+ compound_->SetReplacedBy(replaced_by.str());
}
} else if (item.GetName()==StringRef("one_letter_code", 15)) {
if (item.GetValue().length()==1) {
compound_->SetOneLetterCode(item.GetValue()[0]);
}
- } else if (item.GetName()==StringRef("pdbx_initial_date", 17)) {
- compound_->SetCreationDate(Date::FromString(item.GetValue()));
+ } else if (item.GetName()==StringRef("pdbx_initial_date", 17)) {
+ StringRef pdbx_initial_date = item.GetValue();
+ if (! IsUndefined(pdbx_initial_date)) {
+ compound_->SetCreationDate(Date::FromString(pdbx_initial_date));
+ }
} else if (item.GetName()==StringRef("pdbx_modified_date", 18)) {
- compound_->SetModificationDate(Date::FromString(item.GetValue()));
+ StringRef pdbx_modified_date = item.GetValue();
+ if (! IsUndefined(pdbx_modified_date)) {
+ compound_->SetModificationDate(Date::FromString(pdbx_modified_date));
+ }
}
} else if (item.GetName()==StringRef("atom_id", 7)) {
atom_.name=item.GetValue().str();
+ valid_atom_ = true;
} else if (item.GetName()==StringRef("alt_atom_id", 11)) {
if (compound_->GetID()=="ILE" && item.GetValue()==StringRef("CD1", 3)) {
atom_.alt_name="CD";
@@ -177,14 +195,27 @@ void ChemdictParser::OnEndData()
{
if (compound_)
{
- if (compound_->GetID() != "UNL" &&
+ if (! valid_compound_)
+ {
+ LOG_WARNING("Skipping compound without _chem_comp.id: " << compound_->GetID());
+ }
+ else if (compound_->GetID() != "UNL" &&
! (ignore_reserved_ && IsNameReserved(compound_->GetID())) &&
! (ignore_obsolete_ && compound_->GetObsolete()))
{
if (compound_->GetAtomSpecs().empty()) {
- compound_->AddAtom(atom_);
+ // This happens if we had a single atom
+ if (valid_atom_)
+ {
+ compound_->AddAtom(atom_);
+ }
+ else
+ {
+ LOG_WARNING("Adding compound with no atoms: " << compound_->GetID());
+ }
}
lib_->AddCompound(compound_);
+ ++imported_count_;
}
}
}
diff --git a/modules/io/src/mol/chemdict_parser.hh b/modules/io/src/mol/chemdict_parser.hh
index 7a30cc8fcb434fe825ed887b3159ac5f7ed16180..fe4981c116945e2a1031def2d5c97dc9c435e1d0 100644
--- a/modules/io/src/mol/chemdict_parser.hh
+++ b/modules/io/src/mol/chemdict_parser.hh
@@ -44,7 +44,8 @@ public:
bool ignore_reserved=false, bool ignore_obsolete=false):
StarParser(stream), compound_(new conop::Compound("UNK")),
last_(0), loop_type_(DONT_KNOW), dialect_(dialect),
- ignore_reserved_(ignore_reserved), ignore_obsolete_(ignore_obsolete)
+ ignore_reserved_(ignore_reserved), ignore_obsolete_(ignore_obsolete),
+ imported_count_(0)
{
this->InitTypeMap();
this->InitPDBXTypeMap();
@@ -65,12 +66,19 @@ public:
{
lib_=lib;
}
+
+ /// \brief Get the number of compounds that were successfully parsed
+ int GetImportedCount()
+ {
+ return imported_count_;
+ }
private:
void InitTypeMap();
void InitPDBXTypeMap();
bool IsNameReserved(const String& data_name);
conop::CompoundLibPtr lib_;
conop::CompoundPtr compound_;
+ bool valid_compound_;
typedef enum {
ATOM_NAME=0,
ALT_ATOM_NAME=1,
@@ -93,9 +101,11 @@ private:
std::map<String, int> atom_map_;
LoopType loop_type_;
conop::AtomSpec atom_;
+ bool valid_atom_;
conop::Compound::Dialect dialect_;
bool ignore_reserved_;
bool ignore_obsolete_;
+ int imported_count_;
};
diff --git a/modules/io/src/mol/entity_io_sdf_handler.cc b/modules/io/src/mol/entity_io_sdf_handler.cc
index 129ce2c86be7a7ec56e1cbee52488a231411067a..fed25f86b2b6c745341a51a4fe6fa770feff500f 100644
--- a/modules/io/src/mol/entity_io_sdf_handler.cc
+++ b/modules/io/src/mol/entity_io_sdf_handler.cc
@@ -21,6 +21,7 @@
*/
#include <ost/log.hh>
+#include <ost/profile.hh>
#include <ost/io/mol/sdf_writer.hh>
#include <ost/io/mol/sdf_reader.hh>
@@ -37,14 +38,14 @@ bool EntityIOSDFHandler::RequiresProcessor() const
void EntityIOSDFHandler::Import(mol::EntityHandle& ent,
std::istream& instream)
{
- SDFReader reader(instream);
+ SDFReader reader(instream, IOProfileRegistry::Instance().GetDefault());
reader.Import(ent);
}
void EntityIOSDFHandler::Import(mol::EntityHandle& ent,
const boost::filesystem::path& loc)
{
- SDFReader reader(loc);
+ SDFReader reader(loc, IOProfileRegistry::Instance().GetDefault());
reader.Import(ent);
}
@@ -69,7 +70,8 @@ bool sdf_handler_is_responsible_for(const boost::filesystem::path& loc,
if(type=="auto") {
String match_suf_string=loc.string();
std::transform(match_suf_string.begin(),match_suf_string.end(),match_suf_string.begin(),tolower);
- if(detail::FilenameEndsWith(match_suf_string,".sdf") || detail::FilenameEndsWith(match_suf_string,".sdf.gz")) {
+ if(detail::FilenameEndsWith(match_suf_string,".sdf") || detail::FilenameEndsWith(match_suf_string,".sdf.gz") ||
+ detail::FilenameEndsWith(match_suf_string,".mol") || detail::FilenameEndsWith(match_suf_string,".mol.gz")) {
return true;
}
diff --git a/modules/io/src/mol/mmcif_info.cc b/modules/io/src/mol/mmcif_info.cc
index 77ae604029e8a0b71a00a3c0b48fec6972c5bf61..da292223460a6717e772a3abe60f374e0025bdae 100644
--- a/modules/io/src/mol/mmcif_info.cc
+++ b/modules/io/src/mol/mmcif_info.cc
@@ -196,11 +196,20 @@ void MMCifInfoBioUnit::Merge(MMCifInfoBioUnit& from)
MMCifInfoStructRefSeqPtr
MMCifInfoStructRef::AddAlignedSeq(const String& aid, const String& chain_name,
int seq_begin, int seq_end, int db_begin,
- int db_end)
+ int db_end, bool fault_tolerant)
{
std::map<String, MMCifInfoStructRefSeqPtr>::const_iterator i=seqs_.find(aid);
if (i!=seqs_.end()) {
- throw IOException("duplicate align_id for struct_ref '"+id_+"'");
+ std::stringstream msg;
+ msg << "Duplicate align_id for struct_ref '" << id_ << "'";
+ if (fault_tolerant) {
+ msg << ". Record will be overwritten.";
+ LOG_WARNING(msg.str());
+ }
+ else {
+ throw IOException(msg.str());
+ }
+
}
MMCifInfoStructRefSeqPtr p(new MMCifInfoStructRefSeq(aid, chain_name,
seq_begin, seq_end,
diff --git a/modules/io/src/mol/mmcif_info.hh b/modules/io/src/mol/mmcif_info.hh
index 73d3c25d2b69feafb098ed4f0d1dbd00add0657a..73946c1393d6e285b0b66637ad89ad977a462a26 100644
--- a/modules/io/src/mol/mmcif_info.hh
+++ b/modules/io/src/mol/mmcif_info.hh
@@ -858,7 +858,8 @@ public:
const String& GetDBAccess() const { return db_access_; }
MMCifInfoStructRefSeqPtr AddAlignedSeq(const String& align_id,
const String& chain_name, int seq_begin,
- int seq_end, int db_begin, int db_end);
+ int seq_end, int db_begin, int db_end,
+ bool fault_tolerant);
MMCifInfoStructRefSeqPtr GetAlignedSeq(const String& align_id) const;
MMCifInfoStructRefSeqs GetAlignedSeqs() const
{
diff --git a/modules/io/src/mol/mmcif_reader.cc b/modules/io/src/mol/mmcif_reader.cc
index fd440812ea485a3174501f438ff9e82ac4d4ed62..03abbad79f9b940e99c556c7ab0bb654d20c648c 100644
--- a/modules/io/src/mol/mmcif_reader.cc
+++ b/modules/io/src/mol/mmcif_reader.cc
@@ -31,11 +31,6 @@
namespace ost { namespace io {
-bool is_undef(StringRef value)
-{
- return value.size()==1 && (value[0]=='?' || value[0]=='.');
-}
-
MMCifReader::MMCifReader(std::istream& stream, mol::EntityHandle& ent_handle,
const IOProfile& profile):
StarParser(stream, true), profile_(profile), ent_handle_(ent_handle)
@@ -97,8 +92,20 @@ void MMCifReader::SetRestrictChains(const String& restrict_chains)
bool MMCifReader::OnBeginData(const StringRef& data_name)
{
- LOG_DEBUG("MCIFFReader: " << profile_);
+ LOG_DEBUG("MMCifReader: " << profile_);
Profile profile_import("MMCifReader::OnBeginData");
+ if (chain_count_ > 0) {
+ std::stringstream ss;
+ ss << "Can only read one data block. Found second one 'data_" << data_name << "'.";
+ if (profile_.fault_tolerant) {
+ LOG_WARNING(ss.str());
+ } else {
+ throw IOException(this->FormatDiagnostic(STAR_DIAG_ERROR,
+ ss.str(),
+ this->GetCurrentLinenum()));
+ }
+
+ }
// IDs in mmCIF files can be any string, so no restrictions here
@@ -478,15 +485,21 @@ void MMCifReader::ParseAndAddAtom(const std::vector<StringRef>& columns)
mol::ResNum res_num(0);
bool valid_res_num = false;
if (indices_[PDBX_PDB_MODEL_NUM] != -1) {
+ int model_id = TryGetInt(columns[indices_[PDBX_PDB_MODEL_NUM]],
+ "atom_site.pdbx_PDB_model_num");
if (has_model_) {
- if (curr_model_ != TryGetInt(columns[indices_[PDBX_PDB_MODEL_NUM]],
- "atom_site.pdbx_PDB_model_num")) {
+ if (curr_model_ != model_id) {
+ if (warned_ignored_model_.find(model_id) == warned_ignored_model_.end()) {
+ LOG_WARNING("Ignorning new model " << model_id <<
+ ". Only model " << curr_model_ << " was read.");
+ warned_ignored_model_.insert(model_id);
+ }
return;
}
} else {
+ LOG_INFO("Reading model " << model_id << ".");
has_model_ = true;
- curr_model_ = TryGetInt(columns[indices_[PDBX_PDB_MODEL_NUM]],
- "atom_site.pdbx_PDB_model_num");
+ curr_model_ = model_id;
}
}
@@ -520,10 +533,10 @@ void MMCifReader::ParseAndAddAtom(const std::vector<StringRef>& columns)
occ = this->GetRealOrDefault(columns[indices_[OCCUPANCY]],
"atom_site.occupancy",
1.0,
- is_undef);
+ IsUndefined);
}
if (indices_[B_ISO_OR_EQUIV] != -1) {
- if (!is_undef(columns[indices_[B_ISO_OR_EQUIV]])) {
+ if (!IsUndefined(columns[indices_[B_ISO_OR_EQUIV]])) {
temp = this->TryGetReal(columns[indices_[B_ISO_OR_EQUIV]],
"atom_site.B_iso_or_equiv");
}
@@ -677,6 +690,7 @@ void MMCifReader::ParseAndAddAtom(const std::vector<StringRef>& columns)
"name. Ignoring atoms for everything but the first");
} else {
LOG_WARNING("Residue with number " << res_num
+ << " in chain " << curr_chain_.GetName()
<< " contains a microheterogeneity. Everything but atoms "
"for the residue '" << curr_residue_.GetName()
<< "' will be ignored");
@@ -871,13 +885,13 @@ void MMCifReader::ParseCitation(const std::vector<StringRef>& columns)
}
}
if (indices_[PDBX_DATABASE_ID_PUBMED] != -1) {
- if (!is_undef(columns[indices_[PDBX_DATABASE_ID_PUBMED]])) {
+ if (!IsUndefined(columns[indices_[PDBX_DATABASE_ID_PUBMED]])) {
cit.SetPubMed(this->TryGetInt(columns[indices_[PDBX_DATABASE_ID_PUBMED]],
"citation.pdbx_database_id_PubMed"));
}
}
if (indices_[YEAR] != -1) {
- if (!is_undef(columns[indices_[YEAR]])) {
+ if (!IsUndefined(columns[indices_[YEAR]])) {
cit.SetYear(this->TryGetInt(columns[indices_[YEAR]], "citation.year"));
}
}
@@ -1176,7 +1190,7 @@ void MMCifReader::ParseStruct(const std::vector<StringRef>& columns)
}
if (indices_[PDBX_FORMULA_WEIGHT] != -1) {
- if (!is_undef(columns[indices_[PDBX_FORMULA_WEIGHT]])) {
+ if (!IsUndefined(columns[indices_[PDBX_FORMULA_WEIGHT]])) {
details.SetMass(this->TryGetReal(columns[indices_[PDBX_FORMULA_WEIGHT]],
"struct.pdbx_formula_weight"));
}
@@ -1680,7 +1694,7 @@ void MMCifReader::ParseStructRefSeq(const std::vector<StringRef>& columns)
e=struct_refs_.end(); i!=e; ++i) {
if ((*i)->GetID()==sr_id) {
(*i)->AddAlignedSeq(aln_id, chain_name, entbeg.second, entend.second,
- dbbeg.second, dbend.second);
+ dbbeg.second, dbend.second, profile_.fault_tolerant);
found=true;
break;
}
@@ -1711,7 +1725,7 @@ void MMCifReader::ParseStructRefSeqDif(const std::vector<StringRef>& columns)
std::pair<bool,int> seq_rnum;
if (indices_[SRSD_SEQ_RNUM] != -1) {
- if (!is_undef(columns[indices_[SRSD_SEQ_RNUM]])) {
+ if (!IsUndefined(columns[indices_[SRSD_SEQ_RNUM]])) {
seq_rnum=this->TryGetInt(columns[indices_[SRSD_SEQ_RNUM]],
"_struct_ref_seq_dif.seq_num",
profile_.fault_tolerant);
diff --git a/modules/io/src/mol/mmcif_reader.hh b/modules/io/src/mol/mmcif_reader.hh
index 23cfc2f4bc50c66bc4312680e1fbddeb6d77e879..8521c048d71c94c00d7cd8f63978130073fcd6f0 100644
--- a/modules/io/src/mol/mmcif_reader.hh
+++ b/modules/io/src/mol/mmcif_reader.hh
@@ -691,6 +691,7 @@ private:
String subst_res_id_; ///< work around for missing label_seq_id's
bool has_model_; ///< keep track of models through different atom_sites
int curr_model_; ///< if we have pdbx_PDB_model_num, store no.
+ std::set<int> warned_ignored_model_; // keep track of ignored model warnings
std::vector<std::pair<mol::ChainHandle, String> > chain_id_pairs_;
///< chain and label_entity_id
MMCifEntityDescMap entity_desc_map_; ///< stores entity items
diff --git a/modules/io/src/mol/mmcif_writer.cc b/modules/io/src/mol/mmcif_writer.cc
index 6ea364aa868975fbb5cee1a8956809de576726d6..3e086fa8aa986928021857cdd3488b434f57bc78 100644
--- a/modules/io/src/mol/mmcif_writer.cc
+++ b/modules/io/src/mol/mmcif_writer.cc
@@ -811,6 +811,19 @@ namespace {
String comp_id = res.GetName();
auto at_list = res.GetAtomList();
+
+ // sanity check that we have no duplicates
+ std::set<String> anames;
+ for(auto a: at_list) {
+ if(anames.find(a.GetName()) != anames.end()) {
+ std::stringstream ss;
+ ss << "Duplicate atom \"" << a.GetName() << "\" in residue ";
+ ss << a.GetResidue() << std::endl;
+ throw ost::io::IOException(ss.str());
+ }
+ anames.insert(a.GetName());
+ }
+
String auth_asym_id = res.GetChain().GetName();
if(res.HasProp("pdb_auth_chain_name")) {
auth_asym_id = res.GetStringProp("pdb_auth_chain_name");
@@ -1328,9 +1341,35 @@ namespace {
auto chain_list = ent.GetChainList();
for(auto chain: chain_list) {
String cname = chain.GetName();
+
+ auto res_list = chain.GetResidueList();
+ // sanity check that residue numbers are strictly increasing and have no
+ // insertion codes which is a requirement for mmCIF conform chains
+ // (not only for polymers)
+ for(size_t i = 0; i < res_list.size(); ++i) {
+ if(res_list[i].GetNumber().GetInsCode() != '\0') {
+ std::stringstream ss;
+ ss << "Residue number insertion codes must be empty if ";
+ ss << "mmcif_conform is enabled. Got \"";
+ ss << res_list[i].GetNumber().GetInsCode() << "\" in residue ";
+ ss << res_list[i];
+ throw ost::io::IOException(ss.str());
+ }
+ if(i>0) {
+ if(res_list[i-1].GetNumber().GetNum() >=
+ res_list[i].GetNumber().GetNum()) {
+ std::stringstream ss;
+ ss << "Residue numbers must be strictly increasing in consecutive ";
+ ss << "residues if mmcif_conform is enabled. ";
+ ss << "Got " << res_list[i-1] << " followed by " << res_list[i];
+ ss << ".";
+ throw ost::io::IOException(ss.str());
+ }
+ }
+ }
+
if(preassigned_polymer_chains.find(cname) !=
preassigned_polymer_chains.end()) {
- auto res_list = chain.GetResidueList();
int entity_id = preassigned_polymer_chains[cname];
AddAsymResnum(cname, res_list, entity_info[entity_id], true);
Feed_atom_site(atom_site, cname, entity_id+1, entity_info[entity_id], res_list);
@@ -1338,7 +1377,6 @@ namespace {
entity_info[entity_id], res_list);
} else {
// do automated matching
- auto res_list = chain.GetResidueList();
int entity_id = SetupEntity(cname,
chain.GetType(),
res_list,
diff --git a/modules/io/src/mol/sdf_reader.cc b/modules/io/src/mol/sdf_reader.cc
index 76c7091a297ab99fc422f7bc075e89e98cecb5fb..edafc0ff1cb401af7aa63878a8f5074799aa3bf0 100644
--- a/modules/io/src/mol/sdf_reader.cc
+++ b/modules/io/src/mol/sdf_reader.cc
@@ -37,20 +37,20 @@ namespace ost { namespace io {
using boost::format;
-SDFReader::SDFReader(const String& filename)
- : infile_(filename), instream_(infile_)
+SDFReader::SDFReader(const String& filename, const IOProfile& profile)
+ : infile_(filename), instream_(infile_), profile_(profile)
{
this->ClearState(boost::filesystem::path(filename));
}
-SDFReader::SDFReader(const boost::filesystem::path& loc)
- : infile_(loc), instream_(infile_)
+SDFReader::SDFReader(const boost::filesystem::path& loc, const IOProfile& profile)
+ : infile_(loc), instream_(infile_), profile_(profile)
{
this->ClearState(loc);
}
-SDFReader::SDFReader(std::istream& instream)
- : infile_(), instream_(instream)
+SDFReader::SDFReader(std::istream& instream, const IOProfile& profile)
+ : infile_(), instream_(instream), profile_(profile)
{
this->ClearState(boost::filesystem::path(""));
}
@@ -113,9 +113,11 @@ void SDFReader::Import(mol::EntityHandle& ent)
ProcessV3000Line(line, ent, editor);
}
}
-
- LOG_INFO("imported " << chain_count_ << " chains, " << residue_count_
- << " residues, " << atom_count_ << " atoms");
+ LOG_INFO("imported "
+ << ent.GetChainCount() << " chains, "
+ << ent.GetResidueCount() << " residues, "
+ << ent.GetAtomCount() << " atoms and "
+ << ent.GetBondCount() << " bonds");
}
void SDFReader::ClearState(const boost::filesystem::path& loc)
@@ -204,7 +206,7 @@ void SDFReader::ParseHeader(const String& line, int line_num,
version_=version_str;
}
else {
- String msg="Unsupported SDF version: %s.";
+ String msg="Invalid SDF file or unsupported SDF version: %s.";
throw IOException(str(format(msg) % version_str));
}
// Counts will be overridden in V3000
@@ -353,15 +355,40 @@ void SDFReader::AddBond(const bond_data& bond_tuple, int line_num, mol::EntityHa
try {
type=boost::lexical_cast<int>(boost::trim_copy(s_type));
- if (type<1 || type>8) {
- String msg="Bad bond line %d: Bond type number"
- " '%s' not within accepted range (1-8).";
- throw IOException(str(format(msg) % line_num % s_type));
+ // From SDF spec:
+ // bond type 1 = Single, 2 = Double, [Q] Values 4 through 8 are
+ // 3 = Triple, 4 = Aromatic, for SSS queries oniy.
+ // 5 = Single or Double,
+ // 6 = Single or Aromatic,
+ // 7 = Double or Aromatic,
+ // 8 = Any
+ if (type < 1 || type > 8) {
+ std::stringstream ss;
+ ss << "Bad bond line " << line_num << ": Bond type number "
+ << std::to_string(type) << " not within accepted range (1-8).";
+ if (profile_.fault_tolerant) {
+ LOG_ERROR(ss.str());
+ } else {
+ throw IOException(ss.str());
+ }
+ } else if (type > 3) {
+ std::stringstream ss;
+ ss << "Bad bond line " << line_num << ": Bond type number "
+ << std::to_string(type) << ": values 4-8 are reserved for queries, "
+ << "should not appear in an SDF file.";
+ LOG_WARNING(ss.str());
}
} catch(boost::bad_lexical_cast&) {
- String msg="Bad bond line %d: Can't convert bond type number"
- " '%s' to integral constant.";
- throw IOException(str(format(msg) % line_num % s_type));
+ std::stringstream ss;
+ ss << "Bad bond line " << line_num << ": Can't convert bond type number '"
+ << s_type << "' to integral constant.";
+ if (profile_.fault_tolerant) {
+ ss << " Assuming single bond in fault tolerant mode.";
+ LOG_ERROR(ss.str());
+ type = 1;
+ } else {
+ throw IOException(ss.str());
+ }
}
mol::AtomHandle first,second;
diff --git a/modules/io/src/mol/sdf_reader.hh b/modules/io/src/mol/sdf_reader.hh
index 59e733a937a3ed0dcb7171deec9ff9a32dd808ae..d524a8bf1561323dcbb2b91c10c70097673b352b 100644
--- a/modules/io/src/mol/sdf_reader.hh
+++ b/modules/io/src/mol/sdf_reader.hh
@@ -28,6 +28,7 @@
#include <ost/mol/chain_handle.hh>
#include <ost/mol/residue_handle.hh>
#include <ost/io/module_config.hh>
+#include <ost/io/mol/io_profile.hh>
namespace ost { namespace io {
@@ -36,11 +37,9 @@ namespace ost { namespace io {
class DLLEXPORT_OST_IO SDFReader {
public:
- SDFReader(const String& filename);
- SDFReader(const boost::filesystem::path& loc);
- SDFReader(std::istream& instream);
-
- bool HasNext();
+ SDFReader(const String& filename, const IOProfile& profile);
+ SDFReader(const boost::filesystem::path& loc, const IOProfile& profile);
+ SDFReader(std::istream& instream, const IOProfile& profile);
void Import(mol::EntityHandle& ent);
@@ -97,6 +96,7 @@ private:
boost::filesystem::ifstream infile_;
std::istream& instream_;
boost::iostreams::filtering_stream<boost::iostreams::input> in_;
+ IOProfile profile_;
String version_;
bool v3000_atom_block_;
bool v3000_bond_block_;
diff --git a/modules/io/src/mol/sdf_str.cc b/modules/io/src/mol/sdf_str.cc
index a2977c432bf3ec90647de1b05dc8eb5a2afb3f74..0441b672f0edbc3604888cec50abc35f75e0e54b 100644
--- a/modules/io/src/mol/sdf_str.cc
+++ b/modules/io/src/mol/sdf_str.cc
@@ -37,9 +37,9 @@ String EntityToSDFString(const mol::EntityView& ent) {
return stream.str();
}
-mol::EntityHandle SDFStringToEntity(const String& sdf) {
+mol::EntityHandle SDFStringToEntity(const String& sdf, const IOProfile& profile) {
std::stringstream stream(sdf);
- SDFReader reader(stream);
+ SDFReader reader(stream, profile);
mol::EntityHandle ent = mol::CreateEntity();
reader.Import(ent);
return ent;
diff --git a/modules/io/src/mol/sdf_str.hh b/modules/io/src/mol/sdf_str.hh
index 87987679ed7ad28a6e3023f536ad6c6fb716f7f7..8cc6974593a202be792cc8fcf788fae6f9b780df 100644
--- a/modules/io/src/mol/sdf_str.hh
+++ b/modules/io/src/mol/sdf_str.hh
@@ -22,6 +22,7 @@
#include <ost/io/module_config.hh>
#include <ost/mol/entity_view.hh>
#include <ost/mol/entity_handle.hh>
+#include <ost/io/mol/io_profile.hh>
namespace ost { namespace io {
@@ -33,7 +34,7 @@ String DLLEXPORT_OST_IO
EntityToSDFString(const mol::EntityView& ent);
mol::EntityHandle DLLEXPORT_OST_IO
-SDFStringToEntity(const String& pdb);
+SDFStringToEntity(const String& pdb, const IOProfile& profile);
}}
diff --git a/modules/io/src/mol/sdf_writer.cc b/modules/io/src/mol/sdf_writer.cc
index e9fd0ded4ca1fb85a9a6fd0ea1db3c7a041d0ad4..e270049f4e5f476683f076149211b87078d3dd89 100644
--- a/modules/io/src/mol/sdf_writer.cc
+++ b/modules/io/src/mol/sdf_writer.cc
@@ -22,6 +22,7 @@
#include "sdf_writer.hh"
+#include <time.h>
#include <ost/boost_filesystem_helper.hh>
#include <ost/mol/atom_view.hh>
#include <ost/mol/residue_view.hh>
@@ -251,8 +252,19 @@ bool SDFWriter::VisitChain(const mol::ChainView& chain) {
}
// print header lines
+ // line 1: molecule name
ostr_ << cname << std::endl;
- ostr_ << std::endl;
+ // line 2: program name and time
+ // IIPPPPPPPPMMDDYYHHmmddSSssssssssssEEEEEEEEEEEERRRRRR
+ std::time_t t = std::time(nullptr);
+ std::tm tm = *std::localtime(&t);
+ ostr_ << " " // User's first and last initials
+ << " OST" // program name
+ << std::put_time(&tm, "%m%d%y%H%M") // date/time
+ << "3D" // dimensional codes
+ << " " // scaling factors, energy, internal registry number
+ << std::endl;
+ // line 3: comment (blank)
ostr_ << std::endl;
// print counts line
ostr_ << format("%3d") % chain.GetAtomCount()
diff --git a/modules/io/src/mol/star_parser.hh b/modules/io/src/mol/star_parser.hh
index ba00f43980d4aa6346b74119ceef2dd4b3b1cc73..cee30f04b90d3022649da977d7f81e9da82b1431 100644
--- a/modules/io/src/mol/star_parser.hh
+++ b/modules/io/src/mol/star_parser.hh
@@ -239,6 +239,11 @@ public:
public:
static bool SplitLine(const StringRef& line,
std::vector<StringRef>& parts, bool clear=true);
+ static inline bool IsUndefined(StringRef value)
+{
+ return value.size()==1 && (value[0]=='?' || value[0]=='.');
+}
+
private:
void ParseLoop();
/// \brief Calls the loop parsing functions on the last data item fetched to
diff --git a/modules/io/tests/test_io_sdf.cc b/modules/io/tests/test_io_sdf.cc
index 68009e9c990d1cf11962eaffdcaea489934d7a7e..8217c05c104e0019c980c4b7fb612585d40c880e 100644
--- a/modules/io/tests/test_io_sdf.cc
+++ b/modules/io/tests/test_io_sdf.cc
@@ -149,7 +149,8 @@ BOOST_AUTO_TEST_CASE(write_sdf)
SaveEntity(eh, "testfiles/sdf/compound-out.sdf");
}
BOOST_CHECK(compare_files("testfiles/sdf/compound-out.sdf",
- "testfiles/sdf/compound.sdf"));
+ "testfiles/sdf/compound.sdf",
+ {2, 159, 316, 473}));
}
BOOST_AUTO_TEST_CASE(write_sdf_view)
@@ -164,7 +165,8 @@ BOOST_AUTO_TEST_CASE(write_sdf_view)
SaveEntity(ev, "testfiles/sdf/compound-view-out.sdf");
}
BOOST_CHECK(compare_files("testfiles/sdf/compound-view-out.sdf",
- "testfiles/sdf/compound-view.sdf"));
+ "testfiles/sdf/compound-view.sdf",
+ {2, 111, 220, 329}));
}
BOOST_AUTO_TEST_CASE(nonexisting_file)
diff --git a/modules/io/tests/test_io_sdf.py b/modules/io/tests/test_io_sdf.py
index 7277a399d5a58ca696f50fb19e0775345999313a..40cee207926556fea2916d0a420fcf7f421dee24 100644
--- a/modules/io/tests/test_io_sdf.py
+++ b/modules/io/tests/test_io_sdf.py
@@ -47,7 +47,41 @@ class TestSDF(unittest.TestCase):
# Charge from atom line is ignored
o_at = ent.FindAtom("00001_Simple Ligand", 1, "3")
self.assertEqual(o_at.charge, 0)
-
+
+ def test_fault_tolerant(self):
+ """This file has a "dative" bond (type = 9).
+ This is a non-standard extension from RDKit which should go through only
+ in fault tolerant mode"""
+
+ with self.assertRaises(Exception):
+ ent = io.LoadSDF('testfiles/sdf/dative_bond.sdf')
+
+ # Directly with fault_tolerant
+ PushVerbosityLevel(-1) # Expect message at Error level
+ ent = io.LoadSDF('testfiles/sdf/dative_bond.sdf', fault_tolerant=True)
+ PopVerbosityLevel()
+ self.assertEqual(ent.FindAtom("00001_Simple Ligand", 1, "5").bonds[0].bond_order, 9)
+
+ # Sloppy profile
+ PushVerbosityLevel(-1) # Expect message at Error level
+ ent = io.LoadSDF('testfiles/sdf/dative_bond.sdf', profile="SLOPPY")
+ PopVerbosityLevel()
+ self.assertEqual(ent.FindAtom("00001_Simple Ligand", 1, "5").bonds[0].bond_order, 9)
+
+ # Sloppy profile set as default
+ old_profile = io.profiles['DEFAULT'].Copy()
+ io.profiles['DEFAULT'] = "SLOPPY"
+ PushVerbosityLevel(-1) # Expect message at Error level
+ ent = io.LoadSDF('testfiles/sdf/dative_bond.sdf')
+ PopVerbosityLevel()
+ self.assertEqual(ent.FindAtom("00001_Simple Ligand", 1, "5").bonds[0].bond_order, 9)
+
+ # Test that a restored default profile has fault_tolerant again
+ io.profiles['DEFAULT'] = old_profile
+ with self.assertRaises(Exception):
+ ent = io.LoadSDF('testfiles/sdf/dative_bond.sdf')
+
+
if __name__== '__main__':
from ost import testutils
testutils.RunTests()
diff --git a/modules/io/tests/test_mmcif_writer.cc b/modules/io/tests/test_mmcif_writer.cc
index 2f6f6745bf6106a3cefdea928993ec7014675683..d85ab094906b54fc97db82b3f128ce076f8cace0 100644
--- a/modules/io/tests/test_mmcif_writer.cc
+++ b/modules/io/tests/test_mmcif_writer.cc
@@ -25,28 +25,31 @@
#include <ost/io/mol/mmcif_writer.hh>
#include <ost/mol/mol.hh>
#include <ost/platform.hh>
+#include <ost/log.hh>
using namespace ost;
+using namespace ost::conop;
using namespace ost::io;
BOOST_AUTO_TEST_SUITE( io );
-conop::CompoundLibPtr SetDefaultCompoundLib() {
- // return NULL if not successful, else return newly set default lib
- // REQ: OST_ROOT to be set
- char * ost_root=getenv("OST_ROOT");
- if (!ost_root) return conop::CompoundLibPtr();
- SetPrefixPath(ost_root);
- String lib_path=GetSharedDataPath() + "/compounds.chemlib";
- conop::CompoundLibPtr compound_lib=conop::CompoundLib::Load(lib_path);
- if (compound_lib) {
- conop::Conopology::Instance().SetDefaultLib(compound_lib);
+CompoundLibPtr load_lib()
+{
+ if (!getenv("OST_ROOT")) {
+ LOG_ERROR("OST_ROOT environment variable not set. Can't load "
+ "compound library without a proper OST_ROOT");
+ return CompoundLibPtr();
}
+ SetPrefixPath(getenv("OST_ROOT"));
+ String lib_path=GetSharedDataPath()+"/compounds.chemlib";
+ CompoundLibPtr compound_lib=CompoundLib::Load(lib_path);
return compound_lib;
}
BOOST_AUTO_TEST_CASE(mmcif_writer_force_hetatm)
{
+ CompoundLibPtr lib = load_lib();
+ if (!lib) { return; }
BOOST_TEST_MESSAGE(" Running mmcif_force_hetatm tests...");
/*
Make sure that atoms set to HETATM are written as HETATM. There is some
@@ -75,7 +78,7 @@ BOOST_AUTO_TEST_CASE(mmcif_writer_force_hetatm)
// Create mmCIF stream
MMCifWriter writer;
- writer.SetStructure(ent, SetDefaultCompoundLib(), false);
+ writer.SetStructure(ent, lib, false);
std::stringstream out;
writer.Write("test", out);
@@ -113,7 +116,7 @@ BOOST_AUTO_TEST_CASE(mmcif_writer_force_hetatm)
// Create mmCIF stream
writer=MMCifWriter();
- writer.SetStructure(ent, SetDefaultCompoundLib(), false);
+ writer.SetStructure(ent, lib, false);
out=std::stringstream();
writer.Write("test", out);
@@ -134,6 +137,8 @@ BOOST_AUTO_TEST_CASE(mmcif_writer_force_hetatm)
BOOST_AUTO_TEST_CASE(mmcif_writer_entity1)
{
+ CompoundLibPtr lib = load_lib();
+ if (!lib) { return; }
BOOST_TEST_MESSAGE(" Running mmcif_writer_entity1 tests...");
/*
Make sure molecular entities in mmCIF files written by OST start counting
@@ -151,7 +156,7 @@ BOOST_AUTO_TEST_CASE(mmcif_writer_entity1)
// Create mmCIF stream
MMCifWriter writer;
- writer.SetStructure(ent, SetDefaultCompoundLib(), false);
+ writer.SetStructure(ent, lib, false);
std::stringstream out;
writer.Write("test", out);
@@ -165,6 +170,8 @@ BOOST_AUTO_TEST_CASE(mmcif_writer_entity1)
BOOST_AUTO_TEST_CASE(mmcif_writer_poly_vs_non_poly)
{
+ CompoundLibPtr lib = load_lib();
+ if (!lib) { return; }
BOOST_TEST_MESSAGE(" Running mmcif_writer_poly_vs_non_poly tests...");
/*
Go for small polymers that are not polymer... the story of 2 amino acids (to
@@ -202,7 +209,7 @@ BOOST_AUTO_TEST_CASE(mmcif_writer_poly_vs_non_poly)
// Create mmCIF stream
MMCifWriter writer;
- writer.SetStructure(ent, SetDefaultCompoundLib(), false);
+ writer.SetStructure(ent, lib, false);
std::stringstream out;
writer.Write("test", out);
@@ -296,7 +303,7 @@ BOOST_AUTO_TEST_CASE(mmcif_writer_poly_vs_non_poly)
// Create mmCIF stream
writer=MMCifWriter();
- writer.SetStructure(ent, SetDefaultCompoundLib(), false);
+ writer.SetStructure(ent, lib, false);
out=std::stringstream();
writer.Write("test", out);
@@ -322,6 +329,8 @@ BOOST_AUTO_TEST_CASE(mmcif_writer_poly_vs_non_poly)
BOOST_AUTO_TEST_CASE(mmcif_writer_small_sugars)
{
+ CompoundLibPtr lib = load_lib();
+ if (!lib) { return; }
BOOST_TEST_MESSAGE(" Running mmcif_writer_small_sugars tests...");
/*
While RCSB marks dipeptides and dinucleotides as non-ploymers, sugars are
@@ -367,7 +376,7 @@ BOOST_AUTO_TEST_CASE(mmcif_writer_small_sugars)
// Create mmCIF stream
MMCifWriter writer;
- writer.SetStructure(ent, SetDefaultCompoundLib(), false);
+ writer.SetStructure(ent, lib, false);
std::stringstream out;
writer.Write("test", out);
@@ -407,7 +416,7 @@ BOOST_AUTO_TEST_CASE(mmcif_writer_small_sugars)
// Create mmCIF stream
writer=MMCifWriter();
- writer.SetStructure(ent, SetDefaultCompoundLib(), false);
+ writer.SetStructure(ent, lib, false);
out = std::stringstream();
writer.Write("test", out);
diff --git a/modules/io/tests/testfiles/sdf/dative_bond.sdf b/modules/io/tests/testfiles/sdf/dative_bond.sdf
new file mode 100644
index 0000000000000000000000000000000000000000..30ac15cde8ddd239f13f7c47b33f584df0533591
--- /dev/null
+++ b/modules/io/tests/testfiles/sdf/dative_bond.sdf
@@ -0,0 +1,18 @@
+Simple Ligand
+
+ Teststructure
+ 6 6 0 0 1 0 999 V2000
+ 0.0000 0.0000 0.0000 N 0 3 0 0 0 0
+ 1.0000 0.0000 0.0000 C 0 0 0 0 0 0
+ 0.0000 1.0000 0.0000 O 0 0 0 0 0 0
+ 1.0000 1.0000 0.0000 S 0 0 0 0 0 0
+ 2.0000 2.0000 0.0000 C 0 0 0 0 0 0
+ -1.0000 -1.0000 0.0000 Cl 0 0 0 0 0 0
+ 1 2 2 0 0 0
+ 1 3 1 0 0 0
+ 1 6 1 0 0 0
+ 2 4 1 0 0 0
+ 3 4 1 0 0 0
+ 4 5 9 0 0 0
+M END
+$$$$
diff --git a/modules/mol/alg/doc/qsscore.rst b/modules/mol/alg/doc/qsscore.rst
index 7b6f1711f1441ce8d7631c1bb28fab471d3568c6..4c0a7c7f9589720a6c2f76e3178e8df0e12dafd5 100644
--- a/modules/mol/alg/doc/qsscore.rst
+++ b/modules/mol/alg/doc/qsscore.rst
@@ -2,6 +2,7 @@
--------------------------------------------------------------------------------
.. module:: ost.mol.alg.qsscore
+ :noindex:
:synopsis: New QS score implementation
.. note::
diff --git a/modules/mol/alg/pymod/chain_mapping.py b/modules/mol/alg/pymod/chain_mapping.py
index b19107ba5fc0640ed6a068a27517b0acd257a0aa..d30d2bcec66bd5239de73e0b84c4536c0e154262 100644
--- a/modules/mol/alg/pymod/chain_mapping.py
+++ b/modules/mol/alg/pymod/chain_mapping.py
@@ -190,13 +190,8 @@ class ReprResult:
self._mdl_bb_pos = None
self._ref_full_bb_pos = None
self._mdl_full_bb_pos = None
- self._transform = None
+ self._superposition = None
self._superposed_mdl_bb_pos = None
- self._bb_rmsd = None
- self._gdt_8 = None
- self._gdt_4 = None
- self._gdt_2 = None
- self._gdt_1 = None
self._ost_query = None
self._flat_mapping = None
self._inconsistent_residues = None
@@ -318,24 +313,33 @@ class ReprResult:
self._mdl_full_bb_pos = self._GetFullBBPos(self.mdl_residues)
return self._mdl_full_bb_pos
+
@property
- def transform(self):
- """ Transformation to superpose mdl residues onto ref residues
+ def superposition(self):
+ """ Superposition of mdl residues onto ref residues
Superposition computed as minimal RMSD superposition on
:attr:`ref_bb_pos` and :attr:`mdl_bb_pos`. If number of positions is
smaller 3, the full_bb_pos equivalents are used instead.
- :type: :class:`ost.geom.Mat4`
+ :type: :class:`ost.mol.alg.SuperpositionResult`
"""
- if self._transform is None:
+ if self._superposition is None:
if len(self.mdl_bb_pos) < 3:
- self._transform = _GetTransform(self.mdl_full_bb_pos,
+ self._superposition = _GetSuperposition(self.mdl_full_bb_pos,
self.ref_full_bb_pos, False)
else:
- self._transform = _GetTransform(self.mdl_bb_pos,
+ self._superposition = _GetSuperposition(self.mdl_bb_pos,
self.ref_bb_pos, False)
- return self._transform
+ return self._superposition
+
+ @property
+ def transform(self):
+ """ Transformation to superpose mdl residues onto ref residues
+
+ :type: :class:`ost.geom.Mat4`
+ """
+ return self.superposition.transformation
@property
def superposed_mdl_bb_pos(self):
@@ -350,53 +354,12 @@ class ReprResult:
@property
def bb_rmsd(self):
- """ RMSD between :attr:`ref_bb_pos` and :attr:`superposed_mdl_bb_pos`
-
- :type: :class:`float`
- """
- if self._bb_rmsd is None:
- self._bb_rmsd = self.ref_bb_pos.GetRMSD(self.superposed_mdl_bb_pos)
- return self._bb_rmsd
-
- @property
- def gdt_8(self):
- """ GDT with one single threshold: 8.0
-
- :type: :class:`float`
- """
- if self._gdt_8 is None:
- self._gdt_8 = self.ref_bb_pos.GetGDT(self.superposed_mdl_bb_pos, 8.0)
- return self._gdt_8
-
- @property
- def gdt_4(self):
- """ GDT with one single threshold: 4.0
+ """ RMSD of the binding site backbone atoms after :attr:`superposition`
:type: :class:`float`
"""
- if self._gdt_4 is None:
- self._gdt_4 = self.ref_bb_pos.GetGDT(self.superposed_mdl_bb_pos, 4.0)
- return self._gdt_4
+ return self.superposition.rmsd
- @property
- def gdt_2(self):
- """ GDT with one single threshold: 2.0
-
- :type: :class:`float`
- """
- if self._gdt_2 is None:
- self._gdt_2 = self.ref_bb_pos.GetGDT(self.superposed_mdl_bb_pos, 2.0)
- return self._gdt_2
-
- @property
- def gdt_1(self):
- """ GDT with one single threshold: 1.0
-
- :type: :class:`float`
- """
- if self._gdt_1 is None:
- self._gdt_1 = self.ref_bb_pos.GetGDT(self.superposed_mdl_bb_pos, 1.0)
- return self._gdt_1
@property
def ost_query(self):
@@ -435,10 +398,6 @@ class ReprResult:
self.mdl_residues]
json_dict["transform"] = list(self.transform.data)
json_dict["bb_rmsd"] = self.bb_rmsd
- json_dict["gdt_8"] = self.gdt_8
- json_dict["gdt_4"] = self.gdt_4
- json_dict["gdt_2"] = self.gdt_2
- json_dict["gdt_1"] = self.gdt_1
json_dict["ost_query"] = self.ost_query
json_dict["flat_mapping"] = self.GetFlatChainMapping()
return json_dict
@@ -522,10 +481,9 @@ class ChainMapper:
simply derived from residue numbers (*resnum_alignments* flag).
In case of NW, *pep_subst_mat*, *pep_gap_open* and *pep_gap_ext*
and their nucleotide equivalents are relevant. Two chains are
- considered identical if they fulfill the thresholds given by
- *pep_seqid_thr*, *pep_gap_thr*, their nucleotide equivalents
- respectively. The grouping information is available as
- attributes of this class.
+ considered identical if they fulfill the *pep_seqid_thr* and
+ have at least *min_pep_length* aligned residues. Same logic
+ is applied for nucleotides with respective thresholds.
* Map chains in an input model to these groups. Generating alignments
and the similarity criteria are the same as above. You can either
@@ -551,19 +509,8 @@ class ChainMapper:
identical. 95 percent tolerates the few mutations
crystallographers like to do.
:type pep_seqid_thr: :class:`float`
- :param pep_gap_thr: Additional threshold to avoid gappy alignments with
- high seqid. By default this is disabled (set to 1.0).
- This threshold checks for a maximum allowed fraction
- of gaps in any of the two sequences after stripping
- terminal gaps. The reason for not just normalizing
- seqid by the longer sequence is that one sequence
- might be a perfect subsequence of the other but only
- cover half of it.
- :type pep_gap_thr: :class:`float`
:param nuc_seqid_thr: Nucleotide equivalent for *pep_seqid_thr*
:type nuc_seqid_thr: :class:`float`
- :param nuc_gap_thr: Nucleotide equivalent for *nuc_gap_thr*
- :type nuc_gap_thr: :class:`float`
:param pep_subst_mat: Substitution matrix to align peptide sequences,
irrelevant if *resnum_alignments* is True,
defaults to seq.alg.BLOSUM62
@@ -595,8 +542,7 @@ class ChainMapper:
:type n_max_naive: :class:`int`
"""
def __init__(self, target, resnum_alignments=False,
- pep_seqid_thr = 95., pep_gap_thr = 1.0,
- nuc_seqid_thr = 95., nuc_gap_thr = 1.0,
+ pep_seqid_thr = 95., nuc_seqid_thr = 95.,
pep_subst_mat = seq.alg.BLOSUM62, pep_gap_open = -11,
pep_gap_ext = -1, nuc_subst_mat = seq.alg.NUC44,
nuc_gap_open = -4, nuc_gap_ext = -4,
@@ -606,9 +552,7 @@ class ChainMapper:
# attributes
self.resnum_alignments = resnum_alignments
self.pep_seqid_thr = pep_seqid_thr
- self.pep_gap_thr = pep_gap_thr
self.nuc_seqid_thr = nuc_seqid_thr
- self.nuc_gap_thr = nuc_gap_thr
self.min_pep_length = min_pep_length
self.min_nuc_length = min_nuc_length
self.n_max_naive = n_max_naive
@@ -697,9 +641,9 @@ class ChainMapper:
_GetChemGroupAlignments(self.polypep_seqs, self.polynuc_seqs,
self.aligner,
pep_seqid_thr=self.pep_seqid_thr,
- pep_gap_thr=self.pep_gap_thr,
+ min_pep_length=self.min_pep_length,
nuc_seqid_thr=self.nuc_seqid_thr,
- nuc_gap_thr=self.nuc_gap_thr)
+ min_nuc_length=self.min_nuc_length)
return self._chem_group_alignments
@@ -738,9 +682,9 @@ class ChainMapper:
_GetChemGroupAlignments(self.polypep_seqs, self.polynuc_seqs,
self.aligner,
pep_seqid_thr=self.pep_seqid_thr,
- pep_gap_thr=self.pep_gap_thr,
+ min_pep_length=self.min_pep_length,
nuc_seqid_thr=self.nuc_seqid_thr,
- nuc_gap_thr=self.nuc_gap_thr)
+ min_nuc_length=self.min_nuc_length)
return self._chem_group_types
@@ -1176,7 +1120,8 @@ class ChainMapper:
for t_pos, t in zip(trg_pos, trg_chains):
for m_pos, m in zip(mdl_pos, mdl_chains):
if len(t_pos) >= 3 and len(m_pos) >= 3:
- transform = _GetTransform(m_pos, t_pos, False)
+ transform = _GetSuperposition(m_pos, t_pos,
+ False).transformation
initial_transforms.append(transform)
initial_mappings.append((t,m))
@@ -1438,6 +1383,10 @@ class ChainMapper:
substructure_chem_mapping,
self.n_max_naive))
+ # This step can be slow so give some hints in logs
+ msg = "Computing initial ranking of %d chain mappings" % len(mappings)
+ (ost.LogWarning if len(mappings) > 10000 else ost.LogInfo)(msg)
+
for mapping in mappings:
# chain_mapping and alns as input for lDDT computation
lddt_chain_mapping = dict()
@@ -1756,20 +1705,17 @@ def _GetAlnPropsOne(aln):
:param aln: Alignment to compute properties
:type aln: :class:`seq.AlignmentHandle`
- :returns: Tuple with 3 elements. 1) sequence identify in range [0, 100]
- considering aligned columns 2) Fraction of gaps between
- first and last aligned column in s1 3) same for s2.
+ :returns: Tuple with 2 elements. 1) sequence identify in range [0, 100]
+ considering aligned columns 2) Number of aligned columns.
"""
assert(aln.GetCount() == 2)
- n_gaps_1 = str(aln.GetSequence(0)).strip('-').count('-')
- n_gaps_2 = str(aln.GetSequence(1)).strip('-').count('-')
- gap_frac_1 = float(n_gaps_1)/len(aln.GetSequence(0).GetGaplessString())
- gap_frac_2 = float(n_gaps_2)/len(aln.GetSequence(1).GetGaplessString())
- return (seq.alg.SequenceIdentity(aln), gap_frac_1, gap_frac_2)
+ seqid = seq.alg.SequenceIdentity(aln)
+ n_aligned = sum([1 for col in aln if (col[0] != '-' and col[1] != '-')])
+ return (seqid, n_aligned)
def _GetChemGroupAlignments(pep_seqs, nuc_seqs, aligner, pep_seqid_thr=95.,
- pep_gap_thr=0.1, nuc_seqid_thr=95.,
- nuc_gap_thr=0.1):
+ min_pep_length=6, nuc_seqid_thr=95.,
+ min_nuc_length=4):
"""Returns alignments with groups of chemically equivalent chains
:param pep_seqs: List of polypeptide sequences
@@ -1782,17 +1728,14 @@ def _GetChemGroupAlignments(pep_seqs, nuc_seqs, aligner, pep_seqid_thr=95.,
identical. 95 percent tolerates the few mutations
crystallographers like to do.
:type pep_seqid_thr: :class:`float`
- :param pep_gap_thr: Additional threshold to avoid gappy alignments with high
- seqid. The reason for not just normalizing seqid by the
- longer sequence is that one sequence might be a perfect
- subsequence of the other but only cover half of it. This
- threshold checks for a maximum allowed fraction of gaps
- in any of the two sequences after stripping terminal gaps.
- :type pep_gap_thr: :class:`float`
+ :param min_pep_length: Additional threshold to avoid gappy alignments with high
+ seqid. Number of aligned columns must be at least this
+ number.
+ :type min_pep_length: :class:`int`
:param nuc_seqid_thr: Nucleotide equivalent of *pep_seqid_thr*
:type nuc_seqid_thr: :class:`float`
- :param nuc_gap_thr: Nucleotide equivalent of *nuc_gap_thr*
- :type nuc_gap_thr: :class:`float`
+ :param min_nuc_length: Nucleotide equivalent of *min_pep_length*
+ :type min_nuc_length: :class:`int`
:returns: Tuple with first element being an AlignmentList. Each alignment
represents a group of chemically equivalent chains and the first
sequence is the longest. Second element is a list of equivalent
@@ -1800,9 +1743,9 @@ def _GetChemGroupAlignments(pep_seqs, nuc_seqs, aligner, pep_seqid_thr=95.,
[:class:`ost.ChemType.AMINOACIDS`,
:class:`ost.ChemType.NUCLEOTIDES`]
"""
- pep_groups = _GroupSequences(pep_seqs, pep_seqid_thr, pep_gap_thr, aligner,
+ pep_groups = _GroupSequences(pep_seqs, pep_seqid_thr, min_pep_length, aligner,
mol.ChemType.AMINOACIDS)
- nuc_groups = _GroupSequences(nuc_seqs, nuc_seqid_thr, nuc_gap_thr, aligner,
+ nuc_groups = _GroupSequences(nuc_seqs, nuc_seqid_thr, min_nuc_length, aligner,
mol.ChemType.NUCLEOTIDES)
group_types = [mol.ChemType.AMINOACIDS] * len(pep_groups)
group_types += [mol.ChemType.NUCLEOTIDES] * len(nuc_groups)
@@ -1810,18 +1753,15 @@ def _GetChemGroupAlignments(pep_seqs, nuc_seqs, aligner, pep_seqid_thr=95.,
groups.extend(nuc_groups)
return (groups, group_types)
-def _GroupSequences(seqs, seqid_thr, gap_thr, aligner, chem_type):
+def _GroupSequences(seqs, seqid_thr, min_length, aligner, chem_type):
"""Get list of alignments representing groups of equivalent sequences
:param seqid_thr: Threshold used to decide when two chains are identical.
:type seqid_thr: :class:`float`
:param gap_thr: Additional threshold to avoid gappy alignments with high
- seqid. The reason for not just normalizing seqid by the
- longer sequence is that one sequence might be a perfect
- subsequence of the other but only cover half of it. This
- threshold checks for a maximum allowed fraction of gaps
- in any of the two sequences after stripping terminal gaps.
- :type gap_thr: :class:`float`
+ seqid. Number of aligned columns must be at least this
+ number.
+ :type gap_thr: :class:`int`
:param aligner: Helper class to generate pairwise alignments
:type aligner: :class:`_Aligner`
:param chem_type: ChemType of seqs which is passed to *aligner*, must be in
@@ -1839,8 +1779,8 @@ def _GroupSequences(seqs, seqid_thr, gap_thr, aligner, chem_type):
for g_s_idx in range(len(groups[g_idx])):
aln = aligner.Align(seqs[s_idx], seqs[groups[g_idx][g_s_idx]],
chem_type)
- sid, frac_i, frac_j = _GetAlnPropsOne(aln)
- if sid >= seqid_thr and frac_i < gap_thr and frac_j < gap_thr:
+ sid, n_aligned = _GetAlnPropsOne(aln)
+ if sid >= seqid_thr and n_aligned >= min_length:
matching_group = g_idx
break
if matching_group is not None:
@@ -2970,8 +2910,8 @@ def _IterativeRigidRMSD(initial_transforms, initial_mappings, chem_groups,
trg_chain_groups[best_sc_group_idx].remove(best_sc_mapping[0])
mdl_chain_groups[best_sc_group_idx].remove(best_sc_mapping[1])
- transform = _GetTransform(mapped_mdl_pos, mapped_trg_pos,
- False)
+ transform = _GetSuperposition(mapped_mdl_pos, mapped_trg_pos,
+ False).transformation
# compute overall RMSD for current transform
mapped_mdl_pos.ApplyTransform(transform)
@@ -3286,7 +3226,7 @@ def _ChainMappings(ref_chains, mdl_chains, n_max=None):
return _ConcatIterators(iterators)
-def _GetTransform(pos_one, pos_two, iterative):
+def _GetSuperposition(pos_one, pos_two, iterative):
""" Computes minimal RMSD superposition for pos_one onto pos_two
:param pos_one: Positions that should be superposed onto *pos_two*
@@ -3297,7 +3237,7 @@ def _GetTransform(pos_one, pos_two, iterative):
potentially raises, uses standard superposition as fallback.
:type iterative: :class:`bool`
:returns: Transformation matrix to superpose *pos_one* onto *pos_two*
- :rtype: :class:`geom.Mat4`
+ :rtype: :class:`ost.mol.alg.SuperpositionResult`
"""
res = None
if iterative:
@@ -3307,7 +3247,7 @@ def _GetTransform(pos_one, pos_two, iterative):
pass # triggers fallback below
if res is None:
res = mol.alg.SuperposeSVD(pos_one, pos_two)
- return res.transformation
+ return res
# specify public interface
__all__ = ('ChainMapper', 'ReprResult', 'MappingResult')
diff --git a/modules/mol/alg/pymod/dockq.py b/modules/mol/alg/pymod/dockq.py
index e8ec3f7a5e42d2bfdffa5929266ab9550aabbe14..b448e8b1e29262a097e8ba4f7556f2a6705bd1ad 100644
--- a/modules/mol/alg/pymod/dockq.py
+++ b/modules/mol/alg/pymod/dockq.py
@@ -158,7 +158,9 @@ def _RMSDScores(mdl, ref, mdl_ch1, mdl_ch2, ref_ch1, ref_ch2, dist_thresh=10.0,
cb_mode=False):
# backbone atoms used for superposition
- sup_atoms = ['CA','C','N','O']
+ sup_atoms = ["CA","C","N","O",
+ "P", "OP1", "OP2", "O2'", "O3'", "O4'", "O5'", "C1'", "C2'",
+ "C3'", "C4'", "C5'"]
# make mapped residues accessible by the dockq_idx property
mapped_mdl = mdl.Select(f"cname={mol.QueryQuoteName(mdl_ch1)},{mol.QueryQuoteName(mdl_ch2)} and grdockq_mapped=1")
@@ -215,7 +217,7 @@ def _RMSDScores(mdl, ref, mdl_ch1, mdl_ch2, ref_ch1, ref_ch2, dist_thresh=10.0,
# receptor is by definition the larger chain in ref
n_ch1 = len(ref.FindChain(ref_ch1).residues)
n_ch2 = len(ref.FindChain(ref_ch2).residues)
- if n_ch1 > n_ch2:
+ if n_ch1 >= n_ch2:
ref_receptor_residues = ref_ch1_residues.values()
ref_ligand_residues = ref_ch2_residues.values()
mdl_receptor_residues = \
diff --git a/modules/mol/alg/pymod/lddt.py b/modules/mol/alg/pymod/lddt.py
index 1453db731b9f6c680ca8ca548a3d04af183218fe..ae7fdbc2970a7c8f377bc979f0a14a1e76935899 100644
--- a/modules/mol/alg/pymod/lddt.py
+++ b/modules/mol/alg/pymod/lddt.py
@@ -1,3 +1,4 @@
+import itertools
import numpy as np
from ost import mol
@@ -438,7 +439,7 @@ class lDDTScorer:
no_intrachain=False, penalize_extra_chains=False,
residue_mapping=None, return_dist_test=False,
check_resnames=True, add_mdl_contacts=False,
- interaction_data=None):
+ interaction_data=None, set_atom_props=False):
"""Computes lDDT of *model* - globally and per-residue
:param model: Model to be scored - models are preferably scored upon
@@ -530,6 +531,12 @@ class lDDTScorer:
:type add_mdl_contacts: :class:`bool`
:param interaction_data: Pro param - don't use
:type interaction_data: :class:`tuple`
+ :param set_atom_props: If True, sets generic properties on a per atom
+ level if *local_lddt_prop*/*local_contact_prop*
+ are set as well.
+ In other words: this is the only way you can
+ get per-atom lDDT values.
+ :type set_atom_props: :class:`bool`
:returns: global and per-residue lDDT scores as a tuple -
first element is global lDDT score (None if *target* has no
@@ -610,16 +617,28 @@ class lDDTScorer:
self._ResolveSymmetries(pos, thresholds, symmetries, sym_ref_indices,
sym_ref_distances)
+ atom_indices = list(itertools.chain.from_iterable(res_atom_indices))
+
+ per_atom_exp = np.asarray([self._GetNExp(i, ref_indices)
+ for i in atom_indices], dtype=np.int32)
per_res_exp = np.asarray([self._GetNExp(res_ref_atom_indices[idx],
ref_indices) for idx in range(len(res_indices))], dtype=np.int32)
- per_res_conserved = self._EvalResidues(pos, thresholds,
- res_atom_indices,
- ref_indices, ref_distances)
+
+ per_atom_conserved = self._EvalAtoms(pos, atom_indices, thresholds,
+ ref_indices, ref_distances)
+ per_res_conserved = np.zeros((len(res_atom_indices), len(thresholds)),
+ dtype=np.int32)
+ start_idx = 0
+ for r_idx in range(len(res_atom_indices)):
+ end_idx = start_idx + len(res_atom_indices[r_idx])
+ per_res_conserved[r_idx] = np.sum(per_atom_conserved[start_idx:end_idx,:],
+ axis=0)
+ start_idx = end_idx
n_thresh = len(thresholds)
# do per-residue scores
- per_res_lDDT = [None] * len(model.residues)
+ per_res_lDDT = [None] * model.GetResidueCount()
for idx in range(len(res_indices)):
n_exp = n_thresh * per_res_exp[idx]
if n_exp > 0:
@@ -656,6 +675,40 @@ class lDDTScorer:
residues[r_idx].SetIntProp(conserved_prop,
int(np.sum(per_res_conserved[i,:])))
+ if set_atom_props and (local_lddt_prop or local_contact_prop):
+ atom_list = list()
+ residues = model.residues
+ for i, indices in enumerate(res_atom_indices):
+ r = residues[res_indices[i]]
+ r_idx = ref_res_indices[i]
+ res_start_idx = self.res_start_indices[r_idx]
+ anames = self.compound_anames[self.compound_names[r_idx]]
+ for a_i in indices:
+ a = r.FindAtom(anames[a_i - res_start_idx])
+ assert(a.IsValid())
+ atom_list.append(a)
+
+ summed_per_atom_conserved = per_atom_conserved.sum(axis=1)
+ if local_lddt_prop:
+ # the only place where actually need to compute per-atom lDDT
+ # scores
+ for a_idx in range(len(atom_list)):
+ if per_atom_exp[a_idx] != 0:
+ tmp = summed_per_atom_conserved[a_idx] / per_atom_exp[a_idx]
+ tmp = tmp / n_thresh
+ atom_list[a_idx].SetFloatProp(local_lddt_prop, tmp)
+
+ if local_contact_prop:
+ conserved_prop = local_contact_prop + "_cons"
+ exp_prop = local_contact_prop + "_exp"
+ for a_idx in range(len(atom_list)):
+ # do number of conserved contacts
+ tmp = summed_per_atom_conserved[a_idx]
+ atom_list[a_idx].SetIntProp(conserved_prop, tmp)
+ # do number of expected contacts
+ tmp = per_atom_exp[a_idx] * n_thresh
+ atom_list[a_idx].SetIntProp(exp_prop, tmp)
+
if return_dist_test:
return lDDT, per_res_lDDT, lDDT_tot, lDDT_cons, res_indices, \
per_res_exp, per_res_conserved
diff --git a/modules/mol/alg/pymod/ligand_scoring_base.py b/modules/mol/alg/pymod/ligand_scoring_base.py
index 7a38e4d3baf41fd2a14e9ab8f28b8c04661c6bbd..45376fde5f3d43a7fffb7ce7bbe518b570e8ea9a 100644
--- a/modules/mol/alg/pymod/ligand_scoring_base.py
+++ b/modules/mol/alg/pymod/ligand_scoring_base.py
@@ -1,10 +1,30 @@
+from contextlib import contextmanager
import numpy as np
import networkx
from ost import mol
-from ost import LogWarning, LogScript, LogVerbose, LogDebug
+from ost import LogWarning, LogScript, LogInfo, LogVerbose, LogDebug, GetVerbosityLevel, PushVerbosityLevel, PopVerbosityLevel
from ost.mol.alg import chain_mapping
+
+@contextmanager
+def _SinkVerbosityLevel(n=1):
+ """ Context manager to temporarily reduce the verbosity level by n.
+
+ Example usage:
+ with _SinkVerbosityLevel(2):
+ LogVerbose("Test")
+ Will only write "Test" in script level (2 above)
+ """
+ PushVerbosityLevel(GetVerbosityLevel() - n)
+ try:
+ yield
+ except:
+ raise
+ finally:
+ PopVerbosityLevel()
+
+
class LigandScorer:
""" Scorer to compute various small molecule ligand (non polymer) scores.
@@ -23,9 +43,10 @@ class LigandScorer:
* :class:`ost.mol.alg.ligand_scoring_lddtpli.LDDTPLIScorer`
that assesses the conservation of protein-ligand
- contacts
+ contacts (lDDT-PLI);
* :class:`ost.mol.alg.ligand_scoring_scrmsd.SCRMSDScorer`
- that computes a binding-site superposed, symmetry-corrected RMSD.
+ that computes a binding-site superposed, symmetry-corrected RMSD
+ (BiSyRMSD) and ligand pocket lDDT (lDDT-LP).
All versus all scores are available through the lazily computed
:attr:`score_matrix`. However, many things can go wrong... be it even
@@ -45,7 +66,7 @@ class LigandScorer:
A common use case is to derive a one-to-one mapping between ligands in
the model and the target for which :class:`LigandScorer` provides an
- automated assignment procedure.
+ automated :attr:`assignment` procedure.
By default, only exact matches between target and model ligands are
considered. This is a problem when the target only contains a subset
of the expected atoms (for instance if atoms are missing in an
@@ -65,7 +86,7 @@ class LigandScorer:
:class:`LigandScorer` generally assumes that the
:attr:`~ost.mol.ResidueHandle.is_ligand` property is properly set on all
- the ligand atoms, and only ligand atoms. This is typically the case for
+ the ligand residues, and only ligand atoms. This is typically the case for
entities loaded from mmCIF (tested with mmCIF files from the PDB and
SWISS-MODEL). Legacy PDB files must contain `HET` headers (which is usually
the case for files downloaded from the PDB but not elsewhere).
@@ -284,7 +305,7 @@ class LigandScorer:
1: ("identity", f"Ligands could not be matched (by {iso})"),
2: ("symmetries", "Too many symmetries between ligand atoms were "
"found - increasing max_symmetries might help"),
- 3: ("no_iso", "No fully isomorphic match could be found - enabling "
+ 3: ("no_iso", "No full isomorphic match could be found - enabling "
"substructure_match might allow a match"),
4: ("disconnected", "Ligand graph is disconnected"),
5: ("stoichiometry", "Ligand was already assigned to another ligand "
@@ -301,7 +322,7 @@ class LigandScorer:
respective location in this matrix is 0.
Target ligands are in rows, model ligands in columns. States are encoded
as integers <= 9. Larger numbers encode errors for child classes.
- Use something like ``self.state_decoding[3]`` to get a decscription.
+ Use something like ``self.state_decoding[3]`` to get a decscription.
:rtype: :class:`~numpy.ndarray`
"""
@@ -426,6 +447,7 @@ class LigandScorer:
raise RuntimeError("LigandScorer._score_dir must return one in "
"['+', '-']")
+ LogScript("Computing ligand assignment")
while len(tmp) > 0:
# select high coverage ligand pairs in working array
coverage_thresh = max([x[1] for x in tmp]) - self.coverage_delta
@@ -553,7 +575,27 @@ class LigandScorer:
""" Makes an educated guess why target ligand is not assigned
This either returns actual error states or custom states that are
- derived from them.
+ derived from them. Currently, the following reasons are reported:
+
+ * `no_ligand`: there was no ligand in the model.
+ * `disconnected`: the ligand graph was disconnected.
+ * `identity`: the ligand was not found in the model (by graph
+ isomorphism). Check your ligand connectivity.
+ * `no_iso`: no full isomorphic match could be found. Try enabling
+ `substructure_match=True` if the target ligand is incomplete.
+ * `symmetries`: too many symmetries were found (by graph isomorphisms).
+ Try to increase `max_symmetries`.
+ * `stoichiometry`: there was a possible assignment in the model, but
+ the model ligand was already assigned to a different target ligand.
+ This indicates different stoichiometries.
+ * `no_contact` (lDDT-PLI only): There were no lDDT contacts between
+ the binding site and the ligand, and lDDT-PLI is undefined.
+ * `target_binding_site` (SCRMSD only): no polymer residues were in
+ proximity of the target ligand.
+ * `model_binding_site` (SCRMSD only): the binding site was not found
+ in the model. Either the binding site was not modeled or the model
+ ligand was positioned too far in combination with
+ `full_bs_search=False`.
:param trg_lig_idx: Index of target ligand
:type trg_lig_idx: :class:`int`
@@ -585,7 +627,7 @@ class LigandScorer:
"Ligand was already assigned to an other "
"model ligand (different stoichiometry)")
- # maybe its a symmetry issue?
+ # maybe it's a symmetry issue?
if 2 in tmp:
return self.state_decoding[2]
@@ -593,12 +635,12 @@ class LigandScorer:
# target counterpart.
if 6 in tmp:
mdl_idx = np.where(self.state_matrix[trg_lig_idx,:]==6)[0]
- # we're reporting everything except disconnected error...
- # don't ask...
for i in mdl_idx:
if self.model_ligand_states[i] == 0:
raise RuntimeError("This should never happen")
- if self.model_ligand_states[i] != 4:
+ if self.model_ligand_states[i] != 4 or len(tmp) == 1:
+ # Don't report disconnected if only 1 model ligand is
+ # disconnected, unless that's the only reason
return self.state_decoding[self.model_ligand_states[i]]
# get rid of remaining single ligand issues (only disconnected error)
@@ -617,7 +659,28 @@ class LigandScorer:
""" Makes an educated guess why model ligand is not assigned
This either returns actual error states or custom states that are
- derived from them.
+ derived from them. Currently, the following reasons are reported:
+
+ * `no_ligand`: there was no ligand in the target.
+ * `disconnected`: the ligand graph is disconnected.
+ * `identity`: the ligand was not found in the target (by graph or
+ subgraph isomorphism). Check your ligand connectivity.
+ * `no_iso`: no full isomorphic match could be found. Try enabling
+ `substructure_match=True` if the target ligand is incomplete.
+ * `symmetries`: too many symmetries were found (by graph isomorphisms).
+ Try to increase `max_symmetries`.
+ * `stoichiometry`: there was a possible assignment in the target, but
+ the model target was already assigned to a different model ligand.
+ This indicates different stoichiometries.
+ * `no_contact` (lDDT-PLI only): There were no lDDT contacts between
+ the binding site and the ligand, and lDDT-PLI is undefined.
+ * `target_binding_site` (SCRMSD only): a potential assignment was found
+ in the target, but there were no polymer residues in proximity of the
+ ligand in the target.
+ * `model_binding_site` (SCRMSD only): a potential assignment was
+ found in the target, but no binding site was found in the model.
+ Either the binding site was not modeled or the model ligand was
+ positioned too far in combination with `full_bs_search=False`.
:param mdl_lig_idx: Index of model ligand
:type mdl_lig_idx: :class:`int`
@@ -657,12 +720,12 @@ class LigandScorer:
# target counterpart.
if 6 in tmp:
trg_idx = np.where(self.state_matrix[:,mdl_lig_idx]==6)[0]
- # we're reporting everything except disconnected error...
- # don't ask...
for i in trg_idx:
if self.target_ligand_states[i] == 0:
raise RuntimeError("This should never happen")
- if self.target_ligand_states[i] != 4:
+ if self.target_ligand_states[i] != 4 or len(tmp) == 1:
+ # Don't report disconnected if only 1 target ligand is
+ # disconnected, unless that's the only reason
return self.state_decoding[self.target_ligand_states[i]]
# get rid of remaining single ligand issues (only disconnected error)
@@ -712,10 +775,11 @@ class LigandScorer:
:type: :class:`ost.mol.alg.chain_mapping.ChainMapper`
"""
if self.__chain_mapper is None:
- self.__chain_mapper = \
- chain_mapping.ChainMapper(self.target,
- n_max_naive=1e9,
- resnum_alignments=self.resnum_alignments)
+ with _SinkVerbosityLevel():
+ self.__chain_mapper = \
+ chain_mapping.ChainMapper(self.target,
+ n_max_naive=1e9,
+ resnum_alignments=self.resnum_alignments)
return self.__chain_mapper
@staticmethod
@@ -772,7 +836,7 @@ class LigandScorer:
# Instantiate the editor
if new_editor is None:
- new_editor = new_entity.EditXCS()
+ new_editor = new_entity.EditXCS(mol.BUFFERED_EDIT)
new_chain = new_entity.FindChain(residue.chain.name)
if not new_chain.IsValid():
@@ -796,7 +860,7 @@ class LigandScorer:
new_chain = \
new_editor.InsertChain(new_chain_name)
break
- LogScript("Moved ligand residue %s to new chain %s" % (
+ LogInfo("Moved ligand residue %s to new chain %s" % (
residue.qualified_name, new_chain.name))
else:
msg = \
@@ -887,7 +951,7 @@ class LigandScorer:
self._state_matrix[g_idx,:] = 6
msg = "Disconnected graph observed for target ligand "
msg += str(self.target_ligands[g_idx])
- LogVerbose(msg)
+ LogWarning(msg)
for g_idx, g in enumerate(model_graphs):
if not networkx.is_connected(g):
@@ -895,11 +959,12 @@ class LigandScorer:
self._state_matrix[:,g_idx] = 6
msg = "Disconnected graph observed for model ligand "
msg += str(self.model_ligands[g_idx])
- LogVerbose(msg)
+ LogWarning(msg)
+ LogScript("Computing pairwise scores for all %s x %s ligands" % shape)
for target_id, target_ligand in enumerate(self.target_ligands):
- LogVerbose("Analyzing target ligand %s" % target_ligand)
+ LogInfo("Analyzing target ligand %s" % target_ligand)
if self._target_ligand_states[target_id] == 4:
# Disconnected graph - already updated states and reported
@@ -907,7 +972,7 @@ class LigandScorer:
continue
for model_id, model_ligand in enumerate(self.model_ligands):
- LogVerbose("Compare to model ligand %s" % model_ligand)
+ LogInfo("Comparing to model ligand %s" % model_ligand)
#########################################################
# Compute symmetries for given target/model ligand pair #
@@ -926,30 +991,30 @@ class LigandScorer:
max_symmetries=self.max_symmetries,
model_graph = model_graphs[model_id],
target_graph = target_graphs[target_id])
- LogVerbose("Ligands %s and %s symmetry match" % (
+ LogInfo("Ligands %s and %s symmetry match" % (
str(model_ligand), str(target_ligand)))
except NoSymmetryError:
# Ligands are different - skip
- LogVerbose("No symmetry between %s and %s" % (
+ LogInfo("No symmetry between %s and %s" % (
str(model_ligand), str(target_ligand)))
self._state_matrix[target_id, model_id] = 1
continue
except TooManySymmetriesError:
# Ligands are too symmetrical - skip
- LogVerbose("Too many symmetries between %s and %s" % (
+ LogWarning("Too many symmetries between %s and %s" % (
str(model_ligand), str(target_ligand)))
self._state_matrix[target_id, model_id] = 2
continue
except NoIsomorphicSymmetryError:
# Ligands are different - skip
- LogVerbose("No isomorphic symmetry between %s and %s" % (
+ LogInfo("No isomorphic symmetry between %s and %s" % (
str(model_ligand), str(target_ligand)))
self._state_matrix[target_id, model_id] = 3
continue
except DisconnectedGraphError:
# this should never happen as we guard against
# DisconnectedGraphError when precomputing the graph
- LogVerbose("Disconnected graph observed for %s and %s" % (
+ LogError("Disconnected graph observed for %s and %s" % (
str(model_ligand), str(target_ligand)))
# just set both ligand states to 4
self._model_ligand_states[model_id] = 4
@@ -1161,7 +1226,7 @@ def ComputeSymmetries(model_ligand, target_ligand, substructure_match=False,
symmetries.append((list(isomorphism.values()),
list(isomorphism.keys())))
assert len(symmetries) > 0
- LogDebug("Found %s isomorphic mappings (symmetries)" % len(symmetries))
+ LogVerbose("Found %s isomorphic mappings (symmetries)" % len(symmetries))
elif gm.subgraph_is_isomorphic() and substructure_match:
if not return_symmetries:
return True
@@ -1177,14 +1242,14 @@ def ComputeSymmetries(model_ligand, target_ligand, substructure_match=False,
# Assert that all the atoms in the target are part of the substructure
assert len(symmetries[0][0]) == len(target_ligand.atoms)
n_sym = len(symmetries)
- LogDebug("Found %s subgraph isomorphisms (symmetries)" % n_sym)
+ LogVerbose("Found %s subgraph isomorphisms (symmetries)" % n_sym)
elif gm.subgraph_is_isomorphic():
- LogDebug("Found subgraph isomorphisms (symmetries), but"
+ LogVerbose("Found subgraph isomorphisms (symmetries), but"
" ignoring because substructure_match=False")
raise NoIsomorphicSymmetryError("No symmetry between %s and %s" % (
str(model_ligand), str(target_ligand)))
else:
- LogDebug("Found no isomorphic mappings (symmetries)")
+ LogVerbose("Found no isomorphic mappings (symmetries)")
raise NoSymmetryError("No symmetry between %s and %s" % (
str(model_ligand), str(target_ligand)))
diff --git a/modules/mol/alg/pymod/ligand_scoring_lddtpli.py b/modules/mol/alg/pymod/ligand_scoring_lddtpli.py
index 5710e3efd07d91da727424a59f050e4f39c18ba9..467066471a480227480f98f4be4adae065242b22 100644
--- a/modules/mol/alg/pymod/ligand_scoring_lddtpli.py
+++ b/modules/mol/alg/pymod/ligand_scoring_lddtpli.py
@@ -1,6 +1,6 @@
import numpy as np
-from ost import LogWarning
+from ost import LogWarning, LogInfo
from ost import geom
from ost import mol
from ost import seq
@@ -26,20 +26,16 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
and an lDDT-PLI score is computed. The best possible lDDT-PLI score is
returned.
- By default, classic lDDT is computed. That means, contacts within
- *lddt_pli_radius* are identified in the target and checked if they're
- conserved in the model. Added contacts are not penalized. That means if
- the ligand is nicely placed in the correct pocket, but that pocket now
- suddenly interacts with MORE residues in the model, you still get a high
- score. You can penalize for these added contacts with the
- *add_mdl_contacts* flag. This additionally considers contacts within
- *lddt_pli_radius* in the model but only if the involved atoms can
- be mapped to the target. This is a requirement to 1) extract the respective
- reference distance from the target 2) avoid usage of contacts for which
- we have no experimental evidence. One special case are
- contacts from chains that are NOT mapped to the target binding site. It is
- very well possible that we have experimental evidence for this chain though
- its just too far away from the target binding site.
+ The lDDT-PLI score is a variant of lDDT with a custom inclusion radius
+ (`lddt_pli_radius`), no stereochemistry checks, and which penalizes
+ contacts added in the model within `lddt_pli_radius` by default
+ (can be changed with the `add_mdl_contacts` flag) but only if the involved
+ atoms can be mapped to the target. This is a requirement to
+ 1) extract the respective reference distance from the target
+ 2) avoid usage of contacts for which we have no experimental evidence.
+ One special case are contacts from chains that are not mapped to the target
+ binding site. It is very well possible that we have experimental evidence
+ for this chain though its just too far away from the target binding site.
We therefore try to map these contacts to the chain in the target with
equivalent sequence that is closest to the target binding site. If the
respective atoms can be mapped there, the contact is considered not
@@ -47,7 +43,7 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
Populates :attr:`LigandScorer.aux_data` with following :class:`dict` keys:
- * lddt_pli: The score
+ * lddt_pli: The LDDT-PLI score
* lddt_pli_n_contacts: Number of contacts considered in lDDT computation
* target_ligand: The actual target ligand for which the score was computed
* model_ligand: The actual model ligand for which the score was computed
@@ -84,7 +80,7 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
:type max_symmetries: :class:`int`
:param lddt_pli_radius: lDDT inclusion radius for lDDT-PLI.
:type lddt_pli_radius: :class:`float`
- :param add_mdl_contacts: Whether to add mdl contacts.
+ :param add_mdl_contacts: Whether to penalize added model contacts.
:type add_mdl_contacts: :class:`bool`
:param lddt_pli_thresholds: Distance difference thresholds for lDDT.
:type lddt_pli_thresholds: :class:`list` of :class:`float`
@@ -102,8 +98,8 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
def __init__(self, model, target, model_ligands=None, target_ligands=None,
resnum_alignments=False, rename_ligand_chain=False,
substructure_match=False, coverage_delta=0.2,
- max_symmetries=1e5, lddt_pli_radius=6.0,
- add_mdl_contacts=False,
+ max_symmetries=1e4, lddt_pli_radius=6.0,
+ add_mdl_contacts=True,
lddt_pli_thresholds = [0.5, 1.0, 2.0, 4.0],
lddt_pli_binding_site_radius=None):
@@ -141,10 +137,12 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
""" Implements interface from parent
"""
if self.add_mdl_contacts:
+ LogInfo("Computing lDDT-PLI with added model contacts")
result = self._compute_lddt_pli_add_mdl_contacts(symmetries,
target_ligand,
model_ligand)
else:
+ LogInfo("Computing lDDT-PLI without added model contacts")
result = self._compute_lddt_pli_classic(symmetries,
target_ligand,
model_ligand)
@@ -935,10 +933,11 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
@property
def _chain_mapping_mdl(self):
- if self.__chain_mapping_mdl is None:
- self.__chem_mapping, self.__chem_group_alns, \
- self.__chain_mapping_mdl = \
- self._chain_mapper.GetChemMapping(self.model)
+ if self.__chain_mapping_mdl is None:
+ with ligand_scoring_base._SinkVerbosityLevel():
+ self.__chem_mapping, self.__chem_group_alns, \
+ self.__chain_mapping_mdl = \
+ self._chain_mapper.GetChemMapping(self.model)
return self.__chain_mapping_mdl
# specify public interface
diff --git a/modules/mol/alg/pymod/ligand_scoring_scrmsd.py b/modules/mol/alg/pymod/ligand_scoring_scrmsd.py
index 56ef24a9e9cccf354e2b286aac1794ae40526ae2..15ef917b48eb4bbcd4c0d97fcf3900ceb8c88d40 100644
--- a/modules/mol/alg/pymod/ligand_scoring_scrmsd.py
+++ b/modules/mol/alg/pymod/ligand_scoring_scrmsd.py
@@ -1,20 +1,21 @@
import numpy as np
-from ost import LogWarning
+from ost import LogWarning, LogScript, LogInfo, LogVerbose
from ost import geom
from ost import mol
from ost.mol.alg import ligand_scoring_base
+
class SCRMSDScorer(ligand_scoring_base.LigandScorer):
- """ :class:`LigandScorer` implementing symmetry corrected RMSD.
+ """ :class:`LigandScorer` implementing symmetry corrected RMSD (BiSyRMSD).
:class:`SCRMSDScorer` computes a score for a specific pair of target/model
ligands.
The returned RMSD is based on a binding site superposition.
The binding site of the target structure is defined as all residues with at
- least one atom within *bs_radius* around the target ligand.
+ least one atom within `bs_radius` around the target ligand.
It only contains protein and nucleic acid residues from chains that
pass the criteria for the
:class:`chain mapping <ost.mol.alg.chain_mapping>`. This means ignoring
@@ -23,28 +24,29 @@ class SCRMSDScorer(ligand_scoring_base.LigandScorer):
The respective model binding site for superposition is identified by
naively enumerating all possible mappings of model chains onto their
chemically equivalent target counterparts from the target binding site.
- The *binding_sites_topn* with respect to lDDT score are evaluated and
+ The `binding_sites_topn` with respect to lDDT score are evaluated and
an RMSD is computed.
You can either try to map ALL model chains onto the target binding site by
- enabling *full_bs_search* or restrict the model chains for a specific
+ enabling `full_bs_search` or restrict the model chains for a specific
target/model ligand pair to the chains with at least one atom within
*model_bs_radius* around the model ligand. The latter can be significantly
faster in case of large complexes.
Symmetry correction is achieved by simply computing an RMSD value for
each symmetry, i.e. atom-atom assignments of the ligand as given by
- :class:`LigandScorer`. The lowest RMSD value is returned
+ :class:`LigandScorer`. The lowest RMSD value is returned.
Populates :attr:`LigandScorer.aux_data` with following :class:`dict` keys:
- * rmsd: The score
- * lddt_lp: lDDT of the binding site used for superposition
+ * rmsd: The BiSyRMSD score
+ * lddt_lp: lDDT of the binding pocket used for superposition (lDDT-LP)
* bs_ref_res: :class:`list` of binding site residues in target
* bs_ref_res_mapped: :class:`list` of target binding site residues that
are mapped to model
* bs_mdl_res_mapped: :class:`list` of same length with respective model
residues
- * bb_rmsd: Backbone RMSD (CA, C3' for nucleotides) for mapped residues
- given transform
+ * bb_rmsd: Backbone RMSD (CA, C3' for nucleotides; full backbone for
+ binding sites with fewer than 3 residues) for mapped binding site
+ residues after superposition
* target_ligand: The actual target ligand for which the score was computed
* model_ligand: The actual model ligand for which the score was computed
* chain_mapping: :class:`dict` with a chain mapping of chains involved in
@@ -145,14 +147,13 @@ class SCRMSDScorer(ligand_scoring_base.LigandScorer):
self._get_repr_input = dict()
# update state decoding from parent with subclass specific stuff
- self.state_decoding[10] = ("binding_site",
+ self.state_decoding[10] = ("target_binding_site",
"No residues were in proximity of the "
"target ligand.")
- self.state_decoding[11] = ("model_representation", "No representation "
- "of the reference binding site was found in "
- "the model, i.e. the binding site was not "
- "modeled or the model ligand was positioned "
- "too far in combination with "
+ self.state_decoding[11] = ("model_binding_site", "Binding site was not"
+ " found in the model, i.e. the binding site"
+ " was not modeled or the model ligand was "
+ "positioned too far in combination with "
"full_bs_search=False.")
self.state_decoding[20] = ("unknown",
"Unknown error occured in SCRMSDScorer")
@@ -174,6 +175,9 @@ class SCRMSDScorer(ligand_scoring_base.LigandScorer):
"inconsistent_residues": list()}
representations = self._get_repr(target_ligand, model_ligand)
+ # This step can be slow so give some hints in logs
+ msg = "Computing BiSyRMSD with %d chain mappings" % len(representations)
+ (LogWarning if len(representations) > 10000 else LogInfo)(msg)
for r in representations:
rmsd = _SCRMSD_symmetries(symmetries, model_ligand,
@@ -231,6 +235,8 @@ class SCRMSDScorer(ligand_scoring_base.LigandScorer):
if key not in self._repr:
ref_bs = self._get_target_binding_site(target_ligand)
+ LogVerbose("%d chains are in proximity of the target ligand: %s" % (
+ ref_bs.chain_count, ", ".join([c.name for c in ref_bs.chains])))
if self.full_bs_search:
reprs = self._chain_mapper.GetRepr(
ref_bs, self.model, inclusion_radius=self.lddt_lp_radius,
@@ -262,7 +268,8 @@ class SCRMSDScorer(ligand_scoring_base.LigandScorer):
h = ref_res.handle.GetHashCode()
if h not in ref_residues_hashes and \
h not in ignored_residue_hashes:
- view = self._chain_mapper.target.ViewForHandle(ref_res)
+ with ligand_scoring_base._SinkVerbosityLevel(1):
+ view = self._chain_mapper.target.ViewForHandle(ref_res)
if view.IsValid():
h = ref_res.handle.GetHashCode()
ref_residues_hashes.add(h)
@@ -341,12 +348,15 @@ class SCRMSDScorer(ligand_scoring_base.LigandScorer):
radius = self.model_bs_radius
chains = set()
for at in mdl_ligand.atoms:
- close_atoms = self._chain_mapping_mdl.FindWithin(at.GetPos(),
- radius)
+ with ligand_scoring_base._SinkVerbosityLevel(1):
+ close_atoms = self._chain_mapping_mdl.FindWithin(at.GetPos(),
+ radius)
for close_at in close_atoms:
chains.add(close_at.GetChain().GetName())
if len(chains) > 0:
+ LogVerbose("%d chains are in proximity of the model ligand: %s" % (
+ len(chains), ", ".join(chains)))
# the chain mapping model which only contains close chains
query = "cname="
diff --git a/modules/mol/alg/pymod/qsscore.py b/modules/mol/alg/pymod/qsscore.py
index 2c98d9fa71cf936854cf6838ad779f634372ea5c..e7d2d945411da662d4706bc4a0681de3a170bb09 100644
--- a/modules/mol/alg/pymod/qsscore.py
+++ b/modules/mol/alg/pymod/qsscore.py
@@ -621,7 +621,7 @@ class QSScorer:
mapped_idx_grid_1 = np.ix_(mapped_indices_1_1, mapped_indices_2_1)
mapped_idx_grid_2 = np.ix_(mapped_indices_1_2, mapped_indices_2_2)
- if mapped_idx_grid_1[0].shape[0] == 0 or mapped_idx_grid_2[0].shape[0] == 0:
+ if mapped_indices_1_1.shape[0] == 0 or mapped_indices_2_1.shape[0] == 0:
# dealing with special cases where we have no mapped residues
# we only avoid errors here when using maped_idx_grid_x for indexing
# but run the rest of the algorithm anyways which produces some
@@ -630,7 +630,8 @@ class QSScorer:
shared_mask_d2 = np.full(d2.shape, False, dtype=bool)
mapped_nonshared_mask_d1 = np.full(d1.shape, False, dtype=bool)
mapped_nonshared_mask_d2 = np.full(d2.shape, False, dtype=bool)
- if mapped_idx_grid_1[0].shape[0] == 0:
+ if mapped_indices_1_1.shape[0] == 0 or \
+ mapped_indices_2_1.shape[0] == 0:
# mapped_idx_grid_1 has not a single mapped residue which raises
# an error when calling something like d1[mapped_idx_grid_1]
mapped_d1_contacts = np.full(d1.shape, False, dtype=bool)
@@ -638,7 +639,8 @@ class QSScorer:
mapped_d1_contacts = d1[mapped_idx_grid_1] < contact_d
mapped_nonshared_mask_d1[mapped_idx_grid_1] = mapped_d1_contacts
- if mapped_idx_grid_2[0].shape[0] == 0:
+ if mapped_indices_1_2.shape[0] == 0 or \
+ mapped_indices_2_2.shape[0] == 0:
# mapped_idx_grid_2 has not a single mapped residue which raises
# an error when calling something like d2[mapped_idx_grid_2]
mapped_d2_contacts = np.full(d2.shape, False, dtype=bool)
diff --git a/modules/mol/alg/pymod/scoring.py b/modules/mol/alg/pymod/scoring.py
index 4444e463709703042b6395da3c866804cd971237..0626bff67ad383dccbcbefcdd5f849abd652e69a 100644
--- a/modules/mol/alg/pymod/scoring.py
+++ b/modules/mol/alg/pymod/scoring.py
@@ -96,11 +96,11 @@ class Scorer:
Deals with structure cleanup, chain mapping, interface identification etc.
Intermediate results are available as attributes.
- :param model: Model structure - a deep copy is available as :attr:`model`.
+ :param model: Model structure - a deep copy is available as :attr:`~model`.
Additionally, :func:`ost.mol.alg.Molck` using *molck_settings*
is applied.
:type model: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
- :param target: Target structure - a deep copy is available as :attr:`target`.
+ :param target: Target structure - a deep copy is available as :attr:`~target`.
Additionally, :func:`ost.mol.alg.Molck` using *molck_settings*
is applied.
:type target: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
@@ -201,6 +201,10 @@ class Scorer:
This flag has no influence on patch_dockq
scores.
:type dockq_capri_peptide: :class:`bool`
+ :param lddt_symmetry_settings: Passed as *symmetry_settings* parameter to
+ lDDT scorer. Default: None
+ :type lddt_symmetry_settings: :class:`ost.mol.alg.lddt.SymmetrySettings`
+ :param lddt_inclusion_radius: lDDT inclusion radius.
"""
def __init__(self, model, target, resnum_alignments=False,
molck_settings = None, cad_score_exec = None,
@@ -208,7 +212,8 @@ class Scorer:
usalign_exec = None, lddt_no_stereochecks=False,
n_max_naive=40320, oum=False, min_pep_length = 6,
min_nuc_length = 4, lddt_add_mdl_contacts=False,
- dockq_capri_peptide=False):
+ dockq_capri_peptide=False, lddt_symmetry_settings = None,
+ lddt_inclusion_radius = 15.0):
self._target_orig = target
self._model_orig = model
@@ -235,6 +240,20 @@ class Scorer:
LogScript("Cleaning up input structures")
Molck(self._model, conop.GetDefaultLib(), molck_settings)
Molck(self._target, conop.GetDefaultLib(), molck_settings)
+
+ if resnum_alignments:
+ # If we're dealing with resnum alignments, we ensure that
+ # consecutive peptide and nucleotide residues are connected based
+ # on residue number information. The conop.Processor only connects
+ # these things if the bonds are actually feasible which can lead to
+ # weird behaviour in stereochemistry checks. Let's say N and C are
+ # too close, it's reported as a clash. If they're too far apart,
+ # they're not reported at all. If we're not dealing with resnum
+ # alignments, we're out of luck as we have no direct residue
+ # connectivity information from residue numbers
+ self._resnum_connect(self._model)
+ self._resnum_connect(self._target)
+
self._model = self._model.Select("peptide=True or nucleotide=True")
self._target = self._target.Select("peptide=True or nucleotide=True")
@@ -298,6 +317,8 @@ class Scorer:
self.min_nuc_length = min_nuc_length
self.lddt_add_mdl_contacts = lddt_add_mdl_contacts
self.dockq_capri_peptide = dockq_capri_peptide
+ self.lddt_symmetry_settings = lddt_symmetry_settings
+ self.lddt_inclusion_radius = lddt_inclusion_radius
# lazily evaluated attributes
self._stereochecked_model = None
@@ -308,6 +329,7 @@ class Scorer:
self._target_clashes = None
self._target_bad_bonds = None
self._target_bad_angles = None
+ self._trimmed_model = None
self._chain_mapper = None
self._mapping = None
self._rigid_mapping = None
@@ -316,16 +338,19 @@ class Scorer:
self._aln = None
self._stereochecked_aln = None
self._pepnuc_aln = None
+ self._trimmed_aln = None
# lazily constructed scorer objects
self._lddt_scorer = None
self._bb_lddt_scorer = None
self._qs_scorer = None
self._contact_scorer = None
+ self._trimmed_contact_scorer = None
# lazily computed scores
self._lddt = None
self._local_lddt = None
+ self._aa_local_lddt = None
self._bb_lddt = None
self._bb_local_lddt = None
self._ilddt = None
@@ -342,6 +367,7 @@ class Scorer:
self._contact_model_interfaces = None
self._native_contacts = None
self._model_contacts = None
+ self._trimmed_model_contacts = None
self._ics_precision = None
self._ics_recall = None
self._ics = None
@@ -351,9 +377,25 @@ class Scorer:
self._ips_precision = None
self._ips_recall = None
self._ips = None
- self._per_interface_ics_precision = None
- self._per_interface_ics_recall = None
- self._per_interface_ics = None
+ self._per_interface_ips_precision = None
+ self._per_interface_ips_recall = None
+ self._per_interface_ips = None
+
+ # subset of contact scores that operate on trimmed model
+ # i.e. no contacts from model residues that are not present in
+ # target
+ self._ics_trimmed = None
+ self._ics_precision_trimmed = None
+ self._ics_recall_trimmed = None
+ self._per_interface_ics_precision_trimmed = None
+ self._per_interface_ics_recall_trimmed = None
+ self._per_interface_ics_trimmed = None
+ self._ips_trimmed = None
+ self._ips_precision_trimmed = None
+ self._ips_recall_trimmed = None
+ self._per_interface_ips_precision_trimmed = None
+ self._per_interface_ips_recall_trimmed = None
+ self._per_interface_ips_trimmed = None
self._dockq_target_interfaces = None
self._dockq_interfaces = None
@@ -371,16 +413,21 @@ class Scorer:
self._mapped_target_pos = None
self._mapped_model_pos = None
+ self._mapped_target_pos_full_bb = None
+ self._mapped_model_pos_full_bb = None
self._transformed_mapped_model_pos = None
self._n_target_not_mapped = None
self._transform = None
self._rigid_mapped_target_pos = None
self._rigid_mapped_model_pos = None
+ self._rigid_mapped_target_pos_full_bb = None
+ self._rigid_mapped_model_pos_full_bb = None
self._rigid_transformed_mapped_model_pos = None
self._rigid_n_target_not_mapped = None
self._rigid_transform = None
+ self._gdt_window_sizes = [7, 9, 12, 24, 48]
self._gdt_05 = None
self._gdt_1 = None
self._gdt_2 = None
@@ -441,9 +488,9 @@ class Scorer:
@property
def aln(self):
- """ Alignments of :attr:`model`/:attr:`target` chains
+ """ Alignments of :attr:`~model`/:attr:`~target` chains
- Alignments for each pair of chains mapped in :attr:`mapping`.
+ Alignments for each pair of chains mapped in :attr:`~mapping`.
First sequence is target sequence, second sequence the model sequence.
:type: :class:`list` of :class:`ost.seq.AlignmentHandle`
@@ -454,7 +501,7 @@ class Scorer:
@property
def stereochecked_aln(self):
- """ Stereochecked equivalent of :attr:`aln`
+ """ Stereochecked equivalent of :attr:`~aln`
The alignments may differ, as stereochecks potentially remove residues
@@ -466,7 +513,7 @@ class Scorer:
@property
def pepnuc_aln(self):
- """ Alignments of :attr:`model_orig`/:attr:`target_orig` chains
+ """ Alignments of :attr:`~model_orig`/:attr:`~target_orig` chains
Selects for peptide and nucleotide residues before sequence
extraction. Includes residues that would be removed by molck in
@@ -478,6 +525,19 @@ class Scorer:
self._compute_pepnuc_aln()
return self._pepnuc_aln
+ @property
+ def trimmed_aln(self):
+ """ Alignments of :attr:`~trimmed_model`/:attr:`~target` chains
+
+ Alignments for each pair of chains mapped in :attr:`~mapping`.
+ First sequence is target sequence, second sequence the model sequence.
+
+ :type: :class:`list` of :class:`ost.seq.AlignmentHandle`
+ """
+ if self._trimmed_aln is None:
+ self._trim_model()
+ return self._trimmed_aln
+
@property
def stereochecked_model(self):
""" View of :attr:`~model` that has stereochemistry checks applied
@@ -568,9 +628,23 @@ class Scorer:
self._do_stereochecks()
return self._target_bad_angles
+ @property
+ def trimmed_model(self):
+ """ :attr:`~model` trimmed to target
+
+ Removes residues that are not covered by :class:`target` given
+ :attr:`~mapping`. In other words: no model residues without experimental
+ evidence from :class:`target`.
+
+ :type: :class:`ost.mol.EntityView`
+ """
+ if self._trimmed_model is None:
+ self._trim_model()
+ return self._trimmed_model
+
@property
def chain_mapper(self):
- """ Chain mapper object for given :attr:`target`
+ """ Chain mapper object for given :attr:`~target`
:type: :class:`ost.mol.alg.chain_mapping.ChainMapper`
"""
@@ -584,7 +658,7 @@ class Scorer:
@property
def mapping(self):
- """ Full chain mapping result for :attr:`target`/:attr:`model`
+ """ Full chain mapping result for :attr:`~target`/:attr:`~model`
Computed with :func:`ost.mol.alg.ChainMapper.GetMapping`
@@ -599,7 +673,7 @@ class Scorer:
@property
def rigid_mapping(self):
- """ Full chain mapping result for :attr:`target`/:attr:`model`
+ """ Full chain mapping result for :attr:`~target`/:attr:`~model`
Computed with :func:`ost.mol.alg.ChainMapper.GetRMSDMapping`
@@ -640,15 +714,22 @@ class Scorer:
@property
def lddt_scorer(self):
- """ lDDT scorer for :attr:`~stereochecked_target` (default parameters)
+ """ lDDT scorer for :attr:`~target`/:attr:`~stereochecked_target`
+
+ Depending on :attr:`~lddt_no_stereocheck` and
+ :attr:`~lddt_symmetry_settings`.
:type: :class:`ost.mol.alg.lddt.lDDTScorer`
"""
if self._lddt_scorer is None:
if self.lddt_no_stereochecks:
- self._lddt_scorer = lDDTScorer(self.target)
+ self._lddt_scorer = lDDTScorer(self.target,
+ symmetry_settings = self.lddt_symmetry_settings,
+ inclusion_radius = self.lddt_inclusion_radius)
else:
- self._lddt_scorer = lDDTScorer(self.stereochecked_target)
+ self._lddt_scorer = lDDTScorer(self.stereochecked_target,
+ symmetry_settings = self.lddt_symmetry_settings,
+ inclusion_radius = self.lddt_inclusion_radius)
return self._lddt_scorer
@property
@@ -661,7 +742,9 @@ class Scorer:
:type: :class:`ost.mol.alg.lddt.lDDTScorer`
"""
if self._bb_lddt_scorer is None:
- self._bb_lddt_scorer = lDDTScorer(self.target, bb_only=True)
+ self._bb_lddt_scorer = lDDTScorer(self.target, bb_only=True,
+ symmetry_settings = self.lddt_symmetry_settings,
+ inclusion_radius = self.lddt_inclusion_radius)
return self._bb_lddt_scorer
@property
@@ -683,6 +766,14 @@ class Scorer:
self._contact_scorer = ContactScorer.FromMappingResult(self.mapping)
return self._contact_scorer
+ @property
+ def trimmed_contact_scorer(self):
+ if self._trimmed_contact_scorer is None:
+ self._trimmed_contact_scorer = ContactScorer(self.mapping.target,
+ self.mapping.chem_groups,
+ self.trimmed_model,
+ self.trimmed_aln)
+ return self._trimmed_contact_scorer
@property
def lddt(self):
@@ -714,6 +805,23 @@ class Scorer:
self._compute_lddt()
return self._local_lddt
+ @property
+ def aa_local_lddt(self):
+ """ Per atom lDDT scores in range [0.0, 1.0]
+
+ Computed based on :attr:`~stereochecked_model` but scores for all
+ atoms in :attr:`~model` are reported. If an atom has been removed
+ by stereochemistry checks, the respective score is set to 0.0. If an
+ atom is not covered by the target or is in a chain skipped by the
+ chain mapping procedure (happens for super short chains), the respective
+ score is set to None. In case of oligomers, :attr:`~mapping` is used.
+
+ :type: :class:`dict`
+ """
+ if self._aa_local_lddt is None:
+ self._compute_lddt()
+ return self._aa_local_lddt
+
@property
def bb_lddt(self):
""" Backbone only global lDDT score in range [0.0, 1.0]
@@ -766,7 +874,7 @@ class Scorer:
def qs_global(self):
""" Global QS-score
- Computed based on :attr:`model` using :attr:`mapping`
+ Computed based on :attr:`~model` using :attr:`~mapping`
:type: :class:`float`
"""
@@ -778,7 +886,7 @@ class Scorer:
def qs_best(self):
""" Global QS-score - only computed on aligned residues
- Computed based on :attr:`model` using :attr:`mapping`. The QS-score
+ Computed based on :attr:`~model` using :attr:`~mapping`. The QS-score
computation only considers contacts between residues with a mapping
between target and model. As a result, the score won't be lowered in
case of additional chains/residues in any of the structures.
@@ -844,7 +952,7 @@ class Scorer:
@property
def per_interface_qs_global(self):
- """ QS-score for each interface in :attr:`qs_interfaces`
+ """ QS-score for each interface in :attr:`~qs_interfaces`
:type: :class:`list` of :class:`float`
"""
@@ -854,7 +962,7 @@ class Scorer:
@property
def per_interface_qs_best(self):
- """ QS-score for each interface in :attr:`qs_interfaces`
+ """ QS-score for each interface in :attr:`~qs_interfaces`
Only computed on aligned residues
@@ -881,12 +989,20 @@ class Scorer:
@property
def model_contacts(self):
- """ Same for model
+ """ Same for :attr:`~model`
"""
if self._model_contacts is None:
self._model_contacts = self.contact_scorer.cent2.hr_contacts
return self._model_contacts
+ @property
+ def trimmed_model_contacts(self):
+ """ Same for :attr:`~trimmed_model`
+ """
+ if self._trimmed_model_contacts is None:
+ self._trimmed_model_contacts = self.trimmed_contact_scorer.cent2.hr_contacts
+ return self._trimmed_model_contacts
+
@property
def contact_target_interfaces(self):
""" Interfaces in :class:`target` which have at least one contact
@@ -992,7 +1108,6 @@ class Scorer:
if self._per_interface_ics is None:
self._compute_ics_scores()
return self._per_interface_ics
-
@property
def ips_precision(self):
@@ -1029,6 +1144,131 @@ class Scorer:
self._compute_ips_scores()
return self._ips
+ @property
+ def ics_trimmed(self):
+ """ Same as :attr:`~ics` but with trimmed model
+
+ Model is trimmed to residues which can me mapped to target in order
+ to not penalize contacts in the model for which we have no experimental
+ evidence.
+
+ :type: :class:`float`
+ """
+ if self._ics_trimmed is None:
+ self._compute_ics_scores_trimmed()
+ return self._ics_trimmed
+
+ @property
+ def ics_precision_trimmed(self):
+ """ Same as :attr:`~ics_precision` but with trimmed model
+
+ Model is trimmed to residues which can me mapped to target in order
+ to not penalize contacts in the model for which we have no experimental
+ evidence.
+
+ :type: :class:`float`
+ """
+ if self._ics_precision_trimmed is None:
+ self._compute_ics_scores_trimmed()
+ return self._ics_precision_trimmed
+
+ @property
+ def ics_recall_trimmed(self):
+ """ Same as :attr:`~ics_recall` but with trimmed model
+
+ Model is trimmed to residues which can me mapped to target in order
+ to not penalize contacts in the model for which we have no experimental
+ evidence.
+
+ :type: :class:`float`
+ """
+ if self._ics_recall_trimmed is None:
+ self._compute_ics_scores_trimmed()
+ return self._ics_recall_trimmed
+
+ @property
+ def per_interface_ics_precision_trimmed(self):
+ """ Same as :attr:`~per_interface_ics_precision` but with :attr:`~trimmed_model`
+
+ :attr:`~ics_precision_trimmed` for each interface in
+ :attr:`~contact_target_interfaces`
+
+ :type: :class:`list` of :class:`float`
+ """
+ if self._per_interface_ics_precision_trimmed is None:
+ self._compute_ics_scores_trimmed()
+ return self._per_interface_ics_precision_trimmed
+
+
+ @property
+ def per_interface_ics_recall_trimmed(self):
+ """ Same as :attr:`~per_interface_ics_recall` but with :attr:`~trimmed_model`
+
+ :attr:`~ics_recall_trimmed` for each interface in
+ :attr:`~contact_target_interfaces`
+
+ :type: :class:`list` of :class:`float`
+ """
+ if self._per_interface_ics_recall_trimmed is None:
+ self._compute_ics_scores_trimmed()
+ return self._per_interface_ics_recall_trimmed
+
+ @property
+ def per_interface_ics_trimmed(self):
+ """ Same as :attr:`~per_interface_ics` but with :attr:`~trimmed_model`
+
+ :attr:`~ics` for each interface in
+ :attr:`~contact_target_interfaces`
+
+ :type: :class:`float`
+ """
+
+ if self._per_interface_ics_trimmed is None:
+ self._compute_ics_scores_trimmed()
+ return self._per_interface_ics_trimmed
+
+ @property
+ def ips_trimmed(self):
+ """ Same as :attr:`~ips` but with trimmed model
+
+ Model is trimmed to residues which can me mapped to target in order
+ to not penalize contacts in the model for which we have no experimental
+ evidence.
+
+ :type: :class:`float`
+ """
+ if self._ips_trimmed is None:
+ self._compute_ips_scores_trimmed()
+ return self._ips_trimmed
+
+ @property
+ def ips_precision_trimmed(self):
+ """ Same as :attr:`~ips_precision` but with trimmed model
+
+ Model is trimmed to residues which can me mapped to target in order
+ to not penalize contacts in the model for which we have no experimental
+ evidence.
+
+ :type: :class:`float`
+ """
+ if self._ips_precision_trimmed is None:
+ self._compute_ips_scores_trimmed()
+ return self._ips_precision_trimmed
+
+ @property
+ def ips_recall_trimmed(self):
+ """ Same as :attr:`~ips_recall` but with trimmed model
+
+ Model is trimmed to residues which can me mapped to target in order
+ to not penalize contacts in the model for which we have no experimental
+ evidence.
+
+ :type: :class:`float`
+ """
+ if self._ips_recall_trimmed is None:
+ self._compute_ips_scores_trimmed()
+ return self._ips_recall_trimmed
+
@property
def per_interface_ips_precision(self):
""" Per-interface IPS precision
@@ -1042,7 +1282,6 @@ class Scorer:
self._compute_ips_scores()
return self._per_interface_ips_precision
-
@property
def per_interface_ips_recall(self):
""" Per-interface IPS recall
@@ -1070,12 +1309,55 @@ class Scorer:
self._compute_ips_scores()
return self._per_interface_ips
+ @property
+ def per_interface_ips_precision_trimmed(self):
+ """ Same as :attr:`~per_interface_ips_precision` but with :attr:`~trimmed_model`
+
+ :attr:`~ips_precision_trimmed` for each interface in
+ :attr:`~contact_target_interfaces`
+
+ :type: :class:`list` of :class:`float`
+ """
+ if self._per_interface_ips_precision_trimmed is None:
+ self._compute_ips_scores_trimmed()
+ return self._per_interface_ips_precision_trimmed
+
+
+ @property
+ def per_interface_ips_recall_trimmed(self):
+ """ Same as :attr:`~per_interface_ips_recall` but with :attr:`~trimmed_model`
+
+ :attr:`~ics_recall_trimmed` for each interface in
+ :attr:`~contact_target_interfaces`
+
+ :type: :class:`list` of :class:`float`
+ """
+ if self._per_interface_ips_recall_trimmed is None:
+ self._compute_ips_scores_trimmed()
+ return self._per_interface_ips_recall_trimmed
+
+ @property
+ def per_interface_ips_trimmed(self):
+ """ Same as :attr:`~per_interface_ips` but with :attr:`~trimmed_model`
+
+ :attr:`~ics` for each interface in
+ :attr:`~contact_target_interfaces`
+
+ :type: :class:`float`
+ """
+
+ if self._per_interface_ips_trimmed is None:
+ self._compute_ips_scores_trimmed()
+ return self._per_interface_ips_trimmed
+
@property
def dockq_target_interfaces(self):
- """ Interfaces in :attr:`target` that are relevant for DockQ
+ """ Interfaces in :attr:`~target` that are relevant for DockQ
All interfaces in :attr:`~target` with non-zero contacts that are
- relevant for DockQ. Chain names are lexicographically sorted.
+ relevant for DockQ. Includes protein-protein, protein-nucleotide and
+ nucleotide-nucleotide interfaces. Chain names for each interface are
+ lexicographically sorted.
:type: :class:`list` of :class:`tuple` with 2 elements each:
(trg_ch1, trg_ch2)
@@ -1093,18 +1375,20 @@ class Scorer:
interfaces = cent.interacting_chains
interfaces = [(min(x[0],x[1]), max(x[0],x[1])) for x in interfaces]
- # select the ones with only peptides involved
pep_seqs = set([s.GetName() for s in self.chain_mapper.polypep_seqs])
+ nuc_seqs = set([s.GetName() for s in self.chain_mapper.polynuc_seqs])
+
+ seqs = pep_seqs.union(nuc_seqs)
self._dockq_target_interfaces = list()
for interface in interfaces:
- if interface[0] in pep_seqs and interface[1] in pep_seqs:
+ if interface[0] in seqs and interface[1] in seqs:
self._dockq_target_interfaces.append(interface)
return self._dockq_target_interfaces
@property
def dockq_interfaces(self):
- """ Interfaces in :attr:`dockq_target_interfaces` that can be mapped
+ """ Interfaces in :attr:`~dockq_target_interfaces` that can be mapped
to model
Target chain names are lexicographically sorted
@@ -1180,9 +1464,11 @@ class Scorer:
def irmsd(self):
""" irmsd scores for interfaces in :attr:`~dockq_interfaces`
- irmsd: RMSD of interface (RMSD computed on N, CA, C, O atoms) which
+ irmsd: RMSD of interface (RMSD computed on backbone atoms) which
consists of each residue that has at least one heavy atom within 10A of
- other chain.
+ other chain. Backbone atoms for proteins: "CA","C","N","O", for
+ nucleotides: "P", "OP1", "OP2", "O2'", "O3'", "O4'", "O5'", "C1'",
+ "C2'", "C3'", "C4'", "C5'".
:class:`list` of :class:`float`
"""
@@ -1194,11 +1480,10 @@ class Scorer:
def lrmsd(self):
""" lrmsd scores for interfaces in :attr:`~dockq_interfaces`
- lrmsd: The interfaces are superposed based on the receptor (rigid
- min RMSD superposition) and RMSD for the ligand is reported.
- Superposition and RMSD are based on N, CA, C and O positions,
- receptor is the chain contributing to the interface with more
- residues in total.
+ lrmsd: The two chains involved in the interface are superposed based on
+ the receptor (rigid min RMSD superposition) and the ligand RMSD is
+ reported. Receptor is the chain with more residues. Superposition and
+ RMSD is computed on same backbone atoms as :attr:`~irmsd`.
:class:`list` of :class:`float`
"""
@@ -1208,7 +1493,7 @@ class Scorer:
@property
def dockq_ave(self):
- """ Average of DockQ scores in :attr:`dockq_scores`
+ """ Average of DockQ scores in :attr:`~dockq_scores`
In its original implementation, DockQ only operates on single
interfaces. Thus the requirement to combine scores for higher order
@@ -1222,7 +1507,7 @@ class Scorer:
@property
def dockq_wave(self):
- """ Same as :attr:`dockq_ave`, weighted by native contacts
+ """ Same as :attr:`~dockq_ave`, weighted by native contacts
:type: :class:`float`
"""
@@ -1269,6 +1554,20 @@ class Scorer:
self._extract_mapped_pos()
return self._mapped_target_pos
+ @property
+ def mapped_target_pos_full_bb(self):
+ """ Mapped representative positions in target
+
+ Thats the equivalent of :attr:`~mapped_target_pos` but containing more
+ backbone atoms (N, CA, C for peptide residues and O5', C5', C4', C3', O3
+ for nucleotide residues). mapping is based on :attr:`~mapping`.
+
+ :type: :class:`ost.geom.Vec3List`
+ """
+ if self._mapped_target_pos_full_bb is None:
+ self._extract_mapped_pos_full_bb()
+ return self._mapped_target_pos_full_bb
+
@property
def mapped_model_pos(self):
""" Mapped representative positions in model
@@ -1283,6 +1582,20 @@ class Scorer:
self._extract_mapped_pos()
return self._mapped_model_pos
+ @property
+ def mapped_model_pos_full_bb(self):
+ """ Mapped representative positions in model
+
+ Thats the equivalent of :attr:`~mapped_model_pos` but containing more
+ backbone atoms (N, CA, C for peptide residues and O5', C5', C4', C3', O3
+ for nucleotide residues). mapping is based on :attr:`~mapping`.
+
+ :type: :class:`ost.geom.Vec3List`
+ """
+ if self._mapped_model_pos_full_bb is None:
+ self._extract_mapped_pos_full_bb()
+ return self._mapped_model_pos_full_bb
+
@property
def transformed_mapped_model_pos(self):
""" :attr:`~mapped_model_pos` with :attr:`~transform` applied
@@ -1309,17 +1622,25 @@ class Scorer:
def transform(self):
""" Transform: :attr:`~mapped_model_pos` onto :attr:`~mapped_target_pos`
- Computed using Kabsch minimal rmsd algorithm
+ Computed using Kabsch minimal rmsd algorithm. If number of positions
+ is too small (< 3), :attr:`~mapped_model_pos_full_bb` and
+ :attr:`~mapped_target_pos_full_bb` are used.
:type: :class:`ost.geom.Mat4`
"""
if self._transform is None:
- try:
+ if len(self.mapped_model_pos) < 3:
+ if len(self.mapped_model_pos_full_bb) >=3:
+ res = mol.alg.SuperposeSVD(self.mapped_model_pos_full_bb,
+ self.mapped_target_pos_full_bb)
+ self._transform = res.transformation
+ else:
+ # there is really nothing we can do => set identity matrix
+ self._transform = geom.Mat4()
+ else:
res = mol.alg.SuperposeSVD(self.mapped_model_pos,
self.mapped_target_pos)
self._transform = res.transformation
- except:
- self._transform = geom.Mat4()
return self._transform
@property
@@ -1336,6 +1657,20 @@ class Scorer:
self._extract_rigid_mapped_pos()
return self._rigid_mapped_target_pos
+ @property
+ def rigid_mapped_target_pos_full_bb(self):
+ """ Mapped representative positions in target
+
+ Thats the equivalent of :attr:`~rigid_mapped_target_pos` but containing
+ more backbone atoms (N, CA, C for peptide residues and O5', C5', C4',
+ C3', O3 for nucleotide residues). mapping is based on :attr:`~mapping`.
+
+ :type: :class:`ost.geom.Vec3List`
+ """
+ if self._rigid_mapped_target_pos_full_bb is None:
+ self._extract_rigid_mapped_pos_full_bb()
+ return self._rigid_mapped_target_pos_full_bb
+
@property
def rigid_mapped_model_pos(self):
""" Mapped representative positions in model
@@ -1350,6 +1685,20 @@ class Scorer:
self._extract_rigid_mapped_pos()
return self._rigid_mapped_model_pos
+ @property
+ def rigid_mapped_model_pos_full_bb(self):
+ """ Mapped representative positions in model
+
+ Thats the equivalent of :attr:`~rigid_mapped_model_pos` but containing
+ more backbone atoms (N, CA, C for peptide residues and O5', C5', C4',
+ C3', O3 for nucleotide residues). mapping is based on :attr:`~mapping`.
+
+ :type: :class:`ost.geom.Vec3List`
+ """
+ if self._rigid_mapped_model_pos_full_bb is None:
+ self._extract_rigid_mapped_pos_full_bb()
+ return self._rigid_mapped_model_pos_full_bb
+
@property
def rigid_transformed_mapped_model_pos(self):
""" :attr:`~rigid_mapped_model_pos` with :attr:`~rigid_transform` applied
@@ -1376,17 +1725,25 @@ class Scorer:
def rigid_transform(self):
""" Transform: :attr:`~rigid_mapped_model_pos` onto :attr:`~rigid_mapped_target_pos`
- Computed using Kabsch minimal rmsd algorithm
+ Computed using Kabsch minimal rmsd algorithm. If number of positions
+ is too small (< 3), :attr:`~rigid_mapped_model_pos_full_bb` and
+ :attr:`~rigid_mapped_target_pos_full_bb` are used.
:type: :class:`ost.geom.Mat4`
"""
if self._rigid_transform is None:
- try:
+ if len(self.rigid_mapped_model_pos) < 3:
+ if len(self.rigid_mapped_model_pos_full_bb) >= 3:
+ res = mol.alg.SuperposeSVD(self.rigid_mapped_model_pos_full_bb,
+ self.rigid_mapped_target_pos_full_bb)
+ self._rigid_transform = res.transformation
+ else:
+ # there is really nothing we can do => set identity matrix
+ self._rigid_transform = geom.Mat4()
+ else:
res = mol.alg.SuperposeSVD(self.rigid_mapped_model_pos,
self.rigid_mapped_target_pos)
self._rigid_transform = res.transformation
- except:
- self._rigid_transform = geom.Mat4()
return self._rigid_transform
@property
@@ -1399,7 +1756,14 @@ class Scorer:
:type: :class:`float`
"""
if self._gdt_05 is None:
- n, m = GDT(self.rigid_mapped_model_pos, self.rigid_mapped_target_pos, 7, 1000, 0.5)
+ N = list()
+ wsizes = self._gdt_window_sizes + [len(self.rigid_mapped_model_pos)]
+ for window_size in wsizes:
+ n = GDT(self.rigid_mapped_model_pos,
+ self.rigid_mapped_target_pos,
+ window_size, 1000, 0.5)[0]
+ N.append(n)
+ n = max(N)
n_full = len(self.rigid_mapped_target_pos) + self.rigid_n_target_not_mapped
if n_full > 0:
self._gdt_05 = float(n) / n_full
@@ -1417,7 +1781,14 @@ class Scorer:
:type: :class:`float`
"""
if self._gdt_1 is None:
- n, m = GDT(self.rigid_mapped_model_pos, self.rigid_mapped_target_pos, 7, 1000, 1.0)
+ N = list()
+ wsizes = self._gdt_window_sizes + [len(self.rigid_mapped_model_pos)]
+ for window_size in wsizes:
+ n = GDT(self.rigid_mapped_model_pos,
+ self.rigid_mapped_target_pos,
+ window_size, 1000, 1.0)[0]
+ N.append(n)
+ n = max(N)
n_full = len(self.rigid_mapped_target_pos) + self.rigid_n_target_not_mapped
if n_full > 0:
self._gdt_1 = float(n) / n_full
@@ -1436,7 +1807,14 @@ class Scorer:
:type: :class:`float`
"""
if self._gdt_2 is None:
- n, m = GDT(self.rigid_mapped_model_pos, self.rigid_mapped_target_pos, 7, 1000, 2.0)
+ N = list()
+ wsizes = self._gdt_window_sizes + [len(self.rigid_mapped_model_pos)]
+ for window_size in wsizes:
+ n = GDT(self.rigid_mapped_model_pos,
+ self.rigid_mapped_target_pos,
+ window_size, 1000, 2.0)[0]
+ N.append(n)
+ n = max(N)
n_full = len(self.rigid_mapped_target_pos) + self.rigid_n_target_not_mapped
if n_full > 0:
self._gdt_2 = float(n) / n_full
@@ -1454,7 +1832,14 @@ class Scorer:
:type: :class:`float`
"""
if self._gdt_4 is None:
- n, m = GDT(self.rigid_mapped_model_pos, self.rigid_mapped_target_pos, 7, 1000, 4.0)
+ N = list()
+ wsizes = self._gdt_window_sizes + [len(self.rigid_mapped_model_pos)]
+ for window_size in wsizes:
+ n = GDT(self.rigid_mapped_model_pos,
+ self.rigid_mapped_target_pos,
+ window_size, 1000, 4.0)[0]
+ N.append(n)
+ n = max(N)
n_full = len(self.rigid_mapped_target_pos) + self.rigid_n_target_not_mapped
if n_full > 0:
self._gdt_4 = float(n) / n_full
@@ -1471,7 +1856,14 @@ class Scorer:
:type: :class:`float`
"""
if self._gdt_8 is None:
- n, m = GDT(self.rigid_mapped_model_pos, self.rigid_mapped_target_pos, 7, 1000, 8.0)
+ N = list()
+ wsizes = self._gdt_window_sizes + [len(self.rigid_mapped_model_pos)]
+ for window_size in wsizes:
+ n = GDT(self.rigid_mapped_model_pos,
+ self.rigid_mapped_target_pos,
+ window_size, 1000, 8.0)[0]
+ N.append(n)
+ n = max(N)
n_full = len(self.rigid_mapped_target_pos) + self.rigid_n_target_not_mapped
if n_full > 0:
self._gdt_8 = float(n) / n_full
@@ -1507,7 +1899,7 @@ class Scorer:
""" RMSD
Computed on :attr:`~rigid_transformed_mapped_model_pos` and
- :attr:`rigid_mapped_target_pos`
+ :attr:`~rigid_mapped_target_pos`
:type: :class:`float`
"""
@@ -1545,7 +1937,7 @@ class Scorer:
@property
def patch_qs(self):
- """ Patch QS-scores for each residue in :attr:`model_interface_residues`
+ """ Patch QS-scores for each residue in :attr:`~model_interface_residues`
Representative patches for each residue r in chain c are computed as
follows:
@@ -1561,11 +1953,11 @@ class Scorer:
mdl_patch_two
Results are stored in the same manner as
- :attr:`model_interface_residues`, with corresponding scores instead of
+ :attr:`~model_interface_residues`, with corresponding scores instead of
residue numbers. Scores for residues which are not
:class:`mol.ChemType.AMINOACIDS` are set to None. Additionally,
- interface patches are derived from :attr:`model`. If they contain
- residues which are not covered by :attr:`target`, the score is set to
+ interface patches are derived from :attr:`~model`. If they contain
+ residues which are not covered by :attr:`~target`, the score is set to
None too.
:type: :class:`dict` with chain names as key and and :class:`list`
@@ -1577,7 +1969,7 @@ class Scorer:
@property
def patch_dockq(self):
- """ Same as :attr:`patch_qs` but for DockQ scores
+ """ Same as :attr:`~patch_qs` but for DockQ scores
"""
if self._patch_dockq is None:
self._compute_patchdockq_scores()
@@ -1683,6 +2075,7 @@ class Scorer:
# score variables to be set
lddt_score = None
local_lddt = None
+ aa_local_lddt = None
if self.lddt_no_stereochecks:
lddt_chain_mapping = dict()
@@ -1694,19 +2087,39 @@ class Scorer:
residue_mapping = alns,
check_resnames=False,
local_lddt_prop="lddt",
- add_mdl_contacts = self.lddt_add_mdl_contacts)[0]
+ add_mdl_contacts = self.lddt_add_mdl_contacts,
+ set_atom_props=True)[0]
local_lddt = dict()
+ aa_local_lddt = dict()
for r in self.model.residues:
+
cname = r.GetChain().GetName()
if cname not in local_lddt:
local_lddt[cname] = dict()
+ aa_local_lddt[cname] = dict()
+
+ rnum = r.GetNumber()
+ if rnum not in aa_local_lddt[cname]:
+ aa_local_lddt[cname][rnum] = dict()
+
if r.HasProp("lddt"):
score = round(r.GetFloatProp("lddt"), 3)
- local_lddt[cname][r.GetNumber()] = score
+ local_lddt[cname][rnum] = score
else:
# not covered by trg or skipped in chain mapping procedure
# the latter happens if its part of a super short chain
- local_lddt[cname][r.GetNumber()] = None
+ local_lddt[cname][rnum] = None
+
+ for a in r.atoms:
+ if a.HasProp("lddt"):
+ score = round(a.GetFloatProp("lddt"), 3)
+ aa_local_lddt[cname][rnum][a.GetName()] = score
+ else:
+ # not covered by trg or skipped in chain mapping
+ # procedure the latter happens if its part of a
+ # super short chain
+ aa_local_lddt[cname][rnum][a.GetName()] = None
+
else:
lddt_chain_mapping = dict()
@@ -1718,22 +2131,77 @@ class Scorer:
residue_mapping = stereochecked_alns,
check_resnames=False,
local_lddt_prop="lddt",
- add_mdl_contacts = self.lddt_add_mdl_contacts)[0]
+ add_mdl_contacts = self.lddt_add_mdl_contacts,
+ set_atom_props=True)[0]
local_lddt = dict()
+ aa_local_lddt = dict()
for r in self.model.residues:
cname = r.GetChain().GetName()
if cname not in local_lddt:
local_lddt[cname] = dict()
+ aa_local_lddt[cname] = dict()
+ rnum = r.GetNumber()
+ if rnum not in aa_local_lddt[cname]:
+ aa_local_lddt[cname][rnum] = dict()
+
if r.HasProp("lddt"):
score = round(r.GetFloatProp("lddt"), 3)
- local_lddt[cname][r.GetNumber()] = score
+ local_lddt[cname][rnum] = score
+
+ trg_r = None
+ mdl_r = None
+
+ for a in r.atoms:
+ if a.HasProp("lddt"):
+ score = round(a.GetFloatProp("lddt"), 3)
+ aa_local_lddt[cname][rnum][a.GetName()] = score
+ else:
+ # the target residue is there since we have a score
+ # for the residue.
+ # opt 1: The atom was never there in the
+ # stereochecked target => None
+ # opt 2: The atom has been removed in the model
+ # stereochecks but is there in stereochecked
+ # target => 0.0
+ if trg_r is None:
+ if cname in flat_mapping:
+ for col in alns[cname]:
+ if col[0] != '-' and col[1] != '-':
+ if col.GetResidue(1).number == r.number:
+ trg_r = col.GetResidue(0)
+ break
+ if trg_r is not None:
+ trg_cname = trg_r.GetChain().GetName()
+ trg_rnum = trg_r.GetNumber()
+ tmp = self.stereochecked_target.FindResidue(trg_cname,
+ trg_rnum)
+ if tmp.IsValid():
+ trg_r = tmp
+
+ if mdl_r is None:
+ tmp = self.stereochecked_model.FindResidue(cname, rnum)
+ if tmp.IsValid():
+ mdl_r = tmp
+
+ if trg_r is not None and not trg_r.FindAtom(a.GetName()).IsValid():
+ # opt 1
+ aa_local_lddt[cname][rnum][a.GetName()] = None
+ elif trg_r is not None and trg_r.FindAtom(a.GetName()).IsValid() and \
+ mdl_r is not None and not mdl_r.FindAtom(a.GetName()).IsValid():
+ # opt 2
+ aa_local_lddt[cname][rnum][a.GetName()] = 0.0
+ else:
+ # unknown issue
+ aa_local_lddt[cname][rnum][a.GetName()] = None
+
else:
- # rsc => residue stereo checked...
- mdl_res = self.stereochecked_model.FindResidue(cname, r.GetNumber())
+ mdl_res = self.stereochecked_model.FindResidue(cname, rnum)
if mdl_res.IsValid():
# not covered by trg or skipped in chain mapping procedure
# the latter happens if its part of a super short chain
- local_lddt[cname][r.GetNumber()] = None
+ local_lddt[cname][rnum] = None
+ for a in r.atoms:
+ aa_local_lddt[cname][rnum][a.GetName()] = None
else:
# opt 1: removed by stereochecks => assign 0.0
# opt 2: removed by stereochecks AND not covered by ref
@@ -1746,13 +2214,29 @@ class Scorer:
if col.GetResidue(1).number == r.number:
trg_r = col.GetResidue(0)
break
+ if trg_r is not None:
+ trg_cname = trg_r.GetChain().GetName()
+ trg_rnum = trg_r.GetNumber()
+ tmp = self.stereochecked_target.FindResidue(trg_cname,
+ trg_rnum)
+ if tmp.IsValid():
+ trg_r = tmp
+
if trg_r is None:
- local_lddt[cname][r.GetNumber()] = None
+ local_lddt[cname][rnum] = None
+ for a in r.atoms:
+ aa_local_lddt[cname][rnum][a.GetName()] = None
else:
- local_lddt[cname][r.GetNumber()] = 0.0
+ local_lddt[cname][rnum] = 0.0
+ for a in r.atoms:
+ if trg_r.FindAtom(a.GetName()).IsValid():
+ aa_local_lddt[cname][rnum][a.GetName()] = 0.0
+ else:
+ aa_local_lddt[cname][rnum][a.GetName()] = None
self._lddt = lddt_score
self._local_lddt = local_lddt
+ self._aa_local_lddt = aa_local_lddt
def _compute_bb_lddt(self):
LogScript("Computing backbone lDDT")
@@ -1877,6 +2361,71 @@ class Scorer:
self._per_interface_ics_recall.append(None)
self._per_interface_ics.append(None)
+ def _trim_model(self):
+ trimmed_mdl = mol.CreateEntityFromView(self.mapping.model, True)
+ trimmed_aln = dict()
+
+ for trg_cname, mdl_cname in self.mapping.GetFlatMapping().items():
+ mdl_ch = trimmed_mdl.FindChain(mdl_cname)
+ aln = self.mapping.alns[(trg_cname, mdl_cname)]
+
+ # some limited test that stuff matches
+ assert(mdl_ch.GetResidueCount() == \
+ len(aln.GetSequence(1).GetGaplessString()))
+
+ mdl_residues = mdl_ch.residues
+ mdl_res_idx = 0
+ aligned_mdl_seq = ['-'] * aln.GetLength()
+ for col_idx, col in enumerate(aln):
+ if col[0] != '-' and col[1] != '-':
+ mdl_res = mdl_residues[mdl_res_idx]
+ mdl_res.SetIntProp("aligned", 1)
+ aligned_mdl_seq[col_idx] = col[1]
+ if col[1] != '-':
+ mdl_res_idx += 1
+ aligned_mdl_seq = ''.join(aligned_mdl_seq)
+ aligned_mdl_seq = seq.CreateSequence(mdl_cname, aligned_mdl_seq)
+
+ new_aln = seq.CreateAlignment()
+ new_aln.AddSequence(aln.GetSequence(0))
+ new_aln.AddSequence(aligned_mdl_seq)
+ trimmed_aln[(trg_cname, mdl_cname)] = new_aln
+
+ self._trimmed_model = trimmed_mdl.Select("graligned:0=1")
+ self._trimmed_aln = trimmed_aln
+
+ def _compute_ics_scores_trimmed(self):
+ LogScript("Computing ICS scores trimmed")
+
+ # this is an ugly hack without any efficiency in mind
+ # we're simply taking the entities from mapper and construct
+ # a new contact scorer from scratch
+
+ contact_scorer_res = self.trimmed_contact_scorer.ScoreICS(self.mapping.mapping)
+ self._ics_trimmed = contact_scorer_res.ics
+ self._ics_precision_trimmed = contact_scorer_res.precision
+ self._ics_recall_trimmed = contact_scorer_res.recall
+
+ self._per_interface_ics_precision_trimmed = list()
+ self._per_interface_ics_recall_trimmed = list()
+ self._per_interface_ics_trimmed = list()
+ flat_mapping = self.mapping.GetFlatMapping()
+ for trg_int in self.contact_target_interfaces:
+ trg_ch1 = trg_int[0]
+ trg_ch2 = trg_int[1]
+ if trg_ch1 in flat_mapping and trg_ch2 in flat_mapping:
+ mdl_ch1 = flat_mapping[trg_ch1]
+ mdl_ch2 = flat_mapping[trg_ch2]
+ res = self.trimmed_contact_scorer.ScoreICSInterface(trg_ch1, trg_ch2,
+ mdl_ch1, mdl_ch2)
+ self._per_interface_ics_precision_trimmed.append(res.precision)
+ self._per_interface_ics_recall_trimmed.append(res.recall)
+ self._per_interface_ics_trimmed.append(res.ics)
+ else:
+ self._per_interface_ics_precision_trimmed.append(None)
+ self._per_interface_ics_recall_trimmed.append(None)
+ self._per_interface_ics_trimmed.append(None)
+
def _compute_ips_scores(self):
LogScript("Computing IPS scores")
contact_scorer_res = self.contact_scorer.ScoreIPS(self.mapping.mapping)
@@ -1904,8 +2453,45 @@ class Scorer:
self._per_interface_ips_recall.append(None)
self._per_interface_ips.append(None)
+ def _compute_ips_scores_trimmed(self):
+ LogScript("Computing IPS scores trimmed")
+
+ # this is an ugly hack without any efficiency in mind
+ # we're simply taking the entities from mapper and construct
+ # a new contact scorer from scratch
+ contact_scorer_res = self.trimmed_contact_scorer.ScoreIPS(self.mapping.mapping)
+ self._ips_precision_trimmed = contact_scorer_res.precision
+ self._ips_recall_trimmed = contact_scorer_res.recall
+ self._ips_trimmed = contact_scorer_res.ips
+
+ self._per_interface_ips_precision_trimmed = list()
+ self._per_interface_ips_recall_trimmed = list()
+ self._per_interface_ips_trimmed = list()
+ flat_mapping = self.mapping.GetFlatMapping()
+ for trg_int in self.contact_target_interfaces:
+ trg_ch1 = trg_int[0]
+ trg_ch2 = trg_int[1]
+ if trg_ch1 in flat_mapping and trg_ch2 in flat_mapping:
+ mdl_ch1 = flat_mapping[trg_ch1]
+ mdl_ch2 = flat_mapping[trg_ch2]
+ res = self.trimmed_contact_scorer.ScoreIPSInterface(trg_ch1, trg_ch2,
+ mdl_ch1, mdl_ch2)
+ self._per_interface_ips_precision_trimmed.append(res.precision)
+ self._per_interface_ips_recall_trimmed.append(res.recall)
+ self._per_interface_ips_trimmed.append(res.ips)
+ else:
+ self._per_interface_ips_precision_trimmed.append(None)
+ self._per_interface_ips_recall_trimmed.append(None)
+ self._per_interface_ips_trimmed.append(None)
+
def _compute_dockq_scores(self):
LogScript("Computing DockQ")
+
+ if self.dockq_capri_peptide and len(self.chain_mapper.polynuc_seqs) > 0:
+ raise RuntimeError("Cannot compute DockQ for reference structures "
+ "with nucleotide chains if dockq_capri_peptide "
+ "is enabled.")
+
# lists with values in contact_target_interfaces
self._dockq_scores = list()
self._fnat = list()
@@ -2019,6 +2605,37 @@ class Scorer:
if ch.GetName() not in processed_trg_chains:
self._n_target_not_mapped += len(ch.residues)
+ def _extract_mapped_pos_full_bb(self):
+ self._mapped_target_pos_full_bb = geom.Vec3List()
+ self._mapped_model_pos_full_bb = geom.Vec3List()
+ exp_pep_atoms = ["N", "CA", "C"]
+ exp_nuc_atoms = ["\"O5'\"", "\"C5'\"", "\"C4'\"", "\"C3'\"", "\"O3'\""]
+ trg_pep_chains = [s.GetName() for s in self.chain_mapper.polypep_seqs]
+ trg_nuc_chains = [s.GetName() for s in self.chain_mapper.polynuc_seqs]
+ for trg_ch, mdl_ch in self.mapping.GetFlatMapping().items():
+ aln = self.mapping.alns[(trg_ch, mdl_ch)]
+ trg_ch = aln.GetSequence(0).GetName()
+ if trg_ch in trg_pep_chains:
+ exp_atoms = exp_pep_atoms
+ elif trg_ch in trg_nuc_chains:
+ exp_atoms = exp_nuc_atoms
+ else:
+ # this should be guaranteed by the chain mapper
+ raise RuntimeError("Unexpected error - contact OST developer")
+ for col in aln:
+ if col[0] != '-' and col[1] != '-':
+ trg_res = col.GetResidue(0)
+ mdl_res = col.GetResidue(1)
+ for aname in exp_atoms:
+ trg_at = trg_res.FindAtom(aname)
+ mdl_at = mdl_res.FindAtom(aname)
+ if not (trg_at.IsValid() and mdl_at.IsValid()):
+ # this should be guaranteed by the chain mapper
+ raise RuntimeError("Unexpected error - contact OST "
+ "developer")
+ self._mapped_target_pos_full_bb.append(trg_at.GetPos())
+ self._mapped_model_pos_full_bb.append(mdl_at.GetPos())
+
def _extract_rigid_mapped_pos(self):
self._rigid_mapped_target_pos = geom.Vec3List()
@@ -2047,6 +2664,36 @@ class Scorer:
if ch.GetName() not in processed_trg_chains:
self._rigid_n_target_not_mapped += len(ch.residues)
+ def _extract_rigid_mapped_pos_full_bb(self):
+ self._rigid_mapped_target_pos_full_bb = geom.Vec3List()
+ self._rigid_mapped_model_pos_full_bb = geom.Vec3List()
+ exp_pep_atoms = ["N", "CA", "C"]
+ exp_nuc_atoms = ["\"O5'\"", "\"C5'\"", "\"C4'\"", "\"C3'\"", "\"O3'\""]
+ trg_pep_chains = [s.GetName() for s in self.chain_mapper.polypep_seqs]
+ trg_nuc_chains = [s.GetName() for s in self.chain_mapper.polynuc_seqs]
+ for trg_ch, mdl_ch in self.rigid_mapping.GetFlatMapping().items():
+ aln = self.mapping.alns[(trg_ch, mdl_ch)]
+ trg_ch = aln.GetSequence(0).GetName()
+ if trg_ch in trg_pep_chains:
+ exp_atoms = exp_pep_atoms
+ elif trg_ch in trg_nuc_chains:
+ exp_atoms = exp_nuc_atoms
+ else:
+ # this should be guaranteed by the chain mapper
+ raise RuntimeError("Unexpected error - contact OST developer")
+ for col in aln:
+ if col[0] != '-' and col[1] != '-':
+ trg_res = col.GetResidue(0)
+ mdl_res = col.GetResidue(1)
+ for aname in exp_atoms:
+ trg_at = trg_res.FindAtom(aname)
+ mdl_at = mdl_res.FindAtom(aname)
+ if not (trg_at.IsValid() and mdl_at.IsValid()):
+ # this should be guaranteed by the chain mapper
+ raise RuntimeError("Unexpected error - contact OST "
+ "developer")
+ self._rigid_mapped_target_pos_full_bb.append(trg_at.GetPos())
+ self._rigid_mapped_model_pos_full_bb.append(mdl_at.GetPos())
def _compute_cad_score(self):
if not self.resnum_alignments:
@@ -2477,7 +3124,7 @@ class Scorer:
def _compute_tmscore(self):
res = None
if self.usalign_exec is None:
- LogScript("Computing patch TM-score with USalign exectuable")
+ LogScript("Computing TM-score with built-in USalign")
if self.oum:
flat_mapping = self.rigid_mapping.GetFlatMapping()
LogInfo("Overriding TM-score chain mapping")
@@ -2486,7 +3133,7 @@ class Scorer:
else:
res = bindings.WrappedMMAlign(self.model, self.target)
else:
- LogScript("Computing patch TM-score with built-in USalign")
+ LogScript("Computing TM-score with USalign exectuable")
if self.oum:
LogInfo("Overriding TM-score chain mapping")
flat_mapping = self.rigid_mapping.GetFlatMapping()
@@ -2502,5 +3149,29 @@ class Scorer:
for a,b in zip(res.ent1_mapped_chains, res.ent2_mapped_chains):
self._usalign_mapping[b] = a
+ def _resnum_connect(self, ent):
+ ed = None
+ for ch in ent.chains:
+ res_list = ch.residues
+ for i in range(len(res_list) - 1):
+ ra = res_list[i]
+ rb = res_list[i+1]
+ if ra.GetNumber().GetNum() + 1 == rb.GetNumber().GetNum():
+ if ra.IsPeptideLinking() and rb.IsPeptideLinking():
+ c = ra.FindAtom("C")
+ n = rb.FindAtom("N")
+ if c.IsValid() and n.IsValid() and not mol.BondExists(c, n):
+ if ed is None:
+ ed = ent.EditXCS(mol.BUFFERED_EDIT)
+ ed.Connect(c,n,1)
+ elif ra.IsNucleotideLinking() and rb.IsNucleotideLinking():
+ o = ra.FindAtom("O3'")
+ p = rb.FindAtom("P")
+ if o.IsValid() and p.IsValid()and not mol.BondExists(o, p):
+ if ed is None:
+ ed = ent.EditXCS(mol.BUFFERED_EDIT)
+ ed.Connect(o,p,1)
+
+
# specify public interface
__all__ = ('lDDTBSScorer', 'Scorer',)
diff --git a/modules/mol/alg/src/gdt.cc b/modules/mol/alg/src/gdt.cc
index 93b6e8b173b51b04fb4c0e01c092f661b7507d04..8eedd7796241085c89ec8ed15d8f28ca168afd03 100644
--- a/modules/mol/alg/src/gdt.cc
+++ b/modules/mol/alg/src/gdt.cc
@@ -302,8 +302,126 @@ void GDT(const geom::Vec3List& mdl_pos, const geom::Vec3List& ref_pos,
int n_pos = mdl_pos.size();
+ // deal with special cases that don't produce valid transforms
+ if(n_pos == 1) {
+ transform = geom::Mat4::Identity();
+ transform.PasteTranslation(ref_pos[0] - mdl_pos[0]);
+ n_superposed = 1;
+ return;
+ }
+
+ if(n_pos == 2) {
+ Real mdl_d = geom::Distance(mdl_pos[0], mdl_pos[1]);
+ Real ref_d = geom::Distance(ref_pos[0], ref_pos[1]);
+ Real dd = std::abs(mdl_d - ref_d);
+ if(dd/2 <= distance_thresh) {
+ // the two can be superposed within specified distance threshold
+ // BUT: cannot construct valid transformation from two positions
+ // => Construct matrix with four positions
+ // Two are constructed starting from the center point +- some direction
+ // vector that is orthogonal to the vector connecting the original two
+ // points.
+ geom::Vec3 mdl_center = (mdl_pos[0] + mdl_pos[1])*0.5;
+ geom::Vec3 ref_center = (ref_pos[0] + ref_pos[1])*0.5;
+ Eigen::Matrix<double, 4, 3> eigen_mdl_pos = \
+ Eigen::Matrix<double, 4, 3>::Zero(4, 3);
+ Eigen::Matrix<double, 4, 3> eigen_ref_pos = \
+ Eigen::Matrix<double, 4, 3>::Zero(4, 3);
+ Eigen::Matrix<double,3,3> eigen_rot = \
+ Eigen::Matrix<double,3,3>::Identity();
+
+ geom::Vec3 mdl_dir = geom::Normalize(mdl_pos[1] - mdl_pos[0]);
+ geom::Vec3 ref_dir = geom::Normalize(ref_pos[1] - ref_pos[0]);
+ geom::Vec3 mdl_normal;
+ geom::Vec3 ref_normal;
+
+ // Use cross product to get some normal on mdl_dir
+ // The direction of the second vector doesn't really matter, but shouldnt
+ // be collinear with mdl_dir
+ if(mdl_dir[0] < 0.999) {
+ mdl_normal = geom::Cross(geom::Vec3(1,0,0), mdl_dir);
+ } else {
+ mdl_normal = geom::Cross(geom::Vec3(0,1,0), mdl_dir);
+ }
+
+ // same for ref_dir
+ if(ref_dir[0] < 0.999) {
+ ref_normal = geom::Cross(geom::Vec3(1,0,0), ref_dir);
+ } else {
+ ref_normal = geom::Cross(geom::Vec3(0,1,0), ref_dir);
+ }
+
+ eigen_mdl_pos(0, 0) = mdl_pos[0][0] - mdl_center[0];
+ eigen_mdl_pos(0, 1) = mdl_pos[0][1] - mdl_center[1];
+ eigen_mdl_pos(0, 2) = mdl_pos[0][2] - mdl_center[2];
+ eigen_mdl_pos(1, 0) = mdl_pos[1][0] - mdl_center[0];
+ eigen_mdl_pos(1, 1) = mdl_pos[1][1] - mdl_center[1];
+ eigen_mdl_pos(1, 2) = mdl_pos[1][2] - mdl_center[2];
+ eigen_mdl_pos(2, 0) = mdl_normal[0];
+ eigen_mdl_pos(2, 1) = mdl_normal[1];
+ eigen_mdl_pos(2, 2) = mdl_normal[2];
+ eigen_mdl_pos(3, 0) = -mdl_normal[0];
+ eigen_mdl_pos(3, 1) = -mdl_normal[1];
+ eigen_mdl_pos(3, 2) = -mdl_normal[2];
+ eigen_ref_pos(0, 0) = ref_pos[0][0] - ref_center[0];
+ eigen_ref_pos(0, 1) = ref_pos[0][1] - ref_center[1];
+ eigen_ref_pos(0, 2) = ref_pos[0][2] - ref_center[2];
+ eigen_ref_pos(1, 0) = ref_pos[1][0] - ref_center[0];
+ eigen_ref_pos(1, 1) = ref_pos[1][1] - ref_center[1];
+ eigen_ref_pos(1, 2) = ref_pos[1][2] - ref_center[2];
+ eigen_ref_pos(2, 0) = ref_normal[0];
+ eigen_ref_pos(2, 1) = ref_normal[1];
+ eigen_ref_pos(2, 2) = ref_normal[2];
+ eigen_ref_pos(3, 0) = -ref_normal[0];
+ eigen_ref_pos(3, 1) = -ref_normal[1];
+ eigen_ref_pos(3, 2) = -ref_normal[2];
+
+ Real tmp; // no need to store RMSD
+ TheobaldRMSD(eigen_mdl_pos, eigen_ref_pos, tmp, eigen_rot);
+ transform = geom::Mat4();
+ transform(0, 0) = eigen_rot(0, 0);
+ transform(0, 1) = eigen_rot(0, 1);
+ transform(0, 2) = eigen_rot(0, 2);
+ transform(1, 0) = eigen_rot(1, 0);
+ transform(1, 1) = eigen_rot(1, 1);
+ transform(1, 2) = eigen_rot(1, 2);
+ transform(2, 0) = eigen_rot(2, 0);
+ transform(2, 1) = eigen_rot(2, 1);
+ transform(2, 2) = eigen_rot(2, 2);
+
+ // there are three transformation to be applied to reach ref_pos from
+ // mdl_pos:
+ // 1: shift mdl_pos to center
+ // 2: apply estimated rotation
+ // 3: shift onto average of ref_pos
+ Eigen::Matrix<double,1,3> eigen_avg_mdl = Eigen::Matrix<double,1,3>::Zero();
+ Eigen::Matrix<double,1,3> eigen_avg_ref = Eigen::Matrix<double,1,3>::Zero();
+ eigen_avg_mdl(0,0) = mdl_center[0];
+ eigen_avg_mdl(0,1) = mdl_center[1];
+ eigen_avg_mdl(0,2) = mdl_center[2];
+ eigen_avg_ref(0,0) = ref_center[0];
+ eigen_avg_ref(0,1) = ref_center[1];
+ eigen_avg_ref(0,2) = ref_center[2];
+ Eigen::Matrix<double,1,3> translation = eigen_rot *
+ (-eigen_avg_mdl.transpose()) +
+ eigen_avg_ref.transpose();
+ transform(0, 3) = translation(0, 0);
+ transform(1, 3) = translation(0, 1);
+ transform(2, 3) = translation(0, 2);
+ n_superposed = 2;
+ } else {
+ // the two cannot be superposed within specified distance threshold
+ // => just set n_superposed to 1 and generate the same transformation
+ // as in n_pos == 1
+ transform = geom::Mat4::Identity();
+ transform.PasteTranslation(ref_pos[0] - mdl_pos[0]);
+ n_superposed = 1;
+ }
+ return;
+ }
+
if(window_size > n_pos) {
- throw ost::Error("Window size in GDT algorithm is larger than positions");
+ window_size = n_pos;
}
Eigen::Matrix<double, Eigen::Dynamic, 3> eigen_mdl_pos = \
diff --git a/modules/mol/alg/tests/test_chain_mapping.py b/modules/mol/alg/tests/test_chain_mapping.py
index 6e17a67e600a4d9b5c6b29423608802e89cd8472..b71676b2681d6c2711a7d85789b1023bb1348683 100644
--- a/modules/mol/alg/tests/test_chain_mapping.py
+++ b/modules/mol/alg/tests/test_chain_mapping.py
@@ -133,13 +133,8 @@ class TestChainMapper(unittest.TestCase):
ed.SetResidueNumber(r, mol.ResNum(r.GetNumber().GetNum() + 42))
self.assertRaises(Exception, ChainMapper, tmp_ent, resnum_alignments=True)
- # chain B has a missing Valine... set pep_gap_thr to 0.0 should give an
- # additional chem group
- mapper = ChainMapper(ent, pep_gap_thr=0.0)
- self.assertEqual(len(mapper.chem_groups), 4)
-
# introduce single point mutation in ent... increasing pep_seqid_thr to 100
- # should give an additional chem group too
+ # should give an additional chem group
mapper = ChainMapper(ent, pep_seqid_thr=100.)
self.assertEqual(len(mapper.chem_groups), 3)
tmp_ent = ent.Copy()
@@ -166,9 +161,6 @@ class TestChainMapper(unittest.TestCase):
mapper = ChainMapper(tmp_ent, nuc_seqid_thr=100.)
self.assertEqual(len(mapper.polynuc_seqs), 3)
self.assertEqual(len(mapper.chem_groups), 4)
- mapper = ChainMapper(tmp_ent, nuc_gap_thr=0.0)
- self.assertEqual(len(mapper.polynuc_seqs), 3)
- self.assertEqual(len(mapper.chem_groups), 4)
def test_chem_mapping(self):
diff --git a/modules/mol/alg/tests/test_ligand_scoring.py b/modules/mol/alg/tests/test_ligand_scoring.py
index 64f0b580183c79332f958f7ed049cbd3cda2c3e8..90cb8f81524d7c259bf2d7cbf1891745b8601b30 100644
--- a/modules/mol/alg/tests/test_ligand_scoring.py
+++ b/modules/mol/alg/tests/test_ligand_scoring.py
@@ -646,7 +646,7 @@ class TestLigandScoringFancy(unittest.TestCase):
trg_report, trg_pair_report = sc.get_target_ligand_state_report(0)
exp_lig_report = {'state': 10.0,
- 'short desc': 'binding_site',
+ 'short desc': 'target_binding_site',
'desc': 'No residues were in proximity of the target ligand.'}
exp_pair_report = [{'state': 1, 'short desc': 'identity',
@@ -856,18 +856,18 @@ class TestLigandScoringFancy(unittest.TestCase):
self.assertDictEqual(sc.unassigned_model_ligands_reasons, {
- 'L_ZN': {1: 'model_representation'},
- 'L_NA': {1: 'binding_site'},
+ 'L_ZN': {1: 'model_binding_site'},
+ 'L_NA': {1: 'target_binding_site'},
'L_OXY': {1: 'identity'},
'L_MG_2': {1: 'stoichiometry'},
"L_CMO": {1: 'disconnected'}
})
self.assertDictEqual(sc.unassigned_target_ligands_reasons, {
'G': {1: 'identity'},
- 'H': {1: 'model_representation'},
+ 'H': {1: 'model_binding_site'},
'J': {1: 'stoichiometry'},
'K': {1: 'identity'},
- 'L_NA': {1: 'binding_site'},
+ 'L_NA': {1: 'target_binding_site'},
"L_CMO": {1: 'disconnected'}
})
self.assertTrue("L_OXY" not in sc.score)
@@ -891,23 +891,23 @@ class TestLigandScoringFancy(unittest.TestCase):
# Test with partial bs search (full_bs_search=True)
# Here we expect L_MG_2 to be unassigned because of stoichiometry
- # rather than model_representation, as it no longer matters so far from
+ # rather than model_binding_site, as it no longer matters so far from
# the binding site.
sc = ligand_scoring_scrmsd.SCRMSDScorer(mdl, trg, None, None,
full_bs_search=True)
self.assertEqual(sc.unassigned_model_ligands_reasons, {
- 'L_ZN': {1: 'model_representation'},
- 'L_NA': {1: 'binding_site'},
+ 'L_ZN': {1: 'model_binding_site'},
+ 'L_NA': {1: 'target_binding_site'},
'L_OXY': {1: 'identity'},
'L_MG_2': {1: 'stoichiometry'},
"L_CMO": {1: 'disconnected'}
})
self.assertEqual(sc.unassigned_target_ligands_reasons, {
'G': {1: 'identity'},
- 'H': {1: 'model_representation'},
+ 'H': {1: 'model_binding_site'},
'J': {1: 'stoichiometry'},
'K': {1: 'identity'},
- 'L_NA': {1: 'binding_site'},
+ 'L_NA': {1: 'target_binding_site'},
"L_CMO": {1: 'disconnected'}
})
diff --git a/modules/mol/base/doc/editors.rst b/modules/mol/base/doc/editors.rst
index 01214c4b3a3798d3fedf10a345b04aee56d42c19..848781a7825d67d770c63e132ed2248e75beb540 100644
--- a/modules/mol/base/doc/editors.rst
+++ b/modules/mol/base/doc/editors.rst
@@ -11,16 +11,30 @@ atoms. There are two flavors of editors, one for the internal coordinate system
Edit Modes
--------------------------------------------------------------------------------
-Editors support two modes: An unbuffered edit mode and a buffered edit mode. In
-the unbuffered edit mode, dependent information such as the spatial organizer
-and the internal coordinate system (in case of the XCSEditor) are updated after
-every change. In buffered edit mode, the updates are delayed until one of the
-following happens:
+.. class:: EditMode
- * The last editor goes out of scope.
- * :meth:`XCSEditor.UpdateICS` or :meth:`ICSEditor.UpdateXCS` is called
- explicitly.
+ Editors support two modes: An unbuffered edit mode and a buffered edit mode.
+ .. attribute:: UNBUFFERED_EDIT
+
+ This is the default edit mode in editors. In unbuffered edit mode,
+ dependent information such as the spatial organizer and the internal
+ coordinate system (in case of the :class:`XCSEditor`) are updated after
+ every change.
+
+ .. warning::
+
+ The unbuffered edit mode can be extremely inefficient as every operation
+ can trigger expensive calculations.
+
+ .. attribute:: BUFFERED_EDIT
+
+ In buffered edit mode, the updates are delayed until one of the following
+ happens:
+
+ * The last editor goes out of scope.
+ * :meth:`XCSEditor.UpdateICS` or :meth:`ICSEditor.UpdateXCS` is called
+ explicitly.
The editors follow the RIAA (resource allocation is initialization) principle:
Whenever an editor is requested an internal reference counter is incremented. In
@@ -28,12 +42,41 @@ the destructor, this reference count is decremented. When the count drops to
zero, the dependent information is updated.
In Python, one can not rely on the destructors being called. It is advised to
-always put a call to :meth:`XCSEditor.UpdateICS` or
-:meth:`ICSEditor.UpdateXCS` when the editing is finished. Alternatively,
-starting from Python version 2.6, one can use the \
-`with <http://docs.python.org/reference/compound_stmts.html#with>`_
-statement to make sure the destructor are called and the dependent information
-is updated.
+always put a call to :meth:`XCSEditor.UpdateICS` or
+:meth:`ICSEditor.UpdateXCS` when the editing is finished.
+
+Here is a short example of how to use the editors to correctly create a
+methionine residue:
+
+.. code-block:: python
+
+ entity = ost.mol.CreateEntity()
+
+ # A buffered external coordinate system editor covers most use cases
+ editor = entity.EditXCS(ost.mol.EditMode.BUFFERED_EDIT)
+
+ # Insert a new protein chain
+ chain = editor.InsertChain("A");
+ editor.SetChainType(chain, ost.mol.ChainType.CHAINTYPE_POLY_PEPTIDE_L)
+
+ # Insert a residue
+ res1 = editor.AppendResidue(chain, "MET");
+ res1.SetChemClass(ost.mol.L_PEPTIDE_LINKING)
+ res1.SetChemType(ost.mol.ChemType.AMINOACIDS)
+
+ # Insert atoms
+ editor.InsertAtom(res1, "N", ost.geom.Vec3(21.609, 35.384, 56.705))
+ editor.InsertAtom(res1, "CA", ost.geom.Vec3(20.601, 35.494, 57.793))
+ editor.InsertAtom(res1, "C", ost.geom.Vec3(19.654, 34.300, 57.789))
+ editor.InsertAtom(res1, "O", ost.geom.Vec3(18.447, 34.456, 57.595))
+
+ # Connect them
+ editor.Connect(res1.FindAtom("N"), res1.FindAtom("CA"))
+ editor.Connect(res1.FindAtom("CA"), res1.FindAtom("C"))
+ editor.Connect(res1.FindAtom("C"), res1.FindAtom("O"))
+
+ # Finalize
+ editor.UpdateICS()
Basic Editing Operations
@@ -44,8 +87,8 @@ The basic functionality of editors is implemented in the EditorBase class.
.. note::
To use the editing functions available in :class:`EditorBase`, it is
- recommended to use the external coordinate system :class:`XCSEditor` for
- performance reasons.
+ recommended to use the external coordinate system :class:`XCSEditor` with
+ buffering for performance reasons.
.. class:: EditorBase
@@ -161,7 +204,7 @@ The basic functionality of editors is implemented in the EditorBase class.
:type is_hetatm: bool
:returns: :class:`AtomHandle`
- .. method:: InsertAltAtom(residue, atom_name, alt_group, pos, element="", occupancy=1.0, b_factor=0.0)
+ .. method:: InsertAltAtom(residue, atom_name, alt_group, pos, ele="", occ=1.0, b_factor=0.0)
Insert new atom with alternative position indicator
@@ -178,14 +221,14 @@ The basic functionality of editors is implemented in the EditorBase class.
:type alt_group: string
:param pos: is the position of the atom in global coordinates
:type pos: :class:`~ost.geom.Vec3`
- :param element: is the atom's element. If set to a a valid element,
+ :param ele: is the atom's element. If set to a a valid element,
atom properties such as mass, charge, radius are set
based on default values for that element. If the element
string is empty, or unknown, the properties are filled
with rather meaningless default values.
- :type element: class:`string`
- :param occupancy: The occupancy of the atom. between 0 and 1
- :type occupancy: float
+ :type ele: class:`string`
+ :param occ: The occupancy of the atom. between 0 and 1
+ :type occ: float
:param b_factor: temperature factor.
:type b_factor: float
diff --git a/modules/mol/base/doc/entity.rst b/modules/mol/base/doc/entity.rst
index 6d8f4ea104e66dc0996b1e23cda9188b54bedb31..8afb1a2d3001cfe56322ab8714c173e8d7fb59b1 100644
--- a/modules/mol/base/doc/entity.rst
+++ b/modules/mol/base/doc/entity.rst
@@ -932,6 +932,16 @@ Residue Handle
See :attr:`chain`
+ .. method:: GetChemClass()
+ SetChemClass()
+
+ See :attr:`chem_class`
+
+ .. method:: GetChemType()
+ SetChemType()
+
+ See :attr:`chem_type`
+
.. method:: GetPhiTorsion()
See :attr:`phi_torsion`
@@ -944,14 +954,6 @@ Residue Handle
See :attr:`omega_torsion`
- .. method:: GetChemClass()
-
- See :attr:`chem_class`
-
- .. method:: GetChemType()
-
- See :attr:`chem_type`
-
.. method:: GetSecStructure()
See :attr:`sec_structure`
@@ -2175,7 +2177,9 @@ Residue View
GetPsiTorsion
GetOmegaTorsion
GetChemClass
+ SetChemClass
GetChemType
+ SetChemType
GetSecStructure
IsLigand
SetIsLigand
diff --git a/modules/mol/base/pymod/export_editors.cc b/modules/mol/base/pymod/export_editors.cc
index e2cd89b9b7730a5e7245ed93a2818a57a55e9d06..70872ac8bbed7faa96bd5eac20bc4805facab59c 100644
--- a/modules/mol/base/pymod/export_editors.cc
+++ b/modules/mol/base/pymod/export_editors.cc
@@ -120,6 +120,7 @@ void export_Editors()
.def("InsertAtom", insert_atom_b)
.def("InsertAltAtom", insert_alt_atom_a)
.def("InsertAltAtom", insert_alt_atom_b)
+ .def("AddAltAtomPos", &EditorBase::AddAltAtomPos)
.def("DeleteResidue", &EditorBase::DeleteResidue)
.def("DeleteChain", &EditorBase::DeleteChain)
.def("DeleteAtom", &EditorBase::DeleteAtom)
diff --git a/modules/mol/base/pymod/export_residue.cc b/modules/mol/base/pymod/export_residue.cc
index 7a0883e7fc30d8a938e821a2bf7fbe93a3a762b8..76fd04be7e0ed6e5fa85abdbd3f6cde91d25f2a4 100644
--- a/modules/mol/base/pymod/export_residue.cc
+++ b/modules/mol/base/pymod/export_residue.cc
@@ -225,7 +225,8 @@ void export_Residue()
.add_property("chem_class", &ResidueBase::GetChemClass, set_chemclass)
.def("SetChemClass", set_chemclass)
.def("GetChemType", &ResidueBase::GetChemType)
- .add_property("chem_type", &ResidueBase::GetChemType)
+ .add_property("chem_type", &ResidueBase::GetChemType, &ResidueBase::SetChemType)
+ .def("SetChemType", &ResidueBase::SetChemType)
.add_property("is_ligand", &ResidueBase::IsLigand, &ResidueBase::SetIsLigand)
.def("IsLigand", &ResidueBase::IsLigand)
.def("SetIsLigand", &ResidueBase::SetIsLigand, arg("ligand"))
diff --git a/modules/seq/alg/src/wrap_parasail.cc b/modules/seq/alg/src/wrap_parasail.cc
index bede42a34698d6c8bd72281700e6e5afd265aa04..8938b5ee8e3c5919c003f1f42699eabd03d5d9ee 100644
--- a/modules/seq/alg/src/wrap_parasail.cc
+++ b/modules/seq/alg/src/wrap_parasail.cc
@@ -17,6 +17,7 @@
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
//------------------------------------------------------------------------------
+#include <ost/message.hh>
#include <ost/seq/alg/wrap_parasail.hh>
#if OST_PARASAIL_ENABLED
@@ -131,6 +132,45 @@ ost::seq::AlignmentList Align(const ost::seq::ConstSequenceHandle& s1,
return aln_list;
}
+ost::seq::AlignmentList AlignEmptySeq(const ost::seq::ConstSequenceHandle& s1,
+ const ost::seq::ConstSequenceHandle& s2) {
+ ost::seq::AlignmentHandle aln = ost::seq::CreateAlignment();
+ if(s1.GetLength() == 0 && s2.GetLength() == 0) {
+ ost::seq::SequenceHandle new_s1 = ost::seq::CreateSequence(s1.GetName(), s1.GetString());
+ new_s1.SetOffset(s1.GetOffset());
+ aln.AddSequence(new_s1);
+ ost::seq::SequenceHandle new_s2 = ost::seq::CreateSequence(s2.GetName(), s2.GetString());
+ new_s2.SetOffset(s2.GetOffset());
+ aln.AddSequence(new_s2);
+ } else if(s1.GetLength() == 0) {
+ String x(s2.GetLength(), '-');
+ ost::seq::SequenceHandle new_s1 = ost::seq::CreateSequence(s1.GetName(), x);
+ new_s1.SetOffset(s1.GetOffset());
+ aln.AddSequence(new_s1);
+ ost::seq::SequenceHandle new_s2 = ost::seq::CreateSequence(s2.GetName(), s2.GetString());
+ new_s2.SetOffset(s2.GetOffset());
+ aln.AddSequence(new_s2);
+ } else if(s2.GetLength() == 0) {
+ ost::seq::SequenceHandle new_s1 = ost::seq::CreateSequence(s1.GetName(), s1.GetString());
+ new_s1.SetOffset(s1.GetOffset());
+ aln.AddSequence(new_s1);
+ String x(s1.GetLength(), '-');
+ ost::seq::SequenceHandle new_s2 = ost::seq::CreateSequence(s2.GetName(), x);
+ new_s2.SetOffset(s2.GetOffset());
+ aln.AddSequence(new_s2);
+ } else {
+ throw ost::Error("One seq must be empty in AlignEmptySeq");
+ }
+ if(s1.HasAttachedView()) {
+ aln.AttachView(0, s1.GetAttachedView());
+ }
+ if(s2.HasAttachedView()) {
+ aln.AttachView(1, s2.GetAttachedView());
+ }
+ ost::seq::AlignmentList aln_list = {aln};
+ return aln_list;
+}
+
}
namespace ost{ namespace seq{ namespace alg{
@@ -139,6 +179,9 @@ ost::seq::AlignmentList ParaLocalAlign(const ost::seq::ConstSequenceHandle& s1,
const ost::seq::ConstSequenceHandle& s2,
ost::seq::alg::SubstWeightMatrixPtr& subst,
int gap_open, int gap_ext) {
+ if(s1.GetLength() == 0 || s2.GetLength() == 0) {
+ return ost::seq::AlignmentList();
+ }
return Align(s1, s2, gap_open, gap_ext, subst,
parasail_sw_trace_scan_sat, true);
}
@@ -147,6 +190,9 @@ ost::seq::AlignmentList ParaGlobalAlign(const ost::seq::ConstSequenceHandle& s1,
const ost::seq::ConstSequenceHandle& s2,
ost::seq::alg::SubstWeightMatrixPtr& subst,
int gap_open, int gap_ext) {
+ if(s1.GetLength() == 0 || s2.GetLength() == 0) {
+ return AlignEmptySeq(s1, s2);
+ }
return Align(s1, s2, gap_open, gap_ext, subst,
parasail_nw_trace_scan_sat, false);
}
@@ -155,6 +201,9 @@ ost::seq::AlignmentList ParaSemiGlobalAlign(const ost::seq::ConstSequenceHandle&
const ost::seq::ConstSequenceHandle& s2,
ost::seq::alg::SubstWeightMatrixPtr& subst,
int gap_open, int gap_ext) {
+ if(s1.GetLength() == 0 || s2.GetLength() == 0) {
+ return AlignEmptySeq(s1, s2);
+ }
return Align(s1, s2, gap_open, gap_ext, subst,
parasail_sg_trace_scan_sat, false);
}
diff --git a/singularity/Singularity b/singularity/Singularity
index 395c512ca3289341613932f39e24ae9aed0fc264..8ed1bb886168a62cb735a7528204265376e8b580 100644
--- a/singularity/Singularity
+++ b/singularity/Singularity
@@ -1,5 +1,5 @@
BootStrap: docker
-From: registry.scicore.unibas.ch/schwede/openstructure:2.8.0
+From: registry.scicore.unibas.ch/schwede/openstructure:2.9.0
%post
##############################################################################
# POST