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 single cell's coverage track is almost entirely zero. Pseudobulking pools every fragment belonging to a label into one profile, a dense track resembling a bulk ATAC experiment on that cell type. JBrowse stacks the 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 BigWigs can also be viewed inline from the clustering environment through the Python anywidget interface or JBrowseR.

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 can be compared, whichever tool writes them:

  • Normalization. Groups differ in cell count and 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. 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:

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 genome-wide coverage vector, and the writer dies partway through the groups when memory runs out. groupby picks the rows: the cluster column ("leiden") for one per cluster, or the annotated column ("cell_type") for one 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) writes one BigWig per group. normMethod = "ReadsInTSS" normalizes by signal-in-TSS; "nCells" and "nFrags" are the alternatives
  • A barcoded BAM (Cell Ranger ATAC, or a Signac workflow's input): split it by label with sinto filterbarcodes, then run deepTools bamCoverage on each with --binSize 25 --normalizeUsing CPM --extendReads. RPGC also needs --effectiveGenomeSize
  • A fragments.tsv.gz alone: 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

All the per-cell-type BigWigs go into one MultiQuantitativeTrack whose MultiWiggleAdapter holds one BigWigAdapter per file, each with a name, an optional color, and an optional group. assemblyNames names an assembly already configured (see the assemblies configuration guide). Three cell types:

{
  "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 group them by lineage
  • color. Take each row's from the cluster's color in your analysis, so a cell type matches its UMAP color
  • group. What the sidebar tree branches on, and what Clustering rows reorders

Without per-row names, colors or groups, the bigWigs shorthand takes a plain array of URLs and labels each row 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 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. Percent-encode the + in a cell-type name (T_lymphocyte_2_CD4%2B.bw), the one way those URLs go wrong quietly.

defaultRendering lists every drawing mode, and the track menu switches between them live. multirowxy (the default, and the figures here) compares peak shape; multirowdensity maps score to color and fits more rows. User guide: Multi-quantitative track covers the rest of the menu.

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 from a set of files.

"Add multi-wiggle track", in the "Add track" workflow, takes BigWig URLs one per line, or a JSON array of subadapter objects. Exporting the session gets the JSON config back out. On JBrowse Desktop it reads local .bw files directly.

jbrowse add-track --multiwig takes the whole set of BigWigs in place of the single positional file, labeling rows from the filenames:

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; both drop out for BigWigs already served over HTTP. For per-row names, colors and groups, pass a .json file of subadapter objects 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, what that tool's standard pipeline and cell-type annotation tutorials produce: per-barcode fragments alongside an obs["cell_type"] call. The script's own work:

  • export_coverage(groupby="cell_type", bin_size=25, normalization="RPKM"), 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

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: