The BackMAP Python module: how a simpler Ramachandran number can simplify the life of a protein simulator

Ramachandran numbers are structurally meaningful

In addition to resolving positions of secondary structures (Fig. 3), ℛ corresponds well to structural measures such as radius of gyration (R_g), end-to-end distance (R_e), and chirality (χ). These relationships are shown in Fig. 4. Note that chirality comes in many forms, for example, one could be talking about different stereo-isomers, such as L vs D amino acids, or one could be concerned with left-twisting versus right-twisting backbones, that is, handedness (Mannige, 2017). This report will primarily be focused on chirality in context of backbone twist/handedness.

The trends in Fig. 4 show that as one progresses from low to high ℛ, various structural properties also progress smoothly. Additionally, backbones that display similar ℛ also show little variation in structural properties, as evidenced by the small standard deviation bars. It is also important to note that the standard deviations shown in Fig. 4 were calculated by first populating every possible region of (ϕ, ψ)-space. However, in reality, most regions of (ϕ, ψ)-space are unoccupied due to steric/electrostatic constraints, which means that these error bars are likely to be even smaller for natural protein backbones than those depicted here. Finally, the ℛ number is calculated by taking “sweeps” of the (ϕ, ψ)-space in lines that are parallel to the negatively-sloping diagonal. Interestingly, such “sweeps” encounter only one major (dense) region within (ϕ, ψ)-space (e.g., ℛ’s in the general vicinity of 0.34 represent structures that resemble α-helices). This means that ℛ can also be used to assess the types of secondary structure present in a protein conformation.

Ramachandran codes are stackable

An important aspect of the Ramachandran number (ℛ) lies in its compactness compared to the traditional Ramachandran pair (ϕ, ψ). The value of the conversion from (ϕ, ψ)-space to ℛ-space is that the structure of a protein can be described in various one-dimensional arrays (per-structure “Ramachandran codes” or “ℛ-codes” or multi-angle maps); see, e.g., Fig. 5.

In addition to assuming a small form factor, ℛ-codes may then be stacked side-by-side for visual and computational analysis. There lies its true power.

For example, the one-ℛ-to-one-residue mapping means that the entire residue-by-residue structure of a protein can be shown using a string of ℛ_is (which would show regions of secondary structure and disorder, for starters). Additionally, an entire protein’s backbone makeup can be shown as a histogram in ℛ-space (which may reveal a protein’s topology). The power of this format lies not only in the capacity to distill complex structure into compact spaces, but in its capacity to display many complex structures in this format, side-by-side (stacking).

Peptoid nanosheets (Mannige et al., 2015) will be used here as an example of how multiple structures, in the form of ℛ-codes, may be stacked to provide immediately useful pictograms. Peptoids are stereo-isomers of peptides, where the sidechain is attached to the backbone nitrogen rather than the carbon atom. Since both peptoids and peptides share identical backbone connectivity, the analysis described below could be applied to both peptides and peptoids.

Peptoid nanosheets are a recently discovered peptide-mimic that, in one molecular dynamics simulation (Mannige et al., 2015), were shown to display a novel secondary structure. In the reported model (Mannige et al., 2015), each peptoid within the nanosheet displays backbone conformations that alternate in chirality, causing the backbone to look like a meandering snake that nonetheless maintains an overall linear direction. This secondary structure was discovered by first setting up a nanosheet where all peptoid backbones were restrained to be fully extended, after which the restraints were energetically softened and completely released. Figure 6A shows snapshots of each of these states. As evident in Figs. 6B and 6C, the two types of ℛ-code stacks display salient information at first glance: (1) Fig. 6B shows that the extended backbone first undergoes some rearrangement with softer restraints, and then becomes more binary in arrangement as we look down the backbone; and (2) Fig. 6C shows that lifting restraints on the backbone causes a dramatic change in backbone topology, namely a birth of a bimodal distribution evident in the two parallel horizontal bands.

