Sunday, July 5, 2009

Two web-interfaces to 3DNA, and more

As a nice surprise, I found in the 2009 web-server issue of NAR published on July 1, two articles on web-interface to 3DNA back-to-back:
  1. "3D-DART: a DNA structure modelling server" by van Dijk and Bonvin from Utrecht University, The Netherlands. Excerpt from the abstract:
    As a response to the demand for 3D-structural models reflecting the intrinsic plasticity of DNA we present the 3D-DART server (3DNA-Driven DNA Analysis and Rebuilding Tool). The server provides an easy interface to a powerful collection of tools for the generation of DNA-structural models in custom conformations. The computational engine beyond the server makes use of the 3DNA software suite together with a collection of home-written python scripts.
  2. "Web 3DNA—a web server for the analysis, reconstruction, and visualization of three-dimensional nucleic-acid structures" by Zheng, me and Olson. Excerpt from the abstract:
    The w3DNA (web 3DNA) server is a user-friendly web-based interface to the 3DNA suite of programs for the analysis, reconstruction, and visualization of three-dimensional (3D) nucleic-acid-containing structures, including their complexes with proteins and other ligands.
While I was aware of 3D-DART prior to its publication, I certainly did not expect it to appear in the same issue as w3DNA of which I am a co-author. There is nothing more compelling to illustrate 3DNA's value to the community than a third-party web-interface to it! Combined together, these two web-servers make 3DNA much more accessible to even wider audience. Specifically, they could well serve for education purposes, e.g., to conveniently build a DNA-model of A-, B- or C-form, with user-supplied sequences. Users would be glad to have a choice that better fits their needs. In the long run, the one which provides best user-support will survive.

It is also worthy noting the in the same 2009 NAR web-server issue, another paper titled "SARA: a server for function annotation of RNA structures" by Capriotti and Marti-Renom from Spain also makes use of 3DNA. This serves to emphasize the point that 3DNA is not just for DNA, but for RNA as well -- 3DNA has unique features for RNA that are not found in other currently available software tools that I am aware of.

Friday, June 19, 2009

RasMol 2.6.4 and 3DNA base rectangular block schematic representation

Over the years, I have been using RasMol v2.6.4, the last release in December 1998 by Roger Sayle, the original author of this molecular graphics visualization tool. According to the README file, RasMol v2.6.4 "includes mainly minor bug fixes and portability enhancements" to the version Sayle used internally.

Back in 1990s and early 2000s, RasMol was the most popular tool of its kind. However, even though RasMol v2.7.x series is actively maintained, in recent years JMol, PyMol and other similar tools have surpassed RasMol in popularity. So why still RasMol v2.6.4? To me, RasMol is a light-weight, robust tool for molecular structure visualization, just like ls, grep Unix commands for my daily computing life. Moreover, RasMol v2.6.4 is connected to 3DNA with the visualization of Alchemy format for rectangular base (or base-pair) block representation.

It is straightforward to install RasMol v2.6.4 in a modern computer system with a C compiler. You may also want to customize it, and more importantly, get to know how to use it effectively.
  • RasMol v2.6.4 can be downloaded from http://www.dcs.ed.ac.uk/home/rasmol/v2.6beta/RasMol2.tar.gz
  • In Linux, type "xmkmf; make" to compile, and this gives executable rasmol. By default, this is for 8-bit graphics. I normally renamed it to rasmol8 by "mv rasmol rasmol8" because the 8-bit version is useful for generating gif images (for web graphics, for example). To make a 32-bit version as default, edit Makefile macro DEPTHDEF: "# DEPTHDEF = -DEIGHTBIT" to "DEPTHDEF = -DTHIRTYTWOBIT", then "make clean; make". Move "rasmol8" and "rasmol" to a directory in your command line search path (e.g., ~/bin), copy "rasmol.hlp" to some location and set up environmental variable RASMOLPATH accordingly, you are all set.
  • By default, rasmol binary does not allow for outputting graphics in (automatic) scripting mode, this can be overwritten by changing AllowWrite = False; at line 1163 of "rasmol.c" (in function InitDefaultValues()) to AllowWrite = True;. While you are at it, you may also want to change other defaults. Just remember to "make clean; make".
  • It is worthwhile to study the help file to get family with its contents, especially the syntax for selections. The real power of RasMol comes with its scripting capability. It can be used to automate an otherwise tedious operation.
In relation to 3DNA (and SCHNArP for that matter), RasMol v2.6.4 can be used to visualize the Calladine-Drew style rectangular representation of bases (base-pairs) in Alchemy format. As a matter of fact, I initially picked up Alchemy format in SCHNArP (and later in 3DNA) primarily because it is one of the structure file formats supported by RasMol 2.6.x. JMol does not support 3DNA generated Alchemy format file until recently when I made a request, and PyMol does not support Alchemy format the last time I tried it.

