Codon outliers: genes whose codon usage departs from their own genome

Score every coding sequence of a genome against the genome's own codon usage and find the genes whose choice among synonymous codons that average does not explain, with the false discovery rate controlled across the genes.

What it doesCopy Anchor Link

Choose the file that holds the genome's coding sequences, set the genetic code the genes are read with, and run. Each gene is compared with the pooled codon usage of every gene of its genome that passes the length filter, itself included. Choose several files and each is its own genome, scored against its own usage and, if you name one, read with its own genetic code; the false discovery rate is then controlled once across every gene of every genome. Turn on pooling for a genome split across files by chromosome, and all the chosen files become one genome. The result lists the outlier genes, largest departure first, with the gene name and product from each record's header, and the whole table of genes can be read in the same order.

When you’d use thisCopy Anchor Link

• Finding the genes a genome treats differently. In fast-growing bacteria and yeast the most highly expressed genes use a narrower set of codons than the average gene, so they stand out.
• Asking whether expression is what separates them. The result reports where the ribosomal proteins rank: near the top when expression shapes codon usage, in the middle when it does not.
• Shortlisting genes of unusual composition for a closer look, such as genes in a region of different base composition or of possible foreign origin.

What this computesCopy Anchor Link

For every gene with more than the minimum number of codons:

• B: Karlin's codon bias, the summed difference between the gene's within-amino-acid codon frequencies and the genome's, weighted by each amino acid's share of the gene. 0 is identical usage
• MILC: Supek and Vlahoviček's measure independent of length and composition, a log-likelihood ratio per codon less its expected value, so it does not grow with gene length
• G statistic of the gene's codon choice within each amino acid against the genome's usage, with degrees of freedom the sum over the amino acids the gene uses of (synonymous codons − 1)
• p-value, the upper tail of the χ² distribution on those degrees of freedom
• q-value, adjusted across every tested gene of every genome in the run by Benjamini–Hochberg (independent or positively dependent tests) or Benjamini–Yekutieli (any dependence, less power)
• outlier calls: the genes whose q is at or below the false discovery rate you set
• rank of every gene by MILC within its own genome, and the median rank of the ribosomal proteins when the headers name some
• with several genomes, the outliers per genome and the genetic code each was read with
• the genes left out as too short, with their codon counts

Codons are counted in frame from each gene's first base. A triplet holding any letter other than A, C, G or T is skipped, a trailing partial codon is dropped, and stop codons form a class of their own. Any NCBI genetic code can be chosen, for the run and for any one file. An option counts every alternative start codon of the code as methionine throughout each gene, which is how the reference R implementation reads genes by default; it is off here, because a codon such as an internal TTG is read as its own amino acid.

Every result states, in words generated from the run, why an outlier is not evidence of horizontal transfer, how much longer the called genes are than the rest (the test grows with gene length, so MILC and B say how far a gene departs and q only whether it does), where the ribosomal proteins rank, which genes were left out as too short, which codons the alternative-start option counted as methionine, and, with several genomes, how many outliers each has under the one false discovery rate and which genetic code read each.

Worked exampleCopy Anchor Link

All 524 coding sequences of Mycoplasma genitalium G37 (RefSeq NC_000908.2), genetic code 4, more than 80 codons, Benjamini–Hochberg at q ≤ 0.05:

• 23 genes are left out as too short, and of the 501 scored, 120 are outliers.
• Adhesins lead: MG_RS02985 has MILC 1.046 and B 0.678, and three of the top five genes are adhesins, the P1 adhesin among them.
• The 43 ribosomal proteins rank at a median of 161 of 501: expression is not what separates this genome's genes.

The 89 genes of Saccharomyces cerevisiae chromosome I with more than 80 codons give 59 outliers. The pyruvate kinase CDC19 ranks first and the translation elongation factor EFB1 fourth, both highly expressed, while MILC across all 89 genes is not correlated with expression (Spearman ρ 0.02 against RNA-seq): the measure finds the extreme genes, not an expression scale.

Run both files together, M. genitalium at code 4 and yeast at code 1, and each gene keeps its own genome's reference, so every MILC and p-value is the same as alone. Only the correction changes: across all 590 genes, 193 are outliers, 138 of 501 in M. genitalium and 55 of 89 in yeast. Pooling the two files as one genome instead changes every p-value, because each gene is then measured against the average of two species.

ControlsCopy Anchor Link

List decides how much of the result the table shows: Outliers or Every gene.

LimitationsCopy Anchor Link

• An outlier is not a transferred gene. Codon bias and base composition are poor indicators of horizontal transfer; expression, strand and regional composition produce outliers too.
• The χ² distribution of the G statistic is a large-sample limit, so the p-value of a short gene with many sparsely used amino acids is the least reliable. The length filter is the guard; raise it to be stricter.
• On long genes the test is powerful enough to call small departures, so in a well-sampled genome many genes are called. Rank by MILC or B for size.
• One false discovery rate covers every genome in a run, so adding a genome changes the q-values of the others. Run a genome alone for its own α.
• The reference is the genome's own average. A genome whose genes split into two populations of usage (two replication strands, or two compositional regions) makes the smaller one look deviant.
• Genes are read in frame from their first base: a record that is not a coding sequence in frame gives a meaningless score.

Where the numbers come fromCopy Anchor Link

B and MILC are checked against coRdon, the Bioconductor implementation, on every gene of the two genomes above under both readings of alternative start codons, and the q-values against R's p.adjust, for each genome alone and for both genomes under one correction.

AutomationCopy Anchor Link

The same test is a workflow node and a tool the assistant can call, with the same settings and defaults as the panel; each genome is the name of a coding-sequence file, or a list of records.