What's special about the GpU dinucleotide platform?

Recently, I (together with Drs. Wilma Olson and Harmen Bussemaker – a team with a unique combination of complementary expertise) published a new article in Nucleic Acids Research (NAR): "The RNA backbone plays a crucial role in mediating the intrinsic stability of the GpU dinucleotide platform and the GpUpA/GpA miniduplex". The key findings of this work are summarized in the abstract:

The side-by-side interactions of nucleobases contribute to the organization of RNA, forming the planar building blocks of helices and mediating chain folding. Dinucleotide platforms, formed by side-by-side pairing of adjacent bases, frequently anchor helices against loops. Surprisingly, GpU steps account for over half of the dinucleotide platforms observed in RNA-containing structures. Why GpU should stand out from other dinucleotides in this respect is not clear from the single well-characterized H-bond found between the guanine N2 and the uracil O4 groups. Here, we describe how an RNA-specific H-bond between O2'(G) and O2P(U) adds to the stability of the GpU platform. Moreover, we show how this pair of oxygen atoms forms an out-of-plane backbone ‘edge’ that is specifically recognized by a non-adjacent guanine in over 90% of the cases, leading to the formation of an asymmetric miniduplex consisting of ‘complementary’ GpUpA and GpA subunits. Together, these five nucleotides constitute the conserved core of the well-known loop-E motif. The backbone-mediated intrinsic stabilities of the GpU dinucleotide platform and the GpUpA/GpA miniduplex plausibly underlie observed evolutionary constraints on base identity. We propose that they may also provide a reason for the extreme conservation of GpU observed at most 5'-splice sites.

As a nice surprise, this publication was selected by NAR as a featured article! According to the NAR website:

Featured Articles highlight the best papers published in NAR. These articles are chosen by the Executive Editors on the recommendation of Editorial Board Members and Referees. They represent the top 5% of papers in terms of originality, significance and scientific excellence.

I feel very gratified with the "extra" recognition. From my own perspective, I can easily rank this paper as the top one in my publication list: from the very beginning, I has been struck by the simplicity and elegance of the GpU story. Hopefully, time will verify the validity of this scientific contribution.

Behind the hood, though, there is a long, complex (sometimes perplexing), yet interesting story associated with this work. Here is how it got started. While writing the 3DNA 2008 Nature Protocols (NP) paper, I selected the (previously undocumented) "-p" option of "find_pair" to showcase its capability to identify higher-order base associations, using the large ribosomal subunit (1JJ2) as an example. I noticed the unexpected O2'(G)⋅⋅⋅O2P(U) H-bond within the GpU dinucleotide platform in the pentaplet shown left in Figure A below. I was well aware of Leontis-Westholf's pioneering work on "Geometric nomenclature and classification of RNA base pairs" which involves three distinct edges – the Watson-Crick edge, the Hoogsteen edge, and the Sugar edge, yet without taking into consideration of possible sugar-phosphate backbone interactions (Figure B below). So I decided to double-check, just to be sure that the H-bond was not spurious due to defects in the H-bond detecting scheme of "find_pair", and the results were very surprising.

The following section was re-added into the 3DNA NP paper in the very last revision:

It is also worth noting that the G1971–U1972 platform is stabilized not only by the well-characterized G(N2)⋅⋅⋅U(O4) H-bond interaction, but also by a little-noticed G(O2’)⋅⋅⋅U(O2P) sugar-phosphate backbone interaction (Fig. 6a). Examination of the 50S large ribosomal unit (1JJ2) alone reveals ten such double H-bonded G–U platforms, far more occurrences than those registered by any other dinucleotide platform (including A–A) in this structure. Apparently, the G–U platform is more stable than other platforms with only a single base–base H-bond interaction. We are currently investigating this overrepresented G–U dinucleotide platform in other RNA structures. (p.1226)

What find_pair in 3DNA can do

