Skip to main content

Selection pressure between two genomes (dN/dS)

TL;DR: dN/dS is a per-ortholog-pair measurement, so it is a colour on a synteny track. This builds a human against rhesus macaque ortholog table with jcvi and measures dN and dS on every pair with kaks_from_pairs.py.

Prerequisites

jcvi builds C extensions and will not install against every python. If pip install jcvi fails compiling them, uv venv --python 3.12 followed by uv pip install jcvi biopython gets an interpreter it does build on.

Where the data comes from

The two assemblies are GCA_000001405.29 (human GRCh38) and GCA_003339765.3 (rhesus macaque Mmul_10), with gene models and coding sequence from Ensembl release 116.

What dN/dS says

A coding substitution is synonymous (the codon changes, the amino acid does not) or non-synonymous. Synonymous changes are nearly invisible to selection, so their rate dS approximates the mutation rate; non-synonymous changes are seen, so their rate dN carries what selection did. The ratio is read against 1: below it is purifying selection, where most genes sit, and above it takes positive selection to explain.

Producing the data

dS has to be large enough to estimate and small enough not to saturate. Rhesus macaque sits in that window against human; chimpanzee leaves a denominator near zero on most genes.

Orthologs

The end-to-end script turns each GFF3 into the BED the adapter reads, translates each CDS to a proteome keyed the same way, and runs jcvi:

diamond makedb --in rhesus.pep -d rhesus.pep
diamond blastp --threads 14 --query human.pep --db rhesus.pep \
  --out human.rhesus.last --max-target-seqs 20 --evalue 1e-10 --outfmt 6
python -m jcvi.compara.catalog ortholog --no_strip_names --dbtype prot \
  --align_soft diamond_blastp --no_dotplot human rhesus

Two silent traps:

  • the alignment file has to be query = the first species, subject = the second. Reversed, every id is looked up in the wrong BED and the run ends with A total of 0 anchor was found
  • Ensembl versions transcript ids in its FASTA and not in its GFF3 (ENST00000641515.7 against ENST00000641515), so nothing matches. The script strips the version, and kaks_from_pairs.py takes --strip-version

dN and dS

python3 kaks_from_pairs.py pairs.tsv both.cds.fa.gz \
  --key record --strip-version -o kaks.tsv

kaks_from_pairs.py aligns each pair as protein, back-translates to codons so that nothing shifts frame, and runs Nei-Gojobori.

Filtering paralogs and low-count pairs

True orthologs share one divergence time, so their dS values cluster. A pair an order of magnitude above the cluster is a paralog the aligner preferred, and --max-ds removes it.

The top of a table sorted by dN/dS is the pairs with almost nothing to divide by: HBA1 comes out over 2 off a single synonymous difference. --min-syn-subs is a floor on that count.

Every row also carries that count and a two-sided Fisher exact p, the test MEGA prescribes for small substitution counts. Both are attributeColumns, so clicking a link shows the evidence under its colour.

Loading the blocks table in JBrowse

The output is a pair table with the two rates after the two gene columns, which is the .blocks shape MCScanBlocksAdapter reads:

{
  "type": "SyntenyTrack",
  "trackId": "primate_orthologs",
  "name": "Human / rhesus orthologs (dN/dS)",
  "assemblyNames": ["human", "rhesus"],
  "adapter": {
    "type": "MCScanBlocksAdapter",
    "uri": "primate.blocks.gz",
    "blockAssemblies": ["human", "rhesus"],
    "bedLocations": [{ "uri": "human.bed.gz" }, { "uri": "rhesus.bed.gz" }],
    "assemblyNames": ["human", "rhesus"],
    "attributeColumns": ["dn", "ds", "syn_subs", "fisher_p"]
  }
}
jbrowse add-track-json '{
  "type": "SyntenyTrack",
  "trackId": "primate_orthologs",
  "name": "Human / rhesus orthologs (dN/dS)",
  "assemblyNames": ["human", "rhesus"],
  "adapter": {
    "type": "MCScanBlocksAdapter",
    "uri": "primate.blocks.gz",
    "blockAssemblies": ["human", "rhesus"],
    "bedLocations": [{ "uri": "human.bed.gz" }, { "uri": "rhesus.bed.gz" }],
    "assemblyNames": ["human", "rhesus"],
    "attributeColumns": ["dn", "ds", "syn_subs", "fisher_p"]
  }
}'

