Skip to main content

Single-cell ATAC pseudobulk

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.

Single-cell ATAC-seq measures chromatin accessibility one cell at a time, and pooling the cells of each cluster, called pseudobulking, gives one accessibility profile per cell type. We pool a clustered 10x peripheral blood (PBMC) dataset into one coverage BigWig per cell type outside JBrowse, load the set as one multi-wiggle track with a row per file, and check the rows at T-cell and B-cell marker genes.

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 matches what 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)

Where the data comes from

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

The build script downloads the dataset through SnapATAC2, so there is nothing to download by hand.

Read by URL (no download needed)

Pooling cells into one coverage track per cell type

One ATAC cell contributes only a few thousand fragments, so a single cell's coverage is almost entirely zero. JBrowse stacks the pooled files as rows of one track.1

PBMC marker genes are the control. At a T-cell marker the T-cell rows have signal and the B-cell rows stay flat, and at a B-cell marker the reverse.

Generating one pseudobulk BigWig per cell type

Clustering and cell-type labeling happen upstream, in Cell Ranger ATAC, ArchR, Signac, or SnapATAC2. Whichever tool writes the files, two settings decide whether the rows compare:

  • 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 trades resolution against file size and 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

Too many n_jobs workers exhaust memory and fail the writer partway through the groups. groupby="leiden" gives one BigWig per cluster instead.

Pseudobulk BigWigs from ArchR, a barcoded BAM or a fragments file

Every route ends at one .bw per cell type:

  • 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 GRCh38 with chr-style names

The fragments are aligned to GRCh38, so the BigWigs use its chr1-style chromosome names, and the assembly has to spell them the same way. A BigWig whose names differ from the assembly's, such as Ensembl's 1 against chr1, draws empty unless the assembly has a name-alias table.

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

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 and an optional color and group. Swap each uri for the BigWig your pooling step wrote.

The fence lists three of the twelve cell types the figure draws; the build script writes the whole list, and the hosted demo config, https://jbrowse.org/code/jb2/main/test_data/config_demo.json, has it as pbmc5k_scatac_pseudobulk_hg38:

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

{
  "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://jbrowse.org/demos/scatac_pbmc5k/CD8_Naive.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "CD8 Memory",
        "group": "T cell",
        "color": "#3cb44b",
        "uri": "https://jbrowse.org/demos/scatac_pbmc5k/CD8_Memory.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "Naive B",
        "group": "B cell",
        "color": "#f58231",
        "uri": "https://jbrowse.org/demos/scatac_pbmc5k/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://jbrowse.org/demos/scatac_pbmc5k/CD8_Naive.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "CD8 Memory",
        "group": "T cell",
        "color": "#3cb44b",
        "uri": "https://jbrowse.org/demos/scatac_pbmc5k/CD8_Memory.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "Naive B",
        "group": "B cell",
        "color": "#f58231",
        "uri": "https://jbrowse.org/demos/scatac_pbmc5k/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://jbrowse.org/demos/scatac_pbmc5k/CD8_Naive.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "CD8 Memory",
        "group": "T cell",
        "color": "#3cb44b",
        "uri": "https://jbrowse.org/demos/scatac_pbmc5k/CD8_Memory.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "Naive B",
        "group": "B cell",
        "color": "#f58231",
        "uri": "https://jbrowse.org/demos/scatac_pbmc5k/Naive_B.bw"
      }
    ]
  }
}

The subadapters list has three things to set yourself:

  • the order: subadapters draw in the order given, so group them by lineage
  • color: copy the cluster color from your analysis, so a cell type matches its UMAP color
  • group: the sidebar tree branches on it, and Clustering rows clusters rows within each group

Without per-row names, colors or groups, the bigWigs shorthand takes a plain array of URLs and labels each row from its filename:

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

{
  "type": "MultiQuantitativeTrack",
  "trackId": "scatac_pseudobulk_simple",
  "name": "scATAC pseudobulk",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "bigWigs": [
      "https://jbrowse.org/demos/scatac_pbmc5k/CD8_Naive.bw",
      "https://jbrowse.org/demos/scatac_pbmc5k/CD8_Memory.bw",
      "https://jbrowse.org/demos/scatac_pbmc5k/Naive_B.bw"
    ]
  }
}
jbrowse add-track-json '{
  "type": "MultiQuantitativeTrack",
  "trackId": "scatac_pseudobulk_simple",
  "name": "scATAC pseudobulk",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "bigWigs": [
      "https://jbrowse.org/demos/scatac_pbmc5k/CD8_Naive.bw",
      "https://jbrowse.org/demos/scatac_pbmc5k/CD8_Memory.bw",
      "https://jbrowse.org/demos/scatac_pbmc5k/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://jbrowse.org/demos/scatac_pbmc5k/CD8_Naive.bw",
      "https://jbrowse.org/demos/scatac_pbmc5k/CD8_Memory.bw",
      "https://jbrowse.org/demos/scatac_pbmc5k/Naive_B.bw"
    ]
  }
}

