Skip to main content

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

LinearMarkDisplay is a grammar of graphics over a track: each entry in marks names a shape, a transform list and an encoding from feature fields to channels, the way a BED column feeds a plot in the Alu tutorial. Here the file is a BAM, and the fields are the ones the aligner wrote. A pair straddling a deletion maps with a long insert, and a heterozygous deletion halves the depth, so plotting those two fields finds the deletion without a caller. The mark display is experimental, and its config shape may change.

Prerequisites

  • a JBrowse to open the figures' sessions in (Web or Desktop); every file here is a URL, so nothing needs hosting to read along
  • 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:

The window

30 kb of an EFCAB8 intron on chromosome 20, where the callset says NA12878 carries one copy of a 3.9 kb deletion.

Depth as a coverage step

A bar mark over a coverage transform, which replaces the reads with runs of constant depth.

{
  "type": "AlignmentsTrack",
  "trackId": "na12878_read_depth",
  "name": "NA12878 reads (1000 Genomes, 30x)",
  "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"
    }
  },
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "displayId": "na12878_read_depth-LinearMarkDisplay",
      "marks": [
        {
          "shape": "bar",
          "transform": [{ "type": "coverage" }],
          "encoding": { "y": "coverage", "color": "#c8d8ee" }
        }
      ]
    }
  ]
}
jbrowse add-track-json '{
  "type": "AlignmentsTrack",
  "trackId": "na12878_read_depth",
  "name": "NA12878 reads (1000 Genomes, 30x)",
  "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"
    }
  },
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "displayId": "na12878_read_depth-LinearMarkDisplay",
      "marks": [
        {
          "shape": "bar",
          "transform": [{ "type": "coverage" }],
          "encoding": { "y": "coverage", "color": "#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 reads (1000 Genomes, 30x)",
  "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"
    }
  },
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "displayId": "na12878_read_depth-LinearMarkDisplay",
      "marks": [
        {
          "shape": "bar",
          "transform": [{ "type": "coverage" }],
          "encoding": { "y": "coverage", "color": "#c8d8ee" }
        }
      ]
    }
  ]
}

Open it at chr20:32,925,000-32,955,000.

Thirty kilobases of an EFCAB8 intron in NA12878, the read depth as bars. Between 32,937,500 and 32,941,500 the depth runs at about half of what it is on either side. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

The steps behind the figure above, run against your own data.

  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.
  2. Add your own track: File → Open track..., then paste a URL or choose a local file.This figure uses the “ncbi_refseq_hg38” track.
  3. Drag the bar at the bottom edge of the track to resize it.“ncbi_refseq_hg38”.
  4. Track menu → Show only genes (checked)“ncbi_refseq_hg38”.
  5. Track menu → Gene glyph → Representative transcript“ncbi_refseq_hg38”.
  6. 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 reads (1000 Genomes, 30x)”.
  7. Drag the bar at the bottom edge of the track to resize it.“NA12878 reads (1000 Genomes, 30x)”.
  8. Type your region of interest into the location box and press Enter.This figure is at chr20:32,925,000-32,955,000.

Open this view in JBrowse Desktop ↗

Opens JBrowse Desktop (JBrowse Desktop 5.0+) at this view and saves it as a reopenable session — swap in your own files afterwards.

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%22sessionTracks%22%3A%5B%7B%22type%22%3A%22AlignmentsTrack%22%2C%22trackId%22%3A%22na12878_read_depth%22%2C%22name%22%3A%22NA12878+reads+%281000+Genomes%2C+30x%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%2C%22displays%22%3A%5B%7B%22type%22%3A%22LinearMarkDisplay%22%2C%22displayId%22%3A%22na12878_read_depth-LinearMarkDisplay%22%2C%22marks%22%3A%5B%7B%22shape%22%3A%22bar%22%2C%22transform%22%3A%5B%7B%22type%22%3A%22coverage%22%7D%5D%2C%22encoding%22%3A%7B%22y%22%3A%22coverage%22%2C%22color%22%3A%22%23c8d8ee%22%7D%7D%5D%7D%5D%7D%5D%2C%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

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_read_depth",
  "name": "NA12878 reads (1000 Genomes, 30x)",
  "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"
    }
  },
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "displayId": "na12878_read_depth-LinearMarkDisplay",
      "marks": [
        {
          "shape": "bar",
          "transform": [
            {
              "type": "coverage"
            }
          ],
          "encoding": {
            "y": "coverage",
            "color": "#c8d8ee"
          }
        }
      ]
    }
  ]
}'

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 Session spec tab, or jbrowse set-default-session.

This session spec pastes after &session=spec- on any JBrowse Web instance.

{
  "sessionTracks": [
    {
      "type": "AlignmentsTrack",
      "trackId": "na12878_read_depth",
      "name": "NA12878 reads (1000 Genomes, 30x)",
      "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"
        }
      },
      "displays": [
        {
          "type": "LinearMarkDisplay",
          "displayId": "na12878_read_depth-LinearMarkDisplay",
          "marks": [
            {
              "shape": "bar",
              "transform": [
                {
                  "type": "coverage"
                }
              ],
              "encoding": {
                "y": "coverage",
                "color": "#c8d8ee"
              }
            }
          ]
        }
      ]
    }
  ],
  "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
        }
      ]
    }
  ]
}

Loaded against config: https://jbrowse.org/demos/read_marks/config.json

The same view with jbrowse-anywidget. Point the adapter at your own file.

from jbrowse_anywidget import LinearGenomeView, fetch_hub

view = LinearGenomeView(
    assembly=fetch_hub("hg38"),
    location="chr20:32,925,000-32,955,000",
)

view.add_track({
    "type": "FeatureTrack",
    "trackId": "my_track_1",
    "name": "My track 1",
    "assemblyNames": ["hg38"],
    "adapter": {"type": "BamAdapter", "uri": "https://your-server/your-file"},
})
view.add_track({
    "type": "AlignmentsTrack",
    "trackId": "my_track_2",
    "name": "My track 2",
    "assemblyNames": ["hg38"],
    "adapter": {"type": "CramAdapter", "uri": "https://your-server/your-file"},
})

view  # display the widget

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'
{
  "sessionTracks": [
    {
      "type": "AlignmentsTrack",
      "trackId": "na12878_read_depth",
      "name": "NA12878 reads (1000 Genomes, 30x)",
      "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"
        }
      },
      "displays": [
        {
          "type": "LinearMarkDisplay",
          "displayId": "na12878_read_depth-LinearMarkDisplay",
          "marks": [
            {
              "shape": "bar",
              "transform": [
                {
                  "type": "coverage"
                }
              ],
              "encoding": {
                "y": "coverage",
                "color": "#c8d8ee"
              }
            }
          ]
        }
      ]
    }
  ],
  "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 \
  --session 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.

The CRAM decodes against the assembly the track is added to.

Insert size as a point per pair

A point per pair with template_length on y, coloured by mapping quality on a ramp pinned to 0 to 60. A filter keeps the leftmost mate, where the template length is positive, so each pair counts once. Depth moves to its own axis with "resolve": "independent".

"marks": [
  {
    "shape": "bar",
    "transform": [{ "type": "coverage" }],
    "encoding": {
      "y": { "field": "coverage", "resolve": "independent" },
      "color": "#c8d8ee"
    }
  },
  {
    "shape": "point",
    "transform": [
      {
        "type": "filter",
        "expr": "jexl:feature.template_length > 0 && feature.template_length < 8000"
      }
    ],
    "encoding": {
      "y": "template_length",
      "color": {
        "field": "score",
        "scale": "linear",
        "domain": [0, 60],
        "ramp": ["#bdbdbd", "#1f4e9a"]
      }
    }
  }
]
The same window with each pair's insert size as a point on the left axis and the depth on its own axis on the right. The pairs sit in a band under 1,000 bases, and over the left edge of the dip a second group appears between 4,300 and 4,700 bases, in the full blue of a mapping quality of 60. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

The steps behind the figure above, run against your own data.

  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.
  2. Add your own track: File → Open track..., then paste a URL or choose a local file.This figure uses the “ncbi_refseq_hg38” track.
  3. Drag the bar at the bottom edge of the track to resize it.“ncbi_refseq_hg38”.
  4. Track menu → Show only genes (checked)“ncbi_refseq_hg38”.
  5. Track menu → Gene glyph → Representative transcript“ncbi_refseq_hg38”.
  6. Add your own track: File → Open track..., then paste a URL or choose a local file.This figure uses the “na12878_read_marks” track.
  7. Drag the bar at the bottom edge of the track to resize it.“na12878_read_marks”.
  8. Type your region of interest into the location box and press Enter.This figure is at chr20:32,925,000-32,955,000.

Open this view in JBrowse Desktop ↗

Opens JBrowse Desktop (JBrowse Desktop 5.0+) at this view and saves it as a reopenable session — swap in your own files afterwards.

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%22sessionTracks%22%3A%5B%5D%2C%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_marks%22%2C%22type%22%3A%22LinearMarkDisplay%22%2C%22height%22%3A260%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+read_marks%2Finsert_size

This session spec pastes after &session=spec- on any JBrowse Web instance.

{
  "sessionTracks": [],
  "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_marks",
          "type": "LinearMarkDisplay",
          "height": 260
        }
      ]
    }
  ]
}

Loaded against config: https://jbrowse.org/demos/read_marks/config.json

The same view with jbrowse-anywidget. Point the adapter at your own file.

from jbrowse_anywidget import LinearGenomeView, fetch_hub

view = LinearGenomeView(
    assembly=fetch_hub("hg38"),
    location="chr20:32,925,000-32,955,000",
)

view.add_track({
    "type": "FeatureTrack",
    "trackId": "my_track_1",
    "name": "My track 1",
    "assemblyNames": ["hg38"],
    "adapter": {"type": "BamAdapter", "uri": "https://your-server/your-file"},
})
view.add_track({
    "type": "FeatureTrack",
    "trackId": "my_track_2",
    "name": "My track 2",
    "assemblyNames": ["hg38"],
    "adapter": {"type": "BamAdapter", "uri": "https://your-server/your-file"},
})

view  # display the widget

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'
{
  "sessionTracks": [],
  "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_marks",
          "type": "LinearMarkDisplay",
          "height": 260
        }
      ]
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/demos/read_marks/config.json \
  --session 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. score is the mapping quality on every track type. Hover a point for its values; click it to open the read.

Which reads carry the long inserts

A span over a stack transform is a pileup. A formula step writes the unsigned insert so both mates share a colour, and a ramp pinned at 5 kb paints a spanning pair red.

"marks": [
  {
    "shape": "span",
    "transform": [
      { "type": "formula", "expr": "jexl:abs(feature.template_length)", "as": "insert" },
      { "type": "stack" }
    ],
    "encoding": {
      "row": "row",
      "color": {
        "field": "insert",
        "scale": "linear",
        "domain": [0, 5000],
        "ramp": ["#c8d8ee", "#d62728"]
      }
    }
  }
]

Zoom to the left edge of the dip, chr20:32,936,200-32,939,200.

The left breakpoint at 3 kb, the reads stacked and coloured by their pair's insert. The red reads end together at 32,937,680, where their mates lie 4 kb to the right; the pale reads run across it, and thin out on the far side. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

The steps behind the figure above, run against your own data.

  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.
  2. Add your own track: File → Open track..., then paste a URL or choose a local file.This figure uses the “ncbi_refseq_hg38” track.
  3. Drag the bar at the bottom edge of the track to resize it.“ncbi_refseq_hg38”.
  4. Track menu → Show only genes (checked)“ncbi_refseq_hg38”.
  5. Track menu → Gene glyph → Representative transcript“ncbi_refseq_hg38”.
  6. 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 reads (1000 Genomes, 30x)”.
  7. Drag the bar at the bottom edge of the track to resize it.“NA12878 reads (1000 Genomes, 30x)”.
  8. Type your region of interest into the location box and press Enter.This figure is at chr20:32,936,200-32,939,200.

Open this view in JBrowse Desktop ↗

