Tutorial

The goals of this tutorial are to:

1. Show a standard workflow for analyzing single cell RNA-seq data.
2. Provide a brief overview of how SingleCellProjections.jl and ReproducibleJobs.jl work.
3. Show how to easily project one data set onto another.

We will use an Acute Myeloid Leukemia (AML) data set from the paper:

Henrik Lilljebjörn, Pablo Peña-Martínez, Hanna Thorsson, Rasmus Henningsson, Marianne Rissler, Niklas Landberg, Noelia Puente-Moncada, Sofia von Palffy, Vendela Rissler, Petr Stanek, Jonathan Desponds, Xiangfu Zhong, Gunnar Juliusson, Vladimir Lazarevic, Sören Lehmann, Magnus Fontes, Helena Ågerstam, Carl Sandén, Christina Orsmark-Pietras, Thoas Fioretos. "The AML cellular state space unveils NPM1 immune evasion subtypes with distinct clinical outcomes". Nat Commun 16, 10592 (2025).

You can download the data here:

• Single Cell RNA-seq Data - for this tutorial, we will use some of the .h5 files.
• Cell annotations

The data set contains 38 patient samples and 8 samples from healthy donors. The samples contain ~6000 cells on average, with measurements from 32738 genes.

First, we will load SingleCellProjections.jl and some other useful packages:

• ReproducibleJobs.jl provides a framework used by SingleCellProjections to handle computations, caching etc.
• DataFrames.jl and CSV.jl we use for handling and loading tabular data.
• SparseArrays.jl is used under the hood by SingleCellProjections. We just load it here so we can inspect some of the raw data.

using SingleCellProjections
import SingleCellProjections as SCP
using ReproducibleJobs
using CSV
using DataFrames
using SparseArrays

Loading Data

Let's start by loading two healthy Normal Bone Marrow (NBM) samples. These samples are from the same donor, but one has been enriched for immature (CD34 positive) cells. This was done because AML patients have more immature cells in their bone marrow than a healthy person. Combined, the NBM samples make a good reference when looking at AML samples later.

You'll need to change the paths to wherever you have downloaded your files.

nbm_names = ["NBM10-CD34", "NBM10-MNC"]
nbm_paths = joinpath.("samples", string.(nbm_names, ".h5"))

2-element Vector{String}:
 "samples/NBM10-CD34.h5"
 "samples/NBM10-MNC.h5"

Loading is done using SCP.load_counts.

raw_counts = SCP.load_counts(nbm_paths; sample_names=nbm_names)

load_counts (DataMatrixFunction)
├─ Vector
│  ├─ checksummedfilepath (prefetch) (cached)
│  │  └─ samples/NBM10-CD34.h5@2026-08-03T15:54:02.253
│  └─ checksummedfilepath (prefetch) (cached)
│     └─ samples/NBM10-MNC.h5@2026-08-03T15:54:05.654
├─ extra_id_cols: "feature_type"
├─ prefilter: "feature_type" => Base.Fix2(isequal, "Gene Expression")
└─ sample_names: ["NBM10-CD34", "NBM10-MNC"]

This creates a Job, which is a kind of specification or recipe for what to compute. Jobs are the cornerstone of ReproducibleJobs.jl since it makes it possible to reason about computations without performing them. Importantly, some results are cached, so when we come back another day, we do not need to recompute everything again. This mostly happens under the hood, and it will not be the focus of this tutorial.

To actually retrieve the data, we need to use fetch!.

c = fetch!(raw_counts)

DataMatrix (32738 variables and 13223 observations)
  Block Matrix (32×14)
  Variables: id, name, feature_type, genome, read, pattern, sequence
  Observations: cell_id, sample_name, barcode

The result is a DataMatrix (read more about them here: Data Matrices). Here we have a DataMatrix where the rows (variables) are genes and the columns (observations) are cells. Inside, there is a matrix with the raw counts, and annotation tables for the genes and the cells.

Here are the first few cells and genes:

julia> c.obs[1:6,:]6×3 DataFrame
 Row │ cell_id                        sample_name  barcode            
     │ String                         String       String             
─────┼────────────────────────────────────────────────────────────────
   1 │ NBM10-CD34_AAACCCAAGCGTATGG-1  NBM10-CD34   AAACCCAAGCGTATGG-1
   2 │ NBM10-CD34_AAACCCAAGCGTTACT-1  NBM10-CD34   AAACCCAAGCGTTACT-1
   3 │ NBM10-CD34_AAACCCAGTGACTAAA-1  NBM10-CD34   AAACCCAGTGACTAAA-1
   4 │ NBM10-CD34_AAACCCAGTTTCGACA-1  NBM10-CD34   AAACCCAGTTTCGACA-1
   5 │ NBM10-CD34_AAACGAAAGTCATGCT-1  NBM10-CD34   AAACGAAAGTCATGCT-1
   6 │ NBM10-CD34_AAACGAACACTACCCT-1  NBM10-CD34   AAACGAACACTACCCT-1
julia> c.var[1:6,:]6×7 DataFrame
 Row │ id               name          feature_type     genome  read    pattern  sequence 
     │ String           String        String           String  String  String   String   
─────┼───────────────────────────────────────────────────────────────────────────────────
   1 │ ENSG00000243485  MIR1302-10    Gene Expression  hg19
   2 │ ENSG00000237613  FAM138A       Gene Expression  hg19
   3 │ ENSG00000186092  OR4F5         Gene Expression  hg19
   4 │ ENSG00000238009  RP11-34P13.7  Gene Expression  hg19
   5 │ ENSG00000239945  RP11-34P13.8  Gene Expression  hg19
   6 │ ENSG00000237683  AL627309.1    Gene Expression  hg19

You rarely need to access the raw count matrix directly. It is stored in a blocked format that's efficient for the computations that SingleCellProjections needs. But to take a look, we can convert it to a regular sparse matrix and show a small part of it:

convert(SparseMatrixCSC, c.matrix)[1:100,1:100]

100×100 SparseMatrixCSC{Int64, Int32} with 766 stored entries:
⎡⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠈⠀⠀⠈⠀⠀⠀⠈⠀⠈⠈⠈⠀⠀⠀⠈⠀⠁⠀⠈⠀⠀⠁⠀⎥
⎢⠠⠀⠀⠆⢀⠀⠀⠀⠂⢀⠠⠤⢂⠄⠐⠀⠀⢠⠀⠠⠀⠠⠠⠀⎥
⎢⠈⠀⠀⠀⠨⠁⠂⠀⠀⠀⠀⠁⠑⠀⠀⠀⠁⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠒⠐⠒⠀⠒⠂⠊⠒⠒⠐⠂⠐⠂⠒⠐⠂⠐⠐⠀⠐⠐⠀⠐⎥
⎢⠀⠉⠀⠈⠁⠁⠀⠀⠈⠈⠅⠀⠀⠈⠍⠀⠄⠀⠈⠈⠈⠀⠀⠀⎥
⎢⣔⣦⣶⣖⢲⡦⣚⠤⡂⣴⠖⡶⣆⡔⣴⢒⡐⣔⢀⠄⠐⡶⢆⣆⎥
⎢⠐⡒⠱⠙⡛⠚⡐⠒⠶⠷⢙⡓⢒⡖⠒⠖⠐⠜⠑⠵⠂⡲⠐⠂⎥
⎢⠋⠛⠟⠛⠙⠛⠋⠙⠓⠟⠻⠛⠛⠛⠛⠚⠊⠻⠛⠟⠛⠙⠙⠛⎥
⎢⠤⡦⠤⢤⠤⡦⠆⢥⡦⢶⡶⡦⠴⢶⢄⡠⠠⢔⢤⢴⢤⢴⠢⠀⎥
⎢⡚⡗⣌⣵⢍⣅⡡⠉⣨⢕⢝⣚⡒⠆⠌⣣⠈⠲⢱⣀⠡⡏⠘⡈⎥
⎣⠐⠛⠒⠛⠚⠀⠂⠐⠱⢓⠛⠊⠒⠁⠂⠘⠐⠛⠒⠚⠘⠓⠒⠒⎦

Quality Filtering

Some of the cells in the data set are of poor quality, and we want to get rid of them before continuing the analysis. Commonly used quality measures are the fraction of reads from Mitochondrial genes, the total read count per cell, and the number of genes with non-zero counts. A high fraction of Mitochondrial reads can indicate that a cell is dying. The total read count and number of expressed genes help determine whether there is enough data for a meaningful analysis.

We set up the quality annotations like this:

counts = SCP.var_counts_fraction(raw_counts, "fraction_mt", "name"=>startswith("MT-"))
counts = SCP.var_counts_sum(counts, "total_RNA_count")
counts = SCP.var_counts_sum(!iszero, counts, "nonzero_RNA_count")
counts = SCP.obs_counts_sum(!iszero, counts, "nonzero_cell_count")

obs_counts_sum (DataMatrixFunction)
├─ var_counts_sum (DataMatrixFunction)
│  ├─ var_counts_sum (DataMatrixFunction)
│  │  ├─ var_counts_fraction (DataMatrixFunction)
│  │  │  ├─ load_counts (DataMatrixFunction)
│  │  │  │  ├─ Vector
│  │  │  │  │  ⋮
│  │  │  │  │  
│  │  │  │  ├─ extra_id_cols: "feature_type"
│  │  │  │  ├─ prefilter: "feature_type" => Base.Fix2(isequal, "Gene Expression")
│  │  │  │  └─ sample_names: ["NBM10-CD34", "NBM10-MNC"]
│  │  │  ├─ "fraction_mt"
│  │  │  ├─ "name" => Base.Fix2(startswith, "MT-")
│  │  │  ├─ Returns(true)
│  │  │  └─ project_ids: :intersect
│  │  ├─ "total_RNA_count"
│  │  ├─ Returns(true)
│  │  ├─ f: identity
│  │  └─ project_ids: :intersect
│  ├─ "nonzero_RNA_count"
│  ├─ Returns(true)
│  ├─ f: !iszero
│  └─ project_ids: :intersect
├─ "nonzero_cell_count"
├─ Returns(true)
├─ f: !iszero
└─ project_ids: :no

And apply filtering:

filtered = SCP.filter_obs("fraction_mt" => <(0.15), counts)
filtered = SCP.filter_obs("total_RNA_count" => >(1000), filtered)
filtered = SCP.filter_obs("nonzero_RNA_count" => >(500), filtered)

filter_matrix (Preprocess)
├─ filter_matrix (Preprocess)
│  ├─ filter_matrix (Preprocess)
│  │  ├─ obs_counts_sum (DataMatrixFunction)
│  │  │  ├─ var_counts_sum (DataMatrixFunction)
│  │  │  │  ├─ var_counts_sum (DataMatrixFunction)
│  │  │  │  │  ⋮
│  │  │  │  │  
│  │  │  │  ├─ "nonzero_RNA_count"
│  │  │  │  ├─ Returns(true)
│  │  │  │  ├─ f: !iszero
│  │  │  │  └─ project_ids: :intersect
│  │  │  ├─ "nonzero_cell_count"
│  │  │  ├─ Returns(true)
│  │  │  ├─ f: !iszero
│  │  │  └─ project_ids: :no
│  │  └─ fobs: "fraction_mt" => Base.Fix2(<, 0.15)
│  └─ fobs: "total_RNA_count" => Base.Fix2(>, 1000)
└─ fobs: "nonzero_RNA_count" => Base.Fix2(>, 500)

And if we fetch! the result:

fetch!(filtered)

DataMatrix (32738 variables and 11969 observations)
  Block Matrix (32×14)
  Variables: id, name, feature_type, genome, read, pattern, sequence, nonzero_cell_count
  Observations: cell_id, sample_name, barcode, fraction_mt, total_RNA_count, nonzero_RNA_count

We see that filtering reduced the number of cells from 13223 to 11969, and that the new annotations are present in the data matrix.

Transformation

The raw counts data is not suitable for analyses like PCA, since the data is far from normally distributed. A common strategy to handle this is to transform the data. Here we will use SCTransform (see also original sctransform implementation in R).

transformed = SCP.sctransform(filtered)
fetch!(transformed)

DataMatrix (19371 variables and 11969 observations)
  A+B₁B₂B₃
  Variables: id, name, feature_type, genome, read, pattern, sequence, nonzero_cell_count
  Observations: cell_id, sample_name, barcode, fraction_mt, total_RNA_count, nonzero_RNA_count

From the output, we see that the number of variables has been reduced, since by default, sctransform removes variables present in very few cells.

The matrix is now shown as A+B₁B₂B₃. This is normally not very important from the user's point of view, but it is critical for explaining how SingleCellProjections can be fast and not use too much memory. Instead of storing the SCTransformed matrix as a huge dense matrix, it is stored in memory as a MatrixExpression, in this case a sparse matrix A plus a product of three smaller matrices B₁,B₂ and B₃.

Normalization

After transformation we always want to normalize the data. At the very least, data should be centered for PCA to work properly. This can be achieved by just running SCP.normalize_matrix with the default parameters. Here, we also want to regress out "fraction_mt". You can add more obs annotations (categorical and/or numerical) to regress out if needed.

normalized = SCP.normalize_matrix(transformed, "fraction_mt")
fetch!(normalized)

DataMatrix (19371 variables and 11969 observations)
  A+B₁B₂B₃+(-β)X'
  Variables: id, name, feature_type, genome, read, pattern, sequence, nonzero_cell_count
  Observations: cell_id, sample_name, barcode, fraction_mt, total_RNA_count, nonzero_RNA_count

Now the matrix is shown as A+B₁B₂B₃+(-β)X', i.e. another low-rank term was added to handle the normalization/regression. Since all ReproducibleJobs.jl results are read-only, the first two terms can be reused, ensuring that memory is not wasted.

Principal Component Analysis

Principal Component Analysis (PCA) is commonly used for single cell expression data for two major reasons:

1. It accurately finds a more compact representation of the data. This is important computationally, since it greatly reduces computation time for downstream analyses.
2. It reduces noise by keeping only the common variations in the data set.

reduced = SCP.pca(normalized; nsv=40)
fetch!(reduced)

DataMatrix (40 variables and 11969 observations)
  ReadOnlyArrays.ReadOnlyMatrix{Float64, Matrix{Float64}}
  Variables: PC_id
  Observations: cell_id, sample_name, barcode, fraction_mt, total_RNA_count, nonzero_RNA_count

Now, the data set has been reduced to 40 variables (the top 40 principal components, ranked by variance explained), and is represented as a single dense matrix. The cell (observation) annotations are left untouched.

Annotations

In preparation for the visualization below (or other analyses), we also load some cell-level annotations for our data set.

annots_path = "annotations/scRNA_AML.tsv"
annots = SCP.load_csv(annots_path)
reduced = SCP.annotate_obs(reduced, annots)

Visualization

Analyses should generally be performed on the normalized or PCA-reduced data, not on visualization embeddings. However, for visualization purposes, we want to reduce to 2 or 3 dimensions.

using Colors
using WGLMakie

function scatter_3d(job)
    matrix = fetch!(SCP.get_matrix(job))
    fig = Figure(; size=(768, 768))
    ax = LScene(fig[1, 1])
    scatter!(ax, matrix; color = :black, markersize = 4)
    fig
end

