Aligning PDB structures with Biopython

You can use the Bio.PDB module in Biopython to align PDB files. This is how I did it. The code should be pretty much self-explanatory.

In this example I align the crystal structure of Ubiquitin (PDB code: 1UBQ) to the first structure of a corresponding NMR ensemble (PDB code: 1D3Z, see picture below).

	# The MIT License
	#
	# Copyright (c) 2010-2016 Anders S. Christensen
	#
	# Permission is hereby granted, free of charge, to any person obtaining a copy
	# of this software and associated documentation files (the "Software"), to deal
	# in the Software without restriction, including without limitation the rights
	# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
	# copies of the Software, and to permit persons to whom the Software is
	# furnished to do so, subject to the following conditions:
	#
	# The above copyright notice and this permission notice shall be included in
	# all copies or substantial portions of the Software.
	#
	# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
	# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
	# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
	# # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
	# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
	# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
	# THE SOFTWARE.
	import Bio.PDB
	# Select what residues numbers you wish to align
	# and put them in a list
	start_id = 1
	end_id = 70
	atoms_to_be_aligned = range(start_id, end_id + 1)
	# Start the parser
	pdb_parser = Bio.PDB.PDBParser(QUIET = True)
	# Get the structures
	ref_structure = pdb_parser.get_structure("reference", "1D3Z.pdb")
	sample_structure = pdb_parser.get_structure("samle", "1UBQ.pdb")
	# Use the first model in the pdb-files for alignment
	# Change the number 0 if you want to align to another structure
	ref_model = ref_structure[0]
	sample_model = sample_structure[0]
	# Make a list of the atoms (in the structures) you wish to align.
	# In this case we use CA atoms whose index is in the specified range
	ref_atoms = []
	sample_atoms = []
	# Iterate of all chains in the model in order to find all residues
	for ref_chain in ref_model:
	# Iterate of all residues in each model in order to find proper atoms
	for ref_res in ref_chain:
	# Check if residue number ( .get_id() ) is in the list
	if ref_res.get_id()[1] in atoms_to_be_aligned:
	# Append CA atom to list
	ref_atoms.append(ref_res['CA'])
	# Do the same for the sample structure
	for sample_chain in sample_model:
	for sample_res in sample_chain:
	if sample_res.get_id()[1] in atoms_to_be_aligned:
	sample_atoms.append(sample_res['CA'])
	# Now we initiate the superimposer:
	super_imposer = Bio.PDB.Superimposer()
	super_imposer.set_atoms(ref_atoms, sample_atoms)
	super_imposer.apply(sample_model.get_atoms())
	# Print RMSD:
	print super_imposer.rms
	# Save the aligned version of 1UBQ.pdb
	io = Bio.PDB.PDBIO()
	io.set_structure(sample_structure)
	io.save("1UBQ_aligned.pdb")

view raw align.py hosted with ❤ by GitHub

Numpy vs. cPickles (Python, ofc)

I've been using cPickels for storing data into a Python-friendly format for some time. See my earlier blog post for more on cPickles.
http://combichem.blogspot.dk/2013/02/saving-into-data-into-cpickle-format-in.html

I have also been using Numpy's save function to do the same thing. numpy.save() and numpy.load() is so much simpler, however. I really recommend that people use numpy.save() and numpy.load() over cPickles for most purposes. It is so much more simple.

I always thought a cPickle was much, much faster than Numpy, but I guess I was wrong, according to this stackoverflow I just saw. Below are loading and saving times for a large array. Practically no difference between Numpy and cPickles!

Source: http://stackoverflow.com/questions/16833124/pickle-faster-than-cpickle-with-numeric-data








To save an array, a list or dictionary or whatever called my_array into my_file.npy:

  numpy.save("my_file", my_array)

Note that Numpy appends .npy to the filename automatically.


To load the stored data simply:

  my_array = numpy.load("my_file.npy")

 Really py-fragging-thonicly easy!

You know what really grinds my gears? (In Python)

I can never correctly remember when things are passed as references or copied as local variables inside functions.


Take these two, innocuously looking functions. Because both do the same thing (namely set the contents of a vector, P, to [1, 1]) I call them 1 and a, respectively, since one is not better than the other.


def implementation_1(P):

    P = [1, 1]


def implementation_a(P):

    P[0] = 1
    P[1] = 1


What you would expect is one of the following two options

1. Both functions change P to [1, 1] (permanently).
2. Both functions take a local copy of P and change it to [1, 1], and after the function returns, the local [1, 1] array is forgotten.

A simple test is to do this:


P = [0, 0]
print P

implementation_1(P)
print P

implementation_a(P)
print P

which prints:

[0, 0]
[0, 0]
[1, 1]

So clearly implementation_a() is different from implementation_1(), although they seemingly do the same.

	def implementation_1(P):
	P = [1, 1]
	def implementation_a(P):
	P[0] = 1
	P[1] = 1
	P = [0, 0]
	print P
	implementation_1(P)
	print P
	implementation_a(P)
	print P

view raw list.py hosted with ❤ by GitHub