I once upgraded to v2.7.x series, but noticed immediately a problem with visualization of Alchemy files from 3DNA. Initially, I thought it was a bug in 3DNA introduced in the source code somewhere. After some late night debugging, I re-tried to display the Alchemy in RasMol v2.6.4, and found that the fault was due to v2.7.x! It is actually quite easy to verify this v2.7.x bug: simply load a PDB structure, write to alchemy, then re-load the alchemy file, and you will see it. I reported this bug to the maintainer of v2.7.x, but was told that the bug was in v2.6.4. That's the reason I have stayed with v2.6.4 thereafter, and made the point that v2.7.x does not work with 3DNA generated schematic. As my programming experience grows, I realize more clearly the subtlety of a software product, and would normally trust the original author unless proved otherwise. It is far too easy to introduce a bug somewhere while fixing/improving other parts, unless you know every detail and/or have a through regression test suite for quality control.

In visualizing 3DNA generated Alchemy files with RasMol v2.6.4, it is important to remember to use command line options "-alchemy -noconnect". While the former option is obvious, the later one is to ensure that RasMol respects the connections explicitly specified by 3DNA generated Alchemy files, which are, after all, schematic, not real atomic structures. Sometimes (not always), this option would make a difference.

Saturday, June 13, 2009

blocview: a simple, effective visualization tool for nucleic acid structures

The 'blocview' Perl script in 3DNA was designed as a simple tool that nevertheless illustrates key features of nucleic acid structures effectively. Specifically, the 'bloc' part of the name means 'block', i.e., the base rectangular block in the Calladine-Drew style to show clearly the size (larger purine vs. smaller pyrimidine), identity (by color: red for A, yellow for C, green for G, and blue for T), and the groove (with minor groove edge filled in black); and the 'view' part means the most extended view (as defined by the principal axes of inertia). 'blocview' calls several utility programs of 3DNA and MolScript (for protein ribbons and nucleic acid backbone rods) to prepare the scenes, and then uses Raster3D (more specifically 'render') or PyMol (as of 3DNA v2.0) to generate a PNG image.

It has been my pleasure to notice that 'blocview' generated images are used in the NDB for each and every structure. See for example, an arbitrarily picked link on X-ray Drug-DNA Complexes. A nice thing about these images is that they can be generated automatically: using PDB id 1z8v as an example, the following command
blocview -i=1z8v.png 1z8v.pdb
would produce the image (named 1z8v.png above) shown in the right. In this representation, one can see clearly that there are two unpaired Gs (green) at each of the 5' end of the two DNA chains (red and yellow rods), and a drug molecule (ball-and-stick) binds in the minor groove (black edge of the rectangular blocks). Moreover, the first C-G pair has pronounced negative propeller twist.

It was a nice surprise when I found (several years ago) that such simple images have also be adopted by the PDB, prominently at the summary page, for every nucleic acid containing structure. There are so many sophisticated molecular graphics tools out there, why 'blocview'? After all, it is just a small utility tool in the 3DNA suite of programs, and I have never written a paper on it. Clearly, 'blocview' fills a niche.

On the other hand, given the effectiveness of this simple representation and its adoption by the NDB and PDB, I am surprised to found that 'blocview' images are not used that much in the literature. Maybe the situation will change as 3DNA gets more widely used, and when people pay more attention to its versatile functionalities: certainly, there are more handy features in 3DNA!

More blocview command-line options are available to suit some other common needs: just type 'blocview -h' for details. Furthermore, many aspects of the images can be fine-tuned by configuration files. Astute viewers could have noticed subtle differences of the 1z8v images shown here (above right), and in the NDB and PDB websites. Since 'blocview' is written in Perl script, users can easily check the implementation details to figure out exactly how it works.

Friday, June 12, 2009

Review of scientific manuscripts

While reading an editorial of Maddox [Nature. 1995 Dec 7; 378(6557):521-3], I found the following comment quite interesting:

It's a good joke (which I have often used) that Watson and Crick's paper on the structure of DNA could not be published now. If is only necessary to imagine what people would say if it reached them in the mail: "It's all model-building, just speculation, and such data they have are not theirs but Rosalind Franklin's!" Some would complain that the sentence beginning, "It has not escaped our attention..." is certainly unsubstantiated, and must be an attempt to claim credit for developments in genetics that lies years ahead.

Maddox continued to imagine what could have happened behind the scene that led to the publication of Watson and Crick's paper, together with Franklin's crucial experimental work illustrating the helical signature of crystalline DNA.

