On Aug 30, 2014, at 4:00 AM, Aswani Kumar Kancherla <ashwin2952@gmail.com> wrote:

@Eric:​
​Many thanks for the code that you gave for extracting disulfide bonds and measuring the chiSS. I incorporated it in the code from the link sent by Elaine and executed it on a list of pdb ids. ​It worked fine on most of pdb files. However, for a few pdb files, some of the disulfide bonds are not recognized. Could it be due to the bad geometry in those structures?

If I understand the code right, it  extracts pairs of SG atoms that are bonded together in each model "m" by looking in "m.bonds". Then it goes on to identify the attached beta carbons and then finally calculates the chiSS. 

Yep.

My suspicion is that, in some of the NMR structures (the pdb ids on which the script failed to identify all the disulfide bonds), the disulfide bond geometry is not good and hence SS bond is not listed in "m.bonds". 

Chimera relies on information inside the PDB file itself to "identify" SS bonds.  In particular, LINK and/or CONECT records must be present that identify the SS linkage.

How do I go about, if I want to identify all the pairs of Cysteine SG atoms which are closer than say 2.5A and then calculate the appropriate dihedral angle about CB-SG-SG'-CB'? I will first need to calculate pairwise distances for SG atoms of all cysteines, select those pairs which are much nearer than allowed by their non-bonded distance (sum of van der waal radii) and then calculate the dihedral angle. This exercise will be useful for me to identify pairs of those Cysteines whose SG atoms have optimal geometry for formation of disulfide bond.

Once again many thanks for your valuable time and prompt help. I am pasting the script that I have used below. It worked fine and hence the indentation must be fine.

Well, if there are only a few of these "misses" (and you're certain you've caught all the misses), it might just be easier to add appropriate LINK/CONECT records to the files, or use Chimera's "bond" command to form the missing SS bond.  Code to do what you suggested above would be something like this:

from chimera import openModels, Molecule
for m in openModels.list(modelTypes=[Molecule]):
sgs = [a for a in m.atoms if a.name == "SG"]
for i, sg1 in enumerate(sgs):
for sg2 in sgs[i+1:]:
if sg1.coord().distance(sg2.coord()) > 2.5:
continue
# the code below should look familiar :-)
cbs = []
for sg in [sg1, sg2]:
if len(sg.neighbors) > 2:
break
for nb in sg.neighbors:
if nb.name == "CB":
cbs.append(nb)
break
else:
# not bonded to a CB
break
if len(cbs) != 2:
continue
cb1, cb2 =  cbs
from chimera import dihedral
print cb1, sg1, sg2, cb2, dihedral(cb1.coord(), sg1.coord(), sg2.coord(), cb2.coord())

--Eric


                        Eric Pettersen

                        UCSF Computer Graphics Lab

                        http://www.cgl.ucsf.edu