When you first dive into genomics with R, the sheer volume of data is usually the first thing that hits you. We aren't talking about a few hundred rows of a CSV; we're talking about tens of thousands of genes across multiple samples. To keep this manageable, I want to show you how to handle a basic differential expression analysis. We'll simulate a small dataset representing gene expression in "Healthy" vs "Tumor" samples, as that's the bread and butter of this field.

Setting up the expression matrix

In a real-world scenario, you'd likely be using BiocManager to pull in packages like DESeq2 or edgeR. But before you let a black-box package do the math, you need to understand the underlying data structure. Usually, you're dealing with a matrix where rows are genes and columns are samples.

# Let's simulate some data: 1000 genes, 6 samples (3 Control, 3 Treatment)
set.seed(42)
genes <- paste0("Gene_", 1:1000)
samples <- c("Ctrl1", "Ctrl2", "Ctrl3", "Tumor1", "Tumor2", "Tumor3")

# Create a matrix of random counts
counts <- matrix(rpois(6000, lambda = 100), nrow = 1000, dimnames = list(genes, samples))

# Let's force a few genes to be "upregulated" in the Tumor samples
counts[1:10, 4:6] <- counts[1:10, 4:6] * 5 

expr_df <- as.data.frame(counts)

The trap of calculating raw fold change

Now, the goal is to find which genes are expressed differently between the two groups. My first instinct—and a mistake I made early in my career—is to just take the average of the Tumor group and divide it by the average of the Control group. Let's see what happens when we do that for our first few genes.

# Calculating raw fold change
ctrl_mean <- rowMeans(expr_df[, 1:3])
tumor_mean <- rowMeans(expr_df[, 4:6])
raw_fc <- tumor_mean / ctrl_mean

head(raw_fc)

The problem here is that raw fold change is asymmetric. A gene that goes from 10 to 20 has a fold change of 2.0, but a gene that goes from 20 to 10 has a fold change of 0.5. If you're plotting this or performing statistical tests, this asymmetry creates a massive bias toward upregulated genes. It makes the data incredibly hard to visualize on a linear scale.

Normalization and the Log2 transformation

To fix this, we use a $\log_2$ transformation. This is standard practice in bioinformatics. By taking the $\log_2$ of the ratio (or the difference of the logs), a 2-fold increase becomes 1, and a 2-fold decrease becomes -1. It centers the data around zero, which is much more intuitive.

# The professional way: log2(Tumor / Control)
# We add a small constant (1) to avoid log(0) if we have zero counts
log2_fc <- log2((tumor_mean + 1) / (ctrl_mean + 1))

# Now the data is symmetric. Let's look at those first 10 genes we manipulated.
head(log2_fc, 10)

Now we have a "Log2 Fold Change" (L2FC). This is the number you'll see in almost every genomic paper. If L2FC is 2, the gene expression increased by $2^2$ (4 times). If it's -2, it decreased by 4 times.

Filtering for biologically significant genes

A fold change alone isn't enough; we need to know if the change is statistically significant. For the sake of this basics lesson, we'll use a simple t-test across the rows to generate p-values. In a production pipeline, you'd use a more sophisticated model that accounts for the overdispersion of count data, but the logic remains the same: combine the magnitude of change with the statistical confidence.

# Apply a t-test to every row
p_values <- apply(expr_df, 1, function(x) {
  t.test(x[1:3], x[4:6])$p.value
})

# Combine results into a final analysis table
results <- data.frame(
  Gene = genes,
  L2FC = log2_fc,
  PValue = p_values
)

# Filter for genes that are upregulated (L2FC > 1) AND significant (P < 0.05)
sig_genes <- subset(results, L2FC > 1 & PValue < 0.05)
head(sig_genes)

I usually suggest filtering by both a "fold-change threshold" and a "p-value threshold." If you only filter by p-value, you might find genes that are statistically different but only changed by 5%, which usually isn't biologically meaningful. By combining them, you narrow your list down to the genes that actually matter.

📋 Practical Task

Exercise: Identifying Downregulated Genes in a Simulated Dataset

You have been provided with a gene expression matrix exp_matrix containing 500 genes across 4 samples: Control_A, Control_B, Treat_A, Treat_B. Your goal is to identify the "downregulated" genes.

Write a script that performs the following:

• Calculates the $\log_2$ fold change between the Treatment group and the Control group.
• Calculates a p-value for each gene using a t-test.
• Filters the results to find genes that are downregulated (defined as a $\log_2$ fold change of less than -1) and have a p-value less than 0.05.
• Returns the final list of gene names that meet both criteria.

# Starter data
set.seed(123)
exp_matrix <- matrix(rpois(2000, lambda = 50), nrow = 500, 
                     dimnames = list(paste0("Gene", 1:500), 
                                     c("Control_A", "Control_B", "Treat_A", "Treat_B")))

# Force some genes to be downregulated in Treatment
exp_matrix[100:110, 3:4] <- exp_matrix[100:110, 3:4] / 4

# Your code here: