RNA-seq Analysis Made Simple: From Raw Counts to Biological Insight — My Review

On transforming noise into meaning and the reality of actually getting there

“What if I told you that most RNA-seq mistakes happen before you even run differential expression?

Raw counts look simple, but analyzing them correctly is not.”

The real problem is this: We start with a matrix of raw counts, genes by samples, and assume it reflects biology. But it doesn’t, at least not yet.

The central challenge is transforming messy, biased, high-noise count data into statistically valid and biologically meaningful differential expression results without being misled by sequencing depth, technical variation, low-expression genes, or incorrect statistical assumptions. Because without that transformation step, any downstream conclusion is fragile. You might be detecting artifacts of the experiment rather than true biological signals. In other words, the “biology” you think you’re seeing could simply be a reflection of how the data was generated.

This article works through that problem in full. It replicates a complete RNA-seq workflow using edgeR, limma, and Glimma applied to a mouse mammary gland dataset from Sheridan et al. (2015), covering every step I took, not just as a sequence of commands, but as a conceptual progression from unstable raw counts to statistically defensible biological insight.

A Personal Note on Process!

My replication of this workflow began on my school laptop, which turned out to be a mistake as my institution’s network had blocked R’s ability to communicate with package repositories. Before I could even process any information or set out a plan, I had to find a way to even download the data packs!

Fortunately, I was able to switch to my sister’s laptop, which did not have R or RStudio blocked. The environment had to be built using edgeR, limma, Glimma, the Mus.musculus annotation database, and the RNAseq123 package, and assembling the right set of packages for this specific workflow required a kind of archaeological patience. Version conflicts. Deprecated functions. Packages that seemed correct but weren’t quite what the pipeline expected. The dataset itself, downloaded from GEO, arrived as a tar archive with more files than I needed and no particular guidance on which ones to use.

Eventually, through iteration that felt at the time like failure but in retrospect was simply the process, everything fell into place. And the insight I took from it was this: the analytical pipeline, once the environment is correctly configured, is logical, well-documented, and genuinely elegant.

The Problem Beneath the Data

When raw RNA-seq data arrives, it presents itself as something simple: a table of counts, genes along one axis, samples along the other. The temptation is to treat it as ground truth. It is not.

Several structural problems distort the count matrix before any analysis begins:

Sequencing depth is not uniform across samples. A sample sequenced to greater depth will produce higher counts for every gene, not because the biology differs, but because more reads were generated. A gene recording 500 counts in one sample and 250 in another might be expressing at identical levels; the difference is entirely technical. Without correcting for this, every comparison is potentially misleading.

Lowly expressed genes are not informative; they are noise. Thousands of genes in any given experiment register near-zero counts across all conditions. They cannot be differentially expressed, because there is insufficient signal to detect anything. Including them dilutes statistical power and inflates the number of tests performed.

Variance scales with expression level. High-count genes behave differently from low-count genes in terms of their statistical variability. This heteroscedasticity violates the assumptions of linear models, which are built on the premise of stable variance. Apply a linear model to untransformed count data, and you are solving the wrong problem.

Technical variation can masquerade as biology. Sequencing lane, processing batch, and sample handling are factors that can systematically shift expression profiles in ways indistinguishable from genuine biological signals, unless they are explicitly accounted for.

The consequences of ignoring these issues are not subtle. Genes that appear significant may reflect nothing more than how the experiment was run. Patterns that seem reproducible may be artefacts of the dataset rather than properties of the underlying system. The analysis produces results, but they are not results you can trust.

The edgeR–limma pipeline exists to address each of these distortions in turn.

The Dataset

The experiment profiled three cell populations sorted from the mammary glands of female virgin mice: Basal, Luminal Progenitor (LP), and Mature Luminal (ML), each in biological triplicate, yielding nine samples in total. Reads were sequenced on an Illumina HiSeq 2000 across three lanes, aligned to the mouse reference genome mm10, and summarised to gene-level counts using featureCounts. The data is publicly available on the Gene Expression Omnibus under accession number GSE63310.

