Your First RNA-seq Analysis

The question we're actually answering

Before any code, let's be clear on what RNA-seq is for. Every cell carries the same DNA, but at any moment it's only reading some of those genes. The genes a cell is actively "reading" get transcribed into RNA. So if you measure how much RNA there is for each gene, you get a snapshot of what the cell is actually doing.

RNA-seq is the technology that counts that RNA, gene by gene, across thousands of genes at once. The classic experiment compares two groups - say, untreated cells vs. cells given a drug - and asks a deceptively simple question:

"Which genes change their activity between the two conditions?"

That's differential expression analysis, and it's the workhorse of modern molecular biology - used to understand diseases, find drug targets, and figure out how cells respond to their environment. That's exactly what you're about to do.

The plan

Five steps. Each one is a few lines of code, and we'll explain what every line is doing and - more importantly - why.

1. Set up - install one package.
2. Load the data - a count matrix and a table describing the samples.
3. Run the analysis - let DESeq2 do the heavy statistical lifting.
4. Read the results - find the genes that changed.
5. Visualize - draw a volcano plot, the signature figure of RNA-seq.

1Set up

We'll use DESeq2, the standard R / Bioconductor package for differential expression - the original method the whole field is built on. Make sure your Colab notebook is set to the R runtime (Runtime → Change runtime type → R), then install it once:

# Installing Bioconductor packages takes a few minutes the first time - grab a coffee.
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("DESeq2", update = FALSE, ask = FALSE)

Then load the library:

library(DESeq2)

2Load the data

We'll use a compact, ready-made example count matrix so everything runs in seconds and you can focus on the workflow. (When you're ready for a real published dataset, see the "Level up" box near the end - the steps are identical.)

Two files: a count matrix (genes × samples) and a metadata table that says which condition each sample belongs to.

base <- "https://raw.githubusercontent.com/owkin/PyDESeq2/main/datasets/synthetic/"

# Count matrix: rows = genes, columns = samples - exactly what DESeq2 wants.
counts <- as.matrix(read.csv(paste0(base, "test_counts.csv"), row.names = 1))

# Metadata: one row per sample, with a "condition" column (values A and B).
metadata <- read.csv(paste0(base, "test_metadata.csv"), row.names = 1)
metadata$condition <- factor(metadata$condition)

# Keep the metadata rows in the same order as the count columns.
metadata <- metadata[colnames(counts), , drop = FALSE]

dim(counts)        # genes x samples
head(counts)

Take a moment to actually look at head(counts) and head(metadata). Getting comfortable looking at your data before analyzing it is one of the most underrated habits in bioinformatics - most real-world mistakes are caught just by eyeballing the table.

3Run the analysis

This is where DESeq2 earns its reputation. In two lines it normalizes the data, models the counts statistically, and prepares everything needed to test each gene.

# Build the DESeq2 dataset. design = ~condition is the question we're asking.
dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData   = metadata,
                              design    = ~ condition)

dds <- DESeq(dds)   # fit the model - this does the real work

The design = ~ condition part is you telling DESeq2: "the thing I care about is the difference between conditions A and B." That single choice frames the entire analysis.

4Read the results

Now we ask, for every gene: is the difference between A and B big enough - and consistent enough - to be real?

# contrast = compare condition B vs A
res <- results(dds, contrast = c("condition", "B", "A"))
res <- res[order(res$padj), ]   # sort by adjusted p-value
head(res, 10)                    # the 10 most convincing hits

The results table has a few columns that matter enormously:

• log2FoldChange - how much the gene changed, on a log2 scale. +1 means it doubled in B vs A; −1 means it halved. Sign tells you direction, size tells you magnitude.
• pvalue - the raw probability that you'd see a change this big by pure chance.
• padj - the p-value adjusted for the fact that we tested thousands of genes at once. This is the one you trust, not the raw p-value.

5Visualize: the volcano plot

The volcano plot is the iconic RNA-seq figure. Each dot is a gene. The x-axis is how much it changed (fold change); the y-axis is how confident we are (significance). The genes that matter - big change and high confidence - fly up into the top corners, which is where the name comes from.

# Turn the results into a plain data frame and drop genes with no p-value
df <- as.data.frame(res)
df <- df[!is.na(df$padj), ]
df$sig <- df$padj < 0.05

plot(df$log2FoldChange, -log10(df$padj),
     pch = 20,
     col = ifelse(df$sig, "#a855f7", "#cbd5e1"),
     xlab = "log2 fold change (B vs A)",
     ylab = "-log10 adjusted p-value",
     main = "Volcano plot - your first one!")
abline(h = -log10(0.05), lty = 2, col = "#64748b")   # the padj < 0.05 line

That's it. You just took a count matrix, ran a real statistical workflow, identified the genes that change between conditions, and made the figure you'll see in thousands of published papers. The exact same five steps scale up to a real experiment with 20,000 genes.

# install once
!pip install pydeseq2

import pandas as pd
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

base = "https://raw.githubusercontent.com/owkin/PyDESeq2/main/datasets/synthetic/"
counts = pd.read_csv(base + "test_counts.csv", index_col=0).T   # samples x genes for PyDESeq2
metadata = pd.read_csv(base + "test_metadata.csv", index_col=0)

dds = DeseqDataSet(counts=counts, metadata=metadata, design="~condition")
dds.deseq2()

stat_res = DeseqStats(dds, contrast=["condition", "B", "A"])
stat_res.summary()
res = stat_res.results_df.sort_values("padj")
res.head(10)

Now do it on real data: the airway dataset

The example matrix above kept things fast so you could see the workflow. Here is the same analysis on a real published experiment: the airway study (Himes et al., 2014, GEO accession GSE52778), human airway smooth-muscle cells treated with dexamethasone, a steroid used in asthma. In R, the Bioconductor airway package ships this exact dataset ready to load, so it is fully reproducible. Run this in Colab with an R runtime, or in RStudio.

# install once (a few minutes the first time)
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("airway", "DESeq2"))

library(airway); library(DESeq2)
data("airway")                              # a real RNA-seq dataset, ready to use
airway$dex <- relevel(airway$dex, "untrt")   # make "untreated" the baseline

# note the design: cell line is a batch term, treatment is what we test
dds <- DESeqDataSet(airway, design = ~ cell + dex)
dds <- DESeq(dds)
res <- results(dds)

summary(res)                              # how many genes up/down at padj < 0.1
head(res[order(res$padj), ])              # the strongest hits

Check your understanding