| name | bioconductor-scran |
| description | Implements miscellaneous functions for interpretation of single-cell RNA-seq data. Methods are provided for assignment of cell cycle phase, detection of highly variable and significantly correlated genes, identification of marker genes, and |
| when_to_use | Use when: Pooling-based normalisation (handles zero inflation); Detecting highly variable genes for HVG selection; Graph-based clustering of single cells; Marker gene detection with pairwise tests; Pseudobulk DE analysis with aggregateAcrossCells |
| user-invocable | false |
scran — Comprehensive Skill Guide
Domain: Single-cell RNA-seq
Bioconductor: scran
Paper: Lun ATL, McCarthy DJ, Marioni JC (2016). Genome Biology, 17:75.
Methods for single-cell RNA-seq: pooling-based normalisation, high-variable gene detection, doublet detection, and graph-based clustering.
Dependencies & Environment
Package-intrinsic requirements from the Bioconductor landing page — reproduce in any R environment.
- Version: 1.40.0 · Bioconductor: 3.23 · R: ≥ 4.6
- Depends: SingleCellExperiment, scuttle
- Imports: SummarizedExperiment, S4Vectors, BiocGenerics, BiocParallel, Rcpp, Matrix, edgeR, limma, igraph, statmod, MatrixGenerics, S4Arrays, DelayedArray, BiocSingular, bluster, metapod, dqrng, beachmat
- System requirements: C++11
- Install:
BiocManager::install("scran")
When to Use
- Pooling-based normalisation (handles zero inflation)
- Detecting highly variable genes for HVG selection
- Graph-based clustering of single cells
- Marker gene detection with pairwise tests
- Pseudobulk DE analysis with aggregateAcrossCells
Alternatives: Seurat normalisation, scuttle, sctransform
Do NOT Use When
- Bulk RNA-seq — use DESeq2 or edgeR.
- Very small datasets (<100 cells) where pooling-based size factor estimation fails.
Data Requirements
| Requirement | Details |
|---|
| Min Cells | ≥100 cells (pooling for size factors requires sufficient cells) |
| Input Format | SingleCellExperiment with raw count matrix |
| Notes | Works best after quality filtering with scater/perCellQCMetrics |
Installation
if (!require("BiocManager", quietly=TRUE))
install.packages("BiocManager")
BiocManager::install("scran")
library(scran)
Workflows
scRNA-seq Normalization and HVG Selection
Core OSCA workflow: normalize pooled size factors → model gene variance → select HVGs
Normalize with pooling-based size factors
library(scran)
library(scater)
set.seed(42)
clusters <- quickCluster(sce)
sce <- computeSumFactors(sce, clusters=clusters)
summary(sizeFactors(sce))
sce <- logNormCounts(sce)
Model gene variance and select highly variable genes
dec <- modelGeneVar(sce)
hvgs <- getTopHVGs(dec, n=3000)
dec1 <- modelGeneVar(sce, block=sce$batch)
hvgs_blocked <- getTopHVGs(dec1, n=3000)
Dimensionality reduction and clustering
sce <- runPCA(sce, subset_row=hvgs, ncomponents=30)
sce <- runUMAP(sce, dimred="PCA")
library(bluster)
clusters <- clusterRows(reducedDim(sce,"PCA"), NNGraphParam(k=20))
sce$cluster <- factor(clusters)
plotReducedDim(sce, "UMAP", colour_by="cluster")
Per-cluster marker gene detection
markers <- findMarkers(sce, groups=sce$cluster,
test.type="wilcox", direction="up")
markers[[1]][1:10, c("Top","p.value","FDR")]
Key Functions & Parameters
quickCluster()
Fast pre-clustering for normalisation (required before computeSumFactors)
| Parameter | Description |
|---|
x | SingleCellExperiment or matrix |
min.size | minimum cluster size (default 100) |
method | 'igraph' (default) | 'hclust' |
computeSumFactors()
Pooling-based normalisation (deconvolution); handles near-zero counts
| Parameter | Description |
|---|
x | SingleCellExperiment |
clusters | cluster assignments from quickCluster() |
min.mean | lower bound on average count for HEG selection (default 0.1) |
BPPARAM | parallel backend |
modelGeneVar()
Model per-gene variance to identify HVGs
| Parameter | Description |
|---|
x | log-normalised SCE |
block | factor for blocking by sample/batch |
density.weights | use density weights for mean-variance fit |
getTopHVGs()
Select top HVGs from modelGeneVar output
| Parameter | Description |
|---|
stats | output of modelGeneVar() |
n | number of HVGs (default 2000) |
var.threshold | minimum biological variance threshold |
fdr.threshold | FDR threshold for significant HVGs |
buildSNNGraph()
Build shared nearest-neighbour graph for clustering
| Parameter | Description |
|---|
x | SCE or matrix |
k | number of nearest neighbours (default 10; higher=coarser clusters) |
type | 'rank' (default) | 'number' | 'jaccard' |
d | number of PCs to use |
Scientific Assumptions
The method assumes the following about your data and experiment:
- Pooling reads across cells enables robust size factor estimation even for very sparse single-cell data.
- Most genes are NOT differentially expressed between cells of the same type (for size factor estimation).
- Highly variable genes (HVGs) reflect biological variation, not purely technical noise.
Result Interpretation
sizeFactors(sce): per-cell scaling factors; very high or very low values indicate problematic cells.logcounts: log1p-transformed normalized counts; use for visualization and clustering input.modelGeneVar() result: total variance = bio + tech; only bio > 0 genes are HVGs.pairwiseTTests() / scoreMarkers(): use Cohen's d and AUC for marker specificity, not just p-value.
Best Practices
- Run
quickCluster() before computeSumFactors() — pooling within clusters is critical for accuracy.
- Use
scuttle::logNormCounts() after computeSumFactors() to apply size factors.
- For multi-sample experiments, pass
block=sce$Sample to modelGeneVar() to account for batch.
- Prefer
aggregateAcrossCells() + DESeq2/edgeR for differential abundance (pseudobulk).
- Use
scDblFinder::scDblFinder() (uses scran infrastructure) for doublet detection.
When to Choose This vs Alternatives
| Alternative | Prefer this tool when | Prefer alternative when |
|---|
Seurat SCTransform | You are in the Seurat ecosystem and want regularized NB regression normalization. | You want Bioconductor-native normalization with theoretical underpinning for sparse data. |
Benchmark Evidence
Key studies that have validated or benchmarked this tool:
- Lun ATL, McCarthy DJ, Marioni JC (2016). Genome Biol 17:75. PMID:27122128
Pooling-based size factor estimation outperforms library-size normalization for sparse scRNA-seq data.
Common Errors & Troubleshooting
Error in computeSumFactors: not enough cells in cluster
Cause: Cluster too small for pooling
Fix: Reduce min.size in quickCluster or merge small clusters
Additional Notes from Official Documentation
Extracted from the scran Bioconductor vignette(s)
1 Introduction
Single-cell RNA sequencing (scRNA-seq) is a widely used technique for profiling gene expression in individual cells.
This allows molecular biology to be studied at a resolution that cannot be matched by bulk sequencing of cell populations.
The scran package implements methods to perform low-level processing of scRNA-seq data,
including cell cycle phase assignment, variance modelling and testing for marker genes and gene-gene correlations.
This vignette provides brief descriptions of these methods and some toy examples to demonstrate their use.
Note: A more comprehensive description of the use of scran (along with other packages) in a scRNA-seq analysis workflow is available at https://osca.bioconductor.org .
2 Setting up the data
We start off with a count matrix where each row is a gene and each column is a cell.
These can be obtained by mapping read sequences to a reference genome, and then counting the number of reads mapped to the exons of each gene.
(See, for example, the Rsubread package to do both of these tasks.)
Alternatively, pseudo-alignment methods can be used to quantify the abundance of each transcript in each cell.
For simplicity, we will pull out an existing dataset from the scRNAseq package.
This particular dataset is taken from a study of the human pancreas with the CEL-seq protocol (Grun et al. 2016 ) .
It is provided as a SingleCellExperiment object (from the SingleCellExperiment package), which contains the raw data and various annotations.
We perform some cursory quality control to remove cell
library(scRNAseq)
sce <- GrunPancreasData()
sce
3 Variance modelling
We identify genes that drive biological heterogeneity in the data set by modelling the per-gene variance.
By only using a subset of highly variable genes in downstream analyses like clustering, we improve resolution of biological structure by removing uninteresting genes driven by technical noise.
We decompose the total variance of each gene into its biological and technical components by fitting a trend to the endogenous variances (Lun, McCarthy, and Marioni 2016 ) .
The fitted value of the trend is used as an estimate of the technical component, and we subtract the fitted value from the total variance to obtain the biological component for each gene.
If we have spike-ins, we can use them to fit the trend instead.
This provides a more direct estimate of the technical variance and avoids
dec <- modelGeneVar(sce)
plot(dec$mean, dec$total, xlab="Mean log-expression", ylab="Variance")
curve(metadata(dec)$trend(x), col="blue", add=TRUE)
dec2 <- modelGeneVarWithSpikes(sce, 'ERCC')
plot(dec2$mean, dec2$total, xlab="Mean log-expression", ylab="Variance")
points(metadata(dec2)$mean, metadata(dec2)$var, col="red")
curve(metadata(dec2)$trend(x), col="blue", add=TRUE)
4 Automated PC choice
Principal components analysis is commonly performed to denoise and compact the data prior to downstream analysis.
A common question is how many PCs to retain; more PCs will capture more biological signal at the cost of retaining more noise and requiring more computational work.
One approach to choosing the number of PCs is to use the technical component estimates to determine the proportion of variance that should be retained.
This is implemented in denoisePCA() , which takes the estimates returned by modelGeneVar() or friends.
(For greater accuracy, we use the fit with the spikes; we also subset to only the top HVGs to remove noise.)
Another approach is based on the assumption that each subpopulation should be separated from each other on a different axis of variation.
Thus, we choose th
sced <- denoisePCA(sce, dec2, subset.row=getTopHVGs(dec2, prop=0.1))
ncol(reducedDim(sced, "PCA"))
5 Graph-based clustering
Clustering of scRNA-seq data is commonly performed with graph-based methods due to their relative scalability and robustness. scran provides several graph construction methods based on shared nearest neighbors (Xu and Su 2015 ) through the buildSNNGraph() function.
This is most commonly generated from the selected PCs, after which community detection methods from the igraph package can be used to explicitly identify clusters.
By default, buildSNNGraph() uses the mode of shared neighbor weighting described by Xu and Su ( 2015 ) , but other weighting methods (e.g., the Jaccard index) are also available by setting type= .
An unweighted (k) -nearest neighbor graph can also be constructed with buildKNNGraph() .
We can then use methods from scater to visualize this clustering on a (t) -SNE
g <- buildSNNGraph(sce, use.dimred="PCAsub")
cluster <- igraph::cluster_walktrap(g)$membership
colLabels(sce) <- factor(cluster)
table(colLabels(sce))
6 Identifying marker genes
The scoreMarkers() wrapper function will perform differential expression comparisons between pairs of clusters to identify potential marker genes.
For each pairwise comparison, we compute a variety of effect sizes to quantify the differences between those clusters.
Cohenâs (d) is a standardized log-fold change, representing the number of standard deviations that separate the means of two groups.
This is analogous to the (t) -statistic in Studentâs (t) -test.
The area under the curve (AUC) is the probability that a randomly chosen observation from one group is greater than a random observation from another group.
This is proportional to the U-statistic from the Wilcoxon ranked sum test.
We also compute the log-fold change in the proportion of cells with detectable (i.e., non-zero
markers <- scoreMarkers(sce)
markers
7 Detecting correlated genes
Another useful procedure is to identify significant pairwise correlations between pairs of HVGs.
The idea is to distinguish between HVGs caused by random stochasticity, and those that are driving systematic heterogeneity, e.g., between subpopulations.
Correlations are computed by the correlatePairs method using a slightly modified version of Spearmanâs rho,
tested against the null hypothesis of zero correlation using the same method in cor.test() .
As with variance estimation, if uninteresting substructure is present, this should be blocked on using the block= argument.
This avoids strong correlations due to the blocking factor.
The pairs can be used for choosing marker genes in experimental validation, and to construct gene-gene association networks.
In other situations, the pairs may
of.interest <- top.hvgs[1:200]
cor.pairs <- correlatePairs(sce, subset.row=of.interest)
cor.pairs
8 Converting to other formats
The SingleCellExperiment object can be easily converted into other formats using the convertTo method.
This allows analyses to be performed using other pipelines and packages.
For example, if DE analyses were to be performed using edgeR , the count data in sce could be used to construct a DGEList .
By default, rows corresponding to spike-in transcripts are dropped when get.spikes=FALSE .
As such, the rows of y may not correspond directly to the rows of sce â users should match by row name to ensure correct cross-referencing between objects.
Normalization factors are also automatically computed from the size factors.
The same conversion strategy roughly applies to the other supported formats.
DE analyses can be performed using DESeq2 by converting the object to a DESeqDataSet .
Cells can
y <- convertTo(sce, type="edgeR")
Resources
Run this on BioMate
This skill is the knowledge layer — when, why, and how to use scran. To run this analysis on your own data with managed compute, automated QC, and reproducible outputs, use BioMate — free to start.
▶ Open scran on BioMate →
