Skip to main content

Low-mappability regions (SMN)

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.

A pileup looks the same whether its reads belong at a locus or merely landed there. At the SMN1/SMN2 duplication on chromosome 5, we tell the two apart with tracks that genomes.jbrowse.org publishes for hg38 (reference mappability, gnomAD coverage and problematic-region annotations) and a public 1000 Genomes CRAM colored by mapping quality. We then ask whether the finished T2T-CHM13 reference resolves the duplication, and compare the locus against a control window from the same sample.

Prerequisites

Where the data comes from

The hosted tracks are files already wired into the genomes.jbrowse.org hg38 config. The reads come from the 1000 Genomes high-coverage short-read release and its ONT long-read release (Gustafson et al. 2024).

The build script takes these files from their URLs, so there is nothing to download by hand.

Read by URL (no download needed)

SMN1 and SMN2: a near-identical gene duplication

SMN1 and SMN2 sit about 900 kb apart on chromosome 5 and are roughly 99.9% identical across their ~28 kb. Spinal muscular atrophy turns on the copy number of SMN1, so which copy a read came from is the clinical question. An aligner given a 150 bp read from either copy has two equally good places to put it and reports mapping quality 0 (MAPQ 0). MAPQ is -10 log10 Pr{mapping position is wrong} (SAM specification), so MAPQ 0 means the chosen position is about as likely wrong as right.

Mappability, coverage and read tracks across the SMN1 and SMN2 region

Turn these tracks on from the track selector of the hosted hg38 config at genomes.jbrowse.org, or add them with the configs below to any JBrowse that has the hg38 assembly. UCSC publishes the mappability and problematic-region tracks for hg38 only:

  • Multi-read mappability - Umap M100, the fraction of overlapping 100-mers at each position that are unique in the genome, computed from the reference alone. The file omits positions with no unique 100-mer, so the track goes blank there.
    • Resolution → Bin summary → Minimum, with Show min and max unchecked, draws the worst position in each bin
    • Y axis... → Range pins the axis at 0 to 1
    • In the 2.5 Mb frame, at about a kilobase per pixel, even Minimum sits on the floor everywhere, so read this track in the narrower read view
  • gnomAD v3 Genome Coverage - Mean Coverage, averaged over tens of thousands of genomes. gnomAD drops reads that do not place uniquely before averaging, so this track falls wherever the Umap track is blank.
  • GIAB Problematic Regions - LowMap+SegDup, the sequence GIAB leaves out of its benchmark.
  • Long-read SVs - 1KG Vienna ONT SVs, the long-read callset.

The Umap config holds the Minimum score mode and the 0 to 1 axis:

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

{
  "type": "QuantitativeTrack",
  "trackId": "hg38-umap100Quantitative",
  "name": "Multi-read mappability - Umap M100",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "BigWigAdapter",
    "uri": "https://hgdownload.soe.ucsc.edu/gbdb/hg38/hoffmanMappability/k100.Umap.MultiTrackMappability.bw"
  },
  "displays": [
    {
      "type": "LinearWiggleDisplay",
      "displayId": "hg38-umap100Quantitative-LinearWiggleDisplay",
      "aggregate": "min",
      "extent": "none",
      "scales": { "y": { "domainMin": 0, "domainMax": 1 } }
    }
  ]
}
jbrowse add-track-json '{
  "type": "QuantitativeTrack",
  "trackId": "hg38-umap100Quantitative",
  "name": "Multi-read mappability - Umap M100",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "BigWigAdapter",
    "uri": "https://hgdownload.soe.ucsc.edu/gbdb/hg38/hoffmanMappability/k100.Umap.MultiTrackMappability.bw"
  },
  "displays": [
    {
      "type": "LinearWiggleDisplay",
      "displayId": "hg38-umap100Quantitative-LinearWiggleDisplay",
      "aggregate": "min",
      "extent": "none",
      "scales": { "y": { "domainMin": 0, "domainMax": 1 } }
    }
  ]
}'

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": "QuantitativeTrack",
  "trackId": "hg38-umap100Quantitative",
  "name": "Multi-read mappability - Umap M100",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "BigWigAdapter",
    "uri": "https://hgdownload.soe.ucsc.edu/gbdb/hg38/hoffmanMappability/k100.Umap.MultiTrackMappability.bw"
  },
  "displays": [
    {
      "type": "LinearWiggleDisplay",
      "displayId": "hg38-umap100Quantitative-LinearWiggleDisplay",
      "aggregate": "min",
      "extent": "none",
      "scales": { "y": { "domainMin": 0, "domainMax": 1 } }
    }
  ]
}

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

{
  "trackId": "hg38-gnomad3MeanCoverage",
  "name": "gnomAD v3 Genome Coverage - Mean Coverage",
  "uri": "https://hgdownload.soe.ucsc.edu/gbdb/hg38/gnomAD/coverage/v3-genome/gnomad.coverage.mean.bw",
  "assemblyNames": ["hg38"]
}
jbrowse add-track https://hgdownload.soe.ucsc.edu/gbdb/hg38/gnomAD/coverage/v3-genome/gnomad.coverage.mean.bw \
  --trackId hg38-gnomad3MeanCoverage \
  --name "gnomAD v3 Genome Coverage - Mean Coverage" \
  --assemblyNames hg38

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: https://hgdownload.soe.ucsc.edu/gbdb/hg38/gnomAD/coverage/v3-genome/gnomad.coverage.mean.bw

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

  • Track name: gnomAD v3 Genome Coverage - Mean Coverage
  • Assembly: hg38

Click Add.

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

{
  "trackId": "hg38-alllowmapandsegdupregions",
  "name": "GIAB Problematic Regions - LowMap+SegDup",
  "uri": "https://hgdownload.soe.ucsc.edu/gbdb/hg38/problematic/GIAB/alllowmapandsegdupregions.bb",
  "assemblyNames": ["hg38"]
}
jbrowse add-track https://hgdownload.soe.ucsc.edu/gbdb/hg38/problematic/GIAB/alllowmapandsegdupregions.bb \
  --trackId hg38-alllowmapandsegdupregions \
  --name "GIAB Problematic Regions - LowMap+SegDup" \
  --assemblyNames hg38

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: https://hgdownload.soe.ucsc.edu/gbdb/hg38/problematic/GIAB/alllowmapandsegdupregions.bb

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

  • Track name: GIAB Problematic Regions - LowMap+SegDup
  • Assembly: hg38

Click Add.

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

{
  "trackId": "hg38-lrSv1kgOnt",
  "name": "Long-read SVs - 1KG Vienna ONT SVs",
  "uri": "https://hgdownload.soe.ucsc.edu/gbdb/hg38/lrSv/1kgOnt.bb",
  "assemblyNames": ["hg38"]
}
jbrowse add-track https://hgdownload.soe.ucsc.edu/gbdb/hg38/lrSv/1kgOnt.bb \
  --trackId hg38-lrSv1kgOnt \
  --name "Long-read SVs - 1KG Vienna ONT SVs" \
  --assemblyNames hg38

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: https://hgdownload.soe.ucsc.edu/gbdb/hg38/lrSv/1kgOnt.bb

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

  • Track name: Long-read SVs - 1KG Vienna ONT SVs
  • Assembly: hg38

Click Add.

Then add the NA12878 reads. The CRAM needs its .crai beside it and decodes against hg38; for your own sample, swap the uri. In grammar-of-graphics terms, the track's color is a color scale on the mapq field:

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

