Microarray Data Analysis

TL;DR: Analyze Affymetrix microarray data in R using Bioconductor. Load CEL files with affy/oligo, normalize with RMA, perform differential expression with limma. Still relevant for existing datasets and GEO reanalysis despite RNA-seq dominance.

---

Why Microarrays in 2026?

While RNA-seq has largely replaced microarrays for new studies, you'll encounter microarray data when:

• ✅ Reanalyzing published datasets (GEO has millions of samples)
• ✅ Meta-analysis across studies (standardized platforms)
• ✅ Clinical diagnostics (FDA-approved platforms)
• ✅ Budget constraints (still cheaper for large cohorts)

Microarray vs RNA-seq

Feature	Microarray	RNA-seq
Cost	$50-100/sample	$200-500/sample
Novel transcripts	❌ No	✅ Yes
Dynamic range	Limited (~3 orders)	Wide (>5 orders)
Quantification	Relative	Absolute
Normalization	Cross-sample required	Per-sample possible
Data size	~10 MB/sample	~2-10 GB/sample
Best for	Known genes, large cohorts	Discovery, isoforms, novel genes

---

Workflow Overview

CEL files → RMA normalization → Differential expression → Functional analysis

---

Part 1: Data Structures

ExpressionSet Object

Bioconductor's standard container for expression data:

library(Biobase)

# Structure
ExpressionSet (storageMode: lockedEnvironment)
  assayData: 12625 features, 22 samples      # Expression matrix
  phenoData:                                 # Sample metadata
    sampleNames: CLL11.CEL, CLL12.CEL, ...
    varLabels: SampleID, Disease, Age
  featureData:                               # Probe annotations
  annotation: hgu95av2                       # Chip type

Extracting Components

# Expression matrix
expr_matrix <- exprs(eset)
# Rows = probes/genes, Columns = samples

# Sample metadata
metadata <- pData(eset)
# Columns = clinical variables (disease, age, etc.)

# Probe annotations
features <- fData(eset)
# Map probes to gene symbols

---

Part 2: Loading Data

From GEO (Recommended)

library(GEOquery)

# Download processed data
gse <- getGEO("GSE1009", GSEMatrix = TRUE)
eset <- gse[[1]]  # ExpressionSet object

# Check platform
annotation(eset)

# View sample info
pData(eset)[1:5, 1:5]

From Raw CEL Files (Affymetrix)

Option 1: affy Package (Older Chips)

library(affy)

# Read CEL files from directory
raw_data <- ReadAffy(celfile.path = "CEL_files/")

# Check data
sampleNames(raw_data)
dim(exprs(raw_data))  # Probe-level

# RMA normalization (background correction + quantile norm + summarization)
eset <- rma(raw_data)

Option 2: oligo Package (Newer Chips)

library(oligo)

# List CEL files
cel_files <- list.celfiles("CEL_files/", full.names = TRUE)

# Read
raw_data <- read.celfiles(cel_files)

# RMA normalization
eset <- rma(raw_data)

Creating ExpressionSet from Matrices

library(Biobase)

# If you have separate expression matrix and metadata
expr_matrix <- read.csv("expression.csv", row.names = 1)
metadata <- read.csv("metadata.csv", row.names = 1)

# Create annotated metadata
pheno <- new("AnnotatedDataFrame", data = metadata)

# Create ExpressionSet
eset <- ExpressionSet(
  assayData = as.matrix(expr_matrix),
  phenoData = pheno,
  annotation = "hgu95av2"  # Chip type
)

---

Part 3: Normalization

RMA (Robust Multi-array Average)

Industry standard for Affymetrix data. Combines:

1. Background correction (using PM/MM probes)
2. Quantile normalization (across arrays)
3. Summarization (probe set → gene)

# RMA already includes quantile normalization
eset <- rma(raw_data)

# Verify normalization
boxplot(exprs(eset), main = "Post-RMA", las = 2)

Important: Do NOT run normalizeBetweenArrays() after RMA (already done).

Alternative Methods

# MAS5 (older method)
eset_mas5 <- mas5(raw_data)

# GCRMA (GC-content aware)
library(gcrma)
eset_gcrma <- gcrma(raw_data)

# VSN (variance stabilization)
library(vsn)
eset_vsn <- vsnrma(raw_data)

Recommendation: Stick with RMA unless you have a specific reason.

---

Part 4: Quality Control

Before Normalization

library(affyPLM)

# Fit probe-level model
plm_fit <- fitPLM(raw_data)

# RNA degradation plot
RNAdeg <- AffyRNAdeg(raw_data)
plotAffyRNAdeg(RNAdeg)

# Pseudo-images (spatial artifacts)
image(raw_data[, 1])  # First array

After Normalization

# Boxplots (should be aligned)
boxplot(exprs(eset), main = "Expression Distributions")

# PCA
library(ggplot2)
pca <- prcomp(t(exprs(eset)), scale. = TRUE)
pca_df <- data.frame(PC1 = pca$x[,1], PC2 = pca$x[,2], 
                     Group = pData(eset)$Disease)
ggplot(pca_df, aes(PC1, PC2, color = Group)) + 
  geom_point(size = 3)

# Sample distances heatmap
library(pheatmap)
dists <- dist(t(exprs(eset)))
pheatmap(as.matrix(dists))

---

Part 5: Probe-to-Gene Mapping

Annotate Probes

library(annotate)
library(hgu95av2.db)  # Change to your chip

# Get gene symbols
symbols <- getSYMBOL(featureNames(eset), "hgu95av2.db")
fData(eset)$GeneSymbol <- symbols

# Get Entrez IDs
entrez <- getEG(featureNames(eset), "hgu95av2.db")
fData(eset)$EntrezID <- entrez

# View mappings
head(fData(eset))

Handle Multiple Probes per Gene

library(limma)

# Average expression when multiple probes map to same gene
averaged <- avereps(exprs(eset), ID = fData(eset)$GeneSymbol)

---

Part 6: Differential Expression with limma

Basic Two-Group Comparison

library(limma)

# Design matrix
design <- model.matrix(~ 0 + factor(eset$Disease))
colnames(design) <- c("Control", "Disease")

# Fit linear model
fit <- lmFit(eset, design)

# Define contrast
contrast_matrix <- makeContrasts(
  Disease_vs_Control = Disease - Control,
  levels = design
)

# Fit contrasts
fit2 <- contrasts.fit(fit, contrast_matrix)
fit2 <- eBayes(fit2)

# Extract results
results <- topTable(fit2, coef = "Disease_vs_Control", number = Inf, adjust.method = "BH")
head(results)

Multi-Group Analysis

# Design with multiple groups
design <- model.matrix(~ 0 + eset$CellType)
colnames(design) <- levels(factor(eset$CellType))

fit <- lmFit(eset, design)

# Multiple contrasts
contrast_matrix <- makeContrasts(
  TypeA_vs_TypeB = TypeA - TypeB,
  TypeA_vs_TypeC = TypeA - TypeC,
  TypeB_vs_TypeC = TypeB - TypeC,
  levels = design
)

fit2 <- contrasts.fit(fit, contrast_matrix)
fit2 <- eBayes(fit2)

# View all comparisons
summary(decideTests(fit2))

With Covariates

# Adjust for batch effects or other covariates
design <- model.matrix(~ Disease + Batch + Age, data = pData(eset))

fit <- lmFit(eset, design)
fit <- eBayes(fit)

# Results for Disease effect (adjusted for Batch and Age)
results <- topTable(fit, coef = "DiseaseCase", number = Inf)

---

Part 7: Visualization

Volcano Plot

library(ggplot2)

# Prepare data
results$significant <- results$adj.P.Val < 0.05 & abs(results$logFC) > 1

