Module 6: Statistical Preprocessing & Cancer Signatures: Imputation, ComBat & Volcano Prioritization
1. The Missing Value Dilemma: MCAR vs. MNAR
Mass spectrometry data contains widespread missing values ($10 - 40\%$) arising from two distinct biological and analytical mechanisms:
- Missing Completely at Random (MCAR): Stochastic issues such as chromatographic bubble dropouts or competitive ionization suppression. MCAR values occur uniformly across all abundance ranges and are imputed using $k$-Nearest Neighbors ($k$-NN) or random forest imputation.
- Missing Not at Random (MNAR): Left-censored values caused by peptide or metabolite concentrations falling below the instrument's lower limit of detection (LOD). A biomarker present in glioblastoma but undetectable in healthy controls is genuine biological signal, not analytical error. Left-censored values must be imputed using left-shifted truncated normal distributions (MinProb or QRILC).
2. Normalization & ComBat Batch Effect Adjustment
When patient cohorts are acquired over several days or flow cell column changes, systematic batch drift corrupts biomarker discovery. Empirical Bayes framework (ComBat) adjusts for batch variance while preserving true cancer phenotypic variance:
$$Y_{ijg} = \alpha_g + X_j \beta_g + \gamma_{ig} + \delta_{ig} \epsilon_{ijg}$$Where $\gamma_{ig}$ represents additive batch effects and $\delta_{ig}$ represents multiplicative variance scaling.
3. Differential Analysis & Volcano Plot Visualization
Candidate biomarkers are prioritized based on statistical effect size and significance:
# Complete Limma Differential Expression Pipeline in R
library(limma)
# Log2 transformation and median centering
log_intensities <- log2(feature_matrix + 1)
norm_intensities <- normalizeMedianValues(log_intensities)
# Design matrix with batch and condition covariates
design <- model.matrix(~ 0 + pheno$cohort)
colnames(design) <- c("Control", "Glioblastoma")
# Fit linear model
fit <- lmFit(norm_intensities, design)
contrast_matrix <- makeContrasts(GBM_vs_Ctrl = Glioblastoma - Control, levels = design)
fit_contrast <- contrasts.fit(fit, contrast_matrix)
fit_ebayes <- eBayes(fit_contrast)
# Extract top table with Benjamini-Hochberg FDR correction
results_table <- topTable(fit_ebayes, adjust="BH", number=Inf)
results_table$Significant <- results_table$adj.P.Val < 0.05 & abs(results_table$logFC) >= 1.0
cat("Identified", sum(results_table$Significant), "significant candidate biomarkers (FDR < 0.05, |log2FC| >= 1).\n")
head(results_table[results_table$Significant, ])Top prioritized hits in our WHO CNS5 brain tumor cohort include hyper-accumulated oncometabolite $2\text{-HG}$ ($\log_2\text{FC} = +5.82$), elevated glycolytic enzyme Hexokinase-2 ($HK2$, $\log_2\text{FC} = +2.41$), and suppressed neuroprotective metabolites (N-acetylaspartate, $\log_2\text{FC} = -3.15$).
Common Pitfalls & Troubleshooting Checklist
- Java Heap Exhaustion (GATK / Picard): If processes terminate with
java.lang.OutOfMemoryError, set your environment variable explicitly:export _JAVA_OPTIONS="-Xmx16g -Djava.io.tmpdir=/tmp". - Index Timestamp Mismatch: If tools complain that a BAM index is missing or out of date, regenerate with
samtools index -@ 4 sample.bam. - Conda Channel Conflicts: Avoid ABI incompatibilities between bioconda and conda-forge by setting strict priority:
conda config --set channel_priority strict. - Chromosome Naming Contamination: Always verify whether reference FASTA headers use
chr1or1; divergent naming causes silent alignment drops.