In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add track from pasted JSON, and paste:

{
  "type": "SyntenyTrack",
  "trackId": "primate_orthologs",
  "name": "Human / rhesus orthologs (dN/dS)",
  "assemblyNames": ["human", "rhesus"],
  "adapter": {
    "type": "MCScanBlocksAdapter",
    "uri": "primate.blocks.gz",
    "blockAssemblies": ["human", "rhesus"],
    "bedLocations": [{ "uri": "human.bed.gz" }, { "uri": "rhesus.bed.gz" }],
    "assemblyNames": ["human", "rhesus"],
    "attributeColumns": ["dn", "ds", "syn_subs", "fisher_p"]
  }
}

attributeColumns names the columns after the two gene columns, and each becomes a feature attribute in the detail panel. dn and ds drive Color by... → dN/dS, whose ramp has 1 at its middle and 2 at its top.

Two LinearSyntenyView properties matter for a view this sparse: alpha defaults to 0.2 for whole-genome views where ribbons overlap, and 0.95 shows the colour as it is; drawCurves separates stacked neighbours.

Reading the plot

Human against rhesus macaque across a collinear neighbourhood on human chromosome 12, each ribbon one ortholog pair coloured by dN/dS. Lysozyme (LYZ) is the one gene above the ramp's pivot; its neighbour YEATS4 is at the other end. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

The steps behind the figure above, run against your own data.

  1. Open your genomes: on the JBrowse Desktop start screen click Open new genome (or Show all available genomes to pick a hosted one), then File → Open genome... for each of the rest.This figure's rows, top to bottom: human, rhesus.
  2. Open the view: Add → Linear synteny view, and set the rows to those assemblies in the import form.
  3. Point the import form at your own file: click the arrow between two rows, choose New track, and paste a URL or pick a local file.A synteny track the session already has is quicker: pick it under **Quick start**, which fills in the rows itself.This figure uses the “primate_orthologs” track.
  4. Click Launch.
  5. Synteny view header → palette button → Show color legend (checked)The last item in the same menu the color modes are in.
  6. Synteny view header → palette button → dN/dSThe palette button's tooltip reads "Color by: ...".
  7. Synteny display settings (the sliders button in the view header) → Opacity → drag to 0.95Lower values let dense overlapping ribbons show through each other.
  8. Synteny display settings (the sliders button in the view header) → Curved lines (checked)
  9. Add your own track: open this view's track selector, click + and choose Add track, then paste a URL or choose a local file.Row 1 (human) — This figure uses the “human_genes” track.
  10. Track menu → Show... → Labels → Name onlyRow 1 (human).
  11. Track menu → Set feature height → CompactRow 1 (human).
  12. Type your region of interest into the location box and press Enter.Row 1 (human) — This figure is at 12:67,835,000-70,835,000.
  13. Add your own track: open this view's track selector, click + and choose Add track, then paste a URL or choose a local file.Row 2 (rhesus) — This figure uses the “rhesus_genes” track.
  14. Track menu → Show... → Labels → Name onlyRow 2 (rhesus).
  15. Track menu → Set feature height → CompactRow 2 (rhesus).
  16. Type your region of interest into the location box and press Enter.Row 2 (rhesus) — This figure is at 11:67,401,000-70,319,000.

Open this view in JBrowse Desktop ↗

Opens JBrowse Desktop (JBrowse Desktop 5.0+) at this view and saves it as a reopenable session — swap in your own files afterwards.

Nothing happens? Paste this link into Desktop's Open JBrowse Web link... (start screen, or File → Session):

