Skip to main content

Single-cell ATAC pseudobulk

TL;DR: pseudobulk outside JBrowse, pooling each cluster's cells into one coverage BigWig, then load the whole set as a single MultiWiggle track, which draws one row per file.

Prerequisites

  • cells already clustered and labeled: either a fragments file (or a barcoded BAM) plus a barcode-to-label table, or the project object your analysis tool already holds, an AnnData in SnapATAC2 (Python), an ArchRProject in ArchR, or a Seurat/Signac object in R
  • the pseudobulk tool that follows from whichever of those you have: pip install snapatac2, pip install deeptools sinto, or bedGraphToBigWig for the fragments-file route (ArchR and Signac install from R)
  • a JBrowse instance to load the finished BigWigs into (see the web quickstart, or the desktop quickstart, which opens local .bw files with nothing hosted)

Where the data comes from

SnapATAC2's annotated release of the 10x 5k PBMC scATAC dataset, already clustered and cell-type-labeled by that tool's own pipeline.

Pooling cells into rows

One ATAC cell contributes only a few thousand fragments, so a coverage track of a single cell is almost entirely zero. Pseudobulking pools every fragment belonging to a label into one profile, and each cell type comes out as a dense track resembling a bulk ATAC experiment on that cell type alone. JBrowse stacks the resulting files as rows of one track.

PBMC markers are the check: at a T-cell marker the T-cell rows carry the signal, and at a B-cell marker the B-cell rows light up.

The pseudobulk step runs in the same environment your clustering does, so the BigWigs it writes can also be viewed inline through the Python anywidget interface (or JBrowseR) without leaving the session.

Generating per-group BigWigs

Clustering and cell-type labeling stay upstream, in Cell Ranger ATAC, ArchR, Signac, or SnapATAC2. Two settings decide whether the rows this page draws can be compared to each other, whichever tool writes them:

  • Normalization. Groups differ in cell count and in total fragments, so each track needs normalizing (CPM / RPKM, or per-cell-count) for a peak's height to mean accessibility.
  • Bin size, which trades resolution against file size. Peak shape is the readable part of an ATAC track, so the bin has to stay well inside one peak; export_coverage below uses 25 bp.

SnapATAC2's export_coverage splits cells by a metadata column and writes one normalized BigWig per group in a single call, which covers the pseudobulk step for this dataset:

import snapatac2 as snap

# adata: an AnnData with fragments imported and a cell-type/cluster label in obs
snap.ex.export_coverage(
    adata,
    groupby="cell_type",     # column in adata.obs to split on
    bin_size=25,             # bp per bin
    normalization="RPKM",    # comparable across groups
    out_dir="bw",
    suffix=".bw",
    n_jobs=2,                # each worker holds a genome-wide coverage vector
    # blacklist= takes an ENCODE blacklist BED and drops those intervals from
    # every group. The build script does not pass it, so the figures below are
    # unmasked coverage.
)
# writes bw/<cell_type>.bw, one per group, keyed by group in the returned dict

n_jobs is a memory knob: each worker holds a whole genome-wide coverage vector, and the BigWig writer dies partway through the groups when memory runs out. Two workers fit this dataset.

groupby is the whole decision: pass the cluster column ("leiden") to get one row per cluster, or the annotated column ("cell_type") to get one row per cell type.

Other starting points

Every route ends at one .bw per cell type. The tools are linked under References:

  • An ArchRProject: getGroupBW(groupBy = "CellType", tileSize = 25) groups cells, sums their Tn5 insertions and writes one BigWig per group. normMethod = "ReadsInTSS" normalizes by signal-in-TSS, accounting for depth and data quality together; "nCells" and "nFrags" are the alternatives.
  • A barcoded BAM (Cell Ranger ATAC, or what a Signac workflow starts from): split it by label with sinto filterbarcodes, passing the barcode-to-label table and the barcode tag, then run deepTools bamCoverage on each with --binSize 25 --normalizeUsing CPM --extendReads. RPGC also needs --effectiveGenomeSize; CPM and RPKM do not.
  • A fragments.tsv.gz and nothing else: filter it to each cluster's barcodes, then bedtools genomecov -bg and bedGraphToBigWig per group. This route is unnormalized, so scale each group yourself (1e6 / total fragments for CPM) before the conversion.