ggplot(results, aes(x = logFC, y = -log10(P.Value), color = significant)) +
  geom_point(alpha = 0.5) +
  scale_color_manual(values = c("gray", "red")) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
  labs(title = "Volcano Plot", x = "Log2 Fold Change", y = "-Log10 P-value")

Heatmap of Top Genes

# Top 50 DEGs
top_genes <- rownames(results)[1:50]
heatmap_data <- exprs(eset)[top_genes, ]

pheatmap(heatmap_data,
         scale = "row",
         annotation_col = pData(eset)[, c("Disease"), drop = FALSE],
         show_rownames = TRUE,
         fontsize_row = 6)

MA Plot

limma::plotMA(fit2, coef = 1)

---

Part 8: Functional Enrichment

GO Enrichment

library(clusterProfiler)
library(org.Hs.eg.db)

# Get significant genes
sig_genes <- results$EntrezID[results$adj.P.Val < 0.05]
sig_genes <- na.omit(sig_genes)

# GO enrichment
go_results <- enrichGO(
  gene = sig_genes,
  OrgDb = org.Hs.eg.db,
  ont = "BP",  # Biological Process
  pAdjustMethod = "BH",
  qvalueCutoff = 0.05
)

# Visualize
dotplot(go_results, showCategory = 20)

KEGG Pathway

kegg_results <- enrichKEGG(
  gene = sig_genes,
  organism = "hsa",  # Human
  pAdjustMethod = "BH",
  qvalueCutoff = 0.05
)

dotplot(kegg_results)

---

Complete Analysis Script

# Load libraries
library(GEOquery)
library(limma)
library(pheatmap)

# 1. Download data
gse <- getGEO("GSE1009", GSEMatrix = TRUE)
eset <- gse[[1]]

# 2. Filter low-expressed probes
eset <- eset[rowMeans(exprs(eset)) > 5, ]

# 3. Check normalization
boxplot(exprs(eset), las = 2)

# 4. Design matrix
design <- model.matrix(~ 0 + Disease, data = pData(eset))
colnames(design) <- c("Control", "Disease")

# 5. Differential expression
fit <- lmFit(eset, design)
contrast <- makeContrasts(Disease - Control, levels = design)
fit2 <- contrasts.fit(fit, contrast)
fit2 <- eBayes(fit2)

# 6. Results
results <- topTable(fit2, number = Inf, adjust.method = "BH")

# 7. Export
write.csv(results, "DEG_results.csv")

# 8. Heatmap of top 50
top_genes <- rownames(results)[1:50]
pheatmap(exprs(eset)[top_genes, ], scale = "row", 
         annotation_col = pData(eset)[, "Disease", drop = FALSE])

---

Troubleshooting

"Annotation package not found"

Fix: Install chip-specific annotation:

# List available
BiocManager::available("hgu")

# Install your chip
BiocManager::install("hgu95av2.db")

RMA fails with memory error

Fix: Process in batches:

# Reduce memory footprint
rm(raw_data)
gc()

Batch effects visible in PCA

Fix: Use ComBat:

library(sva)
combat_expr <- ComBat(exprs(eset), batch = pData(eset)$Batch)
exprs(eset) <- combat_expr

---

Tools & Software

R Packages

• Core: affy, oligo, limma, Biobase
• Annotation: annotate, org.Hs.eg.db, hgu95av2.db
• Visualization: pheatmap, ggplot2
• Enrichment: clusterProfiler, DOSE

GUI Alternative

Chipster: chipster.csc.fi
Point-and-click interface for microarray analysis (no coding).

---

Resources

• Bioconductor workflows: bioconductor.org/packages/workflows
• limma user guide: bioconductor.org/packages/limma
• GEO database: ncbi.nlm.nih.gov/geo

---

See also:

• GSEA, GSVA and ssGSEA
• Download public sequence data from SRA
• How to: use data.table for more efficient data management

---

Originally posted on WeChat Official Account: bioinfo_studyroom