9. Single-Cell & Spatial Transcriptomics

Cell Ranger is 10x Genomics’ official processing pipeline for Chromium Single Cell gene expression. It processes raw BCL or FASTQ files, demultiplexes cellular barcodes, maps reads to transcripts, deduplicates UMIs, and generates feature-barcode matrices.

# 1. Run Cell Ranger count on 10x 3' Gene Expression libraries
cellranger count \
  --id=pbmc_10k_sample \
  --transcriptome=/ref/refdata-gex-GRCh38-2024-A \
  --fastqs=/data/pbmc_10k_fastqs/ \
  --sample=pbmc_10k \
  --expect-cells=10000 \
  --localcores=16 \
  --localmem=64

# Filtered gene-barcode matrix:
# pbmc_10k_sample/outs/filtered_feature_bc_matrix.h5
# Interactive web report: pbmc_10k_sample/outs/web_summary.html

STARsolo is built into STAR for mapping single-cell RNA-seq. It emulates 10x Genomics Cell Ranger barcode and UMI algorithms, delivering identical cell-by-gene matrices with up to 10x lower runtime and flexible barcode geometries.

# 1. Run STARsolo on 10x Chromium v3 data
STAR --runThreadN 16 \
  --genomeDir ./star_index/ \
  --readFilesIn sample_R2.fq.gz sample_R1.fq.gz \
  --readFilesCommand zcat \
  --soloType CB_UMI_Simple \
  --soloCBwhitelist 3M-february-2018.txt \
  --soloCBlen 16 --soloUMIlen 12 --soloUMIstart 17 \
  --soloCellFilter EmptyDrops_CR \
  --soloFeatures Gene GeneFull \
  --outSAMtype None \
  --outFileNamePrefix ./starsolo_out/

Alevin and Alevin-fry provide an ultra-fast, memory-lean suite for processing single-cell and single-nucleus RNA-seq data, resolving ambiguous multimapping reads with parsimonious UMI graphs.

# 1. Generate RAD file with Salmon Alevin
salmon alevin -l ISR -1 sample_R1.fq.gz -2 sample_R2.fq.gz \
  --chromiumV3 -i salmon_splici_index \
  -p 16 -o alevin_rad_out --rad

# 2. Permit list generation and collating with alevin-fry
alevin-fry generate-permit-list -d fw -k -i alevin_rad_out -o alevin_fry_work
alevin-fry collate -r alevin_rad_out -t 16 -i alevin_fry_work
alevin-fry quant -m cr-like -i alevin_fry_work -o alevin_counts -t 16 --use-unspliced

Seurat is an R package designed for QC, analysis, and exploration of single-cell RNA-seq and spatial data. It features SCTransform variance-stabilizing normalization, PCA/UMAP embeddings, Louvain/Leiden clustering, and multi-modal CITE-seq integration.

library(Seurat)
library(dplyr)

# 1. Load 10x Genomics matrix and filter low-quality cells
pbmc.data <- Read10X(data.dir = "filtered_feature_bc_matrix/")
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "PBMC", min.cells = 3, min.features = 200)
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 10)

# 2. Run SCTransform normalization
pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt", verbose = FALSE)

# 3. Dimensionality reduction and Clustering
pbmc <- RunPCA(pbmc, verbose = FALSE)
pbmc <- FindNeighbors(pbmc, dims = 1:30)
pbmc <- FindClusters(pbmc, resolution = 0.6)
pbmc <- RunUMAP(pbmc, dims = 1:30)

# 4. Find cluster biomarker markers
markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
head(markers)

Scanpy is a scalable toolkit for analyzing single-cell gene expression data built on AnnData and NumPy/SciPy. It executes preprocessing, highly variable gene selection, PCA, k-NN graph construction, UMAP embeddings, and Leiden community detection across millions of cells.

import scanpy as sc

# 1. Ingest 10x HDF5 matrix
adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")
adata.var_names_make_unique()

# 2. Calculate QC metrics and filter
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)
adata = adata[(adata.obs.n_genes_by_counts < 2500) & (adata.obs.pct_counts_mt < 10)]

# 3. Normalize and identify highly variable genes
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
adata.raw = adata
adata = adata[:, adata.var.highly_variable]

# 4. Dimensionality reduction & Leiden clustering
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, svd_solver="arpack")
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.5)

# 5. Plot and save
sc.pl.umap(adata, color=["leiden"], save="_clusters.png")

scvi-tools is a suite of deep probabilistic models for single-cell omics data. scVI utilizes Variational Autoencoders (VAEs) to learn a non-linear latent space that captures cell identity while explicitly removing batch effects and library size differences.

import scvi
import scanpy as sc

# 1. Setup AnnData with raw counts and batch key
adata = sc.read_h5ad("combined_cohort.h5ad")
scvi.model.SCVI.setup_anndata(adata, layer="counts", batch_key="donor_batch")

