Add calcSignatureCavities/calcSignatureChannels for ensemble signatur… by stoneandbone · Pull Request #2268 · prody/ProDy

stoneandbone · 2026-08-11T23:03:58Z

Added calcSignatureCavities/calcSignatureChannels for ensemble signature detection.

It would superpose homologous structures onto a reference via ProDy's own CE-align, auto-discover homologues via DALI when none are given, and compute a consensus signature cavity/channel across the ensemble mapped onto of the reference. The reference is colored by per-residue Shannon entropy, computed either from an auto-built MSA or a user-supplied one (msa_fasta). DALI's discovery filter cutoffs are overridable via dali_filter_kwargs.

18 unit tests in prody/tests/proteins/test_channels.py
Validated with curated list of PDB structures and/or automatic DaliSearch for non-receptor tyrosine kinases, fibroblast growth factor receptor kinases, and Lox families during development.

…e detection Superpose homologous structures onto a reference via ProDy's own CE-align, auto-discover homologues via DALI when none are given, and compute a consensus signature cavity/channel across the ensemble. The reference is colored by per-residue Shannon entropy, computed either from an auto-built MSA or a user-supplied one (msa_fasta). DALI's discovery filter cutoffs are overridable via dali_filter_kwargs. 18 tests in prody/tests/proteins/test_channels.py, all local-fixture-seeded (no network required).

karolamik13

I tested basic options. I had a problem with MSA data when directly fetched from Pfam.

I used this code:
from prody import *

ref_structure = '1tqnA'

DALI

dali_rec = searchDali('1tqn', 'A')
while not dali_rec.isSuccess:
dali_rec.fetch()
targets = dali_rec.filter(cutoff_len=0.8, cutoff_rmsd=2.0, cutoff_Z=20)
targets = targets[:10]
if ref_structure not in targets:
targets.insert(0, ref_structure)

Pfam

pfam_ID = searchPfam(ref_structure)
pfam_acc = list(pfam_ID.keys())[0]
print("Pfam:", pfam_acc)
pfam_sth = fetchPfamMSA(pfam_acc, alignment='full', outname='CYP3A4')
msa = parseMSA(pfam_sth)
print("MSA:", msa)
print("Number of sequences:", msa.numSequences())
print("Number of columns:", msa.numResidues())

print("P08684:", [x for x in msa.getLabels() if 'P08684' in x])
print("CP3A4_HUMAN:", [x for x in msa.getLabels() if 'CP3A4_HUMAN' in x])
msa_fasta = writeMSA('CYP3A4_full.fasta', msa)

Signature channels

result = calcSignatureChannels(
structures=targets,
ref_structure=ref_structure,
msa_fasta=msa_fasta,
output_path='ensemble_out_ch3b')

When msa_fasta=msa_fasta is not used, then the code works fine.
Am I doing something wrong?

karolamik13 · 2026-08-26T09:01:01Z

@jamesmkrieger maybe you would like to take a look too?

jamesmkrieger

This feels like it has too much that doesn't need to be there. Lots of reinventing the wheel and nested functions. I made some comments on some of these.

There also seems to be a lot of parsePDB calls as well

It's probably better to build an actual PDBEnsemble, although I understand that this may not work if not all the structures have the same atoms.

Try to simplify and modularise and have as many generalisable pieces as possible rather than just one big piece of code that needs maintaining and understanding and is only usable for a specific purpose

jamesmkrieger · 2026-08-27T15:07:53Z


            norm_count = float(count) / float(len(pqr_files))

+            # PDB serials are a fixed 5-column field; wrap past 99999 rather than


no, don't have you're own pdb writer for this. Use existing pdb format strings that use hexadecimal or hybrid36

jamesmkrieger · 2026-08-27T15:12:22Z

    return axes[0] if len(axes) == 1 else axes