{
  "trackId": "na12878_qc_reads",
  "name": "NA12878, 30x Illumina (1000 Genomes)",
  "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram",
  "assemblyNames": ["hg38"],
  "displayDefaults": { "color": { "field": "mapq" } }
}
jbrowse add-track https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram \
  --trackId na12878_qc_reads \
  --name "NA12878, 30x Illumina (1000 Genomes)" \
  --assemblyNames hg38 \
  --displayDefaults '{"color":{"field":"mapq"}}'

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:

{
  "trackId": "na12878_qc_reads",
  "name": "NA12878, 30x Illumina (1000 Genomes)",
  "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram",
  "assemblyNames": ["hg38"],
  "displayDefaults": { "color": { "field": "mapq" } }
}

Open this track in JBrowse Web ↗

Open this track in JBrowse Desktop ↗ (JBrowse Desktop 5.0+)

Open chr5:69,200,000-71,700,000 for the whole region, and a second view at chr5:70,850,000-71,500,000 for the reads.

Two scales of the same place. Top, the region on chr5 with SMN2 and SMN1 banded: RefSeq genes, gnomAD mean coverage, GIAB's low-mappability and segmental-duplication regions, and the 1000 Genomes long-read SV callset. Below it, a second view from SMN1 to where the reads recover, with Umap k100 mappability and NA12878 reads colored by mapping quality.
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%2Fucsc%2Fhg38%2Fconfig.json&session=spec-%7B%22views%22%3A%5B%7B%22type%22%3A%22LinearGenomeView%22%2C%22assembly%22%3A%22hg38%22%2C%22loc%22%3A%22chr5%3A69%2C200%2C000-71%2C700%2C000%22%2C%22highlight%22%3A%5B%7B%22refName%22%3A%22chr5%22%2C%22start%22%3A70924940%2C%22end%22%3A70953012%2C%22assemblyName%22%3A%22hg38%22%7D%2C%7B%22refName%22%3A%22chr5%22%2C%22start%22%3A70049523%2C%22end%22%3A70077595%2C%22assemblyName%22%3A%22hg38%22%7D%5D%2C%22tracks%22%3A%5B%7B%22trackId%22%3A%22hg38-ncbiRefSeqCurated%22%2C%22type%22%3A%22LinearBasicDisplay%22%2C%22height%22%3A60%2C%22showOnlyGenes%22%3Atrue%2C%22geneGlyphMode%22%3A%22longestCoding%22%2C%22displayMode%22%3A%22compact%22%7D%2C%7B%22trackId%22%3A%22hg38-gnomad3MeanCoverage%22%2C%22type%22%3A%22LinearWiggleDisplay%22%2C%22scales%22%3A%7B%22y%22%3A%7B%22domainMin%22%3A0%2C%22domainMax%22%3A40%7D%7D%2C%22height%22%3A100%2C%22extent%22%3A%22none%22%7D%2C%7B%22trackId%22%3A%22hg38-alllowmapandsegdupregions%22%2C%22type%22%3A%22LinearBasicDisplay%22%2C%22height%22%3A40%7D%2C%7B%22trackId%22%3A%22hg38-lrSv1kgOnt%22%2C%22type%22%3A%22LinearBasicDisplay%22%2C%22height%22%3A100%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+qc%2Fsmn_problematic_regions

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 hg38, loaded from a .2bit genome.
  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 one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “NCBI RefSeq - RefSeq Curated”.
    1. Track menu → Show only genes (checked)
    2. Track menu → Gene glyph → Representative transcript
    3. Track menu → Feature height → Compact
  3. Add the next track the same way. This one needs a BigWig file (.bw).This figure uses “gnomAD v3 Genome Coverage - Mean Coverage”.
    1. Track menu → Y axis... → RangePins the axis: Min 0 and Max 40. A field left empty follows the data in view.
    2. Track menu → Resolution → Bin summary → Show min and max (unchecked)
  4. Add the next track the same way. This one needs a BigBed file (.bb).This figure uses “GIAB Problematic Regions - LowMap+SegDup”.
  5. Add the next track the same way. This one needs a BigBed file (.bb).This figure uses “Long-read SVs - 1KG Vienna ONT SVs”.
  6. Type your region of interest into the location box and press Enter.This figure is at chr5:69,200,000-71,700,000.
  7. View menu → Zoom to region / use the location box, then add a highlightThis figure highlights 2 regions; highlights can also be set with the &highlight= URL parameter.

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

{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr5:69,200,000-71,700,000",
      "highlight": [
        {
          "refName": "chr5",
          "start": 70924940,
          "end": 70953012,
          "assemblyName": "hg38"
        },
        {
          "refName": "chr5",
          "start": 70049523,
          "end": 70077595,
          "assemblyName": "hg38"
        }
      ],
      "tracks": [
        {
          "trackId": "hg38-ncbiRefSeqCurated",
          "type": "LinearBasicDisplay",
          "height": 60,
          "showOnlyGenes": true,
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact"
        },
        {
          "trackId": "hg38-gnomad3MeanCoverage",
          "type": "LinearWiggleDisplay",
          "scales": {
            "y": {
              "domainMin": 0,
              "domainMax": 40
            }
          },
          "height": 100,
          "extent": "none"
        },
        {
          "trackId": "hg38-alllowmapandsegdupregions",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "hg38-lrSv1kgOnt",
          "type": "LinearBasicDisplay",
          "height": 100
        }
      ]
    }
  ]
}

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/ucsc/hg38/config.json")
app = JBrowseApp(
    assemblies=config["assemblies"],
    tracks=config.get("tracks", []),
    views=[
        {
            "type": "LinearGenomeView",
            "assembly": "hg38",
            "loc": "chr5:69,200,000-71,700,000",
            "highlight": [
                {
                    "refName": "chr5",
                    "start": 70924940,
                    "end": 70953012,
                    "assemblyName": "hg38"
                },
                {
                    "refName": "chr5",
                    "start": 70049523,
                    "end": 70077595,
                    "assemblyName": "hg38"
                }
            ],
            "tracks": [
                {
                    "trackId": "hg38-ncbiRefSeqCurated",
                    "type": "LinearBasicDisplay",
                    "height": 60,
                    "showOnlyGenes": True,
                    "geneGlyphMode": "longestCoding",
                    "displayMode": "compact"
                },
                {
                    "trackId": "hg38-gnomad3MeanCoverage",
                    "type": "LinearWiggleDisplay",
                    "scales": {
                        "y": {
                            "domainMin": 0,
                            "domainMax": 40
                        }
                    },
                    "height": 100,
                    "extent": "none"
                },
                {
                    "trackId": "hg38-alllowmapandsegdupregions",
                    "type": "LinearBasicDisplay",
                    "height": 40
                },
                {
                    "trackId": "hg38-lrSv1kgOnt",
                    "type": "LinearBasicDisplay",
                    "height": 100
                }
            ]
        }
    ],
)
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/ucsc/hg38/config.json \
  --assembly hg38 --loc chr5:69,200,000-71,700,000 \
  --track hg38-ncbiRefSeqCurated display:LinearBasicDisplay '{"height":60,"showOnlyGenes":true,"geneGlyphMode":"longestCoding","displayMode":"compact"}' \
  --track hg38-gnomad3MeanCoverage display:LinearWiggleDisplay '{"scales":{"y":{"domainMin":0,"domainMax":40}},"height":100,"extent":"none"}' \
  --track hg38-alllowmapandsegdupregions display:LinearBasicDisplay '{"height":40}' \
  --track hg38-lrSv1kgOnt display:LinearBasicDisplay '{"height":100}' \
  --width 1500 --out figure.svg