function scatter_categorical_3d(job, annot_name; bg=nothing, colors=nothing)
    matrix = fetch!(SCP.get_matrix(job))
    annot = fetch!(SCP.value_column_data(SCP.annotation(SCP.get_obs(job), annot_name)))

    fig = Figure(; size=(768, 768))
    ax = LScene(fig[1, 1])

    if bg !== nothing
        bg_matrix = fetch!(SCP.get_matrix(bg))
        scatter!(ax, bg_matrix; color=colorant"#BFCCE6", markersize=2)
    end

    if colors !== nothing
        unique_annotations = unique(annot)
        unique_annotations_set = Set(unique_annotations)
        colors = filter(x->x[1] in unique_annotations_set, colors)

        categories = first.(colors)
        @assert isempty(setdiff(unique_annotations, categories)) # ensure all categories have colors specified

        plots = [scatter!(ax, matrix[:,isequal.(annot, cat)]; markersize=5, color) for (cat,color) in colors]
    else
        categories = unique(annot)
        plots = [scatter!(ax, matrix[:,isequal.(annot, cat)]; markersize=5) for cat in categories]
    end

    axislegend(ax, plots .=> Ref((;markersize=16)), categories) # use a larger markersize in the legend
    fig
end

scatter_categorical_3d (generic function with 1 method)

scatter_3d(data)
scatter_categorical_3d(data, "celltype.aml")

Force Layout

To embed the points in 2 or 3 dimensions using a Force Layout (also known as a SPRING plot), we set it up like this:

fl = SCP.force_layout(reduced; ndim = 3,
                                seed = 4567,
                                k = 100,
                                k_projection = 25)
scatter_3d(fl)

To make a nicer visualization, we use a celltype annotation from the downloaded data to color the plot, and apply a utility function to rotate it such that the most immature cells (Hematopoietic Stem Cells, or HSCs for short) move to the top.

transform = SCP.find_optimal_coord_transform(fl, "celltype.aml"=>isequal("HSC"), "celltype.aml"=>isequal("T-cells"), "celltype.aml"=>isequal("B-cells"))
fl = SCP.transform_coords(fl, transform; keep_var=true)
colors = ["AML Immature" => colorant"#fec44f", "HSC" => colorant"#66b266", "Monocytes" => colorant"#fa9fb5", "GMP" => colorant"#008000", "Megakaryocytic cells" => colorant"#e5687e", "LMPP" => colorant"#756bb1", "Erythroid cells" => colorant"#a7241d", "B-cells" => colorant"#a64ca6", "T-cells" => colorant"#7fb4e5", "NK-cells" => colorant"#196fbe", "Dendritic cells" => colorant"#ff4d00"]
scatter_categorical_3d(fl, "celltype.aml"; colors)

UMAP

using UMAP
umapped = SCP.umap(reduced; ndim=3)
scatter_categorical_3d(umapped, "celltype.aml"; colors)

t-SNE

t-SNE is also supported, just run:

using TSne
tsne_job = SCP.tsne(reduced; ndim=3)

Projections

SingleCellProjections is built to make it very easy to project one dataset onto another. This is useful for comparing new samples against an established reference.

First, we just choose one or more samples to project:

# Name and location of AML file
proj_name = "AML28"
proj_path = joinpath("samples", string(proj_name, ".h5"))
proj_raw_counts = SCP.load_counts(proj_path; sample_names=proj_name)

load_counts (DataMatrixFunction)
├─ Vector
│  └─ checksummedfilepath (prefetch) (cached)
│     └─ samples/AML28.h5@2026-08-03T15:53:06.252
├─ extra_id_cols: "feature_type"
├─ prefilter: "feature_type" => Base.Fix2(isequal, "Gene Expression")
└─ sample_names: ["AML28"]

Then, a single call to SCP.project sets up the entire analysis pipeline and projects the AML sample onto the reference map created from the NBM samples.

proj_fl = SCP.project(fl, raw_counts=>proj_raw_counts)
scatter_categorical_3d(proj_fl, "celltype.aml"; bg=fl, colors)

And here is the projection onto the UMAP embedding of the NBM samples:

proj_umapped = SCP.project(umapped, raw_counts=>proj_raw_counts)
scatter_categorical_3d(proj_umapped, "celltype.aml"; bg=umapped, colors)