clauswilke / PeptideBuilder

A simple Python library to generate model peptides
MIT License
79 stars 30 forks source link

Single-residue structure: angles have zero effect? #49

Open bwllc opened 1 year ago

bwllc commented 1 year ago

Hello,

I think that I understand that the psi_im1 angle will have no effect on the positions of the (non-hydrogen) atoms in residue 1. I believe that's why the angle passed to the PeptideBuilder is called psi_im1, and not psi.

However, I was expecting that the phi angle would cause a rotation around the CA - C bond, and therefore the peptide oxygen atom would move.

This is not the result I am obtaining. Here is minimal example code:

import numpy as np
import PeptideBuilder
from PeptideBuilder import Geometry

np.set_printoptions(precision=2, suppress=True, floatmode="fixed", sign=" ")

for angles in ((-140, 130), (-165, 130), (-140, 155)):
    geo = Geometry.geometry("A")
    geo.phi, geo.psi_im1 = angles
    structure = PeptideBuilder.initialize_res(geo)
    print(f"(phi, psi_im1) = ({geo.phi}, {geo.psi_im1})")
    names = [f"{atom.name:<3}" for atom in structure.get_atoms()]
    coords = [str(atom.coord) for atom in structure.get_atoms()]
    for n, c in zip(names, coords):
        print(n + c)
    print()

Here is the output:

(phi, psi_im1) = (-140, 130)
N  [-0.52  1.36  0.00]
CA [ 0.00  0.00  0.00]
C  [ 1.52  0.00  0.00]
O  [ 2.14  0.52 -0.92]
CB [-0.51 -0.77 -1.21]

(phi, psi_im1) = (-165, 130)
N  [-0.52  1.36  0.00]
CA [ 0.00  0.00  0.00]
C  [ 1.52  0.00  0.00]
O  [ 2.14  0.52 -0.92]
CB [-0.51 -0.77 -1.21]

(phi, psi_im1) = (-140, 155)
N  [-0.52  1.36  0.00]
CA [ 0.00  0.00  0.00]
C  [ 1.52  0.00  0.00]
O  [ 2.14  0.52 -0.92]
CB [-0.51 -0.77 -1.21]

Please help me to understand why the coordinates of all the atoms are the same, regardless of the backbone angles that were selected. Thanks!

ausmeyer commented 1 year ago

I am a bit rusty on this, but hopefully it is useful and not wrong.

This is probably something Claus or Matthew will have to address from a coding perspective, but from a biophysical perspective it doesn't really make sense (I think) to define a psi (or phi obviously) angle without the next nitrogen. I understand we are making a model here, but that would be like a biophysically impossible carbon coordinate.

The angle that should move the first oxygen and second nitrogen in space should be psi, I believe. The phi should move the second C and CB. Here is what happens when you add a second residue and change phi:

for angles in ((-140, 130), (-165, 130)):
    geo = Geometry.geometry("A")
    geo.phi, geo.psi_im1 = (-140, 130)
    structure = PeptideBuilder.initialize_res(geo)
    geo.phi, geo.psi_im1 = angles
    PeptideBuilder.add_residue(structure, geo)
    print(f"(phi, psi_im1) = ({geo.phi}, {geo.psi_im1})")
    names = [f"{atom.name:<3}" for atom in structure.get_atoms()]
    coords = [str(atom.coord) for atom in structure.get_atoms()]
    for n, c in zip(names, coords):
        print(n + c)
    print()

Output:

(phi, psi_im1) = (-140, 130)
N  [-0.52  1.36  0.00]
CA [ 0.00  0.00  0.00]
C  [ 1.52  0.00  0.00]
O  [ 2.14  0.68 -0.81]
CB [-0.51 -0.77 -1.21]
N  [ 2.12 -0.76  0.91]
CA [ 3.57 -0.84  1.00]
C  [ 4.03 -2.27  1.28]
O  [ 5.22 -2.53  1.39]
CB [ 4.09  0.08  2.09]

(phi, psi_im1) = (-165, 130)
N  [-0.52  1.36  0.00]
CA [ 0.00  0.00  0.00]
C  [ 1.52  0.00  0.00]
O  [ 2.14  0.68 -0.81]
CB [-0.51 -0.77 -1.21]
N  [ 2.12 -0.76  0.91]
CA [ 3.57 -0.84  1.00]
C  [ 4.00 -2.03  1.85]
O  [ 5.20 -2.26  2.04]
CB [ 4.14  0.44  1.60]

It looks like to me the first dihedral rotating about the second N-C (CB and C are moving) which is phi.

Then here is what happens when you add a second residue and change psi:

for angles in ((-140, 130), (-140, 155)):
    geo = Geometry.geometry("A")
    geo.phi, geo.psi_im1 = (-140, 130)
    structure = PeptideBuilder.initialize_res(geo)
    geo.phi, geo.psi_im1 = angles
    PeptideBuilder.add_residue(structure, geo)
    print(f"(phi, psi_im1) = ({geo.phi}, {geo.psi_im1})")
    names = [f"{atom.name:<3}" for atom in structure.get_atoms()]
    coords = [str(atom.coord) for atom in structure.get_atoms()]
    for n, c in zip(names, coords):
        print(n + c)
    print()

Output:

(phi, psi_im1) = (-140, 130)
N  [-0.52  1.36  0.00]
CA [ 0.00  0.00  0.00]
C  [ 1.52  0.00  0.00]
O  [ 2.14  0.68 -0.81]
CB [-0.51 -0.77 -1.21]
N  [ 2.12 -0.76  0.91]
CA [ 3.57 -0.84  1.00]
C  [ 4.03 -2.27  1.28]
O  [ 5.22 -2.53  1.39]
CB [ 4.09  0.08  2.09]