Loading the BigWigs as a MultiWiggle track

In JBrowse, all the per-cell-type BigWigs go into one track: a MultiQuantitativeTrack whose MultiWiggleAdapter holds one BigWigAdapter per file. Each subadapter carries a name (the row label), an optional color, and an optional group.

assemblyNames names an assembly already configured in JBrowse, hg38 here, which is what the BigWigs above were built against. See the assemblies configuration guide if it is not set up yet. Minimal three-cell-type example:

{
  "type": "MultiQuantitativeTrack",
  "trackId": "scatac_pseudobulk",
  "name": "scATAC by cell type",
  "category": ["Single cell", "Chromatin accessibility"],
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "subadapters": [
      {
        "type": "BigWigAdapter",
        "name": "CD8 Naive",
        "group": "T cell",
        "color": "#4363d8",
        "uri": "https://example.com/bw/CD8_Naive.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "CD8 Memory",
        "group": "T cell",
        "color": "#3cb44b",
        "uri": "https://example.com/bw/CD8_Memory.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "Naive B",
        "group": "B cell",
        "color": "#f58231",
        "uri": "https://example.com/bw/Naive_B.bw"
      }
    ]
  }
}
jbrowse add-track-json '{
  "type": "MultiQuantitativeTrack",
  "trackId": "scatac_pseudobulk",
  "name": "scATAC by cell type",
  "category": ["Single cell", "Chromatin accessibility"],
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "subadapters": [
      {
        "type": "BigWigAdapter",
        "name": "CD8 Naive",
        "group": "T cell",
        "color": "#4363d8",
        "uri": "https://example.com/bw/CD8_Naive.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "CD8 Memory",
        "group": "T cell",
        "color": "#3cb44b",
        "uri": "https://example.com/bw/CD8_Memory.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "Naive B",
        "group": "B cell",
        "color": "#f58231",
        "uri": "https://example.com/bw/Naive_B.bw"
      }
    ]
  }
}'

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": "MultiQuantitativeTrack",
  "trackId": "scatac_pseudobulk",
  "name": "scATAC by cell type",
  "category": ["Single cell", "Chromatin accessibility"],
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "subadapters": [
      {
        "type": "BigWigAdapter",
        "name": "CD8 Naive",
        "group": "T cell",
        "color": "#4363d8",
        "uri": "https://example.com/bw/CD8_Naive.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "CD8 Memory",
        "group": "T cell",
        "color": "#3cb44b",
        "uri": "https://example.com/bw/CD8_Memory.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "Naive B",
        "group": "B cell",
        "color": "#f58231",
        "uri": "https://example.com/bw/Naive_B.bw"
      }
    ]
  }
}

Three things in that list are worth writing by hand:

  • Order. Subadapters draw in the order given, so list them grouped by lineage.
  • color. Take each row's from the cluster's color in your analysis, so a cell type is the same color in the browser as on the UMAP.
  • group. What the sidebar tree branches on, and what Clustering rows reorders.

If you don't need per-row names, colors or groups, the bigWigs shorthand takes a plain array of URLs and derives each row's label from its filename:

