Synteny visualization (OrthoFinder orthogroups)
TL;DR: an OrthoFinder run groups
genes by homology and says nothing about position, so it produces a synteny
track for genomes no aligner can line up. Orthogroups.tsv converts to the
.blocks table MCScanBlocksAdapter reads, with one BED per genome, and a
duplicated gene becomes several rows.
Prerequisites
- OrthoFinder with DIAMOND
python3- htslib (
bgzip,tabix) wgetnode, for the JBrowse CLI- The
NCBI datasets CLI
(
datasetsanddataformat), for the sets that name a genome's chromosomes from a sequence report:wheat,drosophilaandsolanaceae - gffread, which translates a genome and its annotation into the proteome OrthoFinder compares. Needed for your own genomes; the sets above download a published proteome instead
- A running JBrowse instance (the web quickstart or the desktop quickstart)
OrthoFinder bundles its own DIAMOND: bioconda has
orthofinder, and the davidemms/orthofinder container has the same bundle.
Everything here calls orthofinder on PATH; the
build script header has an
Apptainer shim for the container.
Where the data comes from
Five OrthoFinder sets, one per section below, each built from one Ensembl division's protein FASTA and GFF3 per genome.
vertebrates: human, chicken, frog, gar, zebrafish, Ensembl release 113: https://ftp.ensembl.org/pub/release-113/wheat: Aegilops tauschii, bread wheat, durum, wild emmer, Triticum urartu, T. timopheevii, Ensembl Plants release 63: https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-63/drosophila: Drosophila melanogaster, D. simulans, D. yakuba, D. pseudoobscura, D. virilis, Ensembl Metazoa release 63: https://ftp.ensemblgenomes.ebi.ac.uk/pub/metazoa/release-63/solanaceae: tomato, potato, pepper and Nicotiana attenuata, with coffee as the outgroup, Ensembl Plants release 63: https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-63/grasses: rice, sorghum, maize, brachypodium, setaria, Ensembl Plants release 63: https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-63/- Sequence reports for the six assemblies whose refName aliases are read from NCBI rather than a URL, which the datasets CLI fetches by accession: T. timopheevii GCA_963921465.1, tomato GCA_000188115.5, and the four flies other than D. melanogaster, GCA_016746395.2, GCA_016746365.2, GCA_009870125.2 and GCA_030788295.1: https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/
Orthogroups as a synteny source
OrthoFinder clusters proteins into orthogroups without reference to position, so a table of orthogroups is a synteny track for a human against a zebrafish, past the divergence where minimap2 and MCScan return nothing. Nothing in an orthogroup asserts synteny, so any collinearity in the ribbons is a property of the genomes.
Five sets follow, a section each; the second half of the page builds wheat
from Orthogroups.tsv.
Vertebrates: blocks that survive out to zebrafish
The vertebrates set is human, chicken, frog, spotted gar and zebrafish. Their
orthologs still fall into chromosome-scale blocks, and the teleost genome
duplication shows as a human chromosome answering to more zebrafish chromosomes
than chicken ones.
A band resolves into wedges only where a chromosome's orthologs mostly land on one chromosome of the row below; the build script prints that share per pair.
One locus, one lane per vertebrate
For one locus, a multi-way synteny track draws the same table as a lane per genome in a single linear view, each in its own coordinates.
{
"defaultSession": {
"name": "Vertebrate multi-way synteny track",
"views": [
{
"type": "LinearGenomeView",
"assembly": "human",
"loc": "2:176,090,000-176,290,000",
"tracks": [
{
"trackId": "human_genes",
"type": "LinearBasicDisplay",
"showOnlyGenes": true,
"displayMode": "compact"
},
{
"trackId": "vertebrates_orthogroups",
"type": "MultiWaySyntenyDisplay",
"rowOrder": ["chicken", "frog", "gar", "zebrafish"],
"height": 320
}
]
}
]
}
}jbrowse set-default-session --session - << 'EOF'
{
"name": "Vertebrate multi-way synteny track",
"views": [
{
"type": "LinearGenomeView",
"assembly": "human",
"loc": "2:176,090,000-176,290,000",
"tracks": [
{
"trackId": "human_genes",
"type": "LinearBasicDisplay",
"showOnlyGenes": true,
"displayMode": "compact"
},
{
"trackId": "vertebrates_orthogroups",
"type": "MultiWaySyntenyDisplay",
"rowOrder": ["chicken", "frog", "gar", "zebrafish"],
"height": 320
}
]
}
]
}
EOFThe human HOXD cluster is one block every one of the four keeps.
Wheat: six genomes of one polyploid history
The wheat set is six genomes of wheat's own polyploidy history:
- Aegilops tauschii, the diploid D-genome donor
- bread wheat, hexaploid, genomes A+B+D
- durum, the domesticated tetraploid, A+B
- wild emmer, durum's wild tetraploid ancestor
- Triticum urartu, the diploid A-genome donor
- T. timopheevii, a second, independent tetraploid that also traces to the A-genome donor
Stacked in that order, each adjacent pair is an evolutionary step. Diagonalizing brings each homoeologous group's chromosomes together, with nothing telling the view which chromosomes are homoeologs.
Show all regions - same bp per pixel in the view menu (sameScale in a
session spec) puts the rows on one bp/px, so a row's drawn length is its genome
size.
Reading one chromosome out of the stack
Any pair of the set opens as a two-row view. This one puts Aegilops tauschii's seven chromosomes over bread wheat 4A, with the palette button's Query painting each link by the tauschii chromosome it leaves.
6D and the other three uninvolved chromosomes reach 4A only as scattered singletons. The blocks are the 4AL/5AL and 4AL/7BS translocation pair, RFLP-mapped by Devos et al. and revisited by Dvorak et al.; a group number is shared across subgenomes, so the 5D and 7D bundles are the D-genome counterparts of the 5A and 7B those names refer to.
Bread wheat 4A against Triticum urartu
Triticum urartu, the A-genome donor, is the same two-row view with one assembly swapped.
Urartu's chromosome 4 covers the first two blocks and its chromosome 7 the distal one.
Drosophila: chromosome arms that outlast gene order
The drosophila set is D. melanogaster, its near relatives D. simulans and
D. yakuba, and the distant D. pseudoobscura and D. virilis. Flies keep
their chromosome arms (Muller's elements) across that range and rewrite the gene
order inside them.
The conversion prints, per adjacent pair, the share of a chromosome's links landing on its single best partner in the row below:
chromosome-level correspondence, each row against the next:
melanogaster -> simulans best partner holds 98% of a chromosome's links (5 chromosomes)
simulans -> yakuba best partner holds 89% of a chromosome's links (5 chromosomes)
yakuba -> pseudoobscura best partner holds 86% of a chromosome's links (5 chromosomes)
pseudoobscura -> virilis best partner holds 77% of a chromosome's links (4 chromosomes)Stacked, each melanogaster arm's colour arrives as one bundle in every row.
One locus, one lane per fly
Gene order needs a window. A multi-way synteny track draws a lane per fly in its own coordinates, the header naming the chromosome that fly keeps these orthologs on.
{
"defaultSession": {
"name": "Drosophila multi-way synteny track",
"views": [
{
"type": "LinearGenomeView",
"assembly": "melanogaster",
"loc": "3L:5,789,000-5,931,000",
"tracks": [
{
"trackId": "melanogaster_genes",
"type": "LinearBasicDisplay",
"showOnlyGenes": true,
"displayMode": "compact"
},
{
"trackId": "drosophila_orthogroups",
"type": "MultiWaySyntenyDisplay",
"rowOrder": ["simulans", "yakuba", "pseudoobscura", "virilis"],
"height": 320
}
]
}
]
}
}jbrowse set-default-session --session - << 'EOF'
{
"name": "Drosophila multi-way synteny track",
"views": [
{
"type": "LinearGenomeView",
"assembly": "melanogaster",
"loc": "3L:5,789,000-5,931,000",
"tracks": [
{
"trackId": "melanogaster_genes",
"type": "LinearBasicDisplay",
"showOnlyGenes": true,
"displayMode": "compact"
},
{
"trackId": "drosophila_orthogroups",
"type": "MultiWaySyntenyDisplay",
"rowOrder": ["simulans", "yakuba", "pseudoobscura", "virilis"],
"height": 320
}
]
}
]
}
EOFEvery fly keeps all of the melanogaster genes in this 3L window, from Bre1 to
PXo. The near relatives keep them in order; the distant ones reverse the
block, which the header says as [rev].
Muller element D is melanogaster's 3L, and in the obscura lineage it is fused to the X, so the pseudoobscura lane sits on a chromosome the assembly calls X.
Nightshades: the same genes over four times the DNA
The solanaceae set is tomato, potato, pepper, Nicotiana attenuata and coffee
as the outgroup. Their gene counts are similar and their genome sizes are not,
the difference being repeat sequence between the genes. On one bp/px a row's
length is its genome size.
N. attenuata is still on scaffolds, which the correspondence print shows as a low best-partner share. Its genes are spread over thousands of sequences, of which the build keeps the densest, so the row draws the share that fell on those.
One locus, five lanes, five scales
A multi-way synteny track makes the size difference per-gene: each lane is fitted to the window's orthologs in its own coordinates and says what scale that took.
{
"defaultSession": {
"name": "Nightshade multi-way synteny track",
"views": [
{
"type": "LinearGenomeView",
"assembly": "tomato",
"loc": "SL4.0ch04:62,880,000-63,037,000",
"tracks": [
{
"trackId": "tomato_genes",
"type": "LinearBasicDisplay",
"showOnlyGenes": true,
"displayMode": "compact",
"showLabels": "none"
},
{
"trackId": "solanaceae_orthogroups",
"type": "MultiWaySyntenyDisplay",
"rowOrder": ["potato", "pepper", "tobacco", "coffee"],
"height": 320
}
]
}
]
}
}jbrowse set-default-session --session - << 'EOF'
{
"name": "Nightshade multi-way synteny track",
"views": [
{
"type": "LinearGenomeView",
"assembly": "tomato",
"loc": "SL4.0ch04:62,880,000-63,037,000",
"tracks": [
{
"trackId": "tomato_genes",
"type": "LinearBasicDisplay",
"showOnlyGenes": true,
"displayMode": "compact",
"showLabels": "none"
},
{
"trackId": "solanaceae_orthogroups",
"type": "MultiWaySyntenyDisplay",
"rowOrder": ["potato", "pepper", "tobacco", "coffee"],
"height": 320
}
]
}
]
}
EOFEvery genome keeps the two dozen genes in this tomato window. Pepper and N. attenuata need several times the anchor's span for them, the intergenic expansion arriving as a number in a lane header.
Every lane's genes stay in the anchor's order; the coffee lane is [rev].
Grasses: a whole-genome duplication only maize has
The grasses set is rice, sorghum, maize, brachypodium and foxtail millet.
Maize carries a whole-genome duplication the other four do not.
The ribbon count is a conversion setting; what to do with a duplicated gene is where the build picks it.
One rice window, one lane per grass
One window of the duplication, as a multi-way synteny track under rice's gene track:
{
"defaultSession": {
"name": "Grasses multi-way synteny track",
"views": [
{
"type": "LinearGenomeView",
"assembly": "rice",
"loc": "3:31,590,000-31,775,000",
"tracks": [
{
"trackId": "rice_genes",
"type": "LinearBasicDisplay",
"showOnlyGenes": true,
"displayMode": "compact"
},
{
"trackId": "grasses_orthogroups",
"type": "MultiWaySyntenyDisplay",
"rowOrder": ["sorghum", "brachypodium", "setaria", "maize"],
"height": 320
}
]
}
]
}
}jbrowse set-default-session --session - << 'EOF'
{
"name": "Grasses multi-way synteny track",
"views": [
{
"type": "LinearGenomeView",
"assembly": "rice",
"loc": "3:31,590,000-31,775,000",
"tracks": [
{
"trackId": "rice_genes",
"type": "LinearBasicDisplay",
"showOnlyGenes": true,
"displayMode": "compact"
},
{
"trackId": "grasses_orthogroups",
"type": "MultiWaySyntenyDisplay",
"rowOrder": ["sorghum", "brachypodium", "setaria", "maize"],
"height": 320
}
]
}
]
}
EOFA lane carries one refName, so the maize lane shows one of the two copies.
Both maize copies at once is the stacked view's job. Launch → Linear synteny view (visible region) in the lane track's menu offers a full row per grass over the visible window.
Producing the blocks table
The commands from here on build the wheat set; the other four differ only in
which proteomes go into the directory:
# -og stops after the orthogroups, skipping the gene trees and the species
# tree, which this table does not use and which are most of the runtime.
# -S diamond picks the aligner; -t is threads.
orthofinder -f proteomes -og -S diamond -t "$(getconf _NPROCESSORS_ONLN)"It writes proteomes/OrthoFinder/Results_<date>/Orthogroups/Orthogroups.tsv.
That table needs the header row and the leading id column dropped, each cell
reduced to a gene id, and an empty cell marked .:
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/orthogroups_to_blocks.py
python3 orthogroups_to_blocks.py Orthogroups.tsv -o tauschii.blocks \
--bed tauschii=tauschii.bed --bed wheat=wheat.bed --bed durum=durum.bed \
--bed emmer=emmer.bed --bed urartu=urartu.bed --bed timopheevii=timopheevii.bedThe table is named after one genome by convention only. The script prints the
column order off the header row, which is what blockAssemblies has to be. A
column is named after its proteome file minus extensions;
--assembly COLUMN=NAME renames one, and a key matching no column is an error.
What to do with a duplicated gene
A cell holds every gene of that genome in the orthogroup, and a synteny link runs gene to gene, so a cell holding two genes has no single answer. In the grasses a rice gene commonly has two maize orthologs. Three treatments:
--pick | A rice gene with two maize orthologs | Use when |
|---|---|---|
first | one ribbon, to whichever maize gene OrthoFinder listed first | you want maximum coverage and accept the arbitrary choice |
expand (default) | two ribbons, one per maize copy | the duplication is part of what you are looking at |
single | no ribbon | you want a strictly one-to-one table |
expand writes one row per copy, pairing copies by index across columns, so an
orthogroup costs rows equal to its largest cell. A cell with more genes than
--max-copies is read as a family and contributes nothing; the conversion
counts those. The track draws a gene pair once however many rows name it, so the
extra ribbons stay on the band the duplication is about.
At one locus the two ribbons are countable. Sorghum sits over rice as the control, sharing the grasses' ancestry without maize's duplication.
The genes with a single maize ribbon lost a copy after the duplication.
Making the ids resolve
The BEDs place each gene, and column 4 must match the table's ids byte for byte. OrthoFinder takes an id from the first token of the FASTA header, while a BED from the GFF3 is keyed on the gene, so the build script renames each protein to its gene id when it prepares the proteomes.
With --bed name=file per column, the script reports what share of each
column's ids that BED places and drops the rest. A column placing near none has
an id mismatch, and one placing none stops the conversion. The same output
counts the orthogroups that became several rows.
Loading the orthogroups in JBrowse
One track backs every band of the stack, the same as the MCScan blocks track:
{
"type": "SyntenyTrack",
"trackId": "wheat_orthogroups",
"name": "Wheat orthogroups (OrthoFinder)",
"assemblyNames": [
"tauschii",
"wheat",
"durum",
"emmer",
"urartu",
"timopheevii"
],
"adapter": {
"type": "MCScanBlocksAdapter",
"uri": "tauschii.blocks.gz",
"blockAssemblies": [
"durum",
"emmer",
"tauschii",
"timopheevii",
"urartu",
"wheat"
],
"bedLocations": [
{ "uri": "durum.bed.gz" },
{ "uri": "emmer.bed.gz" },
{ "uri": "tauschii.bed.gz" },
{ "uri": "timopheevii.bed.gz" },
{ "uri": "urartu.bed.gz" },
{ "uri": "wheat.bed.gz" }
],
"assemblyNames": [
"tauschii",
"wheat",
"durum",
"emmer",
"urartu",
"timopheevii"
]
}
}jbrowse add-track-json '{
"type": "SyntenyTrack",
"trackId": "wheat_orthogroups",
"name": "Wheat orthogroups (OrthoFinder)",
"assemblyNames": [
"tauschii",
"wheat",
"durum",
"emmer",
"urartu",
"timopheevii"
],
"adapter": {
"type": "MCScanBlocksAdapter",
"uri": "tauschii.blocks.gz",
"blockAssemblies": [
"durum",
"emmer",
"tauschii",
"timopheevii",
"urartu",
"wheat"
],
"bedLocations": [
{ "uri": "durum.bed.gz" },
{ "uri": "emmer.bed.gz" },
{ "uri": "tauschii.bed.gz" },
{ "uri": "timopheevii.bed.gz" },
{ "uri": "urartu.bed.gz" },
{ "uri": "wheat.bed.gz" }
],
"assemblyNames": [
"tauschii",
"wheat",
"durum",
"emmer",
"urartu",
"timopheevii"
]
}
}'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": "wheat_orthogroups",
"name": "Wheat orthogroups (OrthoFinder)",
"assemblyNames": [
"tauschii",
"wheat",
"durum",
"emmer",
"urartu",
"timopheevii"
],
"adapter": {
"type": "MCScanBlocksAdapter",
"uri": "tauschii.blocks.gz",
"blockAssemblies": [
"durum",
"emmer",
"tauschii",
"timopheevii",
"urartu",
"wheat"
],
"bedLocations": [
{ "uri": "durum.bed.gz" },
{ "uri": "emmer.bed.gz" },
{ "uri": "tauschii.bed.gz" },
{ "uri": "timopheevii.bed.gz" },
{ "uri": "urartu.bed.gz" },
{ "uri": "wheat.bed.gz" }
],
"assemblyNames": [
"tauschii",
"wheat",
"durum",
"emmer",
"urartu",
"timopheevii"
]
}
}blockAssemblies and bedLocations follow the table's columns (OrthoFinder's
proteome scan, alphabetical here), while assemblyNames is the order the stack
draws. Take the column order from what the conversion printed; a mismatch is
reported as a track error naming both lists.
An orthogroup is a set, so any two filled columns are a direct statement about
that pair and row order in the stack is free, unlike a
jcvi .blocks table
anchored on one column.
Assemblies without sequence
A gene-level synteny view never reads a base, so each assembly is a
ChromSizesAdapter built from the
##sequence-region header of its GFF3: a few kilobytes where the sequence is
tens of gigabytes.
jbrowse add-assembly wheat.chrom.sizes --name wheat --load copyEnsembl lists every unplaced scaffold in that header, so the script keeps only
the sequences carrying the most genes and prints what share of each genome's
genes they hold. Raise MAXSEQ where the sequences it would add are
chromosomes.
Reproduce it end to end
build_orthofinder_synteny.sh
runs everything above and writes a config.json with the assemblies, gene
tracks, the synteny track and a stacked default session.
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_orthofinder_synteny.sh
bash build_orthofinder_synteny.sh wheat # or: vertebrates, grasses, drosophila, solanaceae
npx --yes serve orthofinder_wheat_build/jbrowse2 # then open the printed URLThe sets it knows, and what each costs to build:
| Set | Genomes | DIAMOND runs | Annotation source |
|---|---|---|---|
vertebrates | human, chicken, frog, gar, zebrafish | 25 | Ensembl 113 |
grasses | rice, sorghum, maize, brachypodium, setaria | 25 | Ensembl Plants 63 |
wheat | tauschii, wheat, durum, emmer, urartu, timopheevii | 36 | Ensembl Plants 63 |
drosophila | melanogaster, simulans, yakuba, pseudoobscura, virilis | 25 | Ensembl Metazoa 63 |
solanaceae | tomato, potato, pepper, tobacco, coffee | 25 | Ensembl Plants 63 |
Two cuts are environment variables:
| Variable | Default | What it cuts |
|---|---|---|
MAXSEQ | 30 | sequence regions kept per genome, the ones carrying the most genes |
MAXCOPIES | 4 | genes in one orthogroup cell past which it is a gene family rather than a set of copies |
MAXSEQ=60 MAXCOPIES=6 bash build_orthofinder_synteny.sh wheatOrthoFinder searches every proteome against every other, so its DIAMOND count is the square of the set's size. Every step is guarded on its output file, so a re-run picks up where it stopped. Three sets need the NCBI datasets CLI to name chromosomes their GFF3 gives as INSDC accessions, each from its sequence report.
Your own genomes
A set name only tells the script which Ensembl files to download. The same run works on your own genomes given two files per genome, the FASTA and its GFF3. gffread translates each CDS and prints the transcript-to-gene map alongside, so the proteome and the gene rows come from one parse of one file, and reference names and lengths come from the FASTA index it writes.
Column 2 also takes a proteome, which saves translating one Ensembl already
publishes. Its headers then carry a gene:<id> tag that has to match the GFF3's
ID=gene:<id>. The run says which of the two it read, and prints the share of
ids it placed.
Name the files in a manifest, one line per genome:
cat > my_genomes.tsv <<'EOF'
# name genome annotation aliases
speciesA data/speciesA.fa.gz data/speciesA.gff3.gz
speciesB https://host/B.fa.gz https://host/B.gff3.gz GCF_000001405.40
EOF
bash build_orthofinder_synteny.sh my_genomes.tsv
npx --yes serve orthofinder_my_genomes_build/jbrowse2 # then open the printed URLColumn 1 names the assembly. The file columns take a local path or a URL, and two genomes make a valid manifest. Column 4 is optional: an INSDC assembly accession fetches NCBI's sequence report for the submitter's chromosome names, and anything else is read as a two-column alias table you supply.
See also
- Synteny from an ortholog table (grape, peach, cacao)
- Synteny from MCScan anchors (grape, peach)
- Synteny visualization (pairwise minimap2)
- Linear synteny view
- Synteny track
- MCScanBlocksAdapter
References
- Emms and Kelly (2019). OrthoFinder: phylogenetic orthology inference for comparative genomics
- Devos et al. (1995). Structural evolution of wheat chromosomes 4A, 5A, and 7B and its impact on recombination
- Dvorak et al. (2018). Reassessment of the evolution of wheat chromosomes 4A, 5A, and 7B
Feedback on this tutorial is welcome: contact us.