from prody import *
from pylab import *
%matplotlib inline

"""
Explanation: Sequence-Structure Comparison with Evol
The part shows how to compare sequence conservation properties with structural mobility obtained from Gaussian network model (GNM) calculations.
Setting up environment
End of explanation
"""


fetchPfamMSA('PF00074')

"""
Explanation: Entropy Calculation
First, we retrieve MSA for protein for protein family PF00074:
End of explanation
"""


msa = parseMSA('PF00074_full.sth')

"""
Explanation: We parse the MSA file:
End of explanation
"""


msa_refine = refineMSA(msa, label='RNAS1_BOVIN', rowocc=0.8, seqid=0.98)

"""
Explanation: Then, we refine it using refineMSA() based on the sequence of RNAS1_BOVIN:
End of explanation
"""


entropy = calcShannonEntropy(msa_refine)

"""
Explanation: We calculate entropy for refined MSA.
End of explanation
"""


pdb = parsePDB('2W5I', chain='B')

chB_ca = pdb.select('protein and name CA and resid 1 to 119')

"""
Explanation: Mobility calculation
Next, we obtain residue fluctuations or mobility for protein member of the above family. We will use chain B of 2W5I.
End of explanation
"""


gnm = GNM('2W5I')

gnm.buildKirchhoff(chB_ca)

gnm.calcModes(n_modes=None)  # calculate all modes

mobility_1 = calcSqFlucts(gnm[0])

mobility_1to8 = calcSqFlucts(gnm[:8])

mobility_all = calcSqFlucts(gnm[:])

"""
Explanation: We perform GNM as follows:
End of explanation
"""


indices = range(1,120)

bar(indices, entropy, width=1.2, color='grey', hold='True');

xlim(min(indices)-1, max(indices)+1);

plot(indices, mobility_all*(max(entropy)/mean(mobility_all)), color='b',
linewidth=2);

xlabel('residue index')
ylabel('mobility/entropy')

"""
Explanation: Plotting
End of explanation
"""


