microarray-analysis

name	microarray-analysis
description	Analyze DNA microarray data from GEO datasets. Use when asked to analyze microarray data, perform differential expression analysis, or when given a GEO accession (GSE*). Handles Affymetrix (CEL files) and Illumina (BeadChip) platforms with QC, normalization, statistical modeling, and gene ID mapping.

DNA Microarray Analysis Skill

This skill guides analysis of DNA microarray gene expression data following established patterns from this project.

Quick Start

When asked to analyze microarray data:

1. Check if a CLAUDE.md file exists in the study directory - it contains study-specific metadata
2. Determine the platform (Affymetrix vs Illumina) to select appropriate methods
3. Create numbered scripts in scripts/ directory (01_, 02_, etc.)
4. Update analysis.md after each analysis step
5. Generate all plots in results/plots/ with plots/plots.md documentation

Directory Structure

Create this structure for each analysis:

analyses/[FirstAuthor]_[GSE#]/
├── CLAUDE.md           # Study metadata (samples, platform, design)
├── analysis.md         # Progressive documentation (update after each step)
├── scripts/
│   ├── 01_load_and_qc.R
│   ├── 02_normalize.R (or 02_normalize_filtered.R)
│   ├── 03_differential_expression.R
│   ├── 04_annotate_genes.R
│   └── 05_go_enrichment.R
├── data/               # Raw and processed data (gitignored)
└── results/
    ├── tables/         # CSV output files
    └── plots/          # PNG/PDF figures
        └── plots.md    # Plot documentation

Script Portability

Use the here package for all file paths to ensure scripts work regardless of working directory:

library(here)

# Set paths relative to project root
base_dir <- here("analyses", "Author_GSE#####")
data_dir <- file.path(base_dir, "data")
results_dir <- file.path(base_dir, "results")
plots_dir <- file.path(results_dir, "plots")
tables_dir <- file.path(results_dir, "tables")

# Source utility functions from project root
source(here("scripts", "utils", "collapse_to_gene.R"))
source(here("scripts", "utils", "enrichment_dotplot.R"))

Never use hardcoded absolute paths like /Users/name/path/to/file. Always use here() for portability.

Probe Tracking

IMPORTANT: Track probe counts at each processing step and report a summary table at the end of QC and normalization scripts. This provides transparency about data filtering.

Initialize Tracking

# Initialize probe tracking at the start of script
probe_tracking <- data.frame(
  step = character(),
  n_probes = integer(),
  n_removed = integer(),
  description = character(),
  stringsAsFactors = FALSE
)

Track at Each Step

Add tracking after each filtering/processing step:

# After loading raw data
probe_tracking <- rbind(probe_tracking, data.frame(
  step = "1. Raw data loaded",
  n_probes = nrow(raw_data),
  n_removed = 0,
  description = "Total probes in raw data file"
))

# After creating expression matrix
probe_tracking <- rbind(probe_tracking, data.frame(
  step = "2. Expression matrix created",
  n_probes = nrow(expr_mat),
  n_removed = nrow(raw_data) - nrow(expr_mat),
  description = "Probes with valid signal columns"
))

# After detection filtering (Illumina)
probe_tracking <- rbind(probe_tracking, data.frame(
  step = "3. Detection filter",
  n_probes = nrowexpr_filtered
  n_removed  nrowexpr_mat  nrowexpr_filtered
  description  sprintf min_samples_detected



probe_tracking  rbindprobe_tracking data.frame
  step  
  n_probes  nrowexpr_filtered
  n_removed  n_before  nrowexpr_filtered
  description  

Report Summary

At the end of QC script:

# Calculate percentage remaining
probe_tracking$pct_remaining <- round(100 * probe_tracking$n_probes / probe_tracking$n_probes[1], 1)

cat("\n=== Probe Tracking Summary ===\n\n")
print(probe_tracking)

# Save to file
write.csv(probe_tracking, file.path(tables_dir, "probe_tracking.csv"), row.names = FALSE)

Expected Output

Step	Probes	Removed	% Remaining	Description
1. Raw data loaded	29,377	0	100.0%	Total probes in raw data file
2. Expression matrix	29,377	0	100.0%	Probes with valid signal columns
3. Detection filter	20,233	9,144	68.9%	Probes detected in >= 3 samples

NA Value Checking

IMPORTANT: Check for NA values in both signal and detection p-value matrices early in the QC script. This identifies data quality issues before they cause downstream problems.

# Check for NAs in signal data
na_signal <- is.na(expr_mat)
n_na_signal_total <- sum(na_signal)
n_na_signal_by_sample <- colSums(na_signal)
n_na_signal_by_probe <- rowSums(na_signal)

cat(sprintf("Signal data - Total NA values: %d (%.2f%%)\n",
            n_na_signal_total,
            100 * n_na_signal_total / length(expr_mat)))

if (n_na_signal_total > 0) {
  cat("NA counts by sample (signal):\n")
  for (i in seq_along(sample_names)) {
    if (n_na_signal_by_sample[i] > 0) {
      cat(sprintf
                  sample_namesi
                  n_na_signal_by_samplei
                    n_na_signal_by_samplei  nrowexpr_mat
    
  
  catsprintf
              n_na_signal_by_probe  
                n_na_signal_by_probe    nrowexpr_mat



na_detect  detect_mat
n_na_detect_total  na_detect



na_summary  data.frame
  sample  sample_names
  na_signal_count  n_na_signal_by_sample
  na_signal_pct    n_na_signal_by_sample  nrowexpr_mat
  na_detect_count  n_na_detect_by_sample
  na_detect_pct    n_na_detect_by_sample  nrowdetect_mat

write.csvna_summary file.pathtables_dir  row.names  

NA Visualization

If NA values are present, generate diagnostic plots:

if (n_na_signal_total > 0 || n_na_detect_total > 0) {
  # Bar plot comparing signal vs detection NAs by sample
  na_summary_long <- na_summary %>%
    pivot_longer(cols = c(na_signal_count, na_detect_count),
                 names_to = "data_type", values_to = "na_count")

  p_na <- ggplot(na_summary_long, aes(x = sample, y = na_count, fill = data_type)) +
    geom_bar(stat = "identity", position = "dodge") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust  

  ggsavefile.pathplots_dir  p_na width   height  

  
   n_na_signal_by_probe     
    
  

Handling NAs in Downstream Analysis

• Use na.rm = TRUE for summary statistics: colMeans(expr_mat, na.rm = TRUE)
• Clean data before correlation/PCA: expr_clean <- expr_mat[complete.cases(expr_mat), ]
• Document the number of probes removed due to NA values in the probe tracking table

Platform Detection and Data Loading

Affymetrix Arrays (CEL files)

Use affy and affyPLM packages:

library(affy)
library(affyPLM)

# Load CEL files
raw_data <- ReadAffy(filenames = cel_files)

# Basic QC
pm_data <- pm(raw_data)
mm_data <- mm(raw_data)
avg_pm <- colMeans(pm_data)

# MAS5 present/absent calls
calls <- mas5calls(raw_data)
percent_present <- colSums(exprs(calls) == "P") / nrow(exprs(calls)) * 100

# PLM QC (NUSE and RLE)
plm_fit <- fitPLM(raw_data)
nuse_vals <- NUSE(plm_fit, type = "values")
rle_vals <- RLE(plm_fit, type = "values")

# RMA normalization
eset  rmaraw_data

Annotation packages by platform:

• GPL570 (HG-U133_Plus_2): hgu133plus2.db
• GPL96 (HG-U133A): hgu133a.db
• GPL97 (HG-U133B): hgu133b.db

Illumina BeadChip Arrays

Use limma with detection p-values:

library(limma)

# Data typically has AVG_Signal and Detection Pval columns
expr_mat <- as.matrix(data[, signal_cols])
detect_mat <- as.matrix(data[, pval_cols])

# Filter by detection (p < 0.05 in >= 3 samples)
probes_detected <- rowSums(detect_mat < 0.05) >= 3
expr_filtered <- expr_mat[probes_detected, ]

# Log2 transform and quantile normalize
expr_log2 <- log2(expr_filtered + 1)
expr_norm <- normalizeBetweenArrays(expr_log2, method = "quantile")

Annotation packages by platform:

• GPL14951 (HumanHT-12 v4): illuminaHumanv4.db
• GPL10558 (HumanHT-12 v4.0): illuminaHumanv4.db
• GPL6884 (HumanHT-12 v3): illuminaHumanv3.db

Quality Control Plots

Generate these QC plots (numbered sequentially):

1. Raw intensity boxplot - Distribution per sample
2. Raw intensity density - Overlapping density curves
3. Percent present/detected - Bar plot per sample
4. Average intensity - Bar plot per sample
5. Sample correlation heatmap - Hierarchical clustering
6. PCA plot - PC1 vs PC2 colored by condition
7. QC metrics by condition - Boxplots comparing groups

Post-normalization visualization:

• Overall expression heatmap (e.g., 19b) - Top 2000 most variable genes
  • Median-centered, ±3 SD clipping for color utilization
  • 1 - Pearson correlation distance, average linkage for gene clustering
  • Samples ordered by condition (not clustered)
  • Green-black-red color scheme with condition color bars
• Sample dendrogram (e.g., 19c) - Hierarchical clustering of all samples
  • Uses all genes with 1 - Pearson correlation distance
  • Average linkage, horizontal layout
  • Branches and labels colored by sample condition

Affymetrix-specific:

• RNA degradation plot - 5' to 3' bias
• NUSE boxplot - Outliers if median > 1.05
• RLE boxplot - Outliers if |median| > 0.1

Outlier Detection Criteria

Metric	Threshold	Action
NUSE median	> 1.05	Flag as outlier
NUSE IQR	> 0.15	Flag as outlier
RLE median		value
RLE IQR	> 0.3	Flag as outlier
Intensity z-score		z
PCA extreme	Visual outlier	Investigate

After excluding outliers: Re-run normalization (especially RMA which uses quantile normalization).

Statistical Modeling with limma

Paired Design (same subjects, multiple conditions)

For matched samples (e.g., lesional vs non-lesional from same patient):

# Create factors
tissue <- factor(sample_info$condition, levels = c("NonLesional", "Lesional"))
patient <- factor(sample_info$patient_id)

# Design with patient as fixed effect (n >= 10 patients)
design <- model.matrix(~ patient + tissue)

# Fit model
fit <- lmFit(expr_filtered, design)

# Contrasts (tissue effects are relative to baseline)
contrast_matrix <- makeContrasts(
  Lesional_vs_NonLesional = tissueLesional,
  levels = design
)

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

Unpaired Design (independent groups)

For independent samples:

# Create factor
group <- factor(sample_info$condition, levels = c("Control", "Treatment"))

# Cell means model (no intercept)
design <- model.matrix(~ 0 + group)
colnames(design) <- levels(group)

# Fit model
fit <- lmFit(expr_filtered, design)

# Define contrast
contrast_matrix <- makeContrasts(
  Treatment_vs_Control = Treatment - Control,
  levels = design
)

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

Interaction Terms (genotype x treatment)

For factorial designs:

# Create combined factor
group <- factor(paste(sample_info$genotype, sample_info$treatment, sep = "_"))
design <- model.matrix(~ 0 + group)
colnames(design) <- levels(group)

# Define biologically meaningful contrasts
contrast_matrix <- makeContrasts(
  Treatment_in_WT = WT_Treated - WT_Control,
  Treatment_in_Mut = Mut_Treated - Mut_Control,
  Interaction = (Mut_Treated - Mut_Control) - (WT_Treated - WT_Control),
  levels = design
)

Model Diagnostics

Always generate these diagnostic plots:

1. P-value histogram - Should show uniform + spike near 0
2. MA plot - Check for intensity-dependent bias
3. Mean-variance trend (plotSA) - Validates eBayes assumptions

# P-value histogram
png(file.path(plots_dir, "03_pvalue_histogram.png"), width = 800, height = 600)
hist(results$P.Value, breaks = 50,
     main = "P-value Distribution", xlab = "P-value",
     col = "steelblue", border = "white")
abline(h = nrow(results) / 50, col = "red", lty = 2)
legend("topright", "Expected under null", col = "red", lty = 2)
dev.off()


pngfile.pathplots_dir  width   height  
plotMAfit2 coef   main  
dev.off


pngfile.pathplots_dir  width   height  
plotSAfit2 main  
dev.off

Volcano Plot Best Practices

Important: The y-axis should show FDR (adjusted p-value), not raw p-value:

p_volcano <- ggplot(results, aes(x = logFC, y = -log10(adj.P.Val), color = significance)) +
  geom_point(alpha = 0.5, size = 1.5) +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "gray40") +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "gray40") +
  labs(title = "Volcano Plot",
       x  
       y    

Gene ID Mapping

Using the collapse_to_gene Utility Function

Use the shared utility function scripts/utils/collapse_to_gene.R for probe-to-gene mapping. This function:

1. Queries the annotation database for the specified gene ID type
2. Removes probes that map to multiple gene IDs (ambiguous mapping)
3. Keeps the most significant probe per gene (by P.Value, |logFC| tiebreaker)

source(here("scripts", "utils", "collapse_to_gene.R"))

# Basic usage - collapse to gene symbols
de_annotated <- collapse_to_gene(
  de_results,
  annot_pkg = "illuminaHumanv4.db",  # Platform-specific annotation package
  id_column = "SYMBOL",
  probe_column = "ProbeID",
  keytype = "PROBEID",
  collapse = TRUE,    # Set FALSE to keep all probes with annotations
  verbose = TRUE
) %>%
  rename(gene_symbol = SYMBOL) %>%
  relocate(ProbeID, gene_symbol)

# Get probe-level annotations without collapsing
de_probe_level <- collapse_to_gene(
  de_results,
  annot_pkg = "illuminaHumanv4.db",
  id_column = "SYMBOL",
  probe_column  
  collapse    

Annotation Packages by Platform

Platform	GPL ID	Annotation Package
Illumina HumanHT-12 v4	GPL14951, GPL10558	illuminaHumanv4.db
Illumina HumanHT-12 v3	GPL6884	illuminaHumanv3.db
Affymetrix HG-U133 Plus 2	GPL570	hgu133plus2.db
Affymetrix HG-U133A	GPL96	hgu133a.db

Manual Probe-to-Gene Mapping (Alternative)

If not using the utility function:

library(AnnotationDbi)

# Get mapping (platform-specific)
mapping <- AnnotationDbi::select(
  platform.db,
  keys = probe_ids,
  columns = c("SYMBOL", "ENTREZID", "ENSEMBL", "GENENAME"),
  keytype = "PROBEID"
)

# Join with DE results
de_annotated <- de_results |>
  left_join(mapping, by = c("ProbeID" = "PROBEID"))

# Collapse to unique gene symbols (keep most significant probe)
de_clean_symbol <- de_annotated %>%
  filter(!is.na(gene_symbol)) %>%
  arrange(P.Value, -abs(logFC)) 
  filterduplicatedgene_symbol 
  arrangeadj.P.Val

Output Format

Differential Expression Results

Required columns in output CSV:

Column	Description
ensembl_id	Ensembl gene ID (primary identifier)
gene_symbol	HGNC gene symbol
logFC	Log2 fold change
P.Value	Raw p-value
adj.P.Val	FDR-adjusted p-value
ProbeID	Original probe/probeset ID
AveExpr	Average expression
t	Moderated t-statistic
B	Log-odds of differential expression

File Naming Convention

Use numbered prefixes (01_, 02_, etc.) to indicate the script that generated each file:

results/tables/
├── 01_sample_metadata.csv                 # From 01_load_and_qc.R
├── 01_qc_metrics.csv
├── 01_expression_matrix_raw.csv
├── 01_probe_tracking.csv
├── 01_na_value_summary.csv                # NA diagnostics
├── 02_expression_matrix_normalized.csv    # From 02_normalize.R
├── 02_normalization_summary.csv
├── 03_DE_[Comparison].csv                 # From 03_differential_expression.R
├── 04_DE_[Comparison]_annotated.csv       # From 04_annotate_genes.R
├── 04_DE_[Comparison]_clean_symbol.csv    # Unique gene symbols
├── 04_genes_upregulated.csv
├── 04_genes_downregulated.csv
├── 05_GO_BP_upregulated.csv               # From 05_go_enrichment.R
├── 05_GO_BP_downregulated.csv
├── 05_KEGG_upregulated.csv
├── 05_KEGG_downregulated.csv
└── 05_GSEA_GO_BP.csv

results/plots/
├── 01_raw_intensity_boxplot.png           # From 01_load_and_qc.R
├── 01_raw_intensity_density.png
├── 01_PCA_raw.png
├── 01_sample_correlation_heatmap.png
├── 01_na_comparison.png                   # If NAs present
├── 02_normalized_boxplot.png              # From 02_normalize.R
├── 02_PCA_normalized.png
├── 03_pvalue_histogram.png                # From 03_differential_expression.R
├── 03_MA_plot.png
├── 03_mean_variance_trend.png
├── 03_volcano_plot.png
├── 03_heatmap_top_DE.png
├── 05_GO_BP_combined.png                  # From 05_go_enrichment.R
├── 05_KEGG_combined.png
├── 05_GSEA_GO_BP_combined.png
├── 05_GSEA_GO_BP_dotplot.png
└── 05_GSEA_GO_BP_ridgeplot.png

Note: The _clean_symbol.csv file contains one row per gene (collapsed from multiple probesets), selected by lowest P.Value. Use this for downstream analysis.

Progressive Documentation (analysis.md)

Update analysis.md after EACH analysis step with:

# [GSE#] Microarray Analysis Report

## Study Information

| Field | Value |
|-------|-------|
| **GEO Accession** | GSE##### |
| **Platform** | GPL### [Platform Name] |
| **Analysis Date** | YYYY-MM-DD |

## Methods

### Software Environment

| Package | Version | Purpose |
|---------|---------|---------|
| limma | X.X.X | Differential expression |
| affy | X.X.X | CEL file processing |

### Key Parameters
- [List all thresholds and parameters used]

## Results

### Quality Control
- [Summarize QC findings]
- [Document any outliers removed]

### Differential Expression
| Comparison | Total DEGs | Up | Down |
|------------|------------|-----|------|
| ... | ... | ... | ... |

### Key Findings
- [Biological interpretation]

## Output Files
[List all generated files]

Plot Documentation (plots/plots.md)

Create results/plots/plots.md with embedded image references so plots render in markdown viewers:

# Plot Documentation - [GSE#]

## Quality Control

### 01_raw_intensity_boxplot.png
Raw probe intensity distribution by sample before normalization.

![Raw intensity boxplot](./01_raw_intensity_boxplot.png)

### 02_raw_intensity_density.png
Density curves of raw intensities colored by condition.

![Raw intensity density](./02_raw_intensity_density.png)

[Continue for all plots with image references...]

## Differential Expression

### 15_volcano_plot.png
Volcano plot showing -log10(p-value) vs log2 fold change.
- Red: Upregulated (FDR < 0.05, |logFC| > 1)
- Blue: Downregulated (FDR < 0.05, |logFC| > 1)

![Volcano plot](./15_volcano_plot.png)

Important: Use ./ prefix for relative paths (e.g., ![description](./filename.png)) to ensure images render correctly in VS Code and other markdown viewers.

GO and KEGG Enrichment Analysis

Use clusterProfiler for pathway analysis. Run both GO Biological Process and KEGG pathway enrichment for both upregulated AND downregulated genes in each comparison.

Background Gene List

Important: Use all genes that passed the filtering step as background, loaded from the annotation file:

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

# Load background from annotation file (all filtered genes)
annotation_df <- read.csv(file.path(tables_dir, "probeset_annotation.csv"))
background_genes <- annotation_df %>%
  filter(!is.na(entrez_id)) %>%
  pull(entrez_id) %>%
  unique() %>%
  as.character()

GO Biological Process Enrichment (Both Directions)

# Get Entrez IDs for UPREGULATED genes
sig_up <- de_results %>%
  filter(adj.P.Val < 0.05 & logFC > 0 & !is.na(entrez_id)) %>%
  pull(entrez_id) %>%
  unique() %>%
  as.character()

# Get Entrez IDs for DOWNREGULATED genes
sig_down <- de_results %>%
  filter(adj.P.Val < 0.05 & logFC < 0 & !is.na(entrez_id)) %>%
  pull(entrez_id) %>%
  unique() %>%
  as.character()

# GO enrichment - Upregulated
ego_up <- enrichGO(
  gene = sig_up,
  universe  background_genes
  OrgDb  org.Hs.eg.db
  keyType  
  ont  
  pAdjustMethod  
  pvalueCutoff  
  readable  

ego_up_simp  simplifyego_up cutoff   by  


ego_down  enrichGO
  gene  sig_down
  universe  background_genes
  OrgDb  org.Hs.eg.db
  keyType  
  ont  
  pAdjustMethod  
  pvalueCutoff  
  readable  

ego_down_simp  simplifyego_down cutoff   by  

KEGG Pathway Enrichment (Both Directions)

# KEGG - Upregulated
kegg_up <- enrichKEGG(
  gene = sig_up,
  universe = background_genes,
  organism = "hsa",
  keyType = "ncbi-geneid",
  pAdjustMethod = "BH",
  pvalueCutoff = 0.05
)

# KEGG - Downregulated
kegg_down <- enrichKEGG(
  gene = sig_down,
  universe = background_genes,
  organism = "hsa",
  keyType = "ncbi-geneid",
  pAdjustMethod = "BH",
  pvalueCutoff = 0.05
)

Combined Plots with Fold Enrichment

Use the utility functions from scripts/utils/enrichment_dotplot.R for creating combined enrichment plots:

source(here("scripts/utils/enrichment_dotplot.R"))

# Helper function to prepare ORA data with fold enrichment
prepare_ora_data <- function(enrichResult, direction_label) {
  if (is.null(enrichResult) || nrow(as.data.frame(enrichResult)) == 0) return(NULL)

  df <- as.data.frame(enrichResult)
  df <- df %>%
    mutate(
      gene_ratio_num = as.numeric(sub("/.*", "", GeneRatio)),
      gene_ratio_denom = as.numeric(sub(".*/", "", GeneRatio)
      bg_ratio_num  sub  BgRatio
      bg_ratio_denom  sub  BgRatio
      GeneRatioValue  gene_ratio_num  gene_ratio_denom
      BgRatioValue  bg_ratio_num  bg_ratio_denom
      FoldEnrichment  GeneRatioValue  BgRatioValue
      direction  direction_label
    
  df



go_up_ora  prepare_ora_dataego_up_simp 
go_down_ora  prepare_ora_dataego_down_simp 
go_combined_ora  bind_rowsgo_up_ora go_down_ora


p_go_combined  plot_ora_dotplot_combined
  go_combined_ora
  n_show  
  fdr_cutoff  
  title  
  width  
  height  

ggsavefile.pathplots_dir  p_go_combined
       width  p_go_combined 
       height  p_go_combined  dpi  


kegg_up_ora  prepare_ora_datakegg_up 
kegg_down_ora  prepare_ora_datakegg_down 
kegg_combined_ora  bind_rowskegg_up_ora kegg_down_ora

p_kegg_combined  plot_ora_dotplot_combined
  kegg_combined_ora
  n_show  
  fdr_cutoff  
  title  

ggsavefile.pathplots_dir  p_kegg_combined
       width  p_kegg_combined 
       height  p_kegg_combined  dpi  

GSEA Combined Plots

For GSEA results, use plot_gsea_dotplot_combined():

# Run GSEA
gsea_go <- gseGO(geneList = gene_list, OrgDb = org.Hs.eg.db, ...)
gsea_df <- as.data.frame(gsea_go)

# Create combined GSEA plot
p_gsea_combined <- plot_gsea_dotplot_combined(
  gsea_df,
  n_show = 10,
  fdr_cutoff = 0.05,
  title = "GSEA GO BP: Comparison",
  size_label = "Set Size"
)
ggsave(file.path(plots_dir, "05_GSEA_GO_BP_combined.png"), p_gsea_combined,
       width = attr(p_gsea_combined, "width"),
       height = attr(p_gsea_combined, "height"), dpi = 300)

Output Files

For each comparison, produce:

• ORA combined plots:
  • 05_GO_BP_combined.png - Up + down on same figure
  • 05_KEGG_combined.png - Up + down on same figure
• GSEA plots:
  • 05_GSEA_GO_BP_combined.png - Activated + suppressed
  • 05_GSEA_GO_BP_dotplot.png - Standard clusterProfiler dotplot
  • 05_GSEA_GO_BP_ridgeplot.png - Distribution ridge plot
  • 05_GSEA_GO_BP_running.png - Running enrichment for top terms
• CSV tables:
  • 05_GO_BP_upregulated.csv / 05_GO_BP_downregulated.csv
  • 05_KEGG_upregulated.csv / 05_KEGG_downregulated.csv
  • 05_GSEA_GO_BP.csv

Contrast Selection Guidelines

Only extract biologically meaningful contrasts:

Good contrasts:

• Treatment vs Control
• Disease vs Healthy
• Lesional vs Non-Lesional (paired)
• Time point comparisons in time series

Avoid:

• Redundant comparisons (if A-B and B-C, you may not need A-C)
• Comparisons without biological interpretation
• All pairwise when only specific comparisons matter

Common Issues and Solutions

Issue	Solution
Few DEGs with FDR < 0.05	Use nominal P < 0.05 for exploratory GO analysis
Batch effects visible in PCA	Consider ComBat or limma's removeBatchEffect
Low annotation rate	Try alternative ID mapping sources (biomaRt)
High proportion of outlier samples	Review sample preparation notes; may need to exclude

Utility Functions

The project includes reusable utility functions in scripts/utils/:

Function	File	Description
collapse_to_gene()	collapse_to_gene.R	Maps probes to genes, removes ambiguous mappings, collapses to unique genes
plot_ora_dotplot_combined()	enrichment_dotplot.R	Combined ORA dotplot for up/downregulated genes
plot_gsea_dotplot_combined()	enrichment_dotplot.R	Combined GSEA dotplot for activated/suppressed terms

Always source these utilities using here():

source(here("scripts", "utils", "collapse_to_gene.R"))
source(here("scripts", "utils", "enrichment_dotplot.R"))

Reference Files

See these files in the project for complete examples:

• analyses/Rashighi_GSE53146/ - Illumina unpaired design (most current patterns)
• analyses/Regazzetti_GSE65127/ - Affymetrix paired design
• analyses/Natarajan_GSE75819/ - Illumina paired design