+def _calcSignatureEnsemble(structures, ref_structure, output_path, voxel_res,


if this has is specific to channels/cavities/etc then give it a name related to that

This is general to both channels and cavities. They follow a similar analytic process with different sets of functions. So I thought it may worth clumping them together instead of writing the separatly.

jamesmkrieger · 2026-08-27T15:30:26Z

I'd also say don't base everything on CE-align as it performs worse than DALI in many cases. If you have dali mappings then let them be usable and let MSA's drive the mappings too. So much of this machinery already exists in buildPDBEnsemble and SignDy already

jamesmkrieger · 2026-08-27T15:31:38Z

I tested basic options. I had a problem with MSA data when directly fetched from Pfam.

I used this code: from prody import *

ref_structure = '1tqnA'

DALI

dali_rec = searchDali('1tqn', 'A') while not dali_rec.isSuccess: dali_rec.fetch() targets = dali_rec.filter(cutoff_len=0.8, cutoff_rmsd=2.0, cutoff_Z=20) targets = targets[:10] if ref_structure not in targets: targets.insert(0, ref_structure)

Pfam

pfam_ID = searchPfam(ref_structure) pfam_acc = list(pfam_ID.keys())[0] print("Pfam:", pfam_acc) pfam_sth = fetchPfamMSA(pfam_acc, alignment='full', outname='CYP3A4') msa = parseMSA(pfam_sth) print("MSA:", msa) print("Number of sequences:", msa.numSequences()) print("Number of columns:", msa.numResidues())

print("P08684:", [x for x in msa.getLabels() if 'P08684' in x]) print("CP3A4_HUMAN:", [x for x in msa.getLabels() if 'CP3A4_HUMAN' in x]) msa_fasta = writeMSA('CYP3A4_full.fasta', msa)

Signature channels

result = calcSignatureChannels( structures=targets, ref_structure=ref_structure, msa_fasta=msa_fasta, output_path='ensemble_out_ch3b')

When msa_fasta=msa_fasta is not used, then the code works fine. Am I doing something wrong?

This looks reasonable. What error do you get with msa_fasta=msa_fasta?

karolamik13 · 2026-08-27T21:35:10Z

I tested basic options. I had a problem with MSA data when directly fetched from Pfam.
I used this code: from prody import *
ref_structure = '1tqnA'

DALI

dali_rec = searchDali('1tqn', 'A') while not dali_rec.isSuccess: dali_rec.fetch() targets = dali_rec.filter(cutoff_len=0.8, cutoff_rmsd=2.0, cutoff_Z=20) targets = targets[:10] if ref_structure not in targets: targets.insert(0, ref_structure)

Pfam

pfam_ID = searchPfam(ref_structure) pfam_acc = list(pfam_ID.keys())[0] print("Pfam:", pfam_acc) pfam_sth = fetchPfamMSA(pfam_acc, alignment='full', outname='CYP3A4') msa = parseMSA(pfam_sth) print("MSA:", msa) print("Number of sequences:", msa.numSequences()) print("Number of columns:", msa.numResidues())
print("P08684:", [x for x in msa.getLabels() if 'P08684' in x]) print("CP3A4_HUMAN:", [x for x in msa.getLabels() if 'CP3A4_HUMAN' in x]) msa_fasta = writeMSA('CYP3A4_full.fasta', msa)

Signature channels

result = calcSignatureChannels( structures=targets, ref_structure=ref_structure, msa_fasta=msa_fasta, output_path='ensemble_out_ch3b')
When msa_fasta=msa_fasta is not used, then the code works fine. Am I doing something wrong?

This looks reasonable. What error do you get with msa_fasta=msa_fasta?

This is the error:
4832 best_aln = (aligned_ref, aligned_row)
4834 if best_idx is None:
-> 4835 raise ValueError(
4836 "msa_fasta: no row in the supplied MSA could be aligned "
4837 "to the reference structure {0!r}. Review or supply an "
4838 "MSA that includes or aligns to the reference structure."
4839 .format(ref_label))
4841 mismatches = len(ref_seq) - best_matches
4842 mismatch_pct = 100.0 * mismatches / max(len(ref_seq), 1)

ValueError: msa_fasta: no row in the supplied MSA could be aligned to the reference structure '1tqnA'. Review or supply an MSA that includes or aligns to the reference structure.

stoneandbone · 2026-08-27T21:44:39Z

I know what may be the issue. Let me incorporate @jamesmkrieger comments as well as a fix for @karolamik13 .

My concern with buildPDBEnsemble is that it works best with only calphas while the calculation for detection of cavities/channels perform best with full atomic. The resulted MSA from full atomic PDB ensemble (from my prior experience) will be repeated single-letter codes (e.g. a run of Leucine side-chain atoms all reading as L) for every single position. I will work to see if I can get around with this. It may be also nice to return a PDBensemble for downstream SignDy workflow anyhow :).

jamesmkrieger · 2026-08-28T08:47:14Z

It would be nice to have an argument on PDBEnsemble getMSA that returns based on residues rather than atoms.

I suppose the bigger issue here is that unless all your structures have the same sequence and have all been sufficiently well modelled and fixed then you can’t actually build a good ensemble.

We had a similar problem for signature interactions and managed to find a way around it though, so maybe you can reuse and extend the code from there

jamesmkrieger · 2026-08-28T08:50:07Z

Something else that might help is my extension of clustenm to take multiple structures as input and apply pdbfixer to all of them. I think it’s still limited to structures that all have the same residues but that could potentially be generalised further too

- calcChannelSurfaceOverlaps now writes overlap voxels through AtomGroup + writePDB(hybrid36=True) instead of a hand-rolled PDB line writer. - _calcSignatureEnsemble batches all per-structure parsePDB calls into one parsePDB(list_of_labels) call instead of one call per structure. - msa_fasta now goes through prody's own parseMSA/calcShannonEntropy instead of a manual FASTA parse and a hand-built MSA array. - Document threshold, distA, voxel_res, output_path, structures, ref_structure and kwargs on _calcSignatureEnsemble, which only documented kind/detect_records_fn/overlap_fn before. Also fixes the Pfam-MSA anchoring bug karolamik13 reported: a Pfam 'full' alignment labels its rows by UniProt idcode/range (e.g. 'CP3A4_HUMAN/38-493'), not by accession, so accession-only matching never found the reference's row. Try idcode before accession, matching how refineMSA already resolves DBREF labels. Also widen the compiled-parser fallback to catch ValueError, not just OSError/IOError.

writeFasta() (prody/sequence/msaio.c) measured the MSA array's total character buffer with strlen() and stored the result in an int. The buffer is numseq*lenseq raw bytes with no null terminators between rows, so strlen() was never the right way to measure it, and once that product exceeds INT_MAX (as with a real Pfam 'full' family alignment, e.g. PF00067 at ~362k rows x ~6k columns = ~2.19 billion characters) the truncated/garbage result left every count < lenmsa check false, so writeMSA silently produced a FASTA file of headers followed by empty sequences, with no error raised. This was the actual root cause behind a bug report on PR prody#2268: a Pfam full-alignment FASTA produced by fetchPfamMSA + writeMSA had zero-length sequences for every row (including the reference's own), so calcSignatureChannels/calcSignatureCavities could never align a row to the reference. Use numseq*lenseq directly as the buffer length instead of strlen(), and widen the index/length variables to long.

parsePDB already resolves and restricts each 'PDBID'+'ChainID' label structures are now passed straight through to the batch parsePDB() call, and per-structure Chain objects are taken directly from each already chain-restricted AtomGroup's hierview instead of re-indexing by the chain letter. Also tightens _getMappingResnums's docstring.

