Skip to main content

Synteny visualization (a polyploid against itself)

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.

Hexaploid oat has three near-complete copies of its genome, one from each grass that hybridized to make it, so most of its genes exist three times. We draw oat against itself in a dotplot from the gene pairs between those copies, and color each pair by how far it has diverged at silent sites (dS) and at protein-changing ones (dN). The two subgenomes from closely related ancestors should pair at a lower dS than either does with the third. jcvi finds the pairs from a protein self-alignment, and kaks_from_pairs.py measures the rates.

Prerequisites

Selection pressure between two genomes covers installing jcvi and what dN/dS measures.

Where the data comes from

Oat cultivar Williams (Peng et al. 2022), annotated by Ensembl Plants release 63 as GCA_951802345.1.

The build script takes these files from their URLs, so there is nothing to download by hand.

Oat subgenomes and homoeologs

Oat (Avena sativa) is an allohexaploid with 21 chromosomes, seven homoeologous groups of three, one per subgenome (A, C and D). Copies across subgenomes are homoeologs, so a table of them pairs oat with itself and the dotplot puts oat on both axes. A segment moved between groups leaves the diagonal.

Producing the homoeolog table

Gene models, a proteome, and chromosome sizes

jcvi turns the GFF3 into the BED the adapter also reads, one primary transcript per gene:

python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id \
  --primary_only oat.gff3.gz -o oat.all.bed
awk -F'\t' '$1 ~ /^[1-7][ACD]$/' oat.all.bed > oat.bed

The end-to-end script translates the proteome from the CDS, keeping the transcript ids the BED uses.

The assembly needs no sequence, so we'll load it from the chromosome lengths in the GFF3's ##sequence-region header, written to oat.chrom.sizes (a name and a length per line):

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

{ "name": "oat", "uri": "oat.chrom.sizes" }
jbrowse add-assembly oat.chrom.sizes \
  --name oat \
  --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:

oat.chrom.sizes

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

  • Genome name: oat

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

Finding homoeolog pairs from a protein self-alignment

Passing the same prefix twice (oat oat) makes a self-comparison: jcvi drops the gene-against-itself diagonal, then chains the rest into syntenic blocks (co-linear runs of pairs).

diamond makedb --in oat.pep -d oat.pep
diamond blastp --threads 14 --query oat.pep --db oat.pep --out oat.oat.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 --self_remove 100 --no_dotplot oat oat

Two flags matter here:

  • --self_remove defaults to 98 and discards every hit at or above that percent identity. Oat's A-D homoeologs sit above it, so this run sets 100
  • --no_strip_names keeps the ids byte-identical to the BED the adapter joins on

The alignment is the long step. Running it separately keeps DIAMOND at default sensitivity, which finds homoeologs this recent. jcvi picks the file up by name and skips its alignment step.

Three choices decide which pairs survive:

  • Chaining keeps a gene pair (an anchor) only where its neighbours agree, which removes the off-diagonal noise of gene families' best hits.
  • The script keeps the pairs whose chromosomes have different subgenome letters, since a self-comparison also chains tandem and segmental duplicates within a subgenome, and writes them, two transcript ids per line, to oat.pairs.tsv.
  • The script reads oat.oat.anchors and skips oat.oat.lifted.anchors, because liftover recruits extra pairs near an established block, and their dS runs far above the chained ones', which marks them as paralogs.

Measuring dN and dS on each homoeolog pair

Ensembl declares dn and ds in every homology export and leaves both empty in every division, so the script computes them:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/kaks_from_pairs.py
python3 kaks_from_pairs.py oat.pairs.tsv oat.cds.fa.gz \
  --key record --min-syn-subs 3 -o oat.kaks.tsv

--key record reads the CDS by transcript id, and --min-syn-subs drops pairs with too few synonymous differences to trust a ratio. Selection pressure explains the method and the filters.

Loading the homoeolog table as a dotplot track

oat.kaks.tsv lists gene pairs, then dN, dS, the synonymous substitution count and a Fisher exact p, which is the .blocks layout MCScanBlocksAdapter reads. The script copies it to oat.homoeologs.blocks and gzips that and oat.bed for the track. A self-comparison names one assembly twice: in blockAssemblies, in the track's assemblyNames and in both entries of bedLocations:

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

{
  "type": "SyntenyTrack",
  "trackId": "oat_homoeologs",
  "name": "Oat homoeologs (dN/dS)",
  "assemblyNames": ["oat", "oat"],
  "adapter": {
    "type": "MCScanBlocksAdapter",
    "uri": "oat.homoeologs.blocks.gz",
    "blockAssemblies": ["oat", "oat"],
    "bedLocations": ["oat.bed.gz", "oat.bed.gz"],
    "attributeColumns": ["dn", "ds", "syn_subs", "fisher_p"]
  }
}
jbrowse add-track-json '{
  "type": "SyntenyTrack",
  "trackId": "oat_homoeologs",
  "name": "Oat homoeologs (dN/dS)",
  "assemblyNames": ["oat", "oat"],
  "adapter": {
    "type": "MCScanBlocksAdapter",
    "uri": "oat.homoeologs.blocks.gz",
    "blockAssemblies": ["oat", "oat"],
    "bedLocations": ["oat.bed.gz", "oat.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": "oat_homoeologs",
  "name": "Oat homoeologs (dN/dS)",
  "assemblyNames": ["oat", "oat"],
  "adapter": {
    "type": "MCScanBlocksAdapter",
    "uri": "oat.homoeologs.blocks.gz",
    "blockAssemblies": ["oat", "oat"],
    "bedLocations": ["oat.bed.gz", "oat.bed.gz"],
    "attributeColumns": ["dn", "ds", "syn_subs", "fisher_p"]
  }
}

oat.homoeologs.blocks.gz, oat.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.

attributeColumns lists the columns after the two gene columns, and each becomes a feature attribute in the detail panel. dn and ds drive the palette button's dN/dS.

Add → Dotplot view with oat on both axes opens the track as a dotplot, and the session below does the same.

Coloring homoeolog pairs by dS, with A-D pairs as the control

The script ends by printing the median dS for each subgenome pair, the numbers behind the picture below.

Color by value → ds on the palette button paints each pair by its dS. Oat's A and D subgenomes descend from closely related diploid Avena species and its C subgenome from a more distant one, so A-D pairs should come out at a lower dS than A-C or C-D pairs.

In grammar-of-graphics terms, the dotplot's color is a color scale, here the ds field mapped to a ramp. The session pins the ramp's ends with domainMin and domainMax, so a color means one dS wherever the view goes; the menu has no field for them:

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

{
  "defaultSession": {
    "name": "Oat homoeologs by dS",
    "views": [
      {
        "type": "DotplotView",
        "views": [
          {
            "assembly": "oat",
            "displayedRegionNames": [
              "4A",
              "4C",
              "4D",
              "5A",
              "5C",
              "5D",
              "7A",
              "7C",
              "7D"
            ]
          },
          {
            "assembly": "oat",
            "displayedRegionNames": [
              "4A",
              "4C",
              "4D",
              "5A",
              "5C",
              "5D",
              "7A",
              "7C",
              "7D"
            ]
          }
        ],
        "tracks": ["oat_homoeologs"],
        "color": {
          "field": "ds",
          "domainMin": 0,
          "domainMax": 0.25,
          "title": "dS"
        }
      }
    ]
  }
}
jbrowse set-default-session --session - << 'EOF'
{
  "name": "Oat homoeologs by dS",
  "views": [
    {
      "type": "DotplotView",
      "views": [
        {
          "assembly": "oat",
          "displayedRegionNames": [
            "4A",
            "4C",
            "4D",
            "5A",
            "5C",
            "5D",
            "7A",
            "7C",
            "7D"
          ]
        },
        {
          "assembly": "oat",
          "displayedRegionNames": [
            "4A",
            "4C",
            "4D",
            "5A",
            "5C",
            "5D",
            "7A",
            "7C",
            "7D"
          ]
        }
      ],
      "tracks": ["oat_homoeologs"],
      "color": {
        "field": "ds",
        "domainMin": 0,
        "domainMax": 0.25,
        "title": "dS"
      }
    }
  ]
}
EOF

Open this session in JBrowse Web ↗

Oat against itself over groups 4, 5 and 7, each homoeolog pair colored by dS on a pinned ramp. Cells pairing an A chromosome with a D one sit lower on the ramp than cells pairing either with C.
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%2Foat_homoeologs%2Fconfig.json&session=spec-%7B%22views%22%3A%5B%7B%22type%22%3A%22DotplotView%22%2C%22displayName%22%3A%22Oat+self-alignment%2C+dS%22%2C%22views%22%3A%5B%7B%22assembly%22%3A%22oat%22%2C%22displayedRegionNames%22%3A%5B%224A%22%2C%224C%22%2C%224D%22%2C%225A%22%2C%225C%22%2C%225D%22%2C%227A%22%2C%227C%22%2C%227D%22%5D%7D%2C%7B%22assembly%22%3A%22oat%22%2C%22displayedRegionNames%22%3A%5B%224A%22%2C%224C%22%2C%224D%22%2C%225A%22%2C%225C%22%2C%225D%22%2C%227A%22%2C%227C%22%2C%227D%22%5D%7D%5D%2C%22tracks%22%3A%5B%22oat_homoeologs%22%5D%2C%22color%22%3A%7B%22field%22%3A%22ds%22%2C%22domainMin%22%3A0%2C%22domainMax%22%3A0.25%2C%22title%22%3A%22dS%22%7D%7D%5D%7D&sessionName=JBrowse+docs%3A+homoeolog_synteny%2Foat_ds

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 has oat on X and oat on Y.
  2. Open the view: Add → Dotplot view, and set both axes to those assemblies in the import form.
  3. Point the import form at your own file: use the track panel under the two axis selectors, 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 both axes itself.This figure uses “Oat homoeologs (dN/dS)”.
  4. Click Launch.
  5. Click the view's title in its header and enter "Oat self-alignment, dS"
  6. Dotplot header → palette button → Color by value → dsThe palette button's tooltip reads "Color by: ...".

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

{
  "views": [
    {
      "type": "DotplotView",
      "displayName": "Oat self-alignment, dS",
      "views": [
        {
          "assembly": "oat",
          "displayedRegionNames": [
            "4A",
            "4C",
            "4D",
            "5A",
            "5C",
            "5D",
            "7A",
            "7C",
            "7D"
          ]
        },
        {
          "assembly": "oat",
          "displayedRegionNames": [
            "4A",
            "4C",
            "4D",
            "5A",
            "5C",
            "5D",
            "7A",
            "7C",
            "7D"
          ]
        }
      ],
      "tracks": [
        "oat_homoeologs"
      ],
      "color": {
        "field": "ds",
        "domainMin": 0,
        "domainMax": 0.25,
        "title": "dS"
      }
    }
  ]
}

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/oat_homoeologs/config.json")
app = JBrowseApp(
    assemblies=config["assemblies"],
    tracks=config.get("tracks", []),
    views=[
        {
            "type": "DotplotView",
            "displayName": "Oat self-alignment, dS",
            "views": [
                {
                    "assembly": "oat",
                    "displayedRegionNames": [
                        "4A",
                        "4C",
                        "4D",
                        "5A",
                        "5C",
                        "5D",
                        "7A",
                        "7C",
                        "7D"
                    ]
                },
                {
                    "assembly": "oat",
                    "displayedRegionNames": [
                        "4A",
                        "4C",
                        "4D",
                        "5A",
                        "5C",
                        "5D",
                        "7A",
                        "7C",
                        "7D"
                    ]
                }
            ],
            "tracks": [
                "oat_homoeologs"
            ],
            "color": {
                "field": "ds",
                "domainMin": 0,
                "domainMax": 0.25,
                "title": "dS"
            }
        }
    ],
)
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": "DotplotView",
      "displayName": "Oat self-alignment, dS",
      "views": [
        {
          "assembly": "oat",
          "displayedRegionNames": [
            "4A",
            "4C",
            "4D",
            "5A",
            "5C",
            "5D",
            "7A",
            "7C",
            "7D"
          ]
        },
        {
          "assembly": "oat",
          "displayedRegionNames": [
            "4A",
            "4C",
            "4D",
            "5A",
            "5C",
            "5D",
            "7A",
            "7C",
            "7D"
          ]
        }
      ],
      "tracks": [
        "oat_homoeologs"
      ],
      "color": {
        "field": "ds",
        "domainMin": 0,
        "domainMax": 0.25,
        "title": "dS"
      }
    }
  ]
}
JSON

npx @jbrowse/img dotplot --config https://jbrowse.org/demos/oat_homoeologs/config.json \
  --spec session.json \
  --width 900 --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": "DotplotView",
      "displayName": "Oat self-alignment, dS",
      "views": [
        {
          "assembly": "oat",
          "displayedRegionNames": [
            "4A",
            "4C",
            "4D",
            "5A",
            "5C",
            "5D",
            "7A",
            "7C",
            "7D"
          ]
        },
        {
          "assembly": "oat",
          "displayedRegionNames": [
            "4A",
            "4C",
            "4D",
            "5A",
            "5C",
            "5D",
            "7A",
            "7C",
            "7D"
          ]
        }
      ],
      "tracks": [
        "oat_homoeologs"
      ],
      "color": {
        "field": "ds",
        "domainMin": 0,
        "domainMax": 0.25,
        "title": "dS"
      }
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/demos/oat_homoeologs/config.json \
  --width 900 --height 767 \
  --timeout 300000 \
  --spec 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.

With dN/dS chosen on the palette button, almost every pair draws below 1 on the ramp. A ratio over 1 between copies this recently separated rests on few substitutions.

Reproduce it end to end

build_oat_homoeologs.sh runs everything above and writes a config.json with the assembly, the track and a dotplot session:

  1. Turn the GFF3 into a BED of one primary transcript per gene on the 21 chromosomes, and translate the CDS into a proteome keyed on the same ids.
  2. Align the proteome against itself with DIAMOND and chain the hits with jcvi, with --self_remove 100 so the A-D homoeologs survive.
  3. Keep the chained anchors whose two genes sit on different subgenomes, and measure dN and dS on each.
  4. Print the median dS per subgenome pair, the control the figure is read against, and write the config.
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_oat_homoeologs.sh
bash build_oat_homoeologs.sh
npx --yes serve oat_homoeologs_build/jbrowse2  # then open the printed URL

See also

Citations