Opens JBrowse Desktop (JBrowse Desktop 5.0+) at this view and saves it as a reopenable session — swap in your own files afterwards.

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%22sessionTracks%22%3A%5B%7B%22type%22%3A%22AlignmentsTrack%22%2C%22trackId%22%3A%22na12878_read_pileup%22%2C%22name%22%3A%22NA12878+reads+%281000+Genomes%2C+30x%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%2C%22displays%22%3A%5B%7B%22type%22%3A%22LinearMarkDisplay%22%2C%22displayId%22%3A%22na12878_read_pileup-LinearMarkDisplay%22%2C%22marks%22%3A%5B%7B%22shape%22%3A%22span%22%2C%22transform%22%3A%5B%7B%22type%22%3A%22formula%22%2C%22expr%22%3A%22jexl%3Aabs%28feature.template_length%29%22%2C%22as%22%3A%22insert%22%7D%2C%7B%22type%22%3A%22stack%22%7D%5D%2C%22encoding%22%3A%7B%22row%22%3A%22row%22%2C%22color%22%3A%7B%22field%22%3A%22insert%22%2C%22scale%22%3A%22linear%22%2C%22domain%22%3A%5B0%2C5000%5D%2C%22ramp%22%3A%5B%22%23c8d8ee%22%2C%22%23d62728%22%5D%7D%7D%7D%5D%7D%5D%7D%5D%2C%22views%22%3A%5B%7B%22type%22%3A%22LinearGenomeView%22%2C%22assembly%22%3A%22hg38%22%2C%22loc%22%3A%22chr20%3A32%2C936%2C200-32%2C939%2C200%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_pileup%22%2C%22type%22%3A%22LinearMarkDisplay%22%2C%22height%22%3A300%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+read_marks%2Fpileup

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_read_pileup",
  "name": "NA12878 reads (1000 Genomes, 30x)",
  "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"
    }
  },
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "displayId": "na12878_read_pileup-LinearMarkDisplay",
      "marks": [
        {
          "shape": "span",
          "transform": [
            {
              "type": "formula",
              "expr": "jexl:abs(feature.template_length)",
              "as": "insert"
            },
            {
              "type": "stack"
            }
          ],
          "encoding": {
            "row": "row",
            "color": {
              "field": "insert",
              "scale": "linear",
              "domain": [
                0,
                5000
              ],
              "ramp": [
                "#c8d8ee",
                "#d62728"
              ]
            }
          }
        }
      ]
    }
  ]
}'

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 Session spec tab, or jbrowse set-default-session.

This session spec pastes after &session=spec- on any JBrowse Web instance.

{
  "sessionTracks": [
    {
      "type": "AlignmentsTrack",
      "trackId": "na12878_read_pileup",
      "name": "NA12878 reads (1000 Genomes, 30x)",
      "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"
        }
      },
      "displays": [
        {
          "type": "LinearMarkDisplay",
          "displayId": "na12878_read_pileup-LinearMarkDisplay",
          "marks": [
            {
              "shape": "span",
              "transform": [
                {
                  "type": "formula",
                  "expr": "jexl:abs(feature.template_length)",
                  "as": "insert"
                },
                {
                  "type": "stack"
                }
              ],
              "encoding": {
                "row": "row",
                "color": {
                  "field": "insert",
                  "scale": "linear",
                  "domain": [
                    0,
                    5000
                  ],
                  "ramp": [
                    "#c8d8ee",
                    "#d62728"
                  ]
                }
              }
            }
          ]
        }
      ]
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr20:32,936,200-32,939,200",
      "tracks": [
        {
          "trackId": "ncbi_refseq_hg38",
          "type": "LinearBasicDisplay",
          "height": 50,
          "showOnlyGenes": true,
          "geneGlyphMode": "longestCoding"
        },
        {
          "trackId": "na12878_read_pileup",
          "type": "LinearMarkDisplay",
          "height": 300
        }
      ]
    }
  ]
}

Loaded against config: https://jbrowse.org/demos/read_marks/config.json

The same view with jbrowse-anywidget. Point the adapter at your own file.

from jbrowse_anywidget import LinearGenomeView, fetch_hub

view = LinearGenomeView(
    assembly=fetch_hub("hg38"),
    location="chr20:32,936,200-32,939,200",
)

view.add_track({
    "type": "FeatureTrack",
    "trackId": "my_track_1",
    "name": "My track 1",
    "assemblyNames": ["hg38"],
    "adapter": {"type": "BamAdapter", "uri": "https://your-server/your-file"},
})
view.add_track({
    "type": "AlignmentsTrack",
    "trackId": "my_track_2",
    "name": "My track 2",
    "assemblyNames": ["hg38"],
    "adapter": {"type": "CramAdapter", "uri": "https://your-server/your-file"},
})

