Skip to main content

Selection pressure between two genomes (dN/dS)

JBrowse v5 beta

The tutorials target the JBrowse v5 beta, and the v4.3.0 release on the download page lacks some of what they show. To install the beta of JBrowse Web, run npm install -g @jbrowse/cli@next, then jbrowse create jbrowse2 --branch v5.0.0-beta.13. Desktop beta builds are coming soon.

The ratio of non-synonymous to synonymous substitution rates, dN/dS, measures the selection pressure on a protein-coding gene between two species. jcvi builds a human against rhesus macaque ortholog table, kaks_from_pairs.py measures dN and dS on every pair, and JBrowse colors each ortholog link on a synteny track by the ratio. We then read the lysozyme neighbourhood on human chromosome 12.

Prerequisites

jcvi compiles C extensions at install time, and the compile fails on some Python versions. If pip install jcvi fails, uv venv --python 3.12 followed by uv pip install jcvi biopython gives an interpreter it builds on.

Where the data comes from

The 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.

The build script fetches these files, so there is nothing to download by hand.

What dN/dS says

A coding substitution is synonymous (the codon changes, the amino acid does not) or non-synonymous.

  • dS is the synonymous rate. Selection acts weakly on synonymous changes, so dS approximates the mutation rate.
  • dN is the non-synonymous rate, which reflects the selection acting on those changes.
  • dN/dS below 1 indicates purifying selection, where most genes sit. Above 1 needs positive selection to explain.

Producing the human-rhesus ortholog table and its rates

dS has to be large enough to estimate and small enough not to saturate. Rhesus macaque sits in that window against human.1

Calling human-rhesus orthologs with jcvi

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 input mistakes leave jcvi with no orthologs:

  • the alignment file has to be query = the first species, subject = the second. Reversed, jcvi looks up every id 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

Measuring dN and dS on each pair

pairs.tsv is the two gene columns of human.rhesus.anchors. The script skips human.rhesus.lifted.anchors, which recruits extra pairs near an established syntenic block. Their median dS is several times that of the chained pairs, which marks them as paralogs.

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/kaks_from_pairs.py
python3 kaks_from_pairs.py pairs.tsv both.cds.fa.gz \
  --key record --strip-version --min-syn-subs 3 --max-ds 0.3 -o primate.blocks
gzip -kf primate.blocks human.bed rhesus.bed

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 well above the cluster is a paralog the aligner preferred, and --max-ds drops every pair whose dS exceeds its value.

The top of a table sorted by dN/dS is the pairs with almost nothing to divide by: HBA1 lands there off a single synonymous difference. --min-syn-subs is a floor on that count, and drops every pair with fewer synonymous differences than its value.

Every row also has that count and a two-sided Fisher exact p, the test MEGA prescribes for small substitution counts. The track config below exposes both as attributeColumns, so clicking a link shows the evidence under its color.

Loading the human and rhesus genomes and genes

A gene-level synteny view reads no sequence, so each assembly is a .chrom.sizes file (a name and a length per line) that the script writes from the GFF3's ##sequence-region header. We load human here and rhesus the same way under the name rhesus:

Goes in the assemblies array of config.json. See Assemblies.

{ "name": "human", "uri": "human.chrom.sizes" }
jbrowse add-assembly human.chrom.sizes \
  --name human \
  --load copy

In JBrowse Desktop, Open new genome on the start screen (or File → Open genome... in a session), then Open from a URL and paste:

human.chrom.sizes

JBrowse reads the format off the file name. Then fill in:

  • Genome name: human

human.chrom.sizes is relative to a config.json. Replace it with its URL or its path on this computer.

Each panel draws genes from a track under its assembly. The script writes human.genes.gff3.gz and rhesus.genes.gff3.gz, bgzipped and tabix-indexed (prep). For your own annotation the GFF3 needs the same refNames as the assembly:

Goes in the tracks array of config.json. See Tracks.

{
  "type": "FeatureTrack",
  "trackId": "human_genes",
  "name": "human genes",
  "assemblyNames": ["human"],
  "adapter": {
    "type": "Gff3TabixAdapter",
    "uri": "human.genes.gff3.gz"
  }
}
jbrowse add-track human.genes.gff3.gz \
  --trackId human_genes \
  --name "human genes" \
  --assemblyNames human \
  --load copy

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": "FeatureTrack",
  "trackId": "human_genes",
  "name": "human genes",
  "assemblyNames": ["human"],
  "adapter": {
    "type": "Gff3TabixAdapter",
    "uri": "human.genes.gff3.gz"
  }
}

human.genes.gff3.gz is relative to a config.json. Replace it with its URL or its path on this computer.

Add the rhesus track the same way, with rhesus_genes under rhesus.

Loading the ortholog table as a synteny track

