All in One View

Content from Introduction to Single-Cell RNA-Seq

Last updated on 2026-10-27 | Edit this page

Estimated time: 30 minutes

Overview

Questions

What is single-cell RNA sequencing and how does it differ from bulk RNA-seq?
How does the 10x Genomics Chromium platform capture and barcode individual cells?
What does the resulting count matrix look like and why are UMIs important?
What are the main computational steps in a scRNA-seq analysis workflow?
What will this workshop cover and what is out of scope?

Objectives

Explain what single-cell RNA-seq measures and why cellular heterogeneity matters
Describe the 10x Chromium droplet-based workflow from cell capture through sequencing
Interpret the structure of a UMI count matrix and explain why it is sparse
Outline the end-to-end computational pipeline covered in this workshop
Identify the scope and limitations of this workshop

What Is Single-Cell RNA-Seq?

Every tissue in your body is made up of many different cell types, and each cell type carries out its function by expressing a distinct set of genes. Traditional bulk RNA-seq measures gene expression by grinding up a piece of tissue and sequencing the RNA from all cells at once. The result is a single average expression profile that blends signals from every cell type together. This is perfectly fine when the goal is to compare overall gene expression between two conditions, but it hides all of the cell-to-cell variation within the sample.

Think of it like making a smoothie. If you blend strawberries, blueberries, and bananas together, you can tell that the smoothie contains fruit, and you might even detect the flavors of the individual ingredients. But you cannot tell how many strawberries went in, whether one banana was riper than the others, or whether there was a single grape hiding at the bottom of the blender. Bulk RNA-seq gives you the smoothie. Single-cell RNA-seq (scRNA-seq) gives you the fruit bowl: every piece of fruit is kept separate, so you can count, sort, and inspect them one by one.

Animation: Bulk averages, single-cell counts

scRNA-seq works by isolating individual cells, capturing their mRNA, and sequencing each cell’s transcriptome independently. Instead of one expression profile per sample, you get thousands of individual expression profiles, one per cell. This lets you:

Discover cell types that were previously unknown or hard to distinguish
Track how cells change during development, disease progression, or in response to treatment
Identify rare populations such as stem cells or drug-resistant tumor cells that are invisible in bulk data
Map cell-cell communication by looking at which ligands and receptors are expressed in neighboring cell types

Since the first scRNA-seq experiment was published in 2009, the technology has scaled from single cells to millions of cells per experiment. It is now a routine tool in immunology, neuroscience, cancer biology, and developmental biology.

Callout

Bulk vs. Single-Cell

Feature	Bulk RNA-seq	Single-cell RNA-seq
Resolution	One profile per sample	One profile per cell
Input	~1 million cells (pooled)	~1,000 – 20,000+ individual cells
Output	One expression vector per sample	A cells-by-genes matrix
Genes detected per sample	~15,000 – 20,000	~1,000 – 5,000 per cell (due to dropout)
Cost per sample	~$200 – $500	~$2,000 – $5,000 (library + sequencing)
Best for	Comparing conditions, differential expression with biological replicates	Discovering cell types, studying heterogeneity, trajectory analysis

Single-cell experiments detect fewer genes per cell than bulk experiments because each cell contains only a small amount of mRNA. However, when you combine information across thousands of cells, the total number of genes detected across the full dataset is comparable to bulk RNA-seq.

How Does 10x Chromium Work?

The 10x Genomics Chromium platform is the most widely used droplet-based scRNA-seq system. It can profile thousands to tens of thousands of cells in a single run. The core idea is simple: wrap each cell in its own tiny droplet along with a barcoded gel bead, so that every mRNA molecule from that cell gets tagged with the same unique barcode.

Diagram showing the 10x Chromium microfluidic chip where cells, gel beads, and oil are combined to form gel bead-in-emulsion (GEM) droplets. Each GEM contains one gel bead and ideally one cell.

Here is how the process works, step by step:

1. GEM generation. A suspension of single cells is loaded onto a Chromium microfluidic chip along with gel beads and partitioning oil. The chip combines these three inputs into tiny droplets called GEMs (Gel Beads-in-Emulsion). Each GEM is designed to contain exactly one gel bead and at most one cell. In practice, most GEMs are empty (no cell), a fraction contain one cell, and a small number contain two or more cells (these are called doublets and will be addressed during quality control).

2. Cell lysis and barcoding. Inside each GEM, the cell is lysed and its mRNA is released. The gel bead dissolves, releasing millions of oligonucleotide primers that all share the same cell barcode – a 16-nucleotide sequence unique to that bead. Each primer also carries a UMI (Unique Molecular Identifier), a random 12-nucleotide sequence that is different on every primer. The mRNA molecules hybridize to these primers via their poly(A) tails and are reverse-transcribed into cDNA. Because all primers in the same GEM have the same cell barcode, every cDNA molecule from the same cell is tagged with the same barcode. And because each primer has a different UMI, each original mRNA molecule gets a unique tag.

3. Library construction. After reverse transcription, the GEMs are broken and the barcoded cDNA from all cells is pooled together. The cDNA is amplified by PCR, fragmented, and prepared into a sequencing library. Even though the cDNA from all cells is now mixed in one tube, the cell barcodes embedded in each molecule let us computationally sort the reads back to their cell of origin.

Animation: Inside a 10x Chromium GEM

4. Sequencing. The library is sequenced on an Illumina sequencer, producing paired-end reads plus an index read. The read structure is:

Schematic showing three sequencing reads from a 10x Chromium library. Read 1 is 28 bp and contains the 16 bp cell barcode followed by the 12 bp UMI. Read 2 is variable length and contains the cDNA insert that maps to the transcriptome. The I1 index read contains the sample index for demultiplexing.

Read 1 (28 bp): The first 16 bases are the cell barcode (identifies which cell). The next 12 bases are the UMI (identifies which original mRNA molecule). Read 1 does not contain any transcript sequence.
Read 2 (variable length, typically 90–150 bp): This is the actual cDNA insert that gets mapped to the reference transcriptome to determine which gene the mRNA came from.
Index read (I1, 8 bp): The sample index used for demultiplexing when multiple samples are pooled on the same sequencing lane.

The alignment software (Cell Ranger or STARsolo) reads all three components: it uses the index read to assign reads to samples, the cell barcode to assign reads to cells, the UMI to count unique molecules, and Read 2 to identify the gene of origin.

Animation: What the aligner does with each read

Callout

Other Single-Cell Platforms

While this workshop focuses on 10x Genomics Chromium, several other scRNA-seq platforms exist:

Drop-seq: An early droplet-based method developed in the McCarroll lab (Harvard). Uses a similar barcoding approach but with lower per-cell capture efficiency than 10x Chromium.
inDrop: Another droplet-based platform developed in the Bhatt lab (Harvard). Uses hydrogel beads that release barcoded primers upon UV exposure.
Smart-seq2 / Smart-seq3: Plate-based methods that capture full-length transcripts. They provide much deeper coverage per cell but are limited to hundreds of cells per experiment and are significantly more expensive per cell.
Parse Biosciences (Evercode): A combinatorial barcoding method that does not require specialized microfluidic equipment. Cells are split across wells and barcoded in multiple rounds, enabling very high throughput.

Each platform has trade-offs in throughput, cost, sensitivity, and transcript coverage. The 10x Chromium platform is the most common choice for large-scale cell atlas projects and is the standard in most core facilities, which is why we use it here.

The Count Matrix

The end result of raw data processing (covered in the next episode) is a UMI count matrix. This is a table where each row is a gene, each column is a cell, and each entry is the number of unique mRNA molecules (UMIs) detected for that gene in that cell.

Schematic of a count matrix with genes on rows and cells on columns. Most cells in the matrix contain zero, with occasional non-zero integer counts scattered throughout, illustrating the sparsity typical of single-cell RNA-seq data.

A few important properties of this matrix:

It is very sparse. A typical human cell expresses around 2,000 to 5,000 genes out of a total of approximately 30,000 protein-coding genes in the genome. On top of that, scRNA-seq has limited sensitivity – it captures only a fraction of the mRNA molecules actually present in the cell. As a result, roughly 95% of the entries in the count matrix are zeros. This is not an error; it is an inherent property of the technology. Specialized statistical methods have been developed to handle this sparsity.

Counts are integers. Because UMIs deduplicate PCR artifacts, each entry in the matrix represents the number of distinct original mRNA molecules detected, not the number of sequencing reads. This makes UMI counts more quantitative than read counts.

UMIs solve the PCR duplication problem. During library preparation, cDNA is amplified by PCR so there is enough material to sequence. Without UMIs, a single mRNA molecule that was amplified 100 times would produce 100 reads and appear to be highly expressed. With UMIs, all 100 reads share the same UMI sequence and are collapsed into a single count. This means UMI counts are proportional to the actual number of mRNA molecules in the cell, not to PCR amplification efficiency.

Animation: Why UMI counts are integers

For the PBMC 10k v3 dataset we will use in this workshop, the count matrix contains 11,809 cell barcodes and 38,606 genes. After removing genes detected in fewer than 3 cells and low-quality cells, we will work with 11,310 cells and 29,155 genes.

Analysis Workflow Overview

The computational analysis of scRNA-seq data follows a series of well-defined steps. Each step in the pipeline builds on the output of the previous one. This workshop covers the complete workflow from raw sequencing data through biological interpretation.

Flowchart of the scRNA-seq analysis pipeline covered in this workshop. Eight boxes arranged vertically show the progression: Raw Data Processing, Quality Control, Normalization and Feature Selection, Dimensionality Reduction and Clustering, Cell Type Annotation, Multi-Sample Integration, and Differential Expression. Arrows connect each step to the next.

1. Raw Data Processing (Episode 2). Starting from FASTQ files, we align reads to the human reference genome and quantify gene expression per cell. We will use two tools: Cell Ranger (the commercial pipeline from 10x Genomics) and STARsolo (a fast open-source alternative). Both produce the UMI count matrix described above.

2. Quality Control (Episode 3). Not all barcodes in the count matrix represent healthy, real cells. Some correspond to empty droplets, dying cells, or doublets. We filter cells based on three key metrics: the number of detected genes, the total UMI count, and the percentage of reads mapping to mitochondrial genes (a marker of cell stress). This step ensures we only carry high-quality cells forward.

3. Normalization and Feature Selection (Episode 4). Cells are sequenced to different depths, so raw counts are not directly comparable between cells. Normalization corrects for these differences. We will apply two methods: LogNormalize (a simple and widely used approach) and SCTransform (a more sophisticated variance-stabilizing method). We also select highly variable genes – the 2,000 or so genes that vary the most across cells and carry the strongest biological signal.

4. Dimensionality Reduction and Clustering (Episode 5). With 20,000+ genes per cell, we need to reduce the data to a manageable number of dimensions. PCA (Principal Component Analysis) compresses the data into the top principal components. We then use UMAP to create a 2D visualization and apply the Louvain algorithm (Seurat’s default) to group similar cells into clusters.

5. Cell Type Annotation (Episode 6). Clusters are just numbers until we assign biological meaning. We identify cell types by examining marker genes – genes that are specifically expressed in one cluster but not others. We will annotate clusters both manually (using known PBMC markers) and automatically (using SingleR, which compares expression profiles to labeled reference datasets).

6. Multi-Sample Integration (Episode 7). When analyzing cells from multiple samples or experimental conditions, batch effects can cause cells to cluster by sample rather than by cell type. Integration algorithms correct for these technical differences while preserving real biological variation. We will use Seurat’s CCA (Canonical Correlation Analysis) integration on an IFN-beta stimulated PBMC dataset to learn how to combine control and treated samples.

7. Differential Expression (Episode 8). Finally, we ask: which genes change between conditions? We will compare gene expression between control and stimulated cells within each cell type. We will cover both cell-level testing and pseudobulk analysis (the statistically preferred approach for multi-sample experiments). We will also run gene ontology enrichment to connect lists of differentially expressed genes to biological pathways.

Key Points

Single-cell RNA-seq measures gene expression in individual cells, revealing cellular heterogeneity that bulk RNA-seq averages out
The 10x Chromium platform uses droplet-based barcoding: each cell gets a unique 16 bp cell barcode and each mRNA molecule gets a 12 bp UMI
The UMI count matrix is very sparse (~95% zeros) because each cell expresses only a fraction of all genes and capture efficiency is limited
The scRNA-seq analysis pipeline progresses through alignment, QC, normalization, dimensionality reduction, clustering, annotation, integration, and differential expression
This workshop focuses on 10x Chromium 3’ GEX data analyzed with Seurat v5 in R on the Purdue Negishi HPC cluster

Content from Raw Data Processing

Last updated on 2026-10-27 | Edit this page

Estimated time: 60 minutes

Overview

Questions

How do we organize scRNA-seq data files on an HPC cluster?
How do we align 10x Chromium FASTQ files with Cell Ranger?
How do we run STARsolo as an open-source alternative?
How do the outputs from Cell Ranger and STARsolo compare?

Objectives

Set up a working directory and verify FASTQ and reference files on the cluster
Write a SLURM job script for Cell Ranger count and check it with sbatch --test-only
Write a SLURM job script for STARsolo and check it the same way
Interpret Cell Ranger’s web_summary.html and key QC metrics
Compare the filtered count matrices produced by both tools

Prerequisite

Prerequisites

This episode requires an active account on Purdue’s Negishi cluster and familiarity with basic Linux commands (cd, ls, mkdir, cat). You should be comfortable submitting SLURM jobs with sbatch and monitoring them with squeue. This episode requires shell access to the Negishi cluster via SSH. See the SSH Setup section for connection instructions.

Callout

Cell Ranger Licensing

Cell Ranger is proprietary software from 10x Genomics. It is free to download and use, but it requires acceptance of the 10x Genomics End User License Agreement. STARsolo is fully open source (MIT license) and produces equivalent results. If you cannot use Cell Ranger due to licensing constraints, STARsolo is an excellent alternative. This workshop covers both tools so you can choose whichever is appropriate for your situation.

Data Organization

Before running any alignment, we need a clean directory structure. All workshop files live under ${RCAC_SCRATCH}/scrna_workshop/. Let’s create the directory layout and copy the data.

BASH

mkdir -p ${RCAC_SCRATCH}/scrna_workshop/{fastq,reference,cellranger_output,starsolo_output}

This creates the following structure:

${RCAC_SCRATCH}/scrna_workshop/
├── fastq/                  # Raw FASTQ files
├── reference/              # Reference genome and index
├── cellranger_output/      # Cell Ranger results
└── starsolo_output/        # STARsolo results

Now copy the pre-staged workshop data from the shared depot:

BASH

rsync -avP /depot/workshop/data/scrna_workshop/ ${RCAC_SCRATCH}/scrna_workshop/

Extract the 10x barcode whitelist from the Cell Ranger container (needed for STARsolo):

BASH

module load biocontainers
singularity exec ${BIOC_IMAGE_DIR}/cumulusprod_cellranger:10.0.0.sif \
    zcat /software/cellranger-10.0.0/lib/python/cellranger/barcodes/3M-february-2018_TRU.txt.gz \
    > ${RCAC_SCRATCH}/scrna_workshop/reference/3M-february-2018.txt

Verify that the FASTQ files are present:

BASH

ls -1 ${RCAC_SCRATCH}/scrna_workshop/fastq/

Expected output

OUTPUT

pbmc_10k_v3_S1_L001_I1_001.fastq.gz
pbmc_10k_v3_S1_L001_R1_001.fastq.gz
pbmc_10k_v3_S1_L001_R2_001.fastq.gz
pbmc_10k_v3_S1_L002_I1_001.fastq.gz
pbmc_10k_v3_S1_L002_R1_001.fastq.gz
pbmc_10k_v3_S1_L002_R2_001.fastq.gz

Verify that the reference genome is present:

BASH

ls ${RCAC_SCRATCH}/scrna_workshop/reference/

Expected output

OUTPUT

3M-february-2018.txt  refdata-gex-GRCh38-2024-A

Understanding the FASTQ file names

The file names follow the 10x Genomics naming convention:

pbmc_10k_v3_S1_L001_R1_001.fastq.gz
│            │  │    │   └── File number (always 001 for single-file lanes)
│            │  │    └────── Read type: R1 or R2
│            │  └─────────── Lane number: L001 or L002
│            └────────────── Sample index: S1
└─────────────────────────── Sample name

This dataset was sequenced across two lanes (L001, L002), producing six files: two lanes times three reads per lane. The I1 files hold the sample index read, which cellranger count and STARsolo do not need for this dataset.

What is in each read?

Read 1 (R1) contains the cell barcode and UMI. It is 28 bp long: the first 16 bases are the cell barcode (which cell this read came from) and the next 12 bases are the UMI (which original mRNA molecule this read represents). Read 1 does not contain any transcript sequence.

Read 2 (R2) contains the cDNA insert. This is the actual transcript fragment that gets mapped to the reference genome to determine which gene the mRNA came from. It is typically 90–150 bp long.

Let’s inspect the first few reads from each file to see this in practice:

BASH

cd ${RCAC_SCRATCH}/scrna_workshop/fastq
zcat pbmc_10k_v3_S1_L001_R1_001.fastq.gz | head -8

Expected output

You will see two FASTQ records. Each record has four lines: a header line starting with @, the sequence, a + separator, and quality scores. Notice that every R1 sequence is exactly 28 bp – that is the 16 bp cell barcode followed by the 12 bp UMI.

OUTPUT

@A00228:279:HFWFVDMXX:1:1101:3260:1000 1:N:0:NCAAGATG
NACCAACAGTCGAATAGTGTCATCTGCT
+
#FFFFFFFFFFFFFFFFFFFFFFFFFFF
@A00228:279:HFWFVDMXX:1:1101:3821:1000 1:N:0:NCAAGATG
NTGTCTTAGCCTAGGACGGCCTCCGCCA
+
#FFFFFFFFFFFFFFFFFFFFFFFFFFF

BASH

zcat pbmc_10k_v3_S1_L001_R2_001.fastq.gz | head -8

Expected output

R2 reads are longer (91 bp in this dataset) and contain transcript sequence that will be mapped to the genome.

OUTPUT

@A00228:279:HFWFVDMXX:1:1101:3260:1000 2:N:0:NCAAGATG
NTATAAAATCACCACGGTCTTTAGCCATGCACAAACGGTAGTTTTGTGTGTTGGCTGCTCCACTGTCCTCTGCCAGCCTACAGGAGGAAAA
+
#FFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@A00228:279:HFWFVDMXX:1:1101:3821:1000 2:N:0:NCAAGATG
NTTCTATTGGAAACCCGGTCTTTACAAAAAAATACAAAAATCAGCTGGGCGTTGGCCGCGCGTGGTGGCTCACACCTGTAATCTCAGCACT
+
#FFFFFFFFFFFFFFFFFFFFFF:,FFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFF:FFFFFF:FFF:F:FFFFFFFFFFFFFFF

The reference genome

The reference genome is a 10x Genomics pre-built package based on the human GRCh38 assembly. This package contains everything Cell Ranger needs to align reads:

BASH

ls ${RCAC_SCRATCH}/scrna_workshop/reference/refdata-gex-GRCh38-2024-A/

Expected output

OUTPUT

fasta/     genes/     reference.json     star/

The key components are:

fasta/genome.fa – the genome sequence (FASTA)
genes/genes.gtf – the gene annotation (GTF)
star/ – a pre-built STAR genome index (used internally by Cell Ranger)
reference.json – metadata about the reference build

Both Cell Ranger and STARsolo use STAR for alignment under the hood. Cell Ranger ships its own bundled STAR binary and uses the index inside this reference package. For STARsolo, we can either build our own STAR index or reuse the one inside this package.

Processing with Cell Ranger

Cell Ranger is the official analysis pipeline from 10x Genomics. The cellranger count command takes FASTQ files and a reference genome as input and produces a filtered UMI count matrix as output. We will run it as a SLURM batch job on Negishi.

Create the job script:

BASH

cat << 'EOF' > ${RCAC_SCRATCH}/scrna_workshop/run_cellranger.sh
#!/bin/bash
#SBATCH --nodes=1
#SBATCH --ntasks=64
#SBATCH --time=1-00:00:00
#SBATCH --job-name=cellranger
#SBATCH --account=rcac-rnaseq
#SBATCH --partition=cpu
#SBATCH --output=cellranger_%j.out
#SBATCH --error=cellranger_%j.err

# Load modules
ml --force purge
ml biocontainers
ml cellranger/10.0.0

# Move to the output directory
cd ${RCAC_SCRATCH}/scrna_workshop/cellranger_output

# Run Cell Ranger count
cellranger count \
    --id=pbmc10k \
    --transcriptome=${RCAC_SCRATCH}/scrna_workshop/reference/refdata-gex-GRCh38-2024-A \
    --fastqs=${RCAC_SCRATCH}/scrna_workshop/fastq \
    --create-bam true \
    --sample=pbmc_10k_v3 \
    --localcores=${SLURM_CPUS_ON_NODE} \
    --localmem=100
EOF

Check the script without running it:

BASH

cd ${RCAC_SCRATCH}/scrna_workshop
sbatch --test-only run_cellranger.sh

Expected output

OUTPUT

sbatch: Job 44575176 to start at 2026-10-22T18:54:55 a using 64 processors on nodes a324 in partition cpu

Your job number, start time and node will differ.

--test-only makes SLURM check the script’s #SBATCH lines and estimate when the job could start. It submits nothing; the job number it prints is not a real job. A typo in the account, the partition or an option shows up here as an error, before you have waited in a queue.

We do not run the job during the workshop. Cell Ranger needs 64 cores for about 5.5 hours (see the runtime comparison at the end of this episode), and the workshop account runs only a few jobs that size at once. The instructor ran this same script before the workshop and will show its outputs, and the next episode starts from a pre-computed copy of the same matrix.

Callout

After the workshop: run it for real

At the end of the day, first close your RStudio session so its cores are free (Open OnDemand, My Interactive Sessions, Delete). Then submit the job and check that it is queued:

BASH

cd ${RCAC_SCRATCH}/scrna_workshop
sbatch run_cellranger.sh
squeue -u ${USER} # or squeue --me

The workshop account runs only a few 64-core jobs at the same time, and every learner submits one, so your job may stay pending (PD in squeue) for a day or more before its 5.5 hours start. When it has finished, the commands in the rest of this episode work on your own output.

Animation: What happens when you submit

Callout

If your Cell Ranger job dies

A job can be killed by its time limit, a node failure or a cluster outage. Resubmitting the same script then usually stops within seconds, because cellranger_output/pbmc10k/ already exists from the interrupted run. Read the end of cellranger_<jobid>.out and cellranger_<jobid>.err for the reason, then choose one:

Resume. Cell Ranger can continue an interrupted run from its last finished stage. Remove the lock file the killed run left behind and submit the same script again:
BASH
```
rm ${RCAC_SCRATCH}/scrna_workshop/cellranger_output/pbmc10k/_lock
sbatch ${RCAC_SCRATCH}/scrna_workshop/run_cellranger.sh
```

Start over. Move the old run aside and submit again:

BASH

mv ${RCAC_SCRATCH}/scrna_workshop/cellranger_output/pbmc10k \
   ${RCAC_SCRATCH}/scrna_workshop/cellranger_output/pbmc10k_failed
sbatch ${RCAC_SCRATCH}/scrna_workshop/run_cellranger.sh

Either way, the rest of the workshop does not wait for it: episode 3 uses the pre-computed matrix.

Cell Ranger count parameters

Parameter	Value	Description
`--id`	`pbmc10k`	A unique run ID. Cell Ranger creates an output directory with this name.
`--transcriptome`	`.../refdata-gex-GRCh38-2024-A`	Path to the 10x-compatible reference genome package containing the FASTA, GTF, and STAR index.
`--fastqs`	`.../fastq`	Directory containing the FASTQ files. Cell Ranger auto-detects files matching the sample name.
`--create-bam`	`true`	Generate a BAM file with aligned reads. Starting with Cell Ranger v8, BAM output is off by default; set to `true` if you need the BAM for downstream tools such as velocyto or variant calling.
`--sample`	`pbmc_10k_v3`	The sample name prefix in the FASTQ file names. Cell Ranger selects files matching this prefix.
`--localcores`	`${SLURM_CPUS_ON_NODE}`	Number of CPU cores to use. We set this to match the SLURM allocation so Cell Ranger does not try to use more cores than allocated.
`--localmem`	`100`	Maximum memory (in GB) that Cell Ranger is allowed to use. Keep it below the memory SLURM gives the job, which on Negishi follows the number of cores requested.

