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.13. 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 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. The plots use JBrowse's mark display, which maps feature fields to plot channels. In grammar-of-graphics terms, a mark is the shape drawn and scales.y the y scale. 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 and 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.

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

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.

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"}'

hg38 is one of the genomes JBrowse Desktop hosts, with gene tracks already set up: on the start screen click Show all available genomes and pick it. To load the files of this config instead:

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

https://jbrowse.org/genomes/GRCh38/fasta/hg38.prefix.fa.gz
https://jbrowse.org/genomes/GRCh38/fasta/hg38.prefix.fa.gz.fai
https://jbrowse.org/genomes/GRCh38/fasta/hg38.prefix.fa.gz.gzi

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

  • Genome name: hg38
  • refName aliases (under More options): https://s3.amazonaws.com/jbrowse.org/genomes/GRCh38/hg38_aliases.txt
  • cytobands (under More options): 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.

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 this track in JBrowse Web ↗

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

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 ↗
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%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

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 bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  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 genes”.
    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 CRAM file (.cram + .crai).This figure uses “NA12878 depth (1000 Genomes, 30x)”.
  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. It goes after &session=spec- on a JBrowse Web link that loads https://jbrowse.org/demos/read_marks/config.json, which is what the Web button above opens.

{
  "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: 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": "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. Its filter keeps the leftmost mate, where the template length is positive, so each pair counts once, and drops the few over 8 kb. On an alignments track score is the mapping quality, which colors the points.

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 in 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 ↗
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%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

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 bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  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 genes”.
    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 CRAM file (.cram + .crai).This figure uses “NA12878 depth (1000 Genomes, 30x)”.
  4. Add the next track the same way. This one needs a CRAM file (.cram + .crai).This figure uses “NA12878 insert size (1000 Genomes, 30x)”.
  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. It goes after &session=spec- on a JBrowse Web link that loads https://jbrowse.org/demos/read_marks/config.json, which is what the Web button above opens.

{
  "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: 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": "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. 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 we 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, and two tracks read it:

  • the point track's filter under 20 kb keeps the centromere's megabase inserts off the axis
  • the count track's axis, pinned at 60, lets the centromere saturate so 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"
            }
          }
        }
      ]
    }
  ]
}

Open this track in JBrowse Web ↗

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

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" } }
        }
      ]
    }
  ]
}

Open this track in JBrowse Web ↗

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

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 ↗
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%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

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 bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  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 BED (.bed.gz + .tbi).This figure uses “NA12878 chr20, pairs over 1 kb”.
  3. Add the next track the same way. 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”.
  4. Type your region of interest into the location box and press Enter.This figure is at chr20.

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

{
  "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: 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": "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": {
      "trackId": "na12878_chr20_pair_counts",
      "loc": "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": {
      "trackId": "na12878_chr20_pairs",
      "loc": "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.

Checking the chr20 bars against the 1000 Genomes SV callset

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

# -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