JamesK directional PRS 2#2256
Conversation
Extend calcPerturbResponse's directional (return_vectors) branch with two options: force to apply a specific perturbing direction (a 3-vector or an (m, 3) array) instead of random unit forces, and perturb_node to restrict the scan to selected node(s). Both fill only the perturbed rows, so a partial PRS matrix (and its effectiveness/sensitivity) is returned when return_vectors is False; return_matrix returns both the triple and the response ensemble. Response vectors now carry their squared magnitude as the eigenvalue. Add comparePerturbResponses to overlap a response ensemble against an NMA/ModeSet/Mode/Vector (e.g. an ANM or deformation vector) or another response ensemble, reducing per-force overlaps (max/mean/sum) into a per-node profile. Add unit tests for the new options and helper on the adenylate kinase open->closed transition. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
prs-fixes' changes (the LOGGER.update label fix, the response_matrix -> prs_matrix rename, and the matplotlib version-parse fix in showAtomicMatrix) are already present in this branch, so the perturb.py conflict is resolved by keeping this branch's superset version. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
|
For the glutamate receptor ligand-binding domain closure from open dimer from 1fto to glutamate-bound closed dimer AC from 1ftj, we see top contributions from a sequence around L138 that is part of glutamate binding and around K250 in the vicinity of the flip/flop alternative splicing casette that is also related to glutamate response kinetics. So, I think this validates that it captures something useful. |
karolamik13
left a comment
There was a problem hiding this comment.
I tested the code using the code below. As you can see, I used CA atoms of those structures because they were not equal. Something was missing in one somewhere in the side chain. That is why the results might be different. Btw, structures are aligned via VMD (AB and AC chains).
from prody import *
open_dimer = parsePDB('1fto_aligned.pdb').select('protein and chain A B and resid 5 to 260 and name CA')
closed_dimer = parsePDB('1ftj_aligned.pdb').select('protein and chain A C and resid 5 to 260 and name CA')
target = calcDeformVector(open_dimer, closed_dimer)
anm = ANM('1FTO open')
anm.buildHessian(open_dimer)
anm.calcModes()
responses = calcPerturbResponse(anm, atoms=open_dimer, return_vectors=True, repeats=100)
p = comparePerturbResponses(responses, target, stat='max')
import numpy as np
print(' ')
print('Target vector length:', np.linalg.norm(target.getArray()))
print('RMSD open-closed:', calcRMSD(open_dimer, closed_dimer))
indices = np.argsort(p)[::-1][:20]
print('Chain Residue Score')
for i in indices:
print(f'{open_dimer.getChids()[i]:>5} {open_dimer.getResnames()[i]}{open_dimer.getResnums()[i]:<4} {p[i]:.3f}')
Outcome:
@> 3871 atoms and 1 coordinate set(s) were parsed in 0.04s.
@> 5903 atoms and 1 coordinate set(s) were parsed in 0.06s.
@> Hessian was built in 0.23s.
@> 20 modes were calculated in 0.17s.
@> Calculating perturbation response with 100 repeats
@> Calculating perturbation response [ 99%] 1s
Target vector length: 74.39254136000463
RMSD open-closed: 3.287716904084999
Chain Residue Score
A GLN130 0.631
A LEU127 0.630
A ASP126 0.630
A SER164 0.629
A MET162 0.629
A LYS129 0.627
A ARG163 0.626
A TYR161 0.625
A ALA124 0.625
A SER128 0.625
A THR131 0.625
A SER123 0.624
A GLU125 0.622
B ARG163 0.622
B TRP159 0.622
A ALA165 0.620
B MET158 0.620
B LEU215 0.619
B SER164 0.619
A GLU122 0.619
|
This matches what I saw, thanks. I also used CA atoms and did the matching+alignment with alignChains+superpose in prody. |
I've now redone directional PRS, still with plenty of involvement from Claude but starting it myself and putting things in the right place in the code so everything is unified. There is also a new function comparePerturbResponses that allows overlap comparison to ANM modes and deformation vectors, so that most important residues can be identified more easily.