By utilizing ℛ, maps such as those in Fig. 6 provide information about every ϕ and ψ within the backbone. As such, these maps are dubbed MAPs. A Python package called BackMAP created Figs. 6A and 6B, and is provided as a GitHub repository (https://github.com/ranjanmannige/BackMAP). BackMAP takes in a PDB structure file containing a single structure, or multiple structures separated by the code “MODEL.”

Case study: picking out subtle differences from high volume of data

This section expands on the notion that ℛ-numbers—due to their compactness/stackability—can be used to pick out backbone structural trends that would be hard to decipher using any other metric. For example, it is well known that prolines (pro) display unusual backbone behavior: in particular, proline backbones occupy structures that are close to but distinct from α-helical regions. Due to the two-dimensionality of Ramachandran plots (Fig. 7A), such distinctions are hard to visually pick out from Ramachandran plots. However, stacking per-amino-acid ℛ-codes side-by-side make such differences patent (Fig. 7B; see arrow).

It is also known that amino acids preceeding prolines display unusual shifts in backbone twist/chirality. For example, Fig. 8C shows that amino acids appearing before prolines behave differently than they would otherwise (see the upward-facing arrow). Additionally, amino acids following glycines also appear to have their structures modified (Fig. 8D; upward arrow). Note that these results are not new, and it has already been confirmed that, for example, nearest neighbors affect the conformational behavior of an amino acid as witnessed within Ramachandran plots (Ting et al., 2010), and proline changes the backbone conformation of the preceeding residue (Gunasekaran et al., 1998; Ho & Brasseur, 2005). However, Figs. 7 and 8 indicate that such information can be more concisely shown/identified when structures are stacked side-by-side in the form of ℛ-codes. Such subtle changes are often witnessed when protein backbones transition from one state to another.

Installation

BackMAP may either be installed locally by downloading the GitHub repository, or installed directly by running the following line in the command prompt (assuming that pip exists): > pip install backmap.

Usage

The module can either be imported and used within existing scripts, or used as a standalone package using the command “python -m backmap.” First the in-script usage will be discussed.

# Import module                                        1
import backmap                                         2
# Convert (phi, psi) to R                                    3
print backmap.R(phi =0, phi =0) # Expected output : 0.5                  4
print backmap.R( −180, −180) # Expected output : 0.0                     5
print backmap.R(  180,  180) # Expected output : 1.0 (equivalent in meaning to 0) 6

In-script usage II: basic usage for creating multi-angle pictures

The following code shows how MAPs of protein backbones can be generated:

• 1. Select and read a protein PDB file

  Each trajectory frame must be a set of legitimate protein databank (PDB) “ATOM” records separated by “MODEL” keywords (distinct models show up as distinct frames on the x-axis or abscissa).

  
  import backmap                                             1
  import matplotlib.pyplot as plt  # Pyplot is needed for drawing graphs       2
  import numpy as np      # Numpy is needed for data transformations        3
  pdbfn = './tests/pdbs/1xqq.pdb'  # Set pdb name                         4
  data = backmap.read_pdb(pdbfn) # READ PDB in the form of a matrix with columns 5
  

Here, “data” is a 2d array with four columns [“model,” “chain,” “resid,” “R”]. The first row of “data” is the header, with values that follow.

• 2. Select color scheme (color map)

  In addition to custom colormaps listed in the next section, one can also use traditional colormaps available at matplotlib.org (e.g., “Reds” or “Reds_r”).

  
  # setting the name of the colormap         6
  cmap = "SecondaryStructure"                    7
  

• 3. Draw per-chain MAPs

  
  # Grouping by chain                                8
  grouped_data = backmap.group_data_by (data, group_by='chain',      9
                   columns_to_return=['model','resid','R']) 10
  for chain in grouped_data.keys(): # Going through each chain          11
   # Getting the X,Y,Z values for each entry                     12
   models, residues, Rs = grouped_data [chain]                     13
   # Finally, creating (but not showing) the graph                 14
   backmap.draw_xyz(X = models ,   Y = residues ,  Z = Rs       15
         , xlabel = 'Frame #',  ylabel = "Residue #",zlabel = 'ℛ'    16
          ,cmap = cmap  , title = "Chain: ' "+chain+" ' "        17
          ,vmin=0,vmax=1)                                18
   # Now, we display the graph:                              19
   plt.show() # ... one can also use plt.savefig() to save to file      20
  

The code above results in Fig. 9A.

Additionally, by changing how one assigns values to “X” and “Y,” one can easily construct and draw other types of graphs such as time-resolved histograms, per-residue fluctuations when compared to the first (D₁) and previous structure (D₋₁) within the trajectory, etc.

for chain in grouped_data.keys():                                    11
  models, residues, Rs = grouped_data[chain]                             12
                                                            13
  'Begin custom code'                                          14
  X = []; Y=[]; Z=[]; # Will set X=model, Y=R, Z=P(R)                       15
  # Bundling the three lists into one 2d array                            16
  new_data =  np.array(list(zip(models , residues , Rs)))                     17
  # Getting all R values, model by model                               18
  for m in sorted(set(new_data[:,0])): # column 0 is the model column             19
   # Getting all Rs for that model #                                20
   current_rs = new_data[np.where(new_data[:,0]==m)][:,2] # column 2 contains R     21
   # Getting the histogram                                      22
   a,b = np.histogram(current_rs,bins=np.arange(0,1.01,0.01))                   23
   max_count = float(np.max(a))                                   24
   for i in range(len(a)):                                      25
    X.append(m); Y.append((b[i]+b[i+1])/2.0); Z.append(a[i]/float(np.sum(a)));     26
  'End custom code'                                            27
                                                            28
  # Finally, creating (but not showing) the graph                         29
  backmap.draw_xyz(X = X  ,  Y = Y ,     Z = Z   30
    ,xlabel ='Frame #', ylabel ="ℛ",zlabel ="P'(ℛ)"       31
    ,cmap = 'Greys', ylim=[0,1])                                  32
  plt.yticks(np.arange(0,1.00001,0.2))                                33
  # Now, we display the graph:                                     34
  plt.show() # ... one can also use plt.savefig() to save to file               35

In-script usage IV: available color schemes (CMAPs)

Aside from the general color maps (cmaps) that exist in matplotlib (e.g., “Greys,” “Reds,” or, god forbid, “jet”), BackMAP provides two new colormaps: “Chirality” (key: +ve twists: red; −ve twists: blue), and “SecondaryStructure” (key: potential α-helices: red; β-sheets: blue; ppII-helices: cyan). Right twisting backbones are shown in red; left twisting backbones are shown in blue. Figure 10 shows how a single protein ensemble may be described using these schematics. As illustrated in Fig. 10B, cmaps available within the standard matplotlib package do not distinguish between major secondary structures well, while those provided by BackMAP do. In case it is known that the protein backbone accesses non-traditional regions of the Ramachandran plot, a four-color schematic will be needed (see below for more discussions).

Stand alone usage

BackMAP can be used as a stand alone package by running “> python -m backmap -pdb <pdb_dir_or_file>.” The sectons below describe the expected outputs and how they may be interpreted.

Stand alone Example I: a stable protein

Figures 11B–11F below were created by running “> python -m backmap ./tests/pdbs/1xqq.pdb” (Fig. 11A was created using VMD). These graphs indicate that conformational ensemble 1xqq describes a conformationally stable protein, since 1) each conformation shows little change in the ℛ histogram over time (Fig. 11B), and 2) each residue fluctuates little in color (structure) over “time” (Figs. 11C–11F; see Methods). Here and below, it is assumed that discrete models represent distinct states of the protein over “time”.

