PeptideBuilder: A simple Python library to generate model peptides

We present a simple Python library to construct models of polypeptides from scratch. The intended use case is the generation of peptide models with pre-specified backbone angles. For example, using our library, one can generate a model of a set of amino acids in a specific conformation using just a few lines of python code. We do not provide any tools for energy minimization or rotamer packing, since powerful tools are available for these purposes. Instead, we provide a simple Python interface that enables one to add residues to a peptide chain in any desired conformation. Bond angles and bond lengths can be manipulated if so desired, and reasonable values are used by default.


INTRODUCTION
Researchers working in structural biology and related fields frequently have to create, manipulate, or analyze protein crystal structures.To aid this work, many different software tools have been developed.Examples include visualization (Schrödinger, 2013), mutagenesis (Schrödinger, 2013;Leaver-Fay et al., 2011), high-throughput computational analysis (Hamelryck & Manderick, 2003;Grant et al., 2006), ab-initio protein folding and protein design (Leaver-Fay et al., 2011), and homology modeling and threading (Eswar et al., 2006;Zhang, 2008).In comparison, a relatively simple task, the ab-initio creation of a protein structure in a desired conformation, has received little attention.It is possible to perform this task in PyRosetta (Chaudhury, Lyskov & Gray, 2010;Gray et al., 2013), but that approach incurs the overhead of the entire Rosetta protein modeling package (Leaver-Fay et al., 2011).One can also construct peptides manually in some graphical molecular modeling packages, such as Swiss-PdbViewer (Guex & Peitsch, 1997).Finally, the Rose lab has developed Ribosome (Srinivasan, 2013), a small program with the express purpose of creating model peptides.However, Ribosome is implemented in Fortran, an outdated programming language that integrates poorly with modern bioinformatics pipelines.
For a recent analysis by our group, we wanted to systematically enumerate GLY-X-GLY tripeptides in all allowed conformations (Tien et al., 2012).After review of the available software packages, we determined that there was a need for a lightweight library, implemented in a modern programming language, that would allow us to construct arbitrary peptides in any desired conformation.We decided to write this library in the language Python (Python Sofware Foundation, 2013), as this language is widely used in scientific computing.Specifically, many tools suitable for computational biology and bioinformatics are available (Cock et al., 2009), including tools to read, manipulate, and write PDB (Protein Data Bank) files (Hamelryck & Manderick, 2003).This effort resulted in the Python library PeptideBuilder, which we describe here.The library consists of two Python files comprising a total of approximately 2000 lines of code.Both files are provided as Supplemental Information 1.The entire PeptideBuilder package is also available online at https://github.com/mtien/PeptideBuilder.

CONCEPTUAL OVERVIEW
The key function our library provides is to add a residue at the C terminus of an existing polypeptide model, using arbitrary backbone angles.Our library also allows a user to generate an individual amino acid residue and place it into an otherwise empty model.In combination, these two functions enable the construction of arbitrary polypeptide chains.The generated models are stored as structure objects using the PDB module of Biopython (Bio.PDB, Hamelryck & Manderick 2003).The seemless integration with Biopython's PDB module means that we can leverage a wide range of existing functionality, such as writing structures to PDB files or measuring distances between atoms.
Adding a residue to an existing polypeptide chain involves two separate steps.First, we have to establish the desired geometric arrangement of all atoms in the residue to be added.This means we have to determine all bond lengths and angles.In practice, we will usually want to specify the dihedral backbone angles φ and ψ, and possibly the rotamers, whereas all other bond lengths and angles should be set to reasonable defaults for the amino acid under consideration.Once we have determined the desired geometry, we have to calculate the actual position of all atoms in 3-space and then add the atoms to the structure object.The exact calculations required to convert bond lengths and angles into 3D atom coordinates are given in the supporting text of Tien et al. (2012).Our library places all heavy atoms for each residue, but it does not place hydrogens.
We obtained default values for bond lengths and angles by measuring these quantities in a large collection of published crystal structures and recording the average for each quantity, as described (Tien et al., 2012).We set the default for the backbone dihedral angles to the extended conformation (φ = −120 • , ψ = 140 • , ω = 180 • ).We based the default rotamer angles of each individual amino acid on the rotamer library of Shapovalov & Dunbrack (2011).For each amino acid, the rotamer library provided the frequency of each combination of rotamer angles given the backbone conformation.We analyzed this library at the extended backbone conformation (φ = −120 • , ψ = 140 • ) and used the most likely rotamer conformation at that backbone conformation as default.For amino acids with multiple χ 1 or χ 2 angles (such as Ile or Leu) we used the most likely and the second most likely angles from the rotamer library.

EVALUATION
We did extensive testing on our library to verify that we were placing atoms at the correct locations given the bond lengths and angles we specified.First, we collected tri-peptides from published PDB structures, extracted all bond lengths and angles, reconstructed the tri-peptides using our library, and verified that the original tri-peptide and the reconstructed one aligned with an RMSD of zero.Next, we wanted to know how our library would fare in reconstructing longer peptides, in particular when using the default parameter values we used for bond lengths and angles.For this analysis, we focused on the peptide backbone, since the evaluation of tripeptides had shown that our library was capable of placing side-chains correctly if it was given the correct bond lengths and angles.
We selected ten proteins with solved crystal structure.The proteins were chosen to represent a diverse group of common folds.For each protein, we then attempted to reconstruct the backbone of either the first 50 residues in the structure, the first 150 residues in the structure, or the entire structure.In all cases, we extracted backbone bond lengths and angles at each residue, and then reconstructed the protein using four different methods.When placing each residue, we either (i) adjusted only the extracted φ and ψ dihedral angles, (ii) adjusted φ, ψ, and ω dihedral angles, (iii) adjusted all dihedral and planar bond angles, or (iv) adjusted all bond lengths and angles exactly to the values measured in the structure we were reconstructing.In each case, any remaining parameters were left at their default values.
As expected, when we set all bond lengths and angles to exactly the values observed in the reference crystal structure, we could reconstruct the entire backbone with an RMSD close to zero.We did see an accumulation of rounding errors in longer proteins, but these rounding errors amounted to an RMSD of less than 0.01 Å even for a protein of over 600 residues.Hence they are negligible in practice.By contrast, reconstructions relying on just backbone dihedral angles performed poorly.We found that we had to adjust all backbone bond angles, inlcuding planar angles, to obtain accurate reconstructions.Bond lengths, on the other hand, could be left at their default values.Table 1 summarizes our findings for all 10 structures, and Fig. 1 shows the results of the four different methods of reconstruction for one example structure.The python script to generate these reconstructions is provided as part of Supplemental Information 1.
Our results show that the PeptideBuilder software correctly places all atoms at the desired locations.However, they also demonstrate that one needs to be careful when constructing longer peptides.It is not possible to construct an entire protein structure just from backbone dihedral angles and expect the structure to look approximately correct.In particular in tight turns and unstructured loops, small deviations in backbone bond angles can have a major impact on where in 3D space downstream secondary structure elements are located.Hence these angles cannot be neglected when reconstructing backbones.

