Identifiying a conserved residue in multiple PDB structures

Question

I have a few hundred PDB structures of the same protein, and I need to identify a specific conserved residue in all of them.

Originally I wanted to extract the sequences from the PDB files with BioPython, align them and then use that alignment to identify the residue number in each PDB for my conserved residue.

But it turns out that is a bit harder than I expected. Parts of the sequence are often missing in the PDB. I don't see an obvious way to get a proper sequence out of the PDB files. Then there is the issue of mapping the conserved residue number back into the PDB file.

Is there any easier way to do this? It doesn't have to be BioPython, though that would be nice.

Answer 1

My approach to this problem would be to use VMD (Visual Molecular Dynamics), where you can load multiple PDBs, perform structural and/or sequence alignment and analyse residue conservation within one program.

VMD is a powerful molecular visualization program for displaying, animating, and analyzing large biomolecular systems using 3-D graphics and built-in scripting. It is probably an overkill for what you want to do, but I find it rather intuitive to use and powerfull. I also has an amazing documentation and many step-by-step tutorials, starting from this basic introduction.

What you will specifically need for your problem is a MultiSeq 2.0 plugin that you can find under Extensions->Analysis->MultiSeq in new versions of VMD. It is a bioinformatics environment that allows you to load, display, and analyze both sequence and structure data for your proteins within a few mouse clicks. It also has a great set of tutorials to help you start! What you might find interesting for comparing many similar protein structures is an option to create, view and manipulate phylogenetic trees.

Here is step by step tutorial on comparing structures and sequences with MultiSeq.

Here you can find original publications for both VMD and MultiSeq:

VMD:

Humphrey, W., Dalke, A. and Schulten, K., VMD - Visual Molecular Dynamics., J. Molec. Graphics, 1996, 14:33-38.

MultiSeq:

Roberts, E., Eargle, J., Wright, D. and Luthey-Schulten, Z., MultiSeq: Unifying sequence and structure data for evolutionary analysis. BMC Bioinformatics, 2006, 7:382.

Answer 2

Incomplete sequences is a common problem. One way of getting around it is submitting your PDB ID list to an ID mapper. The one at Uniprot works well. Simply copy and paste your PDB ID codes. Make sure you're going from PDB to UniprotKB (see picture below. Hopefully you'll be able to map most if not all of your IDs.

After you have that, download the text file. Then you can extract the sequences from that text file with the following script:

from Bio import SeqIO
filenames = ["YOURINPUTFILE.txt"]
input_format = "swiss" 
output_filename = "YOUROUTPUTFILE.fasta" 

output = open(output_filename, "w")
for filename in filenames:
    for record in SeqIO.parse(filename, input_format):
        sequence = record.seq
        output.write(">%s\n%s\n" % (record.id, sequence))

Then submit those sequences to your favourite MSA tool. I answered a question on how to make an MSA here if you've not made one before.

A minor, and perhaps obvious point on reflection, is that the residue number won't be the same in all your sequences due to indels. You will need to reassess this situation when you get there if you still plan on mapping the conserved residues to 3D structures rather than a sequence.

Answer 3

You can use Ensembl (http://www.ensembl.org/index.html) to retrieve orthologues of the protein, and align them using one of many free, web based programs (such as ClusterOmega)