Structural analysis of nucleic acids used to be a rather tedious process, especially for irregular, complicated RNA structures and nucleic-acid/protein complexes [e.g., the large ribosomal subunit of H. marismortui (1JJ2)]. Without valid base-pairing information arranged properly in a duplex fragment as input, analysis programs such as Curves+ and analyze/cehs in 3DNA would produce meaningless results. The program find_pair in 3DNA was originally created to solve this specific problem, i.e., to generate an input file to 3DNA analysis routines directly from a nucleic-acid containing structure in PDB format. It is what makes nucleic acids structural analysis a routine process — running through thousands of structures from NDB/PDB can be fully automated.

Overall, find_pair has more than fulfilled the goal of its initial design (as stated above). Over the past few years, its functionality has been expanded and continuously refined (kaizen; 改善), making find_pair itself a full-featured application. Now, it is efficient, robust, and its simple command line interface allows for easy integration with other bioinformatics tools. Properly acknowledged or otherwise, find_pair has served (at least) as one of the key components in many other applications (RNAView, BPS, SwS, ARTS, to name just a few). Indeed, find_pair is by far the single program in 3DNA that has received the most questions (as evident from the 3DNA forum).

While I still have to write a method paper to describe the underlying algorithms of find_pair in detail — i.e., for identifying nucleotides, H-bonds, base pairs, high-order base associations, and double helical regions — the basic idea is very intuitive and easy to understand: as summarized in our recent GpU paper, find_pair is purely geometric based (with user adjustable parameters) and allows for the identification of canonical Watson–Crick as well as non-canonical base pairs, made up of normal or modified bases, regardless of tautomeric or protonation state. For example, in the GpU paper, we chose the following set of stringent parameters to ensure that the geometry of each identified base pair is nearly planar and supports at least one inter-base H-bond: (i) a vertical distance (stagger) between base planes ≤ 1.5 Å; (ii) an angle between base normal vectors ≤ 30°; and (iii) a pair of nitrogen and/or oxygen base atoms at a distance ≤ 3.3 Å. Other criteria (documented or otherwise), such as the distance between the origins of the two standard base reference frames, are just filters to speed up the calculations.

In a nutshell, find_pair has the following two core functionalities:

1. The default is to generate input to the analysis routines in 3DNA (analyze/cehs) for double helices. However, there are many more works under the hood than just identifying base pairs: the base pairs must be in proper sequential order, and each strand must be in 5' to 3' direction, for the calculated step parameters (twist, roll etc) to make sense. Moreover, with the "-c" option, one gets an input file to Curves (but not Curves+, yet); with the "-s" or "-1" option, find_pair treats the whole structure as one single strand, and is useful for getting all backbone torsion angles.


2. Detect all base pairs (regardless of in double helical regions or not) and higher-oder (3+) base associations with the "-p" option. This feature (in its preliminary form) was there starting from at least v1.5, which was released at the end of 2002 (just before I left Rutgers), but it was intentionally not documented. The source code of find_pair (as part of 3DNA) was tested and shared within Rutgers (NDB and Dr. Olson's laboratory) before any 3DNA paper was published, and served as the basis for several other projects. We also offered 3DNA (with source code) to a few RNA experts for comments; but we received either no responses or politely-worded negative ones. Things did not work out as (what I thought) they should have been, but that's life and I have learned my lessons. The "-p" option was first explicitly mentioned in the 3DNA 2008 Nature Protocols paper, to illustrate how to identify the two pentaplets in the large ribosomal subunit of H. marismortui (1JJ2).

It is interesting to mention the two papers I've recently come across: the first is on DNA-protein interactions and the second on RNA base-pairing, where new algorithms were developed to detect base pairs and their performances were compared with find_pair. In each of the two cases, it was claimed that find_pair missed certain pairs where the new methods succeeded. As it turned out, however, in the first case, simply relaxing find_pair's default H-bond distance cut-off 4.0 Å to 4.5 Å, as used by the authors, virtually all the missing pairs were recovered. In the second case, the "-p" option, which should have been, was simply not specified.

After nearly a decade of extensive real-world applications and refinements, it is safe to say that find_pair is now a versatile and practical tool for nucleic acids structure analysis. Of course, I will continue to support and further refine find_pair as I see fit. Once in a while, I just cannot stop but to think that find_pair is to nucleic acids what DSSP is to proteins: simple and elegant. As more people become aware of its existence, I would expect find_pair to gain even more widespread usage, especially in RNA-structure related research areas.

Eddy's primer on the BLOSUM62 alignment score matrix

As a Linux/Unix fun, I like its philosophy, as summarized by Doug McIlroy, very much: "Write programs that do one thing and do it well. Write programs to work together." In science, I enjoy more reading an article that focuses on one point and address it clearly and thoroughly. It is only after a complete understanding of the components can it be possible to combine them in unique, purpose-specific ways. Those days, though, such type of simple-and-clean articles is no longer that common as it was in the early days, say 1960s or 70s. This post is the first of a series on such articles I've found useful, or on computational tricks that I have learned over the years.

I came across the primer titled "Where did the BLOSUM62 alignment score matrix come from?" by Sean Eddy [Nat Biotechnol. 2004 Aug;22(8):1035-6] early this month through the BioConductor mailing list where it was recommended by Dr. Philipp Pagel as a "well written article". I read the title and the short abstract from PubMed, and then download the PDF version of the whole article.

As a primer, it is only two pages long (short), and reading twice won't take that much time. It explains the meaning of the BLOSUM62 amino acid score matrix and where is comes from clearly. The number 62, for example, stands for a threshold of 62% identity. Other percentages, such as 80% (more Conservative), 45% (more divergent) are also possible and may be more suitable for specific applications. As noted by Eddy, "Empirically, the BLOSUM matrices have performed very well. BLOSUM62 has become a de facto standard for many protein alignment programs."

With the above background, it is much easier to understand the fundamental difference of the default scoring systems between FASTA/WU-BLASTN vs NCBI BLASTN for alignment of DNA sequences: the former is optimal for alignments at the ‘twilight zone’ (65% identity), while the later (NCBI BLASTN) is optimal for homologous DNA at a much higher 95% identity level. Knowing of such subtlety is very important in avoid making false conclusions.

More significantly (to me), there is a supplemental material -- a well-documented, self-contained "ANSI C program for calculating the implicit target frequencies p_ab of a score matrix": it clarifies every details for those who want to get to the bottom of the topic.

Meaning of nucleotide IUPAC codes

Today, anyone with some basic knowledge of biochemistry should be familiar with A, C, G and T, the four bases of DNA, and probably the A-T and G-C Watson-Crick base-pairs as well. The meaning of the one-letter abbreviations is very clear: A for Adenine, C for Cytosine, G for Guanine, and T for Thymine. Of course, for RNA, there is the U (for Uracil) in place of T of DNA.

In the early days when I entered into the field of DNA structure, I also learned that R stands for puRine, i.e., A and G, and Y for pYrimidine, i.e., C and T (U). Trained as a chemist, I had no difficult at all in understanding and remembering them. To process base sequences in bioinformatics projects, I have come across the IUPAC degeneracy codes of nucleotides, such as S, W, M, K D, V, etc, which I had never been able to really memorize what they represent, except for N (A, C, G, T).

My confusions have been clarified completely, however, due to the web document I happened to find: "Nomenclature for Incompletely Specified Bases in Nucleic Acid Sequences" (1984) by the Nomenclature Committee of the International Union of Biochemistry (NC-IUB). This is the document I wish I could have known of from the very beginning. For completeness of this post, here is a summary table of the whole IUPAC codes. It is based on Table 1 in the above document except for uracil and gap:

Symbol Meaning          Origin of designation
-----------------------------------------------------------
G      G                Guanine
A      A                Adenine
T      T                Thymine
C      C                Cytosine
U      U                Uracil
R    G or A             puRine
Y    T or C             pYrimidine
M    A or C             aMino
K    G or T             Keto
S    G or C             Strong interaction (3 H bonds)
W    A or T             Weak interaction (2 H bonds)
H    A or C or T        not-G, H follows G in the alphabet
B    G or T or C        not-A, B follows A
V    G or C or A        not-T (not-U), V follows U
D    G or A or T        not-C, D follows C
N    G or A or T or C   aNy
    . or -             gap
-----------------------------------------------------------

Errors in PDB entries

In the June 24, 2009 issue of Nature (v459, pp.1038-1039), there is an news item titled "New protein structures replace the old" by Katharine Sanderson, on a 'Dutch software to weed out errors in Protein Data Bank'.

In my experience with software development and using the PDB/NDB, it is certainly not a surprise that there are errors of various types in the macromolecular databases: whenever I apply an algorithm consistently to all the entries in the NDB (which is part of PDB, consisting of only nucleic-acid containing structures), I always notice some inconsistencies. As a more concrete example, blocview, a visualization tool initially developed as a by-product of another project while I was still at Rutgers (and partially involved with the NDB), was once used for correcting errors in the NDB as well.

Pure 're-refinement' of existing structures with software is surely helpful in catching obvious, systemic errors. However, it is impossible to catch all problems, no matter how sophisticated the software could be. Moreover, as put in the comment by yet another phd: "i have hard enough time getting the RCSB to change four atoms in a structure for me." and "scientists must remember that when they click the re-refine button, you read the paper where the structure was reported."

Errors will always be there -- that's just a basic fact of life. It is thus crucial for those who perform structural analysis to draw their conclusions based on not one or just a few purposely selected structures, but on a more objective and extensive ground.

PDB ATOM coordinates record

PDB format is one of the standard formats for biological macromolecular structures (proteins, DNA/RNA, their complexes, etc). It came into existence when the initial Brookhaven PDB was established in 1971. Over the years, PDBML and mmCIF have been proposed as substitutes to allow for more flexibilities, yet PDB format is still the mostly commonly used one. Software dealing with PDB structures each has its own parser, at least for the Coordinate Section (especially the ATOM and HETATM records).

Overall, PDB format is simple and very well documented. The simplicity lies just in its 'rigidity', in FORTRAN 77 style. The ATOM/HETATM record description is excerpted below for easy reference:

COLUMNS        DATA  TYPE    FIELD        DEFINITION
-------------------------------------------------------------------------------------
1 -  6         Record name   "ATOM  "
7 - 11         Integer       serial       Atom  serial number.
13 - 16        Atom          name         Atom name.
17             Character     altLoc       Alternate location indicator.
18 - 20        Residue name  resName      Residue name.
22             Character     chainID      Chain identifier.
23 - 26        Integer       resSeq       Residue sequence number.
27             AChar         iCode        Code for insertion of residues.
31 - 38        Real(8.3)     x            Orthogonal coordinates for X in Angstroms.
39 - 46        Real(8.3)     y            Orthogonal coordinates for Y in Angstroms.
47 - 54        Real(8.3)     z            Orthogonal coordinates for Z in Angstroms.
55 - 60        Real(6.2)     occupancy    Occupancy.
61 - 66        Real(6.2)     tempFactor   Temperature  factor.
77 - 78        LString(2)    element      Element symbol, right-justified.
79 - 80        LString(2)    charge       Charge  on the atom.

It won't take much time/lines in a script language such as Perl/Python/Ruby etc to extract specific information one is interested in, e.g., atomic coordinates. However, there are some subtleties that are beyond simple script parsers. On top of that, one needs to understand that not all self-claimed PDB files are standard compliant.

More specifically, a decent PDB format parser must take the following into considerations:

• The four-character atom name specified in columns 13 to 16. Each biological molecule has a convention in naming atoms. For example, the two H-bonds of the A-T pair are between " N1 " (A) to " N3 " (T), and " N6 " (A) to " O4 " (T). In this regard, it worths noting that babel/openbabel converted PDB files do not follow such naming convention.
• The one-character alternate location indicator (altLoc) in column 17
• The one-character residue insertion code (iCode) in column 27.
• Others details to follow ...

If you are serious about your PDB format parser, there is a very simple thing to do: run it against all the NDB entries if you are working on nucleic acid structures (that's what I did for 3DNA), and the whole PDB entries if you are interested in protein structures. If it does not crash and does what it has been designed for, then congratulate yourself!