In some sense (to my understanding), Watson and Crick were lucky, not only in that they solved the puzzle of DNA structure elegantly by connecting the various experimental pieces together, but also they worked in the Cavendish Laboratory at Cambridge University and had the support from Sir Lawrence Bragg. One could easily imagine what would have happened if they did this work in a not so famous lab, without support from a leading, powerful figure.

In a scientific field (e.g. RNA), what appears to be a common sense, a decade old issue, could have something significant when view from a novel perspective. Such unconventional results, however, are often hard to get across the review process, especially for new comers to a field and without backing of big names. In principle, the reviewers should be fair to the author, being specific with his/her comments, especially in case of negative ones. That would be far more convincing to the authors than some vague, general remarks. As a general rule, I won't accept to review a manuscript for a journal unless I have time to read it thorough carefully, and the expertise to understand the details, in order to provide some concrete comments.

In this regard, it is worth reading the essay by Keith Manchester titled "Historical Opinion: Erwin Chargaff and his 'rules' for the base composition of DNA: why did he fail to see the possibility of complementarity?" [Trends Biochem Sci. 2008 Feb;33(2):65-70]. For completeness, the abstract is quoted below:
Erwin Chargaff was one of the more interesting and colourful figures of the historic decade that heralded the proposal of the double helical structure of DNA by Watson and Crick in 1953. In describing Chargaff's important contribution to the study of DNA, particularly its base composition, this article seeks to suggest why, despite his substantial achievements, he failed to anticipate some of the key features of the Watson-Crick model, particularly complementarity between bases--a failure that left him deeply embittered for the rest of his life.

Saturday, June 6, 2009

ANSI C program memory check with valgrind

Over the years, ANSI C has been my selected language for scientific computations, from SCHNAaP/SCHNArP, 3DNA, to my current research projects. As noted by K&R in the preface of "The C Programming Language" book (2nd edition), C "wears well as one's experience with it grows". To me, ANSI C is a language that the better I know it, the more I love it.

Using gcc with strict compiler options -ansi -pedantic -Werror -Wall -W -Wshadow, my ANSI C programs (e.g., 3DNA) compile cleanly without any change on all the OSes I have access to. Any C programmer should be well aware of the memory leak type of bugs that are so difficult to fix, since it is not checked by gcc. That's where valgrind comes in and I have been using it in 3DNA and my current projects on a regular basis.

The usefulness of the valgrind toolset becomes obvious when a SCHNArP user posted a question about some possible memory leak bugs. Traditionally, SCHNArP works interactively from command line and builds one model (DNA structure) at a time. When the program exists, any meomry leaks are reclaimed by the system so such bugs do not show up. At the time SCHNAaP/SCHNArP was developed, I was not aware of valgrind (or maybe it even did not exist yet!). The user has been using some functions from SCHNArP to build many structures automatically, thus any small amount of memory leaks scales up and eventually aborts out-of-memory. As a test case, the user supplied me with a reproducible example as follows:
    while(1){
GLH_build();
}
Without knowing of valgrind and actually applied it in this situation, fixing such bugs would have been too time consuming than I could manage to spend. However, with the help of valgrind, solving this puzzle did not take that long at all. More specifically, recompiling the sample program with -g option to have debug information, and running it with valgrind --leak-check=full prefixed immediately revealed a couple of locations with memory leaks. Knowing where the problems are, fixing them became straightforward. Now SCHNArP is more robust.

Overall, valgrind is an essential tool for any serious C programers.

Friday, June 5, 2009

3DNA citations reach over 300

When checking Web of Science today, I found, to my great satisfaction, that citations to the 3DNA 2003 NAR paper have reached over 300. While the number is still not comparable to true classics, it is nevertheless a significant one. It is worth noting that the 300+ citations span 77 scientific journals in the broad biomedical and chemical research fields, including Cell, Nature and Science (well, if CNS really means anything). The top six journals (with over double-digit citations) are:
  • 49 -- NUCLEIC ACIDS RESEARCH
  • 27 -- JOURNAL OF MOLECULAR BIOLOGY
  • 23 -- JOURNAL OF THE AMERICAN CHEMICAL SOCIETY
  • 19 -- BIOCHEMISTRY
  • 13 -- JOURNAL OF PHYSICAL CHEMISTRY B
  • 11 -- BIOPHYSICAL JOURNAL
