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

Every coding substitution between two orthologs is either synonymous, changing the codon but not the amino acid, or non-synonymous. Synonymous changes are close to invisible to selection, so their rate dS is roughly the rate at which mutations arrive and fix. Non-synonymous changes are seen, so their rate dN carries whatever selection did.

The ratio is therefore read against 1. Below it, amino acid changes were removed faster than silent ones, which is purifying selection and is where most genes sit most of the time. Above it, amino acid changes fixed faster than silent ones, which takes positive selection to explain.

Producing the data

dS has to be large enough to estimate and small enough not to saturate, and 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 traps, both silent.

The alignment file has to be query = the first species, subject = the second, which is the order jcvi would have used had it run the aligner itself. 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, so the proteome ends up naming ENST00000641515.7 where the BED names ENST00000641515. Nothing matches and nothing says so. The script strips the version; kaks_from_pairs.py takes --strip-version for the same reason.

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

Two species diverged once, so their true orthologs share a divergence time and their dS values cluster. A pair whose dS comes out an order of magnitude above that cluster is a paralog the aligner preferred, and the script's --max-ds removes it.

Sorting the table by dN/dS and reading off the top returns the pairs with almost nothing to divide by: HBA1, about as strongly conserved as a gene gets, comes out over 2 off a single synonymous difference, and so do the others near the top. --min-syn-subs is a floor on that count. dS is per site, so the same rate is much weaker evidence in a short gene than a long one.

Every row also carries that count and a two-sided Fisher exact p, which is the test MEGA prescribes when the numbers of substitutions are small, where the large-sample Z-test over-rejects. They are attributeColumns like the rates, so clicking a link shows how much evidence is 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"]
  }
}'

attributeColumns names the columns after the two gene columns, so each becomes a feature attribute visible in the detail panel, and dn with ds together drive Color by... → dN/dS. syn_subs and fisher_p are what a reader checks a colour against.

That ramp is fixed, with 1 at its middle and 2 at its top, so a pair's colour says which side of 1 it falls on.

Two settings matter for a view this sparse, and both are properties of the LinearSyntenyView rather than of the track. alpha defaults to 0.2, tuned for whole-genome views where thousands of ribbons overlap; at 0.95 the colour is the colour. drawCurves renders the links as beziers, which 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 the ribbons run parallel and colour is the only thing that varies across them. Lysozyme is a good gene to find there: 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.

Click the orange link and the detail panel gives the count and the p behind it: a handful of synonymous differences, and a Fisher p nowhere near significant. One pairwise comparison carries very little power, and the published result rests on codon models across many primate lineages.

The blue is the colour that tests strongly here: a conserved gene accumulates enough synonymous change to measure while holding non-synonymous change near zero. Across the whole table the great majority of pairs sit significantly below 1 and almost none significantly above.

Checking the rates against the raw data

The figure carries its own control. YEATS4 begins just past where LYZ ends, so the two share a locus, a divergence time and a neighbourhood, and they land at opposite ends of the ramp. YEATS4 is also the pair the substitution-count floor keeps: it is conserved and compact, so its dS is low while its synonymous count is adequate.

The script prints the neighbourhood beside the genome-wide distribution as two counts: how many pairs exceed 1, and how many of those survive the Fisher test. Few do the first, and those that survive the second are about what chance alone would give at that many tests.

Reproduce it end to end

build_primate_selection.sh runs everything above: it downloads both annotations from Ensembl, builds the BEDs, proteomes and gene tracks, runs DIAMOND and jcvi, measures dN and dS on every ortholog pair, prints the tables in the section above, downloads JBrowse, 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