Saturday, September 26, 2009

RiboClub 10th Annual Meeting -- it was a good one!

I attended the RiboClub 10th Annual Meeting held during September 21 to 23 at Hotel Cheribourg, Orford, Quebec. Everyday, the schedule was fulfilled from early morning to late night. Indeed, it was so tight and intensive that I could not even find a time to walk around in the beautiful season. Overall, though, the meeting was well-organized and had a fantastic scientific program.

The RiboClub was founded ten years ago by researchers at the University of Sherbrooke, Quebec. Initially, it was local, and then the club was extended to nearby areas, and the whole Canada. As time goes, its influence has also passed the border so at its 10th anniversary, several leading RNA experts (including Phillip Sharp, Jack Szostak, Tim Nilsen, Tom Steitz et al.) from the USA also participated.

I noticed the RiboClub meeting early this year when I was writing up a manuscript on an RNA structural motif. As hinted in another post, "Does 3DNA work for RNA?", I've recently been attracted to the field of RNA structures. Using 3DNA, we have uncovered a simple RNA-specific interaction that is biologically relevant, yet virtually ignored by the community. So attending a RNA meeting would allow me to learn more about RNA and to pass our message across to a wider audience.

The meeting was mostly on various aspects of RNA biology. Of which, 3-dimensional structure is an integral part, yet no specific session was devoted to it. The same applied to bioinformatics tools and applications. Thus, for example, Tom Steitz's talk on the ribosome structures was under the session titled "Translation: targets and impact".

The organizers obviously paid attention to arrange meeting participants at the dinning table. So for Tuesday (Sept. 22) night, I sit next to Dr. Andrew MacMillan from University of Alberta. It was a nice surprise to know that Dr. MacMillan works on "structural and functional characterization of splicesome assembly and activation." I took this opportunity to read the abstracts from his lab and to talk to him about my findings that are related to pre-mRNA splicing. He visited my post the next day (Wednesday, Sept 23) and we discussed it in more details. On the Gala dinner on Wednesday, I was arranged to sit between Dr. Paul Griffiths, "a philosopher of science with a focus on biology and psychology", and Dr. W. Ford Doolittle, a leading scientist in Comparative Genomics. It was a valuable experience to hear them and others around the table talking on politics and science-related issues.

It is worth noting that Wednesday's dinner speaker was Alexander Rich. Wearing his tie from the famous RNA Tie Club, Dr. Rich talked about "The era of RNA awakening: structural biology of RNA in the early years." At the age of 85, he still spoke clearly and logically. His story telling style was very effective and his talk was well-received by the audience. For those who are interested in knowing more about Dr. Rich's work, I would strongly recommend his article "The excitement of discovery."

The "Neo-Traditional Quebec Music" show on Wednesday night was exciting and relaxing, following and in contrast to the three-day long intensive scientific program. Although I did not understand the music that well, I stayed until the very end.

Saturday, September 5, 2009

Double helix groove width parameters from 3DNA

In the 3DNA output (from the analyze program) for a DNA/RNA duplex structure, there is a section on "Minor and major groove widths: direct P-P distances and refined P-P distances which take into account the directions of the sugar-phosphate backbones". The underlying algorithm is that of El Hassan and Calladine (1998). ``Two Distinct Modes of Protein-induced Bending in DNA.'' J. Mol. Biol., v282, pp331-343. Note that the P-P distances need to be subtracted by 5.8 Å to take account of the vdw radii of the phosphate groups (2.9 Å), and for comparisons with NewHelix/FreeHelix and Curves.

Using 3DNA fiber models #1 for A-DNA (calf thymus) and #4 for B-DNA (calf thymus), the groove widths are as follows:
                 Minor Groove        Major Groove
P-P Refined P-P Refined
-----------------------------------------------------
A-DNA (#1) 18.5 16.7 15.2 11.1
B-DNA (#4) 11.7 11.7 17.2 17.2
-----------------------------------------------------
From the above table, it is clearly that for A-DNA, the minor and major groove widths for the refined set are smaller than their corresponding non-refined counterparts (i.e., those based on direct P-P distances). For B-DNA, there are no changes between the two sets. It should be noted that in real structures (i.e., non-perfectly regular, as in X-ray crystal structures in the NDB/PDB), there are nearly always some differences between the refined vs. direct P-P distances. As a general rule, the refined set should be used.

One of the key structural differences between A- and B-DNA is their opposite groove dimensions: for B-DNA, the major groove width (~17 Å) is about 5 Å wider than the minor groove width (~12 Å); whereas for A-DNA, the major groove width (~11 Å) is narrower than the minor groove width (~17 Å) by a similar amount. Since the grooves provide binding sites, the difference between A- and B-DNA grooves has important implications in DNA (groove) recognitions by ligands or proteins.

In retrospect, I implemented the El Hassan and Calladine algorithm for calculating the groove widths mainly because of its simplicity: I can understand clearly how it works visually. The algorithm is described in a two-page appendix of the above cited paper. For those who are interest in DNA structures in general and how groove widths are defined in particular, I would strongly recommend them to read the appendix carefully and try to implement it: there is no substitute for first hand experience. For an idealized cases, as the above for fiber A- and B-DNA, the implementation should be straightforward. To be more realistic, an implementation should account for missing phosphate groups in some structures (for testing purpose, simply delete one P atom from a structure), for example.

As is obvious, 3DNA does not calculate groove depths. Over the years, I have actually been approached with requests/suggestions to provide such parameters to complement groove widths. However, for various reasons, none of the algorithms fits with 3DNA. As a general principle, I do not add new functionality to 3DNA simply for the seek of it. I must understand a new piece clearly in order to integrate it with the rest and to be able to respond concretely to possible questions from users.

Some emacs tricks

As a Linux/Unix fan, I am very familiar with vi and use it for quick and simple text editing purposes. Over the years, however, I have been using emacs the most: I like its color coding and programming language-specific editing mode.

Emacs is well-known for its extensibility, so there are many ways to customize it to suit one's taste. In my experience, I have found the following settings convenient:
  • Highlight the current line with a background color (here "greenyellow")
    (require 'highlight-current-line)
    (highlight-current-line-on t)
    (highlight-current-line-set-bg-color "greenyellow")
  • Set transient mark mode on, so that selected text become more obvious
    (setq transient-mark-mode t)
  • Show line number and column number
    (setq line-number-mode t)
    (setq column-number-mode t)
There are many features in emacs that could be handy and I am always trying to learn more of it.

Sunday, August 30, 2009

JMB celebrates 50 years of protein structure determination

In the September 11, 2009 issue (vol.392, issue 1) of JMB, there are a series of three articles on the determination of the first two protein structures (myoglobin and haemoglobin), an achievement accomplished by Perutz and Kendrew and their colleagues at the Cambridge MRC laboratory in 1950s. Of special interest of this series is the three authors — Bror Strandberg, Richard Dickerson and Michael Rossmann — leading scientists in structure biology, then postdocs actively involved in the late stage of the structure determination.

These reviews are vividly written and provide interesting background information and some technical details on X-ray crystal structure determination (especially on phase angle), and have the following titles:
  1. "Building the Ground for the First Two Protein Structures: Myoglobin and Haemoglobin" by Bror Strandberg
  2. "Myoglobin: A Whale of a Structure!" by Richard Dickerson
  3. "Recollection of the Events Leading to the Discovery of the Structure of haemoglobin" by Michael Rossmann
With limited computing power and software support, protein structure determination was a difficult task at that time. Dickerson and Rossmann had to write their own programs to perform some calculations. However, the firsthand experience in both experiment and code development, on a significant project in a famous lab, may in part accounts for their success in structure biology.

I have known Dickerson's work on nucleic acid structures for a while, firstly through the famous Drew-Dickerson dodecamer (CGCGAATTCGCG), and I am intimately familiar with his NewHelix/FreeHelix programs. Nevertheless, it is only after reading his above article do I become aware of his initial protein experience. I like Dickerson's writing a lot. For example, on commenting the different styles of Kendrew and Perutz, he wrote: "John was the mentor, guide, and organizer. …… In contrast, Max was a hands-on bench biochemist whose center of gravity was always the laboratory itself. …… Both styles had their merits: one learned from John, but one learned with Max."

An interesting point from Rossmann's article is his description of a secret he kept to himself for many decades: "I [Rossmann] had been privileged to work on the haemoglobin project with Max, but it was also a project that Max had given his whole life to develop. In my enthusiasm to look at the results, I stole the final discovery from Max. …… With the realization of what I had done, all desire to explore further was completely gone." While I vaguely remembered this story from reading the book "Max Perutz and the Secret of Life" by Georgina Ferry several months ago, Rossmann's personal account would make it unforgettable.

It is worth noting that Perutz and Rossmann were among those few who initiated the PDB in 1971 at a Cold Spring Harbor meeting. Finally, given the expertise of the three authors, it is not surprising to read in the Epilogue that "Indeed, structural biology has become the unifying factor of just about every aspect of biology."

Saturday, August 22, 2009

Fit a least squares plane to a set of points

In the early days when I first got into the field of DNA structures, I was interested in knowing details on how the "best mean (base-pair) plane" is actually defined. I could not find a solution from the articles I read at that time. For clarification, I sent emails to some authors whose papers mentioned the concept, still the results were less than satisfactory.

The mystery was finally solved after I read the following two articles from early days (traced initially from a comment in the source code of Dickerson's NewHelix, if I recall it correctly):
  • D. M. Blow (1960). ``To Fit a Plane to a Set of Points by Least Squares.'' Acta Cryst., 13, 168.
  • V. Schomaker et al. (1959). ``To Fit a Plane or a Line to a Set of Points by Least Squares.'' Acta Cryst., 12, 600-604.
As is clear from the titles, they focused just on least-squares plane/line fitting algorithms. More specifically, the worked example(s) help make the underlying concept clear.

As a concrete example, let's work out the least-square plane of an adenine using Octave/Matlab. The implementation here is much simpler than those in the literature cited above. Note: this example is mentioned in the 3DNA home page in the Technical Details section, so the order of atoms is kept the same as there.

The x-, y-, and z-coordinates are as follows:
xyz = [
16.461 17.015 14.676 % 91 N9
15.775 18.188 14.459 % 100 C4
14.489 18.449 14.756 % 99 N3
14.171 19.699 14.406 % 98 C2
14.933 20.644 13.839 % 97 N1
16.223 20.352 13.555 % 95 C6
16.984 21.297 12.994 % 96 N6
16.683 19.056 13.875 % 94 C5
17.918 18.439 13.718 % 93 N7
17.734 17.239 14.207 % 92 C8
];

Cmtx = cov(xyz)
[V, LAMBDA] = eig(Cmtx)

The covariance matrix (Cmtx), its corresponding eigenvectors (V) and eigenvalues (LAMBDA) are as follows:
Cmtx =
1.66800 -0.50154 -0.32530
-0.50154 2.06703 -0.58401
-0.32530 -0.58401 0.30605

V =
0.27367 0.83264 -0.48147
0.32243 0.39220 0.86152
0.90617 -0.39102 -0.16113

LAMBDA =
0.00001 0.00000 0.00000
0.00000 1.58452 0.00000
0.00000 0.00000 2.45656
The smallest eigenvalue of Cmtx, 0.00001 is close to zeros, indicating that the base is almost perfectly planar. The corresponding unit eigenvector, i.e., the least squares base normal is: [0.27367 0.32243 0.90617]. Of course, the geometric average of the atoms defines a point the ls-plane pass through.

Referring back to the concept of "best mean (base-pair) plane", it could simply be defined using all base atoms of the pair, or alternatively, as the mean of the two base normals. Obviously, they would lead to (slight) numerical differences. Astute viewers would also notice some other subtleties.

Sunday, August 16, 2009

Curves+ vs 3DNA

While browsing Nucleic Acids Research recently, I noticed the paper titled "Conformational analysis of nucleic acids revisited: Curves+" by Dr. Lavery et al. [advance access published online on July 22, 2009]. I read it through carefully during the weekend and played around with the software. Overall, I was fairly impressed, and also happy to see that "It [Curves+] adopts the generally accepted reference frame for nucleic acid bases and no longer shows any significant difference with analysis programs such as 3DNA for intra- or inter-base pair parameters."

Anyone who has ever worked on nucleic acid structures (especially DNA) should be familiar with Curves, an analysis program that has been widely used over the past twenty years. Only in recent years has 3DNA become popular. By and large, though, it is my opinion that 3DNA and Curves are constructive competitors in nucleic acid structure analysis with complementary functionality. As I put it six years ago, before the 13th Conversation at Albany: "Curves has special features that 3DNA does not want to repeat/compete (e.g. global parameters, groove dimension parameters). Nevertheless, we provide an option in a 3DNA utility program (find_pair) to generate input to Curves directly from a PDB data file" on June 6, 2003, and emphasized again on June 09, 2003: "We also see Curves unique in defining global parameters, bending analysis and groove dimensions." 3DNA's real strength, as demonstrated in our 2008 Nature Protocols paper, lies in its integrated approach that combines nucleic acid structure analysis, rebuilding, and visualization into a single software package (see image to the right and above).

Now the nucleic acid structure community is blessed with the new Curves+, which "is algorithmically simpler and computationally much faster than the earlier Curves approach", yet still provides its 'hallmark' curvilinear axis and "a full analysis of groove widths and depth". When I read the text, I especially liked the INTRODUCTION section, which provides a nice summary of relevant background information on nucleic acid conformational analysis. An important feature of Curves+ is its integration of the analysis of molecular dynamics trajectories. In contrast, 3DNA lacks direct support in this area (even though I know of such applications from questions posted on the 3DNA forum), mostly due to the fact that I am not an 'energetic' person. Of special note is a policy-related advantage Curves+ has over 3DNA: Curves+ is distributed freely, and with source code available. On the other hand, due to Rutgers' license constraints and various other (undocumented) reasons, 3DNA users are still having difficulty in accessing 3DNA v2.0 I compiled several months ago!

It is worth noting that the major differences in slide (+0.47 Å) and X-displacement (+0.77 Å) in Curves+ vs the old Curves (~0.5 Å and ~0.8 Å, respectively) are nearly exactly those uncovered a decade ago in "Resolving the discrepancies among nucleic acid conformational analyses" [Lu and Olson, J Mol Biol. 1999 Jan 29; 285(4):1563-75]:
Except for Curves, which defines the local frame in terms of the canonical B-DNA fiber structure (Leslie et al., 1980), the base origins are roughly coincident in the different schemes, but are significantly displaced (~0.8 Å along the positive x-axis) from the Curves reference. As illustrated below, this offset gives rise to systematic discrepancies of ~0.5 Å in slide and ~0.8 Å in global x-displacement in Curves compared with other programs, and also contributes to differences in rise at kinked steps. (p. 1566)

Please note that Curves+ has introduced new name list variables — most notably, lib= — and other subtle format changes, thus rendering the find_pair generated input files (with option '-c') no longer valid. However, it would be easy to manually edit the input file to make it work for Curves+, since the most significant part — i.e., specifying paired nucleotides — does not change. Given time and upon user request, however, I would consider to write a new script to automate the process.

Overall, it is to the user community's advantage to have both 3DNA and Curves+ or a choice between the two programs, and I am more than willing to build a bridge between them to make users' lives easier.

Saturday, August 15, 2009

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 pab of a score matrix": it clarifies every details for those who want to get to the bottom of the topic.

Sunday, August 9, 2009

On third-party interfaces to 3DNA

One clear sign of 3DNA's acceptance by the nucleic-acid structure-related scientific community is its integration into a wide range of GUI/web-interfaces (e.g., NDB, w3DNA, 3D-DART etc) and many other packages. As the author and maintainer of the 3DNA software package, I am (of course) very happy to see this happen: I can easily imagine that better accessibility would make 3DNA open to an even larger community, including non-experts in nucleic acid structures (e.g., for educational purpose). It should be emphasized, however, that point-and-click interfaces to 3DNA (and any other serious scientific software, for that matter), while convenient to occasional users and for routine jobs, have their limitations.

Overall, a point-and-click interface is only sensible to well-defined, routine tasks. In the context of 3DNA, the following list of tasks (not necessarily complete) is perfectly suitable:
  • Generate one of the 55 fiber models, where everything (model number, sequence or number of repeats) can be unambiguously defined through a web-form
  • Build an arbitrary DNA model with a user-specific parameter set
  • Create a blocview image for a nucleic-acid containing structure specified in PDB format
  • Analyze a regular (i.e., not deviating too much) double helix structure, be it in A-, B-, or Z-form
  • Calculate a list of backbone torsion angles and other parameters (with the "-s" option of find_pair and then analyze)
  • Find a list of all possible (RNA) base-pairs fulfilling a specific set of geometric criteria (with the "-p" option of find_pair).
Other than those (and carefully selected/qualified other functionality), users should be very careful in what they get from running a third-party tool outputting 3DNA-related parameters. I know of the so many subtleties in using 3DNA to solve a specific problem, and most of the time the command-line driven style is the most effective, as it allows for easy automation of try-and-error iterations. Software developers who integrate 3DNA into their packages should make the limitations clear to their users, and direct any 3DNA-specific problems to the 3DNA forum.

In today's informatics world, there are so many "easy-to-use" tools available, claiming to be able to solve all sorts of problems (well, that's understandable -- otherwise, how could one get published, especially in the big journals/magazines?). Any serious scientist, however, should know what he/she is doing: it is easy to get some (magic) numbers by clicking a button, but understand clearly what the numbers mean is yet another story. I cannot emphasize more the importance of knowing one's tools, including their limitations.

Saturday, August 8, 2009

Prefatory articles in Annual Review of Biochemistry make good reading

Over the time, I have read several (mostly structural biology related) prefatory articles published in the Annual Review of Biochemistry and found them quite interesting. Some examples:
  • In Vol. 78 (July 2009), James Wang's article, titled "A Journey in the World of DNA Rings and Beyond", told a story behind the discovery of DNA topoisomerase, an enzyme that converts one form of DNA ring to another.
  • In Vol. 73 (July 2004), Alexander Rich wrote about "The Excitement of Discovery" as a scientist. Rich's research has often been "on the question how molecular structure leads to biological function". I am especially impressed by his vivid description of the work with Crick on the structure of collagen (p.10-12) and its "strong positive effect" on his psyche:
    For one thing, I began to develop some self-assurance in my ability to carry out research and make discoveries. I believe that a form of “scientific maturation” is an important component in developing a confident thrust into research work.
  • In Vol. 40 (July 1971), John Edsall provided interesting insights about his editorial work in his article titled "Some Personal History and Reflections from the Life of a Biochemist".
Overall, such articles are well worth reading -- they provide background information and the context on the discoveries made by those leading scientists.