Callout

What if I have multiple samples?

Q: I have 4 samples (WT_rep1, WT_rep2, KO_rep1, KO_rep2). Do I pass them all to one cellranger count command?

A: No. cellranger count processes one sample per invocation. Use a SLURM array job to run them in parallel. The template below is illustrative: replace the sample names with the prefixes of your own FASTQ files.

BASH

#!/bin/bash
#SBATCH --array=0-3
#SBATCH --nodes=1
#SBATCH --ntasks=16
#SBATCH --time=1-00:00:00
#SBATCH --job-name=cellranger
#SBATCH --account=rcac-rnaseq
#SBATCH --partition=cpu

ml --force purge
ml biocontainers cellranger/10.0.0

SAMPLES=(WT_rep1 WT_rep2 KO_rep1 KO_rep2)
SAMPLE=${SAMPLES[$SLURM_ARRAY_TASK_ID]}

cd ${RCAC_SCRATCH}/scrna_workshop/cellranger_output

cellranger count \
    --id=${SAMPLE} \
    --transcriptome=${RCAC_SCRATCH}/scrna_workshop/reference/refdata-gex-GRCh38-2024-A \
    --fastqs=${RCAC_SCRATCH}/scrna_workshop/fastq \
    --create-bam true \
    --sample=${SAMPLE} \
    --localcores=${SLURM_CPUS_ON_NODE} \
    --localmem=64

For most downstream analyses, merge samples in Seurat (e.g., with merge() or IntegrateLayers()) rather than using cellranger aggr.

Cell Ranger output files

When the job completes, Cell Ranger creates the following directory structure:

BASH

ls ${RCAC_SCRATCH}/scrna_workshop/cellranger_output/pbmc10k/outs/

The key output files are:

pbmc10k/outs/
├── filtered_feature_bc_matrix/     # The count matrix (filtered to real cells)
│   ├── barcodes.tsv.gz             # List of cell barcodes
│   ├── features.tsv.gz             # List of genes
│   └── matrix.mtx.gz               # Sparse count matrix (MatrixMarket format)
├── raw_feature_bc_matrix/          # Count matrix for ALL barcodes (including empty drops)
├── web_summary.html                # Interactive QC report
├── metrics_summary.csv             # Key metrics as a CSV file
├── possorted_genome_bam.bam        # Aligned reads (large file)
└── molecule_info.h5                # Per-molecule info (used for aggregation)

The filtered_feature_bc_matrix/ directory is the most important output. It contains the UMI count matrix with only cell-containing barcodes (empty droplets have been removed by Cell Ranger’s cell calling algorithm). This is the file we will load into Seurat in the next episode.

Animation: Raw matrix to filtered matrix

The three files inside use the MatrixMarket sparse format:

barcodes.tsv.gz – one cell barcode per line
features.tsv.gz – one gene per line (Ensembl ID, gene symbol, feature type)
matrix.mtx.gz – the sparse matrix with row indices (genes), column indices (cells), and UMI counts

Callout

Interpreting web_summary.html

The web_summary.html file is an interactive HTML report that summarizes the quality of your Cell Ranger run. Open it in a web browser (you can copy it to your local machine with scp or view it through Open OnDemand). Here are the key metrics to check:

Metric	Healthy range	What it means
Estimated Number of Cells	Matches expectations (e.g., ~10,000)	The number of barcodes Cell Ranger identified as real cells.
Median Genes per Cell	>1,000 for most tissues	How many genes are detected in a typical cell. Very low values may indicate poor capture.
Reads Mapped Confidently to Transcriptome	>50%	Fraction of reads that align to annotated genes. Low values suggest reference mismatch or sample contamination.
Sequencing Saturation	>50% for most experiments	Fraction of reads that are PCR duplicates. Low saturation means more sequencing would detect additional UMIs. High saturation means you have enough sequencing depth.
Fraction Reads in Cells	>70%	Fraction of reads assigned to cell-containing barcodes. Low values may indicate high ambient RNA.

If any metric appears in red or orange in the web summary, investigate before proceeding to downstream analysis.

Interpreting Cell Ranger output

Screenshot of the Cell Ranger web_summary.html report for the PBMC 10k dataset showing key QC metrics

The web_summary.html file is the first thing to check after every Cell Ranger run. Here are the key QC metrics from our PBMC 10k run:

Metric	Value	What to look for
Estimated Number of Cells	11,809	Close to expected loading (~10k)
Median Genes per Cell	3,371	>1,000 for PBMCs is good
Reads Mapped Confidently to Transcriptome	79.1%	>70% expected for 3’ GEX
Sequencing Saturation	68.2%	>60% adequate; >80% ideal
Fraction Reads in Cells	95.7%	>90% indicates clean cell calling
Valid Barcodes	97.3%	>95% expected

A low Fraction Reads in Cells or low Median Genes per Cell would indicate problems with cell viability or library quality.

Callout

Sequencing saturation of 68.2%: is that enough?

A saturation of 68.2% means that roughly 32% of additional sequencing would yield new UMIs. For most differential expression and cell type identification workflows, this depth is sufficient. Diminishing returns set in above ~80% saturation, so deeper sequencing of this library would provide only marginal gains.

Processing with STARsolo

STARsolo is a single-cell RNA-seq processing mode built into the STAR aligner. It replicates the Cell Ranger pipeline – demultiplexing, alignment, barcode error correction, UMI deduplication, and cell filtering – but runs significantly faster because STAR is a highly optimized aligner. STARsolo is fully open source and produces results that are nearly identical to Cell Ranger.

Building a STAR genome index

STARsolo requires a STAR genome index. You can either build one from scratch or reuse the index that ships inside the Cell Ranger reference package. For this workshop, we will build a fresh index from the genome FASTA and GTF files included in the reference package. This ensures compatibility with the version of STAR installed on the cluster.

Callout

Reusing the Cell Ranger index

If you want to skip the index build, you can point STARsolo directly at the pre-built index inside the Cell Ranger reference: ${RCAC_SCRATCH}/scrna_workshop/reference/refdata-gex-GRCh38-2024-A/star/. However, this only works if the STAR version on the cluster matches the version Cell Ranger used to build the index. If the versions differ, STAR will exit with an error and you will need to build a new index.

Create the index build script:

BASH

cat << 'EOF' > ${RCAC_SCRATCH}/scrna_workshop/build_star_index.sh
#!/bin/bash
#SBATCH --nodes=1
#SBATCH --ntasks=64
#SBATCH --time=1-00:00:00
#SBATCH --job-name=star_index
#SBATCH --account=rcac-rnaseq
#SBATCH --partition=cpu
#SBATCH --output=star_index_%j.out
#SBATCH --error=star_index_%j.err

# Load modules
ml --force purge
ml biocontainers
ml star

REF=${RCAC_SCRATCH}/scrna_workshop/reference/refdata-gex-GRCh38-2024-A
STAR_INDEX=${RCAC_SCRATCH}/scrna_workshop/reference/star_index

mkdir -p ${STAR_INDEX}
gunzip -c ${REF}/genes/genes.gtf.gz > ${REF}/genes/genes.gtf
STAR \
    --runMode genomeGenerate \
    --runThreadN ${SLURM_CPUS_ON_NODE} \
    --genomeDir ${STAR_INDEX} \
    --genomeFastaFiles ${REF}/fasta/genome.fa \
    --sjdbGTFfile ${REF}/genes/genes.gtf
EOF

Check it as you did the Cell Ranger script:

BASH

sbatch --test-only ${RCAC_SCRATCH}/scrna_workshop/build_star_index.sh

After the workshop, submit it for real:

BASH

sbatch ${RCAC_SCRATCH}/scrna_workshop/build_star_index.sh

This step requires approximately 35 GB of RAM and takes about 20 minutes on 64 cores. The resulting index will be written to ${RCAC_SCRATCH}/scrna_workshop/reference/star_index/.

Running the STARsolo alignment

Once the index is ready, create the STARsolo alignment script. Note an important difference from Cell Ranger: STAR expects the cDNA read (R2) first, followed by the barcode read (R1) in the --readFilesIn argument.

BASH

cat << 'EOF' > ${RCAC_SCRATCH}/scrna_workshop/run_starsolo.sh
#!/bin/bash
#SBATCH --nodes=1
#SBATCH --ntasks=64
#SBATCH --time=1-00:00:00
#SBATCH --job-name=starsolo
#SBATCH --account=rcac-rnaseq
#SBATCH --partition=cpu
#SBATCH --output=starsolo_%j.out
#SBATCH --error=starsolo_%j.err

# Load modules
ml --force purge
ml biocontainers
ml star

# Define paths
FASTQ=${RCAC_SCRATCH}/scrna_workshop/fastq
STAR_INDEX=${RCAC_SCRATCH}/scrna_workshop/reference/star_index
WHITELIST=${RCAC_SCRATCH}/scrna_workshop/reference/3M-february-2018.txt

cd ${RCAC_SCRATCH}/scrna_workshop/starsolo_output

STAR \
    --runThreadN ${SLURM_CPUS_ON_NODE} \
    --genomeDir ${STAR_INDEX} \
    --sjdbGTFfile ${RCAC_SCRATCH}/scrna_workshop/reference/refdata-gex-GRCh38-2024-A/genes/genes.gtf \
    --readFilesIn \
        ${FASTQ}/pbmc_10k_v3_S1_L001_R2_001.fastq.gz,${FASTQ}/pbmc_10k_v3_S1_L002_R2_001.fastq.gz \
        ${FASTQ}/pbmc_10k_v3_S1_L001_R1_001.fastq.gz,${FASTQ}/pbmc_10k_v3_S1_L002_R1_001.fastq.gz \
    --readFilesCommand zcat \
    --soloType CB_UMI_Simple \
    --soloCBwhitelist ${WHITELIST} \
    --soloCBstart 1 \
    --soloCBlen 16 \
    --soloUMIstart 17 \
    --soloUMIlen 12 \
    --soloBarcodeMate 0 \
    --soloCBmatchWLtype 1MM_multi_Nbase_pseudocounts \
    --soloUMIfiltering MultiGeneUMI_CR \
    --soloUMIdedup 1MM_CR \
    --clipAdapterType CellRanger4 \
    --soloCellFilter EmptyDrops_CR \
    --outSAMtype BAM SortedByCoordinate \
    --outSAMattributes CR UR CY UY CB UB \
    --outFileNamePrefix starsolo_pbmc10k_
EOF

Check it:

BASH

sbatch --test-only ${RCAC_SCRATCH}/scrna_workshop/run_starsolo.sh

After the workshop, once the index job has finished, submit it for real:

BASH

sbatch ${RCAC_SCRATCH}/scrna_workshop/run_starsolo.sh

STARsolo parameter reference

Parameter	Value	Description
`--runThreadN`	`${SLURM_CPUS_ON_NODE}`	Number of threads. Matches the SLURM allocation.
`--genomeDir`	`.../star_index`	Path to the STAR genome index directory.
`--sjdbGTFfile`	`.../genes/genes.gtf`	Gene annotation GTF file for splice junction detection and gene counting.
`--readFilesIn`	`R2_L001,R2_L002 R1_L001,R1_L002`	Input FASTQ files. Important: the cDNA read (R2) comes first, then the barcode read (R1). Multiple lanes are comma-separated within each read group; a space separates R2 from R1.
`--readFilesCommand`	`zcat`	Decompression command for gzipped FASTQ files.
`--soloType`	`CB_UMI_Simple`	Barcode geometry: a single cell barcode followed by a UMI on the barcode read. This is the standard layout for 10x Chromium libraries.
`--soloCBwhitelist`	`.../3M-february-2018.txt`	The 10x barcode whitelist file. This contains all 3.7 million valid cell barcode sequences for v3 chemistry.
`--soloCBstart`	`1`	Position where the cell barcode starts on the barcode read (1-based).
`--soloCBlen`	`16`	Length of the cell barcode in bases.
`--soloUMIstart`	`17`	Position where the UMI starts on the barcode read (immediately after the cell barcode).
`--soloUMIlen`	`12`	Length of the UMI in bases. For 10x v3 chemistry this is 12; for v2 chemistry it is 10.
`--soloBarcodeMate`	`0`	Which mate carries the barcode. `0` means the barcode is on a separate read (R1), not embedded in the cDNA read.
`--soloCBmatchWLtype`	`1MM_multi_Nbase_pseudocounts`	Barcode error correction strategy. Allows 1 mismatch when matching to the whitelist, handles N bases, and uses pseudocounts to resolve ambiguous corrections. This matches Cell Ranger’s correction algorithm.
`--soloUMIfiltering`	`MultiGeneUMI_CR`	Filters UMIs that map to multiple genes using the Cell Ranger algorithm.
`--soloUMIdedup`	`1MM_CR`	UMI deduplication strategy allowing 1 mismatch, matching Cell Ranger’s approach.
`--clipAdapterType`	`CellRanger4`	Clips adapter sequences using the same method as Cell Ranger 4+.
`--soloCellFilter`	`EmptyDrops_CR`	Cell calling algorithm. `EmptyDrops_CR` replicates Cell Ranger’s EmptyDrops-based method for distinguishing real cells from empty droplets.
`--outSAMtype`	`BAM SortedByCoordinate`	Output a coordinate-sorted BAM file.
`--outSAMattributes`	`CR UR CY UY CB UB`	Include cell barcode and UMI tags in the BAM file (raw and error-corrected versions plus quality scores).
`--outFileNamePrefix`	`starsolo_pbmc10k_`	Prefix for all output file names.

STARsolo output files

STARsolo writes the count matrices under Solo.out/ inside the output directory:

BASH

ls ${RCAC_SCRATCH}/scrna_workshop/starsolo_output/starsolo_pbmc10k_Solo.out/Gene/filtered/

Expected output

OUTPUT

barcodes.tsv  features.tsv  matrix.mtx

The filtered directory contains three files in the same format as Cell Ranger:

barcodes.tsv – cell barcodes that passed the EmptyDrops cell filter
features.tsv – gene identifiers (Ensembl ID and gene symbol)
matrix.mtx – the sparse UMI count matrix in MatrixMarket format

Note that STARsolo’s output files are uncompressed (.tsv and .mtx rather than .tsv.gz and .mtx.gz). Seurat’s Read10X() function handles both compressed and uncompressed formats automatically.

Callout

Runtime comparison

On 64 cores with the PBMC 10k v3 dataset, measured on Negishi (October 2026):

Cell Ranger: about 5.5 hours (peak memory about 75 GB)
STAR index: about 20 minutes (about 35 GB)
STARsolo: about 45 minutes (about 50 GB)

STARsolo is roughly 7x faster because STAR is a highly optimized C++ aligner, while Cell Ranger wraps STAR with additional Python and Java orchestration layers. Both produce nearly identical count matrices.

Interpreting STARsolo output

STARsolo writes alignment statistics to Log.final.out in the output directory. Here are the key metrics from our PBMC 10k run:

Metric	Value	What to look for
Uniquely mapped reads %	89.28%	>80% for human GEX
Multi-mapping reads %	7.14%	<10% typical
Unmapped: too short %	2.56%	>10% suggests trimming or degradation
Mismatch rate per base	0.56%	<1% expected
Mapping speed	945.4M reads/hr	Cell Ranger’s whole run averaged about 120M reads/hr (639M reads in 5.4 h)

Unlike Cell Ranger, STARsolo does not produce an HTML summary report. Check Log.final.out for alignment statistics and Solo.out/Gene/Summary.csv for cell-level metrics such as the number of cells detected and UMIs per cell.

Comparing Outputs

Once both tools have finished, compare their outputs (during the workshop, the instructor shows these numbers from the pre-workshop run). We can count the number of cells (barcodes) and genes (features) detected by each tool directly from the output files.

Cell Ranger:

BASH

echo "Cell Ranger cells:"
zcat ${RCAC_SCRATCH}/scrna_workshop/cellranger_output/pbmc10k/outs/filtered_feature_bc_matrix/barcodes.tsv.gz | wc -l

echo "Cell Ranger genes:"
zcat ${RCAC_SCRATCH}/scrna_workshop/cellranger_output/pbmc10k/outs/filtered_feature_bc_matrix/features.tsv.gz | wc -l

STARsolo:

BASH

echo "STARsolo cells:"
wc -l < ${RCAC_SCRATCH}/scrna_workshop/starsolo_output/starsolo_pbmc10k_Solo.out/Gene/filtered/barcodes.tsv

echo "STARsolo genes:"
wc -l < ${RCAC_SCRATCH}/scrna_workshop/starsolo_output/starsolo_pbmc10k_Solo.out/Gene/filtered/features.tsv

Expected output

Both tools should report very similar numbers:

OUTPUT

Cell Ranger cells:
11809
Cell Ranger genes:
38606

STARsolo cells:
11682
STARsolo genes:
38606

The cell counts may differ by a small amount (typically <1%) due to minor differences in the EmptyDrops cell-calling implementation. The gene count is identical because both tools use the same reference annotation.

When to use each tool

Use Cell Ranger when:

You need results that are directly comparable to 10x Genomics documentation and publications
You want the interactive web_summary.html report
Your institution has Cell Ranger installed and you don’t need to worry about licensing
You are running the standard 10x Chromium workflow with no custom modifications

Use STARsolo when:

You need faster turnaround (especially important for large datasets or parameter sweeps)
You need a fully open-source pipeline
You are working with non-standard barcode configurations or custom protocols
You want fine-grained control over alignment and counting parameters

For this workshop, the pre-computed filtered_feature_bc_matrix/ in the workshop data directory was generated with Cell Ranger. We will use it for all downstream analysis starting in the next episode.

Challenge

Challenge 1: Modify STARsolo for v2 Chemistry

The PBMC 10k dataset uses 10x Chromium v3 chemistry (16 bp barcode + 12 bp UMI). Older experiments used v2 chemistry, which has a different UMI length and barcode whitelist.

What parameters in the STARsolo command would you change to process a v2 chemistry library? (Hint: two parameters change and one additional parameter needs a different file.)

Show me the solution

Three changes are required for v2 chemistry:

--soloUMIlen 10 instead of 12. v2 chemistry uses 10 bp UMIs (compared to 12 bp in v3).
--soloCBwhitelist must point to the v2 whitelist file (737K-august-2016.txt) instead of the v3 whitelist (3M-february-2018.txt). The v2 whitelist contains ~737,000 valid barcodes; the v3 whitelist contains ~3.7 million.
--soloCBlen 16 stays the same – both v2 and v3 use 16 bp cell barcodes. The cell barcode length did not change between chemistry versions.

The updated parameters would look like:

BASH

    --soloCBwhitelist /path/to/737K-august-2016.txt \
    --soloCBlen 16 \
    --soloUMIstart 17 \
    --soloUMIlen 10 \

Always check the chemistry version in your experiment’s documentation or the web_summary.html from a previous Cell Ranger run. Using the wrong whitelist or UMI length will result in very few cells being detected.

Challenge

Challenge 2: QC from Cell Ranger Summary

A colleague shares their Cell Ranger metrics_summary.csv with the following values. Review the metrics and identify which one suggests a potential quality concern.

Metric	Value
Estimated Number of Cells	8,500
Mean Reads per Cell	45,000
Median Genes per Cell	2,100
Fraction Reads in Cells	55.2%
Reads Mapped Confidently to Transcriptome	62.5%
Sequencing Saturation	28.3%

Show me the solution

Two metrics stand out:

Fraction Reads in Cells (55.2%) is low. A healthy value is typically >70%. When nearly half of all reads come from non-cell barcodes, it indicates a high level of ambient RNA (free-floating mRNA from lysed cells in the droplet suspension). This does not necessarily mean the experiment failed, but it means a large fraction of the sequencing budget was spent on background noise. Downstream tools like SoupX or CellBender can help correct for ambient RNA contamination.
Sequencing Saturation (28.3%) is low. This means that 71.7% of the reads are still detecting new UMIs – in other words, additional sequencing would reveal more unique transcripts. A value below 50% suggests the library was under-sequenced. If the budget allows, resequencing the same library would improve gene detection. However, for some applications (e.g., cell type identification), 28% saturation may still be sufficient.

The other metrics look reasonable: 8,500 cells is within expectations for a standard 10x run, median genes per cell >2,000 indicates decent capture quality, and 62.5% reads mapped to transcriptome is acceptable for a human sample.

Key Points

Workshop data should be organized under ${RCAC_SCRATCH}/scrna_workshop/ with separate directories for FASTQ files, reference, and outputs
Cell Ranger count aligns 10x FASTQ files and produces a filtered count matrix, a web summary report, and alignment metrics
STARsolo replicates Cell Ranger’s pipeline but runs roughly 7x faster and is fully open source
sbatch --test-only checks a job script’s resource request without submitting it
Both tools produce a filtered count matrix in the same format: barcodes.tsv, features.tsv, and matrix.mtx
Always inspect the Cell Ranger web_summary.html to check estimated cell count, median genes per cell, sequencing saturation, and fraction of reads in cells before proceeding

Content from Quality Control

Last updated on 2026-10-27 | Edit this page

Estimated time: 60 minutes

Overview

Questions

How do we load a 10x count matrix into R and create a Seurat object?
What quality metrics should we examine for scRNA-seq data?
How do we choose appropriate filtering thresholds?
How does mitochondrial gene percentage help identify dying cells?
What are doublets and how can we detect them?

Objectives

Load a Cell Ranger filtered count matrix into a Seurat v5 object
Calculate and visualize QC metrics: nFeature_RNA, nCount_RNA, and percent.mt
Apply filtering thresholds to remove low-quality cells
Describe what doublets are and when dedicated detection tools are needed

Prerequisite

Prerequisites

This episode requires an RStudio session on the Negishi cluster. Launch RStudio (bioconductor) under Bioinformatics Apps on Open OnDemand as described in the Setup instructions.

Setup

R

library(Seurat)
library(ggplot2)
library(patchwork)

Loading Data into Seurat

The starting point for this episode is the filtered count matrix produced by Cell Ranger. This matrix lives in a filtered_feature_bc_matrix/ directory and consists of three files: barcodes.tsv.gz, features.tsv.gz, and matrix.mtx.gz. Seurat’s Read10X() function reads all three files and assembles them into a sparse matrix in R.

Your Cell Ranger job from the previous episode takes several hours, so we use the copy of the same matrix that you copied from the workshop data directory during setup (scrna_workshop/filtered_feature_bc_matrix/). It was produced by the run_cellranger.sh script from the previous episode.

First, set up the paths so that all file references in this episode use the same base directory:

R

data_dir <- paste0(
    "/scratch/negishi/", Sys.getenv("USER"),
    "/scrna_workshop/filtered_feature_bc_matrix/"
)
work_dir <- paste0(
    "/scratch/negishi/", Sys.getenv("USER"),
    "/scrna_workshop/"
)
setwd(work_dir)

Callout

Using your own Cell Ranger output

If your Cell Ranger job has finished, you can load your own matrix instead by setting data_dir to /scratch/negishi/<username>/scrna_workshop/cellranger_output/pbmc10k/outs/filtered_feature_bc_matrix/. The results in this episode are the same.

Now load the count matrix:

R

pbmc.data <- Read10X(data.dir = data_dir)

Read10X() returns a sparse matrix where rows are genes and columns are cell barcodes. Let’s check its dimensions before creating the Seurat object:

R

dim(pbmc.data)

Expected output

OUTPUT

[1] 38606 11809

This tells us the matrix contains 38,606 genes (rows) and 11,809 cell barcodes (columns).

