How to Cluster Your Own RNA-Seq Data

Oak
rna-seq clustering transcriptomics deseq2 scanpy
Glowing jellyfish gathered in separate colored shoals in a dark tidal basin beneath a branching coral tree of light.

By the end of this you will have two things: a sample-level structure map that tells you how your timepoints or tissues relate to each other, and a set of gene modules, groups of transcripts that move together across your samples, each annotated with the biological processes it is enriched for. You need a counts matrix (genes × samples, integers, not TPM), at least six samples for anything sample-level and ideally twelve or more for gene modules, R 4.3 with DESeq2, cluster, pheatmap, fpc, and fgsea, plus Python 3.11 with scanpy if you have single-cell data. If you only have raw FASTQ files, step 1 covers getting to counts. Everything here is measurement and interpretation. Clusters are hypotheses about your data, not findings about your health, and anything that looks clinically relevant belongs in front of a physician.

1. Get to a counts matrix you trust

Clustering amplifies quantification artifacts, so the upstream choices matter more than the clustering algorithm you pick later. We use Salmon for bulk RNA-seq because selective alignment against the transcriptome is fast, and the bias-correction flags materially change the expression estimates for GC-extreme transcripts.

salmon index -t gencode.v44.transcripts.fa.gz -i idx_v44 --gencode -k 31

salmon quant -i idx_v44 -l A \
  -1 T01_R1.fastq.gz -2 T01_R2.fastq.gz \
  --validateMappings --seqBias --gcBias --posBias \
  --numBootstraps 30 -p 8 -o quants/T01

Then collapse transcripts to genes with tximport, which carries the average transcript length per sample into an offset that DESeq2 uses. This matters when isoform usage shifts between timepoints, because gene-level counts alone would misattribute a length change to an expression change.

library(tximport); library(DESeq2)
tx2gene <- read.csv("tx2gene_gencode_v44.csv")   # tx_id, gene_id
files <- file.path("quants", samples$id, "quant.sf")
names(files) <- samples$id
txi <- tximport(files, type = "salmon", tx2gene = tx2gene)

dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ 1)

Use design = ~ 1 deliberately. For clustering you do not want the model to absorb the structure you are trying to see.

2. Normalize and transform, and understand why this decides your dendrogram

Raw counts are heteroscedastic: a gene with mean count 20,000 has far larger absolute variance than one with mean count 20, so Euclidean distance on untransformed counts is dominated by a few ribosomal and mitochondrial transcripts. Clustering on log2(count + 1) overcorrects in the other direction, inflating the apparent variance of low-count genes where the noise is Poisson. DESeq2’s variance-stabilizing transformation fits the mean-variance trend and applies a transformation that flattens it.

vsd <- vst(dds, blind = TRUE)     # blind = TRUE: ignore the design entirely
mat <- assay(vsd)                 # genes x samples, roughly homoscedastic

For fewer than about 30 samples, rlog(dds, blind = TRUE) is a reasonable alternative and shrinks low-count genes more aggressively, at a large cost in runtime. Either is better than log-CPM for distance-based methods. The general point is that normalization choice is not a preprocessing detail you can defer: systematic comparisons across omics types show that the normalization method changes downstream unsupervised clustering results substantially, and that no single method wins on every dataset1. Run your pipeline with two transformations and check whether the sample groupings survive. If they do not, you have learned something important about how weak the signal is.

3. Filter to genes that carry information

Half the human genome is not expressed in whole blood, and undetected genes contribute nothing but noise to a distance calculation.

keep <- rowSums(counts(dds) >= 10) >= 3      # 10+ reads in at least 3 samples
dds  <- dds[keep, ]

For sample-level clustering, restrict further to the most variable genes. The top 500 to 2,000 by variance across samples is the standard range, and the result is usually stable across that window. If the picture changes drastically between 500 and 2,000 genes, your sample structure is being driven by a small handful of transcripts, which is worth investigating directly.

v   <- rowVars(mat)
top <- order(v, decreasing = TRUE)[1:1000]
matTop <- mat[top, ]

4. Sample-level structure: PCA first, then hierarchical clustering

Start with PCA, because it answers the first question you have, which is whether there is any structure at all. PCA is a rotation into orthogonal axes ordered by variance explained. It does not assign labels, so it is not itself a clustering method, but it is the right first step and a good input to one: clustering on the top 10 to 20 principal components rather than on 1,000 genes removes a large fraction of the gene-level noise while preserving the covariance structure that matters.

pca <- prcomp(t(matTop), scale. = FALSE)
summary(pca)$importance[, 1:5]
plot(pca$x[, 1], pca$x[, 2], pch = 19)

Look at the variance explained by PC1 before anything else. If PC1 carries 60% or more and separates samples by collection date or library prep batch, stop and deal with that (see Common problems). If PC1 carries 15% and the samples are a diffuse cloud, you likely have a stable baseline with no dominant axis of variation, which is itself a meaningful result for a longitudinal personal profile.

For the dendrogram, use correlation distance rather than Euclidean, and Ward linkage:

d  <- as.dist(1 - cor(matTop, method = "spearman"))
hc <- hclust(d, method = "ward.D2")
plot(hc, hang = -1)
pheatmap::pheatmap(matTop, scale = "row",
                   clustering_distance_cols = d,
                   clustering_method = "ward.D2",
                   show_rownames = FALSE)

Hierarchical clustering in plain terms: start with every sample as its own cluster, repeatedly merge the two closest clusters, and record the merge order as a tree. The linkage rule defines “closest” between groups rather than points. Average linkage uses the mean pairwise distance, complete linkage uses the worst-case distance and produces compact balanced clusters, and Ward merges the pair that minimizes the increase in within-cluster sum of squares, which tends to give interpretable, similarly-sized groups on expression data. Spearman correlation is the safer similarity here because it is insensitive to the residual scaling differences that survive normalization, and similarity-measure choice is one of the larger levers on clustering output for omics data generally2.

Whether hierarchical clustering beats k-means depends on what you want. Hierarchical gives you the full nesting structure and does not require choosing k in advance, which is the right property when you have 10 or 20 samples and no prior expectation. K-means is faster, scales to hundreds of thousands of points, and gives you a flat partition, which is what you want for gene modules. The four families worth knowing are partitioning methods (k-means, k-medoids), hierarchical methods, density-based methods (DBSCAN, which can leave points unassigned as noise), and model-based methods (Gaussian mixtures, which give soft probabilistic assignments). For expression data you will mostly use the first two, plus graph-based methods for single cell.

Two distinctions that come up when people move over from supervised modeling: regression predicts a known continuous outcome from labeled examples and is evaluated against held-out truth, while clustering has no labels and no ground truth, so it is evaluated on internal coherence and stability. Segmentation is an application, partitioning a population into actionable groups, and it is usually implemented with clustering algorithms.

5. Gene modules across your timepoints

Now transpose the question: which genes move together over time? Z-score each gene across samples so that the clustering sees shape rather than magnitude, then partition.

z <- t(scale(t(matTop)))            # mean 0, sd 1 per gene
set.seed(42)
km <- kmeans(z, centers = 8, nstart = 50, iter.max = 100)
table(km$cluster)

nstart = 50 is not optional. K-means converges to a local optimum determined by initialization, and a single random start on expression data will give you a different answer each run. With 50 restarts and a fixed seed, the result is reproducible and close to the best partition the algorithm can find at that k.

For time-course data specifically, we prefer soft clustering with Mfuzz over hard k-means. Genes rarely belong cleanly to one temporal program, and fuzzy c-means gives each gene a membership value per cluster, so you can threshold at membership > 0.6 to get the core of each module and discard the ambiguous tail.

library(Mfuzz)
eset <- ExpressionSet(assayData = z)
m    <- mestimate(eset)             # typically 1.15 - 1.30 for expression data
cl   <- mfuzz(eset, c = 8, m = m)
mfuzz.plot2(eset, cl = cl, mfrow = c(2, 4), time.labels = samples$week)

WGCNA is the other standard option and produces modules via a soft-thresholded correlation network. It is worth using, but only if you have at least 15 and preferably 20 or more samples. Below that the correlation estimates are too noisy and it will confidently return modules that do not replicate.

