Skip to main content

Differential transcript usage

TL;DR: a per-transcript statistic goes into the GFF3 attribute column, and a jexl: color callback on the gene track paints it. This page builds that GFF3 from ENCODE quantifications and gives the track configuration that reads it.

Prerequisites

  • curl
  • python3
  • htslib (bgzip, tabix)
  • R with satuRn, SummarizedExperiment, edgeR and limma, all Bioconductor
  • to read along without running anything, the finished analysis is hosted at jbrowse.org/demos/dtu

Where the data comes from

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

Building the GFF3

Four steps take the ENCODE quantifications to a GFF3 the gene glyph can paint. One script runs all four.

Fetch the quantifications. Eight RSEM per-transcript tables from ENCODE's ENTEx panel, skeletal muscle and liver, four donors each, quantified against GENCODE v29. The accessions are written into the script.

Build the matrices. One pass over those tables writes a count matrix and a TPM matrix: counts feed the model, TPM 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, over a counts matrix and a tissue column:

# 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 to know whose proportion the transcript is a share of
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 a
# plain difference between two of them rather than a difference of differences
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 carries the p-value, both FDRs and the model's own estimates per transcript. The isoform fractions the color reads come from the TPM matrix rather than from this table.

Write the statistics into GENCODE. The called genes are subset out of the GENCODE v29 GFF3 and each transcript's numbers are appended to its attribute column. 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 attribute column

This is the part the track configuration depends on. 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

Three properties of that line, each of which fails without an error:

  • keys are lowercase: the GFF parser lowercases them, so dIF= read back as feature.dIF is undefined, and an undefined branch returns the default color
  • values are strings, so numeric comparison requires parseFloat
  • the numbers are on the transcript row and nothing below it: the glyph evaluates the color against the box it paints, so the callback reaches up with feature.parent.dif. feature.dif reads the exon and paints the default

dtu is a flag with the values muscle, liver and ns, set by the same threshold the script reports on. The color branches on it before reading dif, so transcripts the test could not separate stay neutral.

The effect size and the FDR gate

Effect size from TPM, model fit on counts. Isoform fraction is a molar quantity, and read counts scale with abundance times effective length, so a count-based fraction is biased toward long isoforms.

The gate is satuRn's regular FDR. Its empirical FDR assumes most tests are null, which does not hold for this contrast: locfdr reports a misfit, and no transcript passes it. The script prints the minimum empirical FDR beside its count.

Configuring the track

One expression covers the whole transcript: a UTR follows color unless utrColor claims it. labels.name reads GENCODE's transcript_name, which also names the isoform under the cursor, and legend declares what the ramp means. mouseover resolves against the gene, so it summarizes the gene; an isoform's own numbers are in the details panel, one click away.

{
  "type": "FeatureTrack",
  "trackId": "dtu_muscle_vs_liver",
  "name": "Transcript usage: skeletal muscle vs liver (satuRn)",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "Gff3TabixAdapter",
    "uri": "https://jbrowse.org/demos/dtu/dtu_muscle_vs_liver.gff3.gz"
  },
  "displayDefaults": {
    "subfeatureLabels": "below",
    "color": "jexl:feature.parent.dtu=='muscle'?(parseFloat(feature.parent.dif)>0.6?'#901e21':parseFloat(feature.parent.dif)>0.3?'#c63335':'#d5716a'):feature.parent.dtu=='liver'?(parseFloat(feature.parent.dif)<-0.6?'#124f95':parseFloat(feature.parent.dif)<-0.3?'#2370cc':'#6394d5'):'#b2b1ac'",
    "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",
    "legend": [
      { "label": "muscle-preferred, ΔIF > 0.6", "color": "#901e21" },
      { "label": "muscle-preferred, ΔIF 0.3–0.6", "color": "#c63335" },
      { "label": "muscle-preferred, ΔIF 0.1–0.3", "color": "#d5716a" },
      { "label": "no usage shift (FDR ≥ 0.05)", "color": "#b2b1ac" },
      { "label": "liver-preferred, ΔIF 0.1–0.3", "color": "#6394d5" },
      { "label": "liver-preferred, ΔIF 0.3–0.6", "color": "#2370cc" },
      { "label": "liver-preferred, ΔIF > 0.6", "color": "#124f95" }
    ]
  }
}
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\":\"jexl:feature.parent.dtu=='muscle'?(parseFloat(feature.parent.dif)>0.6?'#901e21':parseFloat(feature.parent.dif)>0.3?'#c63335':'#d5716a'):feature.parent.dtu=='liver'?(parseFloat(feature.parent.dif)<-0.6?'#124f95':parseFloat(feature.parent.dif)<-0.3?'#2370cc':'#6394d5'):'#b2b1ac'\",\"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\",\"legend\":[{\"label\":\"muscle-preferred, ΔIF > 0.6\",\"color\":\"#901e21\"},{\"label\":\"muscle-preferred, ΔIF 0.3–0.6\",\"color\":\"#c63335\"},{\"label\":\"muscle-preferred, ΔIF 0.1–0.3\",\"color\":\"#d5716a\"},{\"label\":\"no usage shift (FDR ≥ 0.05)\",\"color\":\"#b2b1ac\"},{\"label\":\"liver-preferred, ΔIF 0.1–0.3\",\"color\":\"#6394d5\"},{\"label\":\"liver-preferred, ΔIF 0.3–0.6\",\"color\":\"#2370cc\"},{\"label\":\"liver-preferred, ΔIF > 0.6\",\"color\":\"#124f95\"}]}"