CATlas serves hg38 per-cell-type coverage, one file per cell type, from https://decoder-genetics.wustl.edu/catlasv1/humanenhancer/data/bw/; its URLs go into the subadapters list like the PBMC files. Percent-encode the + in a cell-type name (T_lymphocyte_2_CD4%2B.bw); left unencoded, the URL breaks and the row loads with no data.

In grammar-of-graphics terms, a mark is the shape drawn, and the track menu switches it between modes: bar (the default, and the figures here) compares peak shape, and span maps score to color and fits more rows.

To check the rows against marker genes, paste two loci into the location box to open them side by side in one view, a T-cell marker (CD8A) and a B-cell marker (MS4A1):

chr2:86,780,000-86,820,000 chr11:60,450,000-60,490,000
Twelve per-cell-type BigWigs from the 10x 5k PBMC scATAC dataset, loaded as one MultiQuantitativeTrack, over CD8A and MS4A1 side by side in one view. The CD8, MAIT and NK rows have signal at CD8A, and only the two B rows have it at MS4A1.
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=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%22LinearWiggleDisplay%22%2C%22height%22%3A300%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+scatac%2Fpbmc5k_marker_swap

This session loads a large dataset and takes a while to open.

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 bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “NCBI RefSeq (GFF3Tabix)”.
    1. Track menu → Feature height → Compact
    2. Track menu → Show only genes (checked)
  3. Add the next track the same way.This figure uses “Single-cell ATAC pseudobulk by cell type (10x 5k PBMC, SnapATAC2)”.
  4. 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.

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

{
  "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": "LinearWiggleDisplay",
          "height": 300
        }
      ]
    }
  ]
}

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/code/jb2/main/test_data/config_demo.json")
app = JBrowseApp(
    assemblies=config["assemblies"],
    tracks=config.get("tracks", []),
    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": "LinearWiggleDisplay",
                    "height": 300
                }
            ]
        }
    ],
)
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/code/jb2/main/test_data/config_demo.json \
  --assembly hg38 --loc 'chr2:86,780,000-86,820,000 chr11:60,450,000-60,490,000' \
  --track ncbi_refseq_109_hg38_latest display:LinearBasicDisplay '{"displayMode":"compact","showOnlyGenes":true,"height":80}' \
  --track pbmc5k_scatac_pseudobulk_hg38 display:LinearWiggleDisplay '{"height":300}' \
  --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": [
    {
      "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": "LinearWiggleDisplay",
          "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 \
  --width 1500 --height 680 \
  --timeout 120000 \
  --spec 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 track from a folder of BigWigs

Two workflows write the list from a set of files:

  • Add multi-row track, in the Add track workflow, takes BigWig URLs one per line, or a JSON array of subadapter objects. On JBrowse Desktop it reads local .bw files directly, and the same form stacks BED or BigBed files, one row of features per file. Exporting the session gets the JSON config back out
  • jbrowse add-track --multiwig takes the whole set of BigWigs in place of the single positional file, labeling rows from the filenames, as in the command below
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; leave both off 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

build_scatac_pseudobulk.sh starts from SnapATAC2's annotated release of the 10x 5k PBMC dataset, whose standard pipeline and cell-type annotation tutorials leave every cell filtered, clustered and labeled, with its fragments beside its cell_type. The script:

  1. pools each cell type's fragments into one BigWig, in 25 bp bins and RPKM, as above
  2. orders the rows by lineage (T cells, NK, B cells, then myeloid) from a map written into the script, because the object's own label order scatters each lineage down the stack. Running it on your own experiment means replacing that map
  3. gives each row its cluster's color from the same object, so a cell type matches its UMAP
  4. loads the set as one track over hg38 and RefSeq genes
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

Rows that stay open everywhere usually mean the normalization step was skipped, since an unnormalized group's height tracks its cell count.

See also

Citations

Notes

  1. The BigWigs can also be viewed inline from the clustering environment through the Python anywidget interface or JBrowseR. ↩