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.17. 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 and write each transcript's statistic into the GFF3 attribute column
  • load the GFF3 as a gene track
  • color each isoform by that attribute

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
# Debian/Ubuntu
sudo apt install tabix
# or conda, which has every tool above
conda install -c bioconda -c conda-forge htslib

The data

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

You don't need these files to follow the page: the track configs below read the hosted results. The build script fetches them to rebuild those results. The test step takes any transcript-by-sample count matrix with a sample table and a transcript-to-gene table.

What it isFile
Eight RSEM per-transcript quantification tables, the four muscle donors then the four liver donors: , , , , , ,ENCFF353NZM.tsv
ENCFF172SLW.tsv
ENCFF140GJI.tsv
ENCFF576DOG.tsv
ENCFF996LRE.tsv
ENCFF641ADT.tsv
ENCFF392VYD.tsv
ENCFF383KWZ.tsv
The GENCODE v29 annotation those quantifications were made againstgencode.v29.annotation.gff3.gz
Read by URL (no download needed)
What it isFile
The four coverage bigWigs the demo's track config loads, one donor per tissue, plus and minus strand: , ,ENCFF007ZBY.bigWig
ENCFF518WGP.bigWig
ENCFF565QRM.bigWig
ENCFF253OSP.bigWig

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 takes three inputs:

  • cnt, a transcript-by-sample count matrix
  • coldata, a table with one tissue row per sample
  • txinfo, a table with an isoform_id and a gene_id for each transcript

The test runs on them:

# 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, ]
txinfo <- txinfo[rownames(cnt), ]
multi <- names(which(table(txinfo$gene_id) > 1))
txinfo <- txinfo[txinfo$gene_id %in% multi, ]
cnt <- cnt[txinfo$isoform_id, ]

# 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.
  • 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).

A transcript is called when its regular FDR and its |ΔIF| pass the cutoffs in the code below. The script keeps each gene with a called transcript that is meaningfully expressed, since a fraction can swing widely on a handful of reads, and appends each transcript's numbers to its attribute column in the GENCODE v29 GFF3.

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 and appends:

  • dif, fdr, dtu and dif_called to each transcript row
  • 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

# satuRn writes NA for a transcript it could not fit
stats = {r['isoform_id']: r
         for r in csv.DictReader(open('results.tsv'), delimiter='\t')
         if 'NA' not in (r['dIF'], r['regular_FDR'])}

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

The statistics in the 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 two more keys, which the track's mouseover reads:

  • dtu_transcripts, the count of called isoforms
  • dtu_top_dif, the largest ΔIF among them

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, encoding.color is a color scale: its field, dif_called (a called transcript's isoform-fraction change), goes through the threshold cuts in scale.domain to the colors in scale.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 legend.title, plus a (no value) row for uncalled transcripts.
  • A UTR follows encoding.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",
    "encoding": {
      "color": {
        "field": "dif_called",
        "scale": {
          "type": "threshold",
          "domain": [-0.6, -0.3, 0, 0.3, 0.6],
          "range": [
            "#124f95",
            "#2370cc",
            "#6394d5",
            "#d5716a",
            "#c63335",
            "#901e21"
          ],
          "names": [
            "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"
          ]
        },
        "legend": { "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\",\"encoding\":{\"color\":{\"field\":\"dif_called\",\"scale\":{\"type\":\"threshold\",\"domain\":[-0.6,-0.3,0,0.3,0.6],\"range\":[\"#124f95\",\"#2370cc\",\"#6394d5\",\"#d5716a\",\"#c63335\",\"#901e21\"],\"names\":[\"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\"]},\"legend\":{\"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",
    "encoding": {
      "color": {
        "field": "dif_called",
        "scale": {
          "type": "threshold",
          "domain": [-0.6, -0.3, 0, 0.3, 0.6],
          "range": [
            "#124f95",
            "#2370cc",
            "#6394d5",
            "#d5716a",
            "#c63335",
            "#901e21"
          ],
          "names": [
            "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"
          ]
        },
        "legend": { "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 ↗

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 hold its reads.
  • satuRn used no genomic coordinates, so the coverage is an independent check on the color.

Each coverage track scales to its own peak by default. To compare the tissues by height, choose Y axis... → Share axis with on one track and tick the other; the two then share one axis that follows the view. 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": {
    "encoding": { "y": { "scale": { "group": "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 '{"encoding":{"y":{"scale":{"group":"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": {
    "encoding": { "y": { "scale": { "group": "coverage" } } }
  }
}

Open this track in JBrowse Web ↗

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": {
    "encoding": { "y": { "scale": { "group": "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 '{"encoding":{"y":{"scale":{"group":"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": {
    "encoding": { "y": { "scale": { "group": "coverage" } } }
  }
}

Open this track in JBrowse Web ↗

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": {
    "encoding": { "y": { "scale": { "group": "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 '{"encoding":{"y":{"scale":{"group":"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": {
    "encoding": { "y": { "scale": { "group": "coverage" } } }
  }
}

Open this track in JBrowse Web ↗

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": {
    "encoding": { "y": { "scale": { "group": "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 '{"encoding":{"y":{"scale":{"group":"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": {
    "encoding": { "y": { "scale": { "group": "coverage" } } }
  }
}

Open this track in JBrowse Web ↗

ATP5F1C with muscle and liver RNA-seq coverage on a shared scale, over GENCODE transcripts colored by satuRn's isoform-fraction change. The marked column is a cassette exon.
Open in JBrowse Web ↗

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",
          "height": 90
        },
        {
          "trackId": "liver_plus",
          "height": 90
        },
        {
          "trackId": "dtu_muscle_vs_liver",
          "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",
                    "height": 90
                },
                {
                    "trackId": "liver_plus",
                    "height": 90
                },
                {
                    "trackId": "dtu_muscle_vs_liver",
                    "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: a data file opens as the track its name says.

npx @jbrowse/img --config https://jbrowse.org/demos/dtu/config.json \
  --assembly hg38 --loc chr10:7,787,600..7,812,400 \
  --track muscle_plus '{"height":90}' \
  --track liver_plus '{"height":90}' \
  --track dtu_muscle_vs_liver '{"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",
          "height": 90
        },
        {
          "trackId": "liver_plus",
          "height": 90
        },
        {
          "trackId": "dtu_muscle_vs_liver",
          "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

build_dtu_demo.sh wraps every step above. The script:

  1. Downloads the ENTEx muscle and liver RSEM quantifications and GENCODE v29.
  2. Tests differential transcript usage with satuRn and computes each transcript's isoform-fraction change from TPM.
  3. Writes the statistics into the GENCODE GFF3, then bgzips and indexes it.
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. ↩