Skip to main content

Differential transcript usage

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.

Differential transcript usage tests whether a gene changes which of its isoforms it expresses between two conditions, here skeletal muscle and liver. We run the test with satuRn, a Bioconductor package, on ENCODE quantifications, write each transcript's statistic into the GFF3 attribute column, and color the gene track by it.

Prerequisites

  • a JBrowse to open the finished GFF3 and coverage tracks in: Desktop takes a local file by path, Web through Add track
  • curl
  • python3
  • htslib (bgzip, tabix)
  • R with satuRn, SummarizedExperiment, edgeR and limma, all Bioconductor

Where the data comes from

ENCODE's ENTEx panel, four skeletal-muscle and four liver donors, quantified with RSEM against GENCODE v29.

The build script fetches these files, so there is nothing to download by hand.

Read by URL (no download needed)

Building the GFF3

The build takes the ENCODE quantifications to a colorable GFF3 in three steps, and one script runs all three.

Build the matrices. The script downloads the eight RSEM tables listed above, with the accessions written into it, and writes a count matrix, which feeds the model, and a TPM matrix, which feeds the effect size.

Test usage. satuRn fits a quasi-binomial model to each transcript's share of its gene's reads and tests that share between the two tissues. It needs a transcript-by-sample count matrix cnt, a coldata table with one tissue row per sample, and a txinfo table with an isoform_id and a gene_id for each transcript:

# after filterByExpr, any gene left with one isoform goes too: usage is a
# within-gene proportion, so a lone isoform is always 100%
keep <- edgeR::filterByExpr(cnt, group = coldata$tissue)
cnt <- cnt[keep, ]
multi <- names(which(table(txinfo$gene_id) > 1))
cnt <- cnt[txinfo$isoform_id[txinfo$gene_id %in% multi], ]

# rowData has to carry isoform_id and gene_id: satuRn reads each transcript's
# gene from there, since usage is a share of the gene's reads
se <- SummarizedExperiment(
  assays = list(counts = cnt),
  colData = coldata,
  rowData = txinfo
)

# the formula names the colData column holding the groups. 0 + tissue drops the
# intercept, so each tissue gets its own coefficient and the contrast below is
# the difference between two of them
se <- satuRn::fitDTU(object = se, formula = ~ 0 + tissue, parallel = FALSE)

design <- model.matrix(~ 0 + tissue, data = coldata)
colnames(design) <- levels(factor(coldata$tissue))
L <- limma::makeContrasts(muscle_vs_liver = muscle - liver, levels = design)

# sort = FALSE leaves the result rows in the order the assay had them, which is
# what lets the isoform fractions be indexed straight into the result
se <- satuRn::testDTU(object = se, contrasts = L, sort = FALSE)
res <- rowData(se)[["fitDTUResult_muscle_vs_liver"]]

res holds, per transcript, the p-value, the regular and empirical FDRs and the model estimates. The script computes each transcript's isoform fraction, its share of the gene's expression, from the TPM matrix.1

Write the statistics into GENCODE. The annotation has to be the release the quantifications were made against. RSEM names each transcript with its version (ENST00000356708.11) and GENCODE raises that version whenever it revises the transcript, so joining against a later release drops every transcript revised since, with no error. ENCODE lists the release on each quantification's file page as its genome annotation (V29 for all eight here).

The script keeps each gene with a called transcript that is meaningfully expressed, since a fraction can swing widely on a handful of reads. It subsets those genes out of the GENCODE v29 GFF3 and appends each transcript's numbers to its attribute column.

To do the same join over your own results, write results.tsv with one row per transcript and the columns isoform_id, regular_FDR and dIF (the isoform-fraction change). This version keeps every GENCODE row, appending dif, fdr, dtu and dif_called to each transcript row and dtu_transcripts and dtu_top_dif to the gene row of each gene with a called transcript. The track reads those keys:

import csv
import gzip

stats = {r['isoform_id']: r
         for r in csv.DictReader(open('results.tsv'), delimiter='\t')}

