Skip to content

Fix crash from altlocs changing between incremental alignments #72

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mwinokan wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Fix crash from altlocs changing between incremental alignments #72

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e6b87fd
Support altlocs changing between alignments
mwinokan Apr 20, 2026
a9fbf6f
remove changes from black formatter
mwinokan Apr 20, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
125 changes: 74 additions & 51 deletions src/xchemalign/aligner.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ def try_make(path):


def read_yaml(path):
with open(path, 'r') as f:
with open(path, "r") as f:
dic = yaml.safe_load(f)

return dic
Expand Down Expand Up @@ -164,7 +164,7 @@ def __init__(self, dir, log_file=None, log_level=0, debug=False):
self._find_version_dir(dir) # sets self.working_dir and self.version_dir
output_path = self.working_dir / "upload-current"
if not log_file:
log_file = output_path / 'aligner.log'
log_file = output_path / "aligner.log"
self.logger = utils.Logger(logfile=log_file, level=log_level)
self.logger.info("Using", self.version_dir, "as current version dir")
self.base_dir = self.version_dir.parent # e.g. path/to
Expand Down Expand Up @@ -199,18 +199,18 @@ def _find_version_dir(self, dir):
exit(1)

self.working_dir = wd1
current_dir = self.working_dir / 'upload-current'
current_dir = self.working_dir / "upload-current"

# check that we at least have an upload_1 dir
if not current_dir.joinpath('upload_1').is_dir():
if not current_dir.joinpath("upload_1").is_dir():
self._log_error("Working dir " + str(dir) + " does not contain and upload_? dirs")
exit(1)

# now find the latest upload_? dir
i = 0
while i < 100:
i += 1
version_dir = current_dir.joinpath(f'upload_{str(i)}')
version_dir = current_dir.joinpath(f"upload_{str(i)}")
if version_dir.is_dir():
self.version_dir = version_dir
else:
Expand Down Expand Up @@ -290,7 +290,7 @@ def _write_output(self, collator_dict, aligner_dict):
)

else:
self.logger.warn('crystal {} not found in input. This is very strange.'.format(k))
self.logger.warn("crystal {} not found in input. This is very strange.".format(k))

