# Synteny from gene symbols (44 E. coli genomes)

**TL;DR:** 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 it names, so
the ortholog table is a join on the gene name over the GFF3 files, and each
genome becomes a lane under the K-12 view carrying its own gene models. The join
sees the core genome and nothing else: at a cluster that differs between strains
the lanes draw their own genes and no ribbons, which is where the page ends.

## Prerequisites

- The
  [NCBI datasets CLI](https://www.ncbi.nlm.nih.gov/datasets/docs/v2/command-line-tools/)
- htslib (`bgzip`, `tabix`)
- `python3`
- A running JBrowse instance (the [web quickstart](https://jbrowse.org/jb2/docs/quickstart_web) or the
  [desktop quickstart](https://jbrowse.org/jb2/docs/quickstart_desktop))

## 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.
Nothing about a genome is typed by hand: the zip `datasets` writes carries an
assembly report, and each lane's name is that report's strain field (MG1655
rather than "K-12 substr. MG1655", `Sflexneri_301` for the Shigella so it does
not read as an E. coli strain), and the report's organism name is a screen of
its own — three accessions the listing handed over turned out to be a Leclercia,
a Salmonella and a suppressed _Shigella_ sp., and the build drops them and says
so. The five strains the [pangenome graph](https://jbrowse.org/jb2/docs/tutorials/pangenome_ecoli) and
[all-vs-all](https://jbrowse.org/jb2/docs/tutorials/allvsall_synteny) pages build from are all here
under the same accessions — MG1655 is the strain those pages call K12 — so the
three pages read one set of genomes three ways. The
[build script](#reproduce-it-end-to-end) pins every accession; the anchor and
the four Shigella are:

- K-12 MG1655:
  https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/
- Shigella flexneri 301:
  https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/006/925/GCF_000006925.2_ASM692v2/
- Shigella dysenteriae Sd197:
  https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/012/005/GCF_000012005.1_ASM1200v1/
- Shigella boydii Sb227:
  https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/012/025/GCF_000012025.1_ASM1202v1/
- Shigella sonnei 53G:
  https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/283/715/GCF_000283715.1_ASM28371v1/
- the finished table, BEDs and config, rehosted so the lanes load without
  rerunning the pipeline: https://jbrowse.org/demos/ecoli_orthologs/config.json

## A join instead of an alignment

The [all-vs-all page](https://jbrowse.org/jb2/docs/tutorials/allvsall_synteny) describes five E. coli
strains to each other by aligning them, which is one minimap2 run per pair and
what holds that demo at five. RefSeq's prokaryotic annotation pipeline names a
gene by its ortholog (_atpA_ is _atpA_ in every strain that carries it), so the
same `.blocks` table can be filled by matching symbols across the GFF3 files,
the route the [primate page](https://jbrowse.org/jb2/docs/tutorials/primate_orthologs_synteny) takes
for eight apes. The download is an annotation and a sequence report per genome:

<!-- from: scripts/build_ecoli_orthologs.sh -->

```bash
# 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 genome's chromosome is the longest sequence in its report, and the plasmids
are dropped: a lane follows one contig at a time. The GFF3 filtered to that
sequence, sorted, bgzipped and tabix-indexed as in the
[web quickstart](https://jbrowse.org/jb2/docs/quickstart_web) is the genome's gene track.

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:

<!-- from: scripts/build_ecoli_orthologs.sh -->

```bash
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/symbols_to_blocks.py
# --anchor names the genome whose genes are the rows; --unnamed is the
# locus-tag shape PGAP falls back to, which joins nothing
python3 symbols_to_blocks.py --anchor MG1655 -o ecoli.blocks --unnamed '_RS[0-9]+$' \
  MG1655=MG1655.gff.gz Sakai=Sakai.gff.gz CFT073=CFT073.gff.gz Sflexneri_301=Sflexneri_301.gff.gz
```

The helper reports how much of each column it filled, and that number is the
screen a strain has to pass. Older PGAP runs named genes by locus tag alone, so
a genome can be complete, current and join nothing; of a hundred accessions
tried for this page, a quarter of those that still downloaded were like that,
and two were phage genomes. Counting the named genes in an annotation says which
before anything is built:

<!-- from: scripts/build_ecoli_orthologs.sh -->

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

## The ortholog track

One `SyntenyTrack` names all forty-four assemblies, each of which is a
`ChromSizesAdapter` over its chromosome's length, since the lanes never read
sequence. `blockAssemblies` and `bedLocations` are positional against the
table's columns, in the order the helper printed; the list below is cut to the
first four for the page. The table and the BEDs are gzipped, which the adapter
undoes itself, since it reads each file whole before the first lane draws. At
the default height the display keeps every lane at a readable pitch and scrolls
the stack inside the track; the `height` here sizes the track to the whole
stack, so the figure below shows every lane at once:

```json addtrack
{
  "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": [
      { "uri": "MG1655.bed.gz" },
      { "uri": "Sakai.bed.gz" },
      { "uri": "CFT073.bed.gz" },
      { "uri": "Sflexneri_301.bed.gz" }
    ],
    "assemblyNames": ["MG1655", "Sakai", "CFT073", "Sflexneri_301"]
  },
  "displays": [
    {
      "type": "MultiWaySyntenyDisplay",
      "displayId": "ecoli_orthologs-MultiWaySyntenyDisplay",
      "color": "jexl:feature.name ? randomColor(feature.name) : '#b0b0b0'",
      "height": 970
    }
  ]
}
```

## One operon, forty-four genomes

Opened on K-12 at the _atp_ operon, the track draws a lane per genome under the
K-12 axis. Every gene is colored by its symbol, so a conserved gene is one color
running down the whole stack, and a lane's header names its chromosome, where it
is looking and `[rev]` where the strain's chromosome reads the other way. Lanes
stack densest first, so the genomes placing the most of the window sit at the
top and the reduced Shigella genomes fall toward the bottom without anything
naming them.

```json session config=https://jbrowse.org/demos/ecoli_orthologs/config.json
{
  "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": 970
          }
        ]
      }
    ]
  }
}
```

<Figure caption="The atp operon on K-12 over forty-three E. coli and Shigella lanes from one gene-symbol ortholog track, each lane drawing its own RefSeq gene models. Every gene's color runs the full stack; the lanes reading the operon reversed are the ones whose chromosome was deposited the other way round." src="/img/multiway_synteny/ecoli_symbol_atp_operon.png" />

## Where the join stops

The O-antigen cluster between _galF_ and _gnd_ is the locus that differs most
between strains: each serotype carries its own set of sugar pathway genes, and
those share no symbol across serotypes. The lanes that keep the whole cluster
sort to the top, and their headers say why: DH10B, HMS174, C3026, MGY, tolC- and
MG1655_TMP32XR1 are all K-12 derivatives, whose cluster is K-12's gene for gene.
Below them the flanking genes chain down every lane, _galF_ on one side and
_wzzB_ and _ugd_ on the other (_gnd_ itself is a locus tag in most PGAP
annotations, so it joins only a handful of lanes), and between the flanks each
lane draws whatever its own annotation holds with no ribbon to it. That is the
accessory genome as the join sees it.

<Figure caption="The O-antigen cluster on K-12 over the same forty-three lanes. The six K-12 derivatives at the top of the stack match the cluster gene for gene; in every lane below, the flanking galF, wzzB and ugd chains run through and the cluster between them is that strain's own, joined to nothing." src="/img/multiway_synteny/ecoli_symbol_oantigen.png" />

For the variable loci the table wants a homology call across the proteomes, an
[OrthoFinder](https://jbrowse.org/jb2/docs/tutorials/orthofinder_synteny) run, or the
[all-vs-all alignment](https://jbrowse.org/jb2/docs/tutorials/allvsall_synteny) that draws the same
locus base by base for five strains.

## Reproduce it end to end

The script fetches the annotations, names and screens the genomes off NCBI's
assembly report, builds the table and writes the config; see
[Prerequisites](#prerequisites).

```bash
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_ecoli_orthologs.sh
bash build_ecoli_orthologs.sh
```

## See also

- [](https://jbrowse.org/jb2/docs/tutorials/primate_orthologs_synteny)
- [](https://jbrowse.org/jb2/docs/tutorials/allvsall_synteny)
- [](https://jbrowse.org/jb2/docs/tutorials/pangenome_ecoli)

## References

- Li W, et al. RefSeq: expanding the Prokaryotic Genome Annotation Pipeline
  reach with protein family model curation. Nucleic Acids Res (2021).
  https://doi.org/10.1093/nar/gkaa1105
- Touchon M, et al. Organised genome dynamics in the Escherichia coli species
  results in highly diverse adaptive paths. PLoS Genet (2009).
  https://doi.org/10.1371/journal.pgen.1000344
- NCBI Datasets: https://www.ncbi.nlm.nih.gov/datasets/