view  # display the widget

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'
{
  "sessionTracks": [
    {
      "type": "AlignmentsTrack",
      "trackId": "na12878_read_pileup",
      "name": "NA12878 reads (1000 Genomes, 30x)",
      "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"
        }
      },
      "displays": [
        {
          "type": "LinearMarkDisplay",
          "displayId": "na12878_read_pileup-LinearMarkDisplay",
          "marks": [
            {
              "shape": "span",
              "transform": [
                {
                  "type": "formula",
                  "expr": "jexl:abs(feature.template_length)",
                  "as": "insert"
                },
                {
                  "type": "stack"
                }
              ],
              "encoding": {
                "row": "row",
                "color": {
                  "field": "insert",
                  "scale": "linear",
                  "domain": [
                    0,
                    5000
                  ],
                  "ramp": [
                    "#c8d8ee",
                    "#d62728"
                  ]
                }
              }
            }
          ]
        }
      ]
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr20:32,936,200-32,939,200",
      "tracks": [
        {
          "trackId": "ncbi_refseq_hg38",
          "type": "LinearBasicDisplay",
          "height": 50,
          "showOnlyGenes": true,
          "geneGlyphMode": "longestCoding"
        },
        {
          "trackId": "na12878_read_pileup",
          "type": "LinearMarkDisplay",
          "height": 300
        }
      ]
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/demos/read_marks/config.json \
  --session 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.

Pin the domain: an unpinned ramp spans the values on screen, so a window with no spanning pair would paint its longest ordinary insert red.

Scanning the chromosome for the same signature

Fetching every read of a chromosome overruns the byte budget, 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

11,327 rows, small enough to fetch whole at any zoom. A FeatureTrack with two marks:

  • a point per pair, tlen on y, coloured by score. x2: start draws the pair at its leftmost read; 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 a right-hand axis pinned at 60 so the centromere saturates and a deletion's ten to fifty stand up.
{
  "type": "FeatureTrack",
  "trackId": "na12878_chr20_pairs",
  "name": "NA12878 chr20, pairs over 1 kb",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "BedTabixAdapter",
    "bedGzLocation": {
      "uri": "https://jbrowse.org/demos/read_marks/NA12878.chr20.discordant_pairs.bed.gz"
    },
    "index": {
      "location": {
        "uri": "https://jbrowse.org/demos/read_marks/NA12878.chr20.discordant_pairs.bed.gz.tbi"
      }
    }
  },
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "displayId": "na12878_chr20_pairs-LinearMarkDisplay",
      "marks": [
        {
          "shape": "point",
          "transform": [
            { "type": "filter", "expr": "jexl:feature.tlen < 20000" }
          ],
          "encoding": {
            "x2": "start",
            "y": "tlen",
            "color": {
              "field": "score",
              "scale": "linear",
              "domain": [0, 60],
              "ramp": ["#bdbdbd", "#1f4e9a"]
            }
          }
        },
        {
          "shape": "bar",
          "transform": [
            {
              "type": "filter",
              "expr": "jexl:feature.tlen > 2000 && feature.tlen < 10000"
            },
            { "type": "bin", "step": "auto" },
            {
              "type": "aggregate",
              "groupby": ["start", "end"],
              "ops": [{ "op": "count" }]
            }
          ],
          "encoding": {
            "y": {
              "field": "count",
              "resolve": "independent",
              "domain": [0, 60]
            },
            "color": "#d62728"
          }
        }
      ]
    }
  ]
}
jbrowse add-track-json '{
  "type": "FeatureTrack",
  "trackId": "na12878_chr20_pairs",
  "name": "NA12878 chr20, pairs over 1 kb",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "BedTabixAdapter",
    "bedGzLocation": {
      "uri": "https://jbrowse.org/demos/read_marks/NA12878.chr20.discordant_pairs.bed.gz"
    },
    "index": {
      "location": {
        "uri": "https://jbrowse.org/demos/read_marks/NA12878.chr20.discordant_pairs.bed.gz.tbi"
      }
    }
  },
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "displayId": "na12878_chr20_pairs-LinearMarkDisplay",
      "marks": [
        {
          "shape": "point",
          "transform": [
            { "type": "filter", "expr": "jexl:feature.tlen < 20000" }
          ],
          "encoding": {
            "x2": "start",
            "y": "tlen",
            "color": {
              "field": "score",
              "scale": "linear",
              "domain": [0, 60],
              "ramp": ["#bdbdbd", "#1f4e9a"]
            }
          }
        },
        {
          "shape": "bar",
          "transform": [
            {
              "type": "filter",
              "expr": "jexl:feature.tlen > 2000 && feature.tlen < 10000"
            },
            { "type": "bin", "step": "auto" },
            {
              "type": "aggregate",
              "groupby": ["start", "end"],
              "ops": [{ "op": "count" }]
            }
          ],
          "encoding": {
            "y": {
              "field": "count",
              "resolve": "independent",
              "domain": [0, 60]
            },
            "color": "#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_pairs",
  "name": "NA12878 chr20, pairs over 1 kb",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "BedTabixAdapter",
    "bedGzLocation": {
      "uri": "https://jbrowse.org/demos/read_marks/NA12878.chr20.discordant_pairs.bed.gz"
    },
    "index": {
      "location": {
        "uri": "https://jbrowse.org/demos/read_marks/NA12878.chr20.discordant_pairs.bed.gz.tbi"
      }
    }
  },
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "displayId": "na12878_chr20_pairs-LinearMarkDisplay",
      "marks": [
        {
          "shape": "point",
          "transform": [
            { "type": "filter", "expr": "jexl:feature.tlen < 20000" }
          ],
          "encoding": {
            "x2": "start",
            "y": "tlen",
            "color": {
              "field": "score",
              "scale": "linear",
              "domain": [0, 60],
              "ramp": ["#bdbdbd", "#1f4e9a"]
            }
          }
        },
        {
          "shape": "bar",
          "transform": [
            {
              "type": "filter",
              "expr": "jexl:feature.tlen > 2000 && feature.tlen < 10000"
            },
            { "type": "bin", "step": "auto" },
            {
              "type": "aggregate",
              "groupby": ["start", "end"],
              "ops": [{ "op": "count" }]
            }
          ],
          "encoding": {
            "y": {
              "field": "count",
              "resolve": "independent",
              "domain": [0, 60]
            },
            "color": "#d62728"
          }
        }
      ]
    }
  ]
}
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 right axis count the pairs between 2 and 10 kb per bin. The centromere, from 26 to 32 Mb, saturates both; outside it the bars stand up in a handful of places, each under a short stack of dark points. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

The steps behind the figure above, run against your own data.

  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.
  2. Add your own track: File → Open track..., then paste a URL or choose a local file.This figure uses the “na12878_chr20_pairs” track.
  3. Drag the bar at the bottom edge of the track to resize it.
  4. Type your region of interest into the location box and press Enter.This figure is at chr20.

Open this view in JBrowse Desktop ↗

Opens JBrowse Desktop (JBrowse Desktop 5.0+) at this view and saves it as a reopenable session — swap in your own files afterwards.

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%3A260%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+read_marks%2Fchromosome

This session spec pastes after &session=spec- on any JBrowse Web instance.

{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr20",
      "tracks": [
        {
          "trackId": "na12878_chr20_pairs",
          "type": "LinearMarkDisplay",
          "height": 260
        }
      ]
    }
  ]
}

Loaded against config: https://jbrowse.org/demos/read_marks/config.json

The same view with jbrowse-anywidget. Point the adapter at your own file.

from jbrowse_anywidget import LinearGenomeView, fetch_hub

view = LinearGenomeView(
    assembly=fetch_hub("hg38"),
    location="chr20",
)

view.add_track({
    "type": "FeatureTrack",
    "trackId": "my_track",
    "name": "My track",
    "assemblyNames": ["hg38"],
    "adapter": {"type": "BamAdapter", "uri": "https://your-server/your-file"},
})

view  # display the widget

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": 260
        }
      ]
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/demos/read_marks/config.json \
  --session 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.

The bar at 34.2 Mb is a homozygous deletion; the one at 32.9 Mb is the intron above.

Checking the bars against the callset

Every deletion over 2 kb 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'

Pairs of 2 to 10 kb in the 100 kb window around each call:

position, chr20sizegenotypepairs 2 to 10 kb in the window
1.58 Mb33.1 kb0/1not in range
32.94 Mb3.9 kb0/118
34.23 Mb3.3 kb1/146
43.64 Mb2.7 kb0/110
43.85 Mb2.6 kb0/125
52.14 Mb2.1 kb1/132
54.03 Mb10.9 kb0/1not in range
55.86 Mb6.0 kb0/116

Every callset deletion in range is a bar, the two homozygous ones tallest. Five 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. Read the first window out of the file directly:

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" }'
windowmean depth
chr20:32,937,680-32,941,58315.6x
chr20:32,930,000-32,937,00034.1x

Around the call, 1,688 pairs sit at the library insert and 41 exceed 2 kb.

Reproduce it end to end

Every step above is wrapped in one script, build_read_marks.sh:

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 it builds the two tracks above 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

References