Skip to main content

Synteny from gene symbols (44 E. coli genomes)

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.

We look at one K-12 operon across forty-three other E. coli and Shigella genomes at once, without aligning any of them. RefSeq's bacterial pipeline gives an orthologous gene the same symbol in every strain that names it, so the ortholog table is a join on the gene name over the GFF3 files. Each genome becomes a lane under the K-12 view, holding the gene models annotated in that genome. The join connects genes that share a symbol, which covers the core genome, and joins a symbol K-12 lacks between the strains that have it. The page ends at a cluster that differs between strains, where many genes in each lane draw grey with no ribbon.

Prerequisites

Where the data comes from

Forty-four RefSeq assemblies, each fetched by accession with the datasets CLI: the classic reference strains across phylogroups A, B1, B2, D and E, four Shigella, and complete genomes picked by striding a datasets summary listing. The five strains the pangenome graph and all-vs-all pages build from are all here under the same accessions (MG1655 is the strain those pages call K12). The build script pins and fetches every accession, so there is nothing to download by hand; the anchor and the four Shigella are:

A join on the gene symbols

RefSeq's Prokaryotic Genome Annotation Pipeline (PGAP) gives a gene the symbol of its ortholog (atpA is atpA in every strain that has it), so matching symbols across the GFF3 files builds an ortholog table with no alignment step: one row per gene, one column per genome, in the .blocks format that MCScanBlocksAdapter reads. Most of the effort goes to PGAP's unnamed and renamed genes. The download is an annotation and a sequence report per genome:

# one RefSeq accession per line; the first is the genome the rows are anchored on
datasets download genome accession --inputfile accessions.txt \
  --include gff3,seq-report --filename genomes.zip
unzip genomes.zip

Each lane takes its name from the strain field of the assembly report: MG1655 for "K-12 substr. MG1655", and Sflexneri_301 for the Shigella so it does not read as an E. coli strain.

PGAP writes a gene's locus tag into its Name when it has no symbol for it, so the join has to be told what an unnamed gene looks like, or every hypothetical protein would look named and match nothing:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/symbols_to_blocks.py
# --anchor names the genome whose genes lead the table; --unnamed is the
# locus-tag shape PGAP falls back to, which joins nothing; --merge-cited
# joins a gene PGAP renamed to the gene it was annotated from; --pseudogenes
# keeps a gene NCBI writes as a pseudogene feature, such as K-12's wbbL
python3 symbols_to_blocks.py --anchor MG1655 -o ecoli.blocks --unnamed '_RS[0-9]+$' --merge-cited --pseudogenes \
  MG1655=MG1655.gff.gz Sakai=Sakai.gff.gz CFT073=CFT073.gff.gz Sflexneri_301=Sflexneri_301.gff.gz

PGAP also renames genes between releases, so most strains here call K-12's gnd gndA. Each CDS lists the protein PGAP annotated the gene from (similar to AA sequence:RefSeq:NP_416533.1, the K-12 gnd protein), and --merge-cited joins the two symbols when that protein is in the table under the other name. Two symbols one genome has side by side, such as K-12's narH and its paralog narY, stay apart.

symbols_to_blocks.py reports how much of each column it filled, which tells whether a strain joins. Older PGAP runs gave genes a locus tag and no symbol, so a genome can be complete, current and join nothing. Counting the named genes in an annotation shows which genomes will fail before anything is built:

gzip -dc strain.gff.gz | awk -F'\t' '$3 == "gene" && $9 ~ /;gene=/' | wc -l

Loading the assemblies, gene tracks and ortholog track

Each strain is an assembly that needs only its chromosome's length, since the lanes never read sequence. The build writes <strain>.chrom.sizes (a name and a length, tab-separated) and loads it as a ChromSizesAdapter. A reader with ten strains repeats this block and the gene track below once per strain:

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

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

MG1655.chrom.sizes

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

  • Genome name: MG1655

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

Each lane draws gene models from the strain's gene track, the sorted, bgzipped and tabix-indexed GFF3 the build made above:

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

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

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

  • Main file: MG1655.gff.gz

Click Next. JBrowse reads the adapter and track type off the file name. Then fill in:

  • Track name: MG1655 genes
  • Assembly: MG1655

Click Add.

MG1655.gff.gz is relative to a config.json. Replace it with its URL or its path on this computer.

One SyntenyTrack names all forty-four assemblies, with blockAssemblies and bedLocations in the order the helper printed. The config below keeps the four genomes the command above joined, and the hosted config.json has all forty-four:

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

{
  "type": "SyntenyTrack",
  "trackId": "ecoli_orthologs",
  "name": "E. coli orthologs by gene symbol (44 genomes, RefSeq)",
  "assemblyNames": ["MG1655", "Sakai", "CFT073", "Sflexneri_301"],
  "adapter": {
    "type": "MCScanBlocksAdapter",
    "mcscanBlocksLocation": { "uri": "ecoli.blocks.gz" },
    "blockAssemblies": ["MG1655", "Sakai", "CFT073", "Sflexneri_301"],
    "bedLocations": [
      "MG1655.bed.gz",
      "Sakai.bed.gz",
      "CFT073.bed.gz",
      "Sflexneri_301.bed.gz"
    ]
  },
  "displays": [
    {
      "type": "MultiWaySyntenyDisplay",
      "color": { "field": "cluster" }
    }
  ]
}
jbrowse add-track-json '{
  "type": "SyntenyTrack",
  "trackId": "ecoli_orthologs",
  "name": "E. coli orthologs by gene symbol (44 genomes, RefSeq)",
  "assemblyNames": ["MG1655", "Sakai", "CFT073", "Sflexneri_301"],
  "adapter": {
    "type": "MCScanBlocksAdapter",
    "mcscanBlocksLocation": { "uri": "ecoli.blocks.gz" },
    "blockAssemblies": ["MG1655", "Sakai", "CFT073", "Sflexneri_301"],
    "bedLocations": [
      "MG1655.bed.gz",
      "Sakai.bed.gz",
      "CFT073.bed.gz",
      "Sflexneri_301.bed.gz"
    ]
  },
  "displays": [
    {
      "type": "MultiWaySyntenyDisplay",
      "color": { "field": "cluster" }
    }
  ]
}'

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": "ecoli_orthologs",
  "name": "E. coli orthologs by gene symbol (44 genomes, RefSeq)",
  "assemblyNames": ["MG1655", "Sakai", "CFT073", "Sflexneri_301"],
  "adapter": {
    "type": "MCScanBlocksAdapter",
    "mcscanBlocksLocation": { "uri": "ecoli.blocks.gz" },
    "blockAssemblies": ["MG1655", "Sakai", "CFT073", "Sflexneri_301"],
    "bedLocations": [
      "MG1655.bed.gz",
      "Sakai.bed.gz",
      "CFT073.bed.gz",
      "Sflexneri_301.bed.gz"
    ]
  },
  "displays": [
    {
      "type": "MultiWaySyntenyDisplay",
      "color": { "field": "cluster" }
    }
  ]
}

ecoli.blocks.gz, MG1655.bed.gz, Sakai.bed.gz, CFT073.bed.gz, Sflexneri_301.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.

Reading the atp operon across forty-four genomes

At the atp operon (the ATP synthase genes) on K-12, the track draws a lane per genome under the K-12 axis. Each lane is drawn in its genome's coordinates, and its header names the strain, contig and span shown, with [rev] where the lane runs against K-12 (lane headers). In grammar-of-graphics terms, the display's color is a color scale, here the cluster field mapped to colors:

  • Every gene is colored by its ortholog group, which the table names after the K-12 gene anchoring it, so a conserved gene is one color running down the whole stack.
  • A gene no group claims is grey, which marks the genes specific to a strain at a glance.
  • Show... → Show legend in the track menu turns on the key. The demo config hides it because the cluster window holds dozens of groups, and the display omits it in any window with more than thirty.

Lanes stack with the genome placing the most of the window's genes first. The default height scrolls the stack inside the track, so the session sets height to fit every lane:

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

{
  "defaultSession": {
    "name": "The atp operon across 44 genomes",
    "views": [
      {
        "type": "LinearGenomeView",
        "assembly": "MG1655",
        "loc": "NC_000913.3:3,910,000-3,925,000",
        "tracks": [
          {
            "trackId": "ecoli_orthologs",
            "type": "MultiWaySyntenyDisplay",
            "height": 1500
          }
        ]
      }
    ]
  }
}
jbrowse set-default-session --session - << 'EOF'
{
  "name": "The atp operon across 44 genomes",
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "MG1655",
      "loc": "NC_000913.3:3,910,000-3,925,000",
      "tracks": [
        {
          "trackId": "ecoli_orthologs",
          "type": "MultiWaySyntenyDisplay",
          "height": 1500
        }
      ]
    }
  ]
}
EOF

