Calculating Torsional Angles

[Christien Williams]

I am trying to calculate the torsional angles across that SSBONDs. When you use:

tor <- torsion.pdb(pdb) it gives you the following result:

Are the values on the far left column the residue/amino acid number in the chain?

[Barry Grant]

There should be one row per residue that torsions could be calculated for. E.g.

> library(bio3d)

> pdb <- read.pdb("5p21")

# Trim to protein only - should be 166 residues (see below)
> pdb.prot <- trim.pdb(pdb, "protein")
> pdb.prot

 Call:  trim.pdb(pdb = pdb, "protein")

   Total Models#: 1
     Total Atoms#: 1323,  XYZs#: 3969  Chains#: 1  (values: A)

     Protein Atoms#: 1323  (residues/Calpha atoms#: 166)
     Nucleic acid Atoms#: 0  (residues/phosphate atoms#: 0)

     Non-protein/nucleic Atoms#: 0  (residues: 0)
     Non-protein/nucleic resid values: [ none ]

   Protein sequence:
      MTEYKLVVVGAGGVGKSALTIQLIQNHFVDEYDPTIEDSYRKQVVIDGETCLLDILDTAG
      QEEYSAMRDQYMRTGEGFLCVFAINNTKSFEDIHQYREQIKRVKDSDDVPMVLVGNKCDL
      AARTVESRQAQDLARSYGIPYIETSAKTRQGVEDAFYTLVREIRQH

+ attr: atom, helix, sheet, seqres, xyz,
        calpha, call

# Calculate our torsions...
> tor.prot <- torsion.pdb(pdb.prot)

# Should be 166 rows, one per residue...
> dim(tor.prot$tbl)
[1] 166   7

# You could label these rows by original chain and residue numbers from your input PDB
> my_residue_labels <- paste(pdb.prot$atom$chain[pdb.prot$calpha], pdb.prot$atom$resno[pdb.prot$calpha], sep="-")
> head(my_residue_labels)
[1] "A-1" "A-2" "A-3" "A-4" "A-5" "A-6"

> rownames(tor.prot$tbl) <-  my_residue_labels

# Now we can more easily look up the angles for e.g. residue 152 from chain A
> tor.prot$tbl["A-152",]
      phi       psi      chi1      chi2      chi3      chi4      chi5
-65.18211 -52.28821 171.36372        NA        NA        NA        NA

Hope this helps! Barry

[Lars Skjærven]

Hi Christien, That's the row number of the matrix ... (check out matrices and data frames in R). No row names have been assigned, but you can do that yourself using rownames

> rownames(tor$tbl) = pdb$atom$resno[ pdb$calpha ]
> head(tor$tbl)
         phi       psi       chi1      chi2      chi3      chi4 chi5
1         NA 133.04002  171.92390  170.1749 175.68456 -171.6513   NA
2 -100.27840 105.17020  178.20639        NA        NA        NA   NA
3  -74.09361 155.83103  -79.79702   78.0313        NA        NA   NA
4  -87.24712 148.52165         NA        NA        NA        NA   NA
5  -45.32540 -56.26063 -172.93132 -176.6713  53.87713   92.3284   NA
6  -72.02024 -34.59925  -64.30896        NA        NA        NA   NA

[Christien Williams]

Thanks to both of you! This helps a lot.

[Christien Williams]

Hi Mr. Grant,

Thank you again for your insight on stripping a pdb file. I’ve been using this code from another answer on the forum:

[1] 230 231 232 [1] "SSBOND 1 CYS A 3 CYS A 40 1555 1555 2.03 " [2] "SSBOND 2 CYS A 4 CYS A 32 1555 1555 2.10 " [3] "SSBOND 3 CYS A 16 CYS A 26 1555 1555 2.09 "

In addition to that code you provided me with in order to strip a PDB file so each row is a single residue, and then using the information above to go to a specific cysteine residue that has a disulfide bond. The ultimate goal is to write a script that that creates a list of the chi1, chi2, and chi3 values for the disulfide bridges of hundreds of proteins. As of right now, I have code that can strip a PDB file, and then prints out the torsional angles of what supposedly is the cysteine that forms a disulfide bridge (for example, per the code above, it goes to “A-308” and takes the torsional angles). The problem is, the result shows that there is no chi2, or chi3 values. However, since the residue at A-308 is a cysteine that has a disulfide bond, both of those chi values should be present. Do you have any idea why this might be happening?

Thank you,

Christien Williams

[Barry Grant]

Best send/upload your example PDB along with the complete snippet of code displaying the problem as a quick test here with a regular PDB file worked fine for me.

[Christien Williams]

Hi Mr. Grant,

A sample protein being used is 4u72 (http://files.rcsb.org/header/4U72.pdb). Here is the code I’ve written:

• {
• i = unlist(strsplit(i, split='.', fixed=TRUE))[1]
• pdbfilesnew <- c(pdbfilesnew, i)
• } Note: Accessing on-line PDB file HEADER OXIDOREDUCTASE 30-JUL-14 4U72 PDB has ALT records, taking A only, rm.alt=TRUE Note: Accessing on-line PDB file HEADER OXIDOREDUCTASE 30-JUL-14 4U72 PDB has ALT records, taking A only, rm.alt=TRUE

[1] "A" "308" [1] "B" "308" phi psi chi1 chi2 chi3 chi4 chi5 -57.29319 -40.69293 168.11632 NA NA NA NA phi psi chi1 chi2 chi3 chi4 chi5

[Christien Williams]

Hi Mr. Grant,

A sample protein being used is 4u72 (http://files.rcsb.org/header/4U72.pdb). Here is the code I’ve written:

• {
• i = unlist(strsplit(i, split='.', fixed=TRUE))[1]
• pdbfilesnew <- c(pdbfilesnew, i)
• } Note: Accessing on-line PDB file HEADER OXIDOREDUCTASE 30-JUL-14 4U72 PDB has ALT records, taking A only, rm.alt=TRUE Note: Accessing on-line PDB file HEADER OXIDOREDUCTASE 30-JUL-14 4U72 PDB has ALT records, taking A only, rm.alt=TRUE

[1] "A" "308" [1] "B" "308" phi psi chi1 chi2 chi3 chi4 chi5 -57.29319 -40.69293 168.11632 NA NA NA NA phi psi chi1 chi2 chi3 chi4 chi5

[Lars Skjærven]

Hi, Chi2 is not calculated for cysteins. Are you looking for CA-CB-SG-SG ? where the latter SG belongs to the opposing residue? You'll have to do this manually (torsion.pdb() does not know that Cys308 contains a disulfide bridge).

# atom selections
> sele1 = atom.select(pdb.prot, elety=c("CA", "CB", "SG"), resno=308, chain="A")
> sele2 = atom.select(pdb.prot, elety=c("SG"), resno=308, chain="B")

# combine xyz indices obtained from the two above selection
> inds.xyz = c(sele1$xyz, sele2$xyz)

# calculate angle
> torsion.xyz(pdb.prot$xyz[, inds.xyz])
[1] 84.08803

[Christien Williams]

Thank you so much Lars! I’m getting really close to reaching my goal. However, I have another question. I used this code and calculated the torsional angles for chi1, chi2, and chi3 angles. The research group I’m working with has also measured chi angles by hand (using WinCoot) to verify the numbers the code is getting. I am getting solid results for the chi1 and chi2 angles using your method, however, for some reason the chi3 values are incorrect. I believe that for a disulfide bond/bridge chi3 values should include the beta carbon and sulfur from both of the chains. So:

sele1 = atom.select(pdb.prot, elety=c("CB", "SG"), resno=308, chain="A")

sele2 = atom.select(pdb.prot, elety=c("SG","CB"), resno=308, chain="B")

Does this look like the right idea?

Thanks again for all of your help in this process

[Lars Skjærven]

Looks right. I guess also the order of the indices here might be important, but you are on the right track. Good luck!