Skip to main content

A grammar of graphics over a BAM (NA12878 insert size)

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.11. Desktop beta builds are coming soon.

A read pair that straddles a deletion maps with a long insert (the distance between its reads), and a heterozygous deletion halves the read depth. We plot those two fields straight from NA12878's reads to find a deletion on chromosome 20 without a variant caller, scan the whole chromosome for the same signature, and check the hits against the 1000 Genomes callset. We draw the plots with JBrowse's mark display, a grammar of graphics over a track: each entry in marks names a mark type, a transform list and an encoding that maps feature fields to channels, as in the Alu tutorial. The mark display is experimental, and its config shape may change.

Prerequisites

  • a JBrowse to open the figures' sessions in (Web or Desktop)
  • samtools and htslib (bgzip, tabix), for cutting the pairs out of the file and for checking a window by hand
  • bcftools, for reading the callset at the end
  • Node.js and the JBrowse CLI, for the build script

Where the data comes from

1000 Genomes high-coverage release (Byrska-Bishop et al. 2022), GRCh38:

Loading hg38

The CRAM decodes against the assembly the track is added to, so the assembly must be the GRCh38 sequence the reads were aligned to. The cytoband table draws each chromosome's banding in the view's overview.

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

{
  "name": "hg38",
  "uri": "https://jbrowse.org/genomes/GRCh38/fasta/hg38.prefix.fa.gz",
  "refNameAliases": {
    "uri": "https://s3.amazonaws.com/jbrowse.org/genomes/GRCh38/hg38_aliases.txt"
  },
  "cytobands": "https://jbrowse.org/genomes/GRCh38/cytoBand.txt"
}
jbrowse add-assembly https://jbrowse.org/genomes/GRCh38/fasta/hg38.prefix.fa.gz \
  --name hg38 \
  --refNameAliases https://s3.amazonaws.com/jbrowse.org/genomes/GRCh38/hg38_aliases.txt \
  --config '{"cytobands":"https://jbrowse.org/genomes/GRCh38/cytoBand.txt"}'

Plotting read depth as bars over the EFCAB8 deletion

The window covers 30 kb of an EFCAB8 intron on chromosome 20, where the callset says NA12878 has one copy of a 3.9 kb deletion. A bar mark over a coverage transform draws the reads as runs of constant depth.

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

{
  "type": "AlignmentsTrack",
  "trackId": "na12878_read_depth",
  "name": "NA12878 depth (1000 Genomes, 30x)",
  "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram",
  "assemblyNames": ["hg38"],
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "scales": { "y": { "title": "Read depth" } },
      "marks": [
        {
          "mark": "bar",
          "transform": [{ "type": "coverage" }],
          "encoding": { "color": { "value": "#c8d8ee" } }
        }
      ]
    }
  ]
}
jbrowse add-track-json '{
  "type": "AlignmentsTrack",
  "trackId": "na12878_read_depth",
  "name": "NA12878 depth (1000 Genomes, 30x)",
  "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram",
  "assemblyNames": ["hg38"],
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "scales": { "y": { "title": "Read depth" } },
      "marks": [
        {
          "mark": "bar",
          "transform": [{ "type": "coverage" }],
          "encoding": { "color": { "value": "#c8d8ee" } }
        }
      ]
    }
  ]
}'

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": "AlignmentsTrack",
  "trackId": "na12878_read_depth",
  "name": "NA12878 depth (1000 Genomes, 30x)",
  "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram",
  "assemblyNames": ["hg38"],
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "scales": { "y": { "title": "Read depth" } },
      "marks": [
        {
          "mark": "bar",
          "transform": [{ "type": "coverage" }],
          "encoding": { "color": { "value": "#c8d8ee" } }
        }
      ]
    }
  ]
}

Open the track at chr20:32,925,000-32,955,000.

Thirty kilobases of an EFCAB8 intron in NA12878, the read depth as bars. Over a few kilobases in the middle of the window the depth runs at about half of what it is on either side.
Open in JBrowse Web ↗Open in JBrowse Desktop ↗

This session loads a large dataset and takes a while to open.

Or rebuild it from your own data:

Opening in Desktop needs JBrowse Desktop 5.0+, and saves the view as a reopenable session.

Nothing happens?

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%2Fdemos%2Fread_marks%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%22chr20%3A32%2C925%2C000-32%2C955%2C000%22%2C%22tracks%22%3A%5B%7B%22trackId%22%3A%22ncbi_refseq_hg38%22%2C%22type%22%3A%22LinearBasicDisplay%22%2C%22height%22%3A50%2C%22showOnlyGenes%22%3Atrue%2C%22geneGlyphMode%22%3A%22longestCoding%22%7D%2C%7B%22trackId%22%3A%22na12878_read_depth%22%2C%22type%22%3A%22LinearMarkDisplay%22%2C%22height%22%3A160%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+read_marks%2Fdepth

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).This figure uses hg38, loaded from a bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  2. Add your own track: File → Open track..., then 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. Drag the bar at the bottom edge of the track to resize it (50px here).
    2. Track menu → Show only genes (checked)
    3. Track menu → Gene glyph → Representative transcript
  3. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a CRAM file (.cram + .crai).This figure uses “NA12878 depth (1000 Genomes, 30x)”.
    1. Drag the bar at the bottom edge of the track to resize it (160px here).
  4. Type your region of interest into the location box and press Enter.This figure is at chr20:32,925,000-32,955,000.

Its genome and tracks are hosted, so File → Open track... there adds your own files beside them.

With your own data, the steps behind the figure:

  1. Put your genome in your JBrowse Web's config.json: jbrowse add-assembly adds it.This figure uses hg38, loaded from a bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  2. Add your own track: File → Open track..., then 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. Drag the bar at the bottom edge of the track to resize it (50px here).
    2. Track menu → Show only genes (checked)
    3. Track menu → Gene glyph → Representative transcript
  3. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a CRAM file (.cram + .crai).This figure uses “NA12878 depth (1000 Genomes, 30x)”.
    1. Drag the bar at the bottom edge of the track to resize it (160px here).
  4. Type your region of interest into the location box and press Enter.This figure is at chr20:32,925,000-32,955,000.

This session spec draws the figure from https://jbrowse.org/demos/read_marks/config.json: paste it after &session=spec- on a JBrowse Web link that loads that config.

{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr20:32,925,000-32,955,000",
      "tracks": [
        {
          "trackId": "ncbi_refseq_hg38",
          "type": "LinearBasicDisplay",
          "height": 50,
          "showOnlyGenes": true,
          "geneGlyphMode": "longestCoding"
        },
        {
          "trackId": "na12878_read_depth",
          "type": "LinearMarkDisplay",
          "height": 160
        }
      ]
    }
  ]
}

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/demos/read_marks/config.json")
app = JBrowseApp(
    assemblies=config["assemblies"],
    tracks=config.get("tracks", []),
    views=[
        {
            "type": "LinearGenomeView",
            "assembly": "hg38",
            "loc": "chr20:32,925,000-32,955,000",
            "tracks": [
                {
                    "trackId": "ncbi_refseq_hg38",
                    "type": "LinearBasicDisplay",
                    "height": 50,
                    "showOnlyGenes": True,
                    "geneGlyphMode": "longestCoding"
                },
                {
                    "trackId": "na12878_read_depth",
                    "type": "LinearMarkDisplay",
                    "height": 160
                }
            ]
        }
    ],
)
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/demos/read_marks/config.json \
  --assembly hg38 --loc chr20:32,925,000-32,955,000 \
  --track ncbi_refseq_hg38 display:LinearBasicDisplay '{"height":50,"showOnlyGenes":true,"geneGlyphMode":"longestCoding"}' \
  --track na12878_read_depth display:LinearMarkDisplay '{"height":160}' \
  --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 — swap an adapter uri for your own file and rerun.