writeFasta() (prody/sequence/msaio.c) measured the MSA array's total character buffer with strlen() and stored the result in an int. The buffer is numseq*lenseq raw bytes with no null terminators between rows, so strlen() was never the right way to measure it, and once that product exceeds INT_MAX (as with a real Pfam 'full' family alignment, e.g. PF00067 at ~362k rows x ~6k columns = ~2.19 billion characters) the truncated/garbage result left every count < lenmsa check false, so writeMSA silently produced a FASTA file of headers followed by empty sequences, with no error raised. This was the actual root cause behind a bug report on PR prody#2268: a Pfam full-alignment FASTA produced by fetchPfamMSA + writeMSA had zero-length sequences for every row (including the reference's own), so calcSignatureChannels/calcSignatureCavities could never align a row to the reference. Use numseq*lenseq directly as the buffer length instead of strlen(), and widen the index/length variables to long. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>

…ity-channel-ensemble # Conflicts: # prody/proteins/channels.py

stoneandbone · 2026-09-25T01:04:54Z

Dear @karolamik13

I tracked this down to two separate bugs, both now fixed on this branch:

Reference row never matched (prody/proteins/channels.py) — the reference structure's row in msa_fasta was looked up by UniProt accession only (e.g. P08684). And the Pfam 'full' Stockholm alignment labels it differently by UniProt idcode/mnemonic + range instead (e.g. CP3A4_HUMAN/38-493)

Another issue about the MSA: corrupted FASTA output (prody/sequence/msaio.c) — writeMSA's compiled FASTA writer measured its output buffer with strlen() and stored the result in a 32-bit int. For a 'full' Pfam alignment like PF00067 (362,663 rows × 6,052 columns ≈ 2.19 billion characters, just over INT_MAX), therefore writeMSA wrote out headers followed by empty sequences.

I reran your exact script against 1tqnA/CYP3A4 and it completes successfully now.

karolamik13 · 2026-09-25T10:04:04Z

Dear @karolamik13

I tracked this down to two separate bugs, both now fixed on this branch:

Reference row never matched (prody/proteins/channels.py) — the reference structure's row in msa_fasta was looked up by UniProt accession only (e.g. P08684). And the Pfam 'full' Stockholm alignment labels it differently by UniProt idcode/mnemonic + range instead (e.g. CP3A4_HUMAN/38-493)

Another issue about the MSA: corrupted FASTA output (prody/sequence/msaio.c) — writeMSA's compiled FASTA writer measured its output buffer with strlen() and stored the result in a 32-bit int. For a 'full' Pfam alignment like PF00067 (362,663 rows × 6,052 columns ≈ 2.19 billion characters, just over INT_MAX), therefore writeMSA wrote out headers followed by empty sequences.

I reran your exact script against 1tqnA/CYP3A4 and it completes successfully now.

Thank you. I will test it today or during the weekend.

karolamik13

I tested those functions, and now I am able to analyze the structure in a correct way.

I would have three suggestions:
(1) Could you include a parameter with selection in your functions? I am analyzing 1TQN, which contains a heme group, but some people might want to analyze "clean" protein or just remove some crystalographic molecules/inhibitors, etc., from their structures. The current CaviTracer code includes everything that is not water, so if you have inhibitors in PDBs or some other small molecules, peptides, etc., in your structures, you will not be able to properly analyze the whole family.

(2) Also, it is advisable to save in a file information about a list of structures, detection parameters, threshold, resolution, MSA source, selected reference sequence, and the status of each structure: aligned, detection successful, missing channels, or error.

(3) Except for the ref.pdb file, where you have conservation in the beta column, I would suggest saving a PNG file and raw data with the results from Shannon Entropy.

Apart from that, I am happy with the code.

@jamesmkrieger Do you have any suggestions?

jamesmkrieger · 2026-10-02T09:47:37Z

I’d say things like saving pngs is for a script running the code, not the core code itself.

Otherwise it also makes sense. I can have a look in more detail late next week when the new changes come in. I am very busy right now