It has no flag for this figure's highlight, so its picture can differ there.

--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": "hg38",
      "loc": "chr5:69,200,000-71,700,000",
      "highlight": [
        {
          "refName": "chr5",
          "start": 70924940,
          "end": 70953012,
          "assemblyName": "hg38"
        },
        {
          "refName": "chr5",
          "start": 70049523,
          "end": 70077595,
          "assemblyName": "hg38"
        }
      ],
      "tracks": [
        {
          "trackId": "hg38-ncbiRefSeqCurated",
          "type": "LinearBasicDisplay",
          "height": 60,
          "showOnlyGenes": true,
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact"
        },
        {
          "trackId": "hg38-gnomad3MeanCoverage",
          "type": "LinearWiggleDisplay",
          "scales": {
            "y": {
              "domainMin": 0,
              "domainMax": 40
            }
          },
          "height": 100,
          "extent": "none"
        },
        {
          "trackId": "hg38-alllowmapandsegdupregions",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "hg38-lrSv1kgOnt",
          "type": "LinearBasicDisplay",
          "height": 100
        }
      ]
    }
  ]
}
JSON

cat > callouts.json <<'JSON'
[
  {
    "type": "text",
    "text": "the dip is mappability: reads that fit in two places are dropped",
    "fontSize": 19,
    "maxWidth": 340,
    "anchor": {
      "trackId": "hg38-gnomad3MeanCoverage",
      "loc": "chr5:70,400,000",
      "fracY": 0.18
    }
  },
  {
    "type": "text",
    "text": "even with long reads, very few calls in this block",
    "fontSize": 18,
    "maxWidth": 340,
    "anchor": {
      "trackId": "hg38-lrSv1kgOnt",
      "loc": "chr5:70,300,000",
      "fracY": 0.55
    }
  }
]
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/ucsc/hg38/config.json \
  --width 1500 --height 622 \
  --spec session.json --annotations callouts.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.

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%2Fucsc%2Fhg38%2Fconfig.json&session=spec-%7B%22sessionTracks%22%3A%5B%7B%22type%22%3A%22AlignmentsTrack%22%2C%22trackId%22%3A%22na12878_qc_reads%22%2C%22name%22%3A%22NA12878%2C+30x+Illumina+%281000+Genomes%29%22%2C%22assemblyNames%22%3A%5B%22hg38%22%5D%2C%22adapter%22%3A%7B%22type%22%3A%22CramAdapter%22%2C%22cramLocation%22%3A%7B%22uri%22%3A%22https%3A%2F%2Fs3.amazonaws.com%2F1000genomes%2F1000G_2504_high_coverage%2Fdata%2FERR3239334%2FNA12878.final.cram%22%7D%2C%22craiLocation%22%3A%7B%22uri%22%3A%22https%3A%2F%2Fs3.amazonaws.com%2F1000genomes%2F1000G_2504_high_coverage%2Fdata%2FERR3239334%2FNA12878.final.cram.crai%22%7D%7D%7D%5D%2C%22views%22%3A%5B%7B%22type%22%3A%22LinearGenomeView%22%2C%22assembly%22%3A%22hg38%22%2C%22loc%22%3A%22chr5%3A70%2C850%2C000-71%2C500%2C000%22%2C%22highlight%22%3A%5B%7B%22refName%22%3A%22chr5%22%2C%22start%22%3A70924940%2C%22end%22%3A70953012%2C%22assemblyName%22%3A%22hg38%22%7D%5D%2C%22tracks%22%3A%5B%7B%22trackId%22%3A%22hg38-ncbiRefSeqCurated%22%2C%22type%22%3A%22LinearBasicDisplay%22%2C%22height%22%3A60%2C%22showOnlyGenes%22%3Atrue%2C%22geneGlyphMode%22%3A%22longestCoding%22%7D%2C%7B%22trackId%22%3A%22hg38-umap100Quantitative%22%2C%22type%22%3A%22LinearWiggleDisplay%22%2C%22scales%22%3A%7B%22y%22%3A%7B%22domainMin%22%3A0%2C%22domainMax%22%3A1%7D%7D%2C%22aggregate%22%3A%22min%22%2C%22extent%22%3A%22none%22%2C%22height%22%3A70%7D%2C%7B%22trackId%22%3A%22hg38-gnomad3MeanCoverage%22%2C%22type%22%3A%22LinearWiggleDisplay%22%2C%22scales%22%3A%7B%22y%22%3A%7B%22domainMin%22%3A0%2C%22domainMax%22%3A40%7D%7D%2C%22height%22%3A90%2C%22extent%22%3A%22none%22%7D%2C%7B%22trackId%22%3A%22na12878_qc_reads%22%2C%22type%22%3A%22LinearAlignmentsDisplay%22%2C%22color%22%3A%7B%22field%22%3A%22mapq%22%7D%2C%22showLegend%22%3Atrue%2C%22forceLoad%22%3Atrue%2C%22scales%22%3A%7B%22y%22%3A%7B%22domainQuantile%22%3A0.99%7D%7D%2C%22height%22%3A260%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+qc%2Fsmn_read_placement

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 hg38, loaded from a .2bit genome.
  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 one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “NCBI RefSeq - RefSeq Curated”.
    1. Track menu → Show only genes (checked)
    2. Track menu → Gene glyph → Representative transcript
  3. Add the next track the same way. This one needs a BigWig file (.bw).This figure uses “Multi-read mappability - Umap M100”.
    1. Track menu → Y axis... → RangePins the axis: Min 0 and Max 1. A field left empty follows the data in view.
    2. Track menu → Resolution → Bin summary → Minimum
    3. Track menu → Resolution → Bin summary → Show min and max (unchecked)
  4. Add the next track the same way. This one needs a BigWig file (.bw).This figure uses “gnomAD v3 Genome Coverage - Mean Coverage”.
    1. Track menu → Y axis... → RangePins the axis: Min 0 and Max 40. A field left empty follows the data in view.
    2. Track menu → Resolution → Bin summary → Show min and max (unchecked)
  5. Add the next track the same way. This one needs a CRAM file (.cram + .crai).This figure uses “NA12878, 30x Illumina (1000 Genomes)”.
    1. Track menu → Color by... → Mapping quality
    2. Track menu → Show... → Show legend (checked)
    3. Click Force load in the track's "Zoom in to see features or force load" message.Loads the region even past the byte-size limit, which can be slow.
    4. Track menu → Coverage axis... → Clip extreme outliers (checked)
  6. Type your region of interest into the location box and press Enter.This figure is at chr5:70,850,000-71,500,000.
  7. View menu → Zoom to region / use the location box, then add a highlightThis figure highlights 1 region; highlights can also be set with the &highlight= URL parameter.

A figure adds its tracks to one session. The jbrowse CLI writes the same tracks into a config.json instead, where every session that opens it has them: add-track where flags cover the whole track, add-track-json where they do not.

jbrowse add-track-json '{
  "type": "AlignmentsTrack",
  "trackId": "na12878_qc_reads",
  "name": "NA12878, 30x Illumina (1000 Genomes)",
  "assemblyNames": [
    "hg38"
  ],
  "adapter": {
    "type": "CramAdapter",
    "cramLocation": {
      "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram"
    },
    "craiLocation": {
      "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram.crai"
    }
  }
}'

Run these where the config.json is, or add --out <dir>, and point each uri at your own file. The config needs hg38 already — jbrowse add-assembly adds yours. The location and the settings the steps carry are session state rather than track config — that half is the Spec tab, or jbrowse set-default-session.

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