genes = {}
called = {}
records = []
with gzip.open('gencode.v29.annotation.gff3.gz', 'rt') as fh:
    for line in fh:
        if line.startswith('#'):
            continue
        cols = line.rstrip('\n').split('\t')
        d = dict(kv.split('=', 1) for kv in cols[8].split(';') if '=' in kv)
        if cols[2] == 'gene':
            genes[d['ID']] = cols
        elif cols[2] == 'transcript':
            r = stats.get(d['transcript_id'])
            if r:
                dif, fdr = float(r['dIF']), float(r['regular_FDR'])
                direction = 'ns'
                if fdr < 0.05 and abs(dif) > 0.1:
                    direction = 'muscle' if dif > 0 else 'liver'
                cols[8] += f';dif={dif:.3f};fdr={fdr:.3g};dtu={direction}'
                if direction != 'ns':
                    cols[8] += f';dif_called={dif:.3f}'
                    called.setdefault(d['Parent'], []).append(dif)
        records.append((cols[0], int(cols[3]), int(cols[4]), cols))

for gene_id, difs in called.items():
    top = max(difs, key=abs)
    genes[gene_id][8] += f';dtu_transcripts={len(difs)};dtu_top_dif={top:.3f}'

records.sort(key=lambda r: r[:3])
with open('dtu_muscle_vs_liver.gff3', 'w') as out:
    out.write('##gff-version 3\n')
    for _, _, _, cols in records:
        out.write('\t'.join(cols) + '\n')

The rows come out in coordinate order, so indexing is the ordinary pair:

bgzip -f dtu_muscle_vs_liver.gff3
tabix -f -p gff dtu_muscle_vs_liver.gff3.gz

satuRn's statistics in the GFF3 attribute column

The track reads each transcript's statistics from the attribute column. A transcript row from the finished file, wrapped:

chr10  HAVANA  transcript  7788129  7807815  .  +  .
  ID=ENST00000356708.11;Parent=ENSG00000165629.19;gene_name=ATP5F1C;
  transcript_name=ATP5F1C-202;...;
  dif=-0.299;fdr=0.0022;if_muscle=0.075;if_liver=0.375;
  tpm_muscle=10.03;tpm_liver=28.88;dtu=liver;dif_called=-0.299

The transcript statistics sit on the transcript row, and the exons, CDS and UTRs take the transcript value. The gene row has dtu_transcripts, the count of called isoforms, and dtu_top_dif, the largest ΔIF among them, which the track's mouseover reads. The keys are lowercase because the GFF parser lowercases them, so a color field named dIF reads nothing and paints every transcript grey.

dtu is muscle, liver or ns, as the script above calls it. dif_called is dif on the called transcripts and absent on the rest, so a transcript the test could not separate has no value to color and stays grey.

Loading GRCh38

GENCODE v29 is a GRCh38 annotation and the coverage is on GRCh38, so we load the hg38 assembly before either track.

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

Coloring each isoform by its usage change

In grammar-of-graphics terms, color is a color scale: its field, dif_called (a called transcript's isoform-fraction change), goes through the threshold cuts in domain to the colors in range, liver-preferred below zero and muscle-preferred above. A value on a cut takes the interval above it.

  • The key lists every interval under the title, plus a (no value) row for uncalled transcripts.
  • A UTR follows color unless utrColor is set.
  • Clicking an isoform opens its numbers in the details panel.

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

{
  "trackId": "dtu_muscle_vs_liver",
  "name": "Transcript usage: skeletal muscle vs liver (satuRn)",
  "uri": "https://jbrowse.org/demos/dtu/dtu_muscle_vs_liver.gff3.gz",
  "assemblyNames": ["hg38"],
  "displayDefaults": {
    "subfeatureLabels": "below",
    "color": {
      "field": "dif_called",
      "scale": "threshold",
      "domain": ["-0.6", "-0.3", "0", "0.3", "0.6"],
      "range": [
        "#124f95",
        "#2370cc",
        "#6394d5",
        "#d5716a",
        "#c63335",
        "#901e21"
      ],
      "labels": [
        "liver, below -0.6",
        "liver, -0.6 to -0.3",
        "liver, -0.3 to 0",
        "muscle, 0 to 0.3",
        "muscle, 0.3 to 0.6",
        "muscle, 0.6 and above"
      ],
      "title": "ΔIF, muscle - liver"
    },
    "labels": {
      "name": "jexl:feature.transcript_name||feature.gene_name||feature.name||feature.id"
    },
    "mouseover": "jexl:feature.gene_name+': '+feature.dtu_transcripts+' isoform(s) with a usage shift, largest ΔIF '+feature.dtu_top_dif"
  }
}
jbrowse add-track https://jbrowse.org/demos/dtu/dtu_muscle_vs_liver.gff3.gz \
  --trackId dtu_muscle_vs_liver \
  --name "Transcript usage: skeletal muscle vs liver (satuRn)" \
  --assemblyNames hg38 \
  --displayDefaults "{\"subfeatureLabels\":\"below\",\"color\":{\"field\":\"dif_called\",\"scale\":\"threshold\",\"domain\":[\"-0.6\",\"-0.3\",\"0\",\"0.3\",\"0.6\"],\"range\":[\"#124f95\",\"#2370cc\",\"#6394d5\",\"#d5716a\",\"#c63335\",\"#901e21\"],\"labels\":[\"liver, below -0.6\",\"liver, -0.6 to -0.3\",\"liver, -0.3 to 0\",\"muscle, 0 to 0.3\",\"muscle, 0.3 to 0.6\",\"muscle, 0.6 and above\"],\"title\":\"ΔIF, muscle - liver\"},\"labels\":{\"name\":\"jexl:feature.transcript_name||feature.gene_name||feature.name||feature.id\"},\"mouseover\":\"jexl:feature.gene_name+': '+feature.dtu_transcripts+' isoform(s) with a usage shift, largest ΔIF '+feature.dtu_top_dif\"}"

In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add track from pasted JSON, and paste:

{
  "trackId": "dtu_muscle_vs_liver",
  "name": "Transcript usage: skeletal muscle vs liver (satuRn)",
  "uri": "https://jbrowse.org/demos/dtu/dtu_muscle_vs_liver.gff3.gz",
  "assemblyNames": ["hg38"],
  "displayDefaults": {
    "subfeatureLabels": "below",
    "color": {
      "field": "dif_called",
      "scale": "threshold",
      "domain": ["-0.6", "-0.3", "0", "0.3", "0.6"],
      "range": [
        "#124f95",
        "#2370cc",
        "#6394d5",
        "#d5716a",
        "#c63335",
        "#901e21"
      ],
      "labels": [
        "liver, below -0.6",
        "liver, -0.6 to -0.3",
        "liver, -0.3 to 0",
        "muscle, 0 to 0.3",
        "muscle, 0.3 to 0.6",
        "muscle, 0.6 and above"
      ],
      "title": "ΔIF, muscle - liver"
    },
    "labels": {
      "name": "jexl:feature.transcript_name||feature.gene_name||feature.name||feature.id"
    },
    "mouseover": "jexl:feature.gene_name+': '+feature.dtu_transcripts+' isoform(s) with a usage shift, largest ΔIF '+feature.dtu_top_dif"
  }
}

Open this track in JBrowse Web ↗

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

Add the four coverage tracks below and open them with the transcript track at chr10:7,787,600-7,812,400, the whole of ATP5F1C, an ATP synthase subunit. The gene is on the plus strand, so the two plus-strand tracks have its reads. satuRn used no genomic coordinates, so the coverage is an independent check on the color.

Each coverage track scales to its peak until the tracks share an axis.

Y axis... → Share axis with on one track, with the other ticked, gives both one axis that follows the view, so the two tissues compare by height. In a config, each track names the same group:

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

{
  "trackId": "liver_plus",
  "name": "Liver RNA-seq, + strand (ENCSR135IAL)",
  "uri": "https://jbrowse.org/demos/dtu/ENCFF565QRM.liver.plus.bigWig",
  "assemblyNames": ["hg38"],
  "displayDefaults": {
    "scales": { "y": { "autoscaleGroup": "coverage" } }
  }
}
jbrowse add-track https://jbrowse.org/demos/dtu/ENCFF565QRM.liver.plus.bigWig \
  --trackId liver_plus \
  --name "Liver RNA-seq, + strand (ENCSR135IAL)" \
  --assemblyNames hg38 \
  --displayDefaults '{"scales":{"y":{"autoscaleGroup":"coverage"}}}'

In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add track from pasted JSON, and paste:

{
  "trackId": "liver_plus",
  "name": "Liver RNA-seq, + strand (ENCSR135IAL)",
  "uri": "https://jbrowse.org/demos/dtu/ENCFF565QRM.liver.plus.bigWig",
  "assemblyNames": ["hg38"],
  "displayDefaults": {
    "scales": { "y": { "autoscaleGroup": "coverage" } }
  }
}

Open this track in JBrowse Web ↗

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

The other three tracks differ in their id, name and file:

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

{
  "trackId": "liver_minus",
  "name": "Liver RNA-seq, - strand (ENCSR135IAL)",
  "uri": "https://jbrowse.org/demos/dtu/ENCFF253OSP.liver.minus.bigWig",
  "assemblyNames": ["hg38"],
  "displayDefaults": {
    "scales": { "y": { "autoscaleGroup": "coverage" } }
  }
}
jbrowse add-track https://jbrowse.org/demos/dtu/ENCFF253OSP.liver.minus.bigWig \
  --trackId liver_minus \
  --name "Liver RNA-seq, - strand (ENCSR135IAL)" \
  --assemblyNames hg38 \
  --displayDefaults '{"scales":{"y":{"autoscaleGroup":"coverage"}}}'

In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add track from pasted JSON, and paste:

{
  "trackId": "liver_minus",
  "name": "Liver RNA-seq, - strand (ENCSR135IAL)",
  "uri": "https://jbrowse.org/demos/dtu/ENCFF253OSP.liver.minus.bigWig",
  "assemblyNames": ["hg38"],
  "displayDefaults": {
    "scales": { "y": { "autoscaleGroup": "coverage" } }
  }
}

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.

{
  "trackId": "muscle_plus",
  "name": "Skeletal muscle RNA-seq, + strand (ENCSR609NZM)",
  "uri": "https://jbrowse.org/demos/dtu/ENCFF007ZBY.muscle.plus.bigWig",
  "assemblyNames": ["hg38"],
  "displayDefaults": {
    "scales": { "y": { "autoscaleGroup": "coverage" } }
  }
}
jbrowse add-track https://jbrowse.org/demos/dtu/ENCFF007ZBY.muscle.plus.bigWig \
  --trackId muscle_plus \
  --name "Skeletal muscle RNA-seq, + strand (ENCSR609NZM)" \
  --assemblyNames hg38 \
  --displayDefaults '{"scales":{"y":{"autoscaleGroup":"coverage"}}}'

In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add track from pasted JSON, and paste:

{
  "trackId": "muscle_plus",
  "name": "Skeletal muscle RNA-seq, + strand (ENCSR609NZM)",
  "uri": "https://jbrowse.org/demos/dtu/ENCFF007ZBY.muscle.plus.bigWig",
  "assemblyNames": ["hg38"],
  "displayDefaults": {
    "scales": { "y": { "autoscaleGroup": "coverage" } }
  }
}

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.

{
  "trackId": "muscle_minus",
  "name": "Skeletal muscle RNA-seq, - strand (ENCSR609NZM)",
  "uri": "https://jbrowse.org/demos/dtu/ENCFF518WGP.muscle.minus.bigWig",
  "assemblyNames": ["hg38"],
  "displayDefaults": {
    "scales": { "y": { "autoscaleGroup": "coverage" } }
  }
}
jbrowse add-track https://jbrowse.org/demos/dtu/ENCFF518WGP.muscle.minus.bigWig \
  --trackId muscle_minus \
  --name "Skeletal muscle RNA-seq, - strand (ENCSR609NZM)" \
  --assemblyNames hg38 \
  --displayDefaults '{"scales":{"y":{"autoscaleGroup":"coverage"}}}'

In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add track from pasted JSON, and paste:

{
  "trackId": "muscle_minus",
  "name": "Skeletal muscle RNA-seq, - strand (ENCSR609NZM)",
  "uri": "https://jbrowse.org/demos/dtu/ENCFF518WGP.muscle.minus.bigWig",
  "assemblyNames": ["hg38"],
  "displayDefaults": {
    "scales": { "y": { "autoscaleGroup": "coverage" } }
  }
}

Open this track in JBrowse Web ↗

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

ATP5F1C on hg38. ENCODE skeletal-muscle and liver RNA-seq coverage on a shared scale, over GENCODE transcripts colored by the isoform-fraction change satuRn measured between the two tissues. The marked column is a cassette exon, included in some isoforms and skipped in others. Muscle coverage is flat there and liver peaks.
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%2Fdtu%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%22chr10%3A7%2C787%2C600..7%2C812%2C400%22%2C%22tracks%22%3A%5B%7B%22trackId%22%3A%22muscle_plus%22%2C%22type%22%3A%22LinearWiggleDisplay%22%2C%22height%22%3A90%7D%2C%7B%22trackId%22%3A%22liver_plus%22%2C%22type%22%3A%22LinearWiggleDisplay%22%2C%22height%22%3A90%7D%2C%7B%22trackId%22%3A%22dtu_muscle_vs_liver%22%2C%22type%22%3A%22LinearBasicDisplay%22%2C%22heightMode%22%3A%22grow%22%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+dtu%2Fdtu_colored_gene_glyph

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 BigWig file (.bw).This figure uses “Skeletal muscle RNA-seq, + strand (ENCSR609NZM)”.
  3. Add the next track the same way. This one needs a BigWig file (.bw).This figure uses “Liver RNA-seq, + strand (ENCSR135IAL)”.
  4. Add the next track the same way. This one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “Transcript usage: skeletal muscle vs liver (satuRn)”.
    1. Track menu → Feature height → Track sizing → Fixed feature height + autogrow track height
  5. Type your region of interest into the location box and press Enter.This figure is at chr10:7,787,600..7,812,400.

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

{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr10:7,787,600..7,812,400",
      "tracks": [
        {
          "trackId": "muscle_plus",
          "type": "LinearWiggleDisplay",
          "height": 90
        },
        {
          "trackId": "liver_plus",
          "type": "LinearWiggleDisplay",
          "height": 90
        },
        {
          "trackId": "dtu_muscle_vs_liver",
          "type": "LinearBasicDisplay",
          "heightMode": "grow"
        }
      ]
    }
  ]
}

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/dtu/config.json")
app = JBrowseApp(
    assemblies=config["assemblies"],
    tracks=config.get("tracks", []),
    views=[
        {
            "type": "LinearGenomeView",
            "assembly": "hg38",
            "loc": "chr10:7,787,600..7,812,400",
            "tracks": [
                {
                    "trackId": "muscle_plus",
                    "type": "LinearWiggleDisplay",
                    "height": 90
                },
                {
                    "trackId": "liver_plus",
                    "type": "LinearWiggleDisplay",
                    "height": 90
                },
                {
                    "trackId": "dtu_muscle_vs_liver",
                    "type": "LinearBasicDisplay",
                    "heightMode": "grow"
                }
            ]
        }
    ],
)
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/dtu/config.json \
  --assembly hg38 --loc chr10:7,787,600..7,812,400 \
  --track muscle_plus display:LinearWiggleDisplay '{"height":90}' \
  --track liver_plus display:LinearWiggleDisplay '{"height":90}' \
  --track dtu_muscle_vs_liver display:LinearBasicDisplay '{"heightMode":"grow"}' \
  --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": "chr10:7,787,600..7,812,400",
      "tracks": [
        {
          "trackId": "muscle_plus",
          "type": "LinearWiggleDisplay",
          "height": 90
        },
        {
          "trackId": "liver_plus",
          "type": "LinearWiggleDisplay",
          "height": 90
        },
        {
          "trackId": "dtu_muscle_vs_liver",
          "type": "LinearBasicDisplay",
          "heightMode": "grow"
        }
      ]
    }
  ]
}
JSON

