How to Parse Codeml Output Without Losing Your Mind

If you’ve ever run codeml (part of the PAML package) you know the feeling: hours of computation, excitement for results, then… a giant text file filled with numbers, log-likelihoods, omega values, and site classes that look like they belong in the Matrix.

Parsing codeml output can be a nightmare. But with a systematic approach, you can turn that wall of text into meaningful insights about molecular evolution. Here’s a step-by-step guide to keep your sanity intact.

---

🧩 Step 1: Know What You’re Looking For

Codeml can run many models (site models, branch models, branch-site models). Depending on the run, the output may contain:

• Log-likelihood (lnL) values → used for likelihood ratio tests (LRTs).

• dN/dS ratios (ω) → telling you about selective pressure.

• Parameter estimates (kappa, frequencies, proportions).

• Bayes Empirical Bayes (BEB) results → list of sites under selection.

👉 Before parsing, ask yourself: Am I comparing models? Extracting site-level selection? Or just summarizing dN/dS?

---

📂 Step 2: Break Down the File Structure

Codeml outputs look messy, but they follow a pattern. For a site model run, for example, you’ll typically see:

1. Header with settings (seqfile, treefile, models used).

2. Log likelihood (lnL) line.

3. Estimated parameters (kappa, omega, proportions, etc.).

4. Codon frequencies.

5. BEB results (sites under selection with posterior probabilities).

A sample lnL line looks like this:

lnL(ntime: 46  np: 52):  -12345.6789 +0.0000

Here:

• lnL = log-likelihood.

• np = number of parameters.

• ntime = number of branch lengths estimated.

---

⚡ Step 3: Automate the Pain Away

Manually copying numbers is risky and slow. Instead, use Python or R scripts to parse automatically.

Example (Python snippet):

import re

def parse_codeml_output(filename):
    results = {}
    with open(filename) as f:
        content = f.read()

    # Extract log-likelihood
    lnL_match = re.search(r'lnL.*?:\s+(-?\d+\.\d+)', content)
    if lnL_match:
        results['lnL'] = float(lnL_match.group(1))

    # Extract omega
    omega_match = re.search(r'omega\s*=\s*([\d\.]+)', content)
    if omega_match:
        results['omega'] = float(omega_match.group(1))

    # Extract BEB sites
    beb_sites = re.findall(r'(\d+)\s+([A-Z])\s+([\d\.]+)\*?', content)
    results['BEB_sites'] = [(int(pos), aa, float(prob)) for pos, aa, prob in beb_sites]

    return results

# Example usage
parsed = parse_codeml_output("codeml_output.txt")
print(parsed)

This will give you a dictionary with lnL, omega, and BEB sites neatly extracted.

---

📊 Step 4: Compare Models

Most codeml analysis relies on likelihood ratio tests (LRTs).

• Compare lnL values from two models (e.g., M1a vs M2a, M7 vs M8).

• Compute 2ΔlnL = 2*(lnL_alt – lnL_null).

• Compare against a chi-square distribution with the correct degrees of freedom.

This step tells you if positive selection is statistically supported.

---

🎨 Step 5: Visualize the Results

Humans are visual creatures — plots help make sense of the numbers.

• Plot sitewise BEB probabilities (bar chart with sites above 0.95 marked).

• Map positively selected sites onto a protein structure (using PyMOL or Chimera).

• Compare dN/dS ratios across branches with a tree heatmap.

Visualization turns codeml’s intimidating numbers into interpretable biology.

---

🛠 Step 6: Troubleshoot Common Pitfalls

• No convergence? → Try different initial omega or kappa values in the control file.

• Weird parameter estimates (e.g., ω > 99)? → Usually means unstable alignment or tree.

• Too many sites under selection? → Check alignment quality (gaps, stop codons).

---

✅ Final Thoughts

Parsing codeml output doesn’t have to be mind-breaking. With the right script and workflow, you can go from raw text file → neat table → clear biological interpretation in minutes.

👉 My golden rule: Don’t trust a single number. Always compare models, validate alignments, and visualize results.

With a few automated parsing tools, you’ll spend less time wrestling with output files — and more time uncovering the evolutionary stories hidden in your data.

---