{
  "sessionTracks": [
    {
      "type": "AlignmentsTrack",
      "trackId": "na12878_qc_reads",
      "name": "NA12878, 30x Illumina (1000 Genomes)",
      "assemblyNames": [
        "hg38"
      ],
      "adapter": {
        "type": "CramAdapter",
        "cramLocation": {
          "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram"
        },
        "craiLocation": {
          "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram.crai"
        }
      }
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr5:70,850,000-71,500,000",
      "highlight": [
        {
          "refName": "chr5",
          "start": 70924940,
          "end": 70953012,
          "assemblyName": "hg38"
        }
      ],
      "tracks": [
        {
          "trackId": "hg38-ncbiRefSeqCurated",
          "type": "LinearBasicDisplay",
          "height": 60,
          "showOnlyGenes": true,
          "geneGlyphMode": "longestCoding"
        },
        {
          "trackId": "hg38-umap100Quantitative",
          "type": "LinearWiggleDisplay",
          "scales": {
            "y": {
              "domainMin": 0,
              "domainMax": 1
            }
          },
          "aggregate": "min",
          "extent": "none",
          "height": 70
        },
        {
          "trackId": "hg38-gnomad3MeanCoverage",
          "type": "LinearWiggleDisplay",
          "scales": {
            "y": {
              "domainMin": 0,
              "domainMax": 40
            }
          },
          "height": 90,
          "extent": "none"
        },
        {
          "trackId": "na12878_qc_reads",
          "type": "LinearAlignmentsDisplay",
          "color": {
            "field": "mapq"
          },
          "showLegend": true,
          "forceLoad": true,
          "scales": {
            "y": {
              "domainQuantile": 0.99
            }
          },
          "height": 260
        }
      ]
    }
  ]
}

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/ucsc/hg38/config.json")
app = JBrowseApp(
    assemblies=config["assemblies"],
    tracks=config.get("tracks", []) + [
        {
            "type": "AlignmentsTrack",
            "trackId": "na12878_qc_reads",
            "name": "NA12878, 30x Illumina (1000 Genomes)",
            "assemblyNames": [
                "hg38"
            ],
            "adapter": {
                "type": "CramAdapter",
                "cramLocation": {
                    "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram"
                },
                "craiLocation": {
                    "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram.crai"
                }
            }
        }
    ],
    views=[
        {
            "type": "LinearGenomeView",
            "assembly": "hg38",
            "loc": "chr5:70,850,000-71,500,000",
            "highlight": [
                {
                    "refName": "chr5",
                    "start": 70924940,
                    "end": 70953012,
                    "assemblyName": "hg38"
                }
            ],
            "tracks": [
                {
                    "trackId": "hg38-ncbiRefSeqCurated",
                    "type": "LinearBasicDisplay",
                    "height": 60,
                    "showOnlyGenes": True,
                    "geneGlyphMode": "longestCoding"
                },
                {
                    "trackId": "hg38-umap100Quantitative",
                    "type": "LinearWiggleDisplay",
                    "scales": {
                        "y": {
                            "domainMin": 0,
                            "domainMax": 1
                        }
                    },
                    "aggregate": "min",
                    "extent": "none",
                    "height": 70
                },
                {
                    "trackId": "hg38-gnomad3MeanCoverage",
                    "type": "LinearWiggleDisplay",
                    "scales": {
                        "y": {
                            "domainMin": 0,
                            "domainMax": 40
                        }
                    },
                    "height": 90,
                    "extent": "none"
                },
                {
                    "trackId": "na12878_qc_reads",
                    "type": "LinearAlignmentsDisplay",
                    "color": {
                        "field": "mapq"
                    },
                    "showLegend": True,
                    "forceLoad": True,
                    "scales": {
                        "y": {
                            "domainQuantile": 0.99
                        }
                    },
                    "height": 260
                }
            ]
        }
    ],
)
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.

cat > tracks.json <<'JSON'
[
  {
    "type": "AlignmentsTrack",
    "trackId": "na12878_qc_reads",
    "name": "NA12878, 30x Illumina (1000 Genomes)",
    "assemblyNames": [
      "hg38"
    ],
    "adapter": {
      "type": "CramAdapter",
      "cramLocation": {
        "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram"
      },
      "craiLocation": {
        "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram.crai"
      }
    }
  }
]
JSON

npx @jbrowse/img --config https://jbrowse.org/ucsc/hg38/config.json \
  --tracks tracks.json \
  --assembly hg38 --loc chr5:70,850,000-71,500,000 \
  --track hg38-ncbiRefSeqCurated display:LinearBasicDisplay '{"height":60,"showOnlyGenes":true,"geneGlyphMode":"longestCoding"}' \
  --track hg38-umap100Quantitative display:LinearWiggleDisplay '{"scales":{"y":{"domainMin":0,"domainMax":1}},"aggregate":"min","extent":"none","height":70}' \
  --track hg38-gnomad3MeanCoverage display:LinearWiggleDisplay '{"scales":{"y":{"domainMin":0,"domainMax":40}},"height":90,"extent":"none"}' \
  --track na12878_qc_reads display:LinearAlignmentsDisplay '{"color":{"field":"mapq"},"showLegend":true,"forceLoad":true,"scales":{"y":{"domainQuantile":0.99}},"height":260}' \
  --width 1500 --out figure.svg

It has no flag for this figure's highlight, so its picture can differ there.

--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'
{
  "sessionTracks": [
    {
      "type": "AlignmentsTrack",
      "trackId": "na12878_qc_reads",
      "name": "NA12878, 30x Illumina (1000 Genomes)",
      "assemblyNames": [
        "hg38"
      ],
      "adapter": {
        "type": "CramAdapter",
        "cramLocation": {
          "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram"
        },
        "craiLocation": {
          "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram.crai"
        }
      }
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr5:70,850,000-71,500,000",
      "highlight": [
        {
          "refName": "chr5",
          "start": 70924940,
          "end": 70953012,
          "assemblyName": "hg38"
        }
      ],
      "tracks": [
        {
          "trackId": "hg38-ncbiRefSeqCurated",
          "type": "LinearBasicDisplay",
          "height": 60,
          "showOnlyGenes": true,
          "geneGlyphMode": "longestCoding"
        },
        {
          "trackId": "hg38-umap100Quantitative",
          "type": "LinearWiggleDisplay",
          "scales": {
            "y": {
              "domainMin": 0,
              "domainMax": 1
            }
          },
          "aggregate": "min",
          "extent": "none",
          "height": 70
        },
        {
          "trackId": "hg38-gnomad3MeanCoverage",
          "type": "LinearWiggleDisplay",
          "scales": {
            "y": {
              "domainMin": 0,
              "domainMax": 40
            }
          },
          "height": 90,
          "extent": "none"
        },
        {
          "trackId": "na12878_qc_reads",
          "type": "LinearAlignmentsDisplay",
          "color": {
            "field": "mapq"
          },
          "showLegend": true,
          "forceLoad": true,
          "scales": {
            "y": {
              "domainQuantile": 0.99
            }
          },
          "height": 260
        }
      ]
    }
  ]
}
JSON