USAGE
The PeptideBuilder software consists of two libraries, Geometry and PeptideBuilder.The Geometry library contains functions that can set up the proper three-dimensional geometry of all 20 amino acids.The PeptideBuilder library contains functions that use this geometry information to construct actual peptides.
PeptideBuilder has one dependency beyond a default python installation, the Biopython package (Cock et al. 2009, http://biopython.org/),which provides the module Bio.
Table 1 Root mean square deviation (RMSD) between reference crystal structures and reconstructions of these structures using varying amounts of modeling detail.Reconstructions using just dihedral backbone angles tend to deviate substantially from the reference structures, whereas reconstructions using all bond angles tend to perform well, even if bond lengths are kept at default values.

Notes.
a PDB ID of the reference structure.In all cases, chain A of the structure was used.b Length of reference structure, in amino acids.c Only dihedral backbone angles φ and ψ have been adjusted to match the reference structure.d Only dihedral backbone angles φ, ψ, and ω have been adjusted to match the reference structure.e RMSD (in Å) over the first 50 residues.f RMSD (in Å) over the first 150 residues.g RMSD (in Å) over the entire structure.

Basic construction of peptides
Let us consider a simple program that places 5 glycines into an extended conformation.First, we need to generate the glycine geometry, which we do using the function Geometry.geometry().This function takes as argument a single character indicating the desired amino acid.Using the resulting geometry, we can then intialize a structure object with this amino acid, using the function PeptideBuilder.initializeres().
If this is the conformation we want to generate, we can simply reuse the previously generated geometry object and write: Creates a new structure containing a single amino acid.make structure Builds an entire peptide with specified amino acid sequence and backbone angles.make extended structure Builds an entire peptide in the extended conformation.make structure from geos Builds an entire peptide from a list of geometry objects.
7 out = Bio .PDB .PDBIO () 8 out .set_structure ( structure ) 9 out .save ( " example .pdb " ) If we want to generate a peptide that is not in the extended conformation, we have to adjust the backbone dihedral angles accordingly.For example, we could place the five glycines into an alpha helix by setting φ = −60 • and ψ = −40 • .We do this by manipulating the phi and psi im1 members of the geometry object.(We are not actually specifying the ψ angle of the residue to be added, but the corresponding angle of the previous residue, ψ i−1 .Hence the member name psi im1.See the Detailed adjustment of residue geometry Section for details.)The code example looks as follows: 1 geo = Geometry .geometry ( " G " ) Several convenience functions exist that simplify common tasks.For example, if we simply want to add a residue at specific backbone angles, we can use an overloaded version of the function add residue() that takes as arguments the structure to which the residue should be added, the amino acid in single-character code, and the φ and ψ i−1 angles: 1 # add an arginine , setting phi = -60 and psi_im1 = -40 2 structure = PeptideBuilder .add_residue ( structure , " R " , -60 , -40) If we want to place an arbitrary sequence of amino acids into an extended structure, we can use the function make extended structure(), which takes as its sole argument a string holding the desired amino-acid sequence: 1 # construct a peptide corresponding to the 2 # sequence " MGGLTR " in extended conformation 3 structure = PeptideBuilder .m ak e _e xt e nd e d_ st r uc t ur e ( " MGGLTR " ) Table 2 summarizes all functions provided by PeptideBuilder.All these functions are documented in the source code using standard Python self-documentation methods.

Detailed adjustment of residue geometry
Geometry objects contain all the bond lengths, bond angles, and dihedral angles necessary to specify a given amino acid.These parameters are stored as member variables, and can be In this function call, the angles are listed in order of standard biochemical convention, χ 1 , χ 2 , χ 3 , and so on, for however many χ angles the amino-acid side chain has.

CONCLUSION
We have developed a Python library to construct model peptides.Our design goals were to make the library simple, lightweight, and easy-to-use.Using our library, one can construct model peptides in only a few lines of Python code, as long as default bond lengths and angles are acceptable.At the same time, all bond-length and bond-angle parameters are user-accessible and can be modified if so desired.We have verified that our library places atoms correctly.As part of this verification effort, we have found that with increasing peptide length it becomes increasingly important to adjust bond angles appropriately to reconstruct biophysically accurate protein structures.

Figure 1
Figure 1 Reconstruction of protein backbone using varying degrees of modeling accuracy.The gray backbone corresponds to chain A of crystal structure 1gfl (green fluorescent protein), and the rainbowcolored backbone corresponds to the reconstructed version thereof.(A) Only φ and ψ dihedral angles are adjusted to match those in the reference structure.(B) All dihedral backbone angles (φ, ψ, and ω) are adjusted to match those in the reference structure.(C) All backbone bond angles are adjusted to match those in the reference structure.(D) All backbone bond lengths and angles are adjusted to match those in the reference structure.RMSD values are given in Table 1.Part (D) shows perfect overlap between the reference and the reconstructed backbone.