In particular, each column in Fig. 11B describes the histogram in Ramachandran number (ℛ) space for a single model/timeframe. These histograms show the presence of both α-helices (at ℛ ≈ 0.34) and β-sheets (at ℛ ≈ 0.52). Additionally, Figs. 11C and 11D describe per-residue conformational plots (colored by two different metrics or CMAPs), which show that most of the protein backbone remains relatively stable over time (e.g., few fluctuations in state or “color” are evident over frame #). Finally, Fig. 11E describes the extent to which a single residue’s state has deviated from the first frame, and Fig. 11F describes the extent to which a single residue’s state has deviated from its state in the previous frame. All these graphs show that this protein is relatively conformationally stable.

Stand alone Example II: an intrinsically disrodered protein

Figure 12 is identical to Fig. 11, except that the panels pertain to an intrinsically disordered protein 2fft whose structural ensemble describes dramatically distinct conformations.

As compared to the conformationally stable protein above, protein 2fft is much more flexible. Figure 11B shows that the states accessed per model are diverse and dramatically fluctuate over the entire range of ℛ (this is especially true when compared to a stable protein, see Fig. 11B).

The diverse states occupied by each residue (Figs. 11C and 11D) confirm the conformational variation displayed by most of the backbone (Figs. 11E and 11F similarly show that most of the residues fluctuate dramatically).

Yet, interestingly, Figs. 11C–11F also show an unusually stable region—residues 15–25—which consistently display the same conformational (α-helical) state at ℛ ≈ 0.34 (interpreted as the color red in Fig. 11C). This trend would be hard to recognize by simply looking at the structural ensemble (Fig. 11A).

A signed Ramachandran number for “misbehaving” backbones

The Ramachandran number increases in value from the bottom left of the Ramachandran plot to the top right in sweeps that are parallel to the negative sloping diagonal. As discussed in Mannige, Kundu & Whitelam (2016), this method of mapping a two-dimensional space into one number is still structurally meaningful and descriptive because (1) most structural features of the protein backbone—for example, radius of gyration (Mannige, Kundu & Whitelam, 2016), end-to-end distance (Mannige, Kundu & Whitelam, 2016), and chirality (Mannige, 2017)—vary little along lines parallel to the negatively-sloping diagonal (this is indicated by relatively small standard deviations in structural metrics for similar ℛ_s; Fig. 4), and (2) most protein backbones display chiral centers and therefore predominantly appear on the top left region of the Ramachandran plot (above the dashed diagonal in Fig. 13A).

However, not all backbones localize in only one half of the Ramachandran plot. Particularly, among biologically relevant amino acids, glycine occupies both regions of the Ramachandran plot (Fig. 13B; of note, the α_L helix region becomes relatively prominent). On the other hand, prolines are known to form polyproline-II helices (ppII in Fig. 13C), which falls on almost the same “sweep” as glycine rich peptides (red dot-dashed line). In situations where both prolines and glycines are abundant, the Ramachandran number (ℛ) would fail to distinguish α_L from ppII (Fig. 13D; regions outlined by rectangles).

To accomodate the situation where achiral backbones are expected (e.g., if peptoids or polygycines are being studied), an additional Ramachandran number—the signed Ramachandran number ℛ_s—is introduced here. ℛ_s is identical to the original number in magnitude, but which changes sign from + to − as you approach ℛ numbers that are to the right (or below) the positively sloped diagonal i.e., (3) $${ℛ}_{\text{S}}=\{\begin{array}{cc}ℛ& ,\text{if\hspace{0.17em}}\psi \ge \varphi \\ ℛ\times -1& ,\text{if\hspace{0.17em}}\psi <\varphi \end{array}$$ Signed Ramachandran numbers are calculated by adding “-signed” within the command line implementation or by adding “signed=True” when making backmap.R() calls.

As an example of the utility of ℛ_s, Fig. 13F shows that ℛ_s easily distinguishes α_D from ppII.

Note that the signed ℛ_s, while useful, would be important in very limited scenarios, as more than 96% of the amino acids in the PDB occupy the upper-left region of the Ramachandran plot (with 3% of “rule breakers” contributed mostly by glycines).

Simplifying the Ramachandran number (ℛ)

This section will derive the simplified Ramachandran number presented in this paper from the more complicated looking Ramachandran number introduced previously (Mannige, Kundu & Whitelam, 2016).

Assuming the bounds ϕ ∈ [ϕ_min, ϕ_max) and ϕ ∈ [ψ_min, ψ_max), the previously described Ramachandran number takes the form (5) $$ℛ(\varphi \mathrm{,}\psi )\equiv \frac{R_{\mathbb{Z}}(\varphi \mathrm{,}\psi )-R_{\mathbb{Z}}({\varphi }_{\text{min}}\mathrm{,}{\varphi }_{\text{min}})}{R_{\mathbb{Z}}({\varphi }_{\text{max}}\mathrm{,}{\varphi }_{\text{max}})-R_{\mathbb{Z}}({\varphi }_{\text{min}}\mathrm{,}{\varphi }_{\text{min}})}\mathrm{,}$$ where, $ℛ(\varphi \mathrm{,}\psi )$ is the Ramachanran number with range (0, 1), and $R_{\mathbb{Z}}(\varphi \mathrm{,}\psi )$ is the unnormalized integer-spaced Ramachandran number whose closed form is (6) $$R_{\mathbb{Z}}(\varphi \mathrm{,}\psi )=\left[(\varphi -\psi +\lambda )\sigma /\sqrt{2}\right]+\left[\sqrt{2}\lambda \sigma \right]\left[(\varphi +\psi +\lambda )\sigma /\sqrt{2}\right]\mathrm{.}$$

Here, [x] rounds x to the closest integer value, σ is a scaling factor, discussed below, and λ is the range of an angle in degrees (i.e., λ = ϕ_max − ϕ_min). Effectively, this equation does the following: (1) It divides up the Ramachandran plot into (360° σ^1/°)² squares, where σ is a user-selected scaling factor that is measured in reciprocal degrees (see Fig. 8B in Mannige, Kundu & Whitelam (2016)); (2) It then assigns integer values to each square by setting the lowest integer value to the bottom left of the Ramachandran plot (ϕ = −180°, ψ = −180°), and proceeding from the bottom left to the top right of the Ramachandran plot by iteratively slicing down −1/2 sloped lines and assigning increasing integer values to each square that one encounters; (3) Finally, the equation assigns any (ϕ, ψ) pair within ϕ, ψ ∈ [−ϕ_min, ϕ_max) to the integer value (R_ℤ) assigned to the divvied-up square that they it exists in. Further discussions on this process are available in the previous publication (Mannige, Kundu & Whitelam, 2016).

Combining the two equations (Eqs. (5) and (6)) results in the following, rather imposing, equation for the Ramachandran number: (7) $$ℛ(\varphi \mathrm{,}\psi )=\frac{\left(\begin{array}{cc}\left[(\varphi -\psi +\lambda )\sigma /\sqrt{2}\right]& +\left[\sqrt{2}\lambda \sigma \right]\left[(\varphi +\psi +\lambda )\sigma /\sqrt{2}\right]\\ -\left[({\varphi }_{\text{min}}-{\psi }_{\text{min}}+\lambda )\sigma /\sqrt{2}\right]& -\left[\sqrt{2}\lambda \sigma \right]\left[({\varphi }_{\text{min}}+{\psi }_{\text{min}}+\lambda )\sigma /\sqrt{2}\right]\end{array}\right)}{\left(\begin{array}{cc}\left[({\varphi }_{\text{max}}-{\psi }_{\text{max}}+\lambda )\sigma /\sqrt{2}\right]& +\left[\sqrt{2}\lambda \sigma \right]\left[({\varphi }_{\text{max}}+{\psi }_{\text{max}}+\lambda )\sigma /\sqrt{2}\right]\\ -\left[({\varphi }_{\text{min}}-{\psi }_{\text{min}}+\lambda )\sigma /\sqrt{2}\right]& -\left[\sqrt{2}\lambda \sigma \right]\left[({\varphi }_{\text{min}}+{\psi }_{\text{min}}+\lambda )\sigma /\sqrt{2}\right]\end{array}\right)}$$

However, useful Eq. (7) is, the complexity of the equation may be a deterrent toward utilizing it. This paper reports a simpler equation that is derived by taking the limit of Eq. (7) as σ tends toward ∞. In particular, when σ→∞, Eq. (7) becomes (8) $$ℛ(\varphi \mathrm{,}\psi )=\underset{\sigma \to \infty }{{\displaystyle \lim }}\overline{ℛ}(\varphi \mathrm{,}\psi )=\frac{\varphi +\psi -({\psi }_{\text{min}}+{\psi }_{\text{min}})}{({\varphi }_{\text{max}}+{\psi }_{\text{max}})-({\varphi }_{\text{min}}\text{+}{\psi }_{\text{min}})}\mathrm{.}$$

Assuming that ϕ, ψ ∈ [−180°, 180°) or [−π, π), (9) $$ℛ(\varphi \mathrm{,}\psi )=\frac{\varphi +\psi +2\pi }{4\pi }\mathrm{.}$$

Conformation of this limit is shown numerically in Fig. A1. Since larger σs indicate higher accuracy, ${{\displaystyle \underset{\sigma \to \infty }{\lim }}}_{}ℛ(\varphi \mathrm{,}\psi )$ represents an exact representation of the Ramachandran number. Using this closed form, this report shows how both static structural features and complex structural transitions may be identified with the help of Ramachandran number-derived plots.

Assuming, a different range (say, ϕ, ψ ∈ [0, 2π)), the Ramachandran number in that frame of reference will be (10) $$ℛ{(\varphi \mathrm{,}\psi )}_{\varphi \mathrm{,}\psi \in [0,2\pi )}=\frac{\varphi +\psi }{4\pi }\mathrm{.}$$ However, in changing the ranges, the meaning of the Ramachandran number will change. This manuscript assumes that all angles (ϕ, ψ, ω) range between −π (−180°) and π (180°).