cat > callouts.json <<'JSON'
[
  {
    "type": "text",
    "text": "SMN1: biallelic loss causes spinal muscular atrophy",
    "fontSize": 18,
    "maxWidth": 330,
    "anchor": {
      "trackId": "hg38-gnomad3MeanCoverage",
      "loc": "chr5:70,925,000",
      "alignX": "left",
      "fracY": 0.15
    }
  }
]
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/ucsc/hg38/config.json \
  --width 1500 --height 805 \
  --timeout 600000 \
  --spec session.json --annotations callouts.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.

GIAB flags two long intervals, much larger than the gene, and the gnomAD track stays low across both. Beyond them, GIAB flags nothing larger than a few kilobases for megabases in either direction. In the read view, the reads stay at MAPQ 0 until well past the end of SMN1, where the Umap track steps up at the same coordinate as the gnomAD coverage.

Zoom the read view to the SMN cassette, chr5:70,889,000-70,989,000.

The SMN cassette, holding SERF1A, SMN1 and NAIP, with the same four tracks and one read per row. Almost every read is dark blue, mapped where it is drawn and fitting somewhere else just as well.
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%2Fucsc%2Fhg38%2Fconfig.json&session=spec-%7B%22sessionTracks%22%3A%5B%7B%22type%22%3A%22AlignmentsTrack%22%2C%22trackId%22%3A%22na12878_qc_reads%22%2C%22name%22%3A%22NA12878%2C+30x+Illumina+%281000+Genomes%29%22%2C%22assemblyNames%22%3A%5B%22hg38%22%5D%2C%22adapter%22%3A%7B%22type%22%3A%22CramAdapter%22%2C%22cramLocation%22%3A%7B%22uri%22%3A%22https%3A%2F%2Fs3.amazonaws.com%2F1000genomes%2F1000G_2504_high_coverage%2Fdata%2FERR3239334%2FNA12878.final.cram%22%7D%2C%22craiLocation%22%3A%7B%22uri%22%3A%22https%3A%2F%2Fs3.amazonaws.com%2F1000genomes%2F1000G_2504_high_coverage%2Fdata%2FERR3239334%2FNA12878.final.cram.crai%22%7D%7D%7D%5D%2C%22views%22%3A%5B%7B%22type%22%3A%22LinearGenomeView%22%2C%22assembly%22%3A%22hg38%22%2C%22loc%22%3A%22chr5%3A70%2C889%2C000-70%2C989%2C000%22%2C%22tracks%22%3A%5B%7B%22trackId%22%3A%22hg38-ncbiRefSeqCurated%22%2C%22type%22%3A%22LinearBasicDisplay%22%2C%22height%22%3A70%2C%22showOnlyGenes%22%3Afalse%2C%22geneGlyphMode%22%3A%22longestCoding%22%7D%2C%7B%22trackId%22%3A%22hg38-umap100Quantitative%22%2C%22type%22%3A%22LinearWiggleDisplay%22%2C%22scales%22%3A%7B%22y%22%3A%7B%22domainMin%22%3A0%2C%22domainMax%22%3A1%7D%7D%2C%22aggregate%22%3A%22min%22%2C%22extent%22%3A%22none%22%2C%22height%22%3A70%7D%2C%7B%22trackId%22%3A%22hg38-gnomad3MeanCoverage%22%2C%22type%22%3A%22LinearWiggleDisplay%22%2C%22scales%22%3A%7B%22y%22%3A%7B%22domainMin%22%3A0%2C%22domainMax%22%3A40%7D%7D%2C%22height%22%3A70%7D%2C%7B%22trackId%22%3A%22na12878_qc_reads%22%2C%22type%22%3A%22LinearAlignmentsDisplay%22%2C%22color%22%3A%7B%22field%22%3A%22mapq%22%7D%2C%22showLegend%22%3Atrue%2C%22featureHeight%22%3A2%2C%22height%22%3A300%2C%22forceLoad%22%3Atrue%2C%22scales%22%3A%7B%22y%22%3A%7B%22domainQuantile%22%3A0.99%7D%7D%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+qc%2Fsmn1_evidence

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 hg38, loaded from a .2bit genome.
  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 one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “NCBI RefSeq - RefSeq Curated”.
    1. Track menu → Show only genes (unchecked)
    2. Track menu → Gene glyph → Representative transcript
  3. Add the next track the same way. This one needs a BigWig file (.bw).This figure uses “Multi-read mappability - Umap M100”.
    1. Track menu → Y axis... → RangePins the axis: Min 0 and Max 1. A field left empty follows the data in view.
    2. Track menu → Resolution → Bin summary → Minimum
    3. Track menu → Resolution → Bin summary → Show min and max (unchecked)
  4. Add the next track the same way. This one needs a BigWig file (.bw).This figure uses “gnomAD v3 Genome Coverage - Mean Coverage”.
    1. Track menu → Y axis... → RangePins the axis: Min 0 and Max 40. A field left empty follows the data in view.
  5. Add the next track the same way. This one needs a CRAM file (.cram + .crai).This figure uses “NA12878, 30x Illumina (1000 Genomes)”.
    1. Track menu → Color by... → Mapping quality
    2. Track menu → Show... → Show legend (checked)
    3. Track menu → Read height → Custom... → 2px
    4. Click Force load in the track's "Zoom in to see features or force load" message.Loads the region even past the byte-size limit, which can be slow.
    5. Track menu → Coverage axis... → Clip extreme outliers (checked)
  6. Type your region of interest into the location box and press Enter.This figure is at chr5:70,889,000-70,989,000.

A figure adds its tracks to one session. The jbrowse CLI writes the same tracks into a config.json instead, where every session that opens it has them: add-track where flags cover the whole track, add-track-json where they do not.

jbrowse add-track-json '{
  "type": "AlignmentsTrack",
  "trackId": "na12878_qc_reads",
  "name": "NA12878, 30x Illumina (1000 Genomes)",
  "assemblyNames": [
    "hg38"
  ],
  "adapter": {
    "type": "CramAdapter",
    "cramLocation": {
      "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram"
    },
    "craiLocation": {
      "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram.crai"
    }
  }
}'

Run these where the config.json is, or add --out <dir>, and point each uri at your own file. The config needs hg38 already — jbrowse add-assembly adds yours. The location and the settings the steps carry are session state rather than track config — that half is the Spec tab, or jbrowse set-default-session.

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