The track loads over the two coverage tracks at ATP5F1C. satuRn used no genomic coordinates, so the coverage lanes are an independent check on the color.

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 the cassette exon, where the muscle lane is flat and the liver lane peaks. Open in JBrowse: Open this view ↗

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 “muscle_plus” track.
  3. Drag the bar at the bottom edge of the track to resize it.“muscle_plus”.
  4. Track menu → Score → Set min/max score...Sets the score-axis minimum (0 here).“muscle_plus”.
  5. Track menu → Score → Set min/max score...Sets the score-axis maximum (16 here).“muscle_plus”.
  6. Add your own track: File → Open track..., then paste a URL or choose a local file.This figure uses the “liver_plus” track.
  7. Drag the bar at the bottom edge of the track to resize it.“liver_plus”.
  8. Track menu → Score → Set min/max score...Sets the score-axis minimum (0 here).“liver_plus”.
  9. Track menu → Score → Set min/max score...Sets the score-axis maximum (16 here).“liver_plus”.
  10. Add your own track: File → Open track..., then paste a URL or choose a local file.This figure uses the “dtu_muscle_vs_liver” track.
  11. Track menu → Set feature height → Track sizing → Fixed feature height + autogrow track height“dtu_muscle_vs_liver”.
  12. Type your region of interest into the location box and press Enter.This figure is at chr10:7,787,600..7,812,400.

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%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%2C%22minScore%22%3A0%2C%22maxScore%22%3A16%7D%2C%7B%22trackId%22%3A%22liver_plus%22%2C%22type%22%3A%22LinearWiggleDisplay%22%2C%22height%22%3A90%2C%22minScore%22%3A0%2C%22maxScore%22%3A16%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

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

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

Loaded against config: https://jbrowse.org/demos/dtu/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="chr10:7,787,600..7,812,400",
)

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.add_track({
    "type": "FeatureTrack",
    "trackId": "my_track_3",
    "name": "My track 3",
    "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": "chr10:7,787,600..7,812,400",
      "tracks": [
        {
          "trackId": "muscle_plus",
          "type": "LinearWiggleDisplay",
          "height": 90,
          "minScore": 0,
          "maxScore": 16
        },
        {
          "trackId": "liver_plus",
          "type": "LinearWiggleDisplay",
          "height": 90,
          "minScore": 0,
          "maxScore": 16
        },
        {
          "trackId": "dtu_muscle_vs_liver",
          "type": "LinearBasicDisplay",
          "heightMode": "grow"
        }
      ]
    }
  ]
}
JSON

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

Reproduce it end to end

Every step above is wrapped in one script, 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/

It fetches the eight RSEM tables and the four coverage bigWigs from ENCODE, downloads the GENCODE v29 GFF3 those quantifications were made against, runs the satuRn fit, and writes dtu_muscle_vs_liver.gff3.gz with its .tbi index: the local build of the file the track configuration above loads from jbrowse.org. Point the adapter's uri at the local copy to open your own run instead. It needs Prerequisites on your PATH.

Along the way it prints the transcript and gene counts at each filtering step, and the minimum empirical FDR beside the regular-FDR count.

See also

References