PCA and Gibbs Free Energy Calculation using Gromacs

[@Zack6ix]

Thank you for this, you're quite inspiring.

[@BioinfoCopilot]

Thank you 😊

[@lucaspaul7208]

Keep doing these are interesting analysis, I am learning a lot from you

[@BioinfoCopilot]

Thank you and yes if time permits 😊

[@AshutoshNayak-d7m]

How did you generate the PCA plot as you have given in the thumbnail image

[@ashutoshnayak5393]

How to get a 3d image of FEL just like the thumbnail, i am only able to generate in 2d ,some other dependencies or packages required.?

[@ashutoshnayak5393]

...

[@Ashrafwazirkhan]

Another masterpiece of work from you. We are learning a lot from you. Sir one request if possible can you upload video on 3D graph of PCA analysis. Further if possible can we correlate the Gibbs free energy with the RMSD. Thanks for wonderful work❤

[@BioinfoCopilot]

Thank you 😊. I will try to do so soon!

[@asif4557]

why we are using Backbone and Ligand. Why not C-Alpha and Ligand or Protein and Ligand, is there any specific reason? or we can consider using these other two options.

[@BioinfoCopilot]

Side chains are subjected to more flexibility meaning giving you a crazy RMSF or RMSD. If you consider backbone then the conformational changes are not that crazy compared to side chains and you get a clear graph indicating the changes in the overall structure of the protein. That’s why side chains are not considered. If you happen to study a particular group of residues then you are welcome to experiment that.

[@LucasPaul-d9q]

Do have any of your videos with docking validation and statistical analysis? If possible I would recommend to do that

[@BioinfoCopilot]

Sure will do that. Thanks for the recommendation.

[@pb_storageBhattacharya]

Hello sir, thank you for the wonderful video. I have a query related to FEL. If we have two FEL plots corresponding to wild type and mutated protein, or lets say an apo-protein and the same protein as holo-protein (bound to some ligand)... how to understand that which system is more stable from their corresponding FEL plots?? If you may kindly elaborate. Thank you

[@BioinfoCopilot]

@@pb_storageBhattacharya Thank you 😊 FEL plots depict conformational stability based on basins, where deeper basins (lower free energy) indicate more stable conformations. Compare the depth and spread of basins in each plot. A deeper, narrow basin often suggests a more stable and well-defined conformation. Check the free energy values of these basins. Generally, a lower minimum energy in one system’s FEL compared to the other suggests higher stability for that system. If one plot shows multiple shallow basins, it may indicate a flexible or less stable structure. A single, deep basin indicates the system consistently favors a particular stable conformation. This can differentiate a stable wild-type (or holo-protein) from a potentially less stable mutant (or apo-protein). For an apo vs. holo comparison, the holo form may exhibit a lower energy basin due to stabilizing interactions with the ligand, which can confirm the structural stabilization upon ligand binding.

[@pb_storageBhattacharya]

Thank you sir for the thorough explanation.🙏 I have two related queries: 1) By depth, do you mean the free energy value which may be plotted along the Z-axis in case of a 3D FEL? 2) if two plots have very close depths, but in one case, the low energy basin (let's assume denoted by deep blue colour code) is spread over a larger cluster of confirmations (larger area) while in the other, the deep blue zone it's narrow and involves a lesser number of clustered conformations (smaller area). Can we say that the one where the low energy basin is spread over a larger cluster is more stable since it involves greater number of stabilized conformations? Let me know if I am wrong. Thank you again for your time and valuable comments.

[@amisupriyo1]

Hello Pritam in my case sh script.sh command not working. IT show bad for loop variable. Please help

[@BioinfoCopilot]

First do chmod +x script.sh Then try to execute it from the parent folder.

[@noopur2604]

Hi Pritam Really insightful explaination! I have a question, sometimes the first two eignvalues are not insufficient to explain the the dynamics of the simlation and does not cover ~90% information on the structural changes happening in the protein. So, my queries are: 1) In this case how do we select the optimum number of eigenvectors 2) Suppose we select a the "-first" and "-last" as 1 and 5 in the gmx anaeig command and proceed with gmx sham, does PC1 and PC2 signify only first and fifth eigenvectors OR the PCA analyses all five eigenvectors and clusters them into two distinct PC1 and PC2?

[@BioinfoCopilot]

After running PCA, you should examine the eigenvalue spectrum (scree plot). This plot shows the eigenvalues in descending order, indicating how much variance each principal component (PC) explains. Look for an "elbow" in the plot, which suggests the point at which additional PCs start to contribute less significantly to the total variance. Check the cumulative variance explained by the eigenvectors. You can sum the variance explained by each successive PC until you reach ~90% of the total variance. This is a common threshold used in PCA to determine how many eigenvectors are sufficient to describe the system's dynamics. If two eigenvectors are not enough, you can increase the number of PCs accordingly. No, PC1 and PC2 do not refer to the first and fifth eigenvectors when you select -first 1 and -last 5 in the gmx anaeig command. gmx anaeig with -first 1 and -last 5: This command selects the first five eigenvectors (1 through 5) for analysis. It does not isolate just the first and fifth eigenvectors. Instead, you are including all five eigenvectors for downstream analysis, and these five are used to generate a reduced description of the system’s dynamics. In the context of gmx sham, when you perform Free Energy Landscape (FEL) analysis, PC1 and PC2 are the two principal components of the selected set of eigenvectors that you provide to the analysis. The PCA doesn't analyze just two specific eigenvectors (like the 1st and 5th); instead, it combines information from all eigenvectors you provided (in this case, eigenvectors 1 through 5) to construct two new PCs (PC1 and PC2) that best describe the overall variance in the dataset. Thanks!

[@noopur2604]

@@BioinfoCopilot Thank you so much for the explaination!! It cleared the concept a lot.

[@BioinfoCopilot]

@noopur2604 My pleasure

[@govind-nf7qg]

plz provide the commands for salt hbond and gmx rama

[@varuntk1183]

Hi, Can we extract the most stable comformation from FEL?

[@BioinfoCopilot]

There is no guarantee of this, and MD simulations are not driven by free energy, which is what is truly necessary to say one has the most “stable” conformation. The total potential energy in the MD simulation is not something that force fields are generally parametrized to ascribe any meaningful value, so there is no way to know if you are in a global energy minimum at any point during the trajectory.

[@varuntk1183]

@@BioinfoCopilot certain publication have done like extract the ensemble having lowest energy minima from fel plot, I was also not able to understand how it can be done

[@BioinfoCopilot]

Use trjconv to dump pdb files at several intervals. Then check the binding orientations. For each pdb file at certain interval check for free energy like I did for 0-100ns.

[@varuntk1183]

Thank you. Let me try it

[@MuhannadDehni]

thank you for great work, I notice from your result Gibbes energy landscape start from 0 Kj/mol. and I try to do that I get the same ,why we didn't get negative value even MMPBSA is = -9 kcal/mol

[@BioinfoCopilot]

The tutorial shown was for demonstration purposes. It might vary from specifics and your simulation environment.