primate.blocks lists each pair's two gene ids, then dN, dS, the synonymous substitution count and the Fisher p, the .blocks shape MCScanBlocksAdapter reads. attributeColumns names the columns after the two gene columns, and the detail panel lists each as a feature attribute:

Goes in the tracks array of config.json. See Tracks.

{
  "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": ["human.bed.gz", "rhesus.bed.gz"],
    "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": ["human.bed.gz", "rhesus.bed.gz"],
    "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": ["human.bed.gz", "rhesus.bed.gz"],
    "attributeColumns": ["dn", "ds", "syn_subs", "fisher_p"]
  }
}

primate.blocks.gz, human.bed.gz, rhesus.bed.gz are relative to a config.json. Replace each with its URL or its path on this computer.

That shows it in the linear view. For the synteny view, open Add → Linear synteny view, pick the track under Quick start, and click Launch.

Color by value → dN/dS in the palette button menu sets the link's color, a color scale that reads dn and ds, on a ramp with 1 at the middle and 2 at the top.

Two view settings matter for a view this sparse:

  • Opacity on the sliders button in the view header sets opacity, which defaults to 0.25 for whole-genome views where ribbons overlap. 0.95 shows the color as it is.
  • Curved lines on the same menu sets drawCurves, which separates stacked neighbours.

Reading lysozyme (LYZ) against its neighbours by dN/dS

The session below opens, from the hosted copy, the stretch around lysozyme (LYZ) on human chromosome 12 where gene order matches rhesus, colored by dN/dS. In your own build, pick Color by value → dN/dS.

Goes at the top level of config.json, replacing any defaultSession there. See Default session.

{
  "defaultSession": {
    "name": "Selection pressure across a primate gene neighbourhood",
    "views": [
      {
        "type": "LinearSyntenyView",
        "views": [
          {
            "assembly": "human",
            "loc": "12:67,835,000-70,835,000",
            "tracks": ["human_genes"]
          },
          {
            "assembly": "rhesus",
            "loc": "11:67,401,000-70,319,000",
            "tracks": ["rhesus_genes"]
          }
        ],
        "tracks": [["primate_orthologs"]],
        "color": { "field": "dnds" },
        "opacity": 0.95,
        "drawCurves": true
      }
    ]
  }
}
jbrowse set-default-session --session - << 'EOF'
{
  "name": "Selection pressure across a primate gene neighbourhood",
  "views": [
    {
      "type": "LinearSyntenyView",
      "views": [
        {
          "assembly": "human",
          "loc": "12:67,835,000-70,835,000",
          "tracks": ["human_genes"]
        },
        {
          "assembly": "rhesus",
          "loc": "11:67,401,000-70,319,000",
          "tracks": ["rhesus_genes"]
        }
      ],
      "tracks": [["primate_orthologs"]],
      "color": { "field": "dnds" },
      "opacity": 0.95,
      "drawCurves": true
    }
  ]
}
EOF

Open this session in JBrowse Web ↗

Human against rhesus macaque across a stretch of conserved gene order on human chromosome 12, each ribbon one ortholog pair colored by dN/dS. Lysozyme (LYZ) is the one gene above the ramp's pivot; its neighbour YEATS4 is at the other end.
Open in JBrowse Web ↗Open in JBrowse Desktop ↗
Desktop link does nothing?

It needs JBrowse Desktop 5.0+. Or 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%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%22color%22%3A%7B%22field%22%3A%22dnds%22%7D%2C%22opacity%22%3A0.95%2C%22drawCurves%22%3Atrue%7D%5D%7D&sessionName=JBrowse+docs%3A+selection_pressure%2Flysozyme

This session loads a large dataset and takes a while to open.

Or rebuild it from your own data:

The short route: open the view in JBrowse Web, where File → Open track... adds your own files beside its hosted genome and tracks.

From scratch with your own data, the steps behind the figure:

  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. For JBrowse Web, run jbrowse add-assembly once for each to put them in your config.json.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 “Human / rhesus orthologs (dN/dS)”.
  4. Click Launch.
  5. Synteny view header → palette button → Color by value → dN/dSThe palette button's tooltip reads "Color by: ...".
  6. Synteny display settings (the sliders button in the view header) → Opacity → drag to 0.95Lower values let dense overlapping ribbons show through each other.
  7. Synteny display settings (the sliders button in the view header) → Curved lines (checked)
  8. Row 1 (human)
    1. Add a track: open this view's track selector and check one your config has, or click + and choose Add track to paste a URL or choose a local file. This one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “Human genes (GRCh38, Ensembl 116)”.
      1. Track menu → Show... → Labels → Name only
      2. Track menu → Feature height → Compact
    2. Type your region of interest into the location box and press Enter.This figure is at 12:67,835,000-70,835,000.
  9. Row 2 (rhesus)
    1. Add a track: open this view's track selector and check one your config has, or click + and choose Add track to paste a URL or choose a local file. This one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “Rhesus macaque genes (Mmul_10, Ensembl 116)”.
      1. Track menu → Show... → Labels → Name only
      2. Track menu → Feature height → Compact
    2. Type your region of interest into the location box and press Enter.This figure is at 11:67,401,000-70,319,000.

This session spec draws the figure. It goes after &session=spec- on a JBrowse Web link that loads https://jbrowse.org/demos/primate_selection/config.json, which is what the Web button above opens.

{
  "views": [
    {
      "type": "LinearSyntenyView",
      "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"
        ]
      ],
      "color": {
        "field": "dnds"
      },
      "opacity": 0.95,
      "drawCurves": true
    }
  ]
}

The same view with jbrowse-anywidget, from the config the figure loads. Swap a track's uri for your own file.

from jbrowse_anywidget import JBrowseApp, fetch_hub

config = fetch_hub("https://jbrowse.org/demos/primate_selection/config.json")
app = JBrowseApp(
    assemblies=config["assemblies"],
    tracks=config.get("tracks", []),
    views=[
        {
            "type": "LinearSyntenyView",
            "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"
                ]
            ],
            "color": {
                "field": "dnds"
            },
            "opacity": 0.95,
            "drawCurves": True
        }
    ],
)
app  # display the widget

jbrowse-img draws the same view as an SVG, from the command line and without a browser. Point --config at your own config, or swap a track for your own file with --bam, --vcfgz and the other file flags.

cat > session.json <<'JSON'
{
  "views": [
    {
      "type": "LinearSyntenyView",
      "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"
        ]
      ],
      "color": {
        "field": "dnds"
      },
      "opacity": 0.95,
      "drawCurves": true
    }
  ]
}
JSON

npx @jbrowse/img synteny --config https://jbrowse.org/demos/primate_selection/config.json \
  --spec session.json \
  --width 1500 --out figure.svg

--out figure.png writes a PNG instead, through rsvg-convert (apt install librsvg2-bin, brew install librsvg).

Hand this to a coding agent. It rebuilds the figure above headlessly, and the session file is then the thing to edit: change the location or a track's settings, or add your own track under sessionTracks, and rerun.

cat > session.json <<'JSON'
{
  "views": [
    {
      "type": "LinearSyntenyView",
      "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"
        ]
      ],
      "color": {
        "field": "dnds"
      },
      "opacity": 0.95,
      "drawCurves": true
    }
  ]
}
JSON

cat > callouts.json <<'JSON'
[
  {
    "type": "text",
    "text": "LYZ, the only pair above 1. Every neighbour is blue",
    "fontSize": 18,
    "maxWidth": 300,
    "anchor": {
      "view": [
        0,
        0
      ],
      "trackId": "human_genes",
      "loc": "12:69,351,000",
      "fracY": 1,
      "dy": 45,
      "alignX": "right",
      "dx": 95
    }
  },
  {
    "type": "arrow",
    "fromAnchor": {
      "view": [
        0,
        0
      ],
      "trackId": "human_genes",
      "loc": "12:69,351,000",
      "fracY": 1,
      "dy": 52,
      "alignX": "right",
      "dx": 88
    },
    "anchor": {
      "view": [
        0,
        0
      ],
      "trackId": "human_genes",
      "loc": "12:69,351,000",
      "fracY": 1,
      "dy": 52,
      "dx": 8
    }
  }
]
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/demos/primate_selection/config.json \
  --width 1500 --height 580 \
  --timeout 300000 \
  --spec session.json --annotations callouts.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.

Gene order is the same in both genomes here, so color is the only thing that varies. Click the orange LYZ ribbon for its synonymous count and Fisher p in the detail panel. Messier and Stewart reported adaptive evolution of primate lysozyme in 1997; the published result rests on codon models across many primate lineages.

YEATS4, a conserved neighbour, as the control for LYZ

YEATS4 begins just past where LYZ ends, so the two share a locus and a divergence time yet sit at opposite ends of the ramp. YEATS4 is conserved and compact, so its dS is low while its synonymous count clears the floor. Blue is the low end of the ramp, where the Fisher test does reach significance, because a conserved gene accumulates measurable synonymous change while its non-synonymous change stays near zero.

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:

  1. Turn each GFF3 into a BED and each CDS into a proteome keyed the same way, stripping Ensembl's transcript versions so the ids match.
  2. Align human against rhesus with DIAMOND and chain the hits with jcvi, keeping the chained anchors and not the lifted ones.
  3. Measure dN and dS on each pair, dropping pairs with too few synonymous differences or a dS far above the cluster.
  4. Print the count of pairs genome-wide with dN/dS above 1, and dN/dS and dS for every pair around LYZ, and write the config.
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

See also

Citations

Notes

  1. Chimpanzee leaves a denominator near zero on most genes. ↩