https://jbrowse.org/code/jb2/main/?config=https%3A%2F%2Fjbrowse.org%2Fdemos%2Fprimate_selection%2Fconfig.json&session=spec-%7B%22views%22%3A%5B%7B%22type%22%3A%22LinearSyntenyView%22%2C%22showColorLegend%22%3Atrue%2C%22views%22%3A%5B%7B%22assembly%22%3A%22human%22%2C%22loc%22%3A%2212%3A67%2C835%2C000-70%2C835%2C000%22%2C%22tracks%22%3A%5B%7B%22trackId%22%3A%22human_genes%22%2C%22type%22%3A%22LinearBasicDisplay%22%2C%22showLabels%22%3A%22name%22%2C%22displayMode%22%3A%22compact%22%7D%5D%7D%2C%7B%22assembly%22%3A%22rhesus%22%2C%22loc%22%3A%2211%3A67%2C401%2C000-70%2C319%2C000%22%2C%22tracks%22%3A%5B%7B%22trackId%22%3A%22rhesus_genes%22%2C%22type%22%3A%22LinearBasicDisplay%22%2C%22showLabels%22%3A%22name%22%2C%22displayMode%22%3A%22compact%22%7D%5D%7D%5D%2C%22tracks%22%3A%5B%5B%22primate_orthologs%22%5D%5D%2C%22colorBy%22%3A%22dnds%22%2C%22alpha%22%3A0.95%2C%22drawCurves%22%3Atrue%7D%5D%7D&sessionName=JBrowse+docs%3A+selection_pressure%2Flysozyme

This session spec pastes after &session=spec- on any JBrowse Web instance.

{
  "views": [
    {
      "type": "LinearSyntenyView",
      "showColorLegend": true,
      "views": [
        {
          "assembly": "human",
          "loc": "12:67,835,000-70,835,000",
          "tracks": [
            {
              "trackId": "human_genes",
              "type": "LinearBasicDisplay",
              "showLabels": "name",
              "displayMode": "compact"
            }
          ]
        },
        {
          "assembly": "rhesus",
          "loc": "11:67,401,000-70,319,000",
          "tracks": [
            {
              "trackId": "rhesus_genes",
              "type": "LinearBasicDisplay",
              "showLabels": "name",
              "displayMode": "compact"
            }
          ]
        }
      ],
      "tracks": [
        [
          "primate_orthologs"
        ]
      ],
      "colorBy": "dnds",
      "alpha": 0.95,
      "drawCurves": true
    }
  ]
}

Loaded against config: https://jbrowse.org/demos/primate_selection/config.json

Hand this to a coding agent. It rebuilds the figure above headlessly, and the session file is then the thing to edit — swap an adapter uri for your own file and rerun.

cat > session.json <<'JSON'
{
  "views": [
    {
      "type": "LinearSyntenyView",
      "showColorLegend": true,
      "views": [
        {
          "assembly": "human",
          "loc": "12:67,835,000-70,835,000",
          "tracks": [
            {
              "trackId": "human_genes",
              "type": "LinearBasicDisplay",
              "showLabels": "name",
              "displayMode": "compact"
            }
          ]
        },
        {
          "assembly": "rhesus",
          "loc": "11:67,401,000-70,319,000",
          "tracks": [
            {
              "trackId": "rhesus_genes",
              "type": "LinearBasicDisplay",
              "showLabels": "name",
              "displayMode": "compact"
            }
          ]
        }
      ],
      "tracks": [
        [
          "primate_orthologs"
        ]
      ],
      "colorBy": "dnds",
      "alpha": 0.95,
      "drawCurves": true
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/demos/primate_selection/config.json \
  --session session.json -o figure.png

Using JBrowse with AI agents covers the rest: hosted genomes to point it at, and why waiting for the render is the part that goes wrong.

The neighbourhood is collinear, so colour is the only thing that varies. Adaptive evolution of primate lysozyme is one of the older results in molecular evolution, the enzyme having been recruited as a digestive protein in foregut fermenters.

Clicking the orange link shows a handful of synonymous differences and a Fisher p nowhere near significant. One pairwise comparison has little power; the published result rests on codon models across many primate lineages. Blue is the colour that tests strongly: a conserved gene accumulates measurable synonymous change while holding non-synonymous change near zero.

Checking the rates against the raw data

YEATS4 begins just past where LYZ ends, so the two share a locus and a divergence time and land at opposite ends of the ramp. It is conserved and compact, so its dS is low while its synonymous count clears the floor.

The script prints two genome-wide counts: pairs exceeding 1, and those surviving the Fisher test. The second is about what chance gives at that many tests.

Reproduce it end to end

build_primate_selection.sh runs everything above and writes a config.json with both assemblies, both gene tracks, the ortholog track and a session opening the locus.

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_primate_selection.sh
bash build_primate_selection.sh
npx --yes serve primate_selection_build/jbrowse2  # then open the printed URL

It needs the tools under Prerequisites on PATH.

See also

References