Now we create a Seurat object, which is the central data structure for all downstream analysis. The CreateSeuratObject() function takes the raw count matrix and applies two initial filters:

min.cells = 3 removes genes detected in fewer than 3 cells (these are too rare to be informative)
min.features = 200 removes cells with fewer than 200 detected genes (these are likely empty droplets or debris)

R

pbmc <- CreateSeuratObject(
    counts = pbmc.data,
    project = "pbmc10k",
    min.cells = 3,
    min.features = 200
)
pbmc

Expected output

OUTPUT

An object of class Seurat
29155 features across 11721 samples within 1 assay
Active assay: RNA (29155 features, 0 variable features)
 1 layer present: counts

After the initial filters, we have 29,155 genes and 11,721 cells. The gene count dropped from 38,606 because many genes were detected in fewer than 3 cells and were removed.

Parameter	Value	Description
`counts`	`pbmc.data`	The raw UMI count matrix (genes x cells).
`project`	`"pbmc10k"`	A name for this project, stored in the object metadata.
`min.cells`	`3`	Keep only genes detected in at least 3 cells. Removes very rare genes.
`min.features`	`200`	Keep only cells with at least 200 detected genes. Removes likely empty droplets.

Understanding the Seurat v5 object

The Seurat object stores all data and analysis results in an organized structure. The key components are:

Assays: containers for expression data. Our object has one assay called RNA. Each assay can hold multiple layers of the data at different processing stages.
Layers: the expression matrices within an assay:
- counts – the raw UMI counts (populated now)
- data – the normalized expression values (populated after normalization)
- scale.data – the scaled/centered values (populated after scaling)
Metadata: a data frame with one row per cell, storing QC metrics, cluster assignments, and any other per-cell annotations.

Diagram of the Seurat v5 object showing the RNA assay with its three layers (counts, data, scale.data) and the per-cell metadata columns

Animation: The Seurat object through the day

You can inspect the metadata that Seurat automatically computed during object creation:

R

head(pbmc@meta.data)

Expected output

OUTPUT

                   orig.ident nCount_RNA nFeature_RNA
AAACCCAAGCGCCCAT-1    pbmc10k       4282         2152
AAACCCAAGGTTCCGC-1    pbmc10k      29509         5990
AAACCCACAGACAAGC-1    pbmc10k        574          358
AAACCCACAGAGTTGG-1    pbmc10k       8400         2881
AAACCCACAGGTATGG-1    pbmc10k       9675         3731
AAACCCACATAGTCAC-1    pbmc10k      10058         3365

Seurat automatically calculated two metrics per cell:

nCount_RNA – the total number of UMIs in that cell
nFeature_RNA – the number of unique genes detected in that cell

Callout

Seurat v5 layer syntax

Seurat v5 introduced a new layer-based architecture. To access expression data, use the $ operator on the assay:

R

# Seurat v5 syntax (use this)
pbmc[["RNA"]]$counts     # raw counts
pbmc[["RNA"]]$data       # normalized data (after NormalizeData)
pbmc[["RNA"]]$scale.data # scaled data (after ScaleData)

If you see older tutorials using the @ slot syntax, that is Seurat v4 and will not work correctly with v5 objects:

R

# Seurat v4 syntax (do NOT use)
pbmc@assays$RNA@counts   # deprecated

QC Metrics

Not every barcode in the count matrix represents a healthy cell. Some barcodes correspond to dying cells, empty droplets that slipped through Cell Ranger’s filter, or doublets (two cells in one droplet). Before we can analyze the data, we need to identify and remove these low-quality barcodes.

We evaluate cell quality using three metrics:

Metric	What it measures	Low-quality signal
`nFeature_RNA`	Number of unique genes detected per cell	Very low: empty droplet or debris. Very high: possible doublet.
`nCount_RNA`	Total UMI counts per cell	Very low: failed capture. Very high: possible doublet.
`percent.mt`	Percentage of UMIs from mitochondrial genes	High (>15–20%): dying or stressed cell.

Let’s calculate the mitochondrial percentage. Human mitochondrial gene names start with MT- (e.g., MT-CO1, MT-ND1, MT-ATP6), so we use a pattern match:

R

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

Now let’s inspect the distribution of all three metrics. First, we look at the summary statistics:

R

summary(pbmc$nFeature_RNA)
summary(pbmc$nCount_RNA)
summary(pbmc$percent.mt)

Expected output

OUTPUT

# nFeature_RNA
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    201    2895    3378    3557    4270    9540

# nCount_RNA
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    505    8047   10357   12280   15166  103741

# percent.mt
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
 0.1752  4.8509  5.7882  6.7695  7.0114 89.8376

Most cells have 2,900-4,300 genes (IQR), 8,000-15,000 UMIs (IQR), and 4.9-7.0% mitochondrial content (IQR). But notice the maximum values: one cell has 9,540 genes (possible doublet), another has over 103,000 UMIs, and one cell has nearly 90% mitochondrial reads (clearly dying). These outliers are what QC filtering removes.

Violin plots

Violin plots show the distribution of each metric across all cells. They are the standard first visualization for scRNA-seq QC.

R

VlnPlot(pbmc,
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
        ncol = 3,
        pt.size = 0,
        layer = "counts")

Violin plots of nFeature_RNA, nCount_RNA, and percent.mt showing the distribution of QC metrics across all cells before filtering

Setting pt.size = 0 hides individual points so the distribution shape is easier to see with 10,000+ cells. If you want to see points, set pt.size = 0.1.

Scatter plots

Scatter plots reveal relationships between metrics. The most informative plot is nCount_RNA vs. nFeature_RNA, which should show a strong positive correlation: cells with more total UMIs generally detect more genes. Points that deviate from this trend are suspicious.

R

ggplot(pbmc@meta.data, aes(x = nCount_RNA, y = nFeature_RNA, color = percent.mt)) +
    geom_point(size = 0.5, alpha = 0.4) +
    scale_color_gradient(low = "grey90", high = "firebrick", name = "% MT") +
    theme_minimal() +
    labs(x = "Total UMI counts", y = "Number of genes")

Scatter plot of total UMI counts versus number of genes per cell colored by mitochondrial percentage, with high-mito cells highlighted in red

What to look for in this plot:

Main cloud: most cells should form a dense cloud along a rough diagonal. These are healthy cells with a proportional relationship between total RNA content and gene diversity.
Bottom-left outliers: cells with very few UMIs and very few genes. These are empty droplets or debris that slipped through the Cell Ranger filter.
Top-right outliers: cells with both very high UMI counts and very high gene counts. These are potential doublets – two cells captured in one droplet, producing double the RNA content.
High mitochondrial cells: points colored in darker shades on the grey-to-red gradient. These often sit below the main cloud because the cell is losing cytoplasmic RNA while retaining mitochondrial RNA.

Callout

What do dying cells look like?

When a cell is damaged or dying, its cell membrane becomes leaky. Cytoplasmic mRNA molecules – which are small and not membrane-bound – escape through the holes in the membrane. However, mitochondrial mRNA is protected inside the double-membraned mitochondria and does not leak out as readily.

The result is a cell with:

High percent.mt: the remaining RNA is disproportionately mitochondrial
Low nFeature_RNA: many cytoplasmic transcripts have leaked out
Low nCount_RNA: total RNA content is reduced

This is why mitochondrial percentage is such a reliable indicator of cell quality. A healthy PBMC typically has 3–8% mitochondrial content. Cells with >15% are likely damaged, and cells with >20% are almost certainly compromised.

Animation: What a dying cell looks like in QC

Setting Thresholds

Now we apply filters to remove low-quality cells. There are no universal thresholds that work for every dataset. The right cutoffs depend on the tissue type, species, experimental protocol, and sequencing depth. The best approach is to inspect the distributions (as we just did) and choose thresholds that remove clear outliers while preserving the bulk of the data.

For this PBMC dataset, we apply the following thresholds:

nFeature_RNA > 200: remove cells with very few genes (likely empty)
nFeature_RNA < 7000: remove cells with an extreme number of genes (possible doublets)
nCount_RNA < 50000: remove cells with unusually high UMI counts (possible doublets)
percent.mt < 15: remove cells with high mitochondrial content (likely dying)

Let’s visualize where these cutoffs fall on the violin plots. Red dashed lines mark the thresholds:

R

p1 <- VlnPlot(pbmc, features = "nFeature_RNA", pt.size = 0.1, layer = "counts") +
    geom_hline(yintercept = c(205,7000), linetype = "dashed", color = "red") +
    NoLegend()
p2 <- VlnPlot(pbmc, features = "nCount_RNA", pt.size = 0.1, layer = "counts") +
    geom_hline(yintercept = 50000, linetype = "dashed", color = "red") +
    NoLegend()
p3 <- VlnPlot(pbmc, features = "percent.mt", pt.size = 0.1, layer = "counts") +
    geom_hline(yintercept = 15, linetype = "dashed", color = "red") +
    NoLegend()
p1 + p2 + p3 + plot_layout(ncol = 3)

Violin plots of nFeature_RNA, nCount_RNA, and percent.mt with red dashed lines showing the filtering thresholds and individual cells visible as points

Cells above or below the red lines will be removed. The small point size (pt.size = 0.1) lets you see how many individual cells fall outside the thresholds.

Animation: Choose the thresholds, watch the cells

Let’s record how many cells we have before filtering:

R

cat("Cells before filtering:", ncol(pbmc), "\n")

Apply the filters:

R

pbmc <- subset(pbmc,
               subset = nFeature_RNA > 200 &
                        nFeature_RNA < 7000 &
                        nCount_RNA < 50000 &
                        percent.mt < 15)
cat("Cells after filtering:", ncol(pbmc), "\n")

Expected output

OUTPUT

Cells before filtering: 11721
Cells after filtering: 11310

Filtering removes ~400 cells (~3.5% of the total). The exact number depends on minor differences in the Cell Ranger version and reference used. You should have roughly 11,300 cells remaining.

Let’s verify the filtering by looking at the QC distributions again:

R

VlnPlot(pbmc,
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
        ncol = 3,
        pt.size = 0,
        layer = "counts")

Violin plots of nFeature_RNA, nCount_RNA, and percent.mt after filtering, showing tighter distributions with extreme outliers removed

The violin plots should now show tighter distributions with the extreme tails removed. The mitochondrial percentage should be capped below 15%.

Let’s also verify with a scatter plot:

R

ggplot(pbmc@meta.data, aes(x = nCount_RNA, y = nFeature_RNA, color = percent.mt)) +
    geom_point(size = 0.5, alpha = 0.4) +
    scale_color_gradient(low = "grey90", high = "firebrick", name = "% MT") +
    theme_minimal() +
    labs(x = "Total UMI counts", y = "Number of genes")

Scatter plot of total UMI counts versus number of genes after filtering, colored by mitochondrial percentage, showing a tighter cloud with outliers removed

Finally, save the filtered object for the next episode:

R

saveRDS(pbmc, file = "pbmc_filtered.rds")

Callout

MAD-based filtering

Instead of choosing fixed thresholds by visual inspection, a more principled approach is MAD-based filtering (Median Absolute Deviation). The idea is to compute the median and MAD of each QC metric, then flag any cell that falls more than a set number of MADs from the median as an outlier.

The formula for an upper threshold is:

threshold = median + N x MAD

and for a lower threshold:

threshold = median - N x MAD

where N is typically 3 or 5. For example, a cell is flagged as high mitochondrial if:

percent.mt > median(percent.mt) + 3 x MAD(percent.mt)

This approach is adaptive: it adjusts automatically to the distribution of each dataset, so you don’t need to guess appropriate cutoffs. It is implemented in the scuttle Bioconductor package via the isOutlier() function.

For this workshop, we use fixed thresholds for simplicity and transparency, but MAD-based filtering is recommended for production analyses where you process many datasets with varying quality.

Doublet Detection

Callout

Doublets in scRNA-seq

A doublet occurs when two cells are captured in the same GEM droplet and tagged with the same cell barcode. The resulting “cell” has roughly double the RNA content of a real single cell, producing a hybrid expression profile that blends two cell types.

How doublets appear in QC metrics: Doublets tend to have high nCount_RNA and high nFeature_RNA because they contain RNA from two cells. Our upper threshold on nFeature_RNA < 7000 catches some of the most extreme doublets, but it cannot identify doublets between similar cell types (e.g., two T cells) because their combined profile looks like a normal high-quality T cell.

Dedicated doublet detection tools such as scDblFinder and DoubletFinder use a simulation-based approach: they create artificial doublets by randomly combining pairs of real cells, then train a classifier to identify real cells that resemble these artificial doublets.

When dedicated detection matters most:

High cell loading (>10,000 cells targeted) where doublet rates exceed 5–8%
Multiplexed experiments (e.g., cell hashing) with cross-sample doublets
Studies where rare intermediate populations could be confused with doublets

For this dataset, 10x Genomics expects about 0.8% doublets per 1,000 cells recovered with v3 chemistry, so about 8–9% for the roughly 11,800 cells here. Many of them combine two cells of the same type (two T cells, say), which barely changes clustering. Doublets of two different types tend to form small, mixed clusters that can be found during cell type annotation (cluster 16 in episode 6, with both monocyte and B cell markers, is one). We will not run a dedicated doublet detection tool in this workshop, but be aware of these tools for your own analyses.

Animation: Why some doublets hide

Challenge

Challenge 1: Propose Your Own Thresholds

Look at the violin plots and scatter plots you generated above. Based on the distributions, propose your own filtering thresholds for nFeature_RNA, nCount_RNA, and percent.mt. Justify your choices. How many cells pass your filters compared to the 11,310 cells that passed the thresholds used above?

Run the following code with your chosen values:

R

# Replace the ??? with your thresholds
pbmc_custom <- subset(pbmc,
                      subset = nFeature_RNA > ??? &
                               nFeature_RNA < ??? &
                               percent.mt < ???)
cat("Cells with custom thresholds:", ncol(pbmc_custom), "\n")

Show me the solution

There is no single correct answer – reasonable thresholds depend on how conservative you want to be. Here is one example of a slightly more conservative approach:

R

pbmc_custom <- subset(pbmc,
                      subset = nFeature_RNA > 300 &
                               nFeature_RNA < 4000 &
                               percent.mt < 10)
cat("Cells with custom thresholds:", ncol(pbmc_custom), "\n")

This retains 7,579 cells compared to 11,310 with the thresholds used above – a loss of about a third of the data. The trade-offs are:

nFeature_RNA > 300 (vs. 200): slightly more aggressive removal of low-complexity cells, but may also remove some small cell types like platelets that naturally express few genes.
nFeature_RNA < 4000 (vs. 7000): this is the biggest driver of cell loss here. The median nFeature_RNA in this dataset is ~3,378, so a 4,000 ceiling cuts into the upper quartile and removes a substantial number of valid cells, not just doublets.
percent.mt < 10 (vs. 15): removes more potentially stressed cells, but some immune cell types (e.g., activated T cells) naturally have moderately elevated mitochondrial content.

The key point is that more aggressive filtering is not always better. Every cell you remove is data you lose. Dropping from 11,310 to 7,579 cells means losing a third of the dataset, which could reduce statistical power and under-represent cell types that naturally have higher gene counts or mitochondrial content. The goal is to remove cells that would distort the analysis, not to enforce an artificially narrow definition of “high quality.”

Challenge

Challenge 2: Aggressive mitochondrial filtering

What happens if you set percent.mt < 5 instead of percent.mt < 15? Run the following code and compare the number of remaining cells. Would this threshold be too aggressive for PBMCs?

R

# pbmc is the filtered object from above (the episode's thresholds, percent.mt < 15)
pbmc_strict <- subset(pbmc, subset = percent.mt < 5)
cat("Cells with percent.mt < 15:", ncol(pbmc), "\n")
cat("Cells with percent.mt < 5:", ncol(pbmc_strict), "\n")

Show me the solution

With percent.mt < 5, only 3,276 cells remain, compared to 11,310 with percent.mt < 15. That is a loss of 71% of the dataset.

This is far too aggressive for PBMCs. Here is why:

The median mitochondrial percentage in this dataset is ~5.8%, meaning a < 5 cutoff removes the majority of cells by definition – not just dying cells, but most of the healthy population.
Different immune cell types have naturally different mitochondrial content, so the loss is uneven. Using the labels from episode 6, a < 5 cutoff keeps about 45% of CD4 T and NK cells but only about 24% of CD14 monocytes and 12% of CD8 T cells and FCGR3A+ monocytes. It biases the downstream analysis toward the cell types with low mitochondrial content.
A threshold of < 5 might be appropriate for some cell lines or tissues where mitochondrial content is uniformly low, but for primary immune cells from blood, 10-15% is a more appropriate cutoff.

As a general rule, look at the distribution of percent.mt in your data. If there is a clear separation between a low-mt peak and a high-mt tail, set your threshold in the valley between them. If the distribution is unimodal (one smooth peak), a fixed cutoff based on MAD (see the callout above) is safer than an arbitrary round number.

Key Points

Read10X() loads a Cell Ranger count matrix and CreateSeuratObject() creates the central data structure for Seurat analysis
The three key QC metrics are nFeature_RNA (genes per cell), nCount_RNA (UMIs per cell), and percent.mt (mitochondrial percentage)
High mitochondrial percentage indicates damaged cells that are losing cytoplasmic RNA through a leaky membrane
Filtering thresholds should be chosen by inspecting the data distributions, not by applying universal fixed cutoffs
Dedicated doublet detection tools (scDblFinder, DoubletFinder) are available for high-loading experiments but are not always necessary at standard loading rates

Content from Normalization and Feature Selection

Last updated on 2026-10-27 | Edit this page

Estimated time: 55 minutes

Overview

Questions

Why do we need to normalize scRNA-seq data before comparing cells?
How does log normalization work and where are the results stored?
What are highly variable features and why do we select them?
What does scaling do and when should we regress out confounders?
How does SCTransform differ from the LogNormalize workflow?

Objectives

Apply log normalization to correct for differences in sequencing depth between cells
Identify highly variable features using variance-stabilizing transformation
Scale the data to prepare for PCA
Describe SCTransform as an alternative to the three-step LogNormalize workflow

Prerequisite

Prerequisites

This episode requires an RStudio session on the Negishi cluster. Launch RStudio (bioconductor) under Bioinformatics Apps on Open OnDemand as described in the Setup instructions. You will also need the pbmc_filtered.rds object saved at the end of the previous episode.

Setup

R

library(Seurat)
library(ggplot2)
library(patchwork)

Loading the Filtered Data

We start from the filtered Seurat object saved at the end of the previous episode. This object contains only cells that passed our QC thresholds.

Set up the working directory to match the previous episode:

R

work_dir <- paste0(
    "/scratch/negishi/", Sys.getenv("USER"),
    "/scrna_workshop/"
)
setwd(work_dir)

Callout

Starting here without the previous episode?

Copy the saved object from the workshop data first: file.copy(paste0(work_dir, "rds/pbmc_filtered.rds"), work_dir)

R

pbmc <- readRDS(paste0(work_dir, "pbmc_filtered.rds"))
pbmc

Expected output

OUTPUT

An object of class Seurat
29155 features across 11310 samples within 1 assay
Active assay: RNA (29155 features, 0 variable features)
 1 layer present: counts

The object has 29,155 genes and 11,310 cells. Only the counts layer is populated at this stage. By the end of this episode we will fill in the data and scale.data layers as well.

Right now, the only expression data we have is raw UMI counts. Before we can compare gene expression between cells, we need to address a fundamental problem: cells are sequenced to different depths.

Consider two cells that are biologically identical – they express the same genes at the same levels. If Cell A was sequenced to 5,000 total UMIs and Cell B to 10,000 total UMIs, every gene in Cell B will appear to have roughly twice the counts of Cell A. This is purely a technical artifact of sequencing depth, not a biological difference. Normalization corrects for this so that expression values are comparable across cells.

Log Normalization

The standard normalization in Seurat is LogNormalize. It applies a simple three-step transformation to each cell independently:

Divide each gene’s count by the cell’s total UMI count
Multiply by a scale factor (default 10,000) so the values are not tiny fractions
Log-transform with log1p (natural log of 1 + x) to compress the dynamic range

The formula for a single gene g in cell c is:

normalized(g, c) = log(1 + count(g, c) / total_counts(c) x 10000)

Animation: Normalization corrects sequencing depth

R

pbmc <- NormalizeData(pbmc,
                      normalization.method = "LogNormalize",
                      scale.factor = 10000)

Parameter	Value	Description
`normalization.method`	`"LogNormalize"`	Divide by total counts, multiply by scale factor, log1p transform.
`scale.factor`	`10000`	Multiplicative factor applied after dividing by total counts. The default 10,000 means normalized values are in “counts per 10,000” before log transformation.

After normalization, the results are stored in the data layer of the RNA assay:

R

pbmc

Expected output

OUTPUT

An object of class Seurat
29155 features across 11310 samples within 1 assay
Active assay: RNA (29155 features, 0 variable features)
 2 layers present: counts, data

Notice that we now have 2 layers: counts (raw) and data (normalized).

You can access each layer directly:

R

# Raw counts for the first 5 genes in the first 3 cells
pbmc[["RNA"]]$counts[1:5, 1:3]

# Normalized values for the same genes and cells
pbmc[["RNA"]]$data[1:5, 1:3]

Expected Output

OUTPUT

# Raw counts for the first 5 genes in the first 3 cells
5 x 3 sparse Matrix of class "dgCMatrix"
                AAACCCAAGCGCCCAT-1 AAACCCAAGGTTCCGC-1 AAACCCACAGACAAGC-1
ENSG00000238009                  .                  .                  .
ENSG00000239945                  .                  .                  .
ENSG00000241860                  .                  .                  .
ENSG00000290385                  .                  .                  .
ENSG00000235146                  .                  .                  .

# Normalized values for the same genes and cells
5 x 3 sparse Matrix of class "dgCMatrix"
                AAACCCAAGCGCCCAT-1 AAACCCAAGGTTCCGC-1 AAACCCACAGACAAGC-1
ENSG00000238009                  .                  .                  .
ENSG00000239945                  .                  .                  .
ENSG00000241860                  .                  .                  .
ENSG00000290385                  .                  .                  .
ENSG00000235146                  .                  .                  .

The first few genes in the matrix are mostly zeros, so the output above is not very illuminating. To see the actual difference between raw and normalized values, we can filter for genes with non-zero counts first:

R

# Raw counts for the first 5 genes in the first 3 cells, with row sum > 0
mat <- pbmc[["RNA"]]$counts[, 1:3]
mat[rowSums(mat) > 0, ][1:5, ]

# Normalized values for the same genes and cells, with row sum > 0
mat2 <- pbmc[["RNA"]]$data[, 1:3]
mat2[rowSums(mat2) > 0, ][1:5, ]

Expected Output

OUTPUT

# Raw counts for the first 5 genes in the first 3 cells, with row sum > 0
5 x 3 sparse Matrix of class "dgCMatrix"
         AAACCCAAGCGCCCAT-1 AAACCCAAGGTTCCGC-1 AAACCCACAGACAAGC-1
ISG15                     1                  2                  .
C1orf159                  .                  1                  .
SDF4                      .                  2                  .
B3GALT6                   .                  1                  .
DVL1                      .                  1                  .

# Normalized values for the same genes and cells, with row sum > 0
5 x 3 sparse Matrix of class "dgCMatrix"
         AAACCCAAGCGCCCAT-1 AAACCCAAGGTTCCGC-1 AAACCCACAGACAAGC-1
ISG15               1.20458          0.5174592                  .
C1orf159            .                0.2918332                  .
SDF4                .                0.5174592                  .
B3GALT6             .                0.2918332                  .
DVL1                .                0.2918332                  .

Visualizing the effect of normalization

To see what normalization does, let’s compare the distribution of a gene before and after normalization. We’ll look at LYZ, a marker for CD14+ monocytes that has high and variable expression across cells.

R

# Extract raw and normalized values for LYZ
lyz_raw <- pbmc[["RNA"]]$counts["LYZ", ]
lyz_norm <- pbmc[["RNA"]]$data["LYZ", ]