In retrospect, the wide acceptance of 3DNA by the community is certainly not by chance. Listed below are some of the major factors:
  1. 3DNA is built upon the best selected features from a thorough understanding of seven popular DNA structural analysis programs, including Curves, NewHelix/FreeHelix, CEHS, CompDNA, and RNA (Running Nucleic Acids by Babcock and Olson). This is summarized back in June 2003, nearly five years ago, in an email I sent to a group of experts before a Workshop at the 13th Conversation at Albany (NY):
    Looking back, 3DNA has never been intended to set up a new standard for the exact definition of various parameters. We selected what made best sense to us from what were already available. This was made possible by knowing the exact technical details of how other programs work. In the meantime, 3DNA also provide a program called "cehs" which gives the original CEHS/SCHNAaP parameters, to which Freehelix parameters would be similar. Believe it or not, "cehs -r" would give the authentic RNA parameters. We also see Curves unique in defining global parameters, bending analysis and groove dimensions."
  2. Active support through emails and the 3DNA forum. Over the years, I have always striven to get back to each 3DNA-related question quickly and concretely. In my memory, I have never ignored users' questions. Instead, I have taken each and every question as a way to improve 3DNA, and refined the code accordingly.
  3. 3DNA's integrated approach that incorporates structural analysis, modeling building, and visualization in a single suite of software package. Our 2008 Nature Protocols paper provides some examples.
  4. 3DNA has been extensively checked against all NDB entries before each major release, so that I can confidently say that it works in real world, not just in some contrived limited example cases.
I still remember that the number of citations to 3DNA was less than 150 nearly two years ago, when I started to wrote the first draft of our 2008 Nature Protocols paper. Now it is more than doubled! I would blog on this topic again when the number reaches 500.

In my opinion, some of 3DNA features are still (heavily) underused. Now that we have a sizable user community, 3DNA could only become better and would be more widely used. I have every reason to believe that in the not-so-distant-future, the citations to 3DNA would reach over 1000.

For current status of 3DNA citations, check Google Scholar.

Monday, May 25, 2009

Examples make it clear

In reading scientific publications, I always quickly get confused with long, complicated equations. Initially, I thought it was because of my lack of domain knowledge: when I get more familiar with a field, things will become easier. Of course, this statement is true. However, it is not exactly the point in the following two senses:
  1. Having been in the field of nucleic acid structures for more than ten years, I can claim I know the field good enough. Yet, I can still easily get lost with large portions of mathematical formula in an article. An example is the Babcock and Olson 1994 JMB paper titled "Nucleic acid structure analysis. Mathematics for local Cartesian and helical structure parameters that are truly comparable between structures". Overall, it is a solid piece of work but I doubt how many people in this field really understand it above the conceptual level.
  2. When I get to the bottom of something, things do become truly simple. But even then the equations still feel complicated. In resolving the discrepancies among nucleic acid conformational analyses, I went through the source code to understanding from the bottom up how the Babcock and Olson method works. It is actually not that complicated, and the implementation with regard to finding the local helical axis and calculating a set of helical parameters can be greatly simplified using a geometric approach, as in SCHNAaP/3DNA. Moreover, some details from the above Babcock and Olson 1994 JMB paper are actually inaccurate.
Another wonderful case I have had is with the CEHS scheme (due to Calladine and El Hassan) to calculate a set of local base-pair parameters. It was only after checking an example did I uncover a subtle discrepancy between the simple CEHS formula (in its original form) and its description, in the draft version of Dr. El Hassan's PhD thesis. This experience led to the joint SCHNAaP/SCHNArP papers with El Hassan, and inspired me later on to work through NewHelix, Curves, CompDNA etc DNA structure analysis programs to reveal "the reasons why the different algorithms come to conflicting structural interpretations".

So, if you are a beginner in a research field, do not be scared off by the fancy mathematical formula. If you want to really understand something, work through an example and/or check through the source code implementations. For the authors, if you would like to get your message across, provide step-by-step examples to the key points you want to make for others to follow.

Saturday, May 16, 2009

Fiber, analyze and rebuild in 3DNA

While browsing the current issue of Nucleic Acids Research (May 2009, Vol. 37, No. 8), I noticed, surprisingly, the article titled Reconstitution of ‘floral quartets’ in vitro involving class B and class E floral homeotic proteins by Rainer Melzer and Günter Theißen from Jena, Germany.

From the title and abstract, I would not expect 3DNA could play a role here. Yet in section Modeling a floral quartet, the authors used fiber to building a B-DNA model of floral quartets, using sequence between the XhoI and XbaI sites of probe A, then analyze the structure, and modified the resultant file 'bp_step.par' accordingly to their specifications, and rebuild a customized DNA model. In the following paragraph, the authors used the analyze/rebuild pair again for the CArG box bound to serum response factor.

This is the intended usage of the 3DNA software package I have had in mind from its initial design. I emphasized such versatile, integrated approach unique with 3DNA in our 2008 Nature Protocols paper. As time goes by, I am confident that more people will begin to appreciate this approach.

Friday, May 15, 2009

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!