cat > callouts.json <<'JSON'
[
  {
    "type": "box",
    "anchor": {
      "loc": "chr10:7,806,974-7,807,010"
    },
    "pad": 9,
    "strokeWidth": 3
  },
  {
    "type": "text",
    "text": "exon skipped",
    "anchor": {
      "loc": "chr10:7,806,974-7,807,010",
      "trackId": "muscle_plus",
      "fracY": 0.5,
      "alignX": "left"
    },
    "textAlign": "end",
    "dx": -14
  },
  {
    "type": "text",
    "text": "exon retained",
    "anchor": {
      "loc": "chr10:7,806,974-7,807,010",
      "trackId": "liver_plus",
      "fracY": 0.35,
      "alignX": "left"
    },
    "textAlign": "end",
    "dx": -14
  }
]
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/demos/dtu/config.json \
  --width 1500 --height 730 \
  --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.

Reproduce it end to end

One script wraps every step above, build_dtu_demo.sh:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_dtu_demo.sh
bash build_dtu_demo.sh dtu_build   # writes ./dtu_build/

The script needs the Prerequisites on your PATH and writes dtu_muscle_vs_liver.gff3.gz with its .tbi index, a local build of the file the track configuration above loads from jbrowse.org. Point the track's uri at the local copy to open your own run.

See also

Citations

Notes

  1. The script takes the isoform fraction from TPM, because read counts scale with effective length and bias a count-based fraction toward long isoforms. It gates on satuRn's regular FDR, because the empirical FDR assumes most tests are null and this contrast breaks that assumption. ↩