{
  "sessionTracks": [
    {
      "type": "AlignmentsTrack",
      "trackId": "na12878_qc_reads",
      "name": "NA12878, 30x Illumina (1000 Genomes)",
      "assemblyNames": [
        "hg38"
      ],
      "adapter": {
        "type": "CramAdapter",
        "cramLocation": {
          "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram"
        },
        "craiLocation": {
          "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram.crai"
        }
      }
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr5:70,889,000-70,989,000",
      "tracks": [
        {
          "trackId": "hg38-ncbiRefSeqCurated",
          "type": "LinearBasicDisplay",
          "height": 70,
          "showOnlyGenes": false,
          "geneGlyphMode": "longestCoding"
        },
        {
          "trackId": "hg38-umap100Quantitative",
          "type": "LinearWiggleDisplay",
          "scales": {
            "y": {
              "domainMin": 0,
              "domainMax": 1
            }
          },
          "aggregate": "min",
          "extent": "none",
          "height": 70
        },
        {
          "trackId": "hg38-gnomad3MeanCoverage",
          "type": "LinearWiggleDisplay",
          "scales": {
            "y": {
              "domainMin": 0,
              "domainMax": 40
            }
          },
          "height": 70
        },
        {
          "trackId": "na12878_qc_reads",
          "type": "LinearAlignmentsDisplay",
          "color": {
            "field": "mapq"
          },
          "showLegend": true,
          "featureHeight": 2,
          "height": 300,
          "forceLoad": true,
          "scales": {
            "y": {
              "domainQuantile": 0.99
            }
          }
        }
      ]
    }
  ]
}

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/ucsc/hg38/config.json")
app = JBrowseApp(
    assemblies=config["assemblies"],
    tracks=config.get("tracks", []) + [
        {
            "type": "AlignmentsTrack",
            "trackId": "na12878_qc_reads",
            "name": "NA12878, 30x Illumina (1000 Genomes)",
            "assemblyNames": [
                "hg38"
            ],
            "adapter": {
                "type": "CramAdapter",
                "cramLocation": {
                    "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram"
                },
                "craiLocation": {
                    "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram.crai"
                }
            }
        }
    ],
    views=[
        {
            "type": "LinearGenomeView",
            "assembly": "hg38",
            "loc": "chr5:70,889,000-70,989,000",
            "tracks": [
                {
                    "trackId": "hg38-ncbiRefSeqCurated",
                    "type": "LinearBasicDisplay",
                    "height": 70,
                    "showOnlyGenes": False,
                    "geneGlyphMode": "longestCoding"
                },
                {
                    "trackId": "hg38-umap100Quantitative",
                    "type": "LinearWiggleDisplay",
                    "scales": {
                        "y": {
                            "domainMin": 0,
                            "domainMax": 1
                        }
                    },
                    "aggregate": "min",
                    "extent": "none",
                    "height": 70
                },
                {
                    "trackId": "hg38-gnomad3MeanCoverage",
                    "type": "LinearWiggleDisplay",
                    "scales": {
                        "y": {
                            "domainMin": 0,
                            "domainMax": 40
                        }
                    },
                    "height": 70
                },
                {
                    "trackId": "na12878_qc_reads",
                    "type": "LinearAlignmentsDisplay",
                    "color": {
                        "field": "mapq"
                    },
                    "showLegend": True,
                    "featureHeight": 2,
                    "height": 300,
                    "forceLoad": True,
                    "scales": {
                        "y": {
                            "domainQuantile": 0.99
                        }
                    }
                }
            ]
        }
    ],
)
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.

cat > tracks.json <<'JSON'
[
  {
    "type": "AlignmentsTrack",
    "trackId": "na12878_qc_reads",
    "name": "NA12878, 30x Illumina (1000 Genomes)",
    "assemblyNames": [
      "hg38"
    ],
    "adapter": {
      "type": "CramAdapter",
      "cramLocation": {
        "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram"
      },
      "craiLocation": {
        "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram.crai"
      }
    }
  }
]
JSON

npx @jbrowse/img --config https://jbrowse.org/ucsc/hg38/config.json \
  --tracks tracks.json \
  --assembly hg38 --loc chr5:70,889,000-70,989,000 \
  --track hg38-ncbiRefSeqCurated display:LinearBasicDisplay '{"height":70,"showOnlyGenes":false,"geneGlyphMode":"longestCoding"}' \
  --track hg38-umap100Quantitative display:LinearWiggleDisplay '{"scales":{"y":{"domainMin":0,"domainMax":1}},"aggregate":"min","extent":"none","height":70}' \
  --track hg38-gnomad3MeanCoverage display:LinearWiggleDisplay '{"scales":{"y":{"domainMin":0,"domainMax":40}},"height":70}' \
  --track na12878_qc_reads display:LinearAlignmentsDisplay '{"color":{"field":"mapq"},"showLegend":true,"featureHeight":2,"height":300,"forceLoad":true,"scales":{"y":{"domainQuantile":0.99}}}' \
  --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'
{
  "sessionTracks": [
    {
      "type": "AlignmentsTrack",
      "trackId": "na12878_qc_reads",
      "name": "NA12878, 30x Illumina (1000 Genomes)",
      "assemblyNames": [
        "hg38"
      ],
      "adapter": {
        "type": "CramAdapter",
        "cramLocation": {
          "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram"
        },
        "craiLocation": {
          "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram.crai"
        }
      }
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr5:70,889,000-70,989,000",
      "tracks": [
        {
          "trackId": "hg38-ncbiRefSeqCurated",
          "type": "LinearBasicDisplay",
          "height": 70,
          "showOnlyGenes": false,
          "geneGlyphMode": "longestCoding"
        },
        {
          "trackId": "hg38-umap100Quantitative",
          "type": "LinearWiggleDisplay",
          "scales": {
            "y": {
              "domainMin": 0,
              "domainMax": 1
            }
          },
          "aggregate": "min",
          "extent": "none",
          "height": 70
        },
        {
          "trackId": "hg38-gnomad3MeanCoverage",
          "type": "LinearWiggleDisplay",
          "scales": {
            "y": {
              "domainMin": 0,
              "domainMax": 40
            }
          },
          "height": 70
        },
        {
          "trackId": "na12878_qc_reads",
          "type": "LinearAlignmentsDisplay",
          "color": {
            "field": "mapq"
          },
          "showLegend": true,
          "featureHeight": 2,
          "height": 300,
          "forceLoad": true,
          "scales": {
            "y": {
              "domainQuantile": 0.99
            }
          }
        }
      ]
    }
  ]
}
JSON

cat > callouts.json <<'JSON'
[
  {
    "type": "text",
    "text": "Dark blue is MAPQ 0: mapped here, but they fit elsewhere just as well",
    "fontSize": 20,
    "maxWidth": 430,
    "anchor": {
      "trackId": "na12878_qc_reads",
      "loc": "chr5:70,895,000",
      "fracY": 0.22
    }
  }
]
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/ucsc/hg38/config.json \
  --width 1500 --height 790 \
  --spec session.json --annotations callouts.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.

Dark blue is MAPQ 0; yellow is MAPQ 30 and above. Advanced → Edit plot... on the reads track shows this color as text and applies edits live.

Does T2T-CHM13 resolve the SMN duplication?

T2T-CHM13, the telomere-to-telomere assembly, is finished on this chromosome. UCSC's hg38-to-CHM13 liftOver forms several long, overlapping chains over the region, some of them reversed:

tabix https://jbrowse.org/ucsc/hg38/liftOver/hg38ToHs1.over.pif.gz \
  tchr5:69200000-71700000

A synteny view between the two assemblies draws each chain as a ribbon. The hosted hg38 config has the track and loads hs1 with it. On another JBrowse, load the hs1 assembly and the chain file, whose query is the genome it lifts to:

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

{
  "name": "hs1",
  "uri": "https://hgdownload.soe.ucsc.edu/goldenPath/hs1/bigZips/hs1.2bit"
}
jbrowse add-assembly https://hgdownload.soe.ucsc.edu/goldenPath/hs1/bigZips/hs1.2bit \
  --name hs1

In JBrowse Desktop, Open new genome on the start screen (or File → Open genome... in a session), then Open from a URL and paste:

https://hgdownload.soe.ucsc.edu/goldenPath/hs1/bigZips/hs1.2bit

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

  • Genome name: hs1

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

{
  "type": "SyntenyTrack",
  "trackId": "hg38_to_hs1_liftOver",
  "name": "hg38 to Human (hs1) liftOver",
  "assemblyNames": ["hs1", "hg38"],
  "adapter": {
    "type": "PairwiseIndexedPAFAdapter",
    "uri": "https://jbrowse.org/ucsc/hg38/liftOver/hg38ToHs1.over.pif.gz",
    "csi": true,
    "queryAssembly": "hs1",
    "targetAssembly": "hg38"
  }
}
jbrowse add-track-json '{
  "type": "SyntenyTrack",
  "trackId": "hg38_to_hs1_liftOver",
  "name": "hg38 to Human (hs1) liftOver",
  "assemblyNames": ["hs1", "hg38"],
  "adapter": {
    "type": "PairwiseIndexedPAFAdapter",
    "uri": "https://jbrowse.org/ucsc/hg38/liftOver/hg38ToHs1.over.pif.gz",
    "csi": true,
    "queryAssembly": "hs1",
    "targetAssembly": "hg38"
  }
}'

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": "hg38_to_hs1_liftOver",
  "name": "hg38 to Human (hs1) liftOver",
  "assemblyNames": ["hs1", "hg38"],
  "adapter": {
    "type": "PairwiseIndexedPAFAdapter",
    "uri": "https://jbrowse.org/ucsc/hg38/liftOver/hg38ToHs1.over.pif.gz",
    "csi": true,
    "queryAssembly": "hs1",
    "targetAssembly": "hg38"
  }
}

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.

GRCh38 above, T2T-CHM13 below, each framed on that assembly's SMN2-to-SMN1 span and its flanks, ribbons from UCSC's liftOver chains and colored by strand. Three chains cross each other.
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=test_data%2Fhg38_hs1_synteny%2Fconfig.json&session=spec-%7B%22views%22%3A%5B%7B%22type%22%3A%22LinearSyntenyView%22%2C%22color%22%3A%7B%22field%22%3A%22strand%22%7D%2C%22drawCurves%22%3Atrue%2C%22cigarMode%22%3A%22matches%22%2C%22minAlignmentLength%22%3A10000%2C%22opacity%22%3A0.5%2C%22levelHeights%22%3A%5B300%5D%2C%22tracks%22%3A%5B%5B%22hg38_hs1_synteny%22%5D%5D%2C%22views%22%3A%5B%7B%22assembly%22%3A%22hg38%22%2C%22loc%22%3A%22chr5%3A69%2C450%2C000-71%2C550%2C000%22%2C%22tracks%22%3A%5B%7B%22trackId%22%3A%22hg38-genes%22%2C%22showOnlyGenes%22%3Atrue%2C%22geneGlyphMode%22%3A%22longestCoding%22%2C%22heightMode%22%3A%22fixed%22%2C%22height%22%3A120%7D%5D%2C%22trackLabels%22%3A%22offset%22%7D%2C%7B%22assembly%22%3A%22hs1%22%2C%22loc%22%3A%22chr5%3A70%2C210%2C000-71%2C980%2C000%22%2C%22tracks%22%3A%5B%7B%22trackId%22%3A%22hs1-genes%22%2C%22showOnlyGenes%22%3Atrue%2C%22geneGlyphMode%22%3A%22longestCoding%22%2C%22showLabels%22%3A%22name%22%2C%22heightMode%22%3A%22fixed%22%2C%22height%22%3A200%7D%5D%2C%22trackLabels%22%3A%22offset%22%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+qc%2Fsmn_vs_t2t

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 genomes: on the JBrowse Desktop start screen click Open new genome (or Show all available genomes to pick a hosted one), then File → Open genome... for each of the rest. For JBrowse Web, run jbrowse add-assembly once for each to put them in your config.json.This figure's rows, top to bottom: hg38, hs1.
  2. Open the view: Add → Linear synteny view, and set the rows to those assemblies in the import form.
  3. Point the import form at your own file: click the arrow between two rows, choose New track, and paste a URL or pick a local file. This one needs a bgzip-indexed PAF (.pif.gz + .tbi).A synteny track the session already has is quicker: pick it under Quick start, which fills in the rows itself.This figure uses “hg38 vs T2T-CHM13 (UCSC liftOver)”.
  4. Click Launch.
  5. Synteny view header → palette button → StrandThe palette button's tooltip reads "Color by: ...".
  6. Synteny display settings (the sliders button in the view header) → Curved lines (checked)
  7. Synteny display settings (the sliders button in the view header) → CIGAR indels → Transparent indelsThe row appears only when the alignments actually carry CIGAR ops — a coarse-tier PIF or a CIGAR-less PAF has none to draw.
  8. Synteny display settings (the sliders button in the view header) → Min length → drag to 10,000bpHides alignments shorter than this, which is what clears the hairball of short spurious chains at whole-genome zoom.
  9. Synteny display settings (the sliders button in the view header) → Opacity → drag to 0.5Lower values let dense overlapping ribbons show through each other.
  10. Drag the bar below a synteny level to resize it.This figure's level is 300px tall.
  11. Row 1 (hg38)
    1. Add a track: open this view's track selector and check one your config has, or click + and choose Add track to paste a URL or choose a local file. This one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “NCBI RefSeq genes”.
      1. Track menu → Show only genes (checked)
      2. Track menu → Gene glyph → Representative transcript
      3. Track menu → Feature height → Track sizing → Fixed feature height + fixed track height
    2. Type your region of interest into the location box and press Enter.This figure is at chr5:69,450,000-71,550,000.
    3. View menu → Show... → Offset
  12. Row 2 (hs1)
    1. Add a track: open this view's track selector and check one your config has, or click + and choose Add track to paste a URL or choose a local file. This one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “NCBI RefSeq genes”.
      1. Track menu → Show only genes (checked)
      2. Track menu → Gene glyph → Representative transcript
      3. Track menu → Show... → Labels → Name only
      4. Track menu → Feature height → Track sizing → Fixed feature height + fixed track height
    2. Type your region of interest into the location box and press Enter.This figure is at chr5:70,210,000-71,980,000.
    3. View menu → Show... → Offset

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

{
  "views": [
    {
      "type": "LinearSyntenyView",
      "color": {
        "field": "strand"
      },
      "drawCurves": true,
      "cigarMode": "matches",
      "minAlignmentLength": 10000,
      "opacity": 0.5,
      "levelHeights": [
        300
      ],
      "tracks": [
        [
          "hg38_hs1_synteny"
        ]
      ],
      "views": [
        {
          "assembly": "hg38",
          "loc": "chr5:69,450,000-71,550,000",
          "tracks": [
            {
              "trackId": "hg38-genes",
              "showOnlyGenes": true,
              "geneGlyphMode": "longestCoding",
              "heightMode": "fixed",
              "height": 120
            }
          ],
          "trackLabels": "offset"
        },
        {
          "assembly": "hs1",
          "loc": "chr5:70,210,000-71,980,000",
          "tracks": [
            {
              "trackId": "hs1-genes",
              "showOnlyGenes": true,
              "geneGlyphMode": "longestCoding",
              "showLabels": "name",
              "heightMode": "fixed",
              "height": 200
            }
          ],
          "trackLabels": "offset"
        }
      ]
    }
  ]
}

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/code/jb2/main/test_data/hg38_hs1_synteny/config.json")
app = JBrowseApp(
    assemblies=config["assemblies"],
    tracks=config.get("tracks", []),
    views=[
        {
            "type": "LinearSyntenyView",
            "color": {
                "field": "strand"
            },
            "drawCurves": True,
            "cigarMode": "matches",
            "minAlignmentLength": 10000,
            "opacity": 0.5,
            "levelHeights": [
                300
            ],
            "tracks": [
                [
                    "hg38_hs1_synteny"
                ]
            ],
            "views": [
                {
                    "assembly": "hg38",
                    "loc": "chr5:69,450,000-71,550,000",
                    "tracks": [
                        {
                            "trackId": "hg38-genes",
                            "showOnlyGenes": True,
                            "geneGlyphMode": "longestCoding",
                            "heightMode": "fixed",
                            "height": 120
                        }
                    ],
                    "trackLabels": "offset"
                },
                {
                    "assembly": "hs1",
                    "loc": "chr5:70,210,000-71,980,000",
                    "tracks": [
                        {
                            "trackId": "hs1-genes",
                            "showOnlyGenes": True,
                            "geneGlyphMode": "longestCoding",
                            "showLabels": "name",
                            "heightMode": "fixed",
                            "height": 200
                        }
                    ],
                    "trackLabels": "offset"
                }
            ]
        }
    ],
)
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.

cat > session.json <<'JSON'
{
  "views": [
    {
      "type": "LinearSyntenyView",
      "color": {
        "field": "strand"
      },
      "drawCurves": true,
      "cigarMode": "matches",
      "minAlignmentLength": 10000,
      "opacity": 0.5,
      "levelHeights": [
        300
      ],
      "tracks": [
        [
          "hg38_hs1_synteny"
        ]
      ],
      "views": [
        {
          "assembly": "hg38",
          "loc": "chr5:69,450,000-71,550,000",
          "tracks": [
            {
              "trackId": "hg38-genes",
              "showOnlyGenes": true,
              "geneGlyphMode": "longestCoding",
              "heightMode": "fixed",
              "height": 120
            }
          ],
          "trackLabels": "offset"
        },
        {
          "assembly": "hs1",
          "loc": "chr5:70,210,000-71,980,000",
          "tracks": [
            {
              "trackId": "hs1-genes",
              "showOnlyGenes": true,
              "geneGlyphMode": "longestCoding",
              "showLabels": "name",
              "heightMode": "fixed",
              "height": 200
            }
          ],
          "trackLabels": "offset"
        }
      ]
    }
  ]
}
JSON

npx @jbrowse/img synteny --config https://jbrowse.org/code/jb2/main/test_data/hg38_hs1_synteny/config.json \
  --spec session.json \
  --width 1200 --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": "LinearSyntenyView",
      "color": {
        "field": "strand"
      },
      "drawCurves": true,
      "cigarMode": "matches",
      "minAlignmentLength": 10000,
      "opacity": 0.5,
      "levelHeights": [
        300
      ],
      "tracks": [
        [
          "hg38_hs1_synteny"
        ]
      ],
      "views": [
        {
          "assembly": "hg38",
          "loc": "chr5:69,450,000-71,550,000",
          "tracks": [
            {
              "trackId": "hg38-genes",
              "showOnlyGenes": true,
              "geneGlyphMode": "longestCoding",
              "heightMode": "fixed",
              "height": 120
            }
          ],
          "trackLabels": "offset"
        },
        {
          "assembly": "hs1",
          "loc": "chr5:70,210,000-71,980,000",
          "tracks": [
            {
              "trackId": "hs1-genes",
              "showOnlyGenes": true,
              "geneGlyphMode": "longestCoding",
              "showLabels": "name",
              "heightMode": "fixed",
              "height": 200
            }
          ],
          "trackLabels": "offset"
        }
      ]
    }
  ]
}
JSON

cat > callouts.json <<'JSON'
[
  {
    "type": "text",
    "fontSize": 20,
    "maxWidth": 420,
    "anchor": {
      "selector": "[data-testid=\"synteny_canvas\"][data-display-drawn=\"true\"]",
      "alignX": "left",
      "alignY": "top",
      "dx": 24,
      "dy": 30
    },
    "text": "The same GRCh38 sequence chains to more than one place in T2T"
  }
]
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/code/jb2/main/test_data/hg38_hs1_synteny/config.json \
  --width 1200 --height 900 \
  --timeout 180000 \
  --spec session.json --annotations callouts.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 gene order is the same in both assemblies, SMN2 first and then SMN1. The crossing chains join each GRCh38 copy to both CHM13 copies, so CHM13 has the same duplication, with the two genes closer together.

Long reads test whether T2T-CHM13 places reads at SMN1 better than GRCh38 does. The 1000 Genomes ONT release aligned GM18501 to both references with one minimap2 pipeline, and scan_mappability_qc.sh prints the share of its records over SMN1 at MAPQ 0 and MAPQ 60 on each. The long reads place better than the short reads at the same gene, and the MAPQ 0 share is the same on both references.

Depth and MAPQ at SMN1 and at a control window

scan_mappability_qc.sh counts the reads in equal windows over SMN1 and over the right-hand end of the read view, from the same library, and finds the same depth in both. A coverage track with no MAPQ filter draws flat across both. The MAPQ 0 share separates them, with most reads at SMN1 at MAPQ 0 and almost none at the control; the gnomAD depth over SMN1 falls to a fraction of the control's for the same reason.

Long-read SV calls across the SMN1 and SMN2 region

The long-read SV track in the wide view is empty across the region. scan_mappability_qc.sh counts calls over the flagged region and an equal-width window on either side, and finds few inside and many on both sides, where the older Database of Genomic Variants (DGV) merged catalogue has records throughout.

Across the whole chromosome, both catalogues put a larger share of their calls inside the flagged regions than those regions' share of chr5. Segmental duplications are copy-number variable, so real variation and artifacts both concentrate there.

Applying the mappability check to another locus

For another hg38 locus, turn on the four tracks above with the same Umap settings, add your reads colored by mapping quality, and put a second window of the same width, from the same sample and outside every flagged interval, beside the first.

In numbers, the comparison is three samtools counts per window, with -q a minimum MAPQ:

samtools view -c "$CRAM" chr5:70,900,000-71,000,000            # every read
samtools view -c -q 1 "$CRAM" chr5:70,900,000-71,000,000       # placed at all
samtools view -c -q 60 "$CRAM" chr5:70,900,000-71,000,000      # placed uniquely

The Umap track's value over a window is one command over the same bigWig:

# positions with no unique 100-mer are absent from the file, so the share of
# the span that has any value is the number to read
bigWigToBedGraph -chrom=chr5 -start=70049000 -end=70077000 \
  k100.Umap.MultiTrackMappability.bw stdout |
  awk -v s=70049000 -v e=70077000 \
    '{cov += $3 - $2} END {printf "%.1f%% has a value\n", 100 * cov / (e - s)}'

Run both on the control window too, at the same width.

Reproduce it end to end

scan_mappability_qc.sh produces every number on this page from the files the figures draw. It asks whether a read can be placed at a locus several ways, and reads each answer against a control:

  1. Over SMN1, SMN2 and two controls, it measures how much of each window the Umap track has a value for, gnomAD's mean depth, and the share of NA12878's reads at MAPQ 0. The near control is a window of the same width at the 5' end of BDP1, past both flagged intervals on the same chromosome, so sample, library and chromosome all match; ACTB on chromosome 7 shows the near control is ordinary.
  2. It lists which problematic-region tracks flag the region, and bins the gnomAD track across it to find where the depth recovers.
  3. It counts GM18501's long-read records over SMN1 at MAPQ 0 and MAPQ 60, once on each reference. The release aligned that sample to both with the same pipeline, so the two counts compare.
  4. It counts each SV catalogue's calls inside the region and in a flank of about the same width on either side, and across chr5 the share of calls whose midpoint falls in a flagged region.

It needs kent tools (bigWigInfo, bigWigToBedGraph, bigBedToBed), bedtools, samtools, curl and awk.

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

The script prints its sections in this page's order, and measures a locus added to its LOCI list the same way.

See also

Citations