3D based on 2D: Calculating helix angles and stacking patterns using forgi 2.0, an RNA Python library centered on secondary structure elements.

Implementation

The forgi library¹⁷ is strongly object-oriented, but takes advantage of module-level Python functions where appropriate. The core object representing the secondary structure is the BulgeGraph object, which holds a Sequence instance for the primary sequence. To include secondary structure based 3D coordinates, the BulgeGraph’s subclass, CoarseGrainRNA is available. For all-atom 3D analysis, the forgi library has a built-in integration with Biopython.

We will briefly describe the three main data structures that hold the sequence, secondary structure and the tertiary structure representation of an RNA.

Primary structure: The Sequence class. Each BulgeGraph holds a Sequence object. Since forgi supports loading of data from PDB files (see below), this Sequence object has to account for many special cases arising from experimental considerations, which will be detailed in the following paragraphs. There are two numbering schemes commonly used for sequences: 1-based indexing and indexing based on an external reference. In particular, many sequences in structural experiments use “1” to dedicate the first residue of a biological macro-molecule, whereas the first residue actually used in the experiment can be upstream of the functional RNA (leading to negative indices) or after the start of the biological unit (leading to an index above 1). The latter is especially common if fragments of larger RNA molecules like the ribosome are studied. Finally, experimenters might decide to insert nucleotides in the middle of a molecule. In order not to affect the numbering of subsequent residues, these inserted residues get the same number as the previous one, followed by a letter (called insertion code).

To handle both kinds of indexing, the Sequence class distinguishes indices by type. Integer indices always refer to 1-based indexing, while tuples compatible to the indices used in Biopython’s PDB module are interpreted as the second kind of indices.