p1 <- ggplot(data.frame(x = as.numeric(lyz_raw)), aes(x = x)) +
    geom_histogram(bins = 50, fill = "steelblue", color = "white") +
    labs(title = "LYZ - Raw counts", x = "UMI counts", y = "Number of cells") +
    theme_minimal()

p2 <- ggplot(data.frame(x = as.numeric(lyz_norm)), aes(x = x)) +
    geom_histogram(bins = 50, fill = "darkred", color = "white") +
    labs(title = "LYZ - Log-normalized", x = "Normalized expression", y = "Number of cells") +
    theme_minimal()

p1 + p2

Two histograms for LYZ. Raw counts show a zero-dominated right-skewed distribution. Log-normalized values reveal a bimodal pattern with peaks near 0.5 and 4.5, separating non-monocytes from LYZ-expressing monocytes

The raw count distribution is heavily right-skewed with a large spike at zero (cells that don’t express LYZ – mostly non-monocytes). After log normalization, the non-zero values are spread more evenly and the influence of sequencing depth differences is reduced.

Identifying Variable Features

The human genome contains roughly 20,000 protein-coding genes, but not all of them are informative for distinguishing cell types. Many genes are:

Housekeeping genes expressed at similar levels in every cell (e.g., ACTB, GAPDH). These carry no information about cell identity.
Low-expression genes detected in only a handful of cells. These are too noisy to contribute meaningfully.
Constant genes that show little variation across cells.

Highly variable features (HVGs) are genes that vary substantially from cell to cell. These are the genes most likely to capture biological differences between cell types. By selecting the top 2,000 variable features, we focus downstream analysis (PCA, clustering) on the genes that matter most, while reducing noise and computation time.

R

pbmc <- FindVariableFeatures(pbmc,
                             selection.method = "vst",
                             nfeatures = 2000)

Parameter	Value	Description
`selection.method`	`"vst"`	Variance-stabilizing transformation. Fits a mean-variance relationship across genes using local regression, then selects genes with the highest standardized variance (i.e., genes that vary more than expected given their expression level).
`nfeatures`	`2000`	Number of variable features to select. 2,000 is the standard default for most scRNA-seq analyses.

Animation: Which genes vary more than expected

Let’s see which genes were selected:

R

head(VariableFeatures(pbmc), 10)

Expected output

OUTPUT

 [1] "PPBP"            "PF4"             "LINC01478"       "LINC01374"       "PTGDS"
 [6] "JCHAIN"          "GP1BB"           "SOX5"            "ENSG00000225885" "EREG"

The top variable features include PPBP and PF4 (platelet markers), GP1BB (also platelets), JCHAIN (B/plasma cells), and PTGDS (dendritic cells). Several lincRNAs and an Ensembl ID also appear due to sparse but extreme expression in a few cells. These genes are highly expressed in specific cell types but absent in others, producing high cell-to-cell variance.

Visualize the variable feature selection with a mean-variance plot:

R

vf_plot <- VariableFeaturePlot(pbmc)
top10 <- head(VariableFeatures(pbmc), 10)
LabelPoints(plot = vf_plot, points = top10, repel = TRUE)

Variable feature plot showing mean expression versus standardized variance for all genes, with the top 2000 variable features highlighted in red and the top 10 labeled

In this plot, each point is a gene. The x-axis shows mean expression across cells and the y-axis shows standardized variance. The red points are the 2,000 selected variable features. Genes in the upper right are both highly expressed and highly variable – these are strong cell-type markers. Genes along the bottom have low variance relative to their expression level and are excluded.

Scaling

The final step before PCA is scaling. ScaleData() applies a z-score transformation to each gene across all cells: it subtracts the mean expression and divides by the standard deviation. After scaling, each gene has mean 0 and standard deviation 1.

Why is this necessary? Without scaling, PCA would be dominated by highly expressed genes simply because they have larger absolute values, not because they carry more biological information. Scaling puts all genes on equal footing so that PCA identifies the axes of true biological variation rather than expression magnitude.

Animation: ScaleData puts genes on one footing

R

pbmc <- ScaleData(pbmc)

By default, ScaleData() only scales the variable features (the 2,000 genes selected above), not all ~29,000 genes. This is faster and sufficient for PCA, which only uses variable features. If you later need scaled values for all genes (e.g., for a heatmap of a non-variable gene), you can re-run ScaleData() with features = rownames(pbmc).

After scaling, a third layer appears in the object:

R

pbmc

Expected output

OUTPUT

An object of class Seurat
29155 features across 11310 samples within 1 assay
Active assay: RNA (29155 features, 2000 variable features)
 3 layers present: counts, data, scale.data

All three layers are now populated: counts (raw), data (log-normalized), and scale.data (z-scored variable features).

Regressing out confounders

ScaleData() can optionally regress out unwanted sources of variation using the vars.to.regress parameter. A common use case is regressing out mitochondrial percentage so that cell damage does not influence clustering:

R

# Optional: regress out mitochondrial percentage
pbmc_regressed <- ScaleData(pbmc, vars.to.regress = "percent.mt")

This fits a linear model for each gene with percent.mt as a covariate and uses the residuals as the scaled expression values. The effect is that variation driven by mitochondrial content is removed before PCA.

When to regress:

Mitochondrial percentage is driving a major axis of variation in PCA (i.e., one of the top PCs strongly correlates with percent.mt)
You see cells grouping by mitochondrial content rather than cell type

When NOT to regress:

Mitochondrial content correlates with real biology. Some cell types (e.g., cardiomyocytes, metabolically active immune cells) naturally have higher mitochondrial gene expression. Regressing it out would remove legitimate biological signal.
The QC filtering in the previous episode already removed high-mt cells, and the remaining variation is modest.

For this workshop, we proceed without regression because our QC filtering already removed cells with percent.mt > 15, and the remaining variation does not dominate the PCA.

Save the normalized object for the next episode:

R

saveRDS(pbmc, file = "pbmc_normalized.rds")

SCTransform as an Alternative

Callout

SCTransform: a one-step alternative

The three-step LogNormalize workflow (NormalizeData + FindVariableFeatures + ScaleData) works well for most datasets, but it has limitations. The log transformation can over-stabilize variance for lowly expressed genes and under-stabilize it for highly expressed genes. This matters most when there is high technical noise in the data.

SCTransform addresses this by using regularized negative binomial regression. It models each gene’s expression as a function of total UMI counts per cell, estimates the expected variance at each expression level, and returns Pearson residuals that serve as normalized, variance-stabilized values. A single function call replaces all three steps:

R

pbmc_sct <- SCTransform(pbmc, vst.flavor = "v2")

Parameter	Value	Description
`vst.flavor`	`"v2"`	Use the improved v2 algorithm (Choudhary & Satija, 2022) with better variance estimation for sparse data.

Key differences from LogNormalize:

Replaces three steps with one: normalization, feature selection, and scaling are all handled internally.
Selects 3,000 variable features by default (vs. 2,000 for FindVariableFeatures), which can improve resolution for complex tissues.
Better variance stabilization across the full range of expression levels.
Stored in a separate assay: results go into an assay called SCT rather than modifying the RNA assay. The RNA assay retains the raw counts for reference.

When to use SCTransform:

For publication-quality analyses, especially with complex tissues or datasets with high technical variance.
When you observe that LogNormalize produces clusters driven by sequencing depth rather than biology.

When LogNormalize is sufficient:

For exploratory analysis and quick overviews.
When the dataset is well-behaved and you want a simpler, faster workflow.
For direct comparability with older analyses and tutorials.

In this workshop, we use LogNormalize throughout for simplicity and because it produces excellent results with this well-characterized PBMC dataset. However, SCTransform is the preferred method for most production analyses.

Challenge

Challenge 1: LogNormalize vs. SCTransform variable features

Run both normalization approaches on the same filtered dataset and compare the top 20 variable features from each method. How many genes appear in both lists?

This challenge is self-paced. SCTransform on the whole dataset needs about 15 GB of memory, so run it in an RStudio session with 12 or more cores (the workshop sessions have 8; start a new one for this). In the new session, reload the object you saved above:

R

library(Seurat)
work_dir <- paste0("/scratch/negishi/", Sys.getenv("USER"), "/scrna_workshop/")
pbmc <- readRDS(paste0(work_dir, "pbmc_normalized.rds"))

R

# LogNormalize pathway (already done above)
lognorm_top20 <- head(VariableFeatures(pbmc), 20)

# SCTransform pathway
pbmc_sct <- SCTransform(pbmc, vst.flavor = "v2")
sct_top20 <- head(VariableFeatures(pbmc_sct), 20)

# Compare
cat("LogNormalize top 20:\n")
print(lognorm_top20)
cat("\nSCTransform top 20:\n")
print(sct_top20)
cat("\nGenes in both lists:\n")
shared <- intersect(lognorm_top20, sct_top20)
print(shared)
cat("\nOverlap:", length(shared), "out of 20\n")

Show me the solution

OUTPUT

LogNormalize top 20:
 [1] "PPBP"            "PF4"             "LINC01478"       "LINC01374"       "PTGDS"
 [6] "JCHAIN"          "GP1BB"           "SOX5"            "ENSG00000225885" "EREG"
[11] "TCF4"            "IGKC"            "ZNF385D"         "GNLY"            "LYPD2"
[16] "CUX2"            "LINGO2"          "CAVIN2"          "CLNK"            "HTR1F"

SCTransform top 20:
 [1] "S100A9" "S100A8" "GNLY"   "IGKC"   "PPBP"   "IGLC2"  "LYZ"    "PF4"    "GP1BB"  "NRGN"   "IGLC3"
[12] "IGHM"   "CAVIN2" "AFF3"   "NKG7"   "GNG11"  "BANK1"  "VCAN"   "TUBB1"  "LINGO2"

Genes in both lists:
[1] "PPBP"   "PF4"    "GP1BB"  "IGKC"   "GNLY"   "LINGO2" "CAVIN2"

The overlap is 7 out of 20 (~35%), which is lower than you might expect. The shared genes are PPBP, PF4, GP1BB, IGKC, GNLY, LINGO2, and CAVIN2.

The two methods rank variable genes quite differently:

LogNormalize + vst picks up several lincRNAs (LINC01478, LINC01374) and even an Ensembl ID (ENSG00000225885) in the top 20. These are genes with sparse but extreme expression in a few cells, which inflates their residual variance in the vst model.
SCTransform surfaces well-known immune markers instead: S100A9, S100A8, LYZ, NKG7, VCAN. Its negative binomial model is more robust to outlier cells, so it favors genes with consistent cell-type-specific expression over genes that are simply noisy.

Despite only 35% overlap in the top 20, the broader variable feature sets (all 2,000-3,000 genes) agree much more. The top-ranked differences rarely change the downstream PCA, clustering, or cell type assignments because both methods capture the same major axes of biological variation. The practical difference is that SCTransform tends to produce cleaner ranked lists with fewer noise-driven genes at the top.

Callout

glmGamPoi speeds up SCTransform

If you see the warning about glmGamPoi not being installed, SCTransform still runs but falls back to a slower implementation. Install it for faster runs:

R

BiocManager::install("glmGamPoi")

Challenge

Challenge 2: Housekeeping vs. Variable Gene Behavior

Plot the expression of ACTB (a housekeeping gene) and LYZ (a variable monocyte marker) before and after normalization. What changes for each gene, and why?

R

# Extract raw and normalized values
actb_raw  <- as.numeric(pbmc[["RNA"]]$counts["ACTB", ])
actb_norm <- as.numeric(pbmc[["RNA"]]$data["ACTB", ])
lyz_raw   <- as.numeric(pbmc[["RNA"]]$counts["LYZ", ])
lyz_norm  <- as.numeric(pbmc[["RNA"]]$data["LYZ", ])

p1 <- ggplot(data.frame(x = actb_raw), aes(x)) +
    geom_histogram(bins = 50, fill = "steelblue") +
    labs(title = "ACTB - Raw counts") + theme_minimal()
p2 <- ggplot(data.frame(x = actb_norm), aes(x)) +
    geom_histogram(bins = 50, fill = "darkred") +
    labs(title = "ACTB - Normalized") + theme_minimal()
p3 <- ggplot(data.frame(x = lyz_raw), aes(x)) +
    geom_histogram(bins = 50, fill = "steelblue") +
    labs(title = "LYZ - Raw counts") + theme_minimal()
p4 <- ggplot(data.frame(x = lyz_norm), aes(x)) +
    geom_histogram(bins = 50, fill = "darkred") +
    labs(title = "LYZ - Normalized") + theme_minimal()

(p1 + p2) / (p3 + p4)

Show me the solution

Four histograms in a 2x2 grid. Top row: ACTB raw counts show a right-skewed distribution peaking near zero; ACTB normalized shows a unimodal bell shape centered around 3.5. Bottom row: LYZ raw counts show a massive zero spike with a flat tail to 500; LYZ normalized reveals a bimodal pattern with peaks near 0.5 and 4.5

ACTB (housekeeping gene):

Before normalization: the raw count distribution is right-skewed because cells with higher sequencing depth have proportionally more ACTB counts. The apparent variation is mostly technical (sequencing depth), not biological.
After normalization: the distribution becomes approximately bell-shaped, centered around 3.5-4 on the log-normalized scale. Because ACTB is expressed at similar levels across all cell types, removing depth effects reveals a tight, unimodal distribution. This is exactly why ACTB is not selected as a variable feature – it doesn’t help distinguish cell types.

LYZ (variable gene):

Before normalization: there is a massive spike at zero (~7,000 cells that don’t express LYZ) and a low, flat spread of non-zero values extending to ~500. The zero-dominated distribution makes it hard to see any structure in the expressing cells.
After normalization: the zero spike persists (you can’t normalize away true absence of expression), but the non-zero cells now reveal a clear bimodal pattern – a small cluster of cells with low expression (~1-2) and a distinct second peak of highly expressing cells (~4-5). This bimodality reflects the biological reality: CD14+ monocytes express LYZ at high levels, while other cell types either don’t express it at all or express it minimally. This separation is why LYZ is selected as a highly variable feature.

The key takeaway: normalization makes expression comparable across cells by removing sequencing depth effects. Housekeeping genes become unimodal (confirming they don’t vary biologically), while cell-type marker genes reveal biologically meaningful structure – like the bimodal split in LYZ – that was obscured by technical noise in the raw counts.

Key Points

Log normalization corrects for sequencing depth by dividing by total counts, scaling, and log-transforming each cell independently
Highly variable features are the ~2,000 genes with the most cell-to-cell variation, capturing biological differences while excluding housekeeping genes and noise
Scaling (z-score transformation) centers and standardizes gene expression so that PCA is not dominated by highly expressed genes
Regressing out confounders like percent.mt during scaling is optional and should only be done when technical variation is obscuring biology
SCTransform is a one-step alternative that uses regularized negative binomial regression and is recommended for publication-quality analyses

Content from Dimensionality Reduction and Clustering

Last updated on 2026-10-27 | Edit this page

Estimated time: 60 minutes

Overview

Questions

Why do we need to reduce the dimensions of scRNA-seq data before clustering?
How do we choose the right number of principal components?
What does UMAP show us and what can it NOT tell us?
How does graph-based clustering work and what does the resolution parameter control?
How do we evaluate whether our clustering resolution is appropriate?

Objectives

Run PCA on variable features and select an appropriate number of PCs using the elbow plot
Generate a UMAP embedding for 2D visualization of cell relationships
Construct a shared nearest neighbor graph and identify clusters with the Louvain algorithm
Explore multiple clustering resolutions and use clustree to assess cluster stability
Visualize known marker genes on the UMAP to preview cell type identity

Prerequisite

Prerequisites

This episode requires an RStudio session on the Negishi cluster. Launch RStudio (bioconductor) under Bioinformatics Apps on Open OnDemand as described in the Setup instructions. You will also need the pbmc_normalized.rds object saved at the end of the previous episode.

Setup

R

library(Seurat)
library(ggplot2)
library(patchwork)
library(clustree)

Loading the Normalized Data

We start from the normalized, scaled Seurat object saved at the end of the previous episode. This object has the counts, data, and scale.data layers populated, with 2,000 variable features selected.

Set up the working directory to match the previous episodes:

R

work_dir <- paste0(
    "/scratch/negishi/", Sys.getenv("USER"),
    "/scrna_workshop/"
)
setwd(work_dir)

Callout

Starting here without the previous episode?

Copy the saved object from the workshop data first: file.copy(paste0(work_dir, "rds/pbmc_normalized.rds"), work_dir)

R

pbmc <- readRDS(paste0(work_dir, "pbmc_normalized.rds"))
pbmc

Expected output

OUTPUT

An object of class Seurat
29155 features across 11310 samples within 1 assay
Active assay: RNA (29155 features, 2000 variable features)
 3 layers present: counts, data, scale.data

Principal Component Analysis (PCA)

Our dataset has ~2,000 variable features, but many of these genes are correlated – they go up and down together across cells because they are co-regulated in the same cell types. PCA (Principal Component Analysis) exploits these correlations to compress the data into a smaller set of principal components (PCs). Each PC is a linear combination of genes that captures a distinct axis of variation in the data.

The key properties of PCA:

PC1 captures the direction of greatest variance across all cells
PC2 captures the next greatest variance, orthogonal to PC1
Each subsequent PC captures decreasing amounts of variance
The first 10–20 PCs typically capture the major biological structure (cell types), while later PCs are dominated by noise

Animation: PCA compresses correlated genes

R

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

Parameter	Value	Description
`features`	`VariableFeatures(object = pbmc)`	Which genes to use for PCA. We use the 2,000 variable features selected in the previous episode.

By default, Seurat computes 50 PCs. Let’s examine which genes drive the first two components:

R

VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")

Dot plots of the top 30 genes by loading magnitude for PC1 and PC2, with PC1 dominated by myeloid markers and PC2 by B cell markers.

Each point represents a gene’s loading (contribution) to that PC. For PC1, you should see monocyte/myeloid markers (CST3, LYZ, S100A9, FCN1) among the top genes. This tells us that PC1 separates myeloid from lymphoid cells – the largest source of variation in PBMCs. PC2 is dominated by B cell markers (BANK1, CD79A, MS4A1, PAX5), indicating that B cells are the next major axis of variation.

We can also visualize PCs as heatmaps to see which genes and cells contribute most:

R

DimHeatmap(pbmc, dims = 1:9, cells = 500, balanced = TRUE)

Parameter	Value	Description
`dims`	`1:9`	Which PCs to plot.
`cells`	`500`	Number of cells to include (randomly sampled from each extreme of the PC).
`balanced`	`TRUE`	Plot equal numbers of cells from the positive and negative ends of each PC.

Each heatmap shows the top genes (rows) ordered by their loading, with cells (columns) ordered by their PC score. Clear blocks of correlated expression indicate PCs that capture real biological structure. PCs that look noisy (no clear pattern) are capturing technical variation or random noise.

Grid of nine heatmaps for PC1 through PC9, each showing top genes ordered by loading and cells ordered by PC score to reveal distinct expression programs.

Choosing the number of PCs

Not all 50 PCs are informative. Later PCs capture decreasing amounts of variance and eventually become dominated by noise. We need to decide how many PCs to retain for downstream steps (UMAP and clustering).

The elbow plot shows the standard deviation explained by each PC:

R

ElbowPlot(pbmc, ndims = 30)

Elbow plot of standard deviation versus PC number showing a steep decline through PC7 and a plateau after PC10, with the elbow around PC 8 to 10.

Look for the elbow – the point where the curve transitions from steep decline to a flatter plateau. PCs before the elbow capture substantial biological signal; PCs after the elbow are mostly noise.

Callout

How many PCs?

For the PBMC 10k dataset, the elbow falls around PC 8–10: the curve drops steeply through PC 7 and flattens after PC 10. We will use 15 PCs for the rest of this episode, a few more than the elbow (see below for why).

In practice, the exact number is not critical. Using 10 vs. 20 PCs usually produces very similar clustering results for well-separated cell types. The choice matters more for subtle distinctions between closely related populations.

When in doubt, err on the side of including more PCs (e.g., 20). Using too few PCs risks losing biological signal from rarer cell types. Using a few extra PCs adds a small amount of noise but rarely changes the major structure.

An alternative to visual inspection is the JackStraw procedure (JackStrawPlot()), which uses a permutation test to identify statistically significant PCs. However, it is computationally expensive for large datasets and is rarely necessary when the elbow is clear.

UMAP Visualization

PCA reduces our data from 2,000 dimensions to 15, but we still cannot visualize 15 dimensions directly. UMAP (Uniform Manifold Approximation and Projection) further reduces the data to just 2 dimensions for plotting. Unlike PCA, UMAP is a non-linear method: it can unfold complex relationships that PCA cannot capture, placing similar cells close together and dissimilar cells far apart in 2D space.

R

pbmc <- RunUMAP(pbmc, dims = 1:15)

Parameter	Value	Description
`dims`	`1:15`	Which PCs to use as input. We use the 15 PCs selected above.

Visualize the UMAP embedding:

R

DimPlot(pbmc, reduction = "umap")

UMAP scatter plot of 11310 cells in a single color before clustering, with cells forming several spatially distinct groups that correspond to different cell types.

At this point, cells are not yet clustered (all are the same color). But you can already see that cells organize into distinct groups on the UMAP. These groups correspond to different cell types.

UMAP parameters

UMAP has two key parameters that control the appearance of the plot:

Parameter	Default	Effect
`n.neighbors`	30	Controls the balance between local and global structure. Low values (5–15) emphasize local neighborhoods, producing tighter clusters. High values (50–200) emphasize global structure, producing smoother layouts.
`min.dist`	0.3	Controls how tightly UMAP packs points. Low values (0.01–0.1) produce dense, compact clusters. High values (0.5–1.0) spread clusters out more evenly.

The defaults work well for most datasets. Adjusting these parameters changes the visual appearance of the plot but does not change the underlying data or the clustering results (which operate on PCA space, not UMAP space).

Callout

UMAP interpretation pitfalls

UMAP plots are powerful visualization tools, but they can be misleading if over-interpreted. Three things you should NOT infer from a UMAP:

Distances between clusters do not reflect true biological similarity. Two clusters that appear far apart on the UMAP may actually be transcriptionally very similar, and vice versa. UMAP distorts global distances to preserve local neighborhoods.
Cluster size on the plot does not reflect population size. A large, spread-out cluster may contain fewer cells than a small, compact one. UMAP adjusts spacing to show local structure, not proportional representation.
Different random seeds produce different layouts. Running UMAP twice with different seeds will give a different arrangement of clusters. The overall groupings will be the same (the same cells will cluster together), but their positions and orientations on the plot will differ. Never compare spatial positions across two separate UMAP runs.

Always verify biological conclusions with quantitative methods (marker gene expression, differential expression tests, cluster composition statistics) rather than relying on visual impressions from UMAP alone.

Animation: What UMAP does not tell you

Callout

t-SNE as an alternative

Before UMAP became standard, t-SNE (t-distributed Stochastic Neighbor Embedding) was the most common visualization method for scRNA-seq. You can run it with RunTSNE(pbmc, dims = 1:15). t-SNE is similar in purpose to UMAP but tends to be slower and produces rounder, more separated clusters. The same interpretation caveats apply. Most current analyses use UMAP because it is faster, better preserves global structure, and scales better to large datasets.

Graph-Based Clustering

With the PCA embedding computed, we can now group cells into clusters. Seurat uses a graph-based clustering approach that works in two steps:

Step 1: Build a neighbor graph

FindNeighbors() constructs a K-nearest neighbor (KNN) graph in PCA space. For each cell, it identifies its k most similar cells based on Euclidean distance in the first 15 PCs. These connections are then refined into a shared nearest neighbor (SNN) graph, where the connection strength between two cells depends on how much their neighborhoods overlap.

R

pbmc <- FindNeighbors(pbmc, dims = 1:15)

