Fixed N-C Peptide Bond Breakage within OpenMM Minimization

[PowersPope]

Whenever amber relaxation was applied to predicted structures, the N-C bond was broken. I manually added in the peptide bond and remove the hydrogen and OXT atoms before applying it to the system.
Also, made some small changes to the README.md to make it a little easier for people to understand how to install the repo.

[baba-hashimoto]

Thank you for your consideration. When I tried the program, I got the following error message: “No template found for residue 30 (PRO). The set of atoms matches CPRO, but the bonds are different. Perhaps the cha in is missing a terminal group? Pro is the last amino acid residue. It seems that I cannot create the input file for the cyclic structure properly.

[PowersPope]

@baba-hashimoto thanks for the info. I am guessing you are predicting the structure of a target and macrocycle together? When I fixed this, I was doing it based on a macrocycle by itself in solution. Though I can fix it for you.
Also, I believe I need to fix my solution as I didn't test it against various peptide sizes and amino acid compositions. Give me a day or two.

[PowersPope]

@baba-hashimoto I just pushed a new fix. Let me know if that fixes your problem. 😄
Made the N-C and atom deletions modular to the input size. Also, made it not atom position dependent, but instead dependent on atom name. This fix did not give me any problems when testing with 8mers or 10mers. However, I didn't test it thoroughly and would like feedback to know if this particular minimize problem persists.

This was only tested for Monomers, don't know how well it functions for Complexes.

[PowersPope]

@baba-hashimoto I fixed the multimer as well now, so that outputs have N-C bonds after minimization. This can all be accessed by using my Repo until they merge this. Please let me now if there are any errors or problems with the minimization.

[baba-hashimoto]

@PowersPope
First of all, I would like to express my gratitude. Thanks to your quick response, I was able to run HighFold to energy minimization without any problems.　It is also working correctly with the complex.
I had no experience using OpenMM, so I was unable to do anything, but thanks to your precise corrections, I was able to run HighFold.

I am afraid to comment on your excellent corrections, but I would like to make two requests for corrections.

･I would like you to add the CONECT information to the end of the relaxed structure PDB file. The number itself should be obtained somewhere in your program, so it would be helpful if you could add it after the energy minimization is complete.

･I would like you to set a threshold value to allow formation before creating the n-c bond.

fasta
SAKIDNLD:
SSPGIWLDCTHLEGKVILVAVHVASGYIEAVIPAETGQETAYFLLLAGRWPVKTHDNGSNFTSTTVKAACWWAGIQEDGIPYNPQSQGVIESMNKELKKIIGQVRDQAEHLKTAVQMAVFIHNHKRKGYSAGERIVIIATDIQ:
SSPGIWLDCTHLEGKVILVAVHVASGYIEAVIPAETGQETAYFLLLAGRWPVKTHDNGSNFTSTTVKAACWWAGIQEDGIPYNPQSQGVIESMNKELKKIIGQVRDQAEHLKTAVQMAVFIHNHKRKGYSAGERIVIIATDIQ

The above is the sequence given as an example by the first author. Please compare the structures of the unrelaxed and relaxed peptides. As far as I have been able to determine, the distance between the ends is quite far apart. I cannot determine whether joining them together will lead to the correct result.

On the other hand, if you run the prediction using HighFold_data/monomer_native/test_fasta/1bh4.fasta, the distance between the two ends is predicted to be close (maybe it's just that the ends are set close together because they are in the training data. I'm sorry if this is a wrong point because I haven't been able to check the training data.

If possible, I think that setting a threshold or alerting the user when the N-C distance is too far would convey the correct information to the user.

The above comment is not related to the essential discussion of the program, so you can ignore it if you like.

[PowersPope]

@baba-hashimoto I am glad that it is functioning for you now!

･I would like you to add the CONECT information to the end of the relaxed structure PDB file. The number itself should be obtained somewhere in your program, so it would be helpful if you could add it after the energy minimization is complete.

I looked through the code and had a hard time trying to find where exactly the PDB writing function/portion is in this script. I found the GoogleDeep Mind portion they are using in alphafold/common/protein.py. However, when I changed stuff there it would only add stuff to the unrelaxed outputs. I could add a print statement during minimization that specifies where I am manually placing the amide bond. However, it may repeat it a bunch I haven't tried it yet.

･I would like you to set a threshold value to allow formation before creating the n-c bond.

I regards to this I believe that you are using --model-type alphafold2 which is fine for monomers. However, if you want to use the fasta input you mentioned above you would need to use one of the multimer models like --model-type alphafold2_multimer_v3 which worked fine for me and the macrocycle was not linear in the unrelaxed output.

[baba-hashimoto]

@PowersPope
Thank you for taking the time to consider my comments.
Regarding the first point, thank you for investigating. I also considered whether it would be possible to add CONECT information in the post-processing. I have attached the file.
makeCONECT.zip

Regarding the second point,
--model-type alphafold2_multimer_v3
I was not aware that it was used.
To run the prediction, type
colabfold_batch --templates --amber --model-type alphafold2 file_input path_output [args]
It was written as above, so I mistakenly recognized it.