diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index 8833e0827..23f29698b 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -14,9 +14,10 @@ import numpy as np from prody import LOGGER, PY3K -from prody.atomic import Atomic +from prody.atomic import Atomic, Chain from prody.utilities import getCoords, isListLike from prody.proteins import writePDB, parsePDB, parsePQR +from prody.proteins.compare import mapChainOntoChain from prody.ensemble import Ensemble from prody.measure import calcCenter, calcTransformation, calcDistance, calcRMSD, superpose @@ -37,7 +38,8 @@ 'calcPoresFromChannelsMultipleFrames', 'getPoreParametersMultipleFrames', 'getPoreResidueNamesMultipleFrames', 'scanChannelParameters', 'scanSurfaceCavityParameters', - 'calcFrequentObjectResidues', 'showFrequentObjectResidues'] + 'calcFrequentObjectResidues', 'showFrequentObjectResidues', + 'calcSignatureCavities', 'calcSignatureChannels'] # Sampling of the enclosure test used to strip the moat (see # ChannelCalculator.calcEnclosure). These are constants, not knobs: the enclosure @@ -3830,8 +3832,11 @@ def calcChannelSurfaceOverlaps(**kwargs): norm_count = float(count) / float(len(pqr_files)) + # PDB serials are a fixed 5-column field; wrap past 99999 rather than + # overflow the column and corrupt every field after it (atom_id is + # never used to identify a FIL voxel downstream). out.write("ATOM {:5d} H FIL T 1 {:8.3f}{:8.3f}{:8.3f}{:6.2f} 1.00\n" - .format(atom_id, x, y, z, norm_count)) + .format(atom_id % 100000, x, y, z, norm_count)) atom_id += 1 @@ -4561,6 +4566,755 @@ def showFrequentObjectResidues(counts_by_chain, top=50): return axes[0] if len(axes) == 1 else axes +def _calcSignatureEnsemble(structures, ref_structure, output_path, voxel_res, + threshold, distA, kwargs, kind, + detect_records_fn, overlap_fn, + dali_filter_kwargs=None, msa_fasta=None): + """Shared alignment/thresholding pipeline behind :func:`calcSignatureCavities` + and :func:`calcSignatureChannels`. The two public functions differ only in + which detection family they dispatch to. + + :arg kind: ``'cavity'`` or ``'channel'`` -- selects output file naming + (``signature__pt{threshold}.pdb``, ``aligned__params.csv``, + ``_overlap.pdb``) and the ``'_records'`` / + ``'signature__residues'`` result keys. + :type kind: str + + :arg detect_records_fn: Called as ``detect_records_fn(prot, det_out, label)`` + once per aligned structure; runs the family's own detection and + parameter/residue extraction and returns a list of record dicts + (each must include a ``'pdb'`` key) to append to the records table. + :type detect_records_fn: callable + + :arg overlap_fn: :func:`calcSurfaceCavityOverlaps` or + :func:`calcChannelSurfaceOverlaps`. + :type overlap_fn: callable + """ + import pandas as pd + import re + from prody.database import searchDali + + if PY3K: + from pathlib import Path + else: + from pathlib2 import Path + + out_dir = Path(output_path).resolve() + align_dir = out_dir / 'aligned' + unaligned_dir = out_dir / 'unaligned' + cav_dir = align_dir / 'cavities' + overlap_dir = out_dir / 'overlaps' + for d in (align_dir, unaligned_dir, cav_dir, overlap_dir): + d.mkdir(parents=True, exist_ok=True) + + def _getMappingResnums(pdb_file, chain=None): + """Return ordered Cα residue numbers for a single chain in a PDB file. + + :arg chain: Chain ID string (``'A'``), or ``None`` to use all protein + chains. + """ + ag = parsePDB(str(pdb_file)) + if ag is None: + return [] + if chain is not None: + ca = ag.select('protein and name CA and chain ' + str(chain).strip()) + else: + ca = ag.select('protein and name CA') + return list(ca.getResnums()) if ca is not None else [] + + def _one_letter(resname): + from prody.atomic.atomic import AAMAP + if resname in ('HSD', 'HSP'): + return AAMAP.get('HIS', resname[0]) + return AAMAP.get(resname, resname[0] if resname else '?') + + + # ── 1. Resolve reference structure string, if given ──────────────────────── + ref_pdb_id = ref_chain = ref_label = None + if ref_structure is not None: + ref_structure_str = str(ref_structure) + if len(ref_structure_str) < 5: + raise ValueError( + "ref_structure must be a 4-character PDB code followed by " + "a chain identifier (e.g. '2srcA'); got {0!r}." + .format(ref_structure_str)) + ref_pdb_id, ref_chain = ref_structure_str[:4], ref_structure_str[4:] + ref_label = ref_structure_str + + # ── 2. Resolve structures, discovering via DALI when not supplied ───────── + if structures is None: + if ref_pdb_id is None: + raise ValueError("Either structures or ref_structure must be provided.") + + dali_rec = searchDali(ref_pdb_id, ref_chain) + while not dali_rec.isSuccess: + LOGGER.info('Waiting on DALI search results for {0}...'.format(ref_label)) + dali_rec.fetch() + + filter_kwargs = {'cutoff_len': 0.5, 'cutoff_rmsd': 2.0, + 'cutoff_Z': 8, 'stringency': True} + if dali_filter_kwargs: + filter_kwargs.update(dali_filter_kwargs) + pdb_ids = dali_rec.filter(**filter_kwargs) + if not pdb_ids: + filter_desc = ', '.join('{0}={1}'.format(k, v) + for k, v in filter_kwargs.items()) + _warn('DALI search for {0} returned no hits passing {1}. ' + 'Provide an explicit structures list instead.' + .format(ref_label, filter_desc)) + return None + structures = list(pdb_ids) + elif not isinstance(structures, list): + raise ValueError("structures must be a list of 'pdbid'+'chain' strings, " + "e.g. ['2srcA', '1bbhA'].") + + entries = [] + for item in structures: + item_str = str(item) + if len(item_str) < 5: + raise ValueError( + "Each element of structures must be a 4-character PDB code " + "followed by a chain identifier (e.g. '2srcA'); got {0!r}." + .format(item_str)) + pdb_id, chain = item_str[:4], item_str[4:] + entries.append((pdb_id, chain, item_str)) + + if not entries: + raise ValueError("No structures found in the provided input.") + + if ref_pdb_id is None: + ref_pdb_id, ref_chain, ref_label = entries[0] + + _ag_ref = parsePDB(ref_pdb_id, chain=ref_chain, folder=str(unaligned_dir)) + if _ag_ref is None: + raise ValueError("Could not parse reference: {0} chain {1}".format(ref_pdb_id, ref_chain)) + _ag_ref_prot = _ag_ref.select('protein') + if _ag_ref_prot is None: + raise ValueError("No protein atoms in reference: {0}".format(ref_pdb_id)) + ref_pdb = align_dir / 'ref.pdb' + writePDB(str(ref_pdb), _ag_ref_prot) + + ref_resnums = _getMappingResnums(ref_pdb, ref_chain) + + # ── 3. CE-align all structures onto the reference ───────────────────────── + report_rows = [] + msa_mappings = {} + + ref_chain_obj = _ag_ref_prot.getHierView()[ref_chain] + if not isinstance(ref_chain_obj, Chain): + raise ValueError("Chain {0} not found in reference: {1}" + .format(ref_chain, ref_pdb_id)) + + for pdb_id, chain, label in entries: + ag = parsePDB(pdb_id, chain=chain, folder=str(unaligned_dir)) + if ag is None: + _warn('No atoms parsed for {0}'.format(label)) + continue + + ag_chain = ag.select('protein') + if ag_chain is None: + _warn('No protein atoms in chain {0} of {1}'.format(chain, pdb_id)) + continue + + mobile_chain_obj = ag_chain.getHierView()[chain] + if not isinstance(mobile_chain_obj, Chain): + _warn('Chain {0} not found in {1}.'.format(chain, pdb_id)) + continue + + mapping = mapChainOntoChain(mobile_chain_obj, ref_chain_obj, + pwalign='cealign', overlap=50.) + if mapping is None: + _warn('CE-align could not map {0} onto the reference.'.format(label)) + continue + + atommap, selection, seqid, cover = mapping + mapped_flags = atommap.getFlags('mapped') + superpose(atommap, selection, weights=mapped_flags) + rmsd = float(calcRMSD(atommap, selection, weights=mapped_flags)) + + # Structures stay in the ensemble regardless of alignment quality -- + # RMSD/coverage/identity are reported, not gated. Only a grossly bad + # superposition (>2 Angstroms) gets flagged. + if rmsd > 2.0: + _warn('{0}: superposition RMSD is {1:.2f} A (>2 Angstroms) -- ' + 'alignment may be less reliable.'.format(label, rmsd)) + + sup_pdb = align_dir / '{0}.pdb'.format(label) + writePDB(str(sup_pdb), ag_chain) + + ca_mask = selection.getNames() == 'CA' + ref_rn = selection.getResnums()[ca_mask] + tgt_rn = atommap.getResnums()[ca_mask] + ca_mapped = mapped_flags[ca_mask] + msa_mappings[label] = { + int(rn): (int(tn) if is_mapped else None) + for rn, tn, is_mapped in zip(ref_rn, tgt_rn, ca_mapped) + } + + report_rows.append({ + 'pdb_chain': label, + 'rmsd_A': round(rmsd, 3), + 'coverage_pct': round(float(cover), 2), + 'identity_pct': round(float(seqid), 2), + 'superposed_pdb': str(sup_pdb), + }) + LOGGER.info('{0} RMSD={1:.3f}A coverage={2:.1f}% identity={3:.1f}%'.format( + label, rmsd, cover, seqid)) + + + df_report = pd.DataFrame(report_rows) + df_report.to_csv(out_dir / 'alignment_report.csv', index=False) + + msa_path = out_dir / 'aligned_resnums.msa' + with open(msa_path, 'w') as f: + f.write('{}\t{}\n'.format(ref_label, ' '.join(str(r) for r in ref_resnums))) + for lbl, mapping in msa_mappings.items(): + vals = ' '.join( + str(mapping[r]) if mapping.get(r) is not None else '-' + for r in ref_resnums + ) + f.write('{}\t{}\n'.format(lbl, vals)) + + # ── 4. Reference entropy coloring + sequence MSA output ─────────────────── + from prody.sequence.msa import MSA + from prody.sequence.analysis import calcShannonEntropy + + msa_fasta_path = None + + if msa_fasta: + # User-supplied MSA replaces the auto-generated MSA from the ensemble of aligned structures. + from Bio import SeqIO + + user_labels = [] + user_seqs = [] + with open(str(msa_fasta)) as f: + for rec in SeqIO.parse(f, 'fasta'): + user_labels.append(rec.id) + user_seqs.append(str(rec.seq)) + + ref_ca_for_msa = _ag_ref_prot.select('protein and name CA') + ref_seq = ref_ca_for_msa.getSequence() if ref_ca_for_msa is not None else '' + + col_to_resnum = {} + + if ref_label in user_labels: + # Exact label match: assumed to correspond 1:1 to the + # reference's own residues, + ref_row_seq = user_seqs[user_labels.index(ref_label)] + ref_idx = 0 + for col, ch in enumerate(ref_row_seq): + if ch not in ('-', '.'): + if ref_idx < len(ref_resnums): + col_to_resnum[col] = ref_resnums[ref_idx] + ref_idx += 1 + else: + # No exact label: pairwise-align the reference's own sequence + # against every row, pick the highest-identity match, to + # recover a ref_resnum -> column-index mapping. + from prody.utilities.seqtools import alignBioPairwise + + best_idx = None + best_aln = None + best_matches = -1 + for idx, seq in enumerate(user_seqs): + ungapped = seq.replace('-', '').replace('.', '') + if not ungapped: + continue + alns = alignBioPairwise(ref_seq, ungapped) + if not alns: + continue + aligned_ref, aligned_row = alns[0][0], alns[0][1] + matches = sum(1 for a, b in zip(aligned_ref, aligned_row) + if a == b and a not in ('-', '.')) + if matches > best_matches: + best_matches = matches + best_idx = idx + best_aln = (aligned_ref, aligned_row) + + if best_idx is None: + raise ValueError( + "msa_fasta: no row in the supplied MSA could be aligned " + "to the reference structure {0!r}. Review or supply an " + "MSA that includes or aligns to the reference structure." + .format(ref_label)) + + mismatches = len(ref_seq) - best_matches + mismatch_pct = 100.0 * mismatches / max(len(ref_seq), 1) + if mismatch_pct > 30.0: + raise ValueError( + "msa_fasta: best-matching row {0!r} for reference {1!r} " + "has {2:.1f}% mismatch (> 30%). Review or supply an MSA " + "that includes or aligns to the reference structure." + .format(user_labels[best_idx], ref_label, mismatch_pct)) + + _warn("msa_fasta: no row labeled {0!r}; using best-matching row " + "{1!r} instead ({2:.1f}% mismatch)." + .format(ref_label, user_labels[best_idx], mismatch_pct)) + + # Map ungapped-position -> column-index within the best-matching + # row's own entry in the user's MSA. + row_col_map = {} + upos = 0 + for col, ch in enumerate(user_seqs[best_idx]): + if ch not in ('-', '.'): + row_col_map[upos] = col + upos += 1 + + aligned_ref, aligned_row = best_aln + ref_idx = 0 + row_upos = 0 + for a, b in zip(aligned_ref, aligned_row): + is_ref_gap = a in ('-', '.') + is_row_gap = b in ('-', '.') + if not is_ref_gap and not is_row_gap and ref_idx < len(ref_resnums): + col = row_col_map.get(row_upos) + if col is not None: + col_to_resnum[col] = ref_resnums[ref_idx] + if not is_ref_gap: + ref_idx += 1 + if not is_row_gap: + row_upos += 1 + + n_cols = max(len(seq) for seq in user_seqs) if user_seqs else 0 + msa_labels = user_labels + msa_array = np.full((len(msa_labels), n_cols), b'-', dtype='S1') + for i, seq in enumerate(user_seqs): + for col, ch in enumerate(seq[:n_cols]): + msa_array[i, col] = ch.encode() + + ensemble_msa = MSA(msa_array, title='signature_{0}_msa'.format(kind), + labels=msa_labels) + entropy = calcShannonEntropy(ensemble_msa, omitgaps=True) + entropy_by_resnum = {} + for col, e in enumerate(entropy): + rn = col_to_resnum.get(col) + if rn is not None: + entropy_by_resnum[rn] = float(e) + else: + homolog_labels = list(msa_mappings.keys()) + msa_labels = [ref_label] + homolog_labels + msa_array = np.full((len(msa_labels), len(ref_resnums)), b'-', dtype='S1') + + ref_ca = _ag_ref_prot.select('name CA') + ref_rn_to_1l = {} + if ref_ca is not None: + ref_rn_to_1l = {int(rn): _one_letter(rname) for rn, rname in + zip(ref_ca.getResnums(), ref_ca.getResnames())} + for j, rn in enumerate(ref_resnums): + if rn in ref_rn_to_1l: + msa_array[0, j] = ref_rn_to_1l[rn].encode() + + for i, label in enumerate(homolog_labels, start=1): + aligned_pdb = align_dir / '{0}.pdb'.format(label) + hom_ag = parsePDB(str(aligned_pdb)) if aligned_pdb.exists() else None + hom_rn_to_1l = {} + if hom_ag is not None: + hom_ca = hom_ag.select('protein and name CA') + if hom_ca is not None: + hom_rn_to_1l = {int(rn): _one_letter(rname) for rn, rname in + zip(hom_ca.getResnums(), hom_ca.getResnames())} + mapping = msa_mappings[label] + for j, rn in enumerate(ref_resnums): + tgt_rn = mapping.get(rn) + if tgt_rn is not None and tgt_rn in hom_rn_to_1l: + msa_array[i, j] = hom_rn_to_1l[tgt_rn].encode() + + msa_fasta_path = out_dir / 'signature_{0}_msa.fasta'.format(kind) + with open(str(msa_fasta_path), 'w') as f: + for lbl, row in zip(msa_labels, msa_array): + seq = ''.join(c.decode() for c in row) + f.write('>{0}\n{1}\n'.format(lbl, seq)) + + ensemble_msa = MSA(msa_array, title='signature_{0}_msa'.format(kind), + labels=msa_labels) + entropy = calcShannonEntropy(ensemble_msa, omitgaps=True) + entropy_by_resnum = {int(rn): float(e) for rn, e in zip(ref_resnums, entropy)} + + ref_resnums_arr = _ag_ref_prot.getResnums() + ref_betas = _ag_ref_prot.getBetas() + new_betas = np.array([entropy_by_resnum.get(int(rn), b) + for rn, b in zip(ref_resnums_arr, ref_betas)]) + _ag_ref_prot.setBetas(new_betas) + writePDB(str(ref_pdb), _ag_ref_prot) + if msa_fasta_path is not None: + LOGGER.info('Reference coloured by Shannon entropy; MSA written: {0}'.format( + msa_fasta_path.name)) + else: + LOGGER.info('Reference coloured by Shannon entropy from user-supplied MSA.') + + # ── 5. Detection ──────────────────────────────────────────────────────── + records = [] + + for row in df_report.itertuples(): + label = row.pdb_chain + sup_pdb = Path(row.superposed_pdb) + if not sup_pdb.exists(): + _warn('Superposed PDB not found: {0}'.format(sup_pdb.name)) + continue + + ag = parsePDB(str(sup_pdb)) + sel = ag.select('protein') if ag is not None else None + prot = sel.copy() if sel is not None else None + if prot is None: + _warn('No protein atoms in {0}'.format(sup_pdb.name)) + continue + + det_out = str(cav_dir / '{0}_{1}'.format(kind, label)) + try: + records.extend(detect_records_fn(prot, det_out, label)) + except Exception as e: + _warn('CaviTracer {0}: {1}'.format(label, e)) + + df_records = pd.DataFrame(records) + df_records.to_csv(out_dir / 'aligned_{0}_params.csv'.format(kind), index=False) + + # ── 6. Surface overlaps ─────────────────────────────────────────────────── + overlap_pdb = str(overlap_dir / '{0}_overlap.pdb'.format(kind)) + if any(cav_dir.glob('*.pqr')): + overlap_fn( + pqr_files=str(cav_dir), + output_file_name=overlap_pdb, resolution=voxel_res + ) + else: + _warn('No PQR files found in {0}; skipping surface overlap.'.format(cav_dir)) + overlap_pdb = None + + # ── 7. Threshold consensus voxels and combine with the reference ────────── + threshold = int(threshold) + if not 0 <= threshold <= 100: + raise ValueError('threshold must be between 0 and 100, got {0}'.format(threshold)) + + sig_coords = np.empty((0, 3)) + sig_occupancies = np.empty((0,)) + if overlap_pdb is not None: + sig_atoms = parsePDB(overlap_pdb) + fil = sig_atoms.select('resname FIL') if sig_atoms is not None else None + if fil is not None: + threshold_frac = threshold / 100.0 + mask = fil.getOccupancies() >= threshold_frac + sig_coords = fil.getCoords()[mask] + sig_occupancies = fil.getOccupancies()[mask] + + sig_pdb_path = out_dir / 'signature_{0}_pt{1}.pdb'.format(kind, threshold) + writePDB(str(sig_pdb_path), _ag_ref_prot) + with open(str(sig_pdb_path)) as f: + ref_lines = f.readlines() + if ref_lines and ref_lines[-1].startswith('END'): + ref_lines = ref_lines[:-1] + + with open(str(sig_pdb_path), 'w') as f: + f.writelines(ref_lines) + start_serial = _ag_ref_prot.numAtoms() + 1 + for i, ((x, y, z), occ) in enumerate(zip(sig_coords, sig_occupancies)): + f.write("ATOM {:5d} H FIL A 1 " + "{:8.3f}{:8.3f}{:8.3f} 1.00{:6.2f}\n".format( + (start_serial + i) % 100000, x, y, z, occ * 100.0)) + f.write('END\n') + + LOGGER.info('Signature {0} written: {1} ({2} voxels, threshold={3}%)'.format( + kind, sig_pdb_path.name, len(sig_coords), threshold)) + + result = { + 'msa_mappings': msa_mappings, + 'alignment_report': df_report, + '{0}_records'.format(kind): df_records, + 'ensemble_voxels': overlap_pdb, + 'signature_pdb': str(sig_pdb_path), + } + if msa_fasta_path is not None: + result['msa_fasta_path'] = str(msa_fasta_path) + + # ── 8. Map reference residues near signature voxels onto each homologue ── + if len(sig_coords) > 0: + ca_coords = _ag_ref_prot.select('name CA').getCoords() + ca_resnums = _ag_ref_prot.select('name CA').getResnums() + ca_resnames = _ag_ref_prot.select('name CA').getResnames() + + dist2_cut = distA ** 2 + near_idx = set() + chunk = 2000 + for i in range(0, len(sig_coords), chunk): + blk = sig_coords[i:i + chunk] + d2 = ((ca_coords[:, None, :] - blk[None, :, :]) ** 2).sum(axis=2) + near_idx.update(np.where(d2.min(axis=1) <= dist2_cut)[0].tolist()) + + sorted_indices = sorted(near_idx, key=lambda i: int(ca_resnums[i])) + ref_near_resnums = [int(ca_resnums[i]) for i in sorted_indices] + ref_str = ', '.join( + '{0}{1}'.format(_one_letter(ca_resnames[i]), int(ca_resnums[i])) + for i in sorted_indices + ) if sorted_indices else 'None' + + signature_residues = {'ref': ref_str} + LOGGER.info('Signature {0}: {1} reference residues within {2} A'.format( + kind, len(near_idx), distA)) + + ref_resnums_set = set(ref_near_resnums) + for label, mapping in msa_mappings.items(): + mapped_resnums = sorted( + tgt_rn for ref_rn, tgt_rn in mapping.items() + if ref_rn in ref_resnums_set and tgt_rn is not None) + if not mapped_resnums: + signature_residues[label] = 'None' + continue + + aligned_pdb = align_dir / '{0}.pdb'.format(label) + hom_ag = parsePDB(str(aligned_pdb)) if aligned_pdb.exists() else None + rn_to_1l = {} + if hom_ag is not None: + hom_ca = hom_ag.select('protein and name CA') + if hom_ca is not None: + rn_to_1l = {int(rn): _one_letter(rname) for rn, rname in + zip(hom_ca.getResnums(), hom_ca.getResnames())} + parts = ['{0}{1}'.format(rn_to_1l[rn], rn) if rn in rn_to_1l else str(rn) + for rn in mapped_resnums] + signature_residues[label] = ', '.join(parts) + + residues_path = out_dir / 'signature_{0}_pt{1}_residues.txt'.format(kind, threshold) + with open(str(residues_path), 'w') as f: + for key, val in signature_residues.items(): + f.write('{0}: {1}\n'.format(key, val)) + + result['signature_{0}_residues'.format(kind)] = signature_residues + result['signature_{0}_residues_path'.format(kind)] = str(residues_path) + LOGGER.info('Signature {0} residues written: {1}'.format(kind, residues_path.name)) + + return result + + +def calcSignatureCavities(structures=None, ref_structure=None, + output_path='output', + voxel_res=0.5, threshold=20, distA=4.5, + dali_filter_kwargs=None, msa_fasta=None, **kwargs): + """Calculate surface cavities for an ensemble of homologous protein structures. + + Structures are superposed onto a reference using ProDy's own CE-align + (:func:`.mapChainOntoChain` with ``pwalign='cealign'``), surface-cavity + detection (:func:`calcSurfaceCavities`) is run on each aligned structure, + and a consensus signature cavity is computed across the ensemble. + + :arg structures: A list of ``'pdbid'+'chain'`` strings (e.g. ``'2srcA'``). + *pdbid* is a 4-character PDB accession code loaded via + :func:`.parsePDB`; *chain* is a single chain identifier. Each string + is used directly as its structure's label. If ``None``, *ref_structure* + must be given, and homologues are discovered automatically via + :func:`.searchDali` (``cutoff_len=0.5, cutoff_rmsd=2.0, cutoff_Z=8, + stringency=True``, overridable via *dali_filter_kwargs*); DALI's own + residue mapping is not used, every homologue is realigned with + CE-align regardless of how it was found. + :type structures: list of str or None + + :arg ref_structure: ``'pdbid'+'chain'`` string for the reference + structure used for CE-align superposition and residue-number mapping. + If ``None`` and *structures* is given, the first entry in *structures* + is used as reference. Required when *structures* is ``None``. + :type ref_structure: str or None + + :arg output_path: Root directory for all output files. + Subdirectories ``aligned/``, ``aligned/cavities/``, and ``overlaps/`` + are created automatically. default is ``'output'`` + :type output_path: str + + :arg voxel_res: Voxel resolution in Å for the consensus cavity map, passed + to :func:`calcSurfaceCavityOverlaps` as ``resolution``. + default is ``0.5`` + :type voxel_res: float + + Additional keyword arguments are forwarded to :func:`calcSurfaceCavities`. + + :arg threshold: Minimum occupancy (integer percent, 0-100) for a voxel to + be retained in the signature cavity. default is ``20`` + :type threshold: int + + :arg dali_filter_kwargs: Overrides merged over the default DALI filter + cutoffs (``cutoff_len=0.5, cutoff_rmsd=2.0, cutoff_Z=8, + stringency=True``) before calling :meth:`.DaliRecord.filter`. Only + used when *structures* is ``None`` (DALI auto-discovery); ignored + otherwise. default is ``None`` + :type dali_filter_kwargs: dict or None + + :arg msa_fasta: Path to a user-supplied FASTA-format MSA. When given, it + replaces the auto-built MSA as the Shannon entropy source. A row anchoring + to *ref_structure* is required: an exact label match is used if present, + otherwise the best-matching row by sequence identity against the reference + structure is used, with a warning. If the best match still has + more than 30% mismatch against the reference (or no row aligns at + all), :exc:`ValueError` is raised. Structural superposition and + mapping (CE-align, ``msa_mappings``) are unaffected either way. + default is ``None`` + :type msa_fasta: str or None + + :arg distA: Distance cutoff in Å for identifying reference residues near + signature voxels. default is ``4.5`` + :type distA: float + + :returns: Dictionary with keys: + + * ``msa_mappings`` — ``{label: {ref_resnum: tgt_resnum}}`` + * ``alignment_report`` — DataFrame with RMSD, coverage, and identity + per structure + * ``cavity_records`` — DataFrame with volume, depth, tetrahedra count, + and lining residues per cavity + * ``ensemble_voxels`` — path to the combined-voxel PDB written by + :func:`calcSurfaceCavityOverlaps` + * ``signature_pdb`` — path to the reference structure combined with + the thresholded consensus voxels, named with the threshold (e.g. + ``signature_cavity_pt20.pdb``) + * ``signature_cavity_residues`` — ``{key: "X123, Y456, ..."}``, + present only when at least one voxel passes *threshold*; ``'ref'`` + holds reference residues within *distA* of a signature voxel, and + each homologue label holds its mapped residues + * ``signature_cavity_residues_path`` — path to the written residue + summary text file (same condition as above) + + :rtype: dict + + Example usage:: + + result = calcSignatureCavities( + structures=targets, + ref_structure='2srcA', + msa_fasta='my_msa.fasta', + output_path='ensemble_out', + ) + print(result['signature_pdb']) + + """ + def _detectRecords(prot, det_out, label): + cavities, surface = calcSurfaceCavities(prot, output_path=det_out, + min_depth=2, separate=False, max_depth=None, **kwargs) + volumes, depths, tetrahedra_counts = getSurfaceCavityParameters(cavities) + cavity_residues = getSurfaceCavityResidueNames( + prot, cavities, surface, distA=4, one_letter_aa=True) + rows = [] + for cav_id, (vol, dep, ntet) in enumerate(zip(volumes, depths, tetrahedra_counts)): + res_str = None + if cavity_residues and cav_id < len(cavity_residues): + parts = cavity_residues[cav_id].split(':', 1) + res_str = parts[1].strip() if len(parts) == 2 else parts[0].strip() + rows.append({ + 'pdb': label, + 'cavity_id': cav_id, + 'volume_A3': round(vol, 3), + 'depth_A': round(dep, 3), + 'tetrahedra_count': ntet, + 'cavity_residues': res_str, + }) + return rows + + return _calcSignatureEnsemble(structures, ref_structure, output_path, + voxel_res, threshold, distA, kwargs, + 'cavity', _detectRecords, + calcSurfaceCavityOverlaps, + dali_filter_kwargs=dali_filter_kwargs, + msa_fasta=msa_fasta) + + +def calcSignatureChannels(structures=None, ref_structure=None, + output_path='output', + voxel_res=0.5, threshold=20, distA=4.5, + dali_filter_kwargs=None, msa_fasta=None, **kwargs): + """Calculate channels for an ensemble of homologous protein structures. + + Identical pipeline to :func:`calcSignatureCavities` -- structures are + superposed onto a reference using ProDy's own CE-align -- but dispatches + to the channel-detection family (:func:`calcChannels`) instead of the + surface-cavity family, producing a consensus signature *channel* (a + through-going tunnel) rather than a signature cavity (an enclosed pocket). + + :arg structures: See :func:`calcSignatureCavities`. + :type structures: list of str or None + + :arg ref_structure: See :func:`calcSignatureCavities`. + :type ref_structure: str or None + + :arg output_path: Root directory for all output files. + Subdirectories ``aligned/``, ``aligned/cavities/``, and ``overlaps/`` + are created automatically. default is ``'output'`` + :type output_path: str + + :arg voxel_res: Voxel resolution in Å for the consensus channel map, + passed to :func:`calcChannelSurfaceOverlaps` as ``resolution``. + default is ``0.5`` + :type voxel_res: float + + Additional keyword arguments are forwarded to :func:`calcChannels`. + + :arg threshold: Minimum occupancy (integer percent, 0-100) for a voxel to + be retained in the signature channel. default is ``20`` + :type threshold: int + + :arg dali_filter_kwargs: See :func:`calcSignatureCavities`. + :type dali_filter_kwargs: dict or None + + :arg msa_fasta: See :func:`calcSignatureCavities`. + :type msa_fasta: str or None + + :arg distA: Distance cutoff in Å for identifying reference residues near + signature voxels. default is ``4.5`` + :type distA: float + + :returns: Dictionary with keys: + + * ``msa_mappings`` — ``{label: {ref_resnum: tgt_resnum}}`` + * ``alignment_report`` — DataFrame with RMSD, coverage, and identity + per structure + * ``channel_records`` — DataFrame with length, bottleneck, volume, and + lining residues per channel + * ``ensemble_voxels`` — path to the combined-voxel PDB written by + :func:`calcChannelSurfaceOverlaps` + * ``signature_pdb`` — path to the reference structure combined with + the thresholded consensus voxels, named with the threshold (e.g. + ``signature_channel_pt20.pdb``) + * ``signature_channel_residues`` — ``{key: "X123, Y456, ..."}``, + present only when at least one voxel passes *threshold*; ``'ref'`` + holds reference residues within *distA* of a signature voxel, and + each homologue label holds its mapped residues + * ``signature_channel_residues_path`` — path to the written residue + summary text file (same condition as above) + + :rtype: dict + + Example usage:: + + result = calcSignatureChannels( + structures=targets, + ref_structure='2srcA', + output_path='ensemble_out', + ) + print(result['signature_pdb']) + + """ + def _detectRecords(prot, det_out, label): + channels, _details = calcChannels(prot, output_path=det_out, + separate=False, **kwargs) + lengths, bottlenecks, volumes = getChannelParameters(channels) + channel_residues = getChannelResidueNames( + prot, channels, distA=4, one_letter_aa=True) + rows = [] + for ch_id, (ln, bn, vol) in enumerate(zip(lengths, bottlenecks, volumes)): + res_str = None + if channel_residues and ch_id < len(channel_residues): + parts = channel_residues[ch_id].split(':', 1) + res_str = parts[1].strip() if len(parts) == 2 else parts[0].strip() + rows.append({ + 'pdb': label, + 'channel_id': ch_id, + 'length_A': round(ln, 3), + 'bottleneck_A': round(bn, 3), + 'volume_A3': round(vol, 3), + 'channel_residues': res_str, + }) + return rows + + return _calcSignatureEnsemble(structures, ref_structure, output_path, + voxel_res, threshold, distA, kwargs, + 'channel', _detectRecords, + calcChannelSurfaceOverlaps, + dali_filter_kwargs=dali_filter_kwargs, + msa_fasta=msa_fasta) + + class Channel: def __init__(self, tetrahedra, centerline_spline, radius_spline, length, bottleneck, volume, cost=None): diff --git a/prody/tests/proteins/test_channels.py b/prody/tests/proteins/test_channels.py new file mode 100644 index 000000000..c7a44a2ff --- /dev/null +++ b/prody/tests/proteins/test_channels.py @@ -0,0 +1,440 @@ +"""This module contains unit tests for signature cavity/channel detection in +:mod:`~prody.proteins.channels`.""" + +import os +import shutil +import tempfile +from unittest.mock import MagicMock, patch + +from prody import LOGGER, parsePDB +from prody.proteins.channels import calcSignatureCavities, calcSignatureChannels +from prody.tests import unittest +from prody.tests.datafiles import pathDatafile + +LOGGER.verbosity = 'none' + + +def _seedLocalStructure(unaligned_dir, fixture_file, pdb_id): + """Copy a local datafile fixture into *unaligned_dir* under the filename + ProDy's local-PDB resolution expects for *pdb_id*, so parsePDB(pdb_id, ...) + resolves it without touching the network.""" + + shutil.copy(pathDatafile(fixture_file), + os.path.join(unaligned_dir, 'pdb' + pdb_id + '.pdb')) + + +class TestCalcSignatureCavitiesInput(unittest.TestCase): + + def testRejectsTooShortStructureString(self): + """A structures entry that is too short to contain a 4-character PDB + code plus a chain identifier should be rejected before any PDB + parsing is attempted.""" + + with self.assertRaises(ValueError) as ctx: + calcSignatureCavities(structures=['54'], ref_structure='5a46A') + + self.assertIn('PDB code', str(ctx.exception)) + self.assertIn('chain', str(ctx.exception)) + + def testRejectsTooShortRefStructureString(self): + """A ref_structure that is too short to contain a 4-character PDB + code plus a chain identifier should be rejected the same way as a + malformed structures entry, before any PDB parsing is attempted.""" + + with self.assertRaises(ValueError) as ctx: + calcSignatureCavities(structures=['5a46A'], ref_structure='54') + + self.assertIn('PDB code', str(ctx.exception)) + self.assertIn('chain', str(ctx.exception)) + + +class TestCalcSignatureCavitiesAlignment(unittest.TestCase): + + def setUp(self): + self.out_dir = tempfile.mkdtemp(prefix='sigcav_') + unaligned_dir = os.path.join(self.out_dir, 'unaligned') + os.makedirs(unaligned_dir, exist_ok=True) + # Reference (1ubiA) and one homolog (1ubjA, same coordinates under a + # different fake code) -- both resolved locally, no network. + _seedLocalStructure(unaligned_dir, 'pdb1ubi.pdb', '1ubi') + _seedLocalStructure(unaligned_dir, 'pdb1ubi.pdb', '1ubj') + + def tearDown(self): + shutil.rmtree(self.out_dir, ignore_errors=True) + + def testUsesCealignNotExternalTMalign(self): + """Superposing a homolog onto the reference must not depend on an + external TMalign binary -- it should succeed using ProDy's own + CE-align, even though TMalign isn't installed on this machine.""" + + result = calcSignatureCavities( + structures=['1ubjA'], ref_structure='1ubiA', + output_path=self.out_dir) + + self.assertIsNotNone(result) + + def testSignaturePdbIsThresholdNamedAndIncludesReference(self): + """calcSignatureCavities collapses thresholding into the same call + (no separate getSignatureCavities step), defaults to a 20% threshold, + and writes a combined PDB naming that threshold, containing the + reference structure's own atoms alongside any consensus voxels.""" + + result = calcSignatureCavities( + structures=['1ubjA'], ref_structure='1ubiA', + output_path=self.out_dir) + + self.assertEqual(os.path.basename(result['signature_pdb']), + 'signature_cavity_pt20.pdb') + self.assertTrue(os.path.exists(result['signature_pdb'])) + + written = parsePDB(result['signature_pdb']) + self.assertIsNotNone(written) + self.assertIsNotNone(written.select('protein'), + 'combined signature PDB must include the reference structure, ' + 'not just voxel pseudo-atoms') + + def testResidueMappingNearSignatureVoxels(self): + """When at least one voxel passes the (default 20%) threshold, the + reference residues within distA=4.5 A of it must be reported, and + mapped onto each homologue via msa_mappings, with a residue-summary + file written alongside the signature PDB.""" + + result = calcSignatureCavities( + structures=['1ubjA'], ref_structure='1ubiA', + output_path=self.out_dir) + + written = parsePDB(result['signature_pdb']) + fil = written.select('resname FIL') + n_voxels = 0 if fil is None else len(fil) + + if n_voxels == 0: + self.skipTest('No signature voxels found for this fixture pair; ' + 'residue-mapping is exercised by the real-data run.') + + self.assertIn('signature_cavity_residues', result) + residues = result['signature_cavity_residues'] + self.assertIn('ref', residues) + self.assertIn('1ubjA', residues) + + self.assertIn('signature_cavity_residues_path', result) + self.assertTrue(os.path.exists(result['signature_cavity_residues_path'])) + + def testEntropyColoringAndMsaOutput(self): + """The reference's B-factor column must carry per-residue Shannon + entropy (not left at the original crystallographic B-factors), and + the ref-anchored alignment used to compute it must be written out as + its own FASTA file alongside the numeric resnum mapping.""" + + result = calcSignatureCavities( + structures=['1ubjA'], ref_structure='1ubiA', + output_path=self.out_dir) + + self.assertIn('msa_fasta_path', result) + self.assertTrue(os.path.exists(result['msa_fasta_path'])) + with open(result['msa_fasta_path']) as f: + fasta = f.read() + self.assertIn('>1ubiA', fasta) + self.assertIn('>1ubjA', fasta) + + self.assertTrue(os.path.exists(os.path.join(self.out_dir, 'aligned_resnums.msa'))) + + ref_written = parsePDB(os.path.join(self.out_dir, 'aligned', 'ref.pdb')) + betas = ref_written.select('protein and name CA').getBetas() + # 1ubjA is an identical copy of the reference -- every aligned column + # has zero variation, so entropy must be uniformly ~0, not left at + # whatever crystallographic B-factors pdb1ubi.pdb originally had. + self.assertTrue((betas < 1e-6).all()) + + +class TestCalcSignatureChannels(unittest.TestCase): + + def setUp(self): + self.out_dir = tempfile.mkdtemp(prefix='sigchan_') + unaligned_dir = os.path.join(self.out_dir, 'unaligned') + os.makedirs(unaligned_dir, exist_ok=True) + _seedLocalStructure(unaligned_dir, 'pdb1ubi.pdb', '1ubi') + _seedLocalStructure(unaligned_dir, 'pdb1ubi.pdb', '1ubj') + + def tearDown(self): + shutil.rmtree(self.out_dir, ignore_errors=True) + + def testDispatchesToChannelFamilyNotCavityFamily(self): + """calcSignatureChannels must reuse the same CE-align pipeline as + calcSignatureCavities but dispatch to the channel-detection family, + producing a distinctly-named signature_channel_pt20.pdb with + 'channel_records' (not 'cavity_records') in the result.""" + + result = calcSignatureChannels( + structures=['1ubjA'], ref_structure='1ubiA', + output_path=self.out_dir) + + self.assertIsNotNone(result) + self.assertIn('channel_records', result) + self.assertNotIn('cavity_records', result) + + self.assertEqual(os.path.basename(result['signature_pdb']), + 'signature_channel_pt20.pdb') + self.assertTrue(os.path.exists(result['signature_pdb'])) + + written = parsePDB(result['signature_pdb']) + self.assertIsNotNone(written) + self.assertIsNotNone(written.select('protein'), + 'combined signature PDB must include the reference structure, ' + 'not just voxel pseudo-atoms') + + def testEntropyColoringAndMsaOutput(self): + """calcSignatureChannels must get entropy coloring and MSA output for + free from the shared pipeline, same as calcSignatureCavities.""" + + result = calcSignatureChannels( + structures=['1ubjA'], ref_structure='1ubiA', + output_path=self.out_dir) + + self.assertIn('msa_fasta_path', result) + self.assertTrue(os.path.exists(result['msa_fasta_path'])) + with open(result['msa_fasta_path']) as f: + fasta = f.read() + self.assertIn('>1ubiA', fasta) + self.assertIn('>1ubjA', fasta) + + ref_written = parsePDB(os.path.join(self.out_dir, 'aligned', 'ref.pdb')) + betas = ref_written.select('protein and name CA').getBetas() + self.assertTrue((betas < 1e-6).all()) + + @patch('prody.database.searchDali') + def testDaliFilterKwargsOverridesDefaultCutoffs(self, mock_search_dali): + """dali_filter_kwargs must be accepted by calcSignatureChannels too, + merged over the default cutoffs and passed through to + DaliRecord.filter, same as calcSignatureCavities.""" + + mock_rec = MagicMock() + mock_rec.isSuccess = True + mock_rec.filter.return_value = ['1ubjA'] + mock_search_dali.return_value = mock_rec + + calcSignatureChannels(ref_structure='1ubiA', + output_path=self.out_dir, + dali_filter_kwargs={'cutoff_rmsd': 1.0}) + + mock_rec.filter.assert_called_once_with( + cutoff_len=0.5, cutoff_rmsd=1.0, cutoff_Z=8, stringency=True) + + def testMsaFastaAcceptedForEntropySource(self): + """calcSignatureChannels must accept msa_fasta and use it as the + entropy source, same as calcSignatureCavities.""" + + ag = parsePDB(pathDatafile('pdb1ubi.pdb')) + ref_seq = ag.select('protein and name CA').getSequence() + msa_path = os.path.join(self.out_dir, 'user.fasta') + with open(msa_path, 'w') as f: + f.write('>1ubiA\n{0}\n'.format(ref_seq)) + f.write('>unrelated_seq\n{0}\n'.format(ref_seq)) + + result = calcSignatureChannels( + structures=['1ubjA'], ref_structure='1ubiA', + output_path=self.out_dir, msa_fasta=msa_path) + + self.assertIsNotNone(result) + self.assertNotIn('msa_fasta_path', result) + + ref_written = parsePDB(os.path.join(self.out_dir, 'aligned', 'ref.pdb')) + betas = ref_written.select('protein and name CA').getBetas() + self.assertTrue((betas < 1e-6).all()) + + +class _FakeDaliRecord(object): + """Minimal stand-in for DaliRecord: succeeds after a fixed number of + fetch() calls, and returns a fixed filter() result.""" + + def __init__(self, success_after=0, filter_result=None): + self._success_after = success_after + self.fetch_calls = 0 + self.filter_result = filter_result if filter_result is not None else [] + self.filter_calls = [] + + @property + def isSuccess(self): + return self.fetch_calls >= self._success_after + + def fetch(self): + self.fetch_calls += 1 + + def filter(self, **kwargs): + self.filter_calls.append(kwargs) + return self.filter_result + + +class TestCalcSignatureCavitiesDaliDiscovery(unittest.TestCase): + + def setUp(self): + self.out_dir = tempfile.mkdtemp(prefix='sigcav_dali_') + unaligned_dir = os.path.join(self.out_dir, 'unaligned') + os.makedirs(unaligned_dir, exist_ok=True) + _seedLocalStructure(unaligned_dir, 'pdb1ubi.pdb', '1ubi') + _seedLocalStructure(unaligned_dir, 'pdb1ubi.pdb', '1ubj') + + def tearDown(self): + shutil.rmtree(self.out_dir, ignore_errors=True) + + @patch('prody.database.searchDali') + def testDaliDiscoveryUsesFilteredHitsAsStructures(self, mock_search_dali): + """Omitting structures must trigger a DALI search against the + reference, filtered with the agreed cutoffs, and feed the surviving + hits into the same pipeline a user-supplied list would use.""" + + mock_rec = MagicMock() + mock_rec.isSuccess = True + mock_rec.filter.return_value = ['1ubjA'] + mock_search_dali.return_value = mock_rec + + result = calcSignatureCavities(ref_structure='1ubiA', + output_path=self.out_dir) + + mock_search_dali.assert_called_once_with('1ubi', 'A') + mock_rec.filter.assert_called_once_with( + cutoff_len=0.5, cutoff_rmsd=2.0, cutoff_Z=8, stringency=True) + self.assertIn('1ubjA', result['msa_mappings']) + + @patch('prody.database.searchDali') + def testDaliFilterKwargsOverridesDefaultCutoffs(self, mock_search_dali): + """dali_filter_kwargs must be merged over the default cutoffs and + passed through to DaliRecord.filter, letting a caller loosen or + tighten the DALI search without editing the library.""" + + mock_rec = MagicMock() + mock_rec.isSuccess = True + mock_rec.filter.return_value = ['1ubjA'] + mock_search_dali.return_value = mock_rec + + calcSignatureCavities(ref_structure='1ubiA', + output_path=self.out_dir, + dali_filter_kwargs={'cutoff_rmsd': 1.0}) + + mock_rec.filter.assert_called_once_with( + cutoff_len=0.5, cutoff_rmsd=1.0, cutoff_Z=8, stringency=True) + + @patch('prody.database.searchDali') + def testDaliFilterKwargsHasNoEffectWithExplicitStructures(self, mock_search_dali): + """dali_filter_kwargs must be ignored -- and DALI never invoked at + all -- when an explicit structures list is given.""" + + calcSignatureCavities(structures=['1ubjA'], ref_structure='1ubiA', + output_path=self.out_dir, + dali_filter_kwargs={'cutoff_rmsd': 1.0}) + + mock_search_dali.assert_not_called() + + def testDaliRetryLoopFetchesUntilSuccess(self): + """The isSuccess/fetch retry loop must keep calling fetch() until the + record reports success before ever calling filter().""" + + fake_rec = _FakeDaliRecord(success_after=2, filter_result=['1ubjA']) + with patch('prody.database.searchDali', return_value=fake_rec): + calcSignatureCavities(ref_structure='1ubiA', output_path=self.out_dir) + + self.assertEqual(fake_rec.fetch_calls, 2) + self.assertEqual(len(fake_rec.filter_calls), 1) + + @patch('prody.proteins.channels._warn') + def testDaliZeroHitsReturnsNoneAndWarnsWithFilterSettings(self, mock_warn): + """When filter() returns no hits, calcSignatureCavities must return + None rather than proceed with a degenerate ensemble, and the warning + must name the actual cutoffs used so the caller can act on it.""" + + fake_rec = _FakeDaliRecord(success_after=0, filter_result=[]) + with patch('prody.database.searchDali', return_value=fake_rec): + result = calcSignatureCavities(ref_structure='1ubiA', + output_path=self.out_dir) + + self.assertIsNone(result) + self.assertTrue(mock_warn.called) + warned_msg = mock_warn.call_args[0][0] + self.assertIn('cutoff_len=0.5', warned_msg) + self.assertIn('cutoff_rmsd=2.0', warned_msg) + self.assertIn('cutoff_Z=8', warned_msg) + + +class TestCalcSignatureCavitiesUserSuppliedMsa(unittest.TestCase): + + def setUp(self): + self.out_dir = tempfile.mkdtemp(prefix='sigcav_msa_') + unaligned_dir = os.path.join(self.out_dir, 'unaligned') + os.makedirs(unaligned_dir, exist_ok=True) + _seedLocalStructure(unaligned_dir, 'pdb1ubi.pdb', '1ubi') + _seedLocalStructure(unaligned_dir, 'pdb1ubi.pdb', '1ubj') + + def tearDown(self): + shutil.rmtree(self.out_dir, ignore_errors=True) + + def _refSequence(self): + ag = parsePDB(pathDatafile('pdb1ubi.pdb')) + ca = ag.select('protein and name CA') + return ca.getSequence() + + def testExactLabelMatchUsesSuppliedMsaForEntropy(self): + """When the user's FASTA has a row exactly labeled ref_structure, + that row anchors entropy directly -- no best-effort warning is + needed, and the auto pseudo-MSA build/write is skipped entirely.""" + + ref_seq = self._refSequence() + msa_path = os.path.join(self.out_dir, 'user.fasta') + with open(msa_path, 'w') as f: + f.write('>1ubiA\n{0}\n'.format(ref_seq)) + f.write('>unrelated_seq\n{0}\n'.format(ref_seq)) + + result = calcSignatureCavities( + structures=['1ubjA'], ref_structure='1ubiA', + output_path=self.out_dir, msa_fasta=msa_path) + + self.assertIsNotNone(result) + self.assertNotIn('msa_fasta_path', result) + self.assertFalse(os.path.exists( + os.path.join(self.out_dir, 'signature_cavity_msa.fasta'))) + + ref_written = parsePDB(os.path.join(self.out_dir, 'aligned', 'ref.pdb')) + betas = ref_written.select('protein and name CA').getBetas() + self.assertTrue((betas < 1e-6).all()) + + @patch('prody.proteins.channels._warn') + def testBestEffortMatchWarnsWhenNoExactLabel(self, mock_warn): + """When no row is labeled exactly ref_structure, the best-matching + row by sequence identity against the reference is used instead, + and a warning is emitted since the match was inferred.""" + + ref_seq = self._refSequence() + msa_path = os.path.join(self.out_dir, 'user.fasta') + with open(msa_path, 'w') as f: + f.write('>some_other_label\n{0}\n'.format(ref_seq)) + + result = calcSignatureCavities( + structures=['1ubjA'], ref_structure='1ubiA', + output_path=self.out_dir, msa_fasta=msa_path) + + self.assertIsNotNone(result) + self.assertTrue(mock_warn.called) + + ref_written = parsePDB(os.path.join(self.out_dir, 'aligned', 'ref.pdb')) + betas = ref_written.select('protein and name CA').getBetas() + self.assertTrue((betas < 1e-6).all()) + + def testBestEffortMatchExceedingMismatchThresholdRaises(self): + """When even the best-matching row has >30% mismatch against the + reference, calcSignatureCavities must raise ValueError rather than + silently colouring from an untrustworthy alignment.""" + + ref_seq = self._refSequence() + unrelated_seq = 'A' * len(ref_seq) + msa_path = os.path.join(self.out_dir, 'user.fasta') + with open(msa_path, 'w') as f: + f.write('>some_other_label\n{0}\n'.format(unrelated_seq)) + + with self.assertRaises(ValueError) as ctx: + calcSignatureCavities( + structures=['1ubjA'], ref_structure='1ubiA', + output_path=self.out_dir, msa_fasta=msa_path) + + self.assertIn('mismatch', str(ctx.exception)) + + +if __name__ == '__main__': + unittest.main()