Parameter	Value	Description
`dims`	`1:15`	Which PCs to use for computing distances. Must match the PCs used for UMAP.
`k.param`	`20` (default)	Number of nearest neighbors per cell. Higher values produce a more connected graph with smoother clustering. Lower values preserve finer local structure.

Step 2: Detect communities

FindClusters() applies a community detection algorithm to the SNN graph. The algorithm finds groups of cells that are densely connected to each other but sparsely connected to the rest of the graph. The resolution parameter controls the granularity of clustering:

R

pbmc <- FindClusters(pbmc, resolution = 0.5)

Parameter	Value	Description
`resolution`	`0.5`	Controls the number of clusters. Lower values (0.1–0.3) produce fewer, larger clusters. Higher values (1.0–2.0) produce more, smaller clusters.
`algorithm`	`1` (default)	Community detection method. `1` = Louvain (default, fast and reliable). `4` = Leiden (more robust against poorly connected clusters, requires the `leidenalg` Python package).

Animation: FindNeighbors and FindClusters

Expected output

At resolution 0.5, you should get approximately 18 clusters, numbered starting from 0. The exact number may vary slightly depending on random seed and software versions.

OUTPUT

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 11310
Number of edges: 410535

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9317
Number of communities: 18
Elapsed time: 1 seconds

The cluster assignments are stored in the object metadata. Let’s check them:

R

head(Idents(pbmc))
table(Idents(pbmc))

Expected output

OUTPUT

AAACCCAAGCGCCCAT-1 AAACCCAAGGTTCCGC-1 AAACCCACAGACAAGC-1 AAACCCACAGAGTTGG-1 AAACCCACAGGTATGG-1
                 0                 12                 13                  1                  7
AAACCCACATAGTCAC-1
                 4
Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17


   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17
2112 1810 1187 1015 1006  716  639  634  517  429  357  234  176  175  122   78   64   39

Cluster 0 is the largest with 2,112 cells (T cells, which are the most abundant PBMC type). Cluster 17 is the smallest with only 39 cells; the next episode shows that it is a rare population of haematopoietic stem and progenitor cells (HSPC).

Now let’s visualize the clusters on the UMAP:

R

DimPlot(pbmc, reduction = "umap", label = TRUE, label.size = 5) +
    NoLegend()

UMAP plot colored by Louvain cluster identity at resolution 0.5, showing 18 clusters numbered 0 through 17 with large clusters in the center and small clusters at the periphery.

Exploring multiple resolutions

The resolution parameter is the most important choice in clustering. Too low and distinct cell types are merged; too high and single cell types are split into artificial sub-clusters. The best way to evaluate the right resolution is to try several and compare.

R

pbmc <- FindClusters(pbmc, resolution = 0.2)
pbmc <- FindClusters(pbmc, resolution = 0.5)
pbmc <- FindClusters(pbmc, resolution = 0.8)
pbmc <- FindClusters(pbmc, resolution = 1.2)

Each call adds a new column to the metadata named RNA_snn_res.X (where X is the resolution value). We can visualize how cells move between clusters at different resolutions using clustree:

R

clustree(pbmc, prefix = "RNA_snn_res.")

Clustree diagram tracing cluster membership from resolution 0.2 to 1.2, with clean branches indicating stable splits and tangled arrows indicating unstable over-splitting.

In the clustree plot, each row is a resolution and each node is a cluster. Arrows show how cells move between clusters as resolution increases. Look for:

Clean splits: one cluster at a lower resolution cleanly divides into two at a higher resolution. This suggests a real biological distinction.
Unstable splits: cells from one cluster scatter across multiple clusters at a higher resolution, or cells bounce between clusters. This suggests over-splitting.

Set the active cluster identity to our chosen resolution:

R

Idents(pbmc) <- "RNA_snn_res.0.5"

Save the clustered object for the next episode:

R

saveRDS(pbmc, file = "pbmc_clustered.rds")

Exploring Clusters

Before formal cell type annotation (next episode), we can already preview cell identity by overlaying known marker genes on the UMAP. If our clustering is good, each marker should light up in a distinct cluster.

R

FeaturePlot(pbmc,
            features = c("CD3D", "MS4A1", "LYZ", "GNLY"),
            ncol = 2)

Four UMAP panels showing expression of CD3D, MS4A1, LYZ, and GNLY on a grey-to-purple gradient, each marking distinct cell populations for T cells, B cells, monocytes, and NK cells.

Gene	Cell type	Expected pattern
`CD3D`	T cells	Lights up the large T cell clusters (the biggest groups on the UMAP)
`MS4A1`	B cells	Lights up a compact cluster separate from T cells
`LYZ`	CD14+ Monocytes	Lights up monocyte clusters, distinct from lymphocytes
`GNLY`	NK cells	Lights up a cluster adjacent to but separate from T cells

We can also view these markers as violin plots to see the expression distribution across clusters:

R

VlnPlot(pbmc,
        features = c("CD3D", "MS4A1", "LYZ", "GNLY"),
        ncol = 2,
        pt.size = 0)

Violin plots of CD3D, MS4A1, LYZ, and GNLY expression across 18 clusters, showing each marker with high expression in specific clusters and near-zero expression elsewhere.

Even without formal annotation, you can see that these markers segregate cleanly into distinct clusters. CD3D marks several clusters (T cell subtypes), MS4A1 is specific to one cluster (B cells), LYZ is strong in one or two clusters (monocytes), and GNLY is specific to another (NK cells). This gives us confidence that the clustering is capturing biologically meaningful groupings. In the next episode, we will systematically identify marker genes for every cluster and assign cell type labels.

Challenge

Challenge 1: Choosing a Resolution with Clustree

Run clustering at resolutions 0.2, 0.5, 0.8, and 1.2 (we already did this above). Use the clustree plot to answer:

At which resolution do clusters start splitting unstably?
Which resolution gives a number of clusters that best matches the known PBMC cell types (~8–10 major types)?
How many clusters do you get at each resolution?

R

# Count clusters at each resolution
cat("Resolution 0.2:", length(unique(pbmc$RNA_snn_res.0.2)), "clusters\n")
cat("Resolution 0.5:", length(unique(pbmc$RNA_snn_res.0.5)), "clusters\n")
cat("Resolution 0.8:", length(unique(pbmc$RNA_snn_res.0.8)), "clusters\n")
cat("Resolution 1.2:", length(unique(pbmc$RNA_snn_res.1.2)), "clusters\n")

# Visualize side by side
p1 <- DimPlot(pbmc, group.by = "RNA_snn_res.0.2", label = TRUE) +
    NoLegend() + ggtitle("Resolution 0.2")
p2 <- DimPlot(pbmc, group.by = "RNA_snn_res.0.5", label = TRUE) +
    NoLegend() + ggtitle("Resolution 0.5")
p3 <- DimPlot(pbmc, group.by = "RNA_snn_res.0.8", label = TRUE) +
    NoLegend() + ggtitle("Resolution 0.8")
p4 <- DimPlot(pbmc, group.by = "RNA_snn_res.1.2", label = TRUE) +
    NoLegend() + ggtitle("Resolution 1.2")

(p1 + p2) / (p3 + p4)

Show me the solution

Side-by-side UMAP plots at resolutions 0.2, 0.5, 0.8, and 1.2 showing progressively more clusters as resolution increases.

Expected cluster counts:

Resolution	Approximate clusters
0.2	16
0.5	18
0.8	21
1.2	25

Resolution 0.5 gives 18 clusters, which captures the 8 major PBMC cell types (CD4+ T, CD8+ T, B cells, CD14+ monocytes, FCGR3A+ monocytes, NK cells, dendritic cells, platelets) along with some sub-clusters (e.g., naive vs. memory T cells), which is biologically reasonable.

Resolution 0.2 under-clusters with 16 clusters: it merges some distinct cell subtypes. You may lose meaningful biological distinctions.

Resolution 0.8 starts to over-split with 21 clusters: the clustree plot shows clusters from resolution 0.5 splitting into sub-clusters, with some cells moving between clusters in ways that are not stable. This does not mean 0.8 is wrong – it may reveal real sub-populations – but the splits become harder to annotate confidently.

Resolution 1.2 clearly over-splits with 25 clusters, many of which are fragments of the same cell type. The clustree plot shows extensive cross-cluster mixing at this resolution.

For this workshop we use resolution 0.5 as a good balance between capturing distinct cell types and avoiding excessive fragmentation.

Challenge

Challenge 2: Mapping Markers to Clusters

Use FeaturePlot() to visualize three markers on the UMAP:

CD3D (T cells)
MS4A1 (B cells)
LYZ (monocytes)

Then use the cluster-labeled UMAP to determine which cluster number(s) correspond to each cell type.

R

FeaturePlot(pbmc, features = c("CD3D", "MS4A1", "LYZ"), ncol = 3)
DimPlot(pbmc, reduction = "umap", label = TRUE, label.size = 5) + NoLegend()

Show me the solution

Three UMAP panels showing expression of CD3D, MS4A1, and LYZ. CD3D expression is concentrated in the upper-left and lower-left regions. MS4A1 lights up a compact group in the upper-right. LYZ is strongest in the lower-left clusters.

UMAP plot with 18 clusters labeled 0 through 17 for cross-referencing with the FeaturePlot above.

By comparing the FeaturePlot panels with the labeled UMAP, you should see the following pattern (exact cluster numbers may differ depending on your run):

Marker	Cell type	Cluster(s)
CD3D	T cells	Clusters 0, 3, 5, 6 and 10 (the large T cell group and its branches). CD3D is expressed in both CD4+ and CD8+ T cells, so it marks multiple clusters.
MS4A1	B cells	Clusters 4, 8 and 14 (one separate group). MS4A1 (also known as CD20) is highly specific to B cells.
LYZ	CD14+ Monocytes	Clusters 1 and 2, with the smaller clusters 11 and 16 (the large monocyte group). LYZ is also high in the dendritic cells (cluster 12) and weaker in the FCGR3A+ monocytes (cluster 9).

The fact that each marker cleanly maps to specific clusters confirms that our clustering at resolution 0.5 is capturing meaningful biological groupings. Notice that CD3D marks multiple clusters – this is expected because there are several T cell subtypes (CD4+ naive, CD4+ memory, CD8+) that form separate clusters but all express the pan-T cell marker CD3D.

In the next episode, we will use FindAllMarkers() to systematically identify the genes that define each cluster and assign proper cell type labels.

Key Points

PCA reduces ~2,000 variable genes to ~15 principal components that capture the major axes of biological variation
The elbow plot helps choose how many PCs to retain; for PBMCs, 15 PCs is typically sufficient
UMAP provides an intuitive 2D visualization but should not be used to infer distances, cluster sizes, or quantitative relationships
Graph-based clustering builds a shared nearest neighbor graph in PCA space and detects communities using the Louvain or Leiden algorithm
The resolution parameter controls clustering granularity; clustree helps identify the resolution where clusters begin to split unstably

Content from Cell Type Annotation

Last updated on 2026-10-27 | Edit this page

Estimated time: 60 minutes

Overview

Questions

How do we systematically find the genes that define each cluster?
How do we assign cell type labels using known marker genes?
How does automated annotation with SingleR work?
When do manual and automated annotations disagree and what should we do?
What are best practices for confident, reproducible annotation?

Objectives

Run FindAllMarkers to identify differentially expressed genes for every cluster
Interpret marker gene statistics: avg_log2FC, pct.1, pct.2, and p_val_adj
Manually annotate clusters using known PBMC marker genes
Use SingleR with the Monaco Immune reference for automated cell type annotation
Compare manual and automated annotations and resolve discrepancies

Prerequisite

Prerequisites

This episode requires an RStudio session on the Negishi cluster. Launch RStudio (bioconductor) under Bioinformatics Apps on Open OnDemand as described in the Setup instructions. You will also need the pbmc_clustered.rds object saved at the end of the previous episode.

Setup

R

library(Seurat)
library(ggplot2)
library(patchwork)
library(SingleR)
library(celldex)
library(dplyr)

Loading the Clustered Data

We start from the clustered Seurat object saved at the end of the previous episode. Make sure the active identity is set to the resolution 0.5 clusters.

Set up the working directory to match the previous episodes:

R

work_dir <- paste0(
    "/scratch/negishi/", Sys.getenv("USER"),
    "/scrna_workshop/"
)
setwd(work_dir)

Callout

Starting here without the previous episode?

Copy the saved object from the workshop data first: file.copy(paste0(work_dir, "rds/pbmc_clustered.rds"), work_dir)

R

pbmc <- readRDS(paste0(work_dir, "pbmc_clustered.rds"))
Idents(pbmc) <- "RNA_snn_res.0.5"
DimPlot(pbmc, reduction = "umap", label = TRUE, label.size = 5) + NoLegend()

Expected output

UMAP plot with 18 clusters labeled 0 through 17 at resolution 0.5, matching the clustering results from the previous episode.

You should see the UMAP with 18 clusters labeled by number (0 to 17), matching the clustering results from the previous episode.

Finding Marker Genes

Up to now, our clusters are just numbers. To give them biological meaning, we need to find the genes that distinguish each cluster from all other clusters. These are called marker genes – genes that are significantly upregulated in one cluster compared to the rest.

FindAllMarkers() runs a differential expression test for every cluster, comparing cells in that cluster against all other cells:

R

pbmc.markers <- FindAllMarkers(pbmc,
                               only.pos = TRUE,
                               min.pct = 0.25,
                               logfc.threshold = 0.25)

Parameter	Value	Description
`only.pos`	`TRUE`	Only return genes that are upregulated in the cluster (positive log fold change). We are interested in markers that define a cluster, not genes that are absent from it.
`min.pct`	`0.25`	Only test genes detected in at least 25% of cells in either the cluster or the rest. This speeds up the computation by skipping very rare genes that cannot be reliable markers.
`logfc.threshold`	`0.25`	Only test genes with at least a 0.25 log2 fold change between the cluster and the rest. This pre-filters out genes with trivially small differences.

This step may take 5 minutes.

Expected output

OUTPUT

Calculating cluster 0
Calculating cluster 1
Calculating cluster 2
Calculating cluster 3
Calculating cluster 4
Calculating cluster 5
Calculating cluster 6
Calculating cluster 7
Calculating cluster 8
Calculating cluster 9
Calculating cluster 10
Calculating cluster 11
Calculating cluster 12
Calculating cluster 13
Calculating cluster 14
Calculating cluster 15
Calculating cluster 16
Calculating cluster 17

The result is a data frame with one row per gene per cluster. Let’s look at the key columns:

R

head(pbmc.markers)

Expected output

OUTPUT

       p_val avg_log2FC pct.1 pct.2 p_val_adj cluster   gene
INPP4B     0   2.607036 0.972 0.345         0       0 INPP4B
CAMK4      0   1.653160 0.950 0.324         0       0  CAMK4
IL7R       0   2.067086 0.935 0.336         0       0   IL7R
CD3D       0   1.523352 0.923 0.334         0       0   CD3D
ITK        0   1.775408 0.920 0.338         0       0    ITK
BCL11B     0   1.586310 0.991 0.413         0       0 BCL11B

The columns mean:

Column	Meaning
`p_val`	Raw p-value from the Wilcoxon rank-sum test (default test).
`avg_log2FC`	Average log2 fold change between the cluster and all other cells. A value of 1.0 means the gene is on average 2x higher in this cluster.
`pct.1`	Fraction of cells in the cluster where the gene is detected.
`pct.2`	Fraction of cells outside the cluster where the gene is detected.
`p_val_adj`	Bonferroni-adjusted p-value. Use this for significance filtering (typically < 0.05).
`cluster`	Which cluster this gene is a marker for.
`gene`	The gene symbol.

A good marker gene has a high avg_log2FC, high pct.1 (expressed in most cells of the cluster), low pct.2 (not expressed in other clusters), and a significant p_val_adj.

Animation: How FindAllMarkers scores a gene

Top markers per cluster

Let’s extract the top 5 markers per cluster ranked by fold change:

R

top5 <- pbmc.markers %>%
    dplyr::group_by(cluster) %>%
    dplyr::slice_max(n = 5, order_by = avg_log2FC)
top5

Expected output

OUTPUT

# A tibble: 90 × 7
# Groups:   cluster [18]
   p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
   <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>
 1     0       5.23 0.258 0.011         0 0       IATPR
 2     0       3.69 0.345 0.036         0 0       ST8SIA1
 3     0       3.17 0.346 0.052         0 0       IL2RA
 4     0       2.88 0.414 0.072         0 0       ICOS
 5     0       2.75 0.521 0.103         0 0       FAAH2
 6     0       4.46 0.29  0.023         0 1       ENSG00000289381
 7     0       4.28 0.5   0.05          0 1       AQP9
 8     0       4.13 0.423 0.043         0 1       MTARC1
 9     0       4.09 1     0.549         0 1       S100A8
10     0       3.98 0.996 0.248         0 1       S100A12

Visualizing markers with a heatmap

A heatmap of top marker genes across all clusters shows the expression patterns at a glance. Each column is a cell (grouped by cluster), each row is a gene, and color intensity represents scaled expression:

R

top3 <- pbmc.markers %>%
    dplyr::group_by(cluster) %>%
    dplyr::slice_max(n = 3, order_by = avg_log2FC)

DoHeatmap(pbmc, features = top3$gene) + NoLegend()

Heatmap showing expression of the top 3 marker genes per cluster across all 18 clusters. Each cluster shows a distinct block of upregulated genes, confirming distinct transcriptional identities.

Each cluster should show a distinct block of highly expressed genes (bright yellow/white) that are low in other clusters (dark purple). Clusters with very similar heatmap profiles may represent subtypes of the same cell type.

Dot plot summary

A dot plot provides a compact summary of marker expression across clusters. The dot size encodes the percentage of cells expressing the gene, and the color intensity encodes the average expression level:

R

DotPlot(pbmc, features = unique(top3$gene)) + coord_flip() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

Dot plot showing expression of top 3 marker genes across all 18 clusters. Dot size indicates the percentage of cells expressing each gene; dot shading indicates average expression level.

Manual Annotation with Known Markers

Now we match the marker gene patterns to known PBMC cell types. The table below lists the canonical markers for each expected cell type:

Cell type	Key markers	Expected %
CD4+ T cells	CD3D, IL7R, CCR7	~35%
CD8+ T cells	CD3D, CD8A, CD8B	~10%
B cells	MS4A1, CD79A	~10%
CD14+ Monocytes	LYZ, CD14, S100A9	~20%
FCGR3A+ Monocytes	FCGR3A, MS4A7	~5%
NK cells	GNLY, NKG7, KLRD1	~10%
Dendritic cells	FCER1A, CST3	~3%
Plasmacytoid DCs (pDC)	LILRA4, IL3RA	~0.5%
Platelets	PPBP, PF4	~2%
Haematopoietic progenitors (HSPC)	CD34	<0.5%

Let’s visualize these markers on the UMAP to see which clusters they map to:

R

FeaturePlot(pbmc,
            features = c("CD3D", "IL7R", "CD8A", "MS4A1",
                         "LYZ", "FCGR3A", "GNLY", "FCER1A"),
            ncol = 4)

Grid of 8 UMAP panels showing expression of CD3D, IL7R, CD8A, MS4A1, LYZ, FCGR3A, GNLY, and FCER1A. Each marker highlights a distinct UMAP region: CD3D and IL7R in the upper-left T cell area, CD8A in a small cluster on the right, MS4A1 in B cells on the right, LYZ in monocytes at the lower-left, FCGR3A in a small cluster at the lower-right, GNLY in NK cells at the upper-center, and FCER1A in a small dendritic cell group.

And as violin plots to see expression across clusters:

R

VlnPlot(pbmc,
        features = c("CD3D", "IL7R", "CD8A", "MS4A1",
                     "LYZ", "FCGR3A", "GNLY", "PPBP"),
        ncol = 4,
        pt.size = 0)

Violin plots of CD3D, IL7R, CD8A, MS4A1, LYZ, FCGR3A, GNLY, and PPBP across all 18 clusters. Each marker shows high expression in the clusters corresponding to its known cell type and low expression elsewhere.

Based on these plots and the FindAllMarkers results, we can build a mapping from cluster numbers to cell type names: find which clusters express each set of markers.

Callout

Check the mapping against your own clusters

The labels below fit the 18 clusters that this lesson’s matrix gives with the package versions in the setup instructions. If your matrix or your package versions differ, the number or the order of the clusters can change, and a label vector written for someone else’s clusters puts the wrong names on yours without any error. Before trusting the labels, plot the canonical markers per cluster and compare each row with the table above and the comment next to each label below:

R

DotPlot(pbmc, features = c("CD3D", "IL7R", "CCR7", "CD8A", "CD8B", "MS4A1",
                           "CD79A", "LYZ", "CD14", "FCGR3A", "MS4A7", "GNLY",
                           "NKG7", "FCER1A", "CLEC10A", "LILRA4", "IL3RA",
                           "PPBP", "PF4", "CD34")) + coord_flip()

The code below also stops if your number of clusters does not match the number of labels.

Animation: Labels are assigned by position

R

# Build the cluster-to-cell-type mapping
# Adjust cluster IDs to match YOUR results
new.cluster.ids <- c(
    "CD4 T",           # 0: CD3D+, IL7R+, CD8A-
    "CD14 Mono",       # 1: LYZ+, CD14+, S100A9+
    "CD14 Mono",       # 2: LYZ+, CD14+, S100A9+
    "CD4 T",           # 3: CD3D+, IL7R+, CCR7+ (naive)
    "B",               # 4: MS4A1+, CD79A+
    "CD8 T",           # 5: CD3D+, CD8A+, CD8B+, NKG7+ (cytotoxic)
    "CD8 T",           # 6: CD3D+, CD8A+, IL7R+, NKG7+, CD8B low
    "NK",              # 7: GNLY+, NKG7+, KLRD1+, CD3D-
    "B",               # 8: MS4A1+, CD79A+
    "FCGR3A+ Mono",    # 9: FCGR3A+, MS4A7+
    "CD8 T",           # 10: CD3D+, CD8A+, CD8B+, CCR7+ (naive)
    "CD14 Mono",       # 11: LYZ+, CD14+
    "DC",              # 12: FCER1A+, CLEC10A+, CST3+
    "Platelet",        # 13: PPBP+, PF4+
    "B",               # 14: MS4A1+, CD79A+, some T cell signal
    "pDC",             # 15: LILRA4+, IL3RA+
    "CD14 Mono",       # 16: LYZ+, CD14+, some MS4A1
    "HSPC"             # 17: CD34+, CD3D-
)
if (nlevels(pbmc) != length(new.cluster.ids)) stop("pbmc has ", nlevels(pbmc), " clusters but new.cluster.ids has ", length(new.cluster.ids), " labels; check the markers and edit new.cluster.ids")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)

Callout

Annotation is iterative

Your first-pass annotation will not always be correct. This is normal. After assigning labels, go back and check:

Do the labels make biological sense? Is the proportion of each cell type reasonable for PBMCs?
Are there clusters where multiple marker sets overlap? This may indicate a mixed cluster that should be sub-clustered at higher resolution.
Are there clusters with no clear markers? These may be low-quality cells that slipped through QC, or a rare cell type you didn’t expect.

Annotation is an iterative process: assign labels, verify with markers, revise, and repeat until you are confident in the assignments.

Visualize the final annotated UMAP:

R

DimPlot(pbmc, reduction = "umap", label = TRUE, repel = TRUE) + NoLegend()

UMAP plot with cells labeled by manually assigned cell type names: CD4 T, CD14 Mono, B, CD8 T, NK, FCGR3A+ Mono, DC, Platelet, pDC, and HSPC.

Store the annotation in the metadata so it persists after saving:

R

pbmc$manual_annotation <- Idents(pbmc)

Automated Annotation with SingleR

Manual annotation requires expert knowledge of marker genes for every expected cell type. SingleR offers an automated alternative: it compares the expression profile of each cell to a labeled reference dataset and assigns the label of the most similar reference sample.

Animation: How SingleR labels a cell

Download the reference