What makes this dataset well-suited to learning is that the biology is not mysterious. Basal cells are known to be transcriptionally distinct from luminal populations. That prior knowledge becomes a continuous sanity check throughout the analysis: at each stage, we can ask whether what the data shows is consistent with what the biology should be doing.

Packaging the Data

The raw data arrives as nine separate text files, one per sample. The readDGE function from edgeR consolidates them into a single DGEList object, a structured container in which the count matrix, sample metadata, and gene annotations are kept together and in relation to one another.

x <- readDGE(files, columns=c(1,3))

After import, the dataset contains 27,179 genes across 9 samples. These numbers are not yet meaningful for comparison; they are shaped by sequencing depth and technical variation, not biology alone. alone. The import step is simply the act of gathering the raw material.

What immediately follows is arguably more important: attaching biological and technical identity to each sample.

group <- as.factor(c(“LP”,”ML”,”Basal”,”Basal”,”Basal”,”ML”,”ML”,”LP”,”LP”))

lane <- as.factor(rep(c(“L004”,”L006",”L008"), c(3,4,2)))

These two vectors, cell type and sequencing lane, are not administrative bookkeeping. They are the variables through which the statistical model will later distinguish signal from noise.

Finally, the gene identifiers in the count matrix are Entrez IDs: numerical codes that are precise but opaque. Linking them to gene symbols and chromosomal locations transforms the data from a table of numbers into something a biologist can actually read.

genes <- select(Mus.musculus, keys=geneid, columns=c(“SYMBOL”,”TXCHROM”), keytype=”ENTREZID”)

genes <- genes[!duplicated(genes$ENTREZID),]

The deduplication step on that second line is easy to overlook and consequential to miss. Some Entrez IDs map to multiple gene symbols in the annotation database. Without removing the duplicates, downstream dimension mismatches emerge silently, corrupting merges without raising any obvious errors.

• Filtering, Normalisation, and Transformation

These three operations happen in sequence, and that sequence is not arbitrary.

Filtering removes genes that cannot contribute to differential expression analysis because they lack sufficient data. The filterByExpr function applies an adaptive threshold approximately 0.2 CPM in this dataset, retaining only genes expressed at a meaningful level in at least some samples.

keep.exprs <- filterByExpr(x, group=group)

x <- x[keep.exprs, keep.lib.sizes=FALSE]

​

The dataset contracts from 27,179 genes to 16,624, nearly 40% removed. This is not loss; it is precision. The density plots make the effect visible: a heavy left tail of near-zero values, present before filtering, simply disappears.

Normalisation addresses compositional bias, the distortion that arises when a few highly expressed genes consume a disproportionate share of sequencing reads, effectively compressing the apparent expression of everything else. TMM normalisation corrects for this by computing scaling factors that align the distributions across samples.

x <- calcNormFactors(x, method=”TMM”)

​

In the real dataset, the normalisation factors were close to 1, a signal that compositional bias was modest. The demonstration with artificially extreme values (samples scaled to 5% and 500% of their true counts) makes the necessity of this step more visible: without it, the distortion is dramatic; with it, the distributions align cleanly.

Log-CPM transformation is then computed for exploratory analysis and visualisation. The log scale compresses the wide dynamic range of expression data into something tractable and stabilises variance enough to make exploratory plots interpretable.

lcpm <- cpm(x, log=TRUE)

​

plotMDS(lcpm)

​

MDS plots reveal whether samples cluster by biological condition or technical factors. This step acts as a diagnostic checkpoint. If expected patterns do not emerge, it signals deeper issues that must be addressed before proceeding.

Only after careful preparation can statistical testing begin!

The three cell types cluster distinctly, with basal cells separated from both luminal populations, LP and ML, which are distinct from each other. Basal’s distance from the others is already suggesting, before a single test has been run, that it is the most transcriptionally divergent population.

Dimensions 3 and 4 reveal the sequencing lane effect, a technical gradient running through the data. Its presence here is not alarming; it is information. It confirms that the lane must be included in the model so that its contribution can be partitioned away from the biological signal.

Building the Statistical Model

The design matrix is where biological intuition becomes mathematical form. It encodes both what we are interested in, cell type, and what we need to control for sequencing lane.

design <- model.matrix(~0+group+lane)

​

Contrasts then define the specific questions we want to answer: which cell types, compared to which others?

design <- model.matrix(~0+group+lane)

co ntr. matrix <- makeContrasts(

BasalvsLP = Basal-LP,

BasalvsML = Basal — ML,

LPvsML = LP — ML,

levels = colnames(design))​

The design matrix encodes the experimental structure, translating biological questions into mathematical form. Contrasts define specific comparisons, turning abstract hypotheses into testable statements.

Removing heteroscedasticity from count data

Variance in RNA-seq data depends on mean expression. The voom method addresses this.

v <- voom(x, design)

(​Figure 4: Means (x-axis) and variances (y-axis) of each gene are plotted to show the dependence between the two before voom is applied to the data (left panel) and how the trend is removed after voom precision weights are applied to the data (right panel) — RNA-seq analysis is easy as 1–2–3 with limma, Glimma, and edgeR )

The plot on the left is created within the voom function, which extracts residual variances from fitting linear models to log-CPM transformed data. Variances are then rescaled to quarter-root variances (or square-root of standard deviations) and plotted against the average log2 count for each gene. The plot on the right is created using plotSA, which plots log2 residual standard deviations against mean log-CPM values. In both plots, each black dot represents a gene. On the left plot, the red curve shows the estimated mean-variance trend used to compute the voom weights. On the right plot, the average log2 residual standard deviation estimated by the empirical Bayes algorithm is marked by a horizontal blue line.

By estimating precision weights, voom transforms the data into a form suitable for linear modeling. This step bridges the gap between count-based data and methods originally developed for continuous data.

What the Data Shows

summary(decideTests(efit))

​

Contrast

Down

Not Significant

Up

Basal vs LP4,6467,1184,860Basal vs ML4,9367,0084,680LP vs ML3,14110,9532,530

Approximately 9,500 differentially expressed genes in both basal comparisons; approximately 5,600 between LP and ML. The pattern is consistent: comparisons involving basal yield the most DE genes, and basal’s separation in the MDS plot now has a quantitative counterpart. The data has been internally coherent from the first exploratory step to the final result.

The Venn diagram of overlapping DE genes across contrasts reinforces this: the majority of genes that differ between basal and LP also differ between basal and ML. Basal does not merely differ from one luminal population; it operates under an entirely different transcriptional programme.

Seeing the Results

MD plots mean-difference plots offer the genome-wide view: log fold change on the vertical axis, average expression on the horizontal, significantly changing genes coloured by direction. Every gene in the dataset is visible at once.

The static plot is useful. Glimma makes it alive:

glMDPlot(tfit, coef=1, status=dt, side.main=”ENTREZID”, counts=lcpm, groups=group, launch=FALSE)

​

Glimma generates an interactive HTML page in which every point in the MD plot is selectable. Choose any gene, and the right panel renders its expression across all nine samples individually. A search bar allows lookup by symbol. What had been a summary becomes an investigation of the difference between a map and the territory it describes.

The heatmap of the top 100 DE genes from the Basal vs LP comparison makes the biology legible at a glance: ML and LP cluster together in their expression patterns, sharing the transcriptional character of luminal cells, while basal stands entirely apart.

Gene set testing with a camera.

load(system.file(“extdata”, “mouse_c2_v5p1.rda”, package = “RNAseq123”))

idx <- ids2indices(Mm.c2,id=rownames(v))

cam.BasalvsLP <- camera(v,idx,design,contrast=contr.matrix[,1])

head(cam.BasalvsLP,5)

• cam.LPvsML <- camera(v,idx,design,contrast=contr.matrix[,3]) head(cam.LPvsML,5)

(RNA-seq analysis is easy as 1–2–3 with limma, Glimma, and edgeR)

Figure 7: Barcode plot of LIM_MAMMARY_LUMINAL_MATURE_UP (red bars, top of plot) and LIM_MAMMARY_LUMINAL_MATURE_DN (blue bars, bottom of plot) gene sets in the LP versus ML contrast

• For each set, an enrichment line that shows the relative enrichment of the vertical bars in each part of the plot is displayed. The experiment of Lim et al. (Lim et al. 2010) is very similar to the current one, with the same sorting strategy used to obtain the different cell populations, except that microarrays were used instead of RNA-seq to profile gene expression. Note that the inverse correlation (the up gene set is down and the down gene set is up) is a result of the way the contrast has been set up (LP versus ML); if reversed, the directionality would agree.

Finally, the gene set testing shifts the focus from individual genes to coordinated biological processes. It asks not just whether genes change, but whether pathways are collectively perturbed.

Three Problems Worth Documenting

I had particular trouble with the duplicate gene annotations. Running select() against the Mus. musculus database returned more rows than the count matrix had genes, a sign that certain Entrez IDs mapped to multiple symbols. The !duplicated(genes$ENTREZID) filter resolves this cleanly, but I realized that its absence produces errors in downstream merges that are difficult to trace back to their origin.

Secondly, the mean-variance trend produced by an early run of voom was erratic rather than smooth; the characteristic curve was disrupted by noise at low expression values. The cause was ordering: I had run voom before filtering. Lowly expressed near-zero-count genes distort the mean-variance model precisely in the region where it matters most. Filter first. Then normalise. Then voom!

On an initial attempt, the primary separation in the first two MDS dimensions ran along the sequencing lane axis rather than the cell type axis, a sign that technical variation was outweighing biological signal. Two errors contributed: the group factor had been assigned in the wrong sample order, and the lane had not been included in the design matrix. The correction required verifying metadata with the table(group, lane) and updating the model.

A Closing Reflection

We began with a deceptively simple matrix of counts and a fundamental question: how do we extract truth from noise? Over the course of this replicate, I learned that the answer is not a single method, but a disciplined sequence of transformations, each addressing a specific distortion in the data.

This workflow showed that reliable RNA-seq analysis depends on a sequence of carefully ordered statistical corrections, each designed to separate biological signal from technical noise.

It also, on occasion, begins with switching to your sister’s laptop. That is part of the process, too.

Footnotes

RNA-seq analysis is easy as 1–2–3 with limma, Glimma, and edgeR — https://bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow.html#normalising-gene-expression-distributions

Charity Law1, Monther Alhamdoosh2, Shian Su3, Xueyi Dong3, Luyi Tian1, Gordon K. Smyth4 and Matthew E. Ritchie5

1The Walter and Eliza Hall Institute of Medical Research, 1G Royal Parade, Parkville, VIC 3052, Melbourne, Australia; Department of Medical Biology, The University of Melbourne, Parkville, VIC 3010, Melbourne, Australia 2CSL Limited, Bio21 Institute, 30 Flemington Road, Parkville, Victoria 3010, Australia 3The Walter and Eliza Hall Institute of Medical Research, 1G Royal Parade, Parkville, VIC 3052, Melbourne, Australia 4The Walter and Eliza Hall Institute of Medical Research, 1G Royal Parade, Parkville, VIC 3052, Melbourne, Australia; School of Mathematics and Statistics, The University of Melbourne, Parkville, VIC 3010, Melbourne, Australia 5The Walter and Eliza Hall Institute of Medical Research, 1G Royal Parade, Parkville, VIC 3052, Melbourne, Australia; Department of Medical Biology, The University of Melbourne, Parkville, VIC 3010, Melbourne, Australia; School of Mathematics and Statistics, The University of Melbourne, Parkville, VIC 3010, Melbourne, Australia

17 December 2018