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 - A running JBrowse instance (the web quickstart or the desktop quickstart)
OrthoFinder bundles its own DIAMOND, so one install covers the first bullet:
bioconda has orthofinder, and the project
publishes a davidemms/orthofinder container with the same bundle. Everything
here calls it as orthofinder on PATH, so without root that container works
too, wrapped as a shim; the build script header has
the Apptainer version.
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 where those genes sit, so a table of orthogroups is a synteny track for a human against a zebrafish. Pairwise minimap2 aligns sequence to sequence and MCScan compares annotations over collinear runs of genes, and both have a divergence past which they 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, and then the second half of the page builds
one of them, wheat, starting from Orthogroups.tsv.
Vertebrates: blocks that survive out to zebrafish
The vertebrates set is human, chicken, frog, spotted gar and zebrafish, whose
common ancestor is a few hundred million years back and whose orthologs still
fall into chromosome-scale blocks. The teleost genome duplication shows up as
counting: a human chromosome answers to one or two chicken chromosomes, and to
more zebrafish ones.
Every band draws one line per ortholog, so it resolves into wedges only where a chromosome's orthologs mostly land on one chromosome of the row below; the build script prints that share for each adjacent pair.
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 this stack lands on the layout wheat figures are conventionally drawn in: each homoeologous group's chromosomes come together, so the hexaploid row reads as groups, with nothing telling the view which chromosomes are homoeologs.
The rows share one bp/px, which is Show all regions - same bp per pixel in
the view menu and sameScale in a session spec. A row's drawn length is then
its genome size: the diploid donors are drawn short and the hexaploid fills the
frame.
Reading one chromosome out of the stack
The track lists every genome in the set, so any pair of them opens as a two-row view with no second file. This one puts Aegilops tauschii's seven chromosomes over bread wheat 4A alone, with the palette button's Query painting each link by the tauschii chromosome it leaves.
The middle bundle leaves the right-hand end of 5D, close by 6D's tick: 6D reaches 4A only as scattered singletons, like the other three uninvolved chromosomes. Those four are the control the blocks are read against.
The blocks themselves are the 4AL/5AL and 4AL/7BS translocation pair, RFLP-mapped by Devos et al. and revisited against the reference assemblies by Dvorak et al. Those names are A and B genome chromosomes, and a group number is shared across the three subgenomes, so the 5D and 7D bundles in the frame are the D-genome counterparts of the 5A and 7B they refer to. The input is orthogroup membership and each gene's position.
Bread wheat 4A against Triticum urartu
Bread wheat 4A is an A-genome chromosome, and the row over it so far has been Aegilops tauschii, the D-genome donor. Triticum urartu, the A-genome donor, is the same two-row view with one assembly swapped.
Urartu's chromosome 4 covers both of the first two blocks, and the distal block is on its chromosome 7. Any of the other four assemblies in the track opens the same way.
Drosophila: chromosome arms that outlast gene order
The drosophila set is five fly genomes. D. simulans and D. yakuba sit
beside D. melanogaster; D. pseudoobscura and D. virilis are roughly 25 and
50 million years out. Flies keep their chromosome arms across that whole range
(Muller's elements, A through F) and rewrite the order of the genes inside them,
so the same orthogroup table answers two different questions depending on how
far you zoom.
The conversion prints the chromosome half as it runs, for each 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, that is a colour per melanogaster arm arriving as one bundle in every row, on a chromosome whose name changes as the lineages rename their own.
One locus, one lane per fly
The gene-order half needs a window, and a multi-way synteny track draws it in a single linear view: a lane per fly, each in its own coordinates, with the lane's header naming the chromosome that fly keeps these orthologs on.
{
"defaultSession": {
"name": "Drosophila multi-way synteny track",
"views": [
{
"type": "LinearGenomeView",
"init": {
"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",
"init": {
"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
}
]
}
}
]
}
EOFTwenty melanogaster genes on 3L, from Bre1 through S6k and mad2 to PXo,
and every one of the four other flies keeps all twenty. The two near relatives
keep them in order too, so their ribbons run parallel. The two distant ones keep
the block and reverse it, which the header says as [rev].
The pseudoobscura lane is the one to read twice. Muller element D is melanogaster's 3L, and in the obscura lineage that element is fused to the X, so a lane fitted to this window's orthologs sits at 59.8 Mb on a chromosome the assembly calls X. Nothing in the table knows that; the lane header is naming the chromosome its own placements landed on.
Nightshades: the same genes over four times the DNA
The solanaceae set is tomato, potato and pepper, Nicotiana attenuata as a
fourth nightshade, and coffee as the outgroup. Their gene counts are within a
factor of 1.5 of each other, 25,574 for coffee to 39,021 for potato, and their
genomes are not: 0.38 Gb of coffee against 2.9 Gb of pepper, most of that
difference being repeat sequence between the genes rather than genes.
Stacked on one bp/px, a row's drawn length is its genome size, so the stack states that difference before any ribbon is read.
N. attenuata is the assembly still on scaffolds here, and the correspondence print says so rather than leaving it to be discovered: its best partner holds 14% of a chromosome's links, against 77% for tomato to potato. Its own genes are spread over thousands of sequences, of which the build keeps the 30 densest, so the row draws the share that fell on those.
One locus, five lanes, five scales
The same table in a multi-way synteny track makes the size difference per-gene rather than per-genome. Each lane is fitted to the orthologs of the window in that genome's own coordinates and then says what scale that took.
{
"defaultSession": {
"name": "Nightshade multi-way synteny track",
"views": [
{
"type": "LinearGenomeView",
"init": {
"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",
"init": {
"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
}
]
}
}
]
}
EOFA 157 kb window on tomato chromosome 4 whose two dozen genes every one of the four other genomes keeps. Potato and coffee draw them at 1.5 times the anchor's span; pepper and N. attenuata need 3 times it for the same genes, which is the intergenic expansion arriving as a number in a lane header.
Every lane's genes stay in the anchor's order, so what changed between them is
the spacing rather than the arrangement. The coffee lane is [rev], the whole
block inverted in the outgroup.
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, and the two
bands it sits between draw visibly more ribbons per gene than the rest of the
stack.
The ribbon count is a conversion setting rather than a property of the genomes, and what to do with a duplicated gene is where the build picks it.
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.
OrthoFinder takes a directory of proteomes, one FASTA per genome, and -og
stops it after the orthogroups, which is all this table needs:
# -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,
naming the directory for the day it ran.
Orthogroups.tsv is already one row per orthogroup and one column per genome.
It needs the header row and the leading Orthogroup 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; nothing in the file makes that column special.
The script prints the column order it read off the header row, which is what
blockAssemblies has to be. A column is named after the proteome file
OrthoFinder read, with the extensions taken off, so a --bed key is the
assembly name only where the proteomes were named that way;
--assembly COLUMN=NAME renames the ones that were not, 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 from one gene to one gene, so a cell holding two genes has no single correct answer. The grasses are where this decides the picture: a rice gene commonly has two maize orthologs, one per copy of the duplication maize carries.
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 |
first is what a one-line awk reduction does, and its link is chosen by file
order. expand writes one row per copy, pairing copies by index across columns,
so an orthogroup costs rows equal to its largest cell and a gene family cannot
blow the table up. A cell with more genes than --max-copies is read as a
family and contributes nothing; the conversion counts those, so a threshold set
below the ploidy in the set shows up in that count.
Each of those rows carries the other genomes' genes too, so a pair that is not duplicated is named on every one of them. The track draws a gene pair once, however many rows name it, which keeps the extra ribbons on the band the duplication is about.
At one locus the default's two ribbons are countable. Sorghum sits over rice here as the control, since it shares the grasses' ancestry without maize's duplication, and the maize row carries both of the regions the duplication left.
The genes with a single maize ribbon lost a copy after the duplication.
Making the ids resolve
The table is coordinate-free. The BEDs place each gene, and column 4 must match the table's ids byte for byte. OrthoFinder takes a sequence's id from the first token of its FASTA header, so the ids are whatever the proteome headers led with, while a BED built from the GFF3 is likely keyed on the gene. The build script settles this by renaming each protein to its gene id when it prepares the proteomes, so both sides speak gene ids.
Pass --bed name=file per column and the script reports what share of each
column's ids that BED places, dropping the rest from the table. Read that share
before loading anything: a column placing near none of its ids has an id
mismatch, and one placing none of them stops the conversion. A column given no
--bed is reported as unchecked. The same output reports how many orthogroups
held a duplicated gene and became several rows.
Loading the orthogroups in JBrowse
One track backs every band of the stack, the same as the
MCScan blocks track,
and jbrowse add-track-json takes it as written, since neither
blockAssemblies nor bedLocations has an add-track flag:
{
"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 are positional against the table's own
columns, which follow OrthoFinder's proteome scan (alphabetical here), while the
track's assemblyNames is the order the stack draws them. Take the column order
from what the conversion printed; a mismatch reads another genome's BED for
every gene lookup, which the adapter reports as a track error naming both lists.
An orthogroup is a set, so any two columns filled on a row are a direct
statement about that pair, and a row with nothing in column 0 is kept like any
other. Row order in the stack is free, and any pair of the six opens as a
two-row view. A
jcvi .blocks table
is anchored on one of its columns.
Assemblies without sequence
A gene-level synteny view never reads a base, so each of these assemblies is a
ChromSizesAdapter built from the
##sequence-region header of its GFF3. The wheat lineage as sequence is tens of
gigabytes to host; as names and lengths it is a few kilobytes.
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. An ortholog on a sequence the assembly
leaves out draws nothing, so the build prints what share of each genome's genes
the kept sequences hold. Read it for a fragmented assembly, and raise MAXSEQ
where the sequences it would add are chromosomes.
Reproduce it end to end
build_orthofinder_synteny.sh
runs everything above. It downloads the proteomes and annotations from Ensembl,
reduces each proteome to one protein per gene, runs OrthoFinder, converts the
orthogroups, 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 |
The two cuts it makes are environment variables, so a set with a different karyotype or ploidy is handled from the command line.
| 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 wheatThe OrthoFinder step is the long one: it searches every proteome against every
other, so the DIAMOND count in the table above is the square of the set's size.
Everything 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 Ensembl GFF3
gives INSDC accessions instead: T. timopheevii in wheat, tomato in
solanaceae, and every fly but melanogaster in drosophila, each from its
sequence report.
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.