diff --git a/src/pdbdepo/pdb_deposition.py b/src/pdbdepo/pdb_deposition.py index e0e7c3b..8fcc2a5 100644 --- a/src/pdbdepo/pdb_deposition.py +++ b/src/pdbdepo/pdb_deposition.py @@ -92,6 +92,70 @@ def read_software_templates(): return d +def substitute_tokens( + template: str, + xtal_name: str, + cmpd_code: str, + ext_codes: list, + pose_id_str: str, +) -> str: + """Replace all known template tokens in a string and return the result. + + $ExternalCodeN (N=2..9) are filled from ext_codes[N-2]; any that remain after substitution + are erased. The numbering starts at 2 because column 1 of the external-codes CSV is the + crystal name and is not exposed as a token. + """ + template = template.replace('$CompoundCode', cmpd_code) + template = template.replace('$CrystalName', xtal_name) + template = template.replace('$PoseID', pose_id_str) + if ext_codes: + for i, code in enumerate(ext_codes): + template = template.replace('$ExternalCode' + str(i + 2), code) + for i in range(1, 9): + template = template.replace('$ExternalCode' + str(i), '') + return template + + +def merge_mmcifgen_into_structure( + structure_block: cif.Block, + mmcifgen_block, + xtal_name: str, + cmpd_code: str, + cmpd_codes_dict: dict, + pose_ids_dict: dict, +) -> None: + """Apply the per-investigation template block from mmcif-gen onto a per-crystal structure block. + + Items that contain template tokens are expanded with crystal-specific values; all other items + are copied across unchanged. + """ + pose_id_str = ', '.join(pose_ids_dict.get(xtal_name, [])) + ext_codes = cmpd_codes_dict.get(xtal_name) + + for item in mmcifgen_block: + to_add = True + if item.loop is not None: + if '_struct.title' in item.loop.tags: + to_add = False + if len(item.loop.values) > 1: + expanded = substitute_tokens(item.loop.values[1], xtal_name, cmpd_code, ext_codes, pose_id_str) + new_loop = structure_block.init_loop('', item.loop.tags) + new_loop.add_row([item.loop.values[0], expanded]) + elif any('$PoseID' in v for v in item.loop.values): + to_add = False + ncols = len(item.loop.tags) + nrows = len(item.loop.values) // ncols + new_loop = structure_block.init_loop('', item.loop.tags) + for r in range(nrows): + loop_row = [item.loop.values[r * ncols + c].replace('$PoseID', pose_id_str) for c in range(ncols)] + new_loop.add_row(loop_row) + elif item.pair is not None and '$PoseID' in item.pair[1]: + structure_block.set_pair(item.pair[0], item.pair[1].replace('$PoseID', pose_id_str)) + to_add = False + if to_add: + structure_block.add_item(item) + + def process_input( base_dir: Path, input_path: Path, @@ -278,47 +342,10 @@ def process_input( mmcifgen_refine_values, ) - # build the comma-separated pose ID string for this crystal - pose_ids = pose_ids_dict.get(xtal_name, []) - pose_id_str = ', '.join(pose_ids) - # add in the common metadata (generated by mmcif-gen) - for item in mmcifgen_block: - to_add = True - if item.loop is not None: - if '_struct.title' in item.loop.tags: - # special case of the structure title loop that needs expanding - to_add = False - if len(item.loop.values) > 1: - template = item.loop.values[1] - template = template.replace('$CompoundCode', cmpd_code) - template = template.replace('$CrystalName', xtal_name) - template = template.replace('$PoseID', pose_id_str) - cmpd_codes = cmpd_codes_dict.get(xtal_name) - if cmpd_codes: - for i, cmpd_code in enumerate(cmpd_codes): - template = template.replace('$ExternalCode' + str(i + 2), cmpd_code) - # erase any non-substituted ExternalCodes - for i in range(1, 9): - template = template.replace('$ExternalCode' + str(i), '') - new_loop = structure_cif_block0.init_loop('', item.loop.tags) - new_loop.add_row([item.loop.values[0], template]) - elif any('$PoseID' in v for v in item.loop.values): - to_add = False - ncols = len(item.loop.tags) - nrows = len(item.loop.values) // ncols - new_loop = structure_cif_block0.init_loop('', item.loop.tags) - for r in range(nrows): - loop_row = [ - item.loop.values[r * ncols + c].replace('$PoseID', pose_id_str) for c in range(ncols) - ] - new_loop.add_row(loop_row) - elif item.pair is not None and '$PoseID' in item.pair[1]: - structure_cif_block0.set_pair(item.pair[0], item.pair[1].replace('$PoseID', pose_id_str)) - to_add = False - if to_add: - # otherwise add the item as it is - added_item = structure_cif_block0.add_item(item) + merge_mmcifgen_into_structure( + structure_cif_block0, mmcifgen_block, xtal_name, cmpd_code, cmpd_codes_dict, pose_ids_dict + ) data_processing_log_file = None data_processing_prog = row.get(Constants.SOAKDB_COL_DATA_PROCESSING_PROGRAM) @@ -1101,7 +1128,14 @@ def read_fragalysis_csv(filename): return d -def run(collator_path, metadata_csv, compound_codes_csv=None, fragalysis_csv=None, debug=False): +def run( + collator_path, + metadata_csv, + compound_codes_csv=None, + fragalysis_csv=None, + mmcifgen_runner=run_mmcifgen, + debug=False, +): info('run on ' + str(datetime.datetime.now())) info('using RDKit version ' + rdBase.rdkitVersion) # info('using InCHI version ' + Chem.GetInchiVersion()) @@ -1127,7 +1161,7 @@ def run(collator_path, metadata_csv, compound_codes_csv=None, fragalysis_csv=Non # run mmcif-gen input_name = 'input' + str(i + 1) - run_mmcifgen(input_name, str(soakdb_file_p), metadata_csv, output_dir_p) + mmcifgen_runner(input_name, str(soakdb_file_p), metadata_csv, output_dir_p) meta_doc = cif.read(str(output_dir_p / (input_name + '_model.cif'))) meta_mmcifgen = meta_doc[0] diff --git a/tests/test_pdb_deposition.py b/tests/test_pdb_deposition.py new file mode 100644 index 0000000..4a2fd63 --- /dev/null +++ b/tests/test_pdb_deposition.py @@ -0,0 +1,274 @@ +import textwrap +from pathlib import Path + +import pytest +from gemmi import cif + +from pdbdepo.pdb_deposition import ( + merge_mmcifgen_into_structure, + read_cmpd_codes, + read_fragalysis_csv, + substitute_tokens, +) + + +# --------------------------------------------------------------------------- +# substitute_tokens — pure unit tests +# --------------------------------------------------------------------------- + + +def test_substitute_tokens_compound_and_crystal(): + result = substitute_tokens('$CompoundCode in $CrystalName', 'XTAL001', 'ABC123', None, '') + assert result == 'ABC123 in XTAL001' + + +def test_substitute_tokens_pose_id_single(): + result = substitute_tokens('keywords, $PoseID', 'XTAL001', 'ABC', None, 'XTAL001a') + assert result == 'keywords, XTAL001a' + + +def test_substitute_tokens_pose_id_multiple(): + result = substitute_tokens('keywords, $PoseID', 'XTAL001', 'ABC', None, 'XTAL001a, XTAL001b') + assert result == 'keywords, XTAL001a, XTAL001b' + + +def test_substitute_tokens_external_code_correct_column(): + # ext_codes[0] -> $ExternalCode2, ext_codes[1] -> $ExternalCode3 + result = substitute_tokens('id $ExternalCode3', 'XTAL001', 'ABC', ['OB-001', 'OB-002'], '') + assert result == 'id OB-002' + + +def test_substitute_tokens_erases_unsubstituted_external_codes(): + result = substitute_tokens('($ExternalCode3)', 'XTAL001', 'ABC', None, '') + assert result == '()' + + +def test_substitute_tokens_no_ext_codes_erases_all(): + result = substitute_tokens('$ExternalCode2 $ExternalCode5', 'X', 'C', None, '') + assert result == ' ' + + +def test_substitute_tokens_no_tokens_unchanged(): + result = substitute_tokens('no tokens here', 'XTAL001', 'ABC', None, '') + assert result == 'no tokens here' + + +def test_substitute_tokens_empty_pose_id_leaves_no_token(): + result = substitute_tokens('keys, $PoseID', 'XTAL001', 'ABC', None, '') + assert '$PoseID' not in result + assert result == 'keys, ' + + +# --------------------------------------------------------------------------- +# read_fragalysis_csv +# --------------------------------------------------------------------------- + + +def test_read_fragalysis_csv_happy_path(tmp_path): + csv = tmp_path / 'meta.csv' + csv.write_text( + textwrap.dedent( + """\ + Code,Long code,Experiment code,Compound code + A71EV2A-x0836a,A71EV2A-x0836_A_301_0_v1,A71EV2A-x0836,Z123 + A71EV2A-x0836b,A71EV2A-x0836_A_302_0_v1,A71EV2A-x0836,Z123 + A71EV2A-x0001a,A71EV2A-x0001_A_1_0_v1,A71EV2A-x0001,Z456 + """ + ) + ) + result = read_fragalysis_csv(str(csv)) + assert result == { + 'A71EV2A-x0836': ['A71EV2A-x0836a', 'A71EV2A-x0836b'], + 'A71EV2A-x0001': ['A71EV2A-x0001a'], + } + + +def test_read_fragalysis_csv_quoted_fields(tmp_path): + """Quoted CSV fields (as produced by the real Fragalysis download) must be parsed correctly.""" + csv = tmp_path / 'meta.csv' + csv.write_text( + textwrap.dedent( + """\ + Code,Long code,Experiment code,Compound code + "A71EV2A-x0836a","A71EV2A-x0836_A_301_0_v1","A71EV2A-x0836","Z123" + "A71EV2A-x0836b","A71EV2A-x0836_A_302_0_v1","A71EV2A-x0836","Z123" + """ + ) + ) + result = read_fragalysis_csv(str(csv)) + assert result == {'A71EV2A-x0836': ['A71EV2A-x0836a', 'A71EV2A-x0836b']} + + +def test_read_fragalysis_csv_none_filename(): + assert read_fragalysis_csv(None) == {} + + +def test_read_fragalysis_csv_skips_short_rows(tmp_path): + csv = tmp_path / 'meta.csv' + csv.write_text('Code,Long code,Experiment code\nA-x0001a,long\n') # only 2 columns on data row + result = read_fragalysis_csv(str(csv)) + assert result == {} + + +# --------------------------------------------------------------------------- +# read_cmpd_codes +# --------------------------------------------------------------------------- + + +def test_read_cmpd_codes_happy_path(tmp_path): + csv = tmp_path / 'codes.csv' + csv.write_text( + textwrap.dedent( + """\ + CrystalName,CompoundCode,OpenBindId + XTAL001,Z001,OB-001 + XTAL002,Z002,OB-002 + """ + ) + ) + result = read_cmpd_codes(str(csv)) + assert result == {'XTAL001': ['Z001', 'OB-001'], 'XTAL002': ['Z002', 'OB-002']} + + +def test_read_cmpd_codes_none_filename(): + assert read_cmpd_codes(None) == {} + + +# --------------------------------------------------------------------------- +# merge_mmcifgen_into_structure — CIF-level tests (in-memory gemmi objects) +# --------------------------------------------------------------------------- + + +def _make_mmcifgen_doc(): + """Return a minimal mmcif-gen CIF document with a _struct.title loop containing tokens.""" + doc = cif.Document() + block = doc.add_new_block('mmcifgen') + loop = block.init_loop('', ['_struct.entry_id', '_struct.title']) + loop.add_row(['INVID', 'structure of $CrystalName with $CompoundCode ($ExternalCode3)']) + return doc + + +def _make_struct_block(): + """Return an empty destination CIF block.""" + doc = cif.Document() + return doc.add_new_block('output') + + +def test_merge_mmcifgen_title_token_substitution(): + mmcifgen_block = _make_mmcifgen_doc()[0] + struct_block = _make_struct_block() + + merge_mmcifgen_into_structure( + struct_block, + mmcifgen_block, + xtal_name='XTAL042', + cmpd_code='Z999', + cmpd_codes_dict={'XTAL042': ['COL2', 'OB-042']}, + pose_ids_dict={}, + ) + + title_item = next(item for item in struct_block if item.loop is not None and '_struct.title' in item.loop.tags) + title = title_item.loop.values[1] + assert 'XTAL042' in title + assert 'Z999' in title + assert 'OB-042' in title + assert '$' not in title + + +def test_merge_mmcifgen_title_erases_unused_external_code(): + mmcifgen_block = _make_mmcifgen_doc()[0] + struct_block = _make_struct_block() + + merge_mmcifgen_into_structure( + struct_block, + mmcifgen_block, + xtal_name='XTAL042', + cmpd_code='Z999', + cmpd_codes_dict={}, # no external codes for this crystal + pose_ids_dict={}, + ) + + title_item = next(item for item in struct_block if item.loop is not None and '_struct.title' in item.loop.tags) + title = title_item.loop.values[1] + assert '$ExternalCode3' not in title + assert '$' not in title + + +def test_merge_mmcifgen_keywords_loop_pose_id(): + """A keywords loop containing $PoseID in any value has it substituted.""" + doc = cif.Document() + block = doc.add_new_block('mmcifgen') + loop = block.init_loop( + '', ['_struct_keywords.entry_id', '_struct_keywords.text', '_struct_keywords.pdbx_keywords'] + ) + loop.add_row(['INVID', 'Diamond Light Source, $PoseID', 'VIRAL PROTEIN']) + + struct_block = _make_struct_block() + merge_mmcifgen_into_structure( + struct_block, + block, + xtal_name='XTAL042', + cmpd_code='Z999', + cmpd_codes_dict={}, + pose_ids_dict={'XTAL042': ['XTAL042a', 'XTAL042b']}, + ) + + kw_item = next( + item for item in struct_block if item.loop is not None and '_struct_keywords.text' in item.loop.tags + ) + idx = kw_item.loop.tags.index('_struct_keywords.text') + text_value = kw_item.loop.values[idx] + assert 'XTAL042a, XTAL042b' in text_value + assert '$PoseID' not in text_value + + +def test_merge_mmcifgen_keywords_pair_pose_id(): + """A pair item whose value contains $PoseID has it substituted.""" + doc = cif.Document() + block = doc.add_new_block('mmcifgen') + block.set_pair('_struct_keywords.text', 'Diamond Light Source, $PoseID') + + struct_block = _make_struct_block() + merge_mmcifgen_into_structure( + struct_block, + block, + xtal_name='XTAL042', + cmpd_code='Z999', + cmpd_codes_dict={}, + pose_ids_dict={'XTAL042': ['XTAL042a']}, + ) + + kw_pair = next(item for item in struct_block if item.pair is not None and '_struct_keywords.text' in item.pair[0]) + assert 'XTAL042a' in kw_pair.pair[1] + assert '$PoseID' not in kw_pair.pair[1] + + +def test_merge_mmcifgen_passthrough_item_added_unchanged(): + """Items with no tokens are added to the destination block as-is.""" + doc = cif.Document() + block = doc.add_new_block('mmcifgen') + block.set_pair('_exptl.method', 'X-RAY DIFFRACTION') + + struct_block = _make_struct_block() + merge_mmcifgen_into_structure(struct_block, block, 'XTAL', 'Z1', {}, {}) + + exptl = next((item for item in struct_block if item.pair is not None and item.pair[0] == '_exptl.method'), None) + assert exptl is not None + assert exptl.pair[1] == 'X-RAY DIFFRACTION' + + +def test_merge_mmcifgen_no_pose_ids_for_crystal(): + """A crystal not present in pose_ids_dict produces an empty substitution — no crash, no token.""" + doc = cif.Document() + block = doc.add_new_block('mmcifgen') + loop = block.init_loop('', ['_struct_keywords.entry_id', '_struct_keywords.text']) + loop.add_row(['INVID', 'keywords, $PoseID']) + + struct_block = _make_struct_block() + merge_mmcifgen_into_structure(struct_block, block, 'MISSING_XTAL', 'Z1', {}, {}) + + kw_item = next( + item for item in struct_block if item.loop is not None and '_struct_keywords.text' in item.loop.tags + ) + idx = kw_item.loop.tags.index('_struct_keywords.text') + assert '$PoseID' not in kw_item.loop.values[idx] diff --git a/tests/test_validate_structure_cif_doc.py b/tests/test_validate_structure_cif_doc.py new file mode 100644 index 0000000..af559f4 --- /dev/null +++ b/tests/test_validate_structure_cif_doc.py @@ -0,0 +1,178 @@ +from pathlib import Path + +import pytest +from gemmi import cif + +from pdbdepo.pdb_deposition import validate_structure_cif_doc + +CIF_DIR = Path(__file__).parent.parent / "test-data" / "cif" + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + + +def _minimal_valid_doc(): + """Build a CIF document that passes all validation checks.""" + doc = cif.Document() + block = doc.add_new_block('TEST') + + block.set_pair('_cell.length_a', '50.0') + block.set_pair('_symmetry.space_group_name_H-M', 'P 21 21 21') + + loop = block.init_loop('', ['_atom_site.id', '_atom_site.label_entity_id']) + loop.add_row(['1', '1']) + + loop = block.init_loop('', ['_entity.id', '_entity.type']) + loop.add_row(['1', 'polymer']) + + loop = block.init_loop( + '', + [ + '_refine.entry_id', + '_refine.ls_R_factor_R_work', + '_refine.ls_R_factor_R_free', + '_refine.ls_d_res_high', + '_refine.ls_d_res_low', + '_refine.B_iso_mean', + ], + ) + loop.add_row(['TEST', '0.20', '0.25', '1.5', '50.0', '20.0']) + + loop = block.init_loop('', ['_software.name', '_software.pdbx_ordinal']) + loop.add_row(['REFMAC', '1']) + + return doc + + +# --------------------------------------------------------------------------- +# File-based tests (existing) +# --------------------------------------------------------------------------- + + +def test_disconnected_loop_raises_on_parse(): + """gemmi rejects a CIF where the same loop category appears in two separate, non-contiguous loops.""" + with pytest.raises(RuntimeError, match="duplicate tag"): + cif.read(str(CIF_DIR / "disconnected_loop.cif")) + + +def test_duplicate_pair_raises_on_parse(): + """gemmi rejects a CIF with a duplicate pair tag.""" + with pytest.raises(RuntimeError, match="duplicate tag"): + cif.read(str(CIF_DIR / "duplicate_pair.cif")) + + +def test_disconnected_pair_no_validation_issues(): + """A CIF with a disconnected pair (same category appearing both before and after a loop) parses without error and passes validation.""" + doc = cif.read(str(CIF_DIR / "disconnected_pair.cif")) + issues = validate_structure_cif_doc(doc) + assert issues == [] + + +# --------------------------------------------------------------------------- +# In-memory validation tests +# --------------------------------------------------------------------------- + + +def test_minimal_valid_doc_passes(): + assert validate_structure_cif_doc(_minimal_valid_doc()) == [] + + +def test_missing_required_loop_category_reported(): + doc = _minimal_valid_doc() + # Remove the _refine loop + block = doc[0] + for item in block: + if item.loop is not None and any(t.startswith('_refine.') for t in item.loop.tags): + item.erase() + break + issues = validate_structure_cif_doc(doc) + assert any('_refine' in i for i in issues) + + +def test_missing_required_pair_category_reported(): + doc = _minimal_valid_doc() + block = doc[0] + for item in block: + if item.pair is not None and item.pair[0].startswith('_cell.'): + item.erase() + break + issues = validate_structure_cif_doc(doc) + assert any('_cell' in i for i in issues) + + +def _doc_with_rfactors(r_work, r_free): + """Build a minimal valid doc with specific R-factor values.""" + doc = cif.Document() + block = doc.add_new_block('TEST') + block.set_pair('_cell.length_a', '50.0') + block.set_pair('_symmetry.space_group_name_H-M', 'P 21 21 21') + loop = block.init_loop('', ['_atom_site.id', '_atom_site.label_entity_id']) + loop.add_row(['1', '1']) + loop = block.init_loop('', ['_entity.id', '_entity.type']) + loop.add_row(['1', 'polymer']) + loop = block.init_loop( + '', + [ + '_refine.entry_id', + '_refine.ls_R_factor_R_work', + '_refine.ls_R_factor_R_free', + '_refine.ls_d_res_high', + '_refine.ls_d_res_low', + '_refine.B_iso_mean', + ], + ) + loop.add_row(['TEST', str(r_work), str(r_free), '1.5', '50.0', '20.0']) + loop = block.init_loop('', ['_software.name', '_software.pdbx_ordinal']) + loop.add_row(['REFMAC', '1']) + return doc + + +def test_rfactor_out_of_range_reported(): + doc = _doc_with_rfactors(r_work=1.5, r_free=0.25) # R_work > 1.0 — invalid + issues = validate_structure_cif_doc(doc) + assert any('R_work' in i for i in issues) + + +def test_rfree_less_than_rwork_reported(): + doc = _doc_with_rfactors(r_work=0.30, r_free=0.20) # R_free < R_work — invalid + issues = validate_structure_cif_doc(doc) + assert any('R_free' in i and 'R_work' in i for i in issues) + + +def test_unsubstituted_placeholder_in_title_reported(): + doc = _minimal_valid_doc() + block = doc[0] + loop = block.init_loop('', ['_struct.entry_id', '_struct.title']) + loop.add_row(['TEST', 'structure of $CompoundCode ($CrystalName)']) + issues = validate_structure_cif_doc(doc) + assert any('placeholder' in i for i in issues) + + +def test_software_ordinal_sequential_passes(): + doc = _minimal_valid_doc() + block = doc[0] + # Replace single-entry software loop with a proper multi-entry sequential one + for item in block: + if item.loop is not None and '_software.pdbx_ordinal' in item.loop.tags: + item.erase() + break + loop = block.init_loop('', ['_software.name', '_software.pdbx_ordinal']) + loop.add_row(['REFMAC', '1']) + loop.add_row(['CCP4', '2']) + assert validate_structure_cif_doc(doc) == [] + + +def test_software_ordinal_non_sequential_reported(): + doc = _minimal_valid_doc() + block = doc[0] + for item in block: + if item.loop is not None and '_software.pdbx_ordinal' in item.loop.tags: + item.erase() + break + loop = block.init_loop('', ['_software.name', '_software.pdbx_ordinal']) + loop.add_row(['REFMAC', '1']) + loop.add_row(['CCP4', '3']) # gap — should be 2 + issues = validate_structure_cif_doc(doc) + assert any('ordinal' in i for i in issues)