… ensembles Address the remaining review requests on PR prody#2268. selection: calcSignatureCavities and calcSignatureChannels take a ProDy selection string (default 'protein', the previous behaviour). The atoms it picks are written to the aligned and signature PDB files and are what detection sees, so a cofactor such as a heme can be kept with 'protein or resname HEM', or an inhibitor dropped. CE-align, the superposition fit, residue mapping, the MSA and the entropy stay on protein atoms only; kept hetero atoms move with the protein but does not affect the superposition. Entropy table: the per-residue Shannon entropy is written to signature_<kind>_entropy.csv (resnum, resname, entropy; NaN where the MSA has no coverage) and returned as result['entropy']. Run summary: signature_<kind>_run_summary.json records the reference label and sequence, the MSA source, path and anchor row, the run parameters, and one entry per input structure with its alignment metrics and detection status (ok, none_found, error or skipped) and message. A structure that cannot be aligned or fails detection is now reported instead of only logged. No images are written. Plotting the entropy is left to the calling script, as suggested in James' comment.

stoneandbone · 2026-10-24T17:22:12Z

Thanks both for the reviews.
Dear @karolamik13, Your three new suggestions:

selection (default 'protein', same as before), for example selection='protein or resname HEM'. Alignment, the fit, the MSA and the entropy stay protein-only.
signature_<kind>_run_summary.json lists the reference, the MSA source and anchor row, the parameters, and each structure’s alignment metrics and detection status (ok, none_found, error, skipped).
Entropy is saved as signature_<kind>_entropy.csv and returned as result['entropy']. Note that, there is no PNG in the library, per James’s comment. A script can plot it with ax.plot(result['entropy']['resnum'], result['entropy']['entropy']).

stoneandbone assigned stoneandbone and karolamik13 and unassigned stoneandbone and karolamik13 Aug 11, 2026

stoneandbone requested a review from karolamik13 August 12, 2026 00:50

karolamik13 requested changes Aug 26, 2026

View reviewed changes

jamesmkrieger reviewed Aug 27, 2026

View reviewed changes

stoneandbone added 3 commits September 24, 2026 17:09

Merge remote-tracking branch 'origin/main' into feature/signature-cav…

48adfb7

…ity-channel-ensemble # Conflicts: # prody/proteins/channels.py

stoneandbone requested a review from MatthewLicht September 25, 2026 01:21

karolamik13 requested changes Oct 2, 2026

View reviewed changes

stoneandbone force-pushed the feature/signature-cavity-channel-ensemble branch from c6760e2 to 48adfb7 Compare October 24, 2026 16:36


		norm_count = float(count) / float(len(pqr_files))

		# PDB serials are a fixed 5-column field; wrap past 99999 rather than

		return axes[0] if len(axes) == 1 else axes


		def _calcSignatureEnsemble(structures, ref_structure, output_path, voxel_res,

Conversation

stoneandbone commented Aug 11, 2026

Uh oh!

karolamik13 left a comment

Choose a reason for hiding this comment

DALI

Pfam

Signature channels

Uh oh!

karolamik13 commented Aug 26, 2026

Uh oh!

jamesmkrieger left a comment

Choose a reason for hiding this comment

Uh oh!

jamesmkrieger Aug 27, 2026

Choose a reason for hiding this comment

Uh oh!

jamesmkrieger Aug 27, 2026

Choose a reason for hiding this comment

Uh oh!

stoneandbone Aug 27, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

jamesmkrieger commented Aug 27, 2026

Uh oh!

jamesmkrieger commented Aug 27, 2026

DALI

Pfam

Signature channels

Uh oh!

karolamik13 commented Aug 27, 2026

DALI

Pfam

Signature channels

Uh oh!

stoneandbone commented Aug 27, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jamesmkrieger commented Aug 28, 2026

Uh oh!

jamesmkrieger commented Aug 28, 2026

Uh oh!

stoneandbone commented Sep 25, 2026

Uh oh!

karolamik13 commented Sep 25, 2026

Uh oh!

karolamik13 left a comment

Choose a reason for hiding this comment

Uh oh!

jamesmkrieger commented Oct 2, 2026

Uh oh!

stoneandbone commented Oct 24, 2026

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

stoneandbone commented Aug 27, 2026 •

edited

Loading