All in One View
Content from Introduction to Single-Cell RNA-Seq
Last updated on 2026-09-08 | Edit this page
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.
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.
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.
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.
4. Sequencing. The library is sequenced on an Illumina sequencer, producing paired-end reads plus an index read. The read structure is:
- 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.
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.
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.
For the PBMC 10k v3 dataset we will use in this workshop, the count matrix contains approximately 11,800 cells and 36,600 genes, but after filtering out lowly-expressed genes and low-quality cells, we will work with roughly 8,000 – 10,000 cells and 15,000 – 20,000 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.
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 Leiden algorithm 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.
- 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-09-08 | Edit this page
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 and submit a SLURM job script for Cell Ranger count
- Write and submit a SLURM job script for STARsolo
- Interpret Cell Ranger’s web_summary.html and key QC metrics
- Compare the filtered count matrices produced by both tools
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.
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.
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 resultsNow copy the pre-staged workshop data from the shared depot:
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.txtVerify that the FASTQ files are present:
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.gzVerify that the reference genome is present:
OUTPUT
3M-february-2018.txt refdata-gex-GRCh38-2024-AUnderstanding 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 nameThis dataset was sequenced across two lanes (L001, L002), producing four files total: two lanes times two reads per lane.
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:
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
+
#FFFFFFFFFFFFFFFFFFFFFFFFFFFR2 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:FFFFFFFFFFFFFFFThe 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:
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=workshop
#SBATCH --partition=cpu
#SBATCH --output=cellranger_%j.out
#SBATCH --error=cellranger_%j.err
# Load modules
ml --force purge
ml biocontainers
ml cellranger
# 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
EOFSubmit the job and monitor its progress:
The job typically takes 30–60 minutes on 16 cores, depending on cluster load.
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 |
64 |
Maximum memory (in GB) that Cell Ranger is allowed to use. |
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:
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=workshop
ml --force purge
ml biocontainers cellranger
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=64For 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:
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.
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
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
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.
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.
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=workshop
#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
EOFThis step requires approximately 32 GB of RAM and takes about 30
minutes on 16 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=workshop
#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_
EOFSubmit the job:
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:
OUTPUT
barcodes.tsv features.tsv matrix.mtxThe 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.
Runtime comparison
On 64 cores with the PBMC 10k v3 dataset, typical runtimes are:
- Cell Ranger: ~7 hours
- STARsolo: ~45 minutes
STARsolo is roughly 9x 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 | 923.7M reads/hr | ~15x faster than Cell Ranger for this run |
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
Now that both tools have finished, let’s compare their outputs. 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 -lSTARsolo:
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.tsvBoth tools should report very similar numbers:
OUTPUT
Cell Ranger cells:
11809
Cell Ranger genes:
38606
STARsolo cells:
11682
STARsolo genes:
38606The 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.htmlreport - 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 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.)
Three changes are required for v2 chemistry:
--soloUMIlen 10instead of12. v2 chemistry uses 10 bp UMIs (compared to 12 bp in v3).--soloCBwhitelistmust 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 16stays 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 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% |
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.
- 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 approximately 9x faster and is fully open source
- 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-09-08 | Edit this page
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
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 in the previous episode. This
matrix lives in the 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.
First, set up the path to the Cell Ranger output so that all file references in this episode use the same base directory:
R
data_dir <- paste0(
"/scratch/negishi/", Sys.getenv("USER"),
"/scrna_workshop/cellranger_output/pbmc10k/outs/filtered_feature_bc_matrix/"
)
work_dir <- paste0(
"/scratch/negishi/", Sys.getenv("USER"),
"/scrna_workshop/"
)
Now load the count matrix:
R
setwd(work_dir)
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)
OUTPUT
[1] 38606 11809This 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 = 3removes genes detected in fewer than 3 cells (these are too rare to be informative) -
min.features = 200removes 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
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: countsAfter 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.
You can inspect the metadata that Seurat automatically computed during object creation:
R
head(pbmc@meta.data)
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 3365Seurat 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
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)
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")
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")
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.
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.
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)
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.
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")
OUTPUT
Cells before filtering: 11721
Cells after filtering: 11310Filtering 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")
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")
Finally, save the filtered object for the next episode:
R
saveRDS(pbmc, file = "pbmc_filtered.rds")
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
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, the 10x Chromium platform at standard loading (~10,000 cells targeted) produces a doublet rate of approximately 3–5%. At this level, most doublets will be removed by our QC filters or will form small insignificant clusters that can be identified and removed during cell type annotation. We will not run a dedicated doublet detection tool in this workshop, but be aware of these tools for your own analyses.
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 default thresholds
(nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 15)?
Run the following code with your chosen values:
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 9,874 with the default thresholds – a loss of roughly 23% 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. 5000): 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 9,874 to 7,579 cells means losing nearly a quarter 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 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
# Starting from the unfiltered object (re-load if needed)
pbmc_strict <- subset(pbmc,
subset = nFeature_RNA > 200 &
nFeature_RNA < 5000 &
percent.mt < 5)
cat("Cells with percent.mt < 5:", ncol(pbmc_strict), "\n")
pbmc_default <- subset(pbmc,
subset = nFeature_RNA > 200 &
nFeature_RNA < 5000 &
percent.mt < 15)
cat("Cells with percent.mt < 15:", ncol(pbmc_default), "\n")
With percent.mt < 5, only 2,709 cells
remain, compared to 9,874 with percent.mt < 15.
That is a loss of over 72% 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
< 5cutoff 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. Activated T cells, monocytes, and NK cells tend to have higher
mitochondrial content (6-12%) because they are metabolically active. A
< 5threshold disproportionately removes these populations, biasing your downstream analysis toward cell types with low metabolic activity. - A threshold of
< 5might 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.
-
Read10X()loads a Cell Ranger count matrix andCreateSeuratObject()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-09-08 | Edit this page
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
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)
R
pbmc <- readRDS(paste0(work_dir, "pbmc_filtered.rds"))
pbmc
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: countsThe 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)
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
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, dataNotice 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]
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, ]
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
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. |
Let’s see which genes were selected:
R
head(VariableFeatures(pbmc), 10)
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)
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.
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
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.dataAll 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
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
SCTrather than modifying theRNAassay. TheRNAassay 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 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?
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")
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.
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 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)
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.
- 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-09-08 | Edit this page
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
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)
R
pbmc <- readRDS(paste0(work_dir, "pbmc_normalized.rds"))
pbmc
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.dataPrincipal 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
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")
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.
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)
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.
How many PCs?
For the PBMC 10k dataset, the elbow typically falls around PC 12–15. We will use 15 PCs for the rest of this episode.
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")
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).
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.
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). |
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 secondsThe cluster assignments are stored in the object metadata. Let’s check them:
R
head(Idents(pbmc))
table(Idents(pbmc))
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 (likely a rare population like platelets or dendritic cells).
Now let’s visualize the clusters on the UMAP:
R
DimPlot(pbmc, reduction = "umap", label = TRUE, label.size = 5) +
NoLegend()
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.")
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)
| 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)
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 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)
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 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()
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, 2, 3, 5, 6 (several groups across the upper-left and lower-left UMAP). CD3D is expressed in both CD4+ and CD8+ T cells, so it marks multiple clusters. |
| MS4A1 | B cells | Cluster 7 (a compact, clearly separated group in the upper-right). MS4A1 (also known as CD20) is highly specific to B cells. |
| LYZ | CD14+ Monocytes | Clusters 1 and 2 (large groups in the lower portion of the UMAP). You may also see weaker LYZ expression in a smaller monocyte cluster (FCGR3A+). |
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.
- 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-09-08 | Edit this page
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
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)
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)
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()
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.
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 17The result is a data frame with one row per gene per cluster. Let’s look at the key columns:
R
head(pbmc.markers)
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 BCL11BThe 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.
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
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()
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))
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% |
| Platelets | PPBP, PF4 | ~2% |
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)
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)
Based on these plots and the FindAllMarkers results, we
can build a mapping from cluster numbers to cell type names. Your exact
cluster numbers will depend on the random seed, but the logic is the
same: find which clusters express each set of markers.
R
# Build the cluster-to-cell-type mapping
# Adjust cluster IDs to match YOUR results
new.cluster.ids <- c(
"CD4 T", # 0: CD3D+, IL7R+, large central cluster
"CD14 Mono", # 1: LYZ+ high, large bottom cluster
"CD14 Mono", # 2: LYZ+ high, adjacent to cluster 1
"CD4 T", # 3: CD3D+, IL7R+, upper-left
"CD8 T", # 4: CD8A+, right side
"CD4 T", # 5: CD3D+, IL7R+, upper-center
"NK", # 6: GNLY+ high, center-right
"NK", # 7: GNLY+, upper-right island
"B", # 8: MS4A1+, right side
"FCGR3A+ Mono", # 9: FCGR3A+ high, lower-right
"CD4 T", # 10: CD3D+, IL7R+, small upper cluster
"CD14 Mono", # 11: LYZ+, between mono clusters
"pDC", # 12: small isolated, low expression of all major markers
"DC", # 13: FCER1A+, small upper-right
"B", # 14: MS4A1+, adjacent to cluster 8
"Platelet", # 15: PPBP+ high, small isolated bottom
"CD14 Mono", # 16: LYZ+, near mono territory
"CD4 T" # 17: CD3D+, small isolated far-left
)
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
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()
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.
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
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.ontThe 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()
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)
OUTPUT
SingleR
Manual B cells Basophils CD4+ T cells CD8+ T cells Dendritic cells Monocytes NK cells Progenitors
CD4 T 3 0 2596 583 0 5 43 22
CD14 Mono 3 0 0 0 32 3259 1 0
CD8 T 1006 0 0 0 0 0 0 0
NK 0 0 3 1 0 0 607 0
B 608 0 0 0 5 0 15 0
FCGR3A+ Mono 0 0 0 0 1 428 0 0
pDC 0 0 0 0 157 19 0 0
DC 4 2 3 2 1 6 0 154
Platelet 0 0 0 0 78 0 0 0
SingleR
Manual T cells
CD4 T 987
CD14 Mono 0
CD8 T 0
NK 662
B 11
FCGR3A+ Mono 0
pDC 0
DC 3
Platelet 0The table reveals a mix of strong agreement and instructive disagreements:
- Strong agreement: CD14+ Monocytes are overwhelmingly called “Monocytes” by SingleR (3,259 of 3,295). B cells agree well (608 of 639). FCGR3A+ Monocytes are correctly called “Monocytes” (428 of 429). These cell types have distinctive transcriptional profiles that both methods recognize.
- CD4 T cell splitting: SingleR splits our CD4 T cells across “CD4+ T cells” (2,596), generic “T cells” (987), and “CD8+ T cells” (583). The generic “T cells” label reflects SingleR’s uncertainty about subtype when activation or memory signatures blur the CD4/CD8 boundary. The 583 cells called “CD8+ T” may be misclassified by SingleR or represent a subset where our manual annotation was too broad.
- CD8 T mislabeled as B cells: this is the largest disagreement. All 1,006 cells we labeled “CD8 T” are called “B cells” by SingleR. This likely means our manual annotation of cluster 4 as CD8 T is wrong and should be reconsidered – check the MS4A1 and CD8A expression in that cluster (visible in the VlnPlot) to determine whether SingleR’s B cell call is more appropriate. Go back and verify marker expression in this cluster.
- NK split with T cells: SingleR calls 607 of our NK cells as “NK cells” but labels 662 as “T cells.” This likely reflects NKT cells or cytotoxic CD8+ T cells that share expression of GNLY and NKG7 with NK cells. Our manual annotation grouped them based on GNLY expression alone, which is insufficient to distinguish NK from NKT populations.
- DC and pDC: SingleR calls our pDC cluster “Dendritic cells” (157), which is correct at a coarser level. However, it calls our DC cluster “Progenitors” (154), suggesting the reference dataset’s progenitor signature overlaps with the markers in that small cluster. This is a case where manual annotation with known markers (FCER1A, CST3) is more trustworthy.
- Platelets misclassified: SingleR calls all 78 platelets “Dendritic cells.” Platelets are often absent from SingleR reference datasets (Monaco Immune Data does not include them), so the classifier assigns the nearest available label. This is a known limitation – always verify rare or non-standard cell types manually.
The key takeaway: neither method is always right. SingleR caught a likely error in our manual CD8 T annotation but failed on platelets and DCs. The best practice is to 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
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
Guidelines for confident annotation
-
Use multiple references. No single reference is perfect. The
celldexpackage 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 1: Identifying T cell subtypes from marker combinations
The CD4 T cell group shows high expression of CD3D, SELL, and IL7R. But look at the CD69 panel – the violin is wide, not uniformly low or high. What does this heterogeneity tell you about the composition of this group?
R
VlnPlot(pbmc, features = c("CD3D", "SELL", "CCR7", "CD69", "IL7R"), ncol = 5, pt.size = 0)
The broad CD69 violin in CD4 T cells reveals that this group is not a single subtype – it contains a mixture of naive and activated T cells that were merged during our annotation.
| Marker | CD4 T pattern | Interpretation |
|---|---|---|
| CD3D | Uniformly high | Confirms T cell identity across the group |
| SELL (CD62L) | High | Lymph node homing receptor; enriched in naive and central memory T cells |
| CCR7 | Moderate, variable | Another homing receptor; high in naive, low in effector memory |
| CD69 | Broad (bimodal) | Early activation marker; low in naive, high in recently activated cells |
| IL7R (CD127) | High | IL-7 receptor for homeostatic survival; expressed on naive and memory T cells |
The CD69-low subset represents naive T cells (SELL-high, CCR7+, CD69-low) that recirculate between blood and lymph nodes but have not yet encountered antigen. The CD69-high subset represents recently activated T cells that have been stimulated.
Also note how these markers behave in other cell types:
- NK cells show the highest CD69 expression of any group, consistent with constitutive activation in circulating NK cells.
- B cells express SELL and IL7R, so these markers alone cannot distinguish T cells from B cells – you need CD3D to confirm T cell identity.
- CD8 T cells show high CD69 but low SELL, suggesting they are predominantly effector or effector memory cells rather than naive.
This illustrates why our coarse “CD4 T” label is a simplification. With higher clustering resolution or sub-clustering, you could separate the naive and activated 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 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)
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 1639The 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), “CMP” (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.” These are the cells we manually annotated as platelets. Monaco does not have a platelet reference profile, so it assigns them to “Progenitors” (the nearest match based on low transcriptional complexity). 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 should not be present in PBMCs. These 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.
- 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-09-08 | Edit this page
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
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
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, dataThe 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)
CTRL STIM
6548 7451Roughly 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
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.STIMThe 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")
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
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")
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")
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"]])
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).
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 annotate cell types and compare between conditions.
Annotating integrated clusters
Find markers for each integrated cluster:
R
ifnb.markers <- FindAllMarkers(ifnb,
only.pos = TRUE,
min.pct = 0.25,
logfc.threshold = 0.25)
Annotate the clusters using known PBMC markers (the same markers from Episode 6):
R
# Adjust cluster IDs to match YOUR results
new.cluster.ids <- c(
"CD14 Mono", # 0
"CD4 T", # 1
"CD4 T", # 2
"CD4 T", # 3
"CD8 T", # 4
"NK", # 5
"B", # 6
"CD14 Mono", # 7
"FCGR3A Mono", # 8
"DC", # 9
"CD14 Mono", # 10
"Mk", # 11
"B", # 12
"Eryth" # 13
)
names(new.cluster.ids) <- levels(ifnb)
ifnb <- RenameIdents(ifnb, new.cluster.ids)
ifnb$celltype <- Idents(ifnb)
R
DimPlot(ifnb, reduction = "umap", label = TRUE, repel = TRUE) + NoLegend()
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 LYZ, S100A9, and CD14
should appear as strong conserved markers for monocytes.
CTRL_p_val CTRL_avg_log2FC CTRL_pct.1 CTRL_pct.2 CTRL_p_val_adj STIM_p_val STIM_avg_log2FC
TYROBP 0 2.356269 0.940 0.227 0 0.000000e+00 2.337569
FCER1G 0 2.197678 0.930 0.233 0 0.000000e+00 2.549589
S100A8 0 4.739502 0.743 0.071 0 0.000000e+00 5.318497
C15orf48 0 2.720354 0.869 0.201 0 0.000000e+00 2.751572
IL8 0 3.651882 0.799 0.144 0 9.967253e-252 3.819957
CD63 0 2.740801 0.958 0.333 0 0.000000e+00 2.889984
S100A9 0 4.606045 0.691 0.069 0 0.000000e+00 5.236085
CTSB 0 3.116819 0.772 0.183 0 0.000000e+00 3.400786
LGALS1 0 2.400262 0.905 0.321 0 0.000000e+00 2.848345
TYMP 0 2.084914 0.860 0.282 0 0.000000e+00 2.267578
STIM_pct.1 STIM_pct.2 STIM_p_val_adj max_pval minimump_p_val
TYROBP 0.952 0.199 0.000000e+00 0.000000e+00 0
FCER1G 0.910 0.197 0.000000e+00 0.000000e+00 0
S100A8 0.445 0.020 0.000000e+00 0.000000e+00 0
C15orf48 0.860 0.201 0.000000e+00 0.000000e+00 0
IL8 0.295 0.028 1.400698e-247 9.967253e-252 0
CD63 0.931 0.276 0.000000e+00 0.000000e+00 0
S100A9 0.584 0.036 0.000000e+00 0.000000e+00 0
CTSB 0.822 0.183 0.000000e+00 0.000000e+00 0
LGALS1 0.867 0.217 0.000000e+00 0.000000e+00 0
TYMP 0.926 0.438 0.000000e+00 0.000000e+00 0Condition-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)
p_val avg_log2FC pct.1 pct.2 p_val_adj
IFIT1 0 7.117626 0.965 0.034 0
IFIT3 0 6.737338 0.971 0.057 0
TNFSF10 0 6.322227 0.959 0.062 0
RSAD2 0 6.619577 0.920 0.042 0
IFIT2 0 6.893956 0.909 0.040 0
MX1 0 4.824046 0.953 0.099 0
CXCL10 0 7.975000 0.873 0.032 0
LY6E 0 4.255152 0.990 0.174 0
CXCL11 0 8.521702 0.800 0.011 0
CCL8 0 9.063100 0.800 0.015 0The top DE genes are all interferon-stimulated genes (ISGs): IFIT1, IFIT3, TNFSF10, RSAD2, IFIT2. 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.
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"))
R
VlnPlot(ifnb,
features = c("IFIT1", "IFIT3", "TNFSF10", "RSAD2"),
split.by = "stim",
idents = "CD14 Mono",
ncol = 4,
pt.size = 0)
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 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?
CellType
Condition CD14 Mono CD4 T CD8 T NK B FCGR3A Mono DC Mk Eryth
CTRL 0.504 0.432 0.477 0.435 0.482 0.530 0.447 0.446 0.458
STIM 0.496 0.568 0.523 0.565 0.518 0.470 0.553 0.554 0.542
Most cell types fall in the 0.43–0.57 range, confirming that integration successfully mixed cells from both conditions. The split UMAP shows the same cell types present in both panels with similar spatial arrangement – each cluster appears in both CTRL and STIM with comparable density.
Some observations worth discussing:
- CD14 Mono is the best-mixed cell type (50.4/49.6), essentially a perfect split. This is expected for the most abundant population.
- CD4 T and NK show the largest deviation (~43/57), with more cells in the stimulated condition. This could reflect a real biological shift – IFN-beta stimulation may promote T cell and NK cell survival or proliferation in culture – rather than a failure of integration.
- Rare populations (Mk, Eryth, DC) show reasonable balance despite small cell numbers, where stochastic variation has more impact.
The key diagnostic: if integration had failed, you would see a cell type that is 90%+ from one condition, or you would see the same cell type forming separate clusters in the two panels. Neither is the case here. The moderate CD4 T skew (43/57) is well within the range of expected biological variation between a control and stimulated sample.
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"))
p_val avg_log2FC pct.1 pct.2 p_val_adj
CCL8 0.000000e+00 9.063 0.800 0.015 0.000000e+00
CXCL11 0.000000e+00 8.522 0.800 0.011 0.000000e+00
CXCL10 0.000000e+00 7.975 0.873 0.032 0.000000e+00
HESX1 1.756212e-217 7.935 0.335 0.002 2.468004e-213
IFIT1 0.000000e+00 7.118 0.965 0.034 0.000000e+00
IFIT2 0.000000e+00 6.894 0.909 0.040 0.000000e+00
IFIT3 0.000000e+00 6.737 0.971 0.057 0.000000e+00
GMPR 1.772282e-239 6.634 0.371 0.005 2.490588e-235
RSAD2 0.000000e+00 6.620 0.920 0.042 0.000000e+00
APOBEC3B 1.261976e-127 6.600 0.273 0.034 1.773454e-123The top 10 genes fall into two functional categories, both consistent with IFN-beta stimulation:
| Gene | log2FC | pct.1 / pct.2 | Function |
|---|---|---|---|
| CCL8 | 9.06 | 0.80 / 0.02 | Chemokine; recruits monocytes and T cells to sites of inflammation |
| CXCL11 | 8.52 | 0.80 / 0.01 | IFN-inducible chemokine; attracts activated T cells via CXCR3 |
| CXCL10 | 7.98 | 0.87 / 0.03 | IFN-inducible chemokine (IP-10); major recruiter of immune cells |
| HESX1 | 7.93 | 0.34 / 0.00 | Transcription factor; emerging role in IFN signaling |
| IFIT1 | 7.12 | 0.97 / 0.03 | Directly induced by type I IFN; inhibits viral translation |
| IFIT2 | 6.89 | 0.91 / 0.04 | Antiviral effector; regulates translation during infection |
| IFIT3 | 6.74 | 0.97 / 0.06 | Antiviral effector; forms complex with IFIT1 and IFIT2 |
| GMPR | 6.63 | 0.37 / 0.01 | Guanosine monophosphate reductase; modulates purine metabolism |
| RSAD2 | 6.62 | 0.92 / 0.04 | Viperin; broad-spectrum antiviral enzyme induced by IFN |
| APOBEC3B | 6.60 | 0.27 / 0.03 | Cytidine deaminase; innate antiviral defense via RNA editing |
Two patterns stand out:
Chemokines dominate the top 3. CCL8, CXCL11, and CXCL10 have the highest fold changes (7.9-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.03) – these chemokines are essentially absent in control cells and massively induced by stimulation.
Classical ISGs fill the rest. IFIT1/2/3, RSAD2, and APOBEC3B are canonical interferon-stimulated genes that directly inhibit viral replication. The IFIT family shows the highest pct.1 values (0.91-0.97), meaning nearly every stimulated monocyte expresses them – a hallmark of a robust, uniform interferon response.
The FeaturePlot confirms this visually: CCL8, CXCL11, CXCL10, and HESX1 are virtually absent in the CTRL panels (gray) and strongly expressed in the STIM panels (red), with the signal concentrated in the monocyte clusters (right side of the UMAP). CXCL10 and CXCL11 also show some expression in the stimulated CD4 T and NK clusters, consistent with these cell types also responding to IFN-beta, though less strongly than monocytes.
- 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-09-08 | Edit this page
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
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)
R
ifnb <- readRDS(paste0(work_dir, "ifnb_annotated.rds"))
Idents(ifnb) <- "celltype"
DimPlot(ifnb, reduction = "umap", label = TRUE, repel = TRUE, split.by = "stim") +
NoLegend()
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")
OUTPUT
p_val avg_log2FC pct.1 pct.2 p_val_adj
IFIT1 0 7.118 0.965 0.034 0
IFIT3 0 6.737 0.971 0.057 0
TNFSF10 0 6.322 0.959 0.062 0
RSAD2 0 6.620 0.920 0.042 0
IFIT2 0 6.894 0.909 0.040 0
MX1 0 4.824 0.953 0.099 0
CXCL10 0 7.975 0.873 0.032 0
LY6E 0 4.255 0.990 0.174 0
CXCL11 0 8.522 0.800 0.011 0
CCL8 0 9.063 0.800 0.015 0
Significant genes (padj < 0.05): 765 You should see approximately 765 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.”
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.
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.
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 2453Now 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. |
OUTPUT
An object of class Seurat
14053 features across 54 samples within 1 assay
Active assay: RNA (14053 features, 0 variable features)
3 layers present: counts, data, scale.dataThe number of samples equals the number of unique cell type x condition x pseudo-replicate combinations. With 9 cell types, 2 conditions, and 3 pseudo-replicates, you get up to 54 pseudobulk samples (some may be fewer if a cell type has very few 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")
OUTPUT
p_val avg_log2FC pct.1 pct.2 p_val_adj
ISG15 0 6.441088 1 1 0
IFI6 0 3.705585 1 1 0
IFI44L 0 3.374846 1 1 0
GBP1 0 2.535916 1 1 0
S100A8 0 -2.276587 1 1 0
MNDA 0 2.661545 1 1 0
RSAD2 0 5.839812 1 1 0
TMSB10 0 1.655318 1 1 0
VAMP5 0 2.977336 1 1 0
SSB 0 2.470712 1 1 0
Significant genes (p_val_adj < 0.05): 1455 The pseudobulk analysis identifies 1455 significant genes – more than the Wilcoxon test (765). 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()
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"))
R
VlnPlot(ifnb,
features = c("IFIT1", "CXCL10", "MX1"),
split.by = "stim",
pt.size = 0,
ncol = 3)
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))
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.
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). |
OUTPUT
Number of significant upregulated genes: 271
SYMBOL ENTREZID
1 IFIT1 3434
2 IFIT3 3437
3 TNFSF10 8743
4 RSAD2 91543
5 IFIT2 3433
6 MX1 4599Running 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")
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)
OUTPUT
category subcategory ID Description GeneRatio
hsa04621 Organismal Systems Immune system hsa04621 NOD-like receptor signaling pathway 22/166
hsa05168 Human Diseases Infectious disease: viral hsa05168 Herpes simplex virus 1 infection 17/166
hsa05169 Human Diseases Infectious disease: viral hsa05169 Epstein-Barr virus infection 18/166
hsa05164 Human Diseases Infectious disease: viral hsa05164 Influenza A 16/166
hsa04612 Organismal Systems Immune system hsa04612 Antigen processing and presentation 11/166
BgRatio RichFactor FoldEnrichment zScore pvalue p.adjust qvalue
hsa04621 187/9382 0.11764706 6.649185 10.472204 1.358486e-12 3.124518e-10 2.702672e-10
hsa05168 180/9382 0.09444444 5.337818 7.886310 1.677668e-08 1.503560e-06 1.300562e-06
hsa05169 205/9382 0.08780488 4.962562 7.698570 1.961165e-08 1.503560e-06 1.300562e-06
hsa05164 173/9382 0.09248555 5.227105 7.531253 6.092210e-08 3.503021e-06 3.030073e-06
hsa04612 82/9382 0.13414634 7.581693 8.033604 1.806443e-07 8.309640e-06 7.187743e-06
geneID
hsa04621 4938/3665/115361/2633/4939/4940/6347/3428/115362/10616/7295/3320/6772/837/2635/79792/114769/3326/2634/834/10010/10628
hsa05168 4938/3665/5610/64135/684/4939/4940/6347/3133/6890/6772/3107/8717/6672/29992/6892/3134
hsa05169 3627/4938/9636/3665/5610/4939/4940/3133/6890/6772/3107/4067/958/8717/9541/4616/6892/3134
hsa05164 8743/91543/4599/3627/4938/4600/3665/5610/64135/4939/4940/6347/6772/8717/103/834
hsa04612 3303/3133/5721/6890/3320/3107/1520/3326/1514/6892/3134
Count
hsa04621 22
hsa05168 17
hsa05169 18
hsa05164 16
hsa04612 11Top KEGG pathways should include entries like “Influenza A”, “Hepatitis C”, “RIG-I-like receptor signaling pathway”, and “NOD-like receptor signaling pathway” – all pathways involving the interferon/antiviral response.
Exporting Results
At the end of an analysis, you should export your results in formats that support reproducibility, sharing, and downstream use.
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:
R
# Export raw count matrix (warning: can be large for big datasets)
counts_matrix <- GetAssayData(pbmc, layer = "counts")
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)
For smaller datasets, you can export as a dense CSV:
R
write.csv(as.matrix(GetAssayData(pbmc, layer = "counts")),
file = "counts_dense.csv")
Format conversion for Python users
If collaborators use the Python/Scanpy ecosystem, you can convert
your Seurat object to AnnData format using the sceasy
package:
R
# Install sceasy if needed: remotes::install_github("cellgeni/sceasy")
library(sceasy)
sceasy::convertFormat(pbmc,
from = "seurat",
to = "anndata",
outFile = "pbmc_annotated.h5ad")
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)
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
renvto 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 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?
You should find that:
- Wilcoxon identifies 765 significant genes
- Pseudobulk identifies 1455 (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 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")
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 | viral process | Genes involved in the viral life cycle, many of which are ISGs that restrict it |
| 4 | viral life cycle | Overlap with viral process; reflects IFIT, OAS, and MX family genes |
| 5 | regulation of viral process | Regulatory genes that modulate antiviral defense (e.g., TRIM, IRF family) |
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.
- 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