For many applications, it is necessary to restrict the RNA alphabet to 4 letters (i.e. 4 residue types), “G”, “C”, “A” and “U”. However, in the cell many RNA molecules are post-transcriptionally modified at certain positions. Many modifications, including the methylation of OH or NH₂ groups and A to I editing, have been implicated with a variety of biological functions¹⁸. During parsing of PDB files, we automatically convert such modified residues to the unmodified parent, but in addition store the modification as an annotation in the Sequence object. The corresponding unmodified parents for 3-letter codes of common modified residues were obtained from PDBeChem¹⁹ (http://www.ebi.ac.uk/pdbe-srv/pdbechem/) and the forgi library has the ability to query this database on the fly if it encounters a new 3-letter code.

Finally, many experimental 3D structures do not contain coordinates for all residues present in the experiment. The forgi Sequence class stores two version of the sequence, with and without missing residues.

The secondary structure of an RNA is internally represented as a graph, the Bulge Graph, where secondary structure elements (stems, single-stranded regions, interior loops and hairpin loops) form the nodes. Whenever these elements are adjacent along the backbone, they are connected by an edge in this graph. During the Bulge Graph creation, each node gets a unique name such as “s0” for the first stem or “h0” for the first hairpin. The concept of the Bulge Graph, illustrated in Figure 1, has been described previously in more detail²⁰ and is related to the independently developed RAG (RNA as Graph) approach²¹.

Figure 1. Illustration of the Bulge Graph representation underlying the forgi library.

Secondary structure elements (stems, bulges and single-stranded regions) are nodes connected by edges. The sequence is shown around the Bulge Graph.

forgi supports element-based transformations of the Bulge Graph, such as condensing the secondary structure to an representation similar to RNAshapes²², which we use to classify pseudo knots, for example (see below).

The BulgeGraph object allows for easy identification, selection and classification of structural domains, such as multi loops, helices consisting of multiple stems and bulges (termed “rods” in forgi), and pseudo knots.

The CoarseGrainRNA class holds 3D structures. 3D structures are loaded from PDB or MMCIF files using Biopython’s PDB module¹³. In order to assign a secondary structure to the RNA, forgi can call either MC-Annotate⁴ (Linux only, no MMCIF-support, binary available at http://major.iric.ca/MajorLabEn/MC-Tools.html) or DSSR⁷ (binary available after free registration at http://forum.x3dna.org). As an alternative, forgi has a built-in heuristic for the detection of canonical base pairs and GU wobble pairs. This heuristic is based on distances along the hydrogen bonds and the coplanarity of the bases, and is not intended to compete with the power of more specialized tools since it will fail in some edge cases, e.g. involving modified residues or residues with missing atoms. This heuristic is a useful fallback, if the above mentioned programs are not available.

During loading of the RNA, a helix axis is assigned to stems as described previously²⁰. For each stem, we store start and end coordinates of the helix axis as well as two twist vectors that point towards the minor groove at the beginning and the end of the stem. Similarly, start and end coordinates for bulges, loops and single-stranded regions are stored and can be accessed using the element’s name.

forgi 2.0 now fully supports co- and multi-fold structures and can load multiple chains that are connected by base pairs into a single CoarseGrainRNA object, while chains not connected by any base pair are loaded into separate objects.

Using Biopython’s KD-Tree implementation (Bio.PDB.NeighborSearch), a list of residues within 6 angstrom from a non-RNA C or N atom is obtained. These residues are considered protein-/ligand interacting. Knowledge of interacting residues is particularly useful to avoid biases in statistics about structural features of bare RNA.

A cleaned version of the PDB as Biopython chains is stored in the CoarseGrainRNA’s chains attribute. This cleaned version has the names of modified residues converted to their canonical parent and non-RNA molecules removed.

Operation

forgi¹⁷ is compatible with Python 2.7, 3.5 and 3.6, should run on all operating systems where its dependencies are available and has been tested on Linux, Mac and Windows. It makes heavy use of the NumPy²³ library to speed-up array-based calculations and also depends on SciPy²⁴, NetworkX²⁵, Biopython^12,13, pandas^26,27 and appdirs, all available via the Python Package Index (PyPi) or Anaconda.

Helpful utility scripts. forgi comes with the two very useful scripts: rnaConvert.py can be used to convert between many common file formats for RNA structures, including the Vienna format and fasta variants, the bpseq format, the connectivity table (ct) format, MMCIF format and PDB files. visualize_rna.py can be used to display a coarse-grained representation of an RNA’s secondary structure in PyMol alongside the all atom structure from PDB files, producing visualizations like those in Figure 3 and Figure 5.

Analysis of stacking geometries. To illustrate how the forgi library makes secondary structure elements usable for 3D structure analysis, we used it to analyze the stacking of adjacent helices in multi loops and pseudo knots in a representative set of RNA 3D structures²⁸ (version 3.36, available at http://rna.bgsu.edu/rna3dhub/nrlist/).

While loading these 3D structures into forgi, DSSR⁷ (Version v1.7.1-2017nov01) was called to obtain the secondary structure and nucleotide level reference stacking annotations. We count a pair of connected helices as stacking, if at least one nucleotide of the first helix’s closing/opening base pair is in a continuous stack with at least one nucleotide of the second helix’s opening/closing base pair. This allows for any number of stacking nucleotides between the stems (not necessarily connected via the backbone) and for bulged out nucleotides which do not contribute to the stack. Our definition of stacking is more relaxed than stricter criteria used elsewhere²⁹.

For each pair of adjacent stems within a junction, we used forgi to calculate a number of properties, such as the angle between the stem vectors, the separation vector between the stems’ ends and an offset calculated as distance between two rays that start at the helix end and extend the helix axis.

The following code is a simplified example how forgi can be used to calculate such properties:

from forgi import load_rna
from forgi.threedee.utilities.vector import vec_angle

def main():
    # Load the PDB into CoarseGrainRNA instances
    # rnas is a list, because a PDB can contain multiple connected components
    # This also loads the DSSR json annotations.
    rnas = load_rna("Path/to/PDB/file.cif", pdb_remove_pk=False,
                    pdb_annotation_tool="DSSR")

    for rna in rnas:
        for junction in rna.junctions:
           print_junction_angles(rna, junction)

def print_junction_angles(rna, junction):
    """
    :param rna: A BulgeGraph Object
    :param junction: A list of element names (strings)
    """

    for ml in junction :
        # rna.edges holds adjacent nodes in the BulgeGraph
        stem1, stem2 = rna.edges[ml]
        # rna.coords holds the 3D coordinates of structure elements
       # rna.coords ["s0"] returns a tuple start –, end–coordinates of stem "s0"
       # Furthermore rna.coords has the get_direction function to give the
       # vector pointing from the start to the end coordinate.
        d1 = rna.coords.get_direction(stem1)
       d2 = rna.coords.get_direction(stem2)
       angle = vec_angle (d1, d2) # imported above
        is_stacking_dssr = (ml in rna.dssr.stacking_loops())
        print("{}\t{}".format(angle, is_stacking_dssr))

We then used pandas^26,27, Matplotlib³⁰ and a custom library (https://github.com/Bernhard10/filterAndView) to analyze and visualize the collected data. The results were collected for different classes of RNA independently. This was necessary, because the representative sets of RNA structures contain, by design, homologs of the same molecule in multiple species.

For the classification of pseudo knots, Reidys’ concept³¹ based on the definition of the mathematical genus is used. Classes of pseudo knots are defined based on their shadow representations, which contain only crossing base pairs and only one base pair each. On the level of these shadows, only four distinct classes exhibit genus 1, two of which are well known: the H-type pseudo knot and the kissing hairpin. pseudo knots with higher genus contain, among others, the case where a genus 1 pseudo knot is nested within another pseudo knot.

In combination with the forgi library, we wrote a tool that is able to convert structures to their shadow representation, identify and classify the pseudo knots within the structure and describe their helix arrangement in the 3D structure. This tool, called pseudoknot_analyzer.py, is distributed with the forgi library in the folder “examples”. We used it to gather statistics about simple H-type and kissing hairpin pseudo knots. Figure 4a, b illustrates how we measured the angle between stems in pseudo knots via vector directions. In kissing hairpins the angles α and β restrict the possible values of the angle γ. We also include the representative structure of an intermolecular kissing hairpin interaction (lacking the green connection in Figure 4b) in our analysis. In this special case α and β are indistinguishable and were assigned arbitrarily.

Results

We used the helix-centered representation of the 3D structure to analyze the geometry of coaxially stacking helices. The relative geometry of two cylinders in 3D space can be described by five parameters: A separation vector (three parameters) and two angles for the relative orientation in 3D space. In Figure 2 and Figure 4, we show the single angle calculated between the vectors along the helix axes, as a proxy for these parameters.

Figure 2. Distribution of angles between adjacent stems in multi loops.

(a) Angles close to 0° mean parallel stems whereas angles close to 180° correspond to potentially stacking stems. Angle distributions in (b) tRNAs and (c) ribosomal RNA are shown as a histogram mapped onto a circular plot illustrating the angles. The (inner) blue bars are instances where DSSR detects stacking on the atom-level scale. The orange bars start at the top of the blue bars and indicate geometries where DSSR does not detect stacking.

Figure 3. Example of a tRNA where stacking between two non-collinear stems occurs (PDB id 4WJ4³²).

Green cylinders were fitted to stems, blue cylinders represent hairpins and the pink cylinder the unpaired nucleotides at the 3’ end. Red cylinders connect stems and indicate single-stranded connections between these stems (multi loop segments). Left: The stems of the anticodon arm (top) and the D-arm stack (according to DSSR) despite being at an angle of 143°. Right: View along the axis of the anticodon arm’s stem. The red nucleotide is the unpaired multi loop segment which mediates stacking to the stem of the D-arm. This illustration was generated using PyMol³³ via the visualize_rna.py wrapper in forgi.

Figure 4. Angles between the stems in pseudo knots.

(a) and (b) show the vector directions used for the angle measurement. The distributions of the angles between adjacent stems building c) an H-type pseudo knot or (d) a kissing hairpin are shown as histograms mapped to a circle. Furthermore, the distribution of angles between the regular stems of kissing hairpin pseudo knots is shown in the second panel of (d). Like mentioned in Figure 2, the angles in c) measured between stacking stems (detected via DSSR) are colored in blue. Geometries without stacking are colored in orange. The color scheme in (d) (shades of orange, blue and green) refer to the respective vectors/measured angles (α, β, γ) in the same color scheme as in (b). Blue bars stand for stacking between the two related stems, whereas orange ones for non-stacking stems. Dark green bars in the second kissing hairpin associated histogram show angles measured between a intramolecular interaction (kissing hairpin pseudo knot), whereas light green bars stand for angles measured between two RNA chains (kissing hairpin interaction).

Figure 5.

Examples of coaxial stacking within (a) an H-type pseudo knot (PDB id 2XD0) and (b-d) the three major structural families of kissing hairpins (PDB ids 5KPY, 4FRN and 1ZCI, from left to right) In all four representations the green cylinders were fitted to stems, the turquoise and pink cylinder represent the unpaired nucleotides at the 5’ and 3’ end respectively. Red cylinders connect stems and indicate single-stranded connections between these stems (pseudo knot or multi loop segments). The light-colored stems in (b)–(d) represent the regular stems of kissing hairpins, whereas the kissing stem is colored yellow. Note that panel (d) (PDB id 1ZCI) shows a kissing hairpin interaction between two RNA chains. Below the 3-dimensional representation of the kissing hairpins, we show a schematic sketch of the structural family’s helix arrangement. These illustrations were generated using PyMol³³ via the visualize_rna.py wrapper in forgi.

The distribution of angles between adjacent stems in tRNAs multi loops (see Figure 2a) shows the expected bimodal distribution with one mode slightly below 90° and another mode between 160° and 170°, which fits to the known helix arrangement in the L-shaped tRNA. As confirmed by comparison to annotations with DSSR, the second peak is almost exclusively due to stacking helices, even at angles as small as 140°. In Figure 3 we show an example of such a coaxial stack that has been strongly bent (possibly by the tRNA synthetase) without completely breaking the stack or affecting the canonical tRNA secondary structure.

The second family of RNA molecules with lots of data available is ribosomal RNA. Here we find that the angle of adjacent stems in multi loops is almost uniformly distributed above 20° until a peak from 135° to 170°, where DSSR annotates roughly half of the geometries as stacking (see Figure 2c). 7% of the stacking geometries and 67% of the non-stacking geometries with an angle above 140° also have an offset above 10Å between the extended stem axes.

Additionally we analyzed the angle distribution between the stems forming simple H-type pseudo knots or kissing hairpins. Most angles measured within an H-type pseudo knot are between 120° and 180° (see Figure 4c). In this range, 35 out of 67 instances correspond to stacking. One example within this range is shown in Figure 5a, which represents a processed non-coding RNA, which regulates a bacterial antiviral system (PDB id 2XD0³⁴).

The distribution illustrated in Figure 4d shows mainly values between 130° and 180° for the angles β and especially α. One additional angle measurement between the two regular stems (see Figure 4d, angle γ) shows one peak at about 20° and one at about 65°. Within the class of kissing hairpins we often find coaxial stacking, but most of the time only one of the two regular stems stacks with the kissing stem in the middle. With the help of the coarse grained representation we were able to divide the class of kissing hairpins into three different structural families (see Figure 5b-d).

The first family is especially common among (A-)riboswitches. Here the two regular stems are oriented almost parallel (γ below to 32°) with the second one (counting from the 5’ end) stacking onto the kissing stem. The angles α and β are above 130°. One example is the structure with PDB id 5KPY³⁵ shown in Figure 5b. Here, the link from the first regular stem stacks onto it and interacts with the kissing stem via base multiplets and A-Minor interactions. This way, the roughly parallel orientation of all 3 stems is stabilized by stacking and base-pairing interactions. Interestingly, the kissing stem is more parallel to the first stem than to the second stem, onto which it stacks directly (α > β).

The second structural family is shown in Figure 5c. Here the two regular stems are close to 90°, whereas the kissing stem is along the arc between them. This conformation is found in several groups of non coding RNA including some 23S ribosomal RNA and the cobalamin riboswitch regulatory element (PDB id 4FRN³⁶) shown in Figure 5c.

The two families as well as some conformations in between are observable within intramolecular pseudo knots. As mentioned, forgi is also able to model multiple RNA chains that interact with each other forming an intermolecular complex. One example is the dimerization initiation site of the HIV type 1 (PDB id 1ZCI³⁷), illustrated in Figure 5d and being the only example of the third family of kissing hairpin pseudo knots. Here the kissing hairpin interaction shows a nearly perfect coaxial stacking between all three involved stems. In this case β is close to 0°, because stacking takes place at the other side of the kissing stem and thus the complementary angle is close to 180°.

Underlying data

The PDB ids analyzed were taken from version 3.36 of the representative sets of RNA 3D structures²⁸ at a cutoff of 4 Å, http://rna.bgsu.edu/rna3dhub/nrlist/release/3.36/4.0A.

Extended data

Open Science Framework: 3D based on 2D: forgi 2.0 Extended Data. https://doi.org/10.17605/OSF.IO/HDJRU⁴³. This project contains the following extended data files:

Extended data are available under the terms of the Creative Commons Zero "No rights reserved" data waiver (CC0 1.0 Public domain dedication).