We use the Monaco Immune Data reference from the celldex package. This reference contains bulk RNA-seq profiles of sorted human immune cell populations, making it well-suited for PBMC annotation.

R

ref <- celldex::MonacoImmuneData()
ref

Expected output

OUTPUT

class: SummarizedExperiment
dim: 46077 114
metadata(0):
assays(1): logcounts
rownames(46077): A1BG A1BG-AS1 ... ZYX ZZEF1
rowData names(0):
colnames(114): DZQV_CD8_naive DZQV_CD8_CM ... G4YW_Neutrophils G4YW_Basophils
colData names(3): label.main label.fine label.ont

The reference has 114 samples across multiple immune cell types. The label.main column has broad cell type labels (e.g., “Monocytes”, “CD4+ T cells”) that match well with our expected PBMC types.

Run SingleR

SingleR works on a SingleCellExperiment object. We convert our Seurat object and then run the annotation:

R

# Convert to SingleCellExperiment
sce <- as.SingleCellExperiment(DietSeurat(pbmc))

# Run SingleR with Monaco Immune reference
results <- SingleR(test = sce,
                   ref = ref,
                   labels = ref$label.main)

Parameter	Value	Description
`test`	`sce`	The query dataset as a SingleCellExperiment. SingleR uses the log-normalized expression values.
`ref`	`ref`	The reference dataset with known labels.
`labels`	`ref$label.main`	Which column of the reference metadata to use as cell type labels. `label.main` gives broad categories; `label.fine` gives more specific subtypes.

This runs for 1–2 minutes. SingleR assigns a label to every individual cell, not to clusters. Let’s transfer the labels to our Seurat object:

R

pbmc$singler_labels <- results$labels

Visualize the SingleR annotations on the UMAP:

R

DimPlot(pbmc, group.by = "singler_labels", reduction = "umap",
        label = TRUE, repel = TRUE) + NoLegend()

UMAP plot with cells colored by SingleR automated annotation labels from the Monaco Immune reference, showing CD4+ T cells, CD8+ T cells, T cells, Monocytes, B cells, NK cells, Dendritic cells, Progenitors, and Basophils.

Comparing manual and automated annotations

The real power of annotation comes from comparing multiple approaches. Let’s create a cross-tabulation of our manual labels versus SingleR labels:

R

table(Manual = pbmc$manual_annotation, SingleR = pbmc$singler_labels)

Expected output

OUTPUT

              SingleR
Manual         B cells Basophils CD4+ T cells CD8+ T cells Dendritic cells Monocytes NK cells Progenitors
  CD4 T              3         0         2586          183               0         5        1           1
  CD14 Mono          3         0            0            0              32      3259        1           0
  B               1614         0            0            0               5         0       15           0
  CD8 T              0         0            9          400               0         0       31           0
  NK                 0         0            3            1               0         0      607           0
  FCGR3A+ Mono       0         0            0            0               1       428        0           0
  DC                 0         0            0            0             157        19        0           0
  Platelet           4         2            3            2               1         6        0         154
  pDC                0         0            0            0              78         0        0           0
  HSPC               0         0            1            0               0         0       11          21
              SingleR
Manual         T cells
  CD4 T            348
  CD14 Mono          0
  B                 11
  CD8 T           1272
  NK                23
  FCGR3A+ Mono       0
  DC                 0
  Platelet           3
  pDC                0
  HSPC               6

The table shows broad agreement, and the disagreements that remain are instructive:

Strong agreement: SingleR calls 3,259 of 3,295 CD14 Monocytes “Monocytes”, 1,614 of 1,645 B cells “B cells”, 607 of 634 NK cells “NK cells”, 428 of 429 FCGR3A+ Monocytes “Monocytes”, and 157 of 176 DCs “Dendritic cells”. These cell types have distinctive transcriptional profiles that both methods recognize.
CD4 T cells: 2,586 of 3,127 are “CD4+ T cells”. Most of the rest (348) get the generic “T cells” label.
CD8 T cells: only 400 of 1,712 are “CD8+ T cells”; 1,272 are generic “T cells”. Nearly all of those come from cluster 6, which expresses CD8A, IL7R and NKG7 but little CD8B. That pattern fits mucosal-associated invariant T (MAIT) cells, which the Monaco reference groups with gamma-delta T cells under its generic “T cells” label. Our coarse “CD8 T” label is defensible; sub-clustering would separate them.
Platelets called “Progenitors”: 154 of our 175 platelets (PPBP+, PF4+) are called “Progenitors”. Monaco Immune Data has no platelets, so SingleR assigns the closest profile it has. Always check rare or non-standard cell types with markers.
pDCs called “Dendritic cells”: all 78 pDCs (LILRA4+, IL3RA+) get the coarser “Dendritic cells” label, which is correct at this level; Monaco’s fine labels separate plasmacytoid from myeloid DCs.
HSPC: the 39 CD34+ cells of cluster 17 are mostly “Progenitors” (21), with some “NK cells” (11). A cluster this small is often a mixture; it is labelled for completeness.

The key takeaway: agreement on the large populations gives confidence in both methods, and every disagreement has a reason you can check with marker genes. A large disagreement (for example, a whole cluster you labelled as T cells called “B cells” by SingleR) is a signal to go back to the markers for that cluster before trusting either label. Use both methods and investigate every disagreement.

Visualize both annotations side by side:

R

p1 <- DimPlot(pbmc, group.by = "manual_annotation", label = TRUE, repel = TRUE) +
    NoLegend() + ggtitle("Manual annotation")
p2 <- DimPlot(pbmc, group.by = "singler_labels", label = TRUE, repel = TRUE) +
    NoLegend() + ggtitle("SingleR annotation")
p1 + p2

Two UMAP plots side by side. Left panel shows manual annotation with CD4 T, CD14 Mono, B, CD8 T, NK, FCGR3A+ Mono, DC, Platelet, pDC, and HSPC labels. Right panel shows SingleR annotation with CD4+ T cells, CD8+ T cells, T cells, Monocytes, B cells, NK cells, Dendritic cells, Progenitors, and Basophils. Most cell groups receive consistent labels between the two methods.

When the two methods agree, you can be very confident in the label. When they disagree, check the marker gene expression for the disputed cluster and decide which label fits better. In general:

Trust manual annotation when the markers are clear and well-established for the tissue type you are studying
Trust SingleR when you are unfamiliar with the tissue or when the reference provides finer subtypes than you would identify manually

Save the final annotated object:

R

saveRDS(pbmc, file = "pbmc_annotated.rds")

Annotation Best Practices

Callout

Guidelines for confident annotation

Use multiple references. No single reference is perfect. The celldex package provides several options for human immune cells:
- MonacoImmuneData() – sorted immune populations, best for PBMCs
- HumanPrimaryCellAtlasData() (HPCA) – broader coverage including non-immune cell types, but less specific for immune subtypes
- DatabaseImmuneCellExpressionData() (DICE) – another immune-focused reference with different experimental conditions
Running SingleR with multiple references and checking for agreement increases confidence.
Cross-validate automated labels with marker genes. Never accept automated labels blindly. Always check that the assigned label is consistent with the expression of established markers for that cell type.
Document your decisions. Record which markers you used, which reference datasets, and any ambiguous cases you resolved manually. This makes your analysis reproducible and allows others to evaluate your annotation choices.
Sub-cluster when needed. If a cluster shows mixed marker expression (e.g., both CD4 and CD8 markers), it may contain two cell types that were merged at the current resolution. Re-run FindClusters() at higher resolution on just that subset of cells to resolve the mixture:
R
```
subset_cells <- subset(pbmc, idents = "mixed_cluster")
subset_cells <- FindNeighbors(subset_cells, dims = 1:15)
subset_cells <- FindClusters(subset_cells, resolution = 0.5)
```

Challenge

Challenge 1: Identifying T cell subtypes from marker combinations

Plot five T cell markers across the annotated groups. In the CD4 T group, CD3D and IL7R are expressed almost everywhere, but SELL and CCR7 only in part of it. What does that tell you about the composition of this group? Then look at CD69, an early activation marker: which groups express it, and can it pick out activated T cells in this dataset?

R

VlnPlot(pbmc, features = c("CD3D", "SELL", "CCR7", "CD69", "IL7R"), ncol = 5, pt.size = 0)

Show me the solution

Violin plots of CD3D, SELL, CCR7, CD69, and IL7R across the ten annotated cell types. CD3D and IL7R are high in CD4 T and CD8 T cells and near zero elsewhere. SELL is expressed in most groups. CCR7 is expressed mainly in CD4 T cells. CD69 is high in CD4 T, CD8 T, B, NK and HSPC cells alike.

Percent of cells in each group with any expression of the marker:

Marker	CD4 T	CD8 T	B	NK	Interpretation
CD3D	92%	91%	9%	13%	Confirms T cell identity
IL7R (CD127)	93%	85%	11%	13%	Naive and memory T cells
SELL (CD62L)	69%	41%	80%	45%	Lymph node homing; naive and central memory T cells
CCR7	55%	26%	34%	6%	Lymph node homing; high in naive, low in effector memory T cells
CD69	88%	87%	90%	87%	Early activation marker

The CD4 T group is a mixture. Nearly every cell is CD3D+ and IL7R+, so the group is T cells throughout. But only about two thirds express SELL and about half express CCR7. The SELL+ CCR7+ cells are naive and central memory T cells, which recirculate through lymph nodes; the SELL- CCR7- cells are effector memory T cells. Most of the CCR7 signal comes from cluster 3, one of the two clusters we merged into “CD4 T”.

CD8 T cells are more differentiated. SELL (41%) and CCR7 (26%) are lower than in CD4 T cells, so the CD8 T group holds more effector and effector memory cells, including the MAIT-like cluster 6.

CD69 does not separate activated cells here. It is expressed in almost 90% of CD4 T, CD8 T, B and NK cells alike. CD69 is an immediate early gene that can also be switched on during sample collection and processing, so in a dataset like this one it does not mark a recently activated subset. Check whether an “activation” marker is specific before interpreting it.

SELL alone does not mark T cells. B cells (80%), pDCs and HSPCs express it too; you need CD3D to confirm T cell identity.

This illustrates why our coarse “CD4 T” label is a simplification. With higher clustering resolution or sub-clustering, you could separate the naive and memory subsets. For many analyses this level of granularity is sufficient, but for studies focused on T cell biology you would want to sub-cluster further.

Challenge

Challenge 2: Comparing two SingleR references

Run SingleR with both MonacoImmuneData and HumanPrimaryCellAtlasData. For which cell types do the two references disagree? Why might they disagree?

R

# Monaco Immune reference (already computed above)
ref_monaco <- celldex::MonacoImmuneData()
results_monaco <- SingleR(test = sce, ref = ref_monaco, labels = ref_monaco$label.main)

# Human Primary Cell Atlas reference
ref_hpca <- celldex::HumanPrimaryCellAtlasData()
results_hpca <- SingleR(test = sce, ref = ref_hpca, labels = ref_hpca$label.main)

# Transfer labels
pbmc$singler_monaco <- results_monaco$labels
pbmc$singler_hpca <- results_hpca$labels

# Compare
table(Monaco = pbmc$singler_monaco, HPCA = pbmc$singler_hpca)

Show me the solution

OUTPUT

                 HPCA
Monaco            B_cell   BM  CMP  GMP HSC_-G-CSF  MEP Monocyte NK_cell Platelets Pre-B_cell_CD34-
  B cells           1575    0    0    0          1    0        3       0         2                1
  Basophils            0    0    0    0          0    0        0       0         2                0
  CD4+ T cells         0    0    0    0          0    0        0       0         0                0
  CD8+ T cells         0    0    0    0          0    0        0       0         1                0
  Dendritic cells     36    0    0    8          1    0      218       0         1                6
  Monocytes            4    0    0    0          7    0     3669       5         5                8
  NK cells             2    0    0    0          0    0        0     649         0                1
  Progenitors          0    1   20    0          0    1        3       0       151                0
  T cells              3    0    0    0          0    0        0      20         1                0
                 HPCA
Monaco            Pro-Myelocyte T_cells
  B cells                     1      41
  Basophils                   0       0
  CD4+ T cells                0    2602
  CD8+ T cells                0     585
  Dendritic cells             0       4
  Monocytes                   0      19
  NK cells                    0      14
  Progenitors                 0       0
  T cells                     0    1639

The two references agree well on some cell types and disagree sharply on others. Here are the key patterns:

1. Strong agreement: Monocytes, B cells, NK cells. Monaco “Monocytes” map almost entirely to HPCA “Monocyte” (3,669 of 3,717). Monaco “B cells” map to HPCA “B_cell” (1,575 of 1,624). Monaco “NK cells” map to HPCA “NK_cell” (649 of 666). These cell types have distinctive profiles that both references recognize.

2. T cell subtyping collapses in HPCA. Monaco distinguishes “CD4+ T cells” (2,602), “CD8+ T cells” (586), and generic “T cells” (1,663). HPCA labels all of these simply as “T_cells.” If your analysis requires T cell subtype resolution, Monaco is the better reference. Interestingly, Monaco assigns 41 B cells to HPCA “T_cells,” suggesting a small population that sits near the T/B boundary (possibly NKT cells or doublets).

3. Dendritic cells split across HPCA labels. Monaco calls 274 cells “Dendritic cells,” but HPCA scatters these across “Monocyte” (218), “GMP” (8), “Pre-B_cell_CD34-” (6), “B_cell” (36), and only 4 as “T_cells.” HPCA lacks a dedicated dendritic cell category with sufficient resolution, so DCs get assigned to the nearest available profile – usually monocytes, which share many myeloid markers.

4. Progenitors become Platelets in HPCA. Monaco labels 176 cells as “Progenitors,” and HPCA calls 151 of those “Platelets.” Most of them (154) are the cells we manually annotated as platelets; another 21 are our HSPC cluster. Monaco does not have a platelet reference profile, so it assigns platelets to “Progenitors”, the nearest match it has. HPCA does include a platelet profile and gets this right. Neither reference is wrong – they just have different coverage.

5. HPCA introduces hematopoietic progenitor labels. HPCA assigns some cells to “BM,” “CMP,” “GMP,” “MEP,” “HSC_-G-CSF,” and “Pro-Myelocyte” – categories that are rare in blood (our 39-cell HSPC cluster is the genuine case). Most of these calls are artifacts of the HPCA reference containing bone marrow profiles. When a cell’s expression sits between two immune types, HPCA may assign a progenitor label rather than committing to either.

Lesson: No single reference is universally best. Monaco gives better immune subtype resolution (especially for T cells and DCs) and avoids spurious non-immune labels, making it the better default for PBMCs. HPCA correctly identifies platelets where Monaco fails. The practical approach: use Monaco as your primary reference for immune tissues, cross-check disagreements with HPCA, and always validate automated labels against known marker genes.

Key Points

FindAllMarkers identifies genes upregulated in each cluster versus all others; key output columns are avg_log2FC, pct.1, pct.2, and p_val_adj
Manual annotation maps cluster marker genes to known cell-type-specific markers and requires domain expertise
SingleR automates annotation by comparing each cell’s expression profile to a labeled reference dataset
Combining manual and automated approaches increases confidence; agreement between methods strongly supports the assigned label
Annotation is iterative and should be cross-validated with multiple references and marker gene visualization

Content from Multi-Sample Integration

Last updated on 2026-10-27 | Edit this page

Estimated time: 60 minutes

Overview

Questions

Why do cells from different samples cluster separately even when they are the same cell type?
How does CCA-based integration correct for batch effects?
How do we evaluate whether integration was successful without removing real biology?
How do we find genes that differ between conditions within a specific cell type?
When should we NOT integrate datasets?

Objectives

Explain the batch effect problem and why simple merging is insufficient
Integrate two conditions using Seurat v5 CCA integration with IntegrateLayers
Evaluate integration by checking condition mixing and preserved marker expression
Identify conserved markers and condition-specific differentially expressed genes
Recognize when integration is inappropriate for the experimental design

Prerequisite

Prerequisites

This episode requires an RStudio session on the Negishi cluster. Launch RStudio (bioconductor) under Bioinformatics Apps on Open OnDemand as described in the Setup instructions.

Setup

R

library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)

Set up the working directory to match the previous episodes:

R

work_dir <- paste0(
    "/scratch/negishi/", Sys.getenv("USER"),
    "/scrna_workshop/"
)
setwd(work_dir)

Why Integration Is Needed

In the previous episodes, we analyzed a single sample of PBMCs. But many experimental designs involve multiple samples – different patients, different time points, or different conditions (treated vs. control). When we combine these samples for joint analysis, a common problem emerges: cells cluster by which sample they came from rather than by what cell type they are.

This happens because of batch effects: technical differences between samples caused by variation in cell handling, library preparation, sequencing depth, or reagent lots. These technical differences can be large enough to dominate over the biological differences between cell types, making it impossible to compare the same cell type across conditions.

Integration algorithms correct for batch effects by finding shared biological structure across samples. The goal is to align the same cell types so they co-cluster, while preserving the real biological differences (like changes in gene expression due to treatment).

Let’s see this problem in practice using a dataset of PBMCs from control and IFN-beta stimulated conditions.

Loading the IFNB Dataset

The IFNB dataset from the SeuratData package contains PBMCs from two conditions:

CTRL: untreated control PBMCs
STIM: PBMCs stimulated with interferon-beta (IFN-beta), a cytokine that activates antiviral immune responses

R

InstallData("ifnb")
ifnb <- LoadData("ifnb")
ifnb

Expected output

An object of class Seurat
14053 features across 13999 samples within 1 assay
Active assay: RNA (14053 features, 0 variable features)
 2 layers present: counts, data

The dataset has approximately 14,000 cells and 14,000 genes. The stim column in the metadata records which condition each cell came from.

R

table(ifnb$stim)

Expected output

CTRL STIM
6548 7451

Roughly 6,500 control cells and 7,500 stimulated cells.

Preparing Data

In Seurat v5, multi-sample data is handled by splitting the assay into layers – one layer per sample. This keeps everything in a single object while allowing integration methods to process each sample separately.

R

ifnb[["RNA"]] <- split(ifnb[["RNA"]], f = ifnb$stim)
ifnb

Expected output

An object of class Seurat
14053 features across 13999 samples within 1 assay
Active assay: RNA (14053 features, 0 variable features)
 4 layers present: counts.CTRL, counts.STIM, data.CTRL, data.STIM

The counts are now stored in two separate layers: counts.CTRL and counts.STIM.

Now run the standard preprocessing steps. When the assay is split, NormalizeData and FindVariableFeatures process each layer independently:

R

ifnb <- NormalizeData(ifnb)
ifnb <- FindVariableFeatures(ifnb, selection.method = "vst", nfeatures = 2000)
ifnb <- ScaleData(ifnb)
ifnb <- RunPCA(ifnb)

Let’s see what happens if we run UMAP and clustering without integration:

R

ifnb <- FindNeighbors(ifnb, dims = 1:30, reduction = "pca")
ifnb <- FindClusters(ifnb, resolution = 0.5)
ifnb <- RunUMAP(ifnb, dims = 1:30, reduction = "pca")

R

DimPlot(ifnb, reduction = "umap", group.by = "stim")

UMAP plot of the IFNB dataset before integration showing CTRL and STIM cells forming separate clusters driven by batch effects rather than cell type identity.

You should see a clear separation between control (CTRL) and stimulated (STIM) cells. Instead of clustering by cell type (T cells with T cells, monocytes with monocytes), the cells cluster by condition. This means we cannot directly compare the same cell type between conditions – which is exactly what we want to do. This is the batch effect problem.

Running Integration

Seurat v5 provides IntegrateLayers() as a unified interface for dataset integration. We will use the CCA (Canonical Correlation Analysis) method, which identifies shared sources of variation between the two conditions and uses them to align the data.

CCA works by:

Finding canonical correlation vectors – directions in gene expression space along which the two datasets are maximally correlated
Identifying anchor pairs – cells from different conditions that are each other’s mutual nearest neighbors in the shared CCA space
Using the anchors to compute a correction vector that aligns matching cell types across conditions

Animation: Why integrate, and what CCA integration does

R

ifnb <- IntegrateLayers(object = ifnb,
                        method = CCAIntegration,
                        orig.reduction = "pca",
                        new.reduction = "integrated.cca")

Parameter	Value	Description
`method`	`CCAIntegration`	Use Canonical Correlation Analysis for integration.
`orig.reduction`	`"pca"`	The unintegrated PCA reduction to use as input.
`new.reduction`	`"integrated.cca"`	Name for the new integrated reduction that will be stored in the object.

Now re-run the downstream steps using the integrated reduction instead of the original PCA:

R

ifnb <- FindNeighbors(ifnb, reduction = "integrated.cca", dims = 1:30)
ifnb <- FindClusters(ifnb, resolution = 0.5)
ifnb <- RunUMAP(ifnb, reduction = "integrated.cca", dims = 1:30)

Visualize the integrated result:

R

DimPlot(ifnb, reduction = "umap", group.by = "stim")

UMAP plot of the IFNB dataset after CCA integration showing control and stimulated cells intermingled within each cluster, indicating successful batch correction.

The control and stimulated cells should now be intermingled within each cluster. Cell types from both conditions co-cluster, which is what we want.

Let’s verify with a split view:

R

DimPlot(ifnb, reduction = "umap", split.by = "stim")

Split UMAP view showing CTRL and STIM conditions side by side after integration, with both panels displaying the same 14 clusters in identical spatial arrangement.

Both panels should show the same overall structure, with the same clusters present in both conditions. This confirms that integration successfully aligned shared cell types.

Before proceeding to analysis, rejoin the layers so that expression data from both conditions is accessible in a single matrix:

R

ifnb[["RNA"]] <- JoinLayers(ifnb[["RNA"]])

Callout

Other integration methods

Seurat v5 supports several integration methods through the same IntegrateLayers() interface. You can swap CCAIntegration for any of these:

Method	Function	Strengths
CCA	`CCAIntegration`	Best for datasets with shared cell types. Robust to differences in cell type composition between batches. The default choice for most experiments.
Harmony	`HarmonyIntegration`	Very fast, scales well to large datasets (100k+ cells). Uses iterative soft clustering. Requires the `harmony` R package. Widely used and well-benchmarked.
RPCA	`RPCAIntegration`	Faster than CCA (uses reciprocal PCA instead of full CCA). Good when datasets are large or when CCA is too slow. May be less robust when cell type composition differs substantially between batches.
scVI	`scVIIntegration`	Deep learning approach. Excellent for complex batch structures (many batches, large composition differences). Requires Python and the `scvi-tools` package. Slower to run but often produces superior results on difficult datasets.

For most standard experiments with two to five batches and shared cell types, CCA or Harmony work well. Use scVI for more complex scenarios (many batches, different tissues, cross-species integration).

Callout

When NOT to integrate

Integration assumes that the same cell types exist in all batches and aligns them. This is the wrong approach when:

All cells are expected to differ between conditions. For example, comparing tumor cells to normal epithelial cells – these are fundamentally different populations and should not be forced to co-cluster.
You want to find condition-specific cell types. Integration pushes cells together, which can mask the appearance of a cell state that exists in only one condition.
The batch effect is minimal. If cells already mix well by condition on the UMAP (as in some well-controlled experiments), integration is unnecessary and may introduce artifacts.

Always visualize the data before integration to assess whether batch effects are actually present. If conditions already mix well, skip integration and analyze directly.

Analysis on Integrated Data

Now that the data is integrated, we can label cell types and compare between conditions.

Cell type labels

The ifnb object comes with cell type labels: the seurat_annotations column in its metadata, assigned by the Seurat developers from marker genes. We use them as our cell types instead of labelling the 14 integrated clusters by hand. Episode 6 showed the manual approach, and Challenge 3 at the end of this episode repeats it on these clusters.

R

# Use the cell type labels that ship with the ifnb dataset
ifnb$celltype <- ifnb$seurat_annotations
Idents(ifnb) <- "celltype"
table(ifnb$celltype)

Expected output