{
  "type": "MultiQuantitativeTrack",
  "trackId": "scatac_pseudobulk_simple",
  "name": "scATAC pseudobulk",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "bigWigs": [
      "https://example.com/bw/CD8_Naive.bw",
      "https://example.com/bw/CD8_Memory.bw",
      "https://example.com/bw/Naive_B.bw"
    ]
  }
}
jbrowse add-track-json '{
  "type": "MultiQuantitativeTrack",
  "trackId": "scatac_pseudobulk_simple",
  "name": "scATAC pseudobulk",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "bigWigs": [
      "https://example.com/bw/CD8_Naive.bw",
      "https://example.com/bw/CD8_Memory.bw",
      "https://example.com/bw/Naive_B.bw"
    ]
  }
}'

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": "MultiQuantitativeTrack",
  "trackId": "scatac_pseudobulk_simple",
  "name": "scATAC pseudobulk",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "bigWigs": [
      "https://example.com/bw/CD8_Naive.bw",
      "https://example.com/bw/CD8_Memory.bw",
      "https://example.com/bw/Naive_B.bw"
    ]
  }
}

A uri reaches anywhere, so a published atlas needs no pipeline at all: CATlas serves hg38 coverage from https://decoder-genetics.wustl.edu/catlasv1/humanenhancer/data/bw/, one file per cell type, and naming the ones you want is the whole track. Percent-encode the + in a cell-type name — T_lymphocyte_2_CD4%2B.bw — which is the one way those URLs go wrong quietly.

The display is a MultiLinearWiggleDisplay, and how the rows are drawn is one slot: defaultRendering lists every mode, and the track menu switches between them live. multirowxy (the default, and the figures on this page) is best for comparing peak shape; multirowdensity maps score to color, which fits more rows in the same space. User guide: Multi-quantitative track covers the rest of the menu.

Loaded, the twelve rows put the marker check in one frame:

Twelve per-cell-type BigWigs from the 10x 5k PBMC scATAC dataset, loaded as one MultiQuantitativeTrack, over CD8A and MS4A1 in one discontinuous view. CD8A is carried by the CD8, MAIT and NK rows; MS4A1 by the two B rows and nothing else. 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, loaded from a bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  2. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “NCBI RefSeq (GFF3Tabix)”.
  3. Track menu → Set feature height → Compact“NCBI RefSeq (GFF3Tabix)”.
  4. Track menu → Show only genes (checked)“NCBI RefSeq (GFF3Tabix)”.
  5. Drag the bar at the bottom edge of the track to resize it.“NCBI RefSeq (GFF3Tabix)”.
  6. Add your own track: File → Open track..., then paste a URL or choose a local file.This figure uses “Single-cell ATAC pseudobulk by cell type (10x 5k PBMC, SnapATAC2)”.
  7. Drag the bar at the bottom edge of the track to resize it.“Single-cell ATAC pseudobulk by cell type (10x 5k PBMC, SnapATAC2)”.
  8. Type your region of interest into the location box and press Enter.This figure is at chr2:86,780,000-86,820,000 chr11:60,450,000-60,490,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=test_data%2Fconfig_demo.json&session=spec-%7B%22views%22%3A%5B%7B%22assembly%22%3A%22hg38%22%2C%22loc%22%3A%22chr2%3A86%2C780%2C000-86%2C820%2C000+chr11%3A60%2C450%2C000-60%2C490%2C000%22%2C%22type%22%3A%22LinearGenomeView%22%2C%22tracks%22%3A%5B%7B%22trackId%22%3A%22ncbi_refseq_109_hg38_latest%22%2C%22type%22%3A%22LinearBasicDisplay%22%2C%22displayMode%22%3A%22compact%22%2C%22showOnlyGenes%22%3Atrue%2C%22height%22%3A80%7D%2C%7B%22trackId%22%3A%22pbmc5k_scatac_pseudobulk_hg38%22%2C%22type%22%3A%22MultiLinearWiggleDisplay%22%2C%22height%22%3A300%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+scatac%2Fpbmc5k_marker_swap

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

