example.py

from Bio.PDB import *


def extract_torsion_angles(structure):
	ic_chain = internal_coords.IC_Chain(structure)
	ic_chain.atom_to_internal_coordinates()

	d = ic_chain.dihedra

	"""
	Update the keys of a dihedra dictionary to make searching easier.

	The keys of the dihedra dictionary (as obtained from
	extract_torsion_angles()) are 4-tuples where each element is
	an AtomKey containing 5 fields (respos, icode, resname, atm,
	altloc, occ). We use the residue position, residue name and
	atom name in order to create a new representation for each
	AtomKey. The final dihedra dictionary will have 4-tuple keys
	where each element of the tuple is one such string.
	"""

	new_d = {}
	for key, value in sorted(d.items(), key=lambda x: x[0]):
		new_key = []
		for atom_key in key:
			# Keep: residue position + resiude name + atom name.
			atom_key_str = (atom_key.akl[0] + '_' +
				atom_key.akl[2] + '_' + atom_key.akl[3])

			new_key.append(atom_key_str)

		new_d[tuple(new_key)] = value

	return new_d


f1 = 'protein1.pdb'
f2 = 'protein2.pdb'

parser = PDBParser(PERMISSIVE=1, QUIET=0)
s1 = parser.get_structure('', f1)
s2 = parser.get_structure('', f2)

c1 = list(s1.get_chains())
c2 = list(s2.get_chains())

print('Chains: ', len(c1), ' ', len(c2))

r1 = list(c1[0].get_residues())
r2 = list(c2[0].get_residues())

print('Residues: ', len(r1), ' ', len(r2))

for i in range(len(r1)):
	resname1 = r1[i].get_resname()
	resname2 = r2[i].get_resname()

	if resname1 != resname2:
		print('Residues not matching')
		break

	atoms1 = list(r1[i].get_atoms())
	atoms2 = list(r2[i].get_atoms())

	atom_names1 = sorted([atom.get_id() for atom in atoms1])
	atom_names2 = sorted([atom.get_id() for atom in atoms2])

	if atom_names1 != atom_names2:
		print('Different atoms!')

d1 = extract_torsion_angles(s1)
d2 = extract_torsion_angles(s2)

print('Torsion angles: ', len(d1), ' ', len(d2))

diff_keys = set(d2.keys()) - set(d1.keys())
print(diff_keys)