On my two data files c7ax_FINAL.pdb and c7eq_FINAL.pdb I wish to perform a root-mean-square-deviation calculation:
My two files are:
< ATOM 1 C ALA A 1 3.953 7.886 1.035 0.00 1.00 ALA
< ATOM 2 O ALA A 1 3.516 8.953 0.733 0.00 0.00 ALA
< ATOM 3 NT ALA A 1 5.207 7.558 0.961 0.00 1.00 ALA
< ATOM 4 HNT ALA A 1 5.464 6.607 1.102 0.00 0.00 ALA
< ATOM 5 CAT ALA A 1 6.224 8.463 0.408 0.00 0.00 ALA
< ATOM 6 HT1 ALA A 1 6.081 8.684 -0.723 0.00 0.00 ALA
< ATOM 7 HT2 ALA A 1 7.281 8.023 0.476 0.00 0.00 ALA
< ATOM 8 HT3 ALA A 1 6.242 9.421 0.971 0.00 0.00 ALA
< ATOM 9 CAY ALA A 1 3.829 3.665 -0.589 0.00 1.00 ALA
< ATOM 10 HY1 ALA A 1 4.149 4.117 -1.530 0.00 0.00 ALA
< ATOM 11 HY2 ALA A 1 4.389 2.789 -0.273 0.00 0.00 ALA
< ATOM 12 HY3 ALA A 1 2.715 3.393 -0.645 0.00 0.00 ALA
< ATOM 13 CY ALA A 1 3.917 4.701 0.520 0.00 0.00 ALA
< ATOM 14 OY ALA A 1 4.828 4.638 1.301 0.00 1.00 ALA
< ATOM 15 N ALA A 1 2.980 5.688 0.542 0.00 0.00 ALA
< ATOM 16 HN ALA A 1 2.165 5.607 0.046 0.00 1.00 ALA
< ATOM 17 CA ALA A 1 3.006 6.799 1.485 0.00 1.00 ALA
< ATOM 18 HA ALA A 1 2.073 7.253 1.419 0.00 0.00 ALA
< ATOM 19 CB ALA A 1 3.241 6.464 2.973 0.00 0.00 ALA
< ATOM 20 HB1 ALA A 1 3.464 7.440 3.349 0.00 0.00 ALA
< ATOM 21 HB2 ALA A 1 4.125 5.884 3.182 0.00 0.00 ALA
< ATOM 22 HB3 ALA A 1 2.252 6.095 3.299 0.00 0.00 ALA
---
> ATOM 1 C ALA A 1 3.976 7.842 1.236 0.00 1.00 ALA
> ATOM 2 O ALA A 1 4.013 8.378 2.325 0.00 0.00 ALA
> ATOM 3 NT ALA A 1 5.027 7.787 0.418 0.00 1.00 ALA
> ATOM 4 HNT ALA A 1 4.956 7.219 -0.411 0.00 0.00 ALA
> ATOM 5 CAT ALA A 1 6.221 8.583 0.554 0.00 0.00 ALA
> ATOM 6 HT1 ALA A 1 6.061 9.618 0.811 0.00 0.00 ALA
> ATOM 7 HT2 ALA A 1 6.891 8.420 -0.307 0.00 0.00 ALA
> ATOM 8 HT3 ALA A 1 6.814 8.138 1.406 0.00 0.00 ALA
> ATOM 9 CAY ALA A 1 3.787 3.479 -0.092 0.00 1.00 ALA
> ATOM 10 HY1 ALA A 1 4.869 3.206 0.041 0.00 0.00 ALA
> ATOM 11 HY2 ALA A 1 3.160 3.112 0.762 0.00 0.00 ALA
> ATOM 12 HY3 ALA A 1 3.421 2.957 -1.005 0.00 0.00 ALA
> ATOM 13 CY ALA A 1 3.693 4.933 -0.222 0.00 0.00 ALA
> ATOM 14 OY ALA A 1 4.325 5.508 -1.108 0.00 1.00 ALA
> ATOM 15 N ALA A 1 2.926 5.567 0.715 0.00 0.00 ALA
> ATOM 16 HN ALA A 1 2.642 5.047 1.543 0.00 1.00 ALA
> ATOM 17 CA ALA A 1 2.783 7.045 0.789 0.00 1.00 ALA
> ATOM 18 HA ALA A 1 2.453 7.408 -0.216 0.00 0.00 ALA
> ATOM 19 CB ALA A 1 1.593 7.521 1.595 0.00 0.00 ALA
> ATOM 20 HB1 ALA A 1 1.687 7.108 2.610 0.00 0.00 ALA
> ATOM 21 HB2 ALA A 1 0.653 7.196 1.126 0.00 0.00 ALA
> ATOM 22 HB3 ALA A 1 1.490 8.595 1.523 0.00 0.00 ALA
The run command is:
match.py -nr c7ax_FINAL.pdb c7eq_FINAL.pdb "[('A:1', 'A:1')]"
I get the following error:
Aligning CA atoms of residues: A:1
No rotations
Traceback (most recent call last):
File "Documents/RMSD_two_PDB/match.py", line 221, in <module>
rmsd = get_raw_rmsd(pdb1, pdb2, segments1, segments2, ['CA'])
File "Documents/RMSD_two_PDB/match.py", line 118, in get_raw_rmsd
return sum_rmsd(atoms1, atoms2)
File "Documents/RMSD_two_PDB/match.py", line 110, in sum_rmsd
return math.sqrt(sum_squared/float(len(atoms1)))
ZeroDivisionError: float division
The program is
#!/usr/bin/env python
"""
Set of routines to calculate the RMSD between two molecular structures.
The module can be run from the command line using PDB files as input.
"""
import math
import numpy
import vector3d, util, molecule, polymer
def rmsd(crds1, crds2):
"""Returns RMSD between 2 sets of [nx3] numpy array"""
assert(crds1.shape[1] == 3)
assert(crds1.shape == crds2.shape)
n_vec = numpy.shape(crds1)[0]
correlation_matrix = numpy.dot(numpy.transpose(crds1), crds2)
v, s, w = numpy.linalg.svd(correlation_matrix)
is_reflection = (numpy.linalg.det(v) * numpy.linalg.det(w)) < 0.0
if is_reflection:
s[-1] = - s[-1]
E0 = sum(sum(crds1 * crds1)) + \
sum(sum(crds2 * crds2))
rmsd_sq = (E0 - 2.0*sum(s)) / float(n_vec)
rmsd_sq = max([rmsd_sq, 0.0])
return numpy.sqrt(rmsd_sq)
def optimal_superposition(crds1, crds2):
"""Returns best-fit rotation matrix as [3x3] numpy matrix"""
assert(crds1.shape[1] == 3)
assert(crds1.shape == crds2.shape)
correlation_matrix = numpy.dot(numpy.transpose(crds1), crds2)
v, s, w = numpy.linalg.svd(correlation_matrix)
is_reflection = (numpy.linalg.det(v) * numpy.linalg.det(w)) < 0.0
if is_reflection:
v[-1,:] = -v[-1,:]
return numpy.dot(v, w)
def get_i_residue(residues, tag):
def get_tag(residue):
tag = ""
if residue.chain_id != " " and residue.chain_id != "":
tag += residue.chain_id + ":"
tag += str(residue.num)
if residue.insert:
tag += residue.insert
return tag
# clean up tag
tag = tag.strip()
if tag[0] == ":":
tag = tag[1:]
if not tag[0].isdigit() and tag[1].isdigit():
tag = tag[0] + ":" + tag[1:]
for i, residue in enumerate(residues):
if tag.lower() == get_tag(residue).lower():
return i
raise "Can't find residue", tag
def get_superposable_atoms(polymer, segments,
atom_types=['CA', 'N', 'C', 'CB']):
result = []
allowed_i = []
residues = polymer.residues()
for res_num_i, res_num_j in segments:
i = get_i_residue(residues, str(res_num_i))
j = get_i_residue(residues, str(res_num_j))
allowed_i.extend(range(i,j))
for i, residue in enumerate(residues):
if i in allowed_i:
result.extend([a for a in residue.atoms()
if a.type in atom_types])
return result
def get_crds(atoms):
crds = numpy.zeros((len(atoms), 3), float)
for i, a in enumerate(atoms):
crds[i,0] = a.pos.x
crds[i,1] = a.pos.y
crds[i,2] = a.pos.z
return crds
def calculate_superposition_matrix(atoms1, atoms2):
def convert_to_matrix3d(numpy_matrix3d):
result = vector3d.Matrix3d()
for i in range(3):
for j in range(3):
result.setElem(i, j, numpy_rotation[j, i])
return result
numpy_rotation = optimal_superposition(get_crds(atoms1), get_crds(atoms2))
return convert_to_matrix3d(numpy_rotation)
def sum_rmsd(atoms1, atoms2):
sum_squared = 0.0
for atom1, atom2 in zip(atoms1, atoms2):
sum_squared += vector3d.pos_distance(atom1.pos, atom2.pos)**2
return math.sqrt(sum_squared/float(len(atoms1)))
def get_raw_rmsd(pdb1, pdb2, segments1, segments2, atom_types):
polymer1 = polymer.Polymer(pdb1)
polymer2 = polymer.Polymer(pdb2)
atoms1 = get_superposable_atoms(polymer1, segments1, atom_types)
atoms2 = get_superposable_atoms(polymer2, segments2, atom_types)
return sum_rmsd(atoms1, atoms2)
def get_best_alignment(pdb1, pdb2, segments1, segments2, atom_types):
"""Returns rmsd and filename of transformed pdb2."""
polymer1 = polymer.Polymer(pdb1)
atoms1 = get_superposable_atoms(polymer1, segments1, atom_types)
polymer2 = polymer.Polymer(pdb2)
atoms2 = get_superposable_atoms(polymer2, segments2, atom_types)
center1 = molecule.get_center(atoms1)
polymer1.transform(vector3d.Translation(-center1))
polymer2.transform(vector3d.Translation(-molecule.get_center(atoms2)))
polymer2.transform(calculate_superposition_matrix(atoms1, atoms2))
rmsd = sum_rmsd(atoms1, atoms2)
temp_pdb2 = util.fname_variant(pdb2)
polymer2.transform(vector3d.Translation(center1))
polymer2.write_pdb(temp_pdb2)
return rmsd, temp_pdb2
def get_rmsd(pdb1, pdb2, segments1, segments2, atom_types):
polymer1 = polymer.Polymer(pdb1)
atoms1 = get_superposable_atoms(polymer1, segments1, atom_types)
polymer2 = polymer.Polymer(pdb2)
atoms2 = get_superposable_atoms(polymer2, segments2, atom_types)
center1 = molecule.get_center(atoms1)
polymer1.transform(vector3d.Translation(-center1))
polymer2.transform(vector3d.Translation(-molecule.get_center(atoms2)))
crds1 = get_crds(atoms1)
crds2 = get_crds(atoms2)
return rmsd(crds1, crds2)
def segments_str(segments):
residues = []
for i, j in segments:
if i == j:
residues.append(str(i))
else:
residues.append("%s-%s" % (i,j))
return ', '.join(residues)
Any ideas? It seems that float(len(atoms1))=0, i.e. the length of the string?