{
  "views": [
    {
      "assembly": "hg38",
      "loc": "chr2:86,780,000-86,820,000 chr11:60,450,000-60,490,000",
      "type": "LinearGenomeView",
      "tracks": [
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "type": "LinearBasicDisplay",
          "displayMode": "compact",
          "showOnlyGenes": true,
          "height": 80
        },
        {
          "trackId": "pbmc5k_scatac_pseudobulk_hg38",
          "type": "MultiLinearWiggleDisplay",
          "height": 300
        }
      ]
    }
  ]
}

Loaded against config: test_data/config_demo.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="chr2:86,780,000-86,820,000 chr11:60,450,000-60,490,000",
)

view.add_track({
    "type": "FeatureTrack",
    "trackId": "my_track_1",
    "name": "My track 1",
    "assemblyNames": ["hg38"],
    "adapter": {"type": "Gff3TabixAdapter", "uri": "https://your-server/your-file"},
})
view.add_track({
    "type": "MultiQuantitativeTrack",
    "trackId": "my_track_2",
    "name": "My track 2",
    "assemblyNames": ["hg38"],
    "adapter": {"type": "MultiWiggleAdapter", "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": [
    {
      "assembly": "hg38",
      "loc": "chr2:86,780,000-86,820,000 chr11:60,450,000-60,490,000",
      "type": "LinearGenomeView",
      "tracks": [
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "type": "LinearBasicDisplay",
          "displayMode": "compact",
          "showOnlyGenes": true,
          "height": 80
        },
        {
          "trackId": "pbmc5k_scatac_pseudobulk_hg38",
          "type": "MultiLinearWiggleDisplay",
          "height": 300
        }
      ]
    }
  ]
}
JSON

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

Building the subadapter list from files

Two workflows write the list for you, from a set of files rather than one entry at a time.

"Add multi-wiggle track", in the "Add track" workflow, takes the BigWig URLs one per line, or a JSON array of subadapter objects, and builds the MultiQuantitativeTrack from them. Exporting the session gets the JSON config back out. On JBrowse Desktop it reads the .bw files straight from local disk with no web server.

jbrowse add-track --multiwig takes the whole set of BigWigs in place of the usual single positional file. The row labels come from the filenames, which the pseudobulk step already named after the groups:

jbrowse add-track --multiwig "$(find bw -name '*.bw' | sort | paste -sd,)" \
  --name "scATAC by cell type" --assemblyNames hg38 \
  --load copy --subDir bw --out /var/www/html/jbrowse2

--load copy --subDir bw copies local files in beside config.json, and both drop out for BigWigs already served over HTTP. To carry per-row names, colors, and groups, pass a .json file of subadapter objects, the same objects as the config above, instead of the comma list.

Reproduce it end to end

One script runs the whole path, build_scatac_pseudobulk.sh:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_scatac_pseudobulk.sh
bash build_scatac_pseudobulk.sh    # builds ./scatac_pseudobulk_build
npx --yes serve scatac_pseudobulk_build/jbrowse2

Its input is SnapATAC2's annotated release of the 10x 5k-PBMC dataset, which is what that tool's standard pipeline and cell-type annotation tutorials already produce: fragments imported, cells QC-filtered and clustered, and each cluster labeled by transferring cell types from a matched multiome reference. That AnnData carries per-barcode fragments alongside an obs["cell_type"] call, which is the pair pseudobulking needs, so the script's own work is short:

  • export_coverage(groupby="cell_type", bin_size=25, normalization="RPKM"), which writes one BigWig per cell type into bw/
  • a sources.json of subadapters, taking each row's color from the same object and its group and position from a lineage map the script states outright. Running it on your own experiment means replacing that map
  • jbrowse create plus add-assembly for hg38 and a RefSeq gene track, then the one MultiQuantitativeTrack those subadapters make up

Navigate the finished instance to the two markers in the figure and read the rows against the labels. Rows that stay open everywhere usually mean the normalization step was skipped, since an unnormalized group's height tracks its cell count.

See also

References

Pseudobulk / coverage tools:

Reference datasets: