Synteny from MCScan anchors (grape, peach)
TL;DR: a pairwise jcvi MCScan run
writes two files that JBrowse loads as separate synteny tracks: .anchors (one
orthologous gene pair per line, via MCScanAnchorsAdapter) and
.anchors.simple (one synteny block per line, via
MCScanSimpleAnchorsAdapter). Both pair genes by name, so each also needs a BED
per genome mapping gene ids to coordinates.
Prerequisites
- jcvi with the LAST aligner
- Or, in place of jcvi, an existing MCScanX run: converting one needs only python3
samtools- htslib (
bgzip,tabix) wgetnode, for the JBrowse CLI- A running JBrowse instance (the web quickstart or the desktop quickstart)
On Debian/Ubuntu, apt install samtools tabix wget last-align covers the
aligner and the file tools; jcvi installs with pip install jcvi and node
comes from nodejs.org.
Where the data comes from
Grape (Jaillon et al. 2007) and peach genomes and gene annotations, Ensembl Plants release 58.
- grape (Vitis_vinifera, PN40024.v4) genome and CDS FASTA: http://ftp.ensemblgenomes.org/pub/plants/release-58/fasta/vitis_vinifera/
- grape gene annotation: http://ftp.ensemblgenomes.org/pub/plants/release-58/gff3/vitis_vinifera/
- peach (Prunus_persica, Prunus_persica_NCBIv2) genome and CDS FASTA: http://ftp.ensemblgenomes.org/pub/plants/release-58/fasta/prunus_persica/
- peach gene annotation: http://ftp.ensemblgenomes.org/pub/plants/release-58/gff3/prunus_persica/
What MCScan compares
MCScan works through the two genomes' gene annotations, which finds synteny between species too divergent to line up base by base with minimap2. An anchor is a gene pair, so there is no CIGAR and nothing to draw below the level of a gene.
For three or more genomes from one MCScan run, see
ortholog tables, which
loads a .blocks table with one track backing every band.
What .anchors and .anchors.simple hold
.anchors is the gene-pair level. Each line is one orthologous pair and its
alignment score, with ### separating synteny blocks:
###
VIT_201s0011g00070.1 Prupe.1G290900.1 1430
VIT_201s0011g00080.1 Prupe.1G290800.1 446
VIT_201s0011g00090.1 Prupe.1G290700.1 147.anchors.simple is the same run reduced to one line per block: the first and
last gene of the block on each side, a score, and the block's orientation:
VIT_201s0011g00070.1 VIT_201s0011g00910.1 Prupe.1G281700.1 Prupe.1G290900.1 149 -
VIT_201s0011g02000.1 VIT_201s0011g02280.2 Prupe.1G345900.1 Prupe.1G348100.1 53 -
VIT_201s0011g02300.1 VIT_201s0011g02530.1 Prupe.1G299800.1 Prupe.1G303200.1 39 +The adapter turns those four gene ids into one feature spanning the block on
each genome, so .anchors.simple draws one ribbon per block where .anchors
draws one per gene pair. Neither file carries coordinates; the BED files supply
them.
BED files
One BED per genome, prepared from its GFF3 before the ortholog run. Only the first six columns are read, and column 4 must match the anchor gene ids byte for byte:
chr1 12836 26777 VIT_201s0011g00010.1 0 +
chr1 33170 35791 VIT_201s0011g00030.1 0 +Column 1 must use the same reference sequence names as the JBrowse assembly.
Which mismatches are loud and which are silent is the
adapters' own gotcha.
The one that bites this pipeline is jcvi stripping isoform suffixes unless run
with --no_strip_names, which the script passes:
strip them on one side only and no row resolves.
Producing the data
One jcvi command writes both anchor files and the BEDs are prepared from each GFF3 beforehand:
python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id \
--primary_only grape.gff3.gz -o grape.bed
python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id \
--primary_only peach.gff3.gz -o peach.bed
python -m jcvi.formats.fasta format grape.cds.fa.gz grape.cds
python -m jcvi.formats.fasta format peach.cds.fa.gz peach.cds
python -m jcvi.compara.catalog ortholog --no_strip_names grape peachThat leaves grape.peach.anchors and grape.peach.anchors.simple in the
working directory. The adapters read anchors and BED files plain or gzipped.
Loading both tracks
Each adapter takes the anchor file plus the two BEDs, and assemblyNames lists
the genomes in the order the anchor columns are in (column 1's genome first):
{
"type": "SyntenyTrack",
"trackId": "grape_peach_anchors",
"name": "Grape peach synteny (MCScan, anchors)",
"assemblyNames": ["grape", "peach"],
"adapter": {
"type": "MCScanAnchorsAdapter",
"uri": "grape.peach.anchors.gz",
"bed1": "grape.bed.gz",
"bed2": "peach.bed.gz",
"assemblyNames": ["grape", "peach"]
}
}jbrowse add-track grape.peach.anchors.gz \
--trackId grape_peach_anchors \
--name "Grape peach synteny (MCScan, anchors)" \
--assemblyNames grape,peach \
--bed1 grape.bed.gz \
--bed2 peach.bed.gz \
--load copyThe simple-anchors track is the same shape with the adapter type and file swapped:
{
"type": "SyntenyTrack",
"trackId": "grape_peach_anchors_simple",
"name": "Grape peach synteny (MCScan, simple anchors)",
"assemblyNames": ["grape", "peach"],
"adapter": {
"type": "MCScanSimpleAnchorsAdapter",
"uri": "grape.peach.anchors.simple.gz",
"bed1": "grape.bed.gz",
"bed2": "peach.bed.gz",
"assemblyNames": ["grape", "peach"]
}
}jbrowse add-track grape.peach.anchors.simple.gz \
--trackId grape_peach_anchors_simple \
--name "Grape peach synteny (MCScan, simple anchors)" \
--assemblyNames grape,peach \
--bed1 grape.bed.gz \
--bed2 peach.bed.gz \
--load copyBoth BEDs are jbrowse add-track flags, which is the CLI tab on each block
above: --bed1 and --bed2 beside the anchors file, and --load copy copies
all three into the config's directory.
Both adapters read the whole file into memory, which suits MCScan's scale. PAF's indexed PIF has no MCScan equivalent.
Both tracks in one synteny view
The two tracks describe the same run at different granularity. Add a linear synteny view (Add → Linear synteny view), pick peach and grape, and turn on both.
The block track is drawn here as an LGVSyntenyDisplay: a synteny track in an
ordinary linear genome view row, drawn as features. No add-track flag covers
the displays array, so this goes in with jbrowse add-track-json, which
copies no data files.
{
"type": "SyntenyTrack",
"trackId": "grape_peach_anchors_simple",
"name": "Grape peach synteny (MCScan, simple anchors)",
"assemblyNames": ["grape", "peach"],
"adapter": {
"type": "MCScanSimpleAnchorsAdapter",
"uri": "grape.peach.anchors.simple.gz",
"bed1": "grape.bed.gz",
"bed2": "peach.bed.gz",
"assemblyNames": ["grape", "peach"]
},
"displays": [
{
"type": "LGVSyntenyDisplay",
"displayId": "grape_peach_anchors_simple-LGVSyntenyDisplay",
"height": 60
}
]
}Read the two together: a bar says a block is there and which way round it runs, and the ribbons above it say whether the genes inside hold their order.
What an anchor looks like up close
Zoom to one block with both gene tracks on and set to Show only genes, and the unit the file is made of is on screen.
MCScan anchors the genes it could pair confidently, so most genes in these windows carry no ribbon. Zooming further widens the ribbons, since the file says only which gene pairs with which.
The same anchors as a dotplot
Either track also loads in a dotplot (Add → Dotplot view), where a gene pair is one point and a block is a run of them. The axes start in each assembly's index order, which scatters the runs; Re-order chromosomes sorts the vertical axis to follow the horizontal one, using the alignments themselves.
Reordering puts each peach chromosome's strongest grape partner on the diagonal,
and its other partners stay off it. The script reads
that off .anchors.simple directly, printing per peach chromosome the grape
chromosomes it shares blocks with and the anchors each pairing carries. Both
genomes descend from the ancestral eudicot hexaploidy (Jaillon et al.) and have
rearranged differently since.
Coming from MCScanX
MCScanX is a different program from jcvi's
MCScan and neither adapter reads its output as it stands: it writes one
.collinearity holding every block it found, self-synteny and cross-species
together, telling the genomes apart only by the two-letter tag it requires on
each chromosome name.
mcscanx_to_anchors.py
splits a run into the same four files jcvi writes, in place of the
jcvi step:
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/mcscanx_to_anchors.py
python3 mcscanx_to_anchors.py --gff xyz.gff --collinearity xyz.collinearity \
--species vv=grape --species pp=peach --strand-gff3 peach=peach.gff3.gzThat writes the four files the track configs above load
unchanged. --species order is the anchors column order, so it has to match the
track's assemblyNames. Two options decide whether the result draws:
--chr-prefix peach=Pp0prepends to the refNames, and--keep-chr-tagleaves MCScanX's tag on them. The tag is stripped by default (vv1becomes1).--strand-gff3 peach=peach.gff3.gzrecovers strand from the annotation the MCScanX input came from, which is what draws an.anchorspair as inverted. Block orientation comes from the collinearity block header.
--fai peach=peach.fa.fai checks the refNames against the assembly; a name the
assembly does not have draws empty in the browser. Scores are converted too, an
anchors score becoming -log10 of MCScanX's e-value where jcvi writes a bit
score.
Naming a third --species writes an ortholog table instead, since one
.collinearity covers every pair. See
ortholog tables.
A genome against itself
MCScanX is as often run on one genome to find its own duplicated blocks. Name a
single --species and the script keeps those blocks:
python3 mcscanx_to_anchors.py --gff grape.gff --collinearity grape.collinearity \
--species vv=grape --strand-gff3 grape=grape.gff3.gzThe anchor files that writes name grape on both sides, so the track lists the assembly twice and both rows of the synteny view are the same genome:
{
"type": "SyntenyTrack",
"trackId": "grape_self_anchors",
"name": "Grape duplicated blocks (MCScanX)",
"assemblyNames": ["grape", "grape"],
"adapter": {
"type": "MCScanAnchorsAdapter",
"uri": "grape.grape.anchors",
"bed1": "grape.bed",
"bed2": "grape.bed",
"assemblyNames": ["grape", "grape"]
}
}jbrowse add-track grape.grape.anchors \
--trackId grape_self_anchors \
--name "Grape duplicated blocks (MCScanX)" \
--assemblyNames grape,grape \
--bed1 grape.bed \
--bed2 grape.bed \
--load copyBoth copies of a block are served, so either one draws its link to the other. A dotplot of the track puts the genome on both axes, each duplicated block a run of points away from the diagonal, and no diagonal itself, since a gene is not its own anchor.
Reproduce it end to end
build_grape_peach_anchors.sh
runs everything above in one shot: the grape and peach genomes from Ensembl
Plants, the jcvi ortholog pipeline, and a config.json with both assemblies,
gene tracks, both MCScan tracks and a default session opening them together.
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_grape_peach_anchors.sh
bash build_grape_peach_anchors.sh
npx --yes serve grape_peach_anchors_build/jbrowse2 # then open the printed URLIts gene ids differ from the samples above, which come from a Phytozome annotation of the same two genomes; the ids are whatever the GFF3 carried.
See also
- Synteny visualization (pairwise minimap2)
- Synteny from an ortholog table (grape, peach, cacao)
- Linear synteny view
- Synteny track
- MCScanAnchorsAdapter
- MCScanSimpleAnchorsAdapter
References
- Tang et al. (2008). Unraveling ancient hexaploidy through multiply-aligned angiosperm gene maps, the MCScan method jcvi implements
- Jaillon et al. (2007). The grapevine genome sequence suggests ancestral hexaploidization in major angiosperm phyla
Feedback on this tutorial is welcome: contact us.