OUTPUT

   CD14 Mono  CD4 Naive T CD4 Memory T    CD16 Mono            B        CD8 T
        4362         2504         1762         1044          978          814
 T activated           NK           DC  B Activated           Mk          pDC
         633          619          472          388          236          132
       Eryth
          55

Thirteen cell types. They are finer than the ten of episode 6 in some places (naive and memory CD4 T cells, resting and activated B cells, an activated T cell group) and include megakaryocytes (Mk, the platelet precursors) and a few erythrocytes (Eryth). CD16 Mono are the FCGR3A+ monocytes of episode 6 (FCGR3A is the gene that encodes CD16).

Checking the integration against the labels

The labels were assigned independently of our integration, so they give a check on it: if integration grouped cells by type, each integrated cluster should consist mostly of one cell type, drawn from both conditions.

R

table(cluster = ifnb$seurat_clusters, celltype = ifnb$celltype)

Expected output

OUTPUT

       celltype
cluster CD14 Mono CD4 Naive T CD4 Memory T CD16 Mono    B CD8 T T activated
     0       4272           0            0        16    1     0           0
     1          0        2332           47         1    1     9          20
     2          5         132         1643         0    0    30          16
     3          0           0            0         0  975     1          15
     4         28           0            0      1027    0     0           0
     5          1           2           51         0    0   765           7
     6          0          38            8         0    0     1         575
     7          0           0            0         0    0     8           0
     8         48           0            0         0    0     0           0
     9          0           0           13         0    1     0           0
     10         0           0            0         0    0     0           0
     11         3           0            0         0    0     0           0
     12         5           0            0         0    0     0           0
     13         0           0            0         0    0     0           0
       celltype
cluster   NK   DC B Activated   Mk  pDC Eryth
     0     0    4           1    2    0     0
     1     0    0           0    6    0     0
     2     1    0           0    1    0     0
     3     0    0         387    0    0     1
     4     0    0           0   13    0     0
     5     9    0           0    0    0     0
     6     0    0           0    0    1     0
     7   608    0           0    0    0     0
     8     0  466           0    0    1     0
     9     1    0           0  175    0     0
     10    0    2           0    0  106     0
     11    0    0           0    0    0    54
     12    0    0           0   39    0     0
     13    0    0           0    0   24     0

Almost every cluster is dominated by one cell type: cluster 0 is CD14 monocytes, 3 is B cells (resting and activated together), 4 is CD16 monocytes, 5 CD8 T cells, 7 NK cells, 8 dendritic cells, 10 pDCs and 11 erythrocytes. Clusters 1 and 2 split the CD4 T cells into mostly naive and mostly memory cells, and the megakaryocytes fall into two small clusters (9 and 12). A cluster that mixed, say, T cells and monocytes would point to an integration or clustering problem. Together with the split view above, which shows every cluster in both conditions, this confirms the integration.

R

DimPlot(ifnb, reduction = "umap", label = TRUE, repel = TRUE) + NoLegend()

Annotated UMAP plot of the integrated IFNB dataset with thirteen labeled cell types: CD14 Mono, CD4 Naive T, CD4 Memory T, CD16 Mono, B, CD8 T, T activated, NK, DC, B Activated, Mk, pDC and Eryth.

Conserved markers

Conserved markers are genes that are markers for a cell type in both conditions. FindConservedMarkers() tests for differential expression within each condition separately and then combines the results:

R

conserved <- FindConservedMarkers(ifnb,
                                  ident.1 = "CD14 Mono",
                                  grouping.var = "stim",
                                  only.pos = TRUE)


head(conserved, 10)

Parameter	Value	Description
`ident.1`	`"CD14 Mono"`	The cell type to find markers for.
`grouping.var`	`"stim"`	The metadata column defining conditions. The test is run within each level of this variable.
`only.pos`	`TRUE`	Only return upregulated genes.

The output has columns for each condition (e.g., CTRL_avg_log2FC, STIM_avg_log2FC) showing that the marker is consistent across both. Genes like S100A8, S100A9, TYROBP and FCER1G appear as strong conserved markers for monocytes.

Expected output

OUTPUT

         CTRL_p_val CTRL_avg_log2FC CTRL_pct.1 CTRL_pct.2 CTRL_p_val_adj
C15orf48          0        3.056281      0.984      0.189              0
S100A8            0        5.231200      0.853      0.063              0
IL8               0        4.008658      0.911      0.133              0
TYROBP            0        2.505484      0.997      0.249              0
FCER1G            0        2.375142      0.992      0.251              0
S100A9            0        5.162259      0.795      0.060              0
CTSB              0        3.378535      0.873      0.174              0
TYMP              0        2.359353      0.961      0.272              0
SOD2              0        2.418683      0.963      0.276              0
S100A4            0        2.704430      0.953      0.270              0
         STIM_p_val STIM_avg_log2FC STIM_pct.1 STIM_pct.2 STIM_p_val_adj
C15orf48          0        3.110432      0.988      0.191              0
S100A8            0        6.080987      0.523      0.016              0
IL8               0        4.025185      0.343      0.025              0
TYROBP            0        2.550902      0.999      0.227              0
FCER1G            0        2.753517      0.991      0.210              0
S100A9            0        5.758266      0.686      0.030              0
CTSB              0        3.791644      0.944      0.174              0
TYMP              0        2.554257      0.998      0.441              0
SOD2              0        3.103627      0.991      0.319              0
S100A4            0        2.986390      0.928      0.186              0
         max_pval minimump_p_val
C15orf48        0              0
S100A8          0              0
IL8             0              0
TYROBP          0              0
FCER1G          0              0
S100A9          0              0
CTSB            0              0
TYMP            0              0
SOD2            0              0
S100A4          0              0

Condition-specific differential expression

The most interesting biological question is: which genes change between control and stimulated cells within the same cell type? To answer this, we subset to a single cell type and compare conditions:

R

# Subset to CD14+ monocytes
mono <- subset(ifnb, idents = "CD14 Mono")

# Set the condition as the active identity
Idents(mono) <- "stim"

# Find DE genes between stimulated and control
mono.de <- FindMarkers(mono,
                       ident.1 = "STIM",
                       ident.2 = "CTRL",
                       min.pct = 0.25,
                       logfc.threshold = 0.25)
head(mono.de, 10)

Expected output

OUTPUT

        p_val avg_log2FC pct.1 pct.2 p_val_adj
IFIT1       0   7.319139 0.985 0.033         0
CXCL10      0   8.036564 0.984 0.035         0
RSAD2       0   6.741673 0.988 0.045         0
TNFSF10     0   6.991279 0.989 0.047         0
IFIT3       0   6.883785 0.992 0.056         0
IFIT2       0   7.179929 0.961 0.039         0
CXCL11      0   8.624208 0.932 0.012         0
CCL8        0   9.134191 0.918 0.017         0
IDO1        0   5.455898 0.965 0.089         0
MX1         0   5.059052 0.960 0.093         0

The top DE genes are interferon-stimulated genes (ISGs) such as IFIT1, RSAD2, TNFSF10, IFIT3, IFIT2 and MX1, together with the IFN-inducible chemokines CXCL10, CXCL11 and CCL8. These are exactly the genes you would expect to be upregulated by IFN-beta stimulation, confirming that our integration preserved the biological signal while correcting the batch effect. The exact fold changes, and the number of genes that pass a significance cutoff, vary slightly between Seurat versions; the leading ISGs do not.

Save the annotated ifnb object for the next episode:

R

saveRDS(ifnb, file = "ifnb_annotated.rds")

Visualizing condition-specific changes

Let’s visualize a few ISGs to see the condition-specific response:

R

FeaturePlot(ifnb,
            features = c("IFIT3", "IFIT1"),
            split.by = "stim",
            cols = c("grey85", "firebrick"))

FeaturePlot of ISG genes split by condition showing near-absent expression in CTRL panels and strong red expression in STIM panels, especially in monocyte clusters.

R

VlnPlot(ifnb,
        features = c("IFIT1", "IFIT3", "TNFSF10", "RSAD2"),
        split.by = "stim",
        idents = "CD14 Mono",
        ncol = 4,
        pt.size = 0)

Violin plots of IFIT1, IFIT3, TNFSF10, and RSAD2 in CD14 monocytes showing flat expression in control and broad high expression in stimulated cells.

IFIT1 and IFIT3 are strongly induced in the stimulated condition across multiple cell types, with particularly strong expression in monocytes. This demonstrates that integration corrected the batch effect (cell types co-cluster) while preserving the biological effect of IFN-beta stimulation (ISGs are still differentially expressed).

Challenge

Challenge 1: Quantifying integration quality

After integration, check whether control and stimulated cells of the same type are well mixed within each cluster. Use DimPlot(split.by = "stim") to visualize, then compute the proportion of each condition per cluster.

R

# Split UMAP by condition
DimPlot(ifnb, reduction = "umap", split.by = "stim", label = TRUE, repel = TRUE) +
    NoLegend()

# Proportion table: condition per cell type
prop.table(table(Condition = ifnb$stim, CellType = Idents(ifnb)), margin = 2)

Are the proportions roughly balanced (close to 50/50) within each cell type?

Show me the solution

OUTPUT

         CellType
Condition CD14 Mono CD4 Naive T CD4 Memory T CD16 Mono         B     CD8 T
     CTRL 0.5077946   0.3905751    0.4875142 0.4856322 0.4161554 0.4324324
     STIM 0.4922054   0.6094249    0.5124858 0.5143678 0.5838446 0.5675676
         CellType
Condition T activated        NK        DC B Activated        Mk       pDC
     CTRL   0.4739336 0.4814216 0.5466102   0.4768041 0.4872881 0.3863636
     STIM   0.5260664 0.5185784 0.4533898   0.5231959 0.5127119 0.6136364
         CellType
Condition     Eryth
     CTRL 0.4181818
     STIM 0.5818182

Split UMAP showing CTRL and STIM panels with cell type labels, demonstrating balanced mixing of conditions within each annotated cell type.

Most cell types fall in the 0.40–0.60 range, and the split UMAP shows every cell type in both panels with similar density. The baseline is not 50/50: STIM has 7,451 of the 13,999 cells (53%), so about 0.47 CTRL and 0.53 STIM is what perfect mixing looks like here.

Some observations worth discussing:

CD14 Mono, the largest population, is close to the baseline (50.8/49.2).
CD4 Naive T and pDC show the largest deviation (about 39/61), with more cells in the stimulated condition. This could reflect a real biological shift in culture rather than a failure of integration; a few percentage points either way are well within what two samples of one experiment give.
Rare populations (Mk, Eryth, pDC) have small cell numbers, so chance has more effect on their proportions.

The key diagnostic: if integration had failed, you would see a cell type that is 90%+ from one condition, or the same cell type forming separate clusters in the two panels. Neither is the case here.

Challenge

Challenge 2: Interferon response genes in monocytes

Find the top 10 differentially expressed genes between stimulated and control CD14+ monocytes (we computed this above as mono.de). Are these genes consistent with an interferon stimulation response?

R

# Top 10 upregulated genes in stimulated vs control CD14+ monocytes
top10_de <- head(mono.de[order(mono.de$avg_log2FC, decreasing = TRUE), ], 10)
print(top10_de)

# Visualize the top 4 genes
FeaturePlot(ifnb,
            features = rownames(top10_de)[1:4],
            split.by = "stim",
            ncol = 4,
            cols = c("grey85", "firebrick"))

Show me the solution

FeaturePlot of CCL8, CXCL11, HESX1, and CXCL10 split by condition showing absent expression in control and strong monocyte-concentrated expression in stimulated cells.

OUTPUT

                p_val avg_log2FC pct.1 pct.2    p_val_adj
CCL8     0.000000e+00   9.134191 0.918 0.017  0.00000e+00
CXCL11   0.000000e+00   8.624208 0.932 0.012  0.00000e+00
HESX1   2.729040e-227   8.265045 0.388 0.001 3.83512e-223
CXCL10   0.000000e+00   8.036564 0.984 0.035  0.00000e+00
IFIT1    0.000000e+00   7.319139 0.985 0.033  0.00000e+00
IFIT2    0.000000e+00   7.179929 0.961 0.039  0.00000e+00
TNFSF10  0.000000e+00   6.991279 0.989 0.047  0.00000e+00
CD38    4.872419e-277   6.974732 0.464 0.005 6.84721e-273
IFIT3    0.000000e+00   6.883785 0.992 0.056  0.00000e+00
GMPR    4.371764e-244   6.765773 0.419 0.005 6.14364e-240

The top 10 genes fall into two functional categories, both consistent with IFN-beta stimulation:

Gene	log2FC	pct.1 / pct.2	Function
CCL8	9.13	0.92 / 0.02	Chemokine; recruits monocytes and T cells to sites of inflammation
CXCL11	8.62	0.93 / 0.01	IFN-inducible chemokine; attracts activated T cells via CXCR3
HESX1	8.27	0.39 / 0.00	Transcription factor; emerging role in IFN signaling
CXCL10	8.04	0.98 / 0.04	IFN-inducible chemokine (IP-10); major recruiter of immune cells
IFIT1	7.32	0.99 / 0.03	Directly induced by type I IFN; inhibits viral translation
IFIT2	7.18	0.96 / 0.04	Antiviral effector; regulates translation during infection
TNFSF10	6.99	0.99 / 0.05	TRAIL; IFN-induced death ligand that kills infected cells
CD38	6.97	0.46 / 0.01	NAD-consuming enzyme induced by type I IFN in monocytes
IFIT3	6.88	0.99 / 0.06	Antiviral effector; forms complex with IFIT1 and IFIT2
GMPR	6.77	0.42 / 0.01	Guanosine monophosphate reductase; modulates purine metabolism

Two patterns stand out:

Chemokines dominate the top. CCL8, CXCL11, and CXCL10 are three of the top four, with fold changes of 8.0-9.1, reflecting the monocyte’s primary role as a sentinel cell: upon IFN stimulation, monocytes broadcast chemokine signals to recruit other immune cells. Note the pct.2 values (0.01-0.04) – these chemokines are essentially absent in control cells and massively induced by stimulation.

Classical ISGs fill the rest. IFIT1/2/3 and TNFSF10 are canonical interferon-stimulated genes; CD38 and GMPR are induced in a subset of cells (pct.1 0.42-0.46). The IFIT family shows the highest pct.1 values (0.96-0.99), meaning nearly every stimulated monocyte expresses them – a hallmark of a robust, uniform interferon response.

The FeaturePlot confirms this visually: CCL8, CXCL11, HESX1, and CXCL10 are virtually absent in the CTRL panels (gray) and strongly expressed in the STIM panels (red), with the signal concentrated in the monocytes (right side of the UMAP); HESX1 is induced in only a subset of them. CXCL10 also reaches the stimulated B cells and is scattered through the T cells, consistent with these cell types also responding to IFN-beta, though less strongly than monocytes.

Challenge

Challenge 3: Label the integrated clusters yourself

We used the labels that ship with the dataset. Label the 14 integrated clusters from marker genes instead, as in episode 6, and compare your labels with seurat_annotations. Plot the canonical markers per cluster first, then write one label per cluster.

R

Idents(ifnb) <- "seurat_clusters"
DotPlot(ifnb, features = c("CD3D", "IL7R", "CCR7", "CD8A", "GNLY", "NKG7",
                           "MS4A1", "CD79A", "CD14", "LYZ", "FCGR3A", "MS4A7",
                           "FCER1A", "LILRA4", "IL3RA", "PPBP", "HBB")) +
    coord_flip()

Callout

A label vector only fits the clusters it was written for

The solution below fits the 14 clusters this lesson’s integration gives with the package versions in the setup instructions. With other versions the number or order of the clusters can change, and the names then land on the wrong cells without any error. The code stops if the number of clusters does not match the number of labels; it cannot detect a change in order, so check your dot plot against the comment beside each label.

Show me the solution

R

# One label per integrated cluster, in cluster order (0 to 13)
my.labels <- c(
    "CD14 Mono",     # 0: CD14+, LYZ+, CST3+
    "CD4 Naive T",   # 1: CD3D+, CCR7 high
    "CD4 Memory T",  # 2: CD3D+, IL7R+, CCR7 lower
    "B",             # 3: MS4A1+, CD79A+
    "CD16 Mono",     # 4: FCGR3A+, MS4A7+
    "CD8 T",         # 5: CD3D+, CD8A+, NKG7+
    "T activated",   # 6: IL7R+, CCR7+, CD3D low
    "NK",            # 7: GNLY+, NKG7+, CD3D-
    "DC",            # 8: LYZ+, CST3+, FCER1A+
    "Mk",            # 9: PPBP+
    "pDC",           # 10: LILRA4+, IL3RA+
    "Eryth",         # 11: HBB+
    "Mk",            # 12: PPBP+, with monocyte genes
    "pDC"            # 13: no clear marker; the shipped labels say pDC
)
if (nlevels(ifnb$seurat_clusters) != length(my.labels)) stop("ifnb has ", nlevels(ifnb$seurat_clusters), " clusters but my.labels has ", length(my.labels), " labels; check the dot plot and edit my.labels")
names(my.labels) <- levels(ifnb$seurat_clusters)
ifnb$manual_celltype <- unname(my.labels[as.character(ifnb$seurat_clusters)])
table(Manual = ifnb$manual_celltype, Shipped = ifnb$celltype)
Idents(ifnb) <- "celltype"

OUTPUT

              Shipped
Manual         CD14 Mono CD4 Naive T CD4 Memory T CD16 Mono    B CD8 T
  B                    0           0            0         0  975     1
  CD14 Mono         4272           0            0        16    1     0
  CD16 Mono           28           0            0      1027    0     0
  CD4 Memory T         5         132         1643         0    0    30
  CD4 Naive T          0        2332           47         1    1     9
  CD8 T                1           2           51         0    0   765
  DC                  48           0            0         0    0     0
  Eryth                3           0            0         0    0     0
  Mk                   5           0           13         0    1     0
  NK                   0           0            0         0    0     8
  pDC                  0           0            0         0    0     0
  T activated          0          38            8         0    0     1
              Shipped
Manual         T activated   NK   DC B Activated   Mk  pDC Eryth
  B                     15    0    0         387    0    0     1
  CD14 Mono              0    0    4           1    2    0     0
  CD16 Mono              0    0    0           0   13    0     0
  CD4 Memory T          16    1    0           0    1    0     0
  CD4 Naive T           20    0    0           0    6    0     0
  CD8 T                  7    9    0           0    0    0     0
  DC                     0    0  466           0    0    1     0
  Eryth                  0    0    0           0    0    0    54
  Mk                     0    1    0           0  214    0     0
  NK                     0  608    0           0    0    0     0
  pDC                    0    0    2           0    0  130     0
  T activated          575    0    0           0    0    1     0

Each cluster gets the shipped label that most of its cells carry, so most cells sit on the diagonal of the table. The off-diagonal cells are the ones a cluster-level label cannot get right: activated B cells inside the B cell cluster, the naive and memory CD4 T cells that the two CD4 clusters share, and a few cells of other types in each cluster. Cluster 13 has no clear marker in the dot plot; with only 24 cells, it is labelled pDC because the shipped annotation says so. Two lessons: a cluster label is a statement about most of its cells, not all of them, and a cluster without markers is a candidate for removal rather than a cell type.

The last line sets the identities back to the shipped labels, so any code you run afterwards uses the same cell types as the rest of the lesson.

Key Points

Batch effects cause cells to cluster by sample rather than by cell type, preventing cross-condition comparisons
CCA integration finds shared correlation structure between conditions and aligns matching cell types using anchor pairs
In Seurat v5, multi-sample data is handled by splitting assay layers, integrating, then rejoining layers for downstream analysis
Always verify integration by checking that conditions mix within clusters AND that known biological differences (like ISG expression) are preserved
Integration is not always appropriate; skip it when conditions have no shared cell types or when batch effects are minimal

Content from Differential Expression Analysis

Last updated on 2026-10-27 | Edit this page

Estimated time: 60 minutes

Overview

Questions

What are the main approaches for differential expression testing in scRNA-seq?
Why is pseudobulk analysis preferred over per-cell testing for multi-sample experiments?
How do we create volcano plots and other visualizations to interpret DE results?
How do we connect lists of DE genes to biological pathways using GO enrichment?
How should we export our results for sharing and reproducibility?

Objectives

Compare Wilcoxon, MAST, and pseudobulk DE testing approaches and choose the right one for your design
Run pseudobulk differential expression using AggregateExpression
Create volcano plots, feature plots, and violin plots to visualize DE results
Perform gene ontology enrichment analysis with clusterProfiler and interpret the results
Export Seurat objects, count matrices, and metadata for reproducibility and sharing

Prerequisite

Prerequisites

This episode requires an RStudio session on the Negishi cluster. Launch RStudio (bioconductor) under Bioinformatics Apps on Open OnDemand as described in the Setup instructions. You will also need the ifnb_annotated.rds object saved at the end of the previous episode.

Setup

R

library(Seurat)
library(ggplot2)
library(patchwork)
library(dplyr)
library(clusterProfiler)
library(org.Hs.eg.db)

Loading the Integrated Data

We continue from the integrated and annotated IFNB dataset produced in the previous episode. Make sure the cell type annotation is set as the active identity.

Set up the working directory to match the previous episodes:

R

work_dir <- paste0(
    "/scratch/negishi/", Sys.getenv("USER"),
    "/scrna_workshop/"
)
setwd(work_dir)

Callout

Starting here without the previous episode?

Copy the saved objects from the workshop data first (pbmc_annotated.rds is used by the export section at the end): file.copy(paste0(work_dir, c("rds/ifnb_annotated.rds", "rds/pbmc_annotated.rds")), work_dir)

R

ifnb <- readRDS(paste0(work_dir, "ifnb_annotated.rds"))
Idents(ifnb) <- "celltype"
DimPlot(ifnb, reduction = "umap", label = TRUE, repel = TRUE, split.by = "stim") +
    NoLegend()

Expected output

UMAP split by CTRL and STIM conditions with thirteen cell types colored and labeled, showing identical spatial arrangement in both panels confirming successful integration.

You should see two UMAP panels (CTRL and STIM) with the same cell type structure in both. Cell types are labeled and colored consistently across panels.

DE Testing Methods

Seurat supports several statistical tests for differential expression. The choice of method depends on your experimental design and whether you have biological replicates.

Method	Seurat parameter	Statistical model	Best for
Wilcoxon rank-sum	`test.use = "wilcox"`	Non-parametric rank test	Fast exploration, comparing clusters within a single sample
MAST	`test.use = "MAST"`	Hurdle model (logistic + Gaussian)	Accounts for dropout and cellular detection rate; per-cell testing with better modeling
DESeq2 (pseudobulk)	`test.use = "DESeq2"`	Negative binomial GLM on aggregated counts	Publication-quality results when biological replicates are available

Let’s start with the Wilcoxon test to quickly identify DE genes between stimulated and control CD14+ monocytes:

R

# Subset to CD14+ monocytes
mono <- subset(ifnb, idents = "CD14 Mono")

# Set condition as active identity
Idents(mono) <- "stim"

# Run Wilcoxon DE test
mono.wilcox <- FindMarkers(mono,
                           ident.1 = "STIM",
                           ident.2 = "CTRL",
                           test.use = "wilcox",
                           min.pct = 0.25,
                           logfc.threshold = 0.25)
head(mono.wilcox, 10)

Parameter	Value	Description
`ident.1`	`"STIM"`	The first group (numerator). Positive `avg_log2FC` means higher in STIM.
`ident.2`	`"CTRL"`	The second group (denominator).
`test.use`	`"wilcox"`	Wilcoxon rank-sum test. Fast, non-parametric, works well for exploration.
`min.pct`	`0.25`	Only test genes detected in at least 25% of cells in either group.
`logfc.threshold`	`0.25`	Only test genes with at least 0.25 log2FC between groups.

R

cat("Significant genes (padj < 0.05):", sum(mono.wilcox$p_val_adj < 0.05), "\n")

Expected output

OUTPUT

        p_val avg_log2FC pct.1 pct.2 p_val_adj