cat > session.json <<'JSON'
{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr20:32,925,000-32,955,000",
      "tracks": [
        {
          "trackId": "ncbi_refseq_hg38",
          "type": "LinearBasicDisplay",
          "height": 50,
          "showOnlyGenes": true,
          "geneGlyphMode": "longestCoding"
        },
        {
          "trackId": "na12878_read_depth",
          "type": "LinearMarkDisplay",
          "height": 160
        }
      ]
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/demos/read_marks/config.json \
  --width 1500 --height 460 \
  --timeout 120000 \
  --spec session.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.

Plotting insert size as one point per read pair

A mark display draws one y axis. Depth runs in the tens and an insert size in the thousands, so the insert goes on a second track over the same file:

  • a point mark plots template_length per pair
  • a filter keeps the leftmost mate, where the template length is positive, so each pair counts once, and drops the few over 8 kb
  • the colour is mapping quality on a ramp pinned to 0 to 60

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

{
  "type": "AlignmentsTrack",
  "trackId": "na12878_read_marks",
  "name": "NA12878 insert size (1000 Genomes, 30x)",
  "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram",
  "assemblyNames": ["hg38"],
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "scales": { "y": { "title": "Insert size (bp)" } },
      "marks": [
        {
          "mark": "point",
          "transform": [
            {
              "type": "filter",
              "expr": "jexl:feature.template_length > 0 && feature.template_length < 8000"
            }
          ],
          "encoding": {
            "y": "template_length",
            "color": {
              "field": "score",
              "scale": "linear",
              "domainMin": 0,
              "domainMax": 60,
              "range": ["#bdbdbd", "#1f4e9a"],
              "title": "Mapping quality"
            }
          }
        }
      ]
    }
  ]
}
jbrowse add-track-json '{
  "type": "AlignmentsTrack",
  "trackId": "na12878_read_marks",
  "name": "NA12878 insert size (1000 Genomes, 30x)",
  "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram",
  "assemblyNames": ["hg38"],
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "scales": { "y": { "title": "Insert size (bp)" } },
      "marks": [
        {
          "mark": "point",
          "transform": [
            {
              "type": "filter",
              "expr": "jexl:feature.template_length > 0 && feature.template_length < 8000"
            }
          ],
          "encoding": {
            "y": "template_length",
            "color": {
              "field": "score",
              "scale": "linear",
              "domainMin": 0,
              "domainMax": 60,
              "range": ["#bdbdbd", "#1f4e9a"],
              "title": "Mapping quality"
            }
          }
        }
      ]
    }
  ]
}'

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": "AlignmentsTrack",
  "trackId": "na12878_read_marks",
  "name": "NA12878 insert size (1000 Genomes, 30x)",
  "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram",
  "assemblyNames": ["hg38"],
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "scales": { "y": { "title": "Insert size (bp)" } },
      "marks": [
        {
          "mark": "point",
          "transform": [
            {
              "type": "filter",
              "expr": "jexl:feature.template_length > 0 && feature.template_length < 8000"
            }
          ],
          "encoding": {
            "y": "template_length",
            "color": {
              "field": "score",
              "scale": "linear",
              "domainMin": 0,
              "domainMax": 60,
              "range": ["#bdbdbd", "#1f4e9a"],
              "title": "Mapping quality"
            }
          }
        }
      ]
    }
  ]
}

Open it under the depth track, on the same window.

The same window, the depth as bars above and each pair's insert size as a point below, each track with a separate y axis. The pairs sit in a low band, and over the left edge of the dip a second group appears well above it, in full blue.
Open in JBrowse Web ↗Open in JBrowse Desktop ↗

This session loads a large dataset and takes a while to open.

Or rebuild it from your own data:

Opening in Desktop needs JBrowse Desktop 5.0+, and saves the view as a reopenable session.

Nothing happens?

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%2Fdemos%2Fread_marks%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%22chr20%3A32%2C925%2C000-32%2C955%2C000%22%2C%22tracks%22%3A%5B%7B%22trackId%22%3A%22ncbi_refseq_hg38%22%2C%22type%22%3A%22LinearBasicDisplay%22%2C%22height%22%3A50%2C%22showOnlyGenes%22%3Atrue%2C%22geneGlyphMode%22%3A%22longestCoding%22%7D%2C%7B%22trackId%22%3A%22na12878_read_depth%22%2C%22type%22%3A%22LinearMarkDisplay%22%2C%22height%22%3A120%7D%2C%7B%22trackId%22%3A%22na12878_read_marks%22%2C%22type%22%3A%22LinearMarkDisplay%22%2C%22height%22%3A200%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+read_marks%2Finsert_size

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).This figure uses hg38, loaded from a bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  2. Add your own track: File → Open track..., then 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. Drag the bar at the bottom edge of the track to resize it (50px here).
    2. Track menu → Show only genes (checked)
    3. Track menu → Gene glyph → Representative transcript
  3. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a CRAM file (.cram + .crai).This figure uses “NA12878 depth (1000 Genomes, 30x)”.
    1. Drag the bar at the bottom edge of the track to resize it (120px here).
  4. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a CRAM file (.cram + .crai).This figure uses “NA12878 insert size (1000 Genomes, 30x)”.
    1. Drag the bar at the bottom edge of the track to resize it (200px here).
  5. Type your region of interest into the location box and press Enter.This figure is at chr20:32,925,000-32,955,000.

Its genome and tracks are hosted, so File → Open track... there adds your own files beside them.

With your own data, the steps behind the figure:

  1. Put your genome in your JBrowse Web's config.json: jbrowse add-assembly adds it.This figure uses hg38, loaded from a bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  2. Add your own track: File → Open track..., then 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. Drag the bar at the bottom edge of the track to resize it (50px here).
    2. Track menu → Show only genes (checked)
    3. Track menu → Gene glyph → Representative transcript
  3. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a CRAM file (.cram + .crai).This figure uses “NA12878 depth (1000 Genomes, 30x)”.
    1. Drag the bar at the bottom edge of the track to resize it (120px here).
  4. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a CRAM file (.cram + .crai).This figure uses “NA12878 insert size (1000 Genomes, 30x)”.
    1. Drag the bar at the bottom edge of the track to resize it (200px here).
  5. Type your region of interest into the location box and press Enter.This figure is at chr20:32,925,000-32,955,000.

This session spec draws the figure from https://jbrowse.org/demos/read_marks/config.json: paste it after &session=spec- on a JBrowse Web link that loads that config.

{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr20:32,925,000-32,955,000",
      "tracks": [
        {
          "trackId": "ncbi_refseq_hg38",
          "type": "LinearBasicDisplay",
          "height": 50,
          "showOnlyGenes": true,
          "geneGlyphMode": "longestCoding"
        },
        {
          "trackId": "na12878_read_depth",
          "type": "LinearMarkDisplay",
          "height": 120
        },
        {
          "trackId": "na12878_read_marks",
          "type": "LinearMarkDisplay",
          "height": 200
        }
      ]
    }
  ]
}

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/demos/read_marks/config.json")
app = JBrowseApp(
    assemblies=config["assemblies"],
    tracks=config.get("tracks", []),
    views=[
        {
            "type": "LinearGenomeView",
            "assembly": "hg38",
            "loc": "chr20:32,925,000-32,955,000",
            "tracks": [
                {
                    "trackId": "ncbi_refseq_hg38",
                    "type": "LinearBasicDisplay",
                    "height": 50,
                    "showOnlyGenes": True,
                    "geneGlyphMode": "longestCoding"
                },
                {
                    "trackId": "na12878_read_depth",
                    "type": "LinearMarkDisplay",
                    "height": 120
                },
                {
                    "trackId": "na12878_read_marks",
                    "type": "LinearMarkDisplay",
                    "height": 200
                }
            ]
        }
    ],
)
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/demos/read_marks/config.json \
  --assembly hg38 --loc chr20:32,925,000-32,955,000 \
  --track ncbi_refseq_hg38 display:LinearBasicDisplay '{"height":50,"showOnlyGenes":true,"geneGlyphMode":"longestCoding"}' \
  --track na12878_read_depth display:LinearMarkDisplay '{"height":120}' \
  --track na12878_read_marks display:LinearMarkDisplay '{"height":200}' \
  --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 — swap an adapter uri for your own file and rerun.