(phi, psi_im1) = (-140, 155)
N  [-0.52  1.36  0.00]
CA [ 0.00  0.00  0.00]
C  [ 1.52  0.00  0.00]
O  [ 2.14  0.96 -0.45]
CB [-0.51 -0.77 -1.21]
N  [ 2.12 -1.08  0.50]
CA [ 3.57 -1.19  0.55]
C  [ 4.03 -2.59  0.20]
O  [ 5.22 -2.88  0.19]
CB [ 4.09 -0.81  1.93]

It looks like the first dihedral rotating is about the CA-CB (the first oxygen moving) which is psi.

It just looks to me like there was a decision made to not adjust the oxygen in space until the second amino acid was added. I didn't make that decision, but it seems reasonable.

ausmeyer commented 1 year ago

As an aside, if you need the psi on the single amino acid, it would be an easy work around to build the dipeptide and then erase everything in the output pdb file after CB. They are just text files after all. Maybe someone else will know how to do it without the work around.

bwllc commented 1 year ago

Thank you very much for your quick replies, @ausmeyer.

I need to study my peptide backbone dihedral angle drawings more closely, but I still think that the oxygen in residue 1 should move if you change phi. Or, if the oxygen does not move, everything attached to CA should spin instead. I'm trying to build a beta sheet geometry from scratch. I need to get the C=O bond pointed at one neighbor strand, while getting the N atom pointed at the other neighbor strand..

You think like a hacker! I had already considered constructing a dipeptide and ignoring residue 1 as a workaround. (I don't have to work with the text file, I also know how to manipulate structural data in Biopython.) I just thought that I was misusing PeptideBuilder somehow. Apparently, I might be using it exactly as it was designed.

I hope we'll hear a bit more from your colleagues regarding this issue.

clauswilke commented 1 year ago

Just wanted to chime in that I have no idea. This is a question for @mtien.

ausmeyer commented 1 year ago

Why would we expect the first oxygen to move when phi is changed?

Phi is angle between the Nitrogen and the C-alpha so the phi of the first residue is symmetric on the nitrogen side and therefore undefined. A phi wouldn’t be defined until the second residue is added.

Psi is the angle between the C-alpha and CB so that would rotate the oxygen around in space when the CB coordination number is correct by adding a second residue.

Both of those expectations play out in the dipeptide above.

bwllc commented 1 year ago

Thanks for the continued replies.

I've looked deeper into the source code.

I can confirm that PeptideBuilder.initialize_res() does not make use of the phi and psi_im1 values of the Geometry that you pass into it. You get the default configuration. In contrast, PeptideBuilder.add_residue() calls PeptideBuilder.add_residue(), and in the latter function we can see that the phi and psi_im1 values are passed to calculateCoordinates(). This explains why the initial residue always has the same backbone orientation.

I'm still of the opinion that there should be at least one degree of freedom in the position of the backbone atoms, even when there is only a single residue, provided that we also have coordinates where we want to place three non-co-linear atoms at known locations: for example, N, Cα and O. That is the situation I have, and I think it is a typical use case.

This might be a misunderstanding on my part, but consider this diagram:

https://en.wikipedia.org/wiki/Ramachandran_plot#/media/File:Protein_backbone_PhiPsiOmega_drawing.svg

If N and Cα are fixed in space, and C(-1) exists and defines the reference angle for ϕ, and I vary ϕ, I expect to see Cβ, C and O move in circles around the N - Cα bond axis.

Now, C(-1) does not exist with the initial residue, but its location is implied by the coordinate system constraints that I described above.

You have suggested a workaround that I was already thinking about. Having an (i-1) residue creates a C(i-1) atom, and so, imposes the same kinds of constraints that I'm describing. I'll try this approach and report back. I won't mark this issue as closed for now.

If anyone wants to continue discussing my other points about specifying protein backbone geometry, we can continue to do that here, or we can move the discussion to another location that people think is more appropriate.

ausmeyer commented 1 year ago

This might be a misunderstanding on my part, but consider this diagram:

https://en.wikipedia.org/wiki/Ramachandran_plot#/media/File:Protein_backbone_PhiPsiOmega_drawing.svg

If N and Cα are fixed in space, and C(-1) exists and defines the reference angle for ϕ, and I vary ϕ, I expect to see Cβ, C and O move in circles around the N - Cα bond axis.

Now, C(-1) does not exist with the initial residue, but its location is implied by the coordinate system constraints that I described above.

This is a misunderstanding. The dihedral angles specify rotation about a bond so biophysically they require 4 backbone atoms to be calculated. Imagine that you are looking down the axis of the N-Ca bond such that N can't be seen because it is immediately behind Ca. Then, you would see Cb sticking off in some direction. No matter what direction it points, the phi is undefined because it requires Cb(-1) to define the coordinate system.

Here is a specified phi the A as Cb(-1), B as N, C as Ca, D as Cb. The angle between A and D on the right is phi.

dihed_specified

Here is an unspecified phi due to missing Cb(-1). Notice, it doesn't matter that B and C are fixed in a coordinate system, you still can't define the position of D with a dihedral angle because you don't have the second half-plane defined.

dihed_unspecified

See the definition of dihedral angles for further details.

https://en.wikipedia.org/wiki/Dihedral_angle