collator_dict[Constants.META_REFERENCE_ALIGNMENTS] = aligner_dict[Constants.META_REFERENCE_ALIGNMENTS]
traverse_dictionary(
Expand Down Expand Up @@ -341,7 +341,7 @@ def _write_output(self, collator_dict, aligner_dict):
del aligner_dict[Constants.META_XTALS][_crystal_name][_aligned_files][_chain][_res][_alt]

# remove this eventually
with open(self.version_dir / 'aligner_tmp.yaml', "w") as stream:
with open(self.version_dir / "aligner_tmp.yaml", "w") as stream:
yaml.dump(aligner_dict, stream, sort_keys=False, default_flow_style=None)

with open(self.version_dir / Constants.METADATA_ALIGN_FILENAME, "w") as stream:
Expand Down Expand Up @@ -405,8 +405,8 @@ def _perform_alignments(self, meta):
self.logger.info(f"Got {len(new_datasets)} new datasets")
if sum([len(dataset.ligand_binding_events) for dataset in datasets.values()]) == 0:
self.logger.error(
'There are no ligand binding events detected! The program will now exit.\n'
'This is most likely because because there are no ligand cif files in the input data.'
"There are no ligand binding events detected! The program will now exit.\n"
"This is most likely because because there are no ligand cif files in the input data."
)
raise Exception

Expand Down Expand Up @@ -461,10 +461,10 @@ def _perform_alignments(self, meta):

if source_fs_model:
self.logger.info(f"Have source fs model at {source_fs_model.ligand_neighbourhood_transforms}!")
ligand_neighbourhood_transforms: dict[
tuple[tuple[str, str, str], tuple[str, str, str]], dt.Transform
] = _load_ligand_neighbourhood_transforms(
source_fs_model.ligand_neighbourhood_transforms, fail_if_not_found=True
ligand_neighbourhood_transforms: dict[tuple[tuple[str, str, str], tuple[str, str, str]], dt.Transform] = (
_load_ligand_neighbourhood_transforms(
source_fs_model.ligand_neighbourhood_transforms, fail_if_not_found=True
)
)
else:
ligand_neighbourhood_transforms = _load_ligand_neighbourhood_transforms(
Expand Down Expand Up @@ -618,7 +618,7 @@ def _perform_alignments(self, meta):

new_meta[Constants.META_XTALS] = {}
for dtag, crystal in crystals.items():
self.logger.info('looking at', dtag)
self.logger.info("looking at", dtag)

new_meta[Constants.META_XTALS][dtag] = {}
crystal_output = new_meta[Constants.META_XTALS][dtag]
Expand All @@ -630,51 +630,75 @@ def _perform_alignments(self, meta):

# Skip if no output for this dataset
if dtag not in updated_fs_model.alignments:
self.logger.warn('skipping {} as aligned structures not found'.format(dtag))
self.logger.warn("skipping {} as aligned structures not found".format(dtag))
continue

# We capture the ligand binding events info as we need to know whether the event map file is present
# This is a bit of a hack as the event map file location is generated by LNA even if there is no event map
# so we need to know whether to actually include it in the metadata.
# It would be better if LNA only included if it actually existed which would make the checking easier.
event_map_dict_list = crystal.get(Constants.META_XTAL_FILES, {}).get(Constants.META_BINDING_EVENT, {})
event_map_dict_list = crystal.get(Constants.META_XTAL_FILES, {}).get(Constants.META_BINDING_EVENT, [])

crystal_output[Constants.META_ALIGNED_FILES] = {}
aligned_output = crystal_output[Constants.META_ALIGNED_FILES]
dataset_output = updated_fs_model.alignments[dtag]
# print(crystal)
# print(dataset_output)
# print(event_map_dict_list)

# Build a lookup from (chain, res, altloc_string) -> event map dict entry so that
# changes in the number of altconfs between versions don't cause a crash.
event_map_lookup = {}
for entry in event_map_dict_list:
key = (
entry.get(Constants.META_PROT_CHAIN),
str(entry.get(Constants.META_PROT_RES)),
str(entry.get(Constants.META_PROT_ALTLOC)),
)
event_map_lookup[key] = entry

for chain_name, chain_output in dataset_output.items():
aligned_chain_output = aligned_output[chain_name] = {}
i = 0
for ligand_residue, ligand_output in chain_output.items():
aligned_ligand_output = aligned_chain_output[ligand_residue] = {}
for altoloc, altloc_output in ligand_output.items():
aligned_altloc_output = aligned_ligand_output[altoloc] = {}
altloc_str = dt.altloc_to_string(altoloc)
lookup_key = (chain_name, ligand_residue, altloc_str)
event_map_entry = event_map_lookup.get(lookup_key)
if event_map_entry is None:
# Try to give a more specific diagnosis: check if the same chain+altloc
# exists under a different residue number (residue renumbering case) vs
# a genuinely new altloc that has no prior event map (altconf added case).
res_changed = any(
k[0] == chain_name and k[2] == altloc_str and k[1] != ligand_residue
for k in event_map_lookup
)
if res_changed:
self.logger.warn(
"No event map metadata found for ligand "
+ ligand_residue
+ " altloc "
+ altloc_str
+ " in crystal "
+ dtag
+ ". The residue number appears to have changed between versions."
+ " Proceeding without event map for this ligand."
)
else:
self.logger.warn(
"No event map metadata found for ligand "
+ ligand_residue
+ " altloc "
+ altloc_str
+ " in crystal "
+ dtag
+ ". This can happen when altconfs are added between versions."
+ " Proceeding without event map for this altloc."
)
for version, version_output in altloc_output.items():
aligned_version_output = aligned_altloc_output[version] = {}
for site_id, aligned_structure_path in version_output.aligned_structures.items():
# Is the event map file present?
# We do this assuming the order is the same as the info in the crystallographic files section
# But first do a check that event_map_dict_list contains the expected entry, as when the
# ligand residue number changes between versions then you get an error, so we try to
# handle this nicely. But there could be other things that cause this error.
if i >= len(event_map_dict_list):
msg = (
'Unexpected number of event maps found for ligand '
+ ligand_residue
+ ' in crystal '
+ dtag
+ '. Possible causes are the ligand residue number'
+ ' or the number of ligands present changing between versions;'
+ ' if it\'s not that, the cause might be even more obscure.'
+ ' Talk to a developer, and good luck.'
+ ' If you can\'t correct this then you can add this crystal to the \`exclude\` list.'
)
self._log_error(msg)
exit(1)
event_map_present = True if Constants.META_FILE in event_map_dict_list[i] else False
event_map_present = (
event_map_entry is not None and Constants.META_FILE in event_map_entry
)

aligned_artefacts_path = version_output.aligned_artefacts[site_id]
aligned_event_map_path = version_output.aligned_event_maps[site_id]
Expand Down Expand Up @@ -707,7 +731,6 @@ def _perform_alignments(self, meta):
aligned_version_output[site_id][
Constants.META_AIGNED_CRYSTALLOGRAPHIC_EVENT_MAP
] = aligned_crystallographic_event_map_path
i += 1

## Add the reference alignments
new_meta[Constants.META_REFERENCE_ALIGNMENTS] = {}
Expand Down Expand Up @@ -753,7 +776,7 @@ def _perform_alignments(self, meta):
if len(list(d.iterdir())) == 0:
empty_dir_count += 1
d.rmdir()
self.logger.info('removing {} empty aligned_files dirs'.format(empty_dir_count))
self.logger.info("removing {} empty aligned_files dirs".format(empty_dir_count))
return new_meta

def _extract_components(self, crystals, aligner_meta):
Expand All @@ -769,7 +792,7 @@ def _extract_components(self, crystals, aligner_meta):
:return:
"""

self.logger.info('extracting components')
self.logger.info("extracting components")

is_covalent = self.config.get(Constants.CONFIG_COVALENT, False)
self.logger.info("Covalent =", is_covalent)
Expand All @@ -778,7 +801,7 @@ def _extract_components(self, crystals, aligner_meta):
num_pdbs = 0
for k1, v1 in aligner_meta.get(Constants.META_XTALS, {}).items(): # k = xtal
if Constants.META_ALIGNED_FILES in v1:
self.logger.info('handling', k1)
self.logger.info("handling", k1)
cif_file = (
crystals.get(k1)
.get(Constants.META_XTAL_FILES, {})
Expand Down Expand Up @@ -825,19 +848,19 @@ def _extract_components(self, crystals, aligner_meta):
)
v6[Constants.META_LIGAND_MOL] = (
str(pdbxtal.ligand_base_file.relative_to(self.base_dir))
+ '.mol'
+ ".mol"
)
v6[Constants.META_LIGAND_SDF] = (
str(pdbxtal.ligand_base_file.relative_to(self.base_dir))
+ '.sdf'
+ ".sdf"
)
v6[Constants.META_LIGAND_PDB] = (
str(pdbxtal.ligand_base_file.relative_to(self.base_dir))
+ '.pdb'
+ ".pdb"
)
v6[Constants.META_LIGAND_SMILES] = (
str(pdbxtal.ligand_base_file.relative_to(self.base_dir))
+ '.smi'
+ ".smi"
)
v6[Constants.META_LIGAND_NAME] = pdbxtal.ligand_name
v6[Constants.META_LIGAND_SMILES_STRING] = pdbxtal.smiles
Expand Down Expand Up @@ -867,9 +890,9 @@ def main():
args = parser.parse_args()

if args.dir:
log = str(Path(args.dir).joinpath('aligner.log'))
log = str(Path(args.dir).joinpath("aligner.log"))
else:
log = 'aligner.log'
log = "aligner.log"

a = Aligner(args.dir, log_file=log, log_level=args.log_level)
logger = a.logger
Expand All @@ -891,7 +914,7 @@ def main():
logger.report()
logger.close()
if logger.logfilename:
to_path = a.version_dir / 'aligner.log'
to_path = a.version_dir / "aligner.log"
print("copying log file", logger.logfilename, "to", to_path)
f = shutil.copy2(logger.logfilename, to_path)
if not f:
Expand Down