cat > session.json <<'JSON'
{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr20:32,925,000-32,955,000",
      "tracks": [
        {
          "trackId": "ncbi_refseq_hg38",
          "type": "LinearBasicDisplay",
          "height": 50,
          "showOnlyGenes": true,
          "geneGlyphMode": "longestCoding"
        },
        {
          "trackId": "na12878_read_depth",
          "type": "LinearMarkDisplay",
          "height": 120
        },
        {
          "trackId": "na12878_read_marks",
          "type": "LinearMarkDisplay",
          "height": 200
        }
      ]
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/demos/read_marks/config.json \
  --width 1500 --height 655 \
  --timeout 120000 \
  --spec session.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.

Each pair in the upper group straddles the missing 3.9 kb. On an alignments track score is the mapping quality; click a point to open its read.

Scanning chromosome 20 for clusters of long-insert pairs

Fetching every read of a chromosome overruns the track's size limit, so cut the long pairs out once, one row per pair, into a BED with a header naming its columns.

# one row per pair with an insert over 1 kb, from the leftmost mate to the
# end of the insert, with the mapping quality in the score column
# -q 20 drops reads the aligner could not place; -F 0x904 drops unmapped,
#   secondary and supplementary records
# REF_PATH lets htslib fetch each reference sequence the CRAM names by MD5
export REF_PATH='https://www.ebi.ac.uk/ena/cram/md5/%s'
samtools view -q 20 -F 0x904 --input-fmt-option required_fields=0x1DF NA12878.final.cram chr20 |
  awk 'BEGIN { OFS = "\t"; print "#chrom", "chromStart", "chromEnd", "name", "score", "strand", "tlen" }
    $7 == "=" && $9 > 1000 { print $3, $4 - 1, $4 - 1 + $9, $1, $5, "+", $9 }' |
  bgzip > NA12878.chr20.discordant_pairs.bed.gz
tabix -p bed NA12878.chr20.discordant_pairs.bed.gz

The BED holds few enough rows to fetch whole at any zoom. Two tracks read it:

  • a point per pair at the middle of its insert, tlen on y, coloured by score; a filter under 20 kb keeps the centromere's megabase inserts off the axis
  • a bar per bin counting pairs of 2 to 10 kb, on an axis pinned at 60 so the centromere saturates and a deletion's few dozen pairs stand up

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

{
  "type": "FeatureTrack",
  "trackId": "na12878_chr20_pairs",
  "name": "NA12878 chr20, pairs over 1 kb",
  "uri": "https://jbrowse.org/demos/read_marks/NA12878.chr20.discordant_pairs.bed.gz",
  "assemblyNames": ["hg38"],
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "scales": { "y": { "title": "Insert size (bp)" } },
      "marks": [
        {
          "mark": "point",
          "transform": [
            { "type": "filter", "expr": "jexl:feature.tlen < 20000" }
          ],
          "encoding": {
            "y": "tlen",
            "color": {
              "field": "score",
              "scale": "linear",
              "domainMin": 0,
              "domainMax": 60,
              "range": ["#bdbdbd", "#1f4e9a"],
              "title": "Mapping quality"
            }
          }
        }
      ]
    }
  ]
}
jbrowse add-track-json '{
  "type": "FeatureTrack",
  "trackId": "na12878_chr20_pairs",
  "name": "NA12878 chr20, pairs over 1 kb",
  "uri": "https://jbrowse.org/demos/read_marks/NA12878.chr20.discordant_pairs.bed.gz",
  "assemblyNames": ["hg38"],
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "scales": { "y": { "title": "Insert size (bp)" } },
      "marks": [
        {
          "mark": "point",
          "transform": [
            { "type": "filter", "expr": "jexl:feature.tlen < 20000" }
          ],
          "encoding": {
            "y": "tlen",
            "color": {
              "field": "score",
              "scale": "linear",
              "domainMin": 0,
              "domainMax": 60,
              "range": ["#bdbdbd", "#1f4e9a"],
              "title": "Mapping quality"
            }
          }
        }
      ]
    }
  ]
}'

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": "FeatureTrack",
  "trackId": "na12878_chr20_pairs",
  "name": "NA12878 chr20, pairs over 1 kb",
  "uri": "https://jbrowse.org/demos/read_marks/NA12878.chr20.discordant_pairs.bed.gz",
  "assemblyNames": ["hg38"],
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "scales": { "y": { "title": "Insert size (bp)" } },
      "marks": [
        {
          "mark": "point",
          "transform": [
            { "type": "filter", "expr": "jexl:feature.tlen < 20000" }
          ],
          "encoding": {
            "y": "tlen",
            "color": {
              "field": "score",
              "scale": "linear",
              "domainMin": 0,
              "domainMax": 60,
              "range": ["#bdbdbd", "#1f4e9a"],
              "title": "Mapping quality"
            }
          }
        }
      ]
    }
  ]
}

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

