P-value based visualization of codon usage data

The standard genetic code of protein coding DNA sequences shows a redundancy, since different triplet codons may be used to code for the same amino acid. In general, codon usages show organism-specific patterns. However, codon usage variation within a single genome can be an important source of information about gene expression levels and events of horizontal gene transfer. In particular, dimensionality reduction methods have widely been used for the analysis of codon usage patterns in microbial genomes. These methods provide a low-dimensional point representation of genes, where the proximity of gene-specific points indicates a similar codon usage of the associated genes. Hence, the resulting two-dimensional scatter plots enable a total view on the genome which may reveal a clustering of genes according to groups of nearby points. These clusters can for instance provide evidence for horizontal gene transfer according to groups of putative alien genes [1, 2] or for translational selection according to groups of highly expressed genes [3, 4].

As a standard method for scatter plot visualization of codon usage data, researchers mostly resort to the so-called correspondence analysis (CA) which has originally been developed for the analysis of contingency tables [5]. From the original formulation it is not completely clear how CA applies to codon counts. Because different preprocessing and normalization schemes have been proposed, the use of CA in codon usage studies has not been without controversy [6]. Nevertheless, CA has been applied for the analysis of many bacterial genomes, including those of Escherichia coli [1, 3], Bacillus subtilis [4, 7, 8], Borrelia burgdorferi [9, 10], Chlamydia trachomatis [11], Mycoplasma genitalium [12], Helicobacter pylori [13] and Pseudomonas aeruginosa [14].

Recently, self-organizing maps [15] have been proposed as an alternative visualization method for codon usage data [2, 16, 17]. Although this method provides a simultaneous clustering of the data which may be useful in certain contexts, it requires to choose the size of a discrete grid on which the genes are mapped in a non-linear way. The grid-size is a critical parameter of the method and directly controls the final clustering in the visualization. Unfortunately, the grid-size of self-organizing maps is a so-called hyperparameter which usually cannot be inferred from the data in an unsupervised manner. Therefore the resulting visualizations bare the risk of being highly subjective.

Here we present a visualization method, which has been tailored to the analysis of codon usage data while not depending on difficult to tune hyperparameters. Our visualization method is based on multidimensional scaling and a new similarity measure for codon usage data. In the following we first introduce our probabilistic similarity measure for codon usage tables and outline the corresponding algorithm for multidimensional scaling based on P-values. Then we provide some visualizations for the analysis of four microbial genomes and discuss our results in comparison with the results obtained from the classical correspondence analysis method.

For the analysis of codon usage tables we developed a special similarity measure which has been derived from the well-known chi-square test for the comparison of two distributions. Unlike the classical chi-square test we do not decide whether two distributions are equal or not, but instead we only use the corresponding P-values to compute a similarity measure for the underlying codon usage tables. For each pair of genes we compare the corresponding codon distributions on the basis of the codon frequencies in the two genes. For a suitable similarity score we average the P-values of the amino acid specific chi-square tests. We start with the counts

for codon of amino acid a _i in the j-th gene. These counts sum up to over the number L _i of different codons for amino acid a _i . Note that n _ij corresponds to the number of occurrences of amino acid a _i in gene j. With these counts we compute the chi-square statistic for each pair (j, k) of genes:

The classical chi-square test for comparison of two distributions is based on the following proposition: under the null hypothesis that the corresponding samples were drawn from the same probability distribution, the variable

is asymptotically chi-square distributed with L _i degrees of freedom. Here we do not perform a chi-square test, but rather calculate the P-value P _ijk associated with the chi-square statistic . The P-values are obtained from the chi-square probability function which is an incomplete gamma function [18]. A small value of P _ijk indicates a significant difference between the codon distributions of gene j and k with respect to amino acid a _i . For a number of M genes in a genome we then assemble the M × M matrix S of similarity scores with non-negative elements

where n _a is the number of amino acids. Note that S has unit diagonal elements, i.e. S _jj = 1, because the P-value for tables with identical counts is one. Consequently all off-diagonal elements are in the range [0, 1].

In order to derive a suitable low-dimensional point representation of genes we apply classical multidimensional scaling (see e.g. [19]) to the above similarities. The objective is to find a two-dimensional point configuration with interpoint distances reflecting the codon usage similarities of the corresponding genes. To perform classical scaling based on similarities we first transform the similarity matrix S into a positive semi-definite matrix C by subtracting the smallest eigenvalue λ_min of S from all of its diagonal elements:

C = S - λ_minI     (3)

where I is the M × M identity matrix. Note that this transformation preserves the equality of diagonal elements. With the M × M centering matrix H with elements

we finally obtain the matrix

B = HCH.     (5)

It can be shown that for a positive semi-definite matrix C the distance matrix D with elements obtained by the standard transformation is Euclidean and B is a centered inner product matrix ([19], pp. 402). Therefore principal components can be obtained from (partial) eigenvalue decomposition of B. Thus, for 2D-visualization we compute the two leading eigenvectors x₁ and x₂ of B associated with the largest and second largest eigenvalue, respectively. The M components of x₁ and x₂ provide the x₁ and x₂ coordinates for the M genes, which are utilized for scatter plot visualization.

We proposed an approach for the visualization of codon usage data, using multidimensional scaling (MDS). In that context we introduced a novel similarity measure for codon usage tables, which has been derived from the classical chi-square test. An important feature of our P-value based similarity measure is that it does not involve any hyperparameters. Therefore a subjective "bias" on the visualization due to user-adjusted parameters is effectively avoided. Our comparisons with the widely-used correspondence analysis (CA) method in most cases showed a slightly better clustering of putative alien genes for our P-value based visualization. In particular the results indicate that our approach is more robust than the CA-based visualization method. The outlier-sensitivity of CA becomes apparent in the plots for all species considered here and has already been observed in previous studies [9]. While in most cases the CA-based visualizations are still useful in terms of a suitable clustering of ribosomal protein and putative alien genes, for T. thermophilus that sensitivity results in an inappropriate plot which complicates interpretation.