IFIT1       0   7.319139 0.985 0.033         0
CXCL10      0   8.036564 0.984 0.035         0
RSAD2       0   6.741673 0.988 0.045         0
TNFSF10     0   6.991279 0.989 0.047         0
IFIT3       0   6.883785 0.992 0.056         0
IFIT2       0   7.179929 0.961 0.039         0
CXCL11      0   8.624208 0.932 0.012         0
CCL8        0   9.134191 0.918 0.017         0
IDO1        0   5.455898 0.965 0.089         0
MX1         0   5.059052 0.960 0.093         0


Significant genes (padj < 0.05): 825

You should see approximately 825 significant genes at padj < 0.05. The top genes include interferon-stimulated genes (IFIT1/2/3, RSAD2, MX1) and IFN-inducible chemokines (CXCL10, CXCL11, CCL8). But notice that the p-values are extremely small (effectively zero). This is a symptom of per-cell testing with thousands of cells – even modest differences become “significant.”

Callout

The pseudobulk revolution

The Wilcoxon and MAST tests treat each cell as an independent observation. But cells from the same donor or the same sample are not independent – they share the same genetic background, the same library preparation, and the same sequencing run. This is called pseudoreplication, and it severely inflates p-values.

For example, if you have 2,000 monocytes from one control sample and 2,000 from one stimulated sample, per-cell testing treats this as n = 4,000. But the true sample size is n = 2 (one sample per condition). Every gene with even a tiny difference will appear significant simply because of the enormous “sample size.”

Pseudobulk analysis solves this by aggregating all cells of the same type from the same sample into a single expression profile. If you have 3 control donors and 3 stimulated donors, you get 3 pseudobulk replicates per condition – the same as a bulk RNA-seq experiment with n = 3 per group. This gives proper type I error control and realistic p-values.

Always use pseudobulk when you have biological replicates. Use per-cell tests only for exploratory analysis or when you have a single sample per condition and have no other option.

Animation: Cells are not replicates

Pseudobulk Analysis

Pseudobulk aggregation sums the raw UMI counts for each gene across all cells of the same cell type from the same sample, producing a bulk-like count matrix with one column per sample-celltype combination.

Callout

The IFNB dataset has no biological replicates

The IFNB dataset has only one orig.ident per condition (IMMUNE_CTRL and IMMUNE_STIM). If we aggregate by celltype + stim + orig.ident, we get one pseudobulk sample per group – not enough for DESeq2, which requires at least 3.

To demonstrate the pseudobulk workflow, we create pseudo-replicates by randomly splitting cells within each condition into 3 groups. This is a common workaround for teaching and exploration, but be aware: pseudo-replicates do not capture true biological variability across donors. For publication-quality DE, you need genuine biological replicates (multiple donors or independent experiments).

R

# Create pseudo-replicates by randomly assigning cells to 3 groups per condition
set.seed(42)
ifnb$pseudo_rep <- paste0("rep", sample(1:3, ncol(ifnb), replace = TRUE))

# Verify: should see 3 roughly equal groups per condition
table(ifnb$stim, ifnb$pseudo_rep)

OUTPUT

       rep1 rep2 rep3
  CTRL 2198 2198 2152
  STIM 2533 2465 2453

Now aggregate using the pseudo-replicate labels:

R

# Create pseudobulk counts
# Group by cell type, condition, and pseudo-replicate
bulk <- AggregateExpression(ifnb,
                            group.by = c("celltype", "stim", "pseudo_rep"),
                            return.seurat = TRUE)
bulk

Parameter	Value	Description
`group.by`	`c("celltype", "stim", "pseudo_rep")`	Metadata columns to group cells by. Each unique combination of cell type, condition, and pseudo-replicate becomes one pseudobulk column.
`return.seurat`	`TRUE`	Return a Seurat object (instead of a raw matrix) so we can use `FindMarkers` on it.

Expected output

OUTPUT

An object of class Seurat
14053 features across 78 samples within 1 assay
Active assay: RNA (14053 features, 0 variable features)
 3 layers present: counts, data, scale.data

The number of samples equals the number of unique cell type x condition x pseudo-replicate combinations. With 13 cell types, 2 conditions, and 3 pseudo-replicates, you get up to 78 pseudobulk samples (fewer if a cell type has no cells in a given replicate).

Now we can run DESeq2-based DE testing on the pseudobulk object. We subset to CD14+ monocytes and compare conditions:

R

# Subset pseudobulk to CD14 Mono samples
bulk.mono <- subset(bulk, celltype == "CD14 Mono")

# Set condition as active identity
Idents(bulk.mono) <- "stim"

# Run DESeq2 on pseudobulk counts
mono.pseudo <- FindMarkers(bulk.mono,
                           ident.1 = "STIM",
                           ident.2 = "CTRL",
                           test.use = "DESeq2")
head(mono.pseudo, 10)

R

cat("Significant genes (padj < 0.05):", sum(mono.pseudo$p_val_adj < 0.05, na.rm = TRUE), "\n")

Expected output

OUTPUT

        p_val avg_log2FC pct.1 pct.2 p_val_adj
ISG15       0   6.526243     1     1         0
IFI6        0   3.627521     1     1         0
GBP1        0   2.583049     1     1         0
S100A8      0  -2.269779     1     1         0
RSAD2       0   5.939389     1     1         0
TMSB10      0   1.705161     1     1         0
VAMP5       0   3.068289     1     1         0
SSB         0   2.631016     1     1         0
PLSCR1      0   2.452130     1     1         0
TNFSF10     0   6.054715     1     1         0


Significant genes (padj < 0.05): 1539

The pseudobulk analysis identifies 1539 significant genes – more than the Wilcoxon test (825). This may seem surprising: we said pseudobulk is more conservative, so why does it find more? Because we used pseudo-replicates, not true biological replicates. The 3 random splits within each condition have very low within-group variance, giving DESeq2 high statistical power.

With genuine biological replicates (different donors with real between-donor variability), the pseudobulk gene count would typically be much lower than the per-cell test. The high count here is an artifact of pseudo-replication and should not be interpreted as pseudobulk being more sensitive in general.

Visualization and Interpretation

Volcano plot

A volcano plot shows the relationship between statistical significance (-log10 p-value, y-axis) and biological effect size (log2 fold change, x-axis). Each point is a gene. Genes in the upper corners are both statistically significant and biologically meaningful.

R

# Use Wilcoxon results for visualization (more genes to plot)
mono.wilcox$gene <- rownames(mono.wilcox)
mono.wilcox$significance <- case_when(
    mono.wilcox$p_val_adj < 0.05 & mono.wilcox$avg_log2FC > 0.5  ~ "Up",
    mono.wilcox$p_val_adj < 0.05 & mono.wilcox$avg_log2FC < -0.5 ~ "Down",
    TRUE ~ "NS"
)

ggplot(mono.wilcox, aes(x = avg_log2FC, y = -log10(p_val_adj), color = significance)) +
    geom_point(size = 0.5, alpha = 0.6) +
    scale_color_manual(values = c("Up" = "firebrick", "Down" = "steelblue", "NS" = "grey70")) +
    geom_vline(xintercept = c(-0.5, 0.5), linetype = "dashed", color = "grey40") +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey40") +
    labs(x = "Log2 fold change (STIM vs CTRL)",
         y = "-log10(adjusted p-value)",
         title = "CD14+ Monocytes: STIM vs CTRL",
         color = "Direction") +
    theme_minimal()

Volcano plot showing upregulated genes in red in the upper right, downregulated genes in blue in the upper left, and non-significant genes in grey near the bottom, with dashed cutoff lines for significance and fold change.

How to read the volcano plot:

Upper right: genes significantly upregulated in stimulated cells (ISGs, chemokines)
Upper left: genes significantly downregulated in stimulated cells
Bottom center: non-significant genes (low fold change or high p-value)
Dashed lines: significance cutoffs (|log2FC| > 0.5 and padj < 0.05)

Feature and violin plots

Visualize specific DE genes on the UMAP and across conditions:

R

FeaturePlot(ifnb,
            features = "IFIT1",
            split.by = "stim",
            cols = c("grey85", "firebrick"))

Feature plot of IFIT1 expression split by condition showing near-absent expression in control cells and strong expression across all cell types in stimulated cells, especially monocytes.

R

VlnPlot(ifnb,
        features = c("IFIT1", "CXCL10", "MX1"),
        split.by = "stim",
        pt.size = 0,
        ncol = 3)

Three violin plot panels showing IFIT1 and MX1 upregulated in stimulated cells of every cell type and CXCL10 mainly in stimulated monocytes, dendritic cells and B cells, with control violins near zero.

Dot plot across cell types

A dot plot can show the top DE genes across all cell types simultaneously, revealing which responses are shared and which are cell-type-specific:

R

top_de_genes <- head(rownames(mono.wilcox), 10)
DotPlot(ifnb,
        features = top_de_genes,
        split.by = "stim",
        cols = c("steelblue", "firebrick")) +
    coord_flip() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

Dot plot showing top 10 DE genes across all cell types split by CTRL and STIM, with dot size indicating percent expressed and color indicating average expression, showing consistently larger and redder dots in stimulated cells.

Functional Enrichment

A list of hundreds of DE genes is difficult to interpret gene by gene. Gene Ontology (GO) enrichment analysis asks: are the DE genes enriched for specific biological processes? This connects individual gene changes to broader biological themes.

Animation: What enrichGO tests

Converting gene symbols to Entrez IDs

clusterProfiler requires Entrez gene IDs. We convert from gene symbols using bitr():

R

# Get significant upregulated genes
sig_genes <- mono.wilcox %>%
    filter(p_val_adj < 0.05, avg_log2FC > 0.5) %>%
    rownames()

cat("Number of significant upregulated genes:", length(sig_genes), "\n")

# Convert gene symbols to Entrez IDs
gene_ids <- bitr(sig_genes,
                 fromType = "SYMBOL",
                 toType = "ENTREZID",
                 OrgDb = org.Hs.eg.db)
head(gene_ids)

Parameter	Value	Description
`fromType`	`"SYMBOL"`	Input ID type (gene symbols like ISG15, IFIT1).
`toType`	`"ENTREZID"`	Output ID type (Entrez gene IDs like 9636, 3434).
`OrgDb`	`org.Hs.eg.db`	Organism annotation database (human).

Expected output

OUTPUT

Number of significant upregulated genes: 293

   SYMBOL ENTREZID
1   IFIT1     3434
2  CXCL10     3627
3   RSAD2    91543
4 TNFSF10     8743
5   IFIT3     3437
6   IFIT2     3433

Running GO enrichment

R

go_results <- enrichGO(gene = gene_ids$ENTREZID,
                       OrgDb = org.Hs.eg.db,
                       ont = "BP",
                       pAdjustMethod = "BH",
                       qvalueCutoff = 0.05,
                       readable = TRUE)

Parameter	Value	Description
`gene`	`gene_ids$ENTREZID`	Vector of Entrez gene IDs to test.
`OrgDb`	`org.Hs.eg.db`	Organism annotation database.
`ont`	`"BP"`	GO ontology category: “BP” (Biological Process), “MF” (Molecular Function), or “CC” (Cellular Component).
`pAdjustMethod`	`"BH"`	Benjamini-Hochberg method for multiple testing correction.
`qvalueCutoff`	`0.05`	Only return terms with q-value below this threshold.
`readable`	`TRUE`	Convert Entrez IDs in the output back to gene symbols for readability.

Visualize the top enriched GO terms with a dot plot:

R

dotplot(go_results, showCategory = 15) +
    ggtitle("GO Biological Process: STIM vs CTRL CD14+ Monocytes")

Expected output

GO enrichment dot plot showing 15 biological process terms with response to virus, defense response to virus, and viral process as the top three by gene ratio, with dot color from red to blue indicating adjusted p-value.

The top enriched GO terms should include:

response to virus
defense response to virus
viral process
viral life cycle
regulation of innate immune response
regulation of viral process
viral genome replication

These terms are entirely consistent with IFN-beta stimulation. IFN-beta is a type I interferon that activates antiviral defense programs, and the GO enrichment confirms that the DE genes are concentrated in exactly these pathways.

KEGG pathway enrichment

For a complementary view, we can also test KEGG pathways:

R

kegg_results <- enrichKEGG(gene = gene_ids$ENTREZID,
                           organism = "hsa",
                           pAdjustMethod = "BH",
                           qvalueCutoff = 0.05)
head(kegg_results, 5)

Expected output

OUTPUT

                   category               subcategory       ID
hsa04621 Organismal Systems             Immune system hsa04621
hsa05168     Human Diseases Infectious disease: viral hsa05168
hsa04612 Organismal Systems             Immune system hsa04612
hsa05169     Human Diseases Infectious disease: viral hsa05169
hsa05164     Human Diseases Infectious disease: viral hsa05164
                                 Description GeneRatio  BgRatio RichFactor
hsa04621 NOD-like receptor signaling pathway    26/179 187/9451 0.13903743
hsa05168    Herpes simplex virus 1 infection    20/179 180/9451 0.11111111
hsa04612 Antigen processing and presentation    13/179  82/9451 0.15853659
hsa05169        Epstein-Barr virus infection    19/179 205/9451 0.09268293
hsa05164                         Influenza A    17/179 173/9451 0.09826590
         FoldEnrichment    zScore       pvalue     p.adjust       qvalue
hsa04621       7.341021 12.168483 9.114111e-16 2.114474e-13 2.114474e-13
hsa05168       5.866543  9.159024 1.551182e-10 1.799371e-08 1.799371e-08
hsa04612       8.370555  9.313570 3.728060e-09 2.883033e-07 2.883033e-07
hsa05169       4.893555  7.830726 9.849154e-09 5.712510e-07 5.712510e-07
hsa05164       5.188330  7.724892 2.584269e-08 1.101704e-06 1.101704e-06
                                                                                                                                           geneID
hsa04621 79792/10010/115362/2634/3320/6773/7295/10628/4939/3428/10616/8767/4938/6347/1508/3665/3326/834/115361/6772/330/4940/2633/114769/837/2635
hsa05168                                      6672/3133/6892/64135/6773/29992/4939/684/4938/6347/6890/3665/3107/8717/6772/330/4940/3134/5610/3105
hsa04612                                                                         3133/6892/3320/1520/3303/1508/6890/3107/3326/5721/3134/3105/1514
hsa05169                                            3133/6892/6773/4939/4616/3627/4938/6890/3665/3107/8717/958/6772/4940/3134/5610/3105/4067/9636
hsa05164                                                     91543/64135/6773/103/4939/3627/4938/6347/4600/3665/4599/8717/834/8743/6772/4940/5610
         Count
hsa04621    26
hsa05168    20
hsa04612    13
hsa05169    19
hsa05164    17

The top KEGG pathways are “NOD-like receptor signaling pathway”, “Herpes simplex virus 1 infection”, “Antigen processing and presentation”, “Epstein-Barr virus infection” and “Influenza A” – all involving the interferon and antiviral response.

Exporting Results

At the end of an analysis, you should export your results in formats that support reproducibility, sharing, and downstream use.

Save the annotated Seurat object

R

saveRDS(ifnb, file = paste0(work_dir, "ifnb_de.rds"))

Export the count matrix and metadata

For the PBMC 10k dataset from earlier episodes, you can export the raw counts and metadata as CSV files. First load the annotated object saved in the cell type annotation episode:

R

pbmc <- readRDS(paste0(work_dir, "pbmc_annotated.rds"))

# Export raw count matrix (warning: can be large for big datasets)
counts_matrix <- GetAssayData(pbmc, layer = "counts")
Matrix::writeMM(counts_matrix, file = "counts.mtx")

# Export metadata
write.csv(pbmc@meta.data, file = "metadata.csv", quote = FALSE)

# Export gene names
write.csv(data.frame(gene = rownames(pbmc)), file = "features.csv",
          row.names = FALSE, quote = FALSE)

A dense CSV stores every zero, so keep it to small datasets or a subset of cells. Here, the first 500 cells:

R

write.csv(as.matrix(GetAssayData(pbmc, layer = "counts")[, 1:500]),
          file = "counts_dense_500cells.csv")

Format conversion for Python users

If collaborators use the Python/Scanpy ecosystem, you can write your Seurat object as an AnnData .h5ad file with the Bioconductor package anndataR. It writes the file directly from R, so no Python installation is needed:

R

# Install anndataR if needed: BiocManager::install("anndataR")
library(anndataR)
write_h5ad(pbmc, "pbmc_annotated.h5ad", mode = "w")  # "w" overwrites an earlier export

The counts, normalized data, cell metadata and the PCA and UMAP embeddings are all carried over. Three warnings are expected. One says the scale.data layer is skipped “with unexpected dimensions”: it holds only the 2,000 variable genes, an AnnData layer must cover every gene, and scaled values are easy to recompute in Python. Two say that matrix column names “cannot be written”: AnnData stores embeddings without column names. To read the file in Python, collaborators need a recent anndata (0.13 reads it; 0.10 does not).

Callout

Why not sceasy?

Older tutorials use sceasy::convertFormat(). Its latest version (0.0.7) calls GetAssayData(slot = ...), which SeuratObject 5 no longer accepts, so it stops with an error on Seurat v5 objects.

GEO submission

For submitting data to NCBI GEO, you typically need:

Raw count matrix (genes x cells) in MTX or CSV format
Cell metadata (barcodes, cell types, conditions) as a TSV/CSV
Processing description (software versions, parameters, filtering criteria)

Callout

Reproducibility checklist

Before sharing or publishing your analysis, ensure you have:

Saved your R session info so others know exactly which package versions you used:
R
```
writeLines(capture.output(sessionInfo()), "session_info.txt")
```
Pinned package versions in your analysis scripts or used renv to create a reproducible environment.
Saved intermediate RDS files at each major step (pbmc_filtered.rds, pbmc_normalized.rds, pbmc_clustered.rds, pbmc_annotated.rds) so others can enter the pipeline at any point.
Documented all parameter choices: QC thresholds, number of PCs, clustering resolution, normalization method, integration method, and DE test used.
Recorded the random seed if you used one (e.g., set.seed(42) before UMAP or clustering).

Challenge

Challenge 1: Wilcoxon vs. Pseudobulk

Run both Wilcoxon (per-cell) and pseudobulk (DESeq2) DE testing for stimulated vs. control CD14+ monocytes. Compare the number of significant genes and the p-value distributions.

R

# Wilcoxon (already computed above)
n_wilcox <- sum(mono.wilcox$p_val_adj < 0.05)
cat("Wilcoxon significant genes:", n_wilcox, "\n")

# Pseudobulk (already computed above)
n_pseudo <- sum(mono.pseudo$p_val_adj < 0.05, na.rm = TRUE)
cat("Pseudobulk significant genes:", n_pseudo, "\n")

# Compare p-value distributions for shared genes
shared_genes <- intersect(rownames(mono.wilcox), rownames(mono.pseudo))
comparison <- data.frame(
    gene = shared_genes,
    wilcox_pval = mono.wilcox[shared_genes, "p_val_adj"],
    pseudo_pval = mono.pseudo[shared_genes, "p_val_adj"]
)
comparison <- comparison[complete.cases(comparison), ]

ggplot(comparison, aes(x = -log10(wilcox_pval), y = -log10(pseudo_pval))) +
    geom_point(size = 0.5, alpha = 0.5) +
    geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed") +
    labs(x = "-log10(p-value) Wilcoxon",
         y = "-log10(p-value) Pseudobulk",
         title = "P-value comparison: Wilcoxon vs Pseudobulk") +
    theme_minimal()

Which method gives more conservative results? Why?

Show me the solution

Scatter plot of negative log10 p-values from Wilcoxon versus pseudobulk for shared genes, with a red dashed diagonal showing equal significance and points broadly scattered indicating different significance rankings between methods.

You should find that:

Wilcoxon identifies 825 significant genes
Pseudobulk identifies 1539 (with pseudo-replicates)

Pseudobulk finds more significant genes here, which may seem to contradict the usual advice that pseudobulk is more conservative. The reason: our pseudo-replicates (random cell splits) have very low within-group variance, so DESeq2 gets very high power. With real biological replicates (different donors), the between-donor variability would inflate variance estimates and the gene count would drop substantially.

The scatter plot shows a complex pattern:

A dense band along the x-axis (Wilcoxon significant, pseudobulk not) – genes with small effects that Wilcoxon detects via per-cell power but pseudobulk misses.
Many points above the diagonal – genes where pseudobulk yields smaller p-values, because DESeq2’s count-based model is more powerful when within-group variance is low.
Some points below the diagonal – genes where Wilcoxon is more significant.

With true biological replicates, you would expect most points to shift below the diagonal and the overall pseudobulk gene count to drop well below the Wilcoxon count.

Rule of thumb: Use Wilcoxon for quick exploration and hypothesis generation. Use pseudobulk for any result you plan to report or publish. And for publication, ensure you have true biological replicates, not pseudo-replicates.

Challenge

Challenge 2: GO Enrichment of IFN-beta Response Genes

Run GO enrichment on the top 200 upregulated genes (by fold change) in stimulated vs. control CD14+ monocytes. What are the top 5 biological processes? Do they make sense for IFN-beta stimulation?

R

# Get top 200 upregulated genes by fold change
top200 <- mono.wilcox %>%
    filter(p_val_adj < 0.05, avg_log2FC > 0) %>%
    arrange(desc(avg_log2FC)) %>%
    head(200)
top200_genes <- rownames(top200)

# Convert to Entrez IDs
top200_ids <- bitr(top200_genes,
                   fromType = "SYMBOL",
                   toType = "ENTREZID",
                   OrgDb = org.Hs.eg.db)

# Run GO enrichment
go_top200 <- enrichGO(gene = top200_ids$ENTREZID,
                      OrgDb = org.Hs.eg.db,
                      ont = "BP",
                      pAdjustMethod = "BH",
                      qvalueCutoff = 0.05,
                      readable = TRUE)

# Show top 5 terms
head(go_top200, 5)

# Visualize
dotplot(go_top200, showCategory = 10) +
    ggtitle("GO BP: Top 200 upregulated genes in STIM CD14+ Monocytes")

Show me the solution

GO enrichment dot plot for top 200 upregulated genes showing 10 biological process terms with response to virus and defense response to virus having the largest gene ratios and deepest red color indicating lowest adjusted p-values.

The top 5 biological processes should include terms like:

Rank	GO Term	Interpretation
1	response to virus	Broad viral response including both detection and effector mechanisms
2	defense response to virus	Direct antiviral effector programs activated by IFN-beta
3	negative regulation of viral genome replication	ISGs that block viral replication directly (IFIT, OAS, MX family genes)
4	regulation of viral process	Regulatory genes that modulate antiviral defense (e.g., TRIM, IRF family)
5	viral genome replication	Overlaps with rank 3; genes acting on the replication step itself

These results make complete biological sense. IFN-beta (interferon beta) is a type I interferon that signals through the IFNAR receptor complex, activating the JAK-STAT pathway, which turns on hundreds of interferon-stimulated genes (ISGs). These ISGs encode proteins that directly inhibit viral replication (MX1, IFIT1/2/3), degrade viral RNA (RSAD2), recruit immune cells via chemokines (CXCL10, CXCL11, CCL8), and alert neighboring cells to the viral threat.

The GO enrichment confirms that our DE analysis is capturing the expected biology: the top upregulated genes in stimulated monocytes are overwhelmingly involved in antiviral and interferon signaling pathways. This serves as important biological validation that the computational pipeline (integration, DE testing, enrichment) is working correctly.

Key Points

Per-cell DE tests (Wilcoxon, MAST) treat each cell as independent, inflating p-values through pseudoreplication
Pseudobulk analysis aggregates cells by sample and cell type, providing statistically valid results with proper type I error control
Volcano plots visualize the relationship between effect size (log2FC) and significance (-log10 p-value) for all tested genes
Gene ontology enrichment connects lists of DE genes to biological pathways, providing functional interpretation of differential expression results
Always export session info, save intermediate objects, and document parameter choices for reproducibility