{
  "type": "FeatureTrack",
  "trackId": "na12878_chr20_pair_counts",
  "name": "NA12878 chr20, pairs of 2 to 10 kb per bin",
  "uri": "https://jbrowse.org/demos/read_marks/NA12878.chr20.discordant_pairs.bed.gz",
  "assemblyNames": ["hg38"],
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "scales": {
        "y": { "domainMin": 0, "domainMax": 60, "title": "Pairs per bin" }
      },
      "marks": [
        {
          "mark": "bar",
          "transform": [
            {
              "type": "filter",
              "expr": "jexl:feature.tlen > 2000 && feature.tlen < 10000"
            },
            { "type": "bin", "step": "auto" },
            { "type": "aggregate", "ops": [{ "op": "count" }] }
          ],
          "encoding": { "color": { "value": "#d62728" } }
        }
      ]
    }
  ]
}
jbrowse add-track-json '{
  "type": "FeatureTrack",
  "trackId": "na12878_chr20_pair_counts",
  "name": "NA12878 chr20, pairs of 2 to 10 kb per bin",
  "uri": "https://jbrowse.org/demos/read_marks/NA12878.chr20.discordant_pairs.bed.gz",
  "assemblyNames": ["hg38"],
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "scales": {
        "y": { "domainMin": 0, "domainMax": 60, "title": "Pairs per bin" }
      },
      "marks": [
        {
          "mark": "bar",
          "transform": [
            {
              "type": "filter",
              "expr": "jexl:feature.tlen > 2000 && feature.tlen < 10000"
            },
            { "type": "bin", "step": "auto" },
            { "type": "aggregate", "ops": [{ "op": "count" }] }
          ],
          "encoding": { "color": { "value": "#d62728" } }
        }
      ]
    }
  ]
}'

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": "FeatureTrack",
  "trackId": "na12878_chr20_pair_counts",
  "name": "NA12878 chr20, pairs of 2 to 10 kb per bin",
  "uri": "https://jbrowse.org/demos/read_marks/NA12878.chr20.discordant_pairs.bed.gz",
  "assemblyNames": ["hg38"],
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "scales": {
        "y": { "domainMin": 0, "domainMax": 60, "title": "Pairs per bin" }
      },
      "marks": [
        {
          "mark": "bar",
          "transform": [
            {
              "type": "filter",
              "expr": "jexl:feature.tlen > 2000 && feature.tlen < 10000"
            },
            { "type": "bin", "step": "auto" },
            { "type": "aggregate", "ops": [{ "op": "count" }] }
          ],
          "encoding": { "color": { "value": "#d62728" } }
        }
      ]
    }
  ]
}