Open this session in JBrowse Web ↗

The atp operon on K-12 over forty-three E. coli and Shigella lanes, each drawing its own RefSeq gene models. A gene's color runs the full stack. Lanes marked [rev] hold this window's genes in the opposite order to K-12 and are mirrored, so their ribbons draw straight.
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%2Fecoli_orthologs%2Fconfig.json&session=spec-%7B%22views%22%3A%5B%7B%22type%22%3A%22LinearGenomeView%22%2C%22assembly%22%3A%22MG1655%22%2C%22loc%22%3A%22NC_000913.3%3A3%2C910%2C000-3%2C925%2C000%22%2C%22tracks%22%3A%5B%7B%22trackId%22%3A%22ecoli_orthologs%22%2C%22type%22%3A%22MultiWaySyntenyDisplay%22%2C%22height%22%3A1500%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+multiway_synteny%2Fecoli_symbol_atp_operon

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 genome: on the JBrowse Desktop start screen click Open new genome (or Show all available genomes to pick a hosted one). For JBrowse Web, jbrowse add-assembly puts it in your config.json.This figure uses MG1655.
  2. Add a track: check one your config has in the track selector, or File → Open track... to paste a URL or choose a local file.This figure uses “E. coli orthologs by gene symbol (44 genomes, RefSeq)”.
  3. Type your region of interest into the location box and press Enter.This figure is at NC_000913.3:3,910,000-3,925,000.

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

{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "MG1655",
      "loc": "NC_000913.3:3,910,000-3,925,000",
      "tracks": [
        {
          "trackId": "ecoli_orthologs",
          "type": "MultiWaySyntenyDisplay",
          "height": 1500
        }
      ]
    }
  ]
}

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/ecoli_orthologs/config.json")
app = JBrowseApp(
    assemblies=config["assemblies"],
    tracks=config.get("tracks", []),
    views=[
        {
            "type": "LinearGenomeView",
            "assembly": "MG1655",
            "loc": "NC_000913.3:3,910,000-3,925,000",
            "tracks": [
                {
                    "trackId": "ecoli_orthologs",
                    "type": "MultiWaySyntenyDisplay",
                    "height": 1500
                }
            ]
        }
    ],
)
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.

npx @jbrowse/img --config https://jbrowse.org/demos/ecoli_orthologs/config.json \
  --assembly MG1655 --loc NC_000913.3:3,910,000-3,925,000 \
  --track ecoli_orthologs display:MultiWaySyntenyDisplay '{"height":1500}' \
  --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": "LinearGenomeView",
      "assembly": "MG1655",
      "loc": "NC_000913.3:3,910,000-3,925,000",
      "tracks": [
        {
          "trackId": "ecoli_orthologs",
          "type": "MultiWaySyntenyDisplay",
          "height": 1500
        }
      ]
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/demos/ecoli_orthologs/config.json \
  --width 1500 --height 1740 \
  --timeout 240000 \
  --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.

The O-antigen cluster, where the gene-symbol join stops

The O-antigen cluster between galF and gnd holds the genes that build a strain's surface sugar chain. Each serotype has a different set, so the locus differs most between strains. Open the cluster:

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

{
  "defaultSession": {
    "name": "The O-antigen cluster across 44 genomes",
    "views": [
      {
        "type": "LinearGenomeView",
        "assembly": "MG1655",
        "loc": "NC_000913.3:2,095,000-2,115,000",
        "tracks": [
          {
            "trackId": "ecoli_orthologs",
            "type": "MultiWaySyntenyDisplay",
            "height": 1500
          }
        ]
      }
    ]
  }
}
jbrowse set-default-session --session - << 'EOF'
{
  "name": "The O-antigen cluster across 44 genomes",
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "MG1655",
      "loc": "NC_000913.3:2,095,000-2,115,000",
      "tracks": [
        {
          "trackId": "ecoli_orthologs",
          "type": "MultiWaySyntenyDisplay",
          "height": 1500
        }
      ]
    }
  ]
}
EOF

Open this session in JBrowse Web ↗

The O-antigen cluster on K-12 over the same forty-three lanes. The flanking ribbons run through every lane; the interior is mostly grey.
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%2Fecoli_orthologs%2Fconfig.json&session=spec-%7B%22views%22%3A%5B%7B%22type%22%3A%22LinearGenomeView%22%2C%22assembly%22%3A%22MG1655%22%2C%22loc%22%3A%22NC_000913.3%3A2%2C095%2C000-2%2C115%2C000%22%2C%22tracks%22%3A%5B%7B%22trackId%22%3A%22ecoli_orthologs%22%2C%22type%22%3A%22MultiWaySyntenyDisplay%22%2C%22height%22%3A1500%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+multiway_synteny%2Fecoli_symbol_oantigen

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 genome: on the JBrowse Desktop start screen click Open new genome (or Show all available genomes to pick a hosted one). For JBrowse Web, jbrowse add-assembly puts it in your config.json.This figure uses MG1655.
  2. Add a track: check one your config has in the track selector, or File → Open track... to paste a URL or choose a local file.This figure uses “E. coli orthologs by gene symbol (44 genomes, RefSeq)”.
  3. Type your region of interest into the location box and press Enter.This figure is at NC_000913.3:2,095,000-2,115,000.

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

{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "MG1655",
      "loc": "NC_000913.3:2,095,000-2,115,000",
      "tracks": [
        {
          "trackId": "ecoli_orthologs",
          "type": "MultiWaySyntenyDisplay",
          "height": 1500
        }
      ]
    }
  ]
}

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/ecoli_orthologs/config.json")
app = JBrowseApp(
    assemblies=config["assemblies"],
    tracks=config.get("tracks", []),
    views=[
        {
            "type": "LinearGenomeView",
            "assembly": "MG1655",
            "loc": "NC_000913.3:2,095,000-2,115,000",
            "tracks": [
                {
                    "trackId": "ecoli_orthologs",
                    "type": "MultiWaySyntenyDisplay",
                    "height": 1500
                }
            ]
        }
    ],
)
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.

npx @jbrowse/img --config https://jbrowse.org/demos/ecoli_orthologs/config.json \
  --assembly MG1655 --loc NC_000913.3:2,095,000-2,115,000 \
  --track ecoli_orthologs display:MultiWaySyntenyDisplay '{"height":1500}' \
  --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": "LinearGenomeView",
      "assembly": "MG1655",
      "loc": "NC_000913.3:2,095,000-2,115,000",
      "tracks": [
        {
          "trackId": "ecoli_orthologs",
          "type": "MultiWaySyntenyDisplay",
          "height": 1500
        }
      ]
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/demos/ecoli_orthologs/config.json \
  --width 1500 --height 1740 \
  --timeout 240000 \
  --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.

Read the lanes from the top:

  • The K-12 derivatives (DH10B, HMS174, C3026, MGY, tolC and MG1655_TMP32XR1) sort to the top and have the cluster gene for gene.
  • Below them the flanks chain down every lane: galF on one side, gnd (gndA in most lanes), ugd and wzzB on the other.
  • The interior is grey, except the rfb genes, wzx and wzy, which join wherever a strain has them.

A row for a symbol K-12 lacks joins only the lanes with that gene, so a serotype gene a few strains share, such as vioB in IAI39 and the BL21 lanes, runs down those lanes and skips the rest. A gene PGAP left under its locus tag, such as ECOLC_RS24020 in the ATCC_8739 lane, has no row.

Most of the grey is each serotype's sugar genes, which arrived by horizontal transfer and have no K-12 counterpart to join. An OrthoFinder run would add genes the annotations named differently, and the all-vs-all alignment draws this locus base by base for five strains.

Reproduce it end to end

The script needs the tools under Prerequisites and works in four steps:

  1. Download each genome's annotation, sequence report and assembly report by accession.
  2. Name each lane after the report's strain field, and drop a genome the report calls neither E. coli nor one of the four Shigella species.
  3. Keep each genome's longest sequence as its chromosome, so a lane follows one contig and no plasmid, and sort and index the GFF3 as in the web quickstart.
  4. Join the annotations on gene symbol, treating PGAP's locus-tag names as unnamed and merging the symbols PGAP renamed, and write the config.
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_ecoli_orthologs.sh
bash build_ecoli_orthologs.sh

See also

Citations