# 2. Initialize and train Variational Autoencoder
vae = scvi.model.SCVI(adata, n_latent=30, gene_likelihood="nb")
vae.train(max_epochs=200, early_stopping=True)

# 3. Extract integrated latent representation
adata.obsm["X_scVI"] = vae.get_latent_representation()

# 4. Cluster on integrated latent space
sc.pp.neighbors(adata, use_rep="X_scVI")
sc.tl.umap(adata)
sc.pl.umap(adata, color=["donor_batch", "cell_type"])

Harmony aligns single-cell transcriptomics data from multiple donors, experimental batches, and technologies. By operating iteratively on low-dimensional PCA space, it projects cells into a shared embedding while preserving real biological variation.

library(harmony)
library(Seurat)

# 1. Merge Seurat objects across batches and run PCA
merged_seurat <- NormalizeData(merged_seurat) %>%
  FindVariableFeatures() %>%
  ScaleData() %>%
  RunPCA(verbose = FALSE)

# 2. Run Harmony integration across donor batches
merged_seurat <- RunHarmony(merged_seurat, group.by.vars = "batch", dims.use = 1:30)

# 3. Build UMAP on Harmony embeddings
merged_seurat <- RunUMAP(merged_seurat, reduction = "harmony", dims = 1:30)
DimPlot(merged_seurat, reduction = "umap", group.by = "batch")

Monocle 3 discovers dynamic biological transitions (stem cell differentiation, immune activation) by learning principal graph trajectories through single-cell transcriptomes and calculating pseudotime along developmental branches.

library(monocle3)

# 1. Construct cell_data_set from Seurat or count matrix
cds <- new_cell_data_set(expression_data, cell_metadata = col_data, gene_metadata = row_data)
cds <- preprocess_cds(cds, num_dim = 50)
cds <- reduce_dimension(cds)

# 2. Cluster cells and learn principal trajectory graph
cds <- cluster_cells(cds)
cds <- learn_graph(cds, use_partition = FALSE)

# 3. Order cells by pseudotime from root node
cds <- order_cells(cds)
plot_cells(cds, color_cells_by = "pseudotime", label_cell_groups = FALSE)

scVelo solves the full transcriptional dynamics of splicing kinetics using a likelihood-based dynamical model. By comparing unspliced pre-mRNA with mature spliced mRNA, it projects directional cell differentiation vectors onto UMAP embeddings.

import scvelo as scv

# 1. Load loom file containing spliced and unspliced layers
adata = scv.read("sample_velocyto.loom", cache=True)

# 2. Preprocess and compute moments
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

# 3. Recover dynamics and compute velocity graph
scv.tl.recover_dynamics(adata)
scv.tl.velocity(adata, mode="dynamical")
scv.tl.velocity_graph(adata)

# 4. Project velocity streams on UMAP
scv.pl.velocity_embedding_stream(adata, basis="umap", color="cell_type", save="velocity_stream.png")

Squidpy brings spatial analysis to AnnData. It constructs spatial cellular neighborhood graphs, computes Moran’s I spatial autocorrelation, tests cellular co-occurrence, and integrates histology images with gene expression in Visium and Xenium assays.

import squidpy as sq

# 1. Ingest 10x Visium spatial dataset
adata = sq.read.visium("visium_data_dir/")

# 2. Build spatial neighborhood graph from physical spot coordinates
sq.gr.spatial_neighbors(adata, coord_type="grid", n_rings=2)

# 3. Calculate spatial autocorrelation (Moran's I) to identify spatially variable genes
sq.gr.spatial_autocorr(adata, mode="moran", genes=adata.var_names[:1000], n_perms=100)

# 4. Compute cellular neighborhood enrichment across cell types
sq.gr.nhood_enrichment(adata, cluster_key="cluster")
sq.pl.nhood_enrichment(adata, cluster_key="cluster", save="nhood_enrichment.png")

CellChat quantitatively infers intercellular communication networks from single-cell gene expression data using a comprehensive database of human and mouse ligand-receptor interactions, co-factors, and multi-subunit complexes.

library(CellChat)

# 1. Initialize CellChat object with expression and cell cluster labels
cellchat <- createCellChat(object = normalized_counts, meta = meta_data, group.by = "cell_type")
CellChatDB <- CellChatDB.human
cellchat@DB <- CellChatDB

# 2. Subset to signaling genes and compute communication probabilities
cellchat <- subsetData(cellchat)
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
cellchat <- computeCommunProb(cellchat, raw.use = TRUE)
cellchat <- aggregateNet(cellchat)

# 3. Visualize intercellular communication network
netVisual_circle(cellchat@net$count, vertex.weight = groupSize, weight.scale = TRUE,
                  label.edge= FALSE, title.name = "Number of Interactions")