In the counts track the aggregate step names no groupby, so it groups by the edges the bin step above it wrote.

Open both tracks on the whole of chr20.

Chromosome 20 end to end. Every pair with an insert under 20 kb is a point at its insert size, and the red bars on the track under it count the pairs between 2 and 10 kb per bin. The centromere, pinched in the banding above the ruler, and the repeats flanking it saturate both; outside them the bars rise in a handful of places, each under a short stack of dark points.
Open in JBrowse Web ↗Open in JBrowse Desktop ↗

This session loads a large dataset and takes a while to open.

Or rebuild it from your own data:

Opening in Desktop needs JBrowse Desktop 5.0+, and saves the view as a reopenable session.

Nothing happens?

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%2Fdemos%2Fread_marks%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%22chr20%22%2C%22tracks%22%3A%5B%7B%22trackId%22%3A%22na12878_chr20_pairs%22%2C%22type%22%3A%22LinearMarkDisplay%22%2C%22height%22%3A160%7D%2C%7B%22trackId%22%3A%22na12878_chr20_pair_counts%22%2C%22type%22%3A%22LinearMarkDisplay%22%2C%22height%22%3A100%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+read_marks%2Fchromosome

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).This figure uses hg38, loaded from a bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  2. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a bgzip-compressed, tabix-indexed BED (.bed.gz + .tbi).This figure uses “NA12878 chr20, pairs over 1 kb”.
    1. Drag the bar at the bottom edge of the track to resize it (160px here).
  3. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a bgzip-compressed, tabix-indexed BED (.bed.gz + .tbi).This figure uses “NA12878 chr20, pairs of 2 to 10 kb per bin”.
    1. Drag the bar at the bottom edge of the track to resize it (100px here).
  4. Type your region of interest into the location box and press Enter.This figure is at chr20.

Its genome and tracks are hosted, so File → Open track... there adds your own files beside them.

With your own data, the steps behind the figure:

  1. Put your genome in your JBrowse Web's config.json: jbrowse add-assembly adds it.This figure uses hg38, loaded from a bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  2. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a bgzip-compressed, tabix-indexed BED (.bed.gz + .tbi).This figure uses “NA12878 chr20, pairs over 1 kb”.
    1. Drag the bar at the bottom edge of the track to resize it (160px here).
  3. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a bgzip-compressed, tabix-indexed BED (.bed.gz + .tbi).This figure uses “NA12878 chr20, pairs of 2 to 10 kb per bin”.
    1. Drag the bar at the bottom edge of the track to resize it (100px here).
  4. Type your region of interest into the location box and press Enter.This figure is at chr20.

This session spec draws the figure from https://jbrowse.org/demos/read_marks/config.json: paste it after &session=spec- on a JBrowse Web link that loads that config.

{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr20",
      "tracks": [
        {
          "trackId": "na12878_chr20_pairs",
          "type": "LinearMarkDisplay",
          "height": 160
        },
        {
          "trackId": "na12878_chr20_pair_counts",
          "type": "LinearMarkDisplay",
          "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/demos/read_marks/config.json")
app = JBrowseApp(
    assemblies=config["assemblies"],
    tracks=config.get("tracks", []),
    views=[
        {
            "type": "LinearGenomeView",
            "assembly": "hg38",
            "loc": "chr20",
            "tracks": [
                {
                    "trackId": "na12878_chr20_pairs",
                    "type": "LinearMarkDisplay",
                    "height": 160
                },
                {
                    "trackId": "na12878_chr20_pair_counts",
                    "type": "LinearMarkDisplay",
                    "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/demos/read_marks/config.json \
  --assembly hg38 --loc chr20 \
  --track na12878_chr20_pairs display:LinearMarkDisplay '{"height":160}' \
  --track na12878_chr20_pair_counts display:LinearMarkDisplay '{"height":100}' \
  --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 — swap an adapter uri for your own file and rerun.

cat > session.json <<'JSON'
{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr20",
      "tracks": [
        {
          "trackId": "na12878_chr20_pairs",
          "type": "LinearMarkDisplay",
          "height": 160
        },
        {
          "trackId": "na12878_chr20_pair_counts",
          "type": "LinearMarkDisplay",
          "height": 100
        }
      ]
    }
  ]
}
JSON

cat > callouts.json <<'JSON'
[
  {
    "type": "text",
    "text": "a tall red bar: a deletion on both copies",
    "maxWidth": 420,
    "fontSize": 18,
    "leader": true,
    "anchor": {
      "track": "na12878_chr20_pair_counts",
      "locus": "chr20:34,250,000",
      "fracY": 0.3
    },
    "dx": 80,
    "dy": -8
  },
  {
    "type": "text",
    "text": "centromere: pairs mis-mapped in repeats",
    "fontSize": 18,
    "leader": true,
    "anchor": {
      "track": "na12878_chr20_pairs",
      "locus": "chr20:26,300,000",
      "fracY": 0.5
    },
    "dx": -80
  }
]
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/demos/read_marks/config.json \
  --width 1500 --height 510 \
  --timeout 120000 \
  --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 bar at 34.2 Mb is a homozygous deletion, and the one at 32.9 Mb is the EFCAB8 deletion from the sections above.

Checking the chr20 bars against the 1000 Genomes SV callset

bcftools lists every deletion over 2 kb that the callset gives NA12878 on the chromosome:

# -s keeps one sample's genotypes; -i then keeps the rows where that sample
#   carries the allele
bcftools view -s NA12878 1KGP_3202.gatksv_svtools_novelins.freeze_V3.wAF.vcf.gz chr20 |
  bcftools query -i 'GT="alt" && INFO/SVTYPE="DEL" && INFO/SVLEN<-2000' \
    -f '%CHROM\t%POS\t%END\t%INFO/SVLEN\t[%GT]\t%INFO/AF\t%INFO/EVIDENCE\n'

Every deletion in that listing between 2 and 10 kb is a bar in the pair-count track, and the two homozygous ones are tallest. Other windows hold ten or more such pairs with no call: the chromosome start and 1.4, 2.8, 32.7 and 48.5 Mb.

To check the EFCAB8 deletion against the reads, compare the depth inside the call with the depth beside it, and count the long pairs around it:

samtools coverage -r chr20:32937680-32941583 NA12878.final.cram | cut -f 1-3,7
samtools coverage -r chr20:32930000-32937000 NA12878.final.cram | cut -f 1-3,7
samtools view -q 20 NA12878.final.cram chr20:32935000-32944000 |
  awk '{ t = $9 < 0 ? -$9 : $9; if (t > 2000) big++; else if (t > 0) norm++ }
    END { print norm " pairs at the library insert, " big " over 2 kb" }'

Depth inside the call comes out about half the depth beside it, as a heterozygous deletion predicts. Around the call, almost every pair sits at the library insert, and a few dozen exceed 2 kb.

Reproduce it end to end

build_read_marks.sh runs the commands above:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_read_marks.sh
bash build_read_marks.sh                     # builds ./read_marks_build/jbrowse2
npx --yes serve read_marks_build/jbrowse2    # then open the printed URL

With no arguments the script builds the depth, insert-size and chromosome-scan tracks over NA12878. Given your own reads, bash build_read_marks.sh reads.cram genome.fa builds them over your file, and CHROM picks the chromosome to scan.

See also

Citations