6. Choose k and test whether the clusters are real

Never accept a k because it produced a pretty heatmap. Use silhouette width, which for each point compares its mean distance to its own cluster against its mean distance to the nearest other cluster, scaled to [-1, 1].

library(cluster)
sil <- sapply(2:15, function(k) {
  cl <- kmeans(z, centers = k, nstart = 25)$cluster
  mean(silhouette(cl, dist(z))[, 3])
})
plot(2:15, sil, type = "b", xlab = "k", ylab = "mean silhouette")

Mean silhouette above 0.5 indicates reasonable separation, 0.25 to 0.5 is weak structure, and below 0.25 means the partition is essentially arbitrary. Expression data often sits in the weak band, which is honest information about the data rather than a reason to keep trying algorithms until one gives a high score.

Then test stability by resampling. fpc::clusterboot bootstraps the samples, reclusters, and reports the mean Jaccard similarity of each recovered cluster to the original.

library(fpc)
cb <- clusterboot(z, B = 100, clustermethod = kmeansCBI, krange = 8, seed = 42)
cb$bootmean      # per-cluster stability

The conventional reading is that a mean Jaccard below 0.6 means the cluster is dissolved and should not be interpreted, 0.6 to 0.75 indicates a pattern with some support, and above 0.85 is a stable cluster. For sample-level subtypes, ConsensusClusterPlus with 1,000 subsamples at 80% of items is the more common approach and gives you consensus matrices to inspect by eye. This kind of stability testing is exactly what separates a subtype claim that replicates from one that does not, and integrative studies that report molecular subgroups build them on repeated resampling across multiple data layers rather than a single partition3.

7. Single-cell data: build a graph, then partition it

If you have single-cell RNA-seq, the distance-based methods above break down. The matrix is sparse, mostly zeros, and dropout (a transcript present in the cell but not captured) makes pairwise distances unreliable, which is the specific problem that low-rank and self-representation methods were developed to address4. The standard pipeline builds a k-nearest-neighbor graph in PC space and partitions the graph.

import scanpy as sc
adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")
adata.var_names_make_unique()

sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True, log1p=False)
adata = adata[adata.obs.pct_counts_mt < 15].copy()

adata.layers["counts"] = adata.X.copy()
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor="seurat_v3",
                            layer="counts")
adata = adata[:, adata.var.highly_variable].copy()
sc.pp.scale(adata, max_value=10)

sc.tl.pca(adata, n_comps=50, svd_solver="arpack")
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)
for r in (0.2, 0.4, 0.6, 1.0):
    sc.tl.leiden(adata, resolution=r, key_added=f"leiden_{r}", flavor="igraph")

Sweep the resolution rather than trusting the default of 1.0, and use sc.tl.dendrogram plus marker-gene checks to decide where clusters stop being distinguishable by biology. The kNN graph is where the biology enters, and the similarity metric used to build it is the main assumption in the whole pipeline: optimal-transport-based similarity, for instance, recovers cell-cell relationships in single-cell omics better than Euclidean or correlation distance in several benchmarks5. If your goal is cell type assignment rather than novel-population discovery, marker-gene-guided methods are the more direct route and avoid the circularity of clustering first and labeling after6.

8. Annotate the modules

A cluster is not a result until you know what is in it. For gene modules, run enrichment against Hallmark or Reactome, using all expressed genes as the background universe rather than the whole genome.

library(fgsea)
paths <- gmtPathways("h.all.v2023.2.Hs.symbols.gmt")
mod   <- names(km$cluster)[km$cluster == 3]
res   <- fora(pathways = paths, genes = mod, universe = rownames(mat))
head(res[order(res$padj), ], 10)

For single-cell, sc.tl.rank_genes_groups(adata, "leiden_0.4", method="wilcoxon") gives per-cluster markers. Treat the resulting p-values as a ranking device and not as inference, because the clusters were defined using the same data that generates the test statistic. If you also have proteomics and blood biomarkers for the same timepoints, integrating layers before clustering rather than clustering each separately tends to recover structure that any single layer misses, which is the motivation behind sequential regularization and deep embedded multi-omics methods78.

Common problems

Batch dominating PC1 is the most frequent failure. If your samples were sequenced in two runs and PC1 separates them cleanly, every downstream cluster is a batch label. Confirm by coloring the PCA plot by run date, library prep kit, and RNA integrity number. If the batch is confounded with your variable of interest, no correction will save it. If it is not confounded, limma::removeBatchEffect on the VST matrix (for visualization and clustering only, never as input to differential testing) or Harmony on the PCs for single cell will usually resolve it.

RNA degradation is the second. A sample with RIN below 7 shows 3’ bias and depleted long transcripts, and it will sit as an outlier that looks biological. Check the 3’/5’ coverage ratio from Picard’s CollectRnaSeqMetrics before believing any outlier.

Whole blood brings two specific hazards. Globin transcripts (HBB, HBA1, HBA2) can consume 50 to 70 percent of reads without depletion, and neutrophil fraction varies enough day to day that it will drive a principal component on its own. Deconvolution with CIBERSORTx or a simple regression on measured complete blood count values tells you whether a cluster reflects cell composition rather than regulatory change.

Sex-linked genes are a giveaway in pooled data. If XIST and RPS4Y1 are among your top variable genes, you are clustering by donor sex, and you should remove chrX and chrY genes before rerunning.

Finally, a personal longitudinal profile has a sample size problem that no method fixes. Eight timepoints give you eight observations, and correlation-based modules built on eight points will include many spurious pairings. Prefer fewer, larger modules, require bootstrap Jaccard above 0.75 before interpreting a module at all, and treat any striking pattern as a question to ask at the next collection rather than a conclusion. If a cluster maps onto something with clinical weight, such as an inflammatory or coagulation signature, that is a conversation with a physician who can order confirmatory testing, not an endpoint you can reach from the matrix.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. Aleesa E. Chua, Leah D. Pfeifer, Emily R. Sekera, et al. Workflow for Evaluating Normalization Tools for Omics Data Using Supervised and Unsupervised Machine Learning. Journal of the American Society for Mass Spectrometry, 2023. https://doi.org/10.1021/jasms.3c00295 ↩

  2. Tianyi Shi, Xiucai Ye, Zeng Zou, et al. PySimi: a unified framework for similarity measure evaluation in spectral clustering with applications to omics data. Frontiers in Genetics, 2026. https://doi.org/10.3389/fgene.2026.1913487 ↩

  3. Matthias Dottermusch, Alice Ryba, Antonia Gocke, et al. Multi-omics integration unravels four molecular subgroups of corticotroph pituitary neuroendocrine tumours with distinct clinicopathological features. Nature Communications, 2026. https://doi.org/10.1038/s41467-026-76292-y ↩

  4. Ye-Sen Sun, Le Ou-Yang, Dao-Qing Dai. LRSK: a low-rank self-representation K -means method for clustering single-cell RNA-sequencing data. Molecular Omics, 2020. https://doi.org/10.1039/d0mo00034e ↩

  5. Geert-Jan Huizing, Gabriel Peyré, Laura Cantini. Optimal transport improves cell–cell similarity inference in single-cell omics data. Bioinformatics, 2022. https://doi.org/10.1093/bioinformatics/btac084 ↩

  6. Shahriar Rahman Niloy, Toushif Muktashid Hasan, Md. Saiduzzaman Apu, et al. Single‐cell marker gene clustering: A unified deep learning framework for marker gene‐based clustering of single‐cell RNA‐sequencing data. Quantitative Biology, 2026. https://doi.org/10.1002/qub2.70047 ↩

  7. Sunghwan Kim, Steffi Oesterreich, Seyoung Kim, et al. Integrative clustering of multi-level omics data for disease subtype discovery using sequential double regularization. Biostatistics, 2016. https://doi.org/10.1093/biostatistics/kxw039 ↩

  8. Jiawei Li, Taoyuan Ye, Yilang Xiao, et al. UDEC-MO: an uncertainty-guided deep embedded clustering framework for bulk and single-cell multi-omics data. Briefings in Bioinformatics, 2026. https://doi.org/10.1093/bib/bbag435 ↩