Xiang-Jun's Corner

3DNA base-mutation functionality for the study of pretein-DNA interactions

2011-10-30T23:12:00.000-04:00

In my last blog post titled "mutate_bases, a new 3DNA tool for in silico base mutations of 3D nucleic acid structures", I outlined the availability of the mutate_bases program and some areas of its possible applications, including to "perform base-pair mutations in DNA-protein complexes".

In the September 6, 2011 issue of PNAS (vol. 108, no. 36), AlQuraishi and McAdamsa from Stanford University published a high-profile article "Direct inference of protein–DNA interactions using compressed sensing methods" (pp. 14819–14824). See also the thoughtful Commentary by Vijay Pande, "(Compressed) sensing and sensibility". Essentially, by combining concepts from compressed sensing and statistical mechanics, AlQuraishi and McAdamsa have developed a novel approach to determine the energy potential of protein–DNA complexes, leading to "an impressive advance in predictive capability" -- ~90% accuracy compared with ~60% for the best-performing alternative computational methods.

I am so glad to find that 3DNA (ref nos. 12 and 13) was used in the work:

In silico mutagenesis was carried out using the 3DNA software package (12, 13), which maintains the backbone atoms of the DNA molecule, but replaces the base pair atoms in a way that is consistent with the backbone orientation in the crystal. (p.14821, middle of the right column)

Hopefully, the new and novel base mutation functionality in 3DNA will find more applications and be in wider use. Surely, I'll respond promptly to users' feedback and refine related 3DNA tools as necessary.

mutate_bases, a new 3DNA tool for in silico base mutations of 3D nucleic acid structures

2011-10-25T21:45:00.002-04:00

In response to repeated requests from 3DNA users over the years, I recently (well, a few months ago) released a new utility program, mutate_bases, to the 3DNA suite of software package. As its name implies, mutate_bases can be used for mutating bases in a nucleic acid structure, with two key and unique features: (1) the sugar-phosphate backbone conformation is untouched; (2) the base reference frame (position and orientation) is reserved, i.e., the mutated structure shares the same base-pair/step parameters as the original one.

mutate_bases is a standalone ANSI C program, on a par with other major 3DNA programs (e.g., find_pair, analyze, rebuild, fiber etc). As seen from the help message below, it can be used for any nucleic-acid-containing structures (DNA, RNA, or their complexes) in PDB format:

------------------------------------------------------------------------
Usage: mutate_bases mutinfo pdbfile outfile
    'mutinfo' can contain upto 5 fields for each mutation:
        [name=residue_name] [icode=insertion_code]
        chain=chain_id seqnum=residue_number
        mutation=residue_name
    Alternatively, 'mutinfo' can be specified with the '-l' option
        followed by a mutation file: '-l mutations.dat'
    o The five fields per mutation can be in any order or CaSe
    o Each field can be abbreviated to its first character
    o Multiple mutations are separated by ';'
    o Fields in [] (i.e., name and icode) are optional
    o Mutation info should be QUOTED to be taken as one entry

Examples:
    mutate_bases 'c=a s=2 m=DA' 355d.pdb 355d_G2A.pdb
    # mutate G2 in chain A of B-DNA 355d to A
    mutate_bases 'c=a s=2 m=DA; c=B s=23 m=DT' 355d.pdb 355d_mutfile.pdb
    # mutate base-pair G-C (C23 in chain B) to A-T
        # also create file 'mutations.dat', see below
    mutate_bases -l mutations.dat 355d.pdb 355d_mutfile.pdb
        # 355d_mutfile.pdb and 355d_mutfile.pdb would be identical
    mutate_bases 'c=A s=74 m=U' 1evv.pdb 1evv_C74U.pdb
    # mutate C74 in chain A of tRNA 1evv to U
------------------------------------------------------------------------

mutate_bases is designed to solve the base mutation problem in a practical sense: robust and efficient, getting its job done and then out of the way. The program can have many possible applications: in addition to perform base-pair mutations in DNA-protein complexes, it should also prove handy RNA modeling and in providing initial structures for QM/MM/MD energy calculations, and in DNA/RNA modeling studies.

Tracking 3DNA citations: Google scholar vs. Web of Science

2011-07-29T23:46:00.001-04:00

Over the years, I have been following citations to 3DNA on a regular basis, using both Google Scholar and Web of Science. From a personal perspective, following the citations has turned out to be an excellent way to keep myself informed of progress related to nucleic acid structures.

Up to early this year, I had found 3DNA citations based on Google Scholar to be significantly larger than those from Web of Science. For example, for the 2003 Nucleic Acids Research (NAR) paper, the difference was well over one hundred. Over the past few months, however, I have noticed that citation numbers to the 2003 NAR paper fluctuating up and down in dozens. Moreover, I have been receiving significantly less email notifications to 3DNA citations from Google Scholar. Apparently, Google has been revising the algorithms for Google Scholar, leading to more conservative citation numbers that are now comparable to those from Web of Science.

Here are the current citation numbers to the three 3DNA publications, from both Google Scholar and Web of Science:

	Google Scholar	Web of Science
2003 NAR	544	473
2008 Nature Protocols	43	41
2009 NAR (w3DNA)	11	13

It is interesting to note that for the w3DNA paper, Web of Science gives a larger number (13) than Google Scholar. A careful check of the citation result from Web of Science shows that all of the 13 citations are legitimate journal articles. Clearly, Web of Science is missing something in this case.

Overall, I trust Web of Science more, which is why I have been using it to compile citations to the 2003 NAR paper. However, I use Google Scholar more frequently for quick reference simply because it is free and easy to access.

GpU dinucleotide platform, the smallest unit with key RNA structural features

2011-07-14T23:18:00.000-04:00

Compared to DNA, RNA has three salient structural features: it contains ribose sugar, uracil, and is normally single-standed. The O2'(G)...O2P(U) H-bond stabilized GpU dinucleotide platform may turn out to be the smallest unit with all those RNA hallmarks (see Figure below).

Firstly, it must have the guanosine ribose to form the O2'(G)...O2P(U) H-bond.

Secondly, the methyl group in position 5 of thymine would cause steric clash with guanosine, thus disrupting the N2(G)...O4(U) base-base H-bond to form the GpU dinucleotide platform.

Thirdly, a dinucleotide, by definition, is single-standed. The two H-bonds, plus the covalent linkage, makes the GpU platform extremely rigid (see Figure 1 of our 2010 NAR paper).

Moreover, the GpU platform is directional: swapping the two bases while keeping the sugar-phosphate backbone fixed does not allow for a base-base H-bond, thus no UpG dinucleotide platform.

PDB v3.3 and partial atom occupancy

2011-06-26T00:06:00.000-04:00

From the PDB mailing list, I know of the recent announcement "PDB Archive Version 4.0 to be Released July 13, 2011". This "ambitious" review of the PDB archive has resulted in a new set of corrected files in ten categories, including biological assemblies, residual B factors, peptide inhibitors/antibiotics, polymers containing nonstandard polymer linkages, and partial occupancy etc. "These data reflect the wwPDB's continuing commitment to providing accurate and detailed data to users worldwide."

I am interested in changes in the PDB format, and read the PDF document "Description of Changes and Corrections for PDB July 2011 Remediation Release":

PDB format files updated in this remediation release comply with PDB Format Version 3.30. PDBx and PDBML data files comply with the PDB Exchange Dictionary v.4.0, and PDBML XML Schema V4.0, respectively.

Specifically, I checked carefully the section on "Partial occupancy", which is quoted in full below:

Problem

In the 2009 remediation, occupancies were corrected in 490 X-ray
and neutron entries. A mistake was made in 104 of these entries:
for atoms with alternate conformer labels and with summed total
occupancy less than 1.0, the occupancies were re-scaled as 1.0/n,
where n is the number of conformers.

Approach

The originally deposited occupancies of the affected atoms were
restored and the remediation was then carried out properly, via:
Atoms with multiple conformations but identical coordinates and B-values were merged and their occupancies were summed.
Atoms which now have (total) occupancies <= 1.0 were left as deposited.
Atoms with (total) occupancies > 1 were rescaled proportionally to a sum of 1.0
Results

The occupancies have been corrected in these entries.

This partial occupancy issue reminds me of an extensive and very informative discussion a few months back in CCP4BB, under the thread "what to do with disordered side chains" and its derivatives, about setting "zero" occupancy and/or high B values for PDB ATOM/HETATM records in disordered regions. Over the course of the discussion, Frances Bernstein made the following comment:

I am absolutely positive that there is software that does its voodoo on ATOM/HETATM records and pays absolutely no attention to anything beyond the x, y, z coordinates (i.e. beyond column 54).

3DNA does pay attention beyond column 54 (up to 80, actually) for ATOM/HETATM records, but internally it does not make use of the occupancy info. In future releases of 3DNA, I am planning to take occupancy/B-factor into consideration, probably through configurable parameters.

Reading through the "Description of Changes and Corrections for PDB July 2011 Remediation Release", I noticed cases of re-corrections of previous corrections in wwPDB remediation efforts. A concrete example is about partial occupancy in the 2009 remediation: among the 490 corrected X-ray and neutron entries, a mistake was made in 104 of them. As a side note, there was a post by Morten Kjeldgaard, titled "PDB changes data in entries?" early this year in CCP4BB, where partial occupancy was used as an example.

Sugar pucker correlates with phosphorus-related distance

2011-06-19T22:30:00.002-04:00

The sugar puckers in DNA/RNA structures are predominately in either C3'-endo or C2'-endo (see Figure below), corresponding to the A- or B-form conformation in a DNA duplex.

Recently, I (re-)read a few articles related to the RNA backbone by Jane Richardson et al., including

"RNA backbone is rotameric" (PNAS 2003)
"RNA backbone: consensus all-angle conformers and modular string nomenclature" (RNA 2008)
"MolProbity: all-atom structure validation for macromolecular crystallography" (Acta Crystallogr D Biol Crystallogr. 2010)
"PHENIX: a comprehensive Python-based system for macromolecular structure solution" (Acta Crystallogr D Biol Crystallogr. 2010)

I somehow became interested in the correlation between sugar pucker and a simple distance parameter, as reported in these papers:

C3'-endo and C2'-endo sugar puckers are highly correlated to the perpendicular distance between the C1'–N1/9 glycosidic bond vector and the following phosphate: > 2.9 Å for C3'-endo and < 2.9 Å for C2'-endo. (p.16 of the MolProbity paper)

Out of curiosity and to get a better understanding of this correlation, I played around with some sample cases both visually in RasMol and numerically. Overall, this is a simple geometric problem, i.e., the shortest distance from a point to a line in 3-dimensional space. Given below is the Octave/Matlab script for calculating the distances for G176 and U176 of PDB entry 1JJ2 (the large ribosomal subunit of Haloarcula marismortui):

function d = get_p3_nc_dist(P3, C1, N)
    N_C1 = N - C1;                 # vector from N to C1'
    nv_N_C1 = N_C1 / norm(N_C1);   # normalized vector
    C1_P3 = P3 - C1;               # vector from C1 to P3
    proj = dot(C1_P3, nv_N_C1);
    d  = norm(C1_P3 - proj * nv_N_C1);
end

## G175 (1jj2)
P3 = [70.104 112.366  44.586];
C1 = [73.017 109.666  45.304];
N = [74.445 109.380  45.288];
d1 = get_p3_nc_dist(P3, C1, N)    # 2.2 Å -- C2'-endo

## U176 (1jj2)
P3 = [66.871 116.402  46.804];
C1 = [68.213 112.454  49.279];
N = [69.678 112.480  49.438];
d2 = get_p3_nc_dist(P3, C1, N)    # 4.6 Å -- C3'-endo

The GpU used in the above example forms a dinucleotide platform, where the sugar of G175 adopts a C2'-endo conformation, and that of U176 has C3'-endo. Indeed, the distance for the G175 nucleotide is 2.2 Å, less than 2.9 Å; whilst the value for U176 is 4.6 Å, greater than 2.9 Å.

It is worth noting the above mentioned articles from Richardson et al. are focused on RNA backbone, without paying attention to the base (pair) geometry. The Zp parameter, which quantifies the z-coordinate of the phosphorus atom in the mean reference frame (see "A-form conformational motifs in ligand-bound DNA structures", JMB 2000), can be easily adapted to the analysis of single stranded RNA structures. For example, the vertical distances of the 3' phosphorus atoms to the G175 and U176 base planes are 1.9 Å and 4.4 Å, respectively.

Since base planes and the phosphorus atoms are the most accurately located entities in a given nucleic acid structure, the nucleotide-based Zp variant presumably would have some advantage over the distance from phosphorus to the glycosidic bond. Naturally, this Zp parameter will be added in future releases of 3DNA.

Conformation of the sugar ring in nucleic acid structures

2011-06-11T00:10:00.000-04:00

The conformation of the five-membered sugar ring in DNA/RNA structure can be characterized using the five corresponding endocyclic torsion angles (see Figure below).

i.e.,

v0: C4'-O4'-C1'-C2'
v1: O4'-C1'-C2'-C3'
v2: C1'-C2'-C3'-C4'
v3: C2'-C3'-C4'-O4'
v4: C3'-C4'-O4'-C1'

Due to the ring constraint, the conformation can be characterized approximately by 5 - 3 = 2 parameters. Using the concept of pseudorotation of the sugar ring, the two parameters are the amplitude (τ_m) and phase angle (P).

One set of widely used formula to convert the five torsion angles to the pseudorotation parameters is due to Altona & Sundaralingam: "Conformational Analysis of the Sugar Ring in Nucleosides and Nucleotides. A New Description Using the Concept of Pseudorotation" [J. Am. Chem. Soc., 1972, 94(23), pp 8205–8212]. As always, the concept is best illustrated with an example. Here I use the sugar ring of G4 (chain A) of the Dickerson-Drew dodecamer (1bna/bdl001), with Matlab/Octave code:

# xyz coordinates of the sugar ring: G4 (chain A), 1bna/bdl001
ATOM     63  C4'  DG A   4      21.393  16.960  18.505  1.00 53.00
ATOM     64  O4'  DG A   4      20.353  17.952  18.496  1.00 38.79
ATOM     65  C3'  DG A   4      21.264  16.229  17.176  1.00 56.72
ATOM     67  C2'  DG A   4      20.793  17.368  16.288  1.00 40.81
ATOM     68  C1'  DG A   4      19.716  17.901  17.218  1.00 30.52

# endocyclic torsion angles:
v0 = -26.7; v1 = 46.3; v2 = -47.1; v3 = 33.4; v4 = -4.4
Pconst = sin(pi/5) + sin(pi/2.5)  # 1.5388
P0 = atan2(v4 + v1 - v3 - v0, 2.0 * v2 * Pconst);  # 2.9034
tm = v2 / cos(P0);  # amplitude: 48.469
P = 180/pi * P0;  # phase angle: 166.35 [P + 360 if P0 < 0]

The Altona & Sundaralingam (1972) pseudorotation parameters are what have been adopted in 3DNA. The Curves+ program, however, uses another set of formula due to Westhof & Sundaralingam: "A Method for the Analysis of Puckering Disorder in Five-Membered Rings: The Relative Mobilities of Furanose and Proline Rings and Their Effects on Polynucleotide and Polypeptide Backbone Flexibility" [J. Am. Chem. Soc., 1983, 105(4), pp 970–976]. The two sets of formula by Altona & Sundaralingam (1972) and Westhof & Sundaralingam (1983) give slightly different numerical values for the two pseudorotation parameters (amplitude τ_mand phase angle P).

Since 3DNA and Curves+ are currently the most commonly used programs for conformational analysis of nucleic acid structures, the subtle differences in pseudorotation parameters may cause confusions for users who use both programs. With the same G4 (chain A, 1bna) sugar ring, here is the Matlab/Octave script showing how Curve+ calculates the pseudorotation parameters:

# xyz coordinates of sugar ring G4 (chain A, 1bna/bdl001)

# endocyclic torsion angles, same as above
v0 = -26.7; v1 = 46.3; v2 = -47.1; v3 = 33.4; v4 = -4.4

v = [v2, v3, v4, v0, v1]; # reorder them into vector v[]
A = 0; B = 0;
for i = 1:5
    t = 0.8 * pi * (i - 1);
    A += v(i) * cos(t);
    B += v(i) * sin(t);
end
A *= 0.4;   # -48.476
B *= -0.4;  # 11.516

tm = sqrt(A * A + B * B);  # 49.825

c = A/tm; s = B/tm;
P = atan2(s, c) * 180 / pi;  # 166.64

For this specific example, i.e., the sugar ring G4 (chain A, 1bna/bdl001), the pseudorotation parameters as calculated by 3DNA following Altona & Sundaralingam (1972) and Curves+ following Westhof & Sundaralingam (1983) are as follows:

         amplitude (τ_m)     phase angle (P)
3DNA        48.469             166.35
Curves+     49.825             166.64

Needless to say, the differences are subtle, and few people will notice/bother at all. For those who do care about such little details, however, hopefully this post will help you understand where the differences actually come from.

Lower case chain identifiers in PDB format

2011-06-05T00:38:00.000-04:00

First formulated in early 1970s, the PDB format is rigid with fixed columns for designated contents in its ATOM/HETATM records. Specificlly, a single column, #22, is assigned for the chain identifier (id). Traditionally, the 26 upper case letters of English alphabet (A-Z), space (i.e., ' '), and the single digits (0-9) have been used as chain ids. Up until the ribosomal structures came up, I guess, those 26 + 1 + 10 = 37 characters had been sufficient for the chain ids.

To the best of my knowledge, for a long time, most PDB parsers assume upper case chain ids. Indeed, 3DNA v1.5 automatically converts each ATOM/HETATM records to upper cases. The first time I became aware of lower case chain ids was when I saw a post in the 3DNA forum, titled "Small bug in find_pair", where a user reported the 'w' vs 'W' chain ids in PDB entry 1VSP. Then I refined 3DNA so that the case of chain ids can be preserved, through an undocumented command line option (as a feature for internal testing purpose).

My view to make 3DNA chain ids case-sensitive has been reinforced when I read the article "Crystal structures of CGG RNA repeats with implications for fragile X-associated tremor ataxia syndrome" recently published in Nucleic Acids Research. The asymmetric unit of the unmodified CGG-repeats-containing duplex (GCGGCGGC)₂, NDB entry NA1017 / PDB entry 3R1C, contains a total of 36 chains: designated as A-Z, plus a-j. Without distinguishing cases of the chain ids, the 3DNA output would become quite confusing.

Thus, in future releases of 3DNA, the default would be switched to preserve the case of chain ids. This chain id 'case' serves as an excellent example that scientific software products, unlike publications per se, are not fixed but need continuous care and maintenance to meet the challenges of an evolving world.

NAR's top ten articles

2011-05-22T16:01:00.000-04:00

Recently, I noticed a new feature in the website of Nucleic Acids Research (NAR), i.e., its selection of top ten articles:

NAR’s Top Ten Articles are updated monthly and show recent articles that have been most often accessed in HTML and PDF formats in the specified month.

In the age of information explosion with flood of scientific journals and articles, it is easy get lost. NAR's pick of top ten and featured articles draws my attention to significant work I may otherwise overlook.

The current top ten articles (March 2011) are all selected from 2009/2010 publications in 'Database', 'Methods online', and 'Survey and Summary'. I am browsing the 2009 article by Thomas LaFramboise, "Single nucleotide polymorphism arrays: a decade of biological, computational and technological advances", to get a better understanding of SNPs.

Posts in the 3DNA forum reach 600

2011-05-15T19:18:00.000-04:00

As of May 6, the total number of posts in the 3DNA forum has reached 600. Created in March 2007, with my debut post titled "Welcome message from Xiang-Jun Lu", the forum is now over four years old. Overall, the forum has served its purpose pretty well. In answering questions, I've been increasingly referring to the posts in the forum. As a concrete example, see the thread of a recent question "Base pair step parameters with a missing base pair".

At less than three posts (about one question) per week on average, I've not felt too much stress in supporting the forum (and maintaining 3DNA) in my spare time. For the most part, I've enjoyed interacting with 3DNA users from everywhere in the world, and with diverse backgrounds. Following the Unix philosophy ("Write programs that do one thing and do it well. Write programs to work together."), 3DNA has proved to be robust and flexible in serving its ever-growing user community. As a matter of fact, few questions I received a couple of years ago were beyond my original consideration of the details while I wrote the code. It is this intimated knowledge of all the underlying algorithms and every bit of their implementations that allows me to answer users' questions quickly and concretely.

As time passes by, however, it has become evident to me that 3DNA needs to be further refined and extended to meet the ever changing needs of its user community. For example, over the past few months, several questions asked in the 3DNA forum are directly relevant but clearly beyond 3DNA's current capabilities. While I'd be interested in implementing some of the requested functionality that make sense to me, doing so is certainly over my spare time limit. On the other hand, my increased understanding of nucleic acid structures and accumulated software expertise make it simply an issue of time and effort to move 3DNA to the next level, far beyond its current application scope and impact.

With posts in the 3DNA forum reaching 600, and citations to 3DNA articles over 600 (Google scholar), I am hopeful something good will happen to the 3DNA project. After all, 6 is a lucky number in traditional Chinese culture.

Fifty years of operon

2011-05-15T18:38:00.002-04:00

In the latest issue of Science, there is a one-page editorial titled "The Birth of the Operon" by François Jacob, who won the Nobel Prize in Physiology or Medicine in 1965:

What is the operon, whose 50th anniversary is being celebrated this week? The word heralded the discovery of how genes are turned on and off, and it launched the now-immense field of gene regulation. ... we cannot presume to know how new ideas will arise and where scientific research will lead.

In the next three paragraphs, Jacob provides an insightful and vivid description of his research related to the discovery of the "operon" – a structural gene-regulatory gene ensemble. In consonant with his comment on scientific discovery, he concludes:

Our breakthrough was the result of “night science”: a stumbling, wandering exploration of the natural world that relies on intuition as much as it does on the cold, orderly logic of “day science.” In today’s vastly expanded scientific enterprise, obsessed with impact factors and competition, we will need much more night science to unveil the many mysteries that remain about the workings of organisms.

It is worth noting that the Journal of Molecular Biology (JMB) has recently published a special issue [Volume 409, Issue 1, Pages 1-88 (27 May 2011)], titled "The Operon Model and its Impact on Modern Molecular Biology" with historical accounts and reviews to celebrate operon's 50th anniversary. It is because of this event that motivated me to read the Jacob and Monod 1961 JMB review article "Genetic regulatory mechanisms in the synthesis of proteins" – I have come across this paper so many times before, and should have definitely read it long ago!

Curves+ web server

2011-05-15T17:34:00.000-04:00

Through Google Scholar, I become aware of the article online in Nucleic Acids Research (NAR), titled "CURVES+ web server for analyzing and visualizing the helical, backbone and groove parameters of nucleic acid structures" by Richard Lavery's group:

Curves+, a revised version of the Curves software for analyzing the conformation of nucleic acid structures, is now available as a web server. This version, which can be freely accessed at http://gbio-pbil.ibcp.fr/cgi/Curves_plus/, allows the user to upload a nucleic acid structure file, choose the nucleotides to be analyzed and after optionally setting a number of input variables, view the numerical and graphic results online or download files containing a set of helical, backbone and groove parameters that fully describe the structure. PDB format files are also provided for offline visualization of the helical axis and groove geometry.

The website looks quite streamlined, with required input information all in a single page, and the test page also ran smoothly. In less than two years following the publication of Curves+, it is nice to see the Curves+ web server version available, making this analysis tool more readily available to the nucleic acids community.

Nowadays, it seems safe (to the best of my knowledge) to say that only 3DNA and Curves+ conform to the 1999 Tsukuba convention for the description of nucleic acid base-pair geometry, and each of them provides a web interface: web 3DNA and web Curves+.

Scientific journals on nucleic acids

2011-05-01T23:49:00.002-04:00

In my knowledge, Nucleic Acids Research (NAR) is a highly respected scientific journal with a broad impact in the field of nucleic acids. Over the years, I have been browsing NAR webpage on a regular basis to keep myself up to date to the latest development in this area. It is thus no surprise that the initial 3DNA paper was submitted to and published in NAR in 2003. Among the 500+ citations to that 3DNA paper, over 1/5 (100+) articles are from NAR itself (as an example, please see my January 22, 2011 blog post titled "Three structural biology papers in the latest issue of NAR cite 3DNA"). My latest contribution to NAR is the GpU story, which was actually selected as a featured article.

Another related journal I am quite familiar with is RNA, a publication of the RNA society. As the "About" section of its webpage succinctly summarizes,

RNA serves as an international forum for publishing original reports on RNA research in the broadest sense. The journal aims to unify this field by cutting across established disciplinary lines and focusing on "RNA-centered" science.

RNA currently has an impact factor (IF) of 5.198 (2009), slightly lower than NAR's 7.479. It is, nevertheless, a very decent journal in RNA-related research, and I frequently visit its website. As a side note, the GpU paper was initially submitted to RNA for its RNA-specific content and as a way to diversify my publication spectrum (as mentioned above, 3DNA was initially published in NAR). Unfortunately, the GpU paper was rejected by the RNA journal after two rounds of review, spanning over 6 months.

Another journal closely related to RNA (name wise) is called RNA Biology, which even has a slightly higher IF of 5.56. Admittedly, I was not familiar with this journal at all. Browsing through its website, I am interested in seeing the journal's explicit policy to reconsider papers "rejected by high impact journals [CNS] for reasons of novelty and impact, rather than the importance of the study or the integrity of the data." By enclosing "the reviewers’ and/or editorial comments" from these high impact journals, "it is possible the article might be accepted [by RNA Biology] based on its previous review. This will allow the urgent and competitive research to be published on the day of submission."

I became aware of the journal DNA Research quite recently through an email. From its website, "DNA Research is an internationally peer-reviewed journal which aims at publishing papers of highest quality in broad aspects of DNA and genome-related research." The journal currently has an IF of 4.917. Browsing a couple of its online issues, I sense that the journal is more on genome- than structure-related research.

While following up 3DNA citations recently, I noticed the paper titled "Insights into the Structures of DNA Damaged by Hydroxyl Radical: Crystal Structures of DNA Duplexes Containing 5-Formyluracil" by Tsunoda and Taknaka. It was published in the Journal of Nucleic Acids, which I have never (but probably should have) heard of before. From its website, "Journal of Nucleic Acids is a peer-reviewed, open access journal that publishes original research articles as well as review articles in all areas of nucleic acids." By virtue of this structure paper and its citation to 3DNA, I think the journal is surely of personal interest, and I have added it into my watch-list.

To sum up, there are currently four scientific journals (I know of) that are devoted to nucleic acids:

Do I still miss something? Please make your suggestion in the comment area.

[revised on May 17, 2011 by adding RNA Biology]

Ebook "Gregory Petsko in Genome Biology: The first 10 years"

2011-04-23T22:10:00.005-04:00

Over the years, I have read some of Gregory Petsko's monthly columns in Genome Biology while browsing the journal online, and I like his sensible and entertaining columns quite a bit. Recently, I became aware of the ebook from BioMed Central, "Gregory Petsko in Genome Biology: The first 10 years":

Structural biologist Gregory Petsko has contributed a thought-provoking and entertaining monthly column to the scientific journal Genome Biology every month since its launch in 2000. To mark the 10th anniversary of Genome Biology this eBook brings together 10 years of Petsko's columns.

I downloaded the epub version of the book, and googled around, trying to find a corresponding ebook reader for my MacBook Pro (Snow Leopard) – even though I have some ebooks in the generic PDF format, I am not that familiar with epub or mobi. I finally settled with NOOK for Mac from B&N. It turns out reading ebooks with specifically-desinged apps such as NOOK is quite a different, yet more enjoyable, experience than through a PDF reader.

Now the ebook has become the top one in casual reading list. I am reading it from the very beginning, one column at a time, to have a historical perspective. So far I found the columns indeed very "thought-provoking and entertaining".

Tips and tricks from "The Geek Stuff"

2011-04-08T23:46:00.001-04:00

As a devoted command line user, I am always interested in learning new tricks to make my life more enjoyable. Recently, I came across Ramesh Natarajan’s blog “The Geek Stuff” which is full of “instruction guides, how-to, troubleshooting tips and tricks on Linux, database, hardware, security and web” to help solve practical problems.

For example, in the section “Best of the Blog”, I recently benefitted quite a bit by reading the following posts:

There are many other helpful tips/tricks as well; since I have bookmarked the site, I will surely come back!

Scripting in Ruby is fun

2011-04-03T21:11:00.000-04:00

Over the years, I have played around with various scripting languages, including awk, bash, Perl, Python and Ruby. By far, I have enjoyed Ruby the most; nowadays, I write scripts nearly exclusively in Ruby.

Created by Yukihiro "Matz" Matsumoto in Japan during the mid-1990s, Ruby became popular worldwide in mid-2000s, with the Rails web application framework. Indeed, I first dug into Ruby through Rails, and by reading David Black's book "Ruby for Rails; Ruby techniques for Rails developers". As an exercise, I implemented the current 3DNA v2.0 website with Rails v1.x. Then I quickly realized that the rapidly evolving Rails framework was beyond my time and interest to follow. However, I did begin to appreciate Ruby's simplicity, consistency and expressiveness. Over the past few years, I have collected over a dozen Ruby-related (e)books, including "The Well-Grounded Rubyist" (David Black, covering v1.9), "The Ruby Programming Language" (David Flanagan and Yukihiro Matsumoto), and "Metaprogramming Ruby: Program Like the Ruby Pros" (Paolo Perrotta). Just as my experience with (ANSI) C, I feel Ruby "wears well as one's experience with it grows" (K&R, in the preface of "The C Programming Language"). The better I know Ruby, the more I enjoy using it.

I recently wrote two Ruby scripts for the analysis of molecular dynamics (MD) simulation trajectories using 3DNA. Honestly, I would not have bothered with Perl for the task (otherwise, it would have been done long time ago), given the sideline nature of my support of 3DNA. Yet, writing and refining the Ruby scripts (with help of git and rake) have turned out to be a pleasant experience. Another reason why scripting in Ruby is fun is due to its large, active and friendly user community; there are many user-contributed libraries (gems) that serve well of common programming needs. As an example, in the 3DNA-MD scripts, I took advantage of the elegant Trollop commandline option parser by William Morgan. I picked Trollop among many other choices because it is self-contained in a single file, simple to use, and "gets out of your way".

In the Ruby community, exciting new developments are happening all the time. Recently, I was drawn to thor, "a simple and efficient tool for building self-documenting command line utilities". Over the past couple of years, I have browsed Sinatra and Sequel – they also look brilliant! Of course, for bioinformatics, there is the BioRuby project.

Overall, in my experience, scripting in Ruby is fun and exciting. Are you a Rubist yet?

DNA fiber models ABC

2011-03-26T23:42:00.000-04:00

Among the 55 fiber models available in 3DNA, the A-, B- and C-DNA types are the most generic – they can be built with bases A, C, G and T in any combination (see table below). Moreover, in addition to the well-known Arnott fiber models (#1, #4 and #7, all from calf thymus), there are newer variants from van Dam & Levitt (#46 and #47) and Premilat & Albiser (#53 to #55).

 1   32.7   2.548  A-DNA (calf thymus)
 4   36.0   3.375  B-DNA (calf thymus)
 7   38.6   3.310  C-DNA (calf thymus)
46   36.0   3.38   B-DNA (BI-type nucleotides)
47   40.0   3.32   C-DNA (BII-type nucleotides)
53  -38.7   3.29   C-DNA (depreciated)
54   32.73  2.56   A-DNA [cf. #1]
55   36.0   3.39   B-DNA [cf. #4]

As shown in Figure 9 of the 3DNA 2003 NAR paper (linked below), the A-, B- and C-DNA fiber models are all right-handed regular straight helices, yet each has distinguished features.

While I could easily envisioned possible applications of the fiber models, especially in connection with analysis and rebuilding routines in 3DNA, it was still a nice surprise to see a recent article by Gossett and Harvey, titled "Computational Screening and Design of DNA-Linked Molecular Nanowires" [Nano Lett., 2011, 11 (2), pp 604–608]. The abstract is quoted below:

DNA can be used as a structural component in the process of making conductive polymers called nanowires. Accurate molecular models could lead to a better understanding of how to prepare these types of materials. Here we present a computational tool that allows potential DNA-linked polymer designs to be screened and evaluated. The approach involves an iterative procedure that adjusts the positions of DNA-linked monomers in order to obtain reasonable molecular geometry compatible with normal DNA conformations and with the properties of the polymer being formed. This procedure has been used to evaluate designs already reported experimentally, as well as to suggest a new design based on pyrrylene vinylene (PV) monomers.

In the article, 3DNA (the web interface version w3DNA) was cited as follows:

The selection of DNA structures is important because the DNA remains fixed throughout the procedure. To reduce the risk of an incorrect result, one should choose a subset of DNA structures that are in some sense representative of DNA conformational space. The DNA structures (A-, B-, and C-form DNA) were obtained using the Web 3DNA web server. We used a poly(dG)-poly(dC) sequence with ideal geometry for each DNA structure. A-DNA was constructed with rise = 2.548 Å and twist = 32.7˚ , B-DNA was constructed with rise = 3.375 Å and twist = 36.0 ˚, and C-DNA was constructed with rise = 3.310 Å and twist = 38.6 ˚.

Indeed, this is a novel application of fiber DNA ABC models!

3DNA citations reach over 500

2011-03-20T19:26:00.000-04:00

On Friday, June 5, 2009, I blogged on the topic titled "3DNA citations reach over 300". At that time, I wrote (towards the end):

I still remember that the number of citations to 3DNA was less than 150 nearly two years ago [~ summer 2007], 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.

When I checked Google scholar for 3DNA citations right now, the citation number is already over 500 for the initial 2003 3DNA NAR paper alone. Combined with the two direct follow-ups – the 2008 Nature Protocols paper and the 2009 NAR web server paper – the three 3DNA publications have been cited a total of 550 times.

Again, as noted in that blog post,

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.

A decade after its initial humber release, 3DNA has been successfully applied to many real-world problems. As spare time permits, I have actively maintained and continuously refined 3DNA based largely on users' feedbacks. Over the time, I also see clearly that 3DNA can be moved to the next level both in functionality and usability to enjoy an even larger/broader impact.

Now more than half-way through, it won't be long when citations to 3DNA reach 1000, and then beyond.

Review article on NMR analysis of protein–DNA interactions by Milon et al

2011-03-13T23:21:00.001-04:00

Through Google scholar, I became aware of a recent review article by Milon et al., titled "Nuclear magnetic resonance analysis of protein–DNA interactions" in the journal J. R. Soc. Interface:

This review focuses on the experimental strategies currently employed to solve structures of protein–DNA complexes and to analyse their dynamics. It highlights how these approaches can help in understanding detailed molecular mechanisms of target recognition.

I browsed through the text to get myself more familiar with NMR the methodology and its applications in protein-DNA recognition. I was surprised that 3DNA was cited in the article, especially with respect to its unique analyze/rebuild complementarity:

In addition, several software programs have been developed to model DNA bending such as the 3DNA program, which allows analysis of DNA structural parameters and enables it to be rebuilt with customized DNA models [76]. Several Web servers have been created recently and provide interesting tools to analyse and rebuild DNA models [77,78].

I am only wishing that 3DNA's neat features could be more widely recognized; hopefully I'd have the opportunity to further refine 3DNA and move it to the next level.

Retraction of scientific publications

2011-03-05T23:24:00.000-05:00

Once in a while, I come across retraction notices of scientific publications in leading journals/magazines. Even for cases not directly related to my research areas, I normally browse through them.

In the March 3, 2011 issue of Nature, there is a retraction of the Letter "Mediation of pathogen resistance by exudation of antimicrobials from roots" [Nature 434, 217–221 (2005)]. I am intrigued by the first sentence of the note:

The authors wish to retract this Letter after a key reference by Walker et al. (ref. 9 in this Letter) was retracted from the scientific literature.

It turns out that the 2003 Walter et al. J. Agric. Food Chem. paper (withdrawn in October 2009) and the 2005 Nature Letter were from the same group. Overall, it took ~6 years each for the two papers to be retracted. As of today, they have been cited 76 and 84 times respectively accordingly to Google scholar.

Evidences for transient Hoogsteen base pairs in canonical DNA duplex

2011-02-27T18:23:00.000-05:00

In the February 24, 2011 issue of Nature, there is an interesting article by Nikolova et al., titled "Transient Hoogsteen base pairs in canonical duplex DNA". Its main discovery is succinctly summarized in the abstract:

By using nuclear magnetic resonance relaxation dispersion spectroscopy in concert with steered molecular dynamics simulations, we have observed transient sequence-specific excursions away from Watson–Crick base-pairing at CA and TA steps inside canonical duplex DNA towards low-populated and short-lived A•T and G•C Hoogsteen base pairs. The observation of Hoogsteen base pairs in DNA duplexes specifically bound to transcription factors and in damaged DNA sites implies that the DNA double helix intrinsically codes for excited state Hoogsteen base pairs as a means of expanding its structural complexity beyond that which can be achieved based on Watson–Crick base-pairing.

Geometrically, the Hoogsteen base pair is related to the Watson-Crick base pair by a 180-degree rotation about the glycosidic bond (N9–C1'). While the A•T Hoogsteen base pair is classic, the similar G•C+ Hoogsteen pair (with protonation of cytosine N3) is equally possible. The A•T and G•C Hoogsteen base pairs have two perfect H-bonds, so they are energetically stable. As for their existence in DNA duplex, the most direct evidence comes from the "trap" experiments (see Fig.3 of the paper). In the News & Views section, Honig and Rohs provide a nice recap of the main point and implications of this work.

As also observed in another recent publication, "Replication infidelity via a mismatch with Watson–Crick geometry", the base sequence has a subtle role in influencing the base-pairing schemes, three-dimensional structures and biological functions of DNA. However, we should not forget that only the Watson-Crick base pairs, and to a less extent, the G-U wobble pair, have the correct symmetry to ensure a "regular" double helical structure.

Canned responses in gmail make it easy to send common messages

2011-02-20T23:35:00.002-05:00

Through Gary Rosenzweig's MacMost Now video #509 "Gmail Labs" (January 28, 2011), I first heard of "Canned Responses" in Gmail Labs:

Email for the truly lazy. Save and then send your common messages using a button next to the compose form.

This is a truly handy feature that I have long been waiting for! Yet even though I am aware of Gmail Labs and enabled quite a few experimental features a while ago, I've not been searching Gmail Labs for new features ever since.

Over the past few weeks, I have found "Canned Responses" increasingly indispensable in my support of the 3DNA forum (as a sideline project):

When I activate a new 3DNA forum registration, I've always included a "standard" message to "make the forum policy upfront and explicit, in order to avoid misunderstandings or surprises." Previously, I had to copy-and-paste, e.g., from a specifically created text file or elsewhere. Surely, this worked, but I had felt intuitively that there must be a better way to get the job done. Well, that's exactly where "Canned Responses" fit in!
Over the past few months, I have been ever more bothered by spam registrations. So as a further filter, I have been sending the following enquiry message to each suspicious registration for activation:
Thanks for your registration at the 3DNA forum. Please tell me a little bit about yourself and elaborate on how 3DNA could be useful to your project; we would like to make the forum spam-free.

See also "Further notes on forum registration and posting" -- you may not need to register.
Here once again, the "Canned Responses" feature makes my life much easier! Moreover, this step turns out to be extremely effective; a large percentage of registrations is filtered out at this final stage.

Are you using gmail? If so, you may also want to give "Canned Responses" a try.

Making data maximally available?

2011-02-13T23:00:00.000-05:00

In the February 11, 2011 issue of Science (Vol. 331 no. 6018 p. 649), there is an editorial, titled "Making Data Maximally Available". Indeed, the issue contains a special section on "Dealing with Data".

Science is driven by data. New technologies have vastly increased the ease of data collection and consequently the amount of data collected, while also enabling data to be independently mined and reanalyzed by others ... It is obvious that making data widely available is an essential element of scientific research.

Especially, I like the following two (proposed) new policies:

To extended data access requirement "to include computer codes involved in the creation or analysis of data." If properly implemented/enforced, this policy could significantly increase the repeatability and assessment of published results. In my experience, I have observed too many times that secrets are hidden in the seemingly "little" subtle details.
"To produce a single list that combines references from the main paper and the SOM" (supporting online material) to "provide credit and reveal data sources more clearly". Potentially, this will also increase the citation of method papers.

Hopefully, other journals will follow Science's lead to make data maximally available, and to present data more transparently.

A G-T mismatch with perfect Watson-Crick geometry

2011-02-06T23:24:00.001-05:00

In the February 1, 2011 issue of PNAS [108(5)], there is an interesting article "Replication infidelity via a mismatch with Watson–Crick geometry" by Bebenek et al. They solved the Pol λ DL ternary complex (PDB id: 3PML) which has a G-T nascent mispair in "perfect" Watson-Crick geometry (see their Fig 4, linked below).

From an H-bonding (energitic) point of view, the G-T mispair (with three "H-bonds") can only be possible if G or T is in the rare enol tautomeric state, or is ionized. The pH dependence of single nucleotide disincorporation seems to be consistent with an ionized base pair. Note the G-T mispair is different from a Wobble pair in which G and T have a relative sheared motion (see also Fig. 4C above).

I am glad to find that 3DNA was used in deriving the parameters. By design, 3DNA should be able to identify such "unusual" mispair as easily as for a normal Watson-Crick pair. As noted in our 2008 3DNA Nature Protocols paper,

By taking advantage of the standard base reference frame and selected geometric features, the find_pair program within 3DNA can identify all possible nucleic-acid base pairs, whether they are canonical Watson–Crick or noncanonical pairs and are made up of normal or modified bases, in any tautomeric form or protonation state. (p1217)

Moreover, 3DNA does notice and signify (wit a * instead of the normal -) the atypical H-bonding feature of the G-T mispair to draw further attentions.

5 T-*---g  [3]  O2 - N2  3.06  N3 * N1  2.97  O4 * O6  2.66

Hopefully, this example helps illustrate some of 3DNA's unique features that would hopefully be more widely recognized and applied.

Trial of scientific papers by blogs and tweets?

2011-01-30T23:50:00.000-05:00

In the January 20, 2011 issue of Nature, there is an interesting News Feature article by Apoorva Mandavilli, titled "Peer review: Trial by Twitter":

Blogs and tweets are ripping papers apart within days of publication, leaving researchers unsure how to react.

Specifically, two widely publicized papers in Science are singled out: one is about the longevity genes identified through genome-wide association study (GWAS), published last July; another is a more recent one on arsenic bacteria, published last December. In each case, while "the popular media was trumpeting the finding, other researchers were taking to the web to criticize the paper’s methodology." Yet, the authors failed to hold up their claims in the papers.

With great interest, I've been following the story on arsenic bacteria. I first noticed this work through Science Podcast, and found the topic of "arsenic" life intriguing. So for general knowledge, I read carefully the abstract and browses through the text. One week later came the Nature editorial "Response required", and from which, I followed the link to Rosie Redfield's blog post "Arsenic-associated bacteria (NASA's claims)". I read the post, and many of the comments therein; while I do not understand many of the technical details, I had no difficulty in following her argument. The lead author of the arsenic bacteria paper, Dr. Felisa Wolfe-Simon, did respond to comments on December 16, 2010. Science also published an interview with Wolfe-Simon, titled "Discoverer Asks for Time, Patience Over Arsenic Bacteria Controversy". On the same day, Redfield was quick to write another blog post "Comments on Dr. Wolfe-Simon's Response", which again has received many comments. So far, the story is still unfolding, and Science has promised to publish technical comments and responses in early 2011.

As a related topic, throughout the CCP4bb, I noticed the letter to editor titled "Is too ‘creative’ language acceptable in crystallography?" by Alexander Wlodawer et al. I agree fully with the authors that "While figures of speech are often useful and even educational, flashy titles combined with hyperbolae and imprecise language can mislead or deceive nonspecialist readers and should therefore be avoided."

In the Internet Age, bloggers and tweeters clearly have an important role to play in the assessment of research findings. In scientific publications, what counts is not how much one claims, but to what extent one can hold up such claims. Solid work holds up over time and the scrutiny of peers.

Three structural biology papers in the latest issue of NAR cite 3DNA

2011-01-22T20:14:00.002-05:00

While browsing the latest 39(2) January 2011 issue of Nucleic Acids Research (NAR), I found, to my great surprise, three papers that cite 3DNA. These papers, all under the "structural biology" section, are of interest to me from their titles and abstracts, so I downloaded the PDF versions and read through each of them.

For this blog post, #100 by incidence, it would be intriguing to look into the context to see how 3DNA is cited.

"Asymmetric DNA recognition by the OkrAI endonuclease, an isoschizomer of BamHI" by Vanamee et al. (Mount Sinai School of Medicine, and New England Biolabs):

Analysis of the stereochemical quality of the protein model and assignment of secondary structure were conducted with PROCHECK (13). DNA analysis was performed with 3DNA (14). Solvent-accessible surface areas were calculated in CNS with the algorithm of Lee and Richards employing a 1.4-Å probe(15). Figures were prepared using PyMOL (www.pymol.org). [p713, from bottom left to middle right]

"DNA intercalation without flipping in the specific ThaI–DNA complex" by Firczuk et al. (Poland, Germany and UK):

An oligoduplex with the correct sequence in standard B-DNA geometry was generated with the program 3DNA (44), and manually adjusted to fit the highly distorted DNA in the structure. ... The programs COOT (45), REFMAC (46) and CNS (47) were used for refinement. [p747, top left]

Analysis with the 3DNA software (44) shows that the intercalation increases the rise between base pairs to about 7 Å or approximately twice its usual value (Figure 5B). Phosphorus–phosphorus (Pn–Pn+1) distances in the DNA backbone are only mildly altered (values range from 5.6 to 7.0 Å). Instead, the extra height of the two CG steps comes at the expense of the twist, which is reduced from its usual value of about 36° (360°/10) to between 10 and 15°. A view toward the major groove shows that the inner base pairs of the recognition sequence are strongly tilted (Figure 5). According to the 3DNA software (44), the first CG step has a negative tilt of about ~12°, which results in the oblique orientation of the following base pairs. The central GC step is characterized by a tilt close to 0°, reflecting the nearly parallel arrangement of the middle bases. Finally, the second CG step has a positive tilt of about 15° which restores the standard orientation of the downstream base pairs. A side view of the DNA indicates a bend at the center of the recognition sequence which is primarily due to the positive ~12° roll of the central GC step into the major groove (Table 1). The 3DNA program also indicates that the propeller twist is positive for the specifically recognized sequence, and (as expected for the standard B-DNA) negative for most of the flanking base pairs. [p749, top right]

Table 1. DNA distortion in complex with ThaI restriction endonuclease: all parameters were calculated with the 3DNA software (44). [p750, middle left]

"On the molecular basis of uracil recognition in DNA: comparative study of T-A versus U-A structure, dynamics and open base pair kinetics" by Fadda and Pomès (Ireland and Canada):

MD simulations were run with versions 3.3.3 up to 4.0.4 of the GROMACS software package (47,48).

Structural parameters were determined with the 3DNA software package (51,52). The pymol (www .pymol.org) software package was used to generate figures. [p769, bottom right]

Established in 1974 and currently with an impact factor of 7.479, NAR has also been chosen by the Special Libraries Association as one of the top 100 most influential journals in medicine and biology over the last 100 years. The citations by the three papers in the latest issue of NAR illustrate unambiguously 3DNA's big impact in structural biology.

Ruby scripts for 3DNA analysis of molecular dynamics simulation trajectories

2011-01-19T23:44:00.001-05:00

Over the years, I've been very pleased to see 3DNA's ever-increasing applications for the analysis of molecular dynamics (MD) simulation trajectories of nucleic acid structures. Among its other features, this illustrates that the command-line driven approach of 3DNA makes it easily integrable into the MD analysis pipeline (with some scripting, of course).

However, the lack of direct support of 3DNA to the ever more popular field of MD simulations has caused several obvious problems:

Repeated efforts – virtually every lab or even MD practitioner could come up with an ad hoc scripting solution.
Hinderance to 3DNA's even wider adoption – new comer to the MD field, or bench scientists interested in dynamics simulations would be scared off.
Known issues with existing approaches – most predominately the unnecessary repetitive run of find_pair to deduce base pairing information for each snapshot (model), which not only takes time, but more seriously some pairs could be missing due to melted out or distortion along the trajectory.

I've been following 3DNA's citations for years and I am well aware of the above issues: in addition to answering relevant questions in the 3DNA forum, I have blogged specifically on the topic a few times:

3DNA for the analysis of molecular dynamics simulations [Saturday, July 24, 2010]
3DNA in the PCCP nucleic acid simulations themed issue [Sunday, December 6, 2009]
3DNA in molecular dynamics simulations [Sunday, October 4, 2009]

Of course, I am in a unique position to help solve the problem. Indeed, for the past couple of years, I've been thinking of writing scripts to make life easier for MD practitioners who care to use 3DNA. However, due to my lack of experience in MD simulations, constraints of "spare" time (plus laziness), and a want of suitable collaborator, I've never found the incentive to get the job done.

I finally decided to write some Ruby scripts to streamline the process of using 3DNA in MD simulations, after a recent question from Aneesh on "script for extracting data from 3DNA output file" in the 3DNA forum. After a few exchanges of views with Aneesh, and especially with Alpay's contribution of Python script and sample dataset, I've finished up two standalone yet connected Ruby scripts to analyze MD simulation trajectories with 3DNA and then extract various structural parameters. The details, including source code and test examples, are available in the 3DNA forum under "Ruby scripts for the analysis of MD simulation trajectories", in a newly created section titled "Molecular dynamics simulations".

The sample file ("sample_md0.pdb") distributed with the current v0.1 of the scripts contains 21 snapshots (models, 0..20), separated by MODEL/ENDMDL pairs. While the sample is based on a trajectory file from AMBER, any MD simulation packages, or NMR ensembles, can be similarly handled as well.

Now the ball is rolling. As time goes by, and with users' feedback, I will refine and expand the functionality of the scripts as necessary. I am confident to see more applications of 3DNA in the "dynamic" molecular simulation field.

Open dictionary with command-control-d in Mac OS X

2011-01-09T22:46:00.000-05:00

In a recent MacMost Newsletter, I came across the following handy trick in Mac OS X: while the cursor is over a word (not necessarily selected), one can press command-control-d to open dictionary which pops up a little window with the word's definition. This works in Safari and Mail, but not in Preview (unfortunately).

Previously, when I need to check the definition of a word, I right-click on it and then follow the link "Look Up in Dictionary". This will launch or pop up Dictionary with detailed information about the word (Dictionary/Thesaurus/wikipedia).

The right-click method seems to integrate better with other Mac OS X applications. For example, it works with Preview as well. However, now that I know it, I sometimes prefer the command-control-d approach better; it is quick and non-obstructive.

IUPAC nucleotide symbols and their complements

2011-01-09T22:07:00.000-05:00

Recently, I was interested in knowing the complements of all the IUPAC nucleotide symbols. As blogged previously, I am quite familiar with the "meaning of nucleotide IUPAC codes" (namely A/C/G/T, and R/Y/N etc). However, when I first check the Gene Infinity website on nucleotide symbols, it still puzzled me for awhile to figure out the meaning of the DNA alphabet (with complements), as except below:

A  C  G  T    M  R  W  S  Y  K    B  D  H  V    N
|  |  |  |    |  |  |  |  |  |    |  |  |  |    |
T  G  C  A    K  Y  W  S  R  M    V  H  D  B    N

For example, some degenerated IUPAC symbols are complemented to themselves (e.g., W–W and S–S), while others are seemingly "hard" to apprehend (e.g., B–V and D–H).

After thinking it for a bit, things begin to become clear. They are based on the complementarity of Watson-Crick base-pairs (A–T and G–C) and the meaning of each degenerated IUPAC nucleotide symbol. For example,

W represents A/T, meaning weak (with only two hydrogen-bonds). The complements of A/T are T/A respectively, which is W again.
B (not A) represents C/G/T, and their complements are G/C/A respectively, which is V (not U/T).

It is easy to verify that all other complementary pairs follow exactly the same basic principle.

Extract images from PDF files using 'pdfimages'

2010-12-13T19:52:00.000-05:00

Once in a while, I need to extract an image (or a portion thereof) from a PDF file. Usually, I open the PDF file using 'preview' or Adobe 'Acrobat Reader', and take a screenshot which is then cropped to the desired slice. This manual method "works", albeit a bit tedious.

Recently, I came across the handy command-line utility program 'pdfimages' which allows for automatic extraction of all images from a PDF file. The basic usage is simply:

pdfimages   yourfile.pdf   prefix

This will extract all images contained in yourfile.pdf to files prefix-001.ppm, prefix-002.ppm etc. With option "-j", images in DCT format are saved as JPEG files.

Nice blog on English writing by Lynn Gaertner-Johnston

2010-12-05T23:53:00.000-05:00

A while back, I was somehow curious to figure out the differences between entitled vs. titled. A colleague referred me to a blog post on this topic by Lynn Gaertner-Johnston. After reading the post, it became clear that over the years I had made the same mistake numerous times. I then corrected nearly all occurrences of "entitled" to "titled" in my blog and 3DNA forum posts. Ever since, I have been following Lynn Gaertner-Johnston's blog, titled "Business Writing", on a regular basis, and found it helpful to improve my English writing skills.

Under-citation of method papers?

2010-11-28T22:23:00.001-05:00

Recently in the CCP4BB, there is an interesting thread with extensive discussions on "Citations in supplementary material". The original poster refers to the recent Acta Cryst D editorial with the same title in which the authors highlights the issue of under-citation to papers published in the International Union of Crystallography (IUCr) journals.

The main point is that method papers are more likely to be cited in the supplementary materials only, which are not indexed by PubMed, Scopus, Web of Science or Google Scholar etc. As a result, they are statistically undercounted, and "Journals and scientists that focus on publishing methodologically oriented papers are particularly affected." Specifically, through a survey of articles on protein or nucleic acid structure determination published in Cell, Nature, Science, and PNAS in 2009, the authors found that "almost half of all references to publications in IUCr journals end up being published in the supplementary material only."

The findings of the editorial resonate with my observations, and I cannot agree more with the authors that "in the end, methods need to be continuously developed and refined in order to ensure progress."

On the other end of the spectrum, some highly influential method papers are heavily cited. As an extreme case, the large number of citations to the 2008 paper "A short history of SHELX" by George Sheldrick helps rocket up the impact factor of Acta Crystallographica A by 20-fold to 49.9 this year!

Belorussian translation of 3DNA webpages

2010-11-21T22:41:00.000-05:00

Recently, I communicated with Paul Bukhovko on the translation of 3DNA webpages http://rutchem.rutgers.edu/~xiangjun/3DNA/ into Belorussian. As the author of the original website, referred to as http://rutchem.rutgers.edu/~olson/3DNA/ (which is simply a soft link to the above URL) in the 2003 3DNA NAR paper, I was very surprisingly pleased when Paul asked for permission to perform the translation, which I gladly granted.

Regarding the process, Paul commented:

Was a pleasure to translate this page! It's kinda fresh and related to my professional interests, so I thought - why not, if the author allows to do so.

The translated page is at URL: http://www.movavi.com/opensource/3DNA-be. Interestingly, when I used Google Translate to convert the Belorussian version back to English, the outcome is pretty readable. In contrast, when the original English version is directly translated to Chinese, the result is beyond recognition!

Proper labeling of O1P and O2P atoms in a phosphate group

2010-11-14T23:26:00.003-05:00

Recently, a question on the 3DNA o1p_o2p utility program in the forum led me to reflect on the proper labeling of O1P/O2P atoms in a phosphate group. As is well-known, in DNA/RNA structures, the phosphate group (see figure below left) connected two neighboring nucleosides.

The two nonbridging oxygen atoms of the phosphate group (the horizontal Os in the O-P=O line, left) are named O1P and O2P in PDB files (see also figure to the right). Stereochemically, O1P and O2P are also designated as pro-R and pro-S oxygens, respectively.

Presumably, the structural files in the PDB and NDB databases should be consistent and follow the standard nomenclature. In practice, however, some entries in the NDB had mislabeled O1P/O2P atoms (e.g., adh026). I first noticed this issue when I superposed the A-DNA adh026 to its 3DNA rebuilt version with the sugar-phosphate backbone. I observed an unreasonably large RMSD only for the octamer adh026, while the RMSDs were much smaller (as expected) for the B-DNA dodecamer bdl084 and the 146-bp nucleosomal DNA in pd0001 (see $X3DNA/examples/analyze_rebuild distributed with 3DNA v2.0). Since the standard building blocks in 3DNA were applied consistently, I traced the cause of the large RMSD problem to the PDB file of adh026 itself, and finally identified it was actually due to the mislabeling of the O1P/O2P atoms.

The utility program o1p_o2p was written specifically for the purpose of checking if the O1P/O2P atoms are properly labeled in a PDB file. In a phosphate group, if O1P/O2P are correctly labeled, then following O1P-->O2P-->O5' in a right-handed sense would point the thumb in the direction of O3' (see the figure up right). As always, how it actually works is best illustrated with an example. Shown below, the second phosphate in adh026 is used (as distributed with 3DNA), with GNU octave script. Here the O1P/O2P atoms are mislabeled since direction has a negative value. In contrast, for a properly labeled phosphate group, direction should be positive.

#ATOM      6  O3*   G A   1       8.396  -3.995  -1.948  1.00 30.86           O  
#ATOM     20  P     G A   2       8.163  -3.069  -0.619  1.00 32.38           P  
#ATOM     21  O1P   G A   2       7.401  -1.917  -1.218  1.00 32.09           O  
#ATOM     22  O2P   G A   2       7.280  -3.934   0.195  1.00 34.05           O  
#ATOM     23  O5*   G A   2       9.600  -2.800  -0.121  1.00 29.41           O  

P   = [8.163  -3.069  -0.619]
O1P = [7.401  -1.917  -1.218]
O2P = [7.280  -3.934   0.195]
O3  = [8.396  -3.995  -1.948]
O5  = [9.600  -2.800  -0.121]

O1P_to_O2P = O2P - O1P                    # -0.12100  -2.01700   1.41300
O2P_to_O5 = O5 - O2P                      #  2.32000   1.13400  -0.31600
O1P_O2P_O5 = cross(O1P_to_O2P, O2P_to_O5) # -0.96497   3.23992   4.54223
P_to_O3 = O3 - P                          #  0.23300  -0.92600  -1.32900

direction = dot(O1P_O2P_O5, P_to_O3)    # -9.2616 < 0: O1P/O2P mislabeld

The O1P/O2P labeling issue is just a little detail I came cross while developing 3DNA. Nevertheless, it serves as an excellent example of the subtleties ~~subtitles~~ that should be taken care of in scientific programming.

Please note that as of 2008, in the remediated PDB/NDB entry adh026, the mislabeled O1P/O2P pair has been correct. More generally, O1P/O2P atoms have now been renamed as OP1/OP2, respectively. 3DNA v2.0 takes care of such naming changes internally; for generated PDB files, however, 3DNA still adopts the conventional O1P/O2P labeling.

Transparency in the peer-review process of scientific papers

2010-11-06T23:44:00.000-04:00

In the Nov. 4, 2010 issue of Nature, there is an interesting Comment, titled "Transparency showcases strength of peer review", by Bernd Pulverer, head of scientific publications at the European Molecular Biology Organization and chief editor of The EMBO Journal. In this article, Pulverer "reflects on his experience at The EMBO Journal of publishing referees’ reports, authors’ responses and editors’ comments alongside papers."

The peer-review process of scientific articles has traditionally been a "black box": (anonymous) reviewers' reports, editors' comments, and authors' responses – extremely valuable information in shaping the final form of published papers – are all hidden from public view. In the Internet era, technology (e.g., online space) is no longer an issue. Now The EMBO Journal has led the way, and "the experience has been overwhelmingly positive." Hopefully, other leading journals (e.g., Nature and Science) would follow the example. Afterall, making the peer-review process transparent is an excellent mean to increase the accountability of science and scientific publications.

Overall, this article is well-written, succinct and logical, and it touches an important topic in scientific publication. Over the past few days, I have read several "Review Process Files" accompanying papers I am interested in, e.g., "Recognition of the amber UAG stop codon by release factor RF1", and found them highly revealing.

Publication of scientific programming code

2010-10-30T22:42:00.001-04:00

Recently on the Nature website, I read with great interest a news article, titled "Publish your computer code: it is good enough", by Nick Barnes, a professional software engineer:

Freely provided working code — whatever its quality — improves programming and enables others to engage with your research

Clearly the author knows the "trade secret" in scientific programming. He lists several common reasons why scientists are reluctant to share their source code, and then provides his responses:

The code is low quality — "software in all trades is written to be good enough for the job intended". All software has bugs. Sharing code would help improve the code itself and advance the research field.
Not a common practice — this is going to change or is already changing.
Demand for support — "Nobody is entitled to demand technical support for freely provided code."
Intellectual property issue — The most value part "lies in your expertise", code not backed by skilled experts is called abandonware. (I cannot agree more with this point.)
Polishing code takes time/effort — not need to, just supply, as supplementary materials in a website, the original code used in your publication.

As is evident from the many comments, this assay is well echoed by the community. As an active computational scientist for over a decade, I share mostly the same opinions. Essentially, the transparency of source code is to ensure repeatability of scientific publications. In the field of computational biology (bioinformatics), it is virtually impossible to reproduce exactly a published figure/table without direct access to details, including the source code.

Chi (χ) torsion angle characterizes base/sugar relative orientation

2010-10-23T00:04:00.004-04:00

Except for pseudouridine, a nucleoside in DNA/RNA contains an N-glycosidic bond that connects the base to the sugar. The chi (χ) torsion angle, which characterizes the relative base/sugar orientation, is defined by O4'-C1'-N1-C2 for pyrimidines (C, T and U), and O4'-C1'-N9-C4 for purines (A and G).

Normally (as in A- and B-form DNA/RNA duplex), χ falls into the ranges of +90° to +180°; –90° to –180° (or 90° to 270°), corresponding to the anti conformation (Figure below, top). Occasionally, χ has values in the range of –90° to +90°, referring to the syn conformation (Figure below, bottom). Note that in left-handed Z-DNA with CG repeating sequence, the purine G is in syn conformation whilst the pyrimidine C is anti.

Presumably, the χ-related anti/syn conformation is a very basic and simple concept. In essence, though, the N-glycosidic bond and the corresponding χ torsion angle illustrate that the base and sugar are two separate entities, i.e. there is an internal degree of freedom between them. In this respect, it is worth noting that the Leontis-Westhod sugar edge for base-pair classification corresponds to the anti form only. When a base is flipped over into the syn conformation, the "sugar edge", defined in connection with the minor (shallow) groove side of a nitrogenous bases, simply does not exist.

Base-flipping (anti/syn conformation switch) is one of the factors associated with the two possible relative orientations in a base pair, characterized explicitly in 3DNA as of type M+N or M–N since the 2003 NAR paper (Figure 2, linked below). I reemphasized this distinction in our 2010 GpU dinucleotide platform paper (in particular, see supplementary Figure S2). Unfortunately, this subtle (but crucial, in my opinion) point has never been taken seriously (or at all) by the RNA community, even with 3DNA's wide adoption. However, as people know 3DNA deeper/better and take RNA base-pair classification more rigorously, I have no doubt they will begin to appreciate the simplicity of this explicit distinction and the resultant full quantification of each and every possible base pair using standard geometric parameters.

On a related issue, current versions of 3DNA (v1.5 and v2.0) output only the χ torsion angle without providing the anti/syn classification. This defect, and many others, will hopefully be rectified in future releases of 3DNA.

Improving the design of existing code by refactoring

2010-10-15T23:08:00.000-04:00

Another software engineering book I read recently is "Refactoring: Improving the Design of Existing Code" by Martin Fowler. According to the author,

Refactoring is the process of changing a software system in such a way that it does not alter the external behavior of the code yet improve its internal structure. It is a disciplined way to clean up code that minimize the chances of introducing bugs. In essence when you refractor you are improving the design of the code after it has been written.

The book is practical in nature; it not just explains the principles but provides a detailed account of over 70 commonly used refactorings. As vividly explained by the author in the first paragraph of Chapter 1, "Refactoring, a First Example", "it is with examples that I [the author] can see what is going on." This approach fits my style perfectly: to really understand a topic, I always find a worked example far more effective than general principles. While the examples in the book are illustrated in Java, the basic ideas can be applied well to other object-oriented or even procedural languages (such as C).

Over the years since I left Dr. Olson's laboratory at Rutgers, I have been maintaining and continuously refining 3DNA. I have taken each user's question as an opportunity to fix bugs and improve its design, thus making the code more robust and efficient. The majority of my efforts, as I now realize, is "refactoring" existing 3DNA code to make it easier to maintain and extend. Reading through this book gives me the chance to put my practices into a broader context. I will surely take advantage of some refactoring examples from the book for further refinements of 3DNA.

NSMB editorial: "Go figure"

2010-10-13T22:08:00.004-04:00

In the October 2010 issue of Nature Structural & Molecular Biology (NSMB, Vol. 17, No. 10) there is another interesting one-page editorial, titled "Go figure", which provides tips on how to make a scientific figure that may worth 1000 words:

A picture may be worth a thousand words, but ensuring that those words make sense is important, especially in the context of a scientific figure. Here are some tips for making your figures count.

A recap of the tips is given below; by and large, they all follow conventional wisdom:

General considerations: Each figure should make just one point and be self-explanatory.
See guidelines. "At all stages, the figures should be clear and legible."
How many figures? The figures should complement the Results section, and be included only necessary.
How many panels? Better only one; multiple panels "should be logically connected."
What’s in a label? Keep it succinct, but make the figure self-explanatory.
Getting colorful. Use color wisely and constantly.
A legendary figure. The figure legend should concise and informative.
A model paper. Better have a figure (at the end) of the final model that conveys "the big picture". Honestly, I do not quite get this point.

As pointed out by the author, "These are just a few guidelines and suggestions for handling figures." Overall, "simplicity rules in scientific figures, as in life." I guess no one would argue with such general advices. However, it would be even more helpful to illustrate such points with concrete examples (I know that seems to be beyond the scope of a one-page editorial).

Further details of the DNA story revealed by Crick's lost correspondence

2010-09-30T23:42:00.001-04:00

In the September 30 issue of Nature (Vol. 467, pp519-524), there is an interesting account of "The lost correspondence of Francis Crick" by Gann & Witkowski. The newly found letters, mostly between Crick and Wilkins, unveil further background information on the exciting DNA story. As the authors put it, "Strained relationships and vivid personalities leap off the pages."

I read Watson's "The Double Helix" book a while ago, and overall I am quite familiar with the DNA story. Still, I found this account fascinating: it provides a "CAST LIST" in "The search for the structure of DNA" with photos (p521); and it succinctly summarizes the relationships among the key players. In science, no other story shows more dramatically the collaborative and competitive nature among scientists working on similar projects.

Franklin’s X-ray diffraction photograph 51 of B-form DNA (see figure above, from Wikipedia), with its unambiguous evidence that DNA was helical, proved crucial for Watson and Crick to determine the structure of DNA. Indeed, the Watson-Crick DNA model corresponds to the B-form DNA, with its base-pairs in the middle, parallel to each other and perpendicular to the linear helical axis.

From this Nature account, however, I noticed for the first time a subtle detail: when B-form DNA photograph 51 was shown to Watson by Wilkins in early 1953, Franklin also already had the A-form DNA diffraction pattern. According to the authors,

It was the A-structure diffraction pattern that had led Franklin away from believing that DNA, in that form at least, was helical, despite her already having produced the most persuasive helical pictures of the B structure — including photograph 51. The crystalline DNA gave better quality diffraction data, more suited to her painstaking, quantita- tive approach, and so she focused on the A form during 1952. It was at this time that she and Gosling made a handwritten, black edged funeral card announcing the death of “DNA Helix(crystalline)”.

When Crick had the opportunity to look the A-form DNA diffraction picture, on 5 June 1953, he wrote (to Wilkins):

This is the first time I have had an opportunity for a detailed study of the picture of Structure A, and I must say I am glad I didn’t see it earlier, as it would have worried me considerably.

I am reading "Blink: The Power of Thinking Without Thinking", a book by Malcolm Gladwell. The above case serves as a vivid example.

New NDB entries have IDs start with the prefix NA?

2010-09-19T22:22:00.001-04:00

Recently, while reading the article titled "Designing Triple Helical Fragments: The Crystal Structure of the Undecamer d(TGGCCTTAAGG) Mimicking T.AT Base Triplets" by Van Hecke (Crystal Growth & Design), I noticed the following:

The atomic coordinates and structure factors have been deposited in the Protein Data Bank⁴¹ and Nucleic Acid Database⁴² (PDB and NDB entry codes 3L1Q and NA0392, respectively).

The NDB id NA0392 reminded me of an email communication I had with Dr. Olson early this year when she told me of the id change of new NDB entries. Occupied with other issues, I did not pay much attention to this point until recently.

Over the years, NDB has established itself prominently as "a repository of three-dimensional structural information about nucleic acids". Traditionally, an NDB id has some associated "meanings", with noticeably exceptions, as discussed one year ago in my blog post "PDB id vs NDB id". Thus, it is not surprising that NDB decided to make changes in naming new ids. However, I cannot find any announcement on the id change in the NDB website; a Google search on "NDB id" did not uncover anything new either. Luckily, NDB provides search by release date ("Released Since"). After a few tries, I traced that the new id policy began to be implemented from around March 2010.

The new NDB id convention appears to be NA (presumably standing for for Nucleic Acid) followed by 4 digits. As more entries available, it is not that hard to imagine that the number of digits must be expanded to 5, 6, or even more. Naturally, PDB id, with a fixed 4-character length (up to now), is (far) more consistent than the NDB id is.

"Code Complete", a practical handbook of software construction

2010-09-12T23:32:00.002-04:00

Over the summer, I read quite a few books on C programming, and more generally on software construction. Among them, the book titled "Code Complete" (2nd edition) by Steve McConnell is the most comprehensive: with over 900 pages, it certainly serves as "a practical handbook of software construction".

I have been writting scientific software applications for over twenty years, using a variety of programming languages. Gradually, writing code per se is taking less time; nowadays by far the most time-consuming parts have increasingly become (1) to understand a scientific topic thoroughly in order to implement it in new code, and (2) to maintain/refine/adapt existing codebase to meet the needs of its user community (e.g., 3DNA). It is for the benefit of the later that I've found the "Code Complete" book very helpful.

As an example, following the book's recommendations (1) "Use boolean variables to simplify complicated tests" (pp301-2) and (2) "Simplify complicated tests with boolean function calls" (p359), I've easily improved the readability of some parts of my code.

Overall, the book is well-written and full of practical advices; I enjoyed reading it.

Identification of C-H...N/O H-honds using 3DNA

2010-09-05T22:20:00.009-04:00

Recently, I came across an interesting article by Kiliszek et al., titled "Atomic resolution structure of CAG RNA repeats: structural insights and implications for the trinucleotide repeat expansion diseases". In addition to its biological implications, this paper uses 3DNA to deal with non-canonical-base-pair-containing helical structures in a sensible way:

The helical parameters were calculated using 3DNA (29). Sequence-independent measures were used, based on vectors connecting the C1' atoms of the paired residues, to avoid computational artefacts arising from non-canonical base pairing.

Another significant point, which is the focus of this post, is the observation that "All the adenosines are in the anti-conformation and the only interaction within each A-A pair is a single C2-H2...N1 hydrogen bond." (Figure below) Given the 0.95 Å ultrahigh resolution of structure 3nj6/na0608, it is likely that this type of A-A is real.

The find_pair program in the currently distributed versions of 3DNA (v2.0 and before), however, does not identify this A-A pair (thus the corresponding NDB list of "Base Pair Step Parameters" is incomplete) for the following two reasons:

No H-bond exists between N/O base atoms – currently a requirement for a base-pair.
By default, only N/O atoms are used in defining H-bonds (see tag hb_atoms in file "misc_3dna.par"). Nevertheless, by adding C as a possible atom in forming H-bond, and manually editing find_pair generated input file to analyze, 3DNA structural parameters can be calculated as usual.

I have updated find_pair in 3DNA to identify such C-H...N/O H-bond automatically (see below for entry 3nj6). Upon further refinements and validations, future releases of 3DNA will have this functionality available.

    1   95  #    1 | ...1>A:...1_:[..G]G-----C[..C]:..10_:A<...2
    2   94  #    2 | ...1>A:...2_:[..G]G-----C[..C]:...9_:A<...2
    3   93  #    3 | ...1>A:...3_:[..C]C-----G[..G]:...8_:A<...2
    4   92  #    4 | ...1>A:...4_:[..A]A-**--A[..A]:...7_:A<...2
    5   91  #    5 | ...1>A:...5_:[..G]G-----C[..C]:...6_:A<...2
    6   90  #    6 | ...1>A:...6_:[..C]C-----G[..G]:...5_:A<...2
    7   89  #    7 | ...1>A:...7_:[..A]A-**--A[..A]:...4_:A<...2
    8   88  #    8 | ...1>A:...8_:[..G]G-----C[..C]:...3_:A<...2
    9   87  #    9 | ...1>A:...9_:[..C]C-----G[..G]:...2_:A<...2
   10   86  #   10 | ...1>A:..10_:[..C]C-----G[..G]:...1_:A<...2
##### Criteria: 4.00  0.00  15.00  2.50  65.00  4.50  7.50   [ O N C]
##### 2 non-Watson-Crick base-pairs, and 1 helix (0 isolated bps)
##### Helix #1 (10): 1 - 10

The Jmol paper: a paradigm shift in crystallographic visualization

2010-09-05T21:19:00.001-04:00

From the Jmol mailing list, I noticed Prof. Robert Hanson's announcement titled "Jmol and crystallography: the paper". Finally, an "official" paper has been published on Jmol, an increasing popular molecular visualization tool.

While the paper is categorized under "teaching and education" and summarized as below:

Recent advances in molecular and crystallographic visualization methods are allowing instructors unprecedented opportunities to enhance student learning using virtual models within a familiar web-browser context. In step with these advances, the latest versions of the Jmol molecular visualization applet offer capabilities that hold potential for revolutionizing the way students learn about symmetry, uncertainty and the overall enterprise of molecular structure determination.

Jmol is certainly as useful for research as other commonly used molecular visualization tools, such as PyMOL and RasMol. Specifically, Jmol applets have become dominated in web-based applications/resources (including RCSB PDB). In fact, the HTML version of the paper itself makes extensive use of Jmol applet for interactive manipulation of over 20 figures.

Over the past couple of years, I have witnessed the many new features added to Jmol, driven by an active and friendly user community. I have been particularly impressed by Prof. Hanson's quick and to-the-point responses to users' requests for new features and bug reports, including to my request for Jmol's "Support of 'alchemy' format for rectangular schematic base-pair geometry". Overall, the paper presents a nice overview of some of Jmol's new functionality, and could become highly cited.

Reading through the paper, I was just a bit surprised that PyMOL is not mentioned. Nowadays, Jmol and PyMOL are apparently the top two molecular graphics software tools, with common and complementary functionality.

Calculation of the chi (χ) torsion angle of pseudouridines in RNA structures

2010-08-21T22:50:00.008-04:00

In RNA structures, chi (χ) is defined as the torsion angle about the N-glycosidic bond between the ribose sugar and the canonical A/C/G/U bases or their modified variants. Specifically, for pyrimidines (C and U), χ is defined by O4'-C1'-N1-C2; and for purines (A and G) by O4'-C1'-N9-C4. See the following figure from the IMB Jena Image Library:

Pseudouridine (5-ribosyluracil, PSU) was the first identified modified nucleoside in RNA and is the most abundant. PSU is unique in that it has a C-glycosidic bond instead of the N-glycosidic bond common to all other nucleosides, canonical or modified. It thus poses a problem as to how to calculate the χ torsion angle: should it be O4'-C1'-C5-C4 reflecting the actual glycosidic bond connection, or should the conventional definition O4'-C1'-N1-C2 still be applied literally? As a concrete example, the figure below shows the (slightly) different numerical values, as given by the two definitions, for PSU 6 on chain A of PDB/NDB entry 3cgp/ar0093.

Needless to say, definition of the χ torsion angle of PSU in RNA structures is a very subtle/minor point, and I am not aware of any discussion on this issue in literature (I'd appreciate your sharing of related information in the comment). In 3DNA, PSU is identified explicitly, and χ is defined by O4'-C1'-C5-C4. In NDB and a couple of other tools (I've played with), χ for PSU is defined by O4'-C1'-N1-C2. Again using 3cgp/ar0093 (figure above) as an example, 3DNA gives -162.7°, whilst NDB gives -163.9°. Hopefully, this post will help clarify a confusion for those who care about such little details.

PDB format, how many variants are there?

2010-08-14T22:13:00.005-04:00

In the field of structural biology of macromolecules (proteins and nucleic acids), the PDB format is still the primary one, even though it has shown its age/limitations, and mmCIF and PDBML have been introduced as alternatives.

In 3DNA, PDB is currently the only accepted format for the analysis routines (analyze/cehs), the primary format for rebuilt structures (rebuild/fiber), and various other output structure files (e.g., bestpairs.pdb etc.). By and large, 3DNA follows the PDB format documented at the RCSB website. In my maintaining and supporting of 3DNA and the forum over the years , I have come across other variants of the PDB format. Software packages customize PDB in subtle, sometimes substantial, ways that cause problems with 3DNA. Some examples follow:

Recently, a user tried to analyze an RNA duplex using w3DNA, but found that "the server seems not to recognize the file". It turned out that the so-called PDB format does not have residue name/number and chain id aligned in columns properly. For example, among other issues, the chain id in the problematic PDB is not at column #22 as specified by the standard.

A couple of months back, another user asked how to handle AMBER-introduced ribonucleotide abreviations (RA/RC/RG/RU). Compared to the above "PDB" variant, the AMBER PDB is well-behavored, and can be easily recognized by 3DNA after introducing corresponding one letter matches in file baselist.dat. Interestingly, a couple of years ago, DA/DC/DG/DT was introduced into the PDB format to distinguish deoxyribonucleotide (DNA) from ribonucleotide A/C/G/U (RNA). The AMBER convention, by using RA/RC/RG/RU for RNA, appears to be more symmetric.

From my (limited) experience in helping a 3DNA user several ago, I know that CHARMM uses three-letter residue names for DNA (e.g., ADE/CYT/GUA/THY). There is actually a short Perl utility program (x3dna2charmm_pdb) in 3DNA for converting 3DNA generated PDB file to that recognized by CHARMM.

As noted in another post, "PDB ATOM coordinates record", Babel/Open Babel converted PDB files (mostly of small molecules) do not follow atom naming conventions.

As widely used and standard as PDB format is, there are still many variants in common use. Thus, when you meet a problem related to a PDB file, it helps to check if it conforms to the standard. With 3DNA, I've always tried to incorporate as many PDB variants as possible, so users can run 3DNA on them transparently.

Calculation of DNA bending angle

2010-08-07T19:34:00.008-04:00

Once in a while, I receive the question (via email or on the 3DNA forum) on how to calculate DNA bending angle, mostly associated with DNA-protein complexes. This question certainly qualifies as an FAQ; Indeed, I did consider writing an entry on it in the 3DNA website. However, I've hesitated to do so because the topic is more subtle and complicated that it appears.

On its face, an angle is defined by two vectors; let's call them a and b, and if each is normalized, then the angle (in degrees) between them is: acos(dot(a, b)) * 180/pi. Geometrically, after moving the tails of the vectors into the same position (e.g., origin), the heads would normally define a plane, unless a and b are strictly parallel (0°) or anti-parallel (180°).

DNA structures are three-dimensional, normally far more complicated than a single number can quantify. The concept of DNA bending angle, as I understand it, is only applicable to DNA structures with two relatively straight fragments (as in CAP-DNA complexes). Under such situations, one can fit a least-squares (ls) linear helical axis to each of the two fragments, and calculate the angle between them. Towards this end, 3DNA outputs the following section when it judges that the input structure is not strongly curved. Using 355d/bdl084, which is distributed with 3DNA, as an example:

Global linear helical axis defined by equivalent C1' and RN9/YN1 atom pairs                                    
Deviation from regular linear helix: 3.30(0.52)
Helix:    -0.127  -0.275  -0.953
HETATM 9998  XS    X X 999      17.536  25.713  25.665
HETATM 9999  XE    X X 999      12.911  15.677  -9.080
Average and standard deviation of helix radius:
P: 9.42(0.82), O4': 6.37(0.85),  C1': 5.85(0.86)

Where the Helix: line gives the normalized vector along the "best-fit" helical axis. The two HETATM records provides the two end points on the helix, and they are directly related to the Helix: line by a simple equation. Following the above example, we have (Octave/Matlab code):


XE = [12.911  15.677  -9.080];
XS = [17.536  25.713  25.665];

dd = XE - XS
%   -4.6250  -10.0360  -34.7450

Helix = dd / norm(dd)
%  -0.12685  -0.27526  -0.95296  ==> [-0.127  -0.275  -0.953]

With the two HETATM records, one can easily add them into the original PDB file to display the helical axis using a molecular graphics programs (e.g., RasMol, Jmol, PyMol). Moreover, the two helix vectors can be used to reorient the original PDB structure into a view so that one helical fragment lies along the x-axis, and the other in the xy-plane. As documented in detail in recipes #4 on "Automatic identification of double-helical regions in a DNA–RNA junction" (2008 3DNA Nature Protocols paper), "The chosen view allows for easy visualization and protractor measurement of the overall bending angle between the two relatively straight helices."

The following points are well worth noting:

The ls fitting procedure used in 3DNA follows SCHNArP, which was based on the algorithm in the well-known NewHelix program, maintained by Dr. Richard Dickerson unto the 1990s. While fitting a global linear helical axis to strongly curved DNA structures clearly makes no much sense with derived parameters (NewHelix itself has been replaced by FreeHelix, also from Dickerson), I do believe it is meaningful to fit a linear helix with relatively straight DNA fragment. That's why I have kept this functionality in SCHNArP and 3DNA; 3DNA bending angle calculation serves as an example illustrating the point – it provides an "intuitive" way for biologists to understand how the bending angle is calculated; it can actually be measured directly.
Instead of directly ls-fitting a linear helical axis with 3DNA, one can alternatively superimpose a regular fiber model into the DNA fragment, and then derive the straight helical axis from the fitted coordinates. The two approaches normally gives slightly different numerical values, as would be expected.
Overall, bending angle is (at most) an approximate measure of DNA curvature. In my opinion, the concept is only applicable for comparing a set of structures, each with two relatively straight fragments. Even in such cases, the relative spatial relationship between two helical fragments are more complicated than a simple (bending) angle could quantify. Specially, it is not meaningful to make a selling point in slight differences of DNA bending angles.

Genomic code for nucleosome positioning?

2010-08-06T21:08:00.004-04:00

Over the past few years, the research on nucleosome positioning has been a hot topic: Segal et al., through a series of high-profile publications, have proposed that "the intrinsic nucleosome sequence preferences contribute substantially to nucleosome organization and chromatin function in vivo" (i.e., there exists a "genomic code") and developed bioinformatics tools to predict patterns of nucleosome positions based on genomic DNA sequence. However, as is clear from the title of a paper by Zhang et al., "Intrinsic histone-DNA interactions are not the major determinant of nucleosome positions in vivo", it is highly debatable if the nucleosome sequence preferences is strong enough to be called a "code", in the canonical sense as for the genetic code.

While I am (generally) aware of the controversy about the "genomic code", my understanding of the issue has become much clearer after reading the correspondences between both sides published in the 17(8), August 2010 issue of Nature Structural & Molecular Biology [NSMB]:

"Nucleosome sequence preferences influence in vivo nucleosome organization" (pp918 - 920) by Kaplan … and Segal
"Evidence against a genomic code for nucleosome positioning" Reply to “Nucleosome sequence preferences influence in vivo nucleosome organization” (pp920 - 923) by Zhang … and Struhl
"A preoccupied position on nucleosomes" (p923) by Pugh

Reading carefully through the text, it is clear how important the little details could be: 20-bp vs 40-bp windows? How significant the preference should be to be claimed as a "code"? Pugh clarifies another source of confusion, in the terminology of nucleosome ‘occupancy’ vs ‘positioning’. In my understanding, the dispute on such details itself speaks strongly against the existence of "genomic code".

It should be noted that early work by Drew & Travers (JMB85) and Satchwell et al. (JMB86) provided evidence that the rotational positioning of DNA about the histone octamer is regulated by DNA sequence-dependent bending preferences. Goodsell & Dickerson quantified the Satchwell trinucleotide model by assigning a set of roll angles to predict DNA bending from base sequence. This bending model, along with five others, was implemented by SCHNArP for building approximate 3D DNA structures (with fixed values for slide, twist etc other parameters). However, this functionality is dropped out from 3DNA, simply because of the arbitrariness/subtleties in picking up the parameters.

Principle axes of inertia as a mean to reorient a molecule

2010-08-01T21:22:00.007-04:00

In structural biology, it is often helpful to reorient a molecule in a specific view for easy visualization and comparison. With no a priori knowledge of a feature, the (objective) principle axes of inertia can be used as a mean to automatically reset a structure into its most extended orientation (and two other orthogonal views). In 3DNA, rotate_mol has such a functionality which is explored by blocview to generate the characteristic, simplified molecular images adopted by the RCSB PDB and NDB.

Recently, the thread (initiated by Dr. Pascal Auffinger) at the 3DNA forum titled "rotate_mol ambiguous orientations" has led me to think the topic deeper. As always, the issue is best illustrated with an example. The file 355d_A6_T19.pdb contains the coordinates of the base pair A6-T19 of PDB entry 355d. Let xyz be a variable in Octave (Matlab) containing the coordinates of all atoms in file 355d_A6_T19.pdb, the following m-file provides the detailed procedure:

xyz = [
   13.954  12.080  15.071
   13.092  11.246  15.953
   14.769  11.453  14.020
   13.056  13.202  14.406
   12.184  13.965  15.222
   11.308  14.846  14.376
   12.033  15.991  13.896
   10.654  14.188  13.168
    9.262  14.480  13.256
   11.373  14.820  11.974
   11.810  16.185  12.511
   13.051  16.742  11.953
   14.253  16.097  11.802
   15.197  16.853  11.305
   14.585  18.084  11.106
   15.066  19.322  10.605
   16.329  19.520  10.209
   14.198  20.357  10.535
   12.938  20.162  10.946
   12.370  19.048  11.444
   13.259  18.034  11.497
   14.865  29.619   5.710
   14.172  30.844   5.281
   16.107  29.695   6.515
   13.846  28.701   6.510
   12.684  28.174   5.864
   11.909  27.311   6.826
   12.626  26.090   7.129
   11.584  27.960   8.173
   10.189  27.833   8.395
   12.409  27.166   9.184
   12.544  25.808   8.521
   13.716  24.979   8.867
   13.496  23.688   9.310
   12.381  23.202   9.467
   14.631  22.971   9.551
   15.943  23.388   9.393
   16.852  22.613   9.632
   16.107  24.759   8.938
   17.488  25.309   8.745
   14.995  25.475   8.711 ]

cov_mtx = cov(xyz)

% cov_mtx =
%     3.4094    2.3860   -1.4222
%     2.3860   33.7089  -15.7175
%    -1.4222  -15.7175    7.8786

[v, lambda] = eig(cov_mtx)

% v =
%    0.094830   0.992835  -0.072709
%    0.419520  -0.106092  -0.901525
%    0.902779  -0.054989   0.426575

% lambda =
%     0.42535          0          0
%           0    3.23319          0
%           0          0   41.33836

x = v(:, 2);
y = v(:, 3);
z = cross(x, y);

% z =
%   -0.094830
%   -0.419520
%   -0.902779

num = size(xyz, 1);
xyz2 = (xyz - ones(num, 1) * mean(xyz)) * [x y z]

% xyz2 =
%     1.148136   10.032849   -0.514731
%     0.332293   11.223635   -0.879359
%     2.081609   10.090517    0.619843
%    .................................
%     2.277813   -4.170312   -0.501254
%     3.601180   -4.848890   -0.688714
%     1.110301   -4.831784   -0.491249

Among other subtitles, note that:

lambda 3 has the largest value (41.33836), i.e., with the biggest variance, and its corresponding eigenvector is set as the vertical y-axis.
lambda 2 shows intermediate variation (3.23319), and its corresponding eigenvector is set as the horizontal x-axis.
The z-axis is derived as the cross-product of x- and y-axes to ensure a right-handed frame. Note that it is the opposite of eigenvector 1.
The original xyz coordinates are moved to its geometric center [xyz - ones(num, 1) * mean(xyz)], and then projected onto the preferred frame to get the transformed coordinates (xyz2) with the molecule in the most extended view (vertical then horizontal).

3DNA for the analysis of molecular dynamics simulations

2010-07-24T22:06:00.006-04:00

To my surprise and delight, 3DNA has been increasingly used in the analysis of trajectories of molecular dynamics (MD) simulations of nucleic acid structures. While I have heard of CHARM/AMBER for a long time and know of the basic principles of MD simulations, I have never performed any simulations (yet). Over the years, I have been mostly interested in methods development and algorithms implementation for better understanding of nucleic acids structures, and have always used NDB/PDB entries for my tests/applications. Thus unlike Curves+, 3DNA lacks direct support for the analysis of MD trajectories.

From the MD-related questions I have received over emails and on the 3DNA forum, I sense that users do not have that much of a problem or inconvenience in applying 3DNA for the task. There is one issue that does stand out, however, i.e., the use of find_pair for each and every structure in the trajectories to prepare an input file that contains base-pairing information for analyze to calculate various parameters. As an example, a recent article by Ricci et al., "Molecular Dynamics of DNA: Comparison of Force Fields and Terminal Nucleotide Definitions" (J Phys Chem B. 2010 Jul 8. [Epub ahead of print]), says that:

It is noteworthy that, because the DNA structure endured strong deformations in these simulations, 3DNA could not perform the analysis of the angles for all the 5000 snapshots. Actually, for system 3 the calculation of the backbone parameters could be performed in only 59% of the snapshots, whereas in system 4 only 18% of the configurations allowed this calculation.

Dr. Netz, the corresponding author, confirmed to me that:

I don't think this is a 3DNA problem, but is rather an effect of the distortion caused by the problems in parameterization of the GROMOS forcefield.

Essentially, though, the problem lies in the coupled use of find_pair/analyze; for strongly deformed DNA structures along the simulation trajectories, some base pairs melt out, and are thus beyond the recognition of find_pair. On the other hand, one does not need to run find_pair for each structure; simply prepare a template input file for analyze, with the help of find_pair and probably some manual editing, then all the structures can be analyzed, no matter how distorted they are. Furthermore, saving the find_pair step should speed up the analysis of thousands of structures, typical of MD simulations.

In fact, 3DNA (from v1.5) has a Perl script named nmr_strs that implements the above idea in a preliminary form. The script is intended to serve as a starting point for the analysis of multiple similar structures, in an NMR ensemble or from MD simulations. See 3DNA wiki page on nmr_strs for details.

In the long run, a better integration between 3DNA and AMBER or other MD packages would solve this problem from the upstream, thus making users' life (downstream) much easier. Until that happens, I hope this post could help clarify the issue somewhat in this important application area of 3DNA. As always, I greatly value your comments.

PDB fingerprints make salient features of PDB entries obvious at a glance

2010-07-16T22:07:00.011-04:00

In the past couple of days, I have followed with great interest the thread "PDBprints - salient, at-a-glance info about PDB entries" in CCP4BB. The idea behind PDBprints (short for PDB finger prints) is very simple: now with over 66 thousands entries, no one could possibly be intimately familiar with many of the structures in PDB. Yet "it would be very useful to have a way of obtaining at-a-glance information about the salient features of one or more PDB entries." In this scheme, key features (e.g., experimental method, species, protein/NA/ligand components, published or not) are represented by stylized icons called PDBlogos; a set of which organized in a specific order gives PDBprints.

The initial announcement, posted on Thursday by Gerard Kleywegt (the current Director of PDBe), has been well-received by the community. I am very impressed by the quality of the (over a dozen) feedbacks on icons coloring for presence/absence of a feature, and the proper icons for species etc. On the other hand, I am surprised by the quick and concrete response Gerard posted today, addressing the various issues.

Over the years, I have taken PDB/NDB mostly as data repositories; I regularly download (update) nucleic-acid containing structures for local analyses. Among the three wwPDB sites (RCSB, PDBe and PDBj), I've only visited the RCSB site on a regular basis. Now I find more features are available from PDBe, which I will browse more in the future.

In 3DNA, analyze and rebuild are two sides of the same coin

2010-07-03T23:46:00.003-04:00

Recently, I noticed the following two papers:

"A measure of bending in nucleic acids structures applied to A-tract DNA" by Lankas et al. (2010), Nucleic Acids Res.,38(10), 3414–3422.
"The chirality of DNA: elasticity cross-terms at base-pair level including A-tracts and the influence of ionic strength" by Noy and Golestanian (2010), J. Phys. Chem. B, 114, 8022–8031.

They both cite 3DNA, where the combined usage of its analyze/rebuild components has played a significant role.

When I first had access to the CEHS scheme, I was immediately attracted by its mathematical rigor which allows for complete reversibility in DNA structural analysis and model rebuilding. Historically, CEHS was the initial seed of the SCHNAaP/SCHNArP programs, and analyze/rebuild in 3DNA were directly derived from them.

Frequently, I think of analyze/rebuild as two sides of the same coin: starting from a nucleic acid structure (e.g., a DNA double helix), the analyze program gives a set of base-pair (propeller, buckle etc) and step (roll, slide etc) parameters. The structural parameters can be used to rebuild the structure, which is virtually identical in base geometry (i.e., without taking consideration of the sugar-phosphate backbone) to the original structure. Conversely, analyzing the rebuilt structure again gives exactly the same set of structural parameters describing the relative base geometry.

Reversibility is essentially a simple concept, nevertheless a very powerful one. Over the years, I am glad to see more people are taking advantage of this 3DNA unique feature in addressing real-world problems related to nucleic acid structures. I am confident to see more such applications.

Get all torsion angles in a nucleic acid structure

2010-06-27T23:20:00.002-04:00

In the field of nucleic acid structure analysis, a commonly calculated set of parameters is the torsion angles. Specifically, they include the main chain and chi (χ) torsion angles defined as follows:

α: O3'(i-1)-P-O5'-C5'
β: P-O5'-C5'-C4'
γ: O5'-C5'-C4'-C3'
δ: C5'-C4'-C3'-O3'
ε: C4'-C3'-O3'-P(i+1)
ζ: C3'-O3'-P(i+1)-O5'(i+1)
χ: for pyrimidines (Y, i.e., T, U, C): O4'-C1'-N1-C2; for purines (R, i.e., A, G): O4'-C1'-N9-C4

A related set of parameters characterizes the sugar conformation:

ν0: C4'-O4'-C1'-C2'
ν1: O4'-C1'-C2'-C3'
ν2: C1'-C2'-C3'-C4'
ν3: C2'-C3'-C4'-O4'
ν4: C3'-C4'-O4'-C1'
tm: amplitude of pseudorotation of the sugar ring
P: phase angle of pseudorotation of the sugar ring

These torsion angles are clearly defined and are readily available from various informatics programs. Not surprisingly, 3DNA also provides a complete set of DNA/RNA backbone torsions, calculated robustly and efficiently. The key is the "-s" option of find_pair" program, which is nevertheless little used, mostly because it is not the default.

Using the Haloarcula marismortui 50S large ribosomal subunit as an example (1jj2), the output file 1jj2.outs from the following command contains all the above mentioned main chain and sugar conformational parameters:

find_pair -s 1jj2.pdb stdout | analyze

Please see the Jena "Nucleic acid backbone parameters" website for a diagram (based on Saenger's book) defining the various backbone torsion angles.

Subscribing to mailing lists in daily digest mode

2010-06-20T15:48:00.004-04:00

I am currently on quite a few mailing lists that are relevant to my research topics in a broad sense, including Jmol, pdb-l and CCP4. Over the years, I've found it very handy to keep informed of a field by following in its (major) mailing list. However, I could easily get lost with so many posts each day on some active lists. So subscribing to the daily digest mode, provided by many mailing list management software, has become a norm. Thus, I would receive only one (or a few) email(s) from a mailing list. I could then easily browse through the subject lines to decide if to read further of a specific topic.

Among the three lists mentioned above, Jmol is quite active with several core contributors, most impressively from Prof. Robert Hanson. The PDB mailing list is of unbelievably low volume, with less than a post per day. CCP4 bullet board is overall the most well-organized; I am surprised by the knowledge-base from the posts there – it is a great resource in structural biology.

AMBER mailing list and Computational Chemistry List (CCL) are the other two lists I subscribed to before. However, I could not get into the daily digest mode; soon my mail box was flooded by many unrelated posts so I had to un-subscribe from them.

Journal impact factor and individual researcher h-index

2010-06-20T13:07:00.005-04:00

The June 17, 2010 issue of Nature (Vol. 465, No. 7300) has extensive discussions on assessing the impact (influence/significance) of a journal or an individual researcher using quantitative metrics. Not surprisingly, there is no consensus. Over the past 20 years, the field of bibliometrics (scientometrics) has seen a ten-fold explosion in publications. In fact, the title of the editorial is "Assessing assessment"; it boils down to whether we should use a quantitative metric, a combination of several such metrics, or which metrics to choose among the so many possibilities.

Among the various view points, I agree more with David Pendlebury, Citation Analyst of Thomson Reuters. "No one enjoys being measured — unless he or she comes out on top." "Importantly, publication-based metrics provide an objective counterweight ... to bias of many kinds." "there are dangers [to] put too much faith in them. A quantitative profile should always be used to foster discussion, rather than to end it."

Among the many currently available metrics, it seems fair to say:

impact factor (IF), introduced in 1963, is the most important one to measure the impact of a journal. Nowadays, it is common to see the IF of a journal at its home page. For example, currently Nucleic Acids Research has an IF of 6.878. However, as emphasized by Anthony van Raan, director of the Centre for Science and Technology Studies at Leiden University in the Netherlands: "You should never use the journal impact factor to evaluate research performance for an article or for an individual — that is a mortal sin." Indeed, "In 2005, 89% of Nature’s impact factor was
generated by 25% of the articles."
h-index, introduced in 2005 by Hirsch, is currently the most influential metric to quantify the productivity and impact of an individual researcher. An h-index is "defined as the number of papers with citation number ≥h". It "gives an estimate of the importance, significance, and broad impact of a scientist’s cumulative research contributions." I found it very informative to read the original, 4-page long PNAS paper to understand clearly why h-index is defined that way and to be aware of some of its caveats.
number of citations is undoubtedly the most objective metric to measure the impact of a publication.

Overall, it is easy to come up with a quantitative measure; the point is what that number really means. Along the line, it is of crucial importance to know exactly how a metric is calculated. No metric could be perfect; however, a metric, defined transparently and applied consistently, is more objective and convincing than other means.

3DNA citations from mainland Chinese researchers

2010-06-11T23:11:00.008-04:00

Over the past couple of weeks, I have been quite surprised to notice the following four papers, by researchers from mainland China, that cite 3DNA:

Feng et al. "Crystal structure of the crenarchaeal conserved chromatin protein Cren7 and double-stranded DNA complex." Protein Sci. 2010 Jun; 19(6):1253-7.
Zhang et al. "Structural insights into the interaction of the crenarchaeal chromatin protein Cren7 with DNA." Mol Microbiol. 2010 Mar 16. [Epub ahead of print]
Li et al. "Molecular dynamics studies of the 3D structure and planar ligand binding of a quadruplex dimer." J Mol Model. 2010 May 29. [Epub ahead of print]
Chen et al. "Theoretical studies on the stacking interactions between methylated and unmethylated DNA bases." ACTA CHIMICA SINICA (in Chinese) 68(8):739-46. Published: Apr. 28, 2010.

For more ten years, this is the first time (to the best of my knowledge) that 3DNA has been cited by scientists from China, and it so happens that four come in a row. I am so happy to see that structural biology is shaping up in China. Surely, I would expect more 3DNA-citing papers from China to come in the following years.

block_atom, a little used 3DNA utility script

2010-06-11T22:25:00.006-04:00

Following my previous post, "How is a base-pair rectangular block defined in 3DNA?", I feel it is in order to mention a tiny Perl script named block_atom that seems to have been little used by the 3DNA user community.

The basic idea of block_atom is to get an ALCHEMY data file which combines both atomic and block representations of nucleic-acid-containing structure in PDB format. The ALCHEMY file can then be displayed interactively using RasMol (v2.6.4, by its original author) (or Jmol).

Again, the point is best illustrated with an example. Here I am using 355d, the famous Drew-Dickerson B-DNA dodecamer solved at a high resolution by Williams and colleagues. Let the input PDB file be 355d.pdb, then the output file is automatically named 355d.alc, which can be displayed using RasMol with the options '-alchemy -noconnect'.

block_atom 355d.pdb
rasmol -alchemy -noconnect 355d.alc

One view of RasMol rendered image is as below, which shows more clearly the minor-groove (colored blue) of the DNA duplex and the various intra-base-pair deformations (e.g., propeller and buckle):
I initially wrote this utility script mostly for verification purpose, i.e., to make sure that the base reference frames are defined properly. For those who want to know more about 3DNA, having a look of the script (which is tiny) and understanding how it works could be a good exercise.

How is a base-pair rectangular block defined in 3DNA?

2010-06-05T22:27:00.009-04:00

One of 3DNA's unique features is the base-pair (and base) rectangular blocks, as shown in the figure below. Since such a schematic was initialized by Calladine and Drew in their Understanding DNA book, I normally refer it as the Calladine-Drew style representation.

By default, a base-pair [BP, (a)] has a dimension of 10x4.5x0.5 Å; a purine [R, (b) left] 4.5x4.5x0.5 Å; a pyrimidine [Y, (b) right] 3x4.5x0.5 Å; and a mean base [M, (c)], which is exactly half of the base-pair, 5x4.5x0.5 Å.

The blocks are put into specific files: Block_BP.alc, Block_R.alc, and Block_Y.alc respectively. To use M for R and Y, one just needs to copy file Block_M.alc to overwrite Block_R.alc and Block_Y.alc in the current directory or the system directory ($X3DNA/config/). These blocks are used in the rebuilding and visualization components of 3DNA (see my blog post "blocview: a simple, effective visualization tool for nucleic acid structures").

The blocks are stored in the ALCHEMY format for easy specification of the nodes and edges, and since ALCHEMY is a format supported by RasMol (and Jmol). As an example, Block_BP.alc has the following content:

   12 ATOMS,    12 BONDS
    1 N      -2.2500   5.0000   0.2500
    2 N      -2.2500   0.5000   0.2500
    3 N      -2.2500   0.5000  -0.2500
    4 N      -2.2500   5.0000  -0.2500
    5 C       2.2500   3.5000   0.2500
    6 C       2.2500   0.5000   0.2500
    7 C       2.2500   0.5000  -0.2500
    8 C       2.2500   3.5000  -0.2500
    9 C      -2.2500   5.0000   0.2500
   10 C      -2.2500   0.5000   0.2500
   11 C      -2.2500   0.5000  -0.2500
   12 C      -2.2500   5.0000  -0.2500
    1     1     2
    2     2     3
    3     3     4
    4     4     1
    5     5     6
    6     6     7
    7     7     8
    8     5     8
    9     9     5
   10    10     6
   11    11     7
   12    12     8

Astute viewers may observe that nodes 1-4 are identified as N (nitrogen) and have exactly the same coordinates as 9-12 (C, carbon). This is a trick so RasMol can display the the minor groove edge in a different color than the other five sides of the rectangular, as shown in the following figure:

Note that the rectangular is preset in the standard base reference frame. Thus the nodes have y-coordinates of +5 Å and -5 Å along the long edge of the base pair, and x-coordinates of +2.25 Å and -2.25 Å along the short edge.

As an extra bonus of storing the blocks in external ALCHEMY text files, the dimensions of the blocks are be readily changed. For example, the depth of a block (z-coordinates) can be easily increased from 0.5 to 1.0 Å to make it thicker (see $X3DNA/config/Block_BP1.alc). Moreover, the blocks do not need to be rectangular either – see $X3DNA/config/block/Block_R_nr.alc for an example.

Naming conventions: RasMol, PyMOL and Jmol

2010-05-28T22:22:00.003-04:00

Over the years, I've used the following three molecular graphics programs the most: RasMol, PyMOL and Jmol. It is interesting to note their different naming conventions: while their names all end with m-o-l (for molecule, I believe), the cases (used in convention) are different: Mol as in RasMol, MOL as in PyMOL, and mol as in Jmol.

For Jmol, there is actually an FAQ item "How do you write Jmol?" that reads "Capital J, lower case mol. Please do not write it any other way." For PyMOL, I am not aware of a specification on its name – I know that's the "official name" from its home page. Regarding RasMol, that's the most common way to name it, even though the title of its 1995 publication is "RASMOL: biomolecular graphics for all."

As for the prefix in the names, J in Jmol stands for Java; Py in PyMOL for Python; and Ras in RasMol could well be the initials of Roger A. Sayle, the original author of the program.

MODEL/ENDMDL records in PDB X-ray crystal structures

2010-05-23T23:08:00.007-04:00

Over the years, I have always been using the PDB format for nucleic acid structures. Naturally I've thought that I know the format well, especially the "Coordinate Section." According to the PDB documentation,

MODEL/ENDMDL records are used only when more than one structure is presented in the entry, as is often the case with NMR entries.

In my experience, I have always connected the MODEL/ENDMDL pair only with the different models in an NMR entry. However, I was recently bitten by a subtlety in PDB format that is related to the MODEL/ENDMDL records in X-ray crystal structures which contain more than one asymmetric unit in their biological assembly.

As always, the point is best illustrated with an example – here I am using the four-stranded DNA Holliday junction X-ray crystal structure (1zf4/ud0061) solved by Ho and colleagues [PNAS 2005 May 17;102(20):7157-62]. As shown in the NDB website, the asymmetric unit of 1zf4/ud0061 contains only two chains; it is the biological assembly that has the four-stranded DNA junction.

I downloaded the biological assembly coordinates (in PDB format) from the NDB (named it '1zf4.pdb'), and ran blocview on it: blocview -i=1zf4.png 1zf4.pdb. However, the generated image [see Figure (A) below] has only half of what expected – as if the downloaded file contains only coordinates of an asymmetric unit. The mystery was gone when I checked the PDB file and realized that the MODEL/ENDMDL records now also apply to X-ray crystal structures to delineate symmetric-related units. Since 3DNA is designed to handle only one structure – it stops processing whenever an END or ENDMDL record is encountered. Simply change each 'ENDMDL' record to ' ENDMDL' (i.e., adding a space, or any character, for that matter) and run blocview again will get the expected image [Figure (B)].

Note that by default, RasMol also only displays the first model in a PDB file. To see multiple structures, the option '-nmrpdb' must be specified in the command line.

3DNA forum registration -- your id

2010-05-15T17:59:00.003-04:00

Over the past couple of months, I have received a few emails with a subject line, as the title of this post, "3DNA forum registration -- your id". Presumably, such an email is from a (potential) 3DNA user, requesting for activation after 3DNA forum registration.

Clearly, "your id" should have been a specified user id. In a couple of cases, I tried to figure out the corresponding 3DNA forum ids based on email addresses and activated them. More recently, I've switched to ask the sender a very simple question: "So what is 'your id' specifically?" In at least two cases, I have never heard back from the original senders with 'your id' specified. To a certain extent, this is a surprising result: providing such information won't take more than a couple of minutes, and this will help me (and others) to help them better with 3DNA-related questions. I can certainly understand occasional negligence, but users must follow some common-sense rules. Thus those registrations are deleted, as junk ones, after some grace period.

The 3DNA forum is open to public for browsing, without registration. At its current settings, no emails are sent from the forum; do not confuse an online forum with a mailing list. You only need to register if you want (and are certainly more than welcome) to get more actively involved in, e.g. asking/answering questions, sharing tips/tricks with other users. These policies are enforced to make the 3DNA forum free from spams, as much as possible.

Some key combinations in MacBook Pro (Snow Leopard)

2010-05-08T08:40:00.005-04:00

One of the inconveniences I experienced when first switching to a MacBook Pro running Mac OS X (Snow Leopard) was its missing of the PgUp, PgDn, Home and End keys. Over the past few months, I have learned that the functionality of such 'convenience' keys can be achieved via a combination of the arrow-keys (bottom right), with the 'fn', 'option', or 'command' keys (bottom left), as follows:

fn + right-arrow – end of a document
fn + left-arrow – beginning of a document
fn + up-arrow – Page up
fn + down-arrow – Page down
command + right-arrow – beginning of a line (Home)
command + left-arrow – end of a line (End)
option + right-arrow – one word to the right
option + left-arrow – one word to the left

According to Wikipedia, the 'fn' key "is a modifier key on many keyboards, especially on laptops, used in a compact layout to combine keys which are usually kept separate."

In addition to modifying cursor movements as noted above, the 'fn' key in Mac OS X can also be used to switch the default functionalities of F1 to F12 (e.g., F1 and F2 for screen brightness control) to the standard function keys. As an example, in MS Word, the keyboard shortcut for toggling case is "shift + f3", a trick I recently learned. In the default setting, one must press "fn + shift + f3" to achieve the desired effect.

Any tricks to share? I'd like to hear them!

One year of blogging

2010-05-01T14:01:00.005-04:00

When I checked the date of my first blog post today, I was a bit surprised to find that it is exactly one year since I begin to blog on May 2, 2009. Altogether, I have written over 60 posts, slightly more than one per week. At this time, I feel it appropriate to summarize my thought on blogging in general, to provide a perspective to those who care to visit here and make comments.

Why? The initial motivation was to use blog as a platform to express my personal views on issues I am interested in. As made clear in my first post, "this is Xiang-Jun's Corner on the Internet: all views are mine, and I am opinionated." Over the time, the blog posts have served as a convenient notebook (searchable and archived) either for my personal reference, or to refer others to a particular post (e.g., "On maintaining the 3DNA forum" when being asked a 3DNA-related question via email).
What? "Random thoughts, mostly on scientific issues". I write only on issues I am familiar with and feel comfortable to say something, to the limit that I can respond quickly and concretely to users comments. I am always open to suggestions and will be prompt in acknowledging errors and making corrections. So far, the largest portion of the posts has been devoted to nucleic acid structures in general, and 3DNA-related topics in particular.
How often? Due to time constraints, I will try to write one post per week at the minimum, maybe two, or in rare occasion three. Exceptions are possible, but by and large, I will aim for ~100 posts per year.
Comment? Currently, the policy is set such that "Anyone - includes Anonymous Users" can make a comment, and the comments are moderated. So far, I have always approved all the comments as soon as I see them in my gmail alert, and follow up where appropriate. Note that due to global time difference, commenters may experience some lag in time.
Does it work? Not unexpectedly, my blog has gradually attracted attentions from quite a broad audience, especially those interested in nucleic acid structures (3DNA), including leading scientists in the field (I know from emails I have received).
Hot posts? According to Google Analytics, the most frequently visited nine posts are as follows (with posting date in parentheses):
1. Curves+ vs 3DNA (Sunday, August 16, 2009)
2. Does 3DNA work for RNA? (Friday, July 10, 2009)
3. Two web-interfaces to 3DNA, and more (Sunday, July 5, 2009)
4. Fit a least squares plane to a set of points (Saturday, August 22, 2009)
5. Two 3DNA figures made into a textbook on structural biology (Sunday, May 3, 2009)
6. Chemical diagram of Watson-Crick base-pairs (Saturday, January 23, 2010)
7. What's special about the GpU dinucleotide platform? (Friday, April 2, 2010)
8. Double helix groove width parameters from 3DNA (Saturday, September 5, 2009)
9. How to calculate torsion angle? (Saturday, October 31, 2009)
While some of the posts are well-expected to be in the list, a few of them (e.g., ls-plane fitting, calculation of torsion angle) could look a bit surprising. It does, however, verify an observation based on my personal experience and intuitive feeling about a technical niche that my expertise can make a difference.

Again, as I wrote in my first blog post, "Now the ball is rolling, and only time can tell where the destination will be -- but surely it will no longer stand where it was!" One year later, I can confidently say that the ball in rolling in the right direction, as I'd have hoped for. Of course, I know for sure more time and efforts are need to move to the next level, and I value your feedback!

C2B2 annual retreat; DNA shape matters

2010-04-30T23:17:00.012-04:00

As in previous years, today's C2B2 annual retreat provided a good opportunity to know the projects going on, and to meet the people performing the researches in the broad area of computational biology and bioinformatics at Columbia. The depth and breadth of the topics covered, and the close collaborations among different labs were amazing; revealing the multidisciplinary nature of bioinformatics and systems biology. After lunch, there was only slightly more than 30 minutes left for poster session. I quickly walked around the posters and I only wished to have more time to read more carefully some of them and discuss with the presenters.

From a DNA structural perspective, of special note is the plenary talk by Dr. Tom Tullius, titled "Function follow form: the role of DNA shape in the working of human genome". I am aware of Dr. Tullius work on hydroxyl radical footprinting, and his 2009 Science paper on the correlation between local DNA topography (the shape of the backbone and grooves) with functional noncoding regions in the human genome. His talk connects naturally to the work by Rohs et al. (Nature, 2009) on the "The role of DNA shape in protein-DNA recognition" which shows that narrow DNA minor groove (associated with A-tracts and negative electrostatic potential) interacts preferably with arginine residues. The last talk by Dr. Richard Mann, titled "New insights into transcript factors specificity" echoed nicely with Dr. Tullius's presentation, emphasizing the importance of DNA shape.

Along the line, it is very interesting to note that Dr. Tullius mentioned a recently published paper titled "The shape of the DNA minor groove directs binding by the DNA-bending protein Fis" by Stella et al., Genes Dev. 2010 Apr 15; 24(8):814-26. I happened to read this paper just a coupe of days ago and was impressed by the work. Additionally, I am happy to see that 3DNA's analyze/rebuild/fiber components combined together helped to pinpoint the fact that only three parameters – roll, twist, and slide – are needed to explain observed DNA shapes in Fis-DNA complexes. See also a related blog post, titled "Fiber, analyze and rebuild in 3DNA".

Overall, it is nice to see that the community has gradually realized that DNA structural features, not just (or in addition to) its base sequence, are a major player in protein/ligand-DNA recognition. From my (certainly not that extensive) understanding of the field, the impressive list of publications over the past few years by Dr. Remo Rohs (Barry Honig group at Columbia University) has played a significant role in "shaping" DNA.

Life is complicated -- is there a way to make it simpler?

2010-04-23T23:18:00.007-04:00

To mark the 10th anniversary on completion of the draft sequence of human genome, the April 01, 2010 issue of Nature [464 (7289)] published a series of interesting and revealing articles, including an Editorial and (historical) accounts from Francis Collins and Craig Venter. I browsed through the whole list, and I especially liked the News article titled "Life is complicated" by Erika Check Hayden, a senior reporter for Nature. I was attracted by its catchy title and brief summary: "The more biologists look, the more complexity there seems to be. Erika Check Hayden asks if there's a way to make life simpler." Obviously, the title of this blog post was inspired by the sources.

Over the past decade, the Human Genome Project has helped to clarify the number of genes from previously assumed ~100,000 to the "true" number of only ~21,000. The dramatically reduced number of genes illustrates the crucial importance of non-coding (used to be called "junk") DNA to biology, yet what non-coding DNA does is still befuddling. Nowadays, we are faced with data deluge from sequencing, gene expression, protein (transcription factor) binding, and other new technologies. The complexity of biology has grown significantly, instead of simplified, even for the most extensively studied protein p53. As put by Jennifer Doudna, “The more we know, the more we realize there is to know.”

The community has gradually realized that information gathering does not always bring corresponding increase of meaningful biological insights. We are facing an age of “drowning in information, starved for knowledge.” For example, systems biology, a new discipline “supposed to help scientists make sense of the complexity”, has turned out that “In many cases, the models themselves quickly become so complex that they are unlikely to reveal insights about the system, degenerating instead into mazes of interactions that are simply exercises in cataloguing.” I cannot agree more with the comment by Leonid Kruglyak: it is naive to think that “you can simply take very large amounts of data and run a data-mining program and understand what is going on in a generic way.”

Are we lost in the sea of biological data? Not necessarily. Eric Davidson's work is an excellent example in “taking smarter systems approaches” to reveal overarching biological rules. Instead of a “machine learning” type top-down approach, the “insights have come when scientists systematically analyse the components of processes that are easily manipulated in the laboratory — largely in model organisms. They’re still using a systems approach, but focusing it through a more traditional, bottom–up lens.” Through this systemic bottom-up approach, Davidson's group has deciphered the mechanism of how gene expressions are controlled through regulatory interactions and specify the construction of sea-urchin’s skeleton.

Interestingly, Eric Davidson gave a seminar at C2B2 on Thursday, April 22, titled "Causal Systems Biology: the Sea Urchin Embryo Gene Regulatory Network". I was very impressed by his talk, and the points he made in the Nature News article.

3DNA in the June 2010 issue of JBSD on "Current Perspectives on Nucleosome Positioning"

2010-04-16T22:16:00.007-04:00

While updating 3DNA citations this week, I noticed five of them are from the same June 2010 issue of JOURNAL OF BIOMOLECULAR STRUCTURE & DYNAMICS (JBSD), which is focused on "Current Perspectives on Nucleosome Positioning". Most of the papers in this issue are from well-known laboratories in computational structural biology. It is my pleasure to see 3DNA being widely used in the important research area of nucleosome positioning.

I browsed through the abstracts of all the papers in this JBSD issue to refresh my knowledge of this field. While DNA sequence surely plays some role in nucleosome positioning, I remain to be convinced of the existence of a nucleosome "code" (yet) in the sense of the generally applicable "genetic code". Overall, DNA is so flexible and the signal is so week, thus allowing for tailored data fitting to specific analysis, which is not transferable to other situations. Clearly, the area is hot, yet still wide open.

NSMB editorial: "Making your point-by-point"

2010-04-08T20:36:00.008-04:00

In the April 2010 issue of Nature Structural & Molecular Biology [NSMB, 17(4)], there is another interesting editorial, titled "Making your point-by-point". This editorial addresses an important issue in the process of publishing papers in peer-reviewed journals, that is: how to make effective point-by-point response to "those ever-demanding editors and reviewers"?

Overall, it can be helpful to put yourself in the reviewer’s shoes and compose a response s/he would find appropriate, where the concerns raised are considered and fully addressed. In its ideal state, the review process is a positive and constructive back and forth, an intellectual discussion in which the manuscript is the ultimate beneficiary.

Here is my re-cap of the main points, as I understand it. I am also taking this opportunity to read this one-page editorial one more time.

What to do?

Keep to the point – "makes a series of [succinct] points in response [directly] to each point raised by the reviewers."
Keep it objective – be diplomatic in your point-by-point response to the reviewers, "even if the reviewer’s wording might have seemed overly strong." You could be forthright in your cover letter to the editors, though.
Keep things under control – "Know when to go to the bench and when to argue."
The scope of things – "Say clearly and succinctly" when "some requests might genuinely be beyond the scope of the manuscript or might simply be unfeasible." "Try not to salami-slice", one strong and solid paper is (much) better than two weak ones!

Some don'ts, especially:

Mentioning celebrity endorsements. "you never know—they could be moonlighting as your most critical anonymous reviewer."
Trying to guess who the reviewers are when communicating to the editors – it does not help. Additionally, you could be plain wrong in your guess (again, you never know) – they are anonymous, literally.

Generally speaking, I think authors should be appreciative of the work of the reviewers and editors. Occasionally, I serve as a reviewer and I know the time and efforts it takes to make a fair and thorough assessment of a manuscript.

It is certainly not just because of politeness that in our 2008 3DNA Nature Protocols paper, we acknowledged:

We also thank the editor and the anonymous reviewers whose comments helped to clarify the presentation of the protocols.

More recently, in our 2010 NAR GpU paper, we acknowledged:

They also thank the anonymous reviewers, whose comments helped clarify the presentation of the manuscript.

What's special about the GpU dinucleotide platform?

2010-04-02T20:29:00.004-04:00

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

2010-03-26T23:45:00.010-04:00

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:

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.

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.

One computer, three operating systems

2010-03-20T17:27:00.005-04:00

While so far I have been quite happy with my new MacBook Pro, running Mac OS X 10.6 (Snow Leopard), I still feel more comfortable with the Ubuntu Linux programming environment I have been using for the past few years. Moreover, to make sure that my software (e.g., 3DNA) is strictly ANSI C compliant, and compiles without changes on the most commonly used operating systems (OSes), I need to have direct access to Linux and Windows. Luckily, the Intel-based hardware architecture of MacBook Pro and the free VirtualBox software make it possible to have the three OSes – Mac OS X, Ubuntu Linux, and Windows – in one computer.

Installing VirtualBox on Mac OS X was a snap. Specifically, I added the following two guest OSes:

Windows XP, with 1 GB RAM and 70 GB (virtual) hard disk
Ubuntu 9.10, with 2 GB RAM and 90 GM disk space

For seamless integration between each of the two guest OSes and the host Mac OS X, and for improved performance, I also created shared folders and installed guest additions for Windows and Linux. For Windows XP, the process had been quite straight forward. For Linux guest addition, however, I had some problems and solved them by following the instructions on "How To Install VirtualBox Guest Additions in Linux".

Now in Fullscreen Mode (command-F), I can run Ubuntu Linux or Windows XP as if it is native for each. Very cool!

Hoogsteen base-pair

2010-03-13T11:27:00.012-05:00

The A·U (or A·T) Hoogsteen pair is a well-known base pair (bp), named after the scientist who discovered it. As shown in the Figure below (left), in the Hoogsteen bp scheme, adenine uses its N7 and N6 atoms (at the major groove edge) to form two H-bonds with the N3 and O4 atoms from uracil, respectively. Interestingly, if the uracil base ring is flipped around the N7(A)…N3(U) H-bond by 180 degrees, N6(A) can also form an H-bond with O2(U), i.e., N6(A)…O2(U): this pairing scheme is called the reverse Hoogsteen bp (right).

I first came to know about the Hoogsteen bp from Saenger's book ("Principles of Nucleic Acid Structure"). Over the years, I have read many articles mentioning the Hoogsteen bp and touched this topic myself in the 2003 3DNA NAR publication. However, I have never read Hoogsteen's two original publications on this topic until recently:

The two-page long preliminary report, titled "The structure of crystals containing a hydrogen-bonded complex of 1-methylthymine and 9-methyladenine", was published in Acta Cryst. (1959). 12, 822-3. The paper contained only a single reference to the Watson-Crick DNA structure paper, published in Nature in 1953. I found it very revealing to understand why Hoogsteen used the methyl-ed derivatives of thymine and adenine, and how the failed initial interpretation of the experimental "vector-density map" using the Watson-Crick A-T bp led to the discovery of the new base-pairing scheme:
The fact that the first trial structure could not be refined led to a more critical scrutiny of the generalized projection and a greater emphasis on the significance of certain spurious peaks and on relatively large variations in the heights of peaks that were assumed to represent atoms. The correct structure was finally discovered by changing the positions of a few atoms in the 9-methyladenine portion of the asymmetric unit.
The more extensive account of the Hoogsteen bp story, titled "The Crystal and Molecular Structure of a Hydrogen-Bonded Complex Between 1-Methylthymine and 9-Methyladenine", published in Acta Cryst. (1963) 16, 907-16.

I like these two papers, and more generally those focused-articles, where authors get directly to a point and addressed it thoroughly and clearly. Most publications nowadays are very ambitious, trying to solve "big problems": the papers are generally far more complicated and often have "reproducibility" problems.

As a side note, the term Hoogsteen "edge" appears quite frequently in today's publications of RNA structures: in the Leontis-Westhof bp classification scheme, the term simply means the major groove edge in what would be a Watson-Crick bp geometry.

MacMost, a valuable resource to help you get the most of your Mac OS X

2010-03-06T18:16:00.005-05:00

Upon receiving my new Mac OS X Snow Leopard, I googled around, trying to find some tutorials on the web. Somehow, I came across a video clip by Gary Rosenzweig. I then visited MacMost.com and watched more video clips over there during the past week.

Overall, I like the videos quite a bit: ~5 minutes long each, these podcasts show various Mac-related tips and tricks in an easy to follow fashion. Specifically, I like the following:

#347: "Quick Look" – a functionality of looking at the contents of a file without opening it. "Quick Look" seems to be unique to Mac OS X since I am not familiar with it in Linux and Windows. So far, I have found it especially handy in Mail for quickly checking contents of attachments.
#357: "Do Macs Need Anti-Virus Software?" – it is assuring to know that "There are currently no active Mac viruses", and helpful to be aware that "anti-virus software could cause unexpected problems."
#363: "Learning to Program with Scratch" – it is from this video clip that I came to know the Scratch programming language from the MIT Media Lab. Unlike professional computer languages such as C, C++, Java, Ruby, Perl etc, Scratch targets the general public, especially for kids to learn mathematical and computational ideas by programing using a simple drag-and-drop interface. Using Scratch, it is really easy and cool create interactive stories and animations, and to share them on the web.

I will certainly keep visiting back MacMost.com and watch more videos as they become available. Little by little, I will learn new tricks to make my Mac life more enjoyable.

Mac OS X Snow Leopard -- I'm loving it (mostly)!

2010-02-26T21:03:00.004-05:00

Recently, when it was time for a new laptop, I decided to buy a MacBook Pro (Intel-based with Mac OS X 10.6.2 -- Snow Leopard). Over the past few days, I have been playing around with it, migrating files from my Ubuntu Linux box. So far, things have gone through smoothly, thus by and large, I am enjoying my new Mac.

Over the years, I have been using Ubuntu Linux and I have been very happy with it, especially for software development. Lyx and OpenOffice are handy for writing technical documents. However, I have realized that when it comes to write a manuscript for publication, and to communicate effectively with non-Linux collaborators, MS Word (with EndNote) is the standard. So I set up a Windows XP virtual machine via VirtualBox on my Ubuntu Linux box, which avoids the problem of dual booting and allows for easy file sharing between Linux and Windows.

Mac OS X is Unix/Linux based but has native support for MS Office and Adobe Suite of programs, so it seems an ideal choice for a new laptop. Mac OS X 10.6 (Snow Leopard) is claimed to be "The world's most advanced operating system. Finely tuned." Other things aside, I do appreciate the fact that 10.6 (Snow Leopard) is a refinement of 10.5 (Leopard) from installation to shutdown -- "In ways big and small, Mac OS X Snow Leopard makes your Mac faster, more reliable, and easier to use."

So far, I have configured Mail to access my Columbia emails. I must say that Mail is way better than Columbia's CubMail web-interface, and I like Mail's native integration with iCal and Address Book. Safari still needs some getting used to, from my mostly Firefox experience. However, it is nice to find that some websites, which does not work in Firefox but IE, display properly with Safari. Preview appears to be powerful for PDF and image viewing and manipulations. I have installed Xcode, and may explore it more, if nothing but to see what an IDE has to offer. Of course, it is nice to have direct access to MS Office (mostly for Word and PowerPoint, so no need to play around with OpenOffice), EndNote, Adobe (Acrobat, Photoshop, Illustrator), etc.

Some nuisances up to this point:

Keyboard missing numeric keypad and Home/End/PgUp/PgDn
Ctrl-C/V etc keyboard shortcuts I am used to now become Command-C/V etc
File and directory names are not case-sensitive -- most surprising!

Overall, my new MacBook Pro is a very nice toy to play. As I become more familiar with it, I may like it more, hopefully.

NSMB editorial: "Scientific writing 101"

2010-02-06T11:03:00.004-05:00

In the February 2010 issue of Nature Structural & Molecular Biology [NSMB, 17(2), p.139], there is a nice Editorial titled "Scientific writing 101". This short one-page essay is a good example of a (scientific) writing that is "a pleasures of reading".

"Less is more when it comes to writing a good scientific paper. Tell a story in clear, simple language and keep in mind the importance of the ‘big picture’."

Specifically, the editorial makes the following points:

Tell a story. A scientific paper is not a chronology; the data should be presented and interpreted in context.
Be clear. "Clear, simple language allows the data and their interpretation to come through."
Provide an informative title and abstract. "Make the abstract clear and try to get the ‘big picture’ across."
Make the introduction short and concise.
Clearly distinguish Results from Discussion. "Discussion should put those results in a broader context." It "should be an interpretation of those results..."
Cover letter is important. You should spell check your manuscript, and number the pages, etc.

In this blog post, I am just recapping the key points of the editorial, and taking the opportunity to re-read it. Following the simple principles outlined in the editorial would be beneficial to everyone in the scientific community.

Chinese new year greeting card

2010-01-30T19:27:00.004-05:00

Recently, I received the following nice new year greeting card from Dr. Pascal Auffinger (IBMC-CNRS, France). While never meeting each other in person, we have communicated via numerous emails: over the years, we have discussed extensively on 3DNA-related topics, and on nucleic acid structures in general. Pascal is one of my respected scientists, and his greeting card makes my support of 3DNA so gratifying.

Somehow, the picture reminds me Mountain Tai I climbed while at college, and the famous poem "会当凌绝顶，一览众山小".

According to the Chinese Zodiac, the Year of 2010 is the Year of Tiger (Feb. 14, 2010 to Feb. 2, 2011). Interestingly, this Chinese New Year's day coincides with Valentine's Day, and it is on Sunday.

Chemical diagram of Watson-Crick base-pairs

2010-01-23T21:55:00.006-05:00

Once in awhile, I need to refresh my memory about the chemical identities of the most common nucleobases: A, C, G, T/U. Sometimes, it is also necessary to explain to non-(bio)chemists about the concept of H-bond donor vs acceptor, and the major- vs minor-groove of the DNA double helix. In such cases, I use the following (customized) chemical diagram of Waton-Crick base-pairs (WC-bp):

Before taking the effort to create my own version of the Waton-Crick base-pair diagram, I googled around and found many illustrations (like the one in wikipedia). However, none of them suits my needs perfectly:

Trained as a chemist, I would like to see chemical bond types (double vs. single bond);
Working extensively with PDB format, I want to have the atom numbering information as well.

So I ended up to (re)create my own version of the WC-bp diagram: I used Chemtool to sketch the framework, and Xfig to fine-tune it. Overall, the diagram serves my purpose quite well, and hopefully others would find it useful as well.

Requests for SCHNAaP/SCHNArP source code

2010-01-07T21:38:00.003-05:00

Recently, I received several requests for the source code of the SCHNAaP/SCHNArP, a software package for the analysis and rebuilding of double helical nucleic acid structures. This suite of programs was developed ten years ago during my PhD work on DNA base-stacking interactions with Dr. Chris Hunter at the University of Sheffield, England.

Users become interested in SCHNAaP/SCHNArP mostly because of 3DNA, which can be taken as its superseded, more popular version. Due to Rutgers' policy of not releasing the source code of 3DNA, users who would like to understand details of the underlying algorithms thus turn back to SCHNAaP/SCHNArP. The interface is a bit aged, but the mathematics is still valid: it could serves well as a start point for those who really want to get into the world of nucleic acid structures.

Overall, though, I have a mixed feeling in this situation. On one hand, I am happy to see people becoming interested in my (previous) work. On the other hand, however, it also becomes clear that Rutgers' current licensing policy has blocked 3DNA's further circulation and adoption by the scientific community. Given the current trend of open-source software development, I see no reason to continue keeping 3DNA closed source. Making 3DNA open source (under a proper license term, of course) would allow for interested users to get more directly involved in the project, and thus to move the software to the next level.

ORCID -- an international research identification system?

2009-12-19T10:42:00.003-05:00

From the Nature news article titled "Credit where credit is due" in (462:7275, p. 825 on December 17, 2009), I came cross the ORCID initiative:

Name ambiguity and attribution are persistent, critical problems imbedded in the scholarly research ecosystem. The ORCID Initiative represents a community effort to establish an open, independent registry that is adopted and embraced as the industry’s de facto standard. Our mission is to resolve the systemic name ambiguity, by means of assigning unique identifiers linkable to an individual's research output, to enhance the scientific discovery process and improve the efficiency of funding and collaboration.

Overall, I think it is a good idea. If properly implemented and widely adopted, ORCID could help solve lots of issues associated with various ways of spelling a person's name due to, e.g., cultural differences. For example, put Chinese way, one's family name comes before one's given name, just the opposite of the western convention. Additionally, when a given name has two characters (quite common), there are could be a space or a hyphen (as I normally put in Xiang-Jun) or nothing in between. Combined with possible first name initials, there are already many ways to spell out a Chinese name.

The above Nature article, “Credit Where Credit is Due”, helps introduce the ORCID Initiative. As an specific example, it points to another article on page 843, where Nature profiles a research group trying to "complete the reference human genome sequence, which is still full of errors nearly a decade after the first draft was announced in 2000." Nature acknowledges that "It is essential work", "But it is also work that offers few academic rewards beyond the satisfaction of a job well done — it is unlikely to result in a high profile publication." Hopefully, by adopting the ORCID system, contributions of such types (e.g., software support and maintenance) would be more properly acknowledged by the scientific community.

Given the high profile of the founding parties, I am hopeful that the ORCID initiative would move forward as promised. I will keep an eye on it and see how it evolves.

Ribosomal structure: it helps to know some background information

2009-12-18T19:07:00.003-05:00

This year's Nobel Prize in Chemistry has been awarded to Venki Ramakrishnan, Tom Steitz, and Ada Yonath "for studies of the structure and function of the ribosome."

My connection with ribosomal structure began with the 50S large subunit of Haloarcula marismortui solved by Tom Steitz's group. Ever since the fully refined crystal structure at 2.4 Å resolution was published in 2001 (PDB entry: 1jj2; NDB code: rr0033), I have been using it to check 3DNA's applicability. In the two 3DNA papers (2003 NAR and 2008 NP), 1jj2 was used as an example to illustrate how find_pair can identify higher-order base-associations in complicated RNA containing structures. At the time, though, my understanding of the ribosomal RNA structure was purely geometrical: for quite a while, I got overwhelmed by the various biological terminologies, including the various S-es: 50S large ribosomal subunit vs. the 23S and 5S rRNA; and of course, the 30S small subunit vs. 16S rRNA.

Over the past year or so, I have become more interested in RNA structures. After reading a lot of related articles, gradually I feel things are becoming clearer than before. Nevertheless, there is something still missing, since my focus has (mostly) been on recent X-ray crystal structure-related work. My understanding of the ribosomal structure was finally put into context, thanks to following two recent publications:

One in Cell by James Williamson, titled "The Ribosome at Atomic Resolution".
Another one in Mol. Cell (in parallel and at the same time) by Joseph Puglisi, titled "Resolving the Elegant Architecture of the Ribosome".

These two papers not only summarized the significance of work of the three Nobel laureates — "the atomic resolution structures of the ribosomal subunits provide an extraordinary context for understanding one of the most fundamental aspects of cellular function: protein synthesis" — but also provided background information of decades of work from other players, including Harry Noller, Peter Moore, and Joachim Frank. Solving the ribosomal structure serves as a good example of how the fact that scientific research is both cooperative and competitive in nature.

Not all PDB entries are reliable; some could be plain fake

2009-12-11T21:28:00.012-05:00

With interest, I have browsed the recent thread in the PDB mailing list (pdb-l), "Retraction of 12 Structures" posted by Michael Sadowski and followed-up by Kevin Karplus et al. The story is about Krishna Murthy, a former scientist at the University of Alabama at Birmingham (UAB), who has been alleged to fabricate protein structures and published papers on them. Here is an informative comment by firebug36 from the above link:

I am a protein crystallographer myself, so just trust me - the results this gentleman [Murthy] published were falsified, and not in a smart way. The structures [for C3b] deposited in the Protein Data Bank made no physical sense.
Allegations against UAB group were first brought to light by several prominent people in the field, and not UAB officials:

http://www.nature.com/nature/journal/v448/n7154/full/nature06102.html

Accordingly to the post of Kevin Karplus, "several of the PDB files by Krishna Murthy's group were identified as problematic in the RosettaHoles paper". Naturally, then, comes the question, "should we remove ALL the PDB files from Krishna Murthy's group as suspect?"

The way Murthy's case coming to spotlight may represent an exception rather than norm. Imagine the scenario that he did not publish his C3b structure in Nature which caught the attention from leading crystallographers (Bert Janssen1, Randy Read2, Axel Brünger and Piet Gros), maybe Murthy is still publishing on protein structures today. In a sense, it is a hard to believe how Murthy could falsify 12 protein structures and published 9 papers in prestigious journals (including Nature, Cell, PNAS, JMB, Biochemistry, JBC etc) which have been cited 449 times.

PDB contains the state-of-the-art experimental data of bio-macromolecular structures. Yet, the archive is certainly full of inconsistencies/errors of various types. It would be helpful to know how many PDB entries are largely or partially wrong, and which can be taken as "gold standard" as far data quality is concerned.

This case gives an excellent lesson for those performing data-mining on macromolecular structures. Nowadays, PDB structures are many and keep increasing rapidly, but they are clearly of varying quality. Structural bioinformatics is about solving biology problems using informatics tools. Thus knowing the caveats of your data (how reliable are they?) and tools (what are their limitations?) is a prerequisite to draw sound scientific conclusions.

3DNA in the PCCP nucleic acid simulations themed issue

2009-12-06T17:21:00.005-05:00

While checking 3DNA-related citations through Web of Science for this past week, I found a total of nine times, as follows:

Five times to the 3DNA 2003 NAR paper
Once to the 3DNA 2008 NP paper
Three times to the 2001 standard base reference frame paper

Most interestingly, all the citations are from the same nucleic acid simulations themed issue of Physical Chemistry Chemical Physics 11 (45). Honestly, I was quite a bit (nicely) surprised by the fact, so I browsed the articles online. Edited by Charles Laughton and Modesto Orozco, the 2009 PCCP "themed issue exemplifies the rich diversity of cutting-edge research in the field of nucleic acids simulation." Indeed, quite a few well-known experts are among the authors of the two perspectives and 16 papers.

While not an "energetic" person myself, over the years I have been keeping an eye on MD simulations and MM calculations of nucleic acid structures. It is my pleasure to see that the 3DNA is being widely used (certainly more than I originally expected) by the nucleic acid simulations community. Given time, and with a suitable collaborator, I am open to consider adapting 3DNA to currently available MD simulation packages to make life easier for practitioners in this "dynamic" field.

See the effect of C preprocessing with the -E option of gcc

2009-12-05T23:18:00.009-05:00

Recently, I was interested in understanding better of one part of a C program, which is very generic, covering a lot of grounds with preprocessing options (#if ... #endif). However, I would like to see the minimum that covers the section I cared about. I vaguely remembered there is an option in the gcc compiler, from reading the book "An Introduction to GCC" a while ago, that can stop the process right after the preprocessing step (just like -c option stops after the compilation step). A quick check of "man gcc" revealed that it is -E:

-E Stop after the preprocessing stage; do not run the compiler proper. The output is in the form of preprocessed source code, which is sent to the standard output.

Input files which don't require preprocessing are ignored.

After setting the proper macros and running gcc with -E, I could immediately focus on the component to get what I wanted.

I have been using gcc for more than ten years, still there are more handy tricks to uncover.

Here is a simplified example showing the effect of the -E option. The following C function is saved in file "check_gcc_E.c".

--------------------------------------------------------------------------
#define GREETING "Hello World"

void check_gcc_E(void) {
#ifdef VERBOSE
   printf("Macro VERBOSE is defined.\n");
   printf("So you see: '%s'\n", GREETING);
#endif
   printf("Hello everyone!\n");
}
--------------------------------------------------------------------------

gcc -E -DVERBOSE check_gcc_E.c

--------------------------------------------------------------------------
# 1 "check_gcc_E.c"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "check_gcc_E.c"

void check_gcc_E(void) {

   printf("Macro VERBOSE defined\n");
   printf("So you see: '%s'\n", "Hello World");

   printf("Hello everyone!\n");
}
--------------------------------------------------------------------------

gcc -E check_gcc_E.c

--------------------------------------------------------------------------
# 1 "check_gcc_E.c"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "check_gcc_E.c"

void check_gcc_E(void) {
   printf("Hello everyone!\n");
}
--------------------------------------------------------------------------

Registered COPPA Users in phpBB3

2009-11-20T22:39:00.004-05:00

Recently, I received an email from a 3DNA user who registered at the 3DNA forum, but could not see anything at all once logged in. 3DNA forum is based on phpBB3, and has been running for over three years now. So at the very beginning, I thought how could it be? I had never heard of any such problem/complain from 3DNA forum registers before. I even created a temporary test login account and found no problem. So I communicated with the user and asked her/him to log in using my test account, and again everything was fine!

To reproduce the problem, I logged in as the user, and found one thing spurious: the user was in the group of "Registered COPPA Users", not the normal "Registered Users". I did not know what that COPPA stands for. So I googled the phase "Registered COPPA users" and the top hit led me into the phpBB3 document on Group Management, and the section I am interested in reads as follows:

Registered COPPA users are basically the same as registered users, except that they fall under the COPPA, or Child Online Privacy Protection Act, law, meaning that they are under the age of 13 in the U.S.A. Managing the permissions this usergroup has is important in protecting these users. COPPA doesn't apply to users living outside of the U.S.A. and can be disabled altogether.

So a registered COPPA user is, by definition, under the age of 13. By default phpBB3 does not even allow such a child to read any content! In the context of 3DNA forum, this policy simply does not make any sense — the contents (in the public section) are viewable by any one without registration.

It turned out that at registration stage, the first question is: "To continue with the registration procedure please tell us when you were born." Two dynamically generated dates are given, one is "Before" a date defining an age over 13, and the other "On or after" it for below 13. Obviously the 3DNA user mentioned above clicked the wrong button.

After knowing where the problem was and how it was created, fixing it was straightforward. Interestingly, when I then checked the 3DNA forum registered users, I found five of them were in the "Registered COPPA Users" group. Obviously, the previous (wrong) registers did not complain — possibly lost interest in pursuing further, so this issue did not surfaced until recently.

In a real world as we live, what seems simple may not be. Nothing should be taken for granted.

3DNA in PDB

2009-11-20T20:36:00.005-05:00

As mentioned previously, PDB makes use of blocview (part of 3DNA) to generate the simple yet effective images for nucleic-acid-containing structures. That's the connection I knew of between 3DNA and PDB. By pure chance, however, I recently noticed the 3DNA entry in PDB — it is actually a protein structure, completely unrelated to the 3DNA software package!

Just out of curiosity, I browsed the abstract of the Liu et al. article, titled "Halogenated benzenes bound within a non-polar cavity in T4 lysozyme provide examples of I...S and I...Se halogen-bonding" [J Mol Biol. 2009 Jan 16;385(2):595-605]. I then downloaded the full PDF version of the paper and read it carefully through. This work studied binding interactions of benzenes with the internal cavity of L99A mutated T4 lysozyme. The authors demonstrated that the center of the phenyl ring can be shifted by more than one angstrom due to different halogen-substitutions (where the 3DNA entry corresponds to C6H5I), and (further) proved the concept that "the protein is flexible and adapts to the size and shape of the ligand". At better than 2.0 Å resolution, they also observed the I...S and I...Se halogen-bonds.

I became interested in this paper not just because of the name of 3DNA, for which a quick browsing over the abstract would be sufficient. This paper also reminded me of an early article I published with the title "Influence of fluorine on aromatic interactions":

Non-covalent interactions between aromatic ligands influence the conformations of metal complexes, and the system [M(OAr)₂L₂] has been used to investigate the difference between phenyl–phenyl, phenyl–pentafluorophenyl and pentafluorophenyl–pentafluorophenyl interactions. X-Ray crystal structures show that pentafluorophenyl groups adopt partially stacked orientations with the two aromatic rings close to parallel and with significant π overlap. In contrast, phenyl groups are skewed away from each other with only edge-to-face contacts. Phenyl–pentafluorophenyl interactions adopt a coplanar fully stacked geometry. These results have been rationalised on the basis of energy calculations (carried out blind) using a variety of empirical models for treating weak non-covalent interactions. The major cause of the different behaviour of the three systems lies in the electrostatic interactions between the π systems.

Knowing of the pattern of a PDB id, — 4 characters long: the first character is a numeral in the range 0-9, while the rest can be either numerals or letters — I played around with some other possible ids with my name initials in it. Indeed I found one, 1XJL, a protein structure of human annexin A2 in the presence of calcium ions. If you are bimolecular structure-oriented, why not have a try with some ids of special meaning to you — you might be related to PDB in some unexpected way!

How shear affects twist angle of a dinucleotide step?

2009-11-14T22:20:00.009-05:00

A recent post in the 3DNA forum, titled "NUPARM vs X3DNA twist values", made me to rethink the issue of how or why shear affects twist angle of a dinucleotide step.

To me, this problem has long been solved as demonstrated by the following two well-cited publications:

The Tsukuba report, a.k.a., "A Standard Reference Frame for the Description of Nucleic Acid Base-pair Geometry". When Dr. Olson and I were drafting this report, I felt clearly the need to caution the community of the intrinsic correlations between base-pair parameters and the associated step parameters (Figure 3 there) to avoid possible mis-interpretations in structural analysis. This is specially the case for the effect of shear on twist, since the G–U wobble base-pair is common in RNA and it has a ~2.0 A shear.

The 3DNA 2003 NAR paper. There is a subsection on the "Treatment of non-Watson–Crick base pairing motifs", and Figure 3 addressed specially on the issue:
"Large Shear of the G–U wobble base pair influences the calculated but not the ‘observed’ Twist. The 3DNA numerical values of Twist, 20° (top) and 43° (bottom), differ from the visualization of nearly equivalent Twist suggested by the angle between successive C1'···C1' vectors (finely dotted lines)."

It was thus a bit surprising that such question still popping up. On second thought, however, it is quite understandable: one cannot expect everyone to read that two papers; not to mention remembering such details. So I am glad that this question was brought up to my attention, and it made me thinking possible ways to document more thoroughly the many 3DNA-related "technical details" that are crucial for better understanding of nucleic acid structures.

Coming back to the shear on twist angle issue, the figure at the left shows a G–U wobble pair example (top), and a simple rationale: the base-pair is approximately of 10Å-by-5Å (as defined in SCHNArP/3DNA), so a 2Å shift will lead to an angle:

atan2(2, 10) * 180 / pi = 11.3 degrees

(i.e., the red dotted line relative to the bottom horizontal line).

To a first order approximation, that is the difference between RC8–YC6 (or C1'–C1') vs. the base-centered mean y-axis of the pair for calculating twist angle. So whenever one has a G–U wobble pair next to a normal Watson-Crick pair, there would be ~11 degrees difference in "calculated" twist angle between the two approaches (NewHelix/CEHS/SCHNAaP/NUPARM vs 3DNA/Curve+). Moreover, when a G–U wobble is next to a U–G wobble pair, the difference would be doubled to ~23 degrees!

It is worth mentioning that the issue here (as in other similar cases) is not which number is "correct" or which is "wrong": a number is a number. It is its interpretation that matters, and it is here that "details" do count.

It's sad to hear that Warren DeLano, author of PyMOL, passed away

2009-11-08T16:35:00.003-05:00

From a couple of mailing lists, I heard the sad news that Warren DeLano, author of PyMOL, passed away on Tuesday morning, November 3rd. He was only 37!

I have never met Dr. DeLano personally, nor even I communicated with him by email, but I am very aware of PyMol, the de facto standard nowadays for molecular graphics. In writing 3DNA Nature Protocols paper, I dug more deeply into PyMol. I was impressed by its interactive interface to .r3d files (Raster3D) and the high quality ray-traced images it produced. So I came up with a Perl script (x3dna_r3d2png) to convert automatically from a 3DNA generated .r3d file to a PNG image through the PyMol engine.

Through his seminal contributions to PyMol, Dr. DeLano achieved something very few others in computational chemistry/biology can match: he successfully mobilized literately thousands of software programmers and ordinary users from multi-disciplines to join him to produce phenomenal pictures, each of which is worth a thousand words!

It was due to Dr. DeLano's vision that he made PyMol open source so the community now has the possibility/opportunity to continue support and further improve the software. At this stage, however, no one is likely to knows PyMol code to the depth Dr. DeLano did, not to mention the leadership and enthusiasm that he brought to the project. Whatever the case, the community undoubtedly would appreciate Dr. DeLano's valuable contributions.

Thanks, Dr. DeLano, for bring PyMol to the world!

How to calculate torsion angle?

2009-10-31T19:49:00.009-04:00

Given the x-, y-, and z-coordinates of four points in 3-dimensional (3D) space, how to calculate torsion angle? Overall, this is a well-solved problem in structural biology, and one can find detailed description of the algorithm in text books and on-line documents. The algorithm for calculating torsion angle is implementated in virtually every software package in structural chemistry or biology.

Basic as it is, however, in my experience, it is very important to have a detailed appreciation of how the method works in order to really get into the 3D world. Here is a worked example using Octave/Matlab of my simplified, geometry-based implementation of how to calculate torsion angle, including how to determine its sign. No theory or (complicated) mathematical formula, just a step-by-step illustration of how I solve this problem.

Coordinates of four points in variable abcd:

abcd = [ 21.350  31.325  22.681
22.409  31.286  21.483
22.840  29.751  21.498
23.543  29.175  22.594 ];

Two auxiliary functions: norm_vec() to normalize a vector; get_orth_norm_vec() to get the orthogonal component (normalized) of a vector with reference to another vector, which should have been normalized.
```
function ovec = norm_vec(vec)
  ovec = vec / norm(vec);
endfunction

function ovec = get_orth_norm_vec(vec, vref)
  temp = vec - vref * dot(vec, vref);
  ovec = norm_vec(temp);
endfunction
```

Get three vectors: b_c is the normalized vector b→c; b_a_orth is the orthogonal component (normalized) of vector b→a with reference to b→c; c_d_orth is similarly defined, as the orthogonal component (normalized) of vector c→d with reference to b→c.

b_c = norm_vec(abcd(3, :) - abcd(2, :))
  % [0.2703158  -0.9627257   0.0094077]
b_a_orth = get_orth_norm_vec(abcd(1, :) - abcd(2, :), b_c)
  % [-0.62126  -0.16696   0.76561]
c_d_orth = get_orth_norm_vec(abcd(4, :) - abcd(3, :), b_c)
  % [0.41330   0.12486   0.90199]

Now the torsion angle is the angle between the two vectors, b_a_orth and c_d_orth, and can be easily calculated by their dot product. The sign of the torsion angle is determined by the relative orientation of the cross product of the same two vectors with reference to the middle vector b→c. Here they are in opposite direction, thus the torsion angle is negative.
```
angle_deg = acos(dot(b_a_orth, c_d_orth)) * 180 / pi
  % 65.609
sign = dot(cross(b_a_orth, c_d_orth), b_c)
  % -0.91075
if (sign < 0)
  ang_deg = -angle_deg  % -65.609
endif
```

A related concept is the so-called dihedral angle, or more generally the angle between two planes. As long as the normal vectors to the two corresponding planes are defined, the angle between them is easy to work out.

Moreover, the method to calculate twist angles of helical nucleic acid structures in SCHNAaP and 3DNA is essentially the same.

Upgrade to Ubuntu 9.10

2009-10-30T22:04:00.002-04:00

I am now upgraded to Ubuntu 9.10 Karmic Koala, released On October 29 by Canonical Ltd. The system is now up and running, even though there are still some (minor) issues to be resolved. Overall, it was an exciting exploration, and Internet and Google search were essential for solving most of the problems.

Normally, I am not that quick to catch up with a new software release, but wait a few more weeks when most initial bugs have been fixed. I was quick this time mainly because I had been trapped into an earlier 9.10 alpha release (in development branch), when I tried to fix a printing problem (without success). Worse yet, suddenly, my VirtualBox to run Windows XP did not start, and SCIM Chinese input stopped functioning, etc. It caused me quite some time and trouble to get my Linux box back to work. So over the past few months, I did not perform any update, but eagerly waiting for the stable release. It is a relief that the new official 9.10 release does solve my problems.

With Ubuntu 9.10, now I have access to OpenOffice 3.1.1 and Lyx 1.6.4 for text processing, GNU Emacs 22.2.1, and GCC 4.4.1, among other things. It is nice to stay current, not only in science, but also in IT.

It is worth mentioning that Ubuntu is "an ethical concept of African origin emphasizing community, sharing and generosity." It is amazing how useful and robust the free, open-source software can be!

Sharing published data?

2009-10-23T23:45:00.007-04:00

A letter in the recent issue of Science [VOL 326 23 OCTOBER 2009] titled "The Antidote to Bias in Research" by Allison mentioned the Nature Genetics article on the issue of wide-spread non-repeatability of published microarray gene expression studies. It also touched on another recent study titled "Empirical Study of Data Sharing by Authors Publishing in PLoS Journals" by Savage and Vickers [PLoS ONE 4(9): e7078]:

In another study, refusal to share data despite policies requiring sharing was nearly ubiquitous among authors publishing in Public Library of Science journals.

I then read the Savage and Vickers article, and was surprised to find that only one author out of 10 sent them an original data set. Therefore, the authors wrote:

In conclusion, our findings suggest that explicit journal policies requiring data sharing do not lead to authors making their data sets available to independent investigators.

I can understand why the repeatability rate is so low from the Nature Genetics paper given the so many details associated with a published table or figure. However, it is a bit hard to explain why it is so difficult to share just the original dataset associated with a published work. I would imagine that the authors would be honored that others care about their publications, and contact them directly.

In my experience, I have also been simply ignored frequently for clarifications of some details in published work, or requests for some datasets. The highest successful rate is asking for PDF represents from corresponding authors.

If the basic principles of scientific publications are strictly enforced, even just by the big journals (so many, nowadays), a lot of unfounded big claims would be gone or can be easily seen through. Well, I know that's only possible in an ideal world.

How to quantify the relative geometry between two base-pairs? — Part 2

2009-10-23T21:07:00.007-04:00

In part I of this series, and indeed other occasions in my blog, 3DNA home page and forum, I have attributed to several related programs where 3DNA has benefited from. Before going into details on how the various parameters are calculated in 3DNA, however, I feel it is in order to make it clear my philosophy on creating scientific software in general, and 3DNA in particular.

Some underlying considerations

Whenever possible/applicable, I always try to get a thorough understanding of related software available. By thorough I mean at the source code level to see how an algorithm is implemented. The benefits of doing this are several folds:
1. Not to recreate the wheel, but to build upon previous work.
2. It is the most effective way to learn — many math formulae or text descriptions of algorithms, while useful for getting general ideas, are lack in details or vague in nature. In bioinformatics, most of the fundamental mathematics are not new. It is more about a novel combination of various known-parts, applying to a specific problem. Reading the source code is the only way to see unambiguously the implementation details.
3. With a clear understanding of how the parameters are actually calculated, one can make better use of a software tool, know its limitations and thus avoid misinterpretations.
4. To make objective, convincing (even to the original authors) comparisons/comments on existing tools, and importantly, to create something better if needed.

Specially, 3DNA has benefited the most from my SCHNAaP/SCHNArP programs: the underlying algorithms in 3DNA for calculating the various structural parameters (propeller, buckle, shear, roll, slide, x-displacement, inclination etc) are exactly the same as there — analyze and rebuild are direct derivatives from SCHNAaP (a for analysis) and SCHNArP (r for rebuilding/reconstruction), respectively; the ideas of standard stacking diagrams, the Zp parameter, building atomic models with sugar-phosphate backbone, and the base/bp rectangular block representations are also from SCHNAaP/SCHNArP. I also borrowed ideas from Babcock's RNA, Bansal's NUPARM, Dickerson's NewHelix/FreeHelix etc., which will be mentioned explicitly in following sessions.
Code wise, 3DNA was created from scratch in strict ANSI C (over 25K lines). In addition to implementing a unique combination of existing algorithms, I have added many significant novel methods including the find_pair program. However, feature-rich has never been my goal. Instead, I only consider to add new functions that I understand clearly and find useful. On the other hand, I am always quick to fix bugs.

Overall, it seems that I am an outlier rather than norm in scientific programming. The reasons could be as follows:

It is hard to understand other people's implementations of algorithms since it needs to read between the lines. This is especially the case with a unfamiliar computer language; undocumented, or bad code. It is easier to create one's own program than to understand and modify third-party software.
For publication purpose, it is more attractive to work on something "new" instead of building on others.

Thus, in bioinformatics field, there are far more "new" methods papers published than refinements of previous ones. However, what claimed/appeared to be novel may not be the case under the hook, due to different terminologies used in different fields. Moreover, while it is easy to create a self-claimed "new" program, others could simply do the same. In my experience, it is well worth the effort to really understand leading programs if one is serious about getting into an informatics field: it was really revealing and exciting when I went thorough CEHS, NewHelix/FreeHelix, RNA, Curves etc. Then I know clearly what was "the start-of-the-art" of the field and how I can do better: that's exactly how 3DNA came about!

In the following sessions, I will provide details on how 3DNA calculates the two sets of parameters to quantify the relative geometry of a dinucleotide step.

Duplicate tab function in Adobe READER 9.1 is very handy

2009-10-18T19:21:00.003-04:00

In reading a scientific paper, one often needs to jump back and forth for referred to tables, figures and references, etc. Those days, PDF is the standard way to share an e-document, and Acroread is the norm to view its content. I used to open two copies of the same document in two instances of Acroread (or one Acroread, one xpdf) so I can read continuously in one while moving around in the other, mostly at the end for references. This works, but obviously less than ideal.

One day, purely by chance, I (mouse) right-clicked the tab for the document I was viewing, and noticed that it popped up with two options: Detach Tab and Duplicate Tab. Clicking Duplicate Tab led to duplication of the same document in another tab. Actually, this process can be repeated more than once, thus allowing for multiple views of the same document simultaneously. Very neat!

Nowadays, whenever I read a scientific publication in Acroread, I often duplicate tab to have two views. It is only a click away to switch between the two. Thus when a citation is referred to in the main text, I can immediately see at the reference section what it is about.

A word of caution: I am using Ubuntu 9.04, with Acroread v9.1.2 05/25/2009. I have no idea from which version and on what platform such function was added.

How to quantify the relative geometry between two base-pairs? — Part 1

2009-10-18T12:10:00.008-04:00

In the field of double helical DNA (and RNA) structures, the following two sets of parameters to normally used to quantify the relative geometry of the two base pairs in a dinucleotide step:

shift, slide, rise, tilt, roll, and twist — which I call them stacking or simple step parameters
x-displacement, y-displacement, helical rise, inclination, tip, and helical twist — which I call them helical parameters

3DNA calculates all of these parameters. Over the years, I have been approached quite a few times on the local helical parameters interpretation question. In literature, I've noticed numerous times of confusions the community still has over the two set of parameters, especially with regard to the issue of "twist vs helical twist" and "rise vs. helical rise". This series of blog posts is aimed to clarify the problem by providing some background information and step-by-step worked examples illustrating how the parameters are actually calculated in 3DNA (which follows SCHNAaP/SCHNArP). In my experience, knowing how the parameters are derived is the key to understand what they mean, and thus to avoid misinterpretations.

Part 1 — background information

The Calladine and El-Hassan Scheme (CEHS) calculates only the six step parameters, i.e., shift, slide, rise, tilt, roll, and twist. The original CEHS implementation was a few FORTRAN subroutines taking advantage the code base of Dickerson's well-known NewHelix program. In developing SCHNAaP/SCHNArP, a collaborative project between the then Sheffield and Cambridge groups, I also went through the source code of NewHelix to get a thorough knowledge of its internals in order to integrate its nice features into SCHNAaP.
One thing I noticed was x-displacement, a parameter characterizing the "hole" of A-DNA in top view. The global helical axis from NewHelix was (still is) appealing to me, especially for a short, relatively straight fragment. Even for a clearly non-straight duplex, e.g., in super helical nucleosome core particle DNA, or a severely kinked DNA, the deviation from regular linear helix serves as a good parameter to quantify its overall non-linearity. As a side note, a single so-called bending angle is often misleading — for example, the bending angle as commonly cited in literature (mostly calculated from Curves), is strongly influenced by the two terminal base-pairs.
So I devised a parallel set of global helical parameters following the CEHS scheme of angle combination: here instead of roll-tilt, I used tip-inclination. The definition of global helical axis and the point the helix passes through are based on my simplified implementation in ANSI C of the algorithms used by NewHelix.
Thus SCHNAaP calculates both a set of local CEHS step parameters, and a set of global helical parameters — a unique combination of CEHS and NewHelix, made possible only after a thorough understanding of both methods. SCHNArP was developed to complete the circle, i.e., to rebuild a structure given a set of parameters, either the local step parameters: shift, slide, rise, tilt, roll, and twist, or the global helical parameters: x-displacement, y-displacement, helical rise, inclination, tip, and helical twist.
The RNA (Running Nucleic Acids) program by Babcock et al. calculates a set of local helical parameters. While working at Rutgers in Dr. Olson's group, I was interested in knowing how the local helical axis and the point its passes through were defined in the RNA program, since it would naturally substitute for the global one in SCHNAaP/SCHNArP to make two sets of purely local parameters. As it turned out, however, the algorithm was finally implemented in the 3DNA software package, as summarized in the email I sent out before the 13th Conversation at Albany:
The way to calculate the local helical axis and helical parameters in 3DNA is *essentially the same* as in RNA, but expressed in a much simpler way. First, the local helical axis is calculated as dx-times-dy, following Bansal, which gives the *same* result as the "single rotation axis" detailed in the RNA paper. I still could not figure out WHY, but verified this numerically. The location where the helix passes through is based directly on the RNA paper (i.e., following Chasles's theorem). The procedure to calculate helical parameters, i.e., tip/inclination/x-disp/y-disp/etc, following the SCHNAaP paper (i.e, with tip-inclination combination, which is consistent with the roll-tilt, propeller-buckle combinations in 3DNA). What amazed us is that this much simpler tip-inclination implementation in 3DNA gives *exactly* the same numerical values as the original RNA algorithm.
In summary, the two sets of local parameters as implemented in 3DNA are based on a unique combination of nice features from several programs well-known in the community, including SCHNAaP/SCHNArP (CEHS), NewHelix, RNA, and NUPARM. However, the 3DNA calculated parameters are numerically different from any of them, partially because 3DNA adopts the standard base reference frame. To the best of my knowledge, 3DNA is the only software that allows for rigorous conversion between the two sets of local parameters — a utility program in 3DNA called step_hel illustrates this simple fact. Again, as I wrote before the 13th Conversation at Albany:
Take each base-pair as a rigid block, only 6 parameters are required to relate one bp to the other. Two sets are commonly used: one set is "shift, slide, rise, tilt, roll and twist", and the other set is "x-displacement, y-displacement, helical rise, inclination, tip, and helical twist". Obviously, these two sets of parameters should be directly convertible, as demonstrated by Calladine and Drew in their 1984 JMB B-to-A transition paper, and illustrated in 3DNA manuscript.
Links to other programs referred to in this section:
- SCHNAaP/SCHNArP by Lu et al.
- The RNA program by Babcock et al.
- FreeHelix98 by Dickerson
- NUPARM from Banal et al.
- Standards, especially the email exchanges before the 13th Conversation at Albany

Blogger's duo-editing modes allow for flexibility and convenience

2009-10-10T21:06:00.007-04:00

In my experience using Blogger over the past several months, I have begun to appreciate more its double editing mode: Compose and Edit HTML. The Compose mode is a simple WYSIWYG editor, convenient for most common tasks. Once in a while, however, I get stuck with some nasty formatting issues that could drive one crazy to fix. This is where the Edit HTML mode comes in handy, which allows for full flexibility in editing raw HTML.

In principle, I am pretty competent with HTML and could use the Edit HTML mode directly. However, raw HTML is verbose and thus not that convenient. So I normally start a blog post with the Compose mode for most of the content, and switch to the Edit HTML mode only when necessary. Blogger makes the switch between the two modes a simple button click. This is in contrast to the single editing mode in phpBB3 (BBCode) used by the 3DNA forum, which does not allow for direct access to HTML.

Ideally, a software tool should be both flexible and convenient. In reality, however, not that many software could strike a balance between the two factors. Blogger is a nice example. Interesting, in composing this post, I switched back and forth between the two editing modes a couple of occasions: one after copy-and-pasting Edit HTML to stop red coloring of the following text, and the other to qualify the 3DNA forum link text.

Fiber models in 3DNA make it easy to build regular DNA helices

2009-10-09T20:47:00.008-04:00

3DNA contains 55 fiber models compiled from literature. To the best of my knowledge, this is the most comprehensive collection of its kind (detailed list given below), including:

Chandrasekaran & Arnott (from #1 to #43) — the most well-known set of fiber models
Alexeev et al. (#44-#45)
van Dam & Levitt (#46-#47)
Premilat & Albiser (#48-#55)

The utility program fiber puts the generation of all these fiber models in a simple, consistent interface, and can generate structures in either standard PDB or PDBML format. Of those models, some can be built with an arbitrary sequence of A, C, G and T (e.g., A-/B-/C-DNA from calf thymus), while others are of fixed sequences (e.g., Z-DNA with GC repeats). The sequence can be specified either from command-line or a plain text file.

Once 3DNA in properly installed, the command-line interface is the most versatile and convenient way to generate, e.g., a regular double strand DNA (mostly, B-DNA) of arbitrary sequence. Moreover, the w3DNA and 3DART web-interface to 3DNA makes it easy to generate a regular DNA model, especially for occasional use or educational purposes.

Theoretically, there is nothing to worth showing off in 3DNA's fiber model generation functionality. However, it serves as a clear example of the differences between a "proof of concept" and a practical software application. I initially decided to work on this issue simply for my own convenience. At that time, I had access to A-DNA and B-DNA fiber model generators, each as a separate program. Moreover, the constructed models did not comply to the PDB format in atom naming and other subtitles.

I started with the Chandrasekaran & Arnott fiber models which I had a copy of data files. Nevertheless, there are many details to work out, typos to correct, etc to put them in a consistent framework. For other models, I read each original publication, and typed into computer (with help from Dr. Olson's secretary) raw atomic cylindrical coordinates for each model. Again, quite a few inconsistencies popped up between the different publications with a time span over decades.

Overall, it was a quite tedious undertaking, requiring great attention to details. I am glad that I did that: I learned so much from the process, and more importantly, others can benefit from my effort. As I put in the 3DNA Nature Protocol paper (BOX 6 | FIBER-DIFFRACTION MODELS),

In preparing this set of fiber models, we have taken great care to ensure the accuracy and consistency of the models. For completeness and user verification, 3DNA includes, in addition to 3DNA-processed files, the original coordinates collected from the literature.

For those who want to understand what's going on under the hood, there is no better way than to try to reproduce the process using, e.g., fiber B-DNA as an example.

From the very beginning, I had expected the fiber functionality to be easily reachable and thus beneficial to all those who are interested in nucleic acid structures, especially to build a regular DNA duplex of chosen sequence. In my sense, the fiber program has not been widely used as it should have been, probably due to the fact that many people are not aware of its existence or capability. Hopefully, this blog post would help to get my message across.

PS. Given below is the content of the README file for fiber models in 3DNA


1. The repeating units of each fiber structure are mostly based on the
   work of Chandrasekaran & Arnott (from #1 to #43). More recent fiber
   models are based on Alexeev et al. (#44-#45), van Dam & Levitt (#46
   -#47) and Premilat & Albiser (#48-#55).

2. Clean up of each residue
   a. currently ignore hydrogen atoms [can be easily added]
   b. change ME/C7 group of thymine to C5M
   c. re-assign O3' atom to be attached with C3'
   d. change distance unit from nm to A [most of the entries]
   e. re-ordering atoms according to the NDB convention

3. Fix up of problem structures.
   a. str#8 has no N9 atom for guanine
   b. str#10 is not available from the disk, manually input
   c. str#14 C5M atom was named C5 for Thymine, resulting two C5 atoms
   d. str#17 has wrong assignment of O3' atom on Guanine
   e. str#33 has wrong C6 position in U3
   f. str#37 to #str41 were typed in manually following Arnott's
        new list as given in "Oxford Handbook of Nucleic Acid Structure"
        edited by S. Neidle (Oxford Press, 1999)
   g. str#38 coordinates for N6(A) and N3(T) are WRONG as given in the
        original literature
   h. str#39 and #40 have the same O3' coordinates for the 2nd strand

4. str#44 & 45 have fixed strand II residues (T)

5. str#46 & 47 have +z-axis upwards (based on BI.pdb & BII.pdb)

6. str#48 to 55 have +z-axis upwards

List of 55 fiber structures

id#  Twist   Rise        Structure description
    (dgrees)  (A)
-------------------------------------------------------------------------------
 1   32.7   2.548  A-DNA  (calf thymus; generic sequence: A, C, G and T)
 2   65.5   5.095  A-DNA  poly d(ABr5U) : poly d(ABr5U)
 3    0.0  28.030  A-DNA  (calf thymus) poly d(A1T2C3G4G5A6A7T8G9G10T11) :
                                        poly d(A1C2C3A4T5T6C7C8G9A10T11)
 4   36.0   3.375  B-DNA  (calf thymus; generic sequence: A, C, G and T)
 5   72.0   6.720  B-DNA  poly d(CG) : poly d(CG)
 6  180.0  16.864  B-DNA  (calf thymus) poly d(C1C2C3C4C5) : poly d(G6G7G8G9G10)
 7   38.6   3.310  C-DNA  (calf thymus; generic sequence: A, C, G and T)
 8   40.0   3.312  C-DNA  poly d(GGT) : poly d(ACC)
 9  120.0   9.937  C-DNA  poly d(G1G2T3) : poly d(A4C5C6)
10   80.0   6.467  C-DNA  poly d(AG) : poly d(CT)
11   80.0   6.467  C-DNA  poly d(A1G2) : poly d(C3T4)
12   45.0   3.013  D-DNA  poly d(AAT) : poly d(ATT)
13   90.0   6.125  D-DNA  poly d(CI) : poly d(CI)
14  -90.0  18.500  D-DNA  poly d(A1T2A3T4A5T6) : poly d(A1T2A3T4A5T6)
15  -60.0   7.250  Z-DNA  poly d(GC) : poly d(GC)
16  -51.4   7.571  Z-DNA  poly d(As4T) : poly d(As4T)
17    0.0  10.200  L-DNA  (calf thymus) poly d(GC) : poly d(GC)
18   36.0   3.230  B'-DNA alpha poly d(A) : poly d(T) (H-DNA)
19   36.0   3.233  B'-DNA beta2 poly d(A) : poly d(T) (H-DNA  beta)
20   32.7   2.812  A-RNA  poly (A) : poly (U)
21   30.0   3.000  A'-RNA poly (I) : poly (C)
22   32.7   2.560  Hybrid poly (A) : poly d(T)
23   32.0   2.780  Hybrid poly d(G) : poly (C)
24   36.0   3.130  Hybrid poly d(I) : poly (C)
25   32.7   3.060  Hybrid poly d(A) : poly (U)
26   36.0   3.010  10-fold poly (X) : poly (X)
27   32.7   2.518  11-fold poly (X) : poly (X)
28   32.7   2.596  Poly (s2U) : poly (s2U) (symmetric base-pair)
29   32.7   2.596  Poly (s2U) : poly (s2U) (asymmetric base-pair)
30   32.7   3.160  Poly d(C) : poly d(I) : poly d(C)
31   30.0   3.260  Poly d(T) : poly d(A) : poly d(T)
32   32.7   3.040  Poly (U) : poly (A) : poly(U) (11-fold)
33   30.0   3.040  Poly (U) : poly (A) : poly(U) (12-fold)
34   30.0   3.290  Poly (I) : poly (A) : poly(I)
35   31.3   3.410  Poly (I) : poly (I) : poly(I) : poly(I)
36   60.0   3.155  Poly (C) or poly (mC) or poly (eC)
37   36.0   3.200  B'-DNA beta2  Poly d(A) : poly d(U)
38   36.0   3.240  B'-DNA beta1  Poly d(A) : poly d(T)
39   72.0   6.480  B'-DNA beta2  Poly d(AI) : poly d(CT)
40   72.0   6.460  B'-DNA beta1  Poly d(AI) : poly d(CT)
41  144.0  13.540  B'-DNA  Poly d(AATT) : poly d(AATT)
42   32.7   3.040  Poly(U) : poly d(A) : poly(U) [cf. #32]
43   36.0   3.200  Beta Poly d(A) : Poly d(U) [cf. #37]
44   36.0   3.233  Poly d(A) : poly d(T) (Ca salt)
45   36.0   3.233  Poly d(A) : poly d(T) (Na salt)
46   36.0   3.38   B-DNA (BI-type nucleotides; generic sequence: A, C, G and T)
47   40.0   3.32   C-DNA (BII-type nucleotides; generic sequence: A, C, G and T)
48   87.8   6.02   D(A)-DNA  ploy d(AT) : ploy d(AT) (right-handed)
49   60.0   7.20   S-DNA  ploy d(CG) : poly d(CG) (C_BG_A, right-handed)
50   60.0   7.20   S-DNA  ploy d(GC) : poly d(GC) (C_AG_B, right-handed)
51   31.6   3.22   B*-DNA  poly d(A) : poly d(T)
52   90.0   6.06   D(B)-DNA  poly d(AT) : poly d(AT) [cf. #48]
53  -38.7   3.29   C-DNA (generic sequence: A, C, G and T) (depreciated)
54   32.73  2.56   A-DNA (generic sequence: A, C, G and T) [cf. #1]
55   36.0   3.39   B-DNA (generic sequence: A, C, G and T) [cf. #4]
-------------------------------------------------------------------------------
List 1-41 based on Struther Arnott: ``Polynucleotide secondary structures:
     an historical perspective'', pp. 1-38 in ``Oxford Handbook of Nucleic
     Acid Structure'' edited by Stephen Neidle (Oxford Press, 1999).
     
     #42 and #43 are from Chandrasekaran & Arnott: "The Structures of DNA
     and RNA Helices in Oriented Fibers", pp 31-170 in "Landolt-Bornstein
     Numerical Data and Functional Relationships in Science and Technology"
     edited by W. Saenger (Springer-Verlag, 1990).

#44-#45 based on Alexeev et al., ``The structure of poly(dA) . poly(dT)
     as revealed by an X-ray fiber diffraction''. J. Biomol. Str. Dyn, 4,
     pp. 989-1011, 1987.

#46-#47 based on van Dam & Levitt, ``BII nucleotides in the B and C forms
     of natural-sequence polymeric DNA: a new model for the C form of DNA''.
     J. Mol. Biol., 304, pp. 541-561, 2000.

#48-#55 based on Premilat & Albiser, ``A new D-DNA form of poly(dA-dT) .
     poly(dA-dT): an A-DNA type structure with reversed Hoogsteen Pairing''.
     Eur. Biophys. J., 30, pp. 404-410, 2001 (and several other publications).

3DNA in molecular dynamics simulations

2009-10-04T22:13:00.008-04:00

While updating 3DNA citations last Friday, I came across the paper "Flexibility of Short-Strand RNA in Aqueous Solution as Revealed by Molecular Dynamics Simulation: Are A-RNA and A′-RNA Distinct Conformational Structures?" by Ouyang et al. [Aust. J. Chem. 2009, 62, 1054–1061]. Through molecular dynamics (MD) simulations over a 30 ns period, the authors found that "the identification of distinct A-RNA and A′-RNA structures ... may not be generally relevant in the context of RNA in the aqueous phase." Overall, the paper is nicely written.

I have never performed MD simulations in my research experience, so normally I only read abstracts of such publications just to get a general idea of the main conclusions. What had attracted my attention to this work was its simultaneous citations to three earlier papers:

Resolving the discrepancies among nucleic acid conformational analyses.

Overview of nucleic acid analysis programs.

A standard reference frame for the description of nucleic acid base-pair geometry.

This is quite unusual, since nowadays it is far more common to only cite 3DNA itself (mostly 2003 NAR and/or 2008 NP). So I decided to have a look of the whole paper. It turned out that the authors used the extra citations to justify their choice of using 3DNA instead of Curves to calculate the helical parameters, including "the three main structural descriptors commonly used to differentiate between the two forms of RNA – namely major groove width, inclination and the number of base pairs in a helical twist [turn]".

I communicated to the authors about the availability of Curves+. Specifically, using one case (413D, one of the three structures in their Table 1, "Comparison of different results by CURVES and 3DNA programs"), I've tried to illustrate the point that Curves+ and 3DNA now give directly comparable parameters.

Of course, I am glad to see 3DNA being applied to molecular dynamics simulations of nucleic acid structures. Hopefully, more such applications will show up in the future, and I am willing to offer my help in ways that make sense to me.

Whenever in doubt, check with the author

2009-10-03T22:15:00.004-04:00

Once in a while, I send emails to authors of papers I am interested in, sometimes simply to ask for PDF reprints, mostly to request for clarifications of points I cannot understand fully. Of course, the responses I have received vary significantly: some authors are responsive and are able to answer my questions concretely; while others respond less professionally; in no small percentage, I get no feedback at all. Whatever the case, though, sending querying emails is convenient, and the responses I get (even no response at all) are informative. Naturally, I would take more seriously the papers whose authors are responsive. On the other hand, in my memory, I have never ignored a reader's question of my publications.

Seeking clarification on a scientific software from author(s) or maintainer(s) is even more important due to inherent subtleties of (undocumented) details, as is common in (bio)informatics. In supporting 3DNA over the years, I've experienced quite a few cases where authors of some articles are misinformed in making judgment about 3DNA's functionality. In one case, I read in a paper claiming 3DNA cannot handle Hoogsteen base-pairs while Curves can. A few email exchanges with the corresponding author (who was very responsive and professional) turned out that an internally modified version of Curves was used. More recently, I found a paper claiming that find_pair from 3DNA failed to identify some base pairs in DNA-protein complexes where a new method succeeded. I asked for the missing list, and immediately noticed that simply relaxing some criteria recovered virtually all of the pairs. Thus, to make a convincing comparison of scientific software, it is crucial to check with the original authors to avoid misunderstandings. Serious scientific software developers and maintainers always welcome users' feedback. Why not ask for clarifications if one really wants to make a (strong) point in comparison? Of course, it is another story for unsupported software.

The Internet age has provided unprecedented convenience for scientific communication. It would be a pity not to take full advantage of it. One simple and important thing to do is: whenever in doubt, ask for clarification from corresponding author of a publication or maintainer of a software.

On reproducibility of scientific publications

2009-09-27T22:50:00.010-04:00

In the September 25, 2009 issue of Science (Vol. 325, pp.1622-3), I read with interest the letter from Osterweil et al. "Forecast for Reproducible Data: Partly Cloudy" and the response from Nelson. This exchange of views highlights the difficulty/importance for one research team to precisely reproduce results from another when elaborate computation is involved. As is well known, subtle differences in computer hardware and software, different versions of the same software, or even different options of the same version, could all play a role. Without specifying those details, it is virtually impossible to repeat a publication exactly.

This reminds me of a recent paper "Repeatability of published microarray gene expression analyses" by Ioannidis et al. [Nat Genet. 2009, 41(2):149-55]. In the abstract, the authors summarized their findings:

Here we evaluated the replication of data analyses in 18 articles on microarray-based gene expression profiling published in Nature Genetics in 2005-2006. One table or figure from each article was independently evaluated by two teams of analysts. We reproduced two analyses in principle and six partially or with some discrepancies; ten could not be reproduced. The main reason for failure to reproduce was data unavailability, and discrepancies were mostly due to incomplete data annotation or specification of data processing and analysis.

Specifically, please note that:

The authors are experts on microarray analysis, not occasional application software users.
These 18 articles surveyed were published in Nature Genetics, one of the top journals of its field.
Not a single analysis could be reproduced exactly: two were reproduced in principle, six only partially, and the other ten not at all.

Without being able to reproduce exactly the results from others, it is hard to build upon previous work and move forward. The various reasons for lack of reproducibility listed by Ioannidis et al. are certainly not limited to microarray analysis. As far as I can tell, they also exist in the fields of RNA structure analysis and predictions, energetics of protein-DNA interactions, quantum mechanics calculations, and molecular dynamics simulations etc.

In my experience and understanding, the methods section in journal articles is not, and should not aim to be, detailed enough for exact duplication by a qualified reader. Instead, most such reproducibility issues would be gone if journals require that authors provide raw data, detailed procedures used to process the data, software version and options used to generate the figures and tables reported in the publication. Such information could be made available in (journal or authors) websites. This is an effective way to solve the problem, especially for computational, informatics-related articles. Over the years, for papers I am the first author or I have made major contributions, I've always kept a folder for each article to include every detail (data files, scripts etc) so that the published tables and figures can be repeated precisely. This has turned out to be extremely helpful when I want to refer back to early publications, or when I was asked by readers for further details.

As noted by Osterweil et al., "repeatability, reproducibility, and transparency are the hallmarks of the scientific enterprise." To really achieve the goal, every scientist needs to pay more attention to details and be responsive. Do not be fool around by the impressive introduction or extensive discussions (which are important, of course) in a paper: to get the bottom of something, it is usually the details that count.

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

2009-09-26T22:53:00.008-04:00

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.

Double helix groove width parameters from 3DNA

2009-09-05T21:45:00.008-04:00

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

2009-09-05T19:46:00.009-04:00

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.