Adjust primary DE model for sex covariate (~ condition + gender)

The primary DESeq2 model now controls for sex as an additive covariate,
addressing the main biological weakness flagged in audit: unadjusted
~ condition model despite having sex metadata available.

Changes:
- 01_qc.R: filter samples with NA gender before balanced sampling,
  update VST design formula
- 03_deseq2.R: ~ condition → ~ condition + gender
- 09_sensitivity_analysis.R: match model, filter NA genders from
  full cohort for consistency

apeglm shrinkage continues to target the condition coefficient only.
The sex interaction model (script 11) remains as a separate complementary
analysis testing condition × gender interaction effects.

Note: committed results were generated with ~ condition and need
regeneration. CI pipeline-rebuild will regenerate from scratch.

Author: Ekin-Kahraman

README.md

@@ -33,7 +33,7 @@ GSE152075 (n=484, GEO)
  02 PCA ─────────── VST (blind=TRUE) → PCA for sample-level exploratory analysis
     │
     ▼
- 03 DE ──────────── Balanced subset (n=60) → DESeq2 (~ condition) → apeglm shrinkage
+ 03 DE ──────────── Balanced subset (n=60) → DESeq2 (~ condition + gender) → apeglm shrinkage
     │
     ├──→ 04 Sensitivity ─── Full cohort (n=484) DE → concordance check (99.7% sign agreement)
     ├──→ 05 Diagnostics ─── Cook's distance, dispersion, MA, volcano, scree
@@ -48,7 +48,7 @@ GSE152075 (n=484, GEO)
 ## Methods Overview
 
 - Bulk RNA-seq preprocessing and quality control
-- Differential expression modelling with DESeq2 plus `apeglm` log2 fold-change shrinkage
+- Differential expression modelling with DESeq2 (`~ condition + gender`) plus `apeglm` log2 fold-change shrinkage
 - Variance-stabilising transformation (VST) for visualisation
 - Functional enrichment analysis (GO/KEGG)
 - Full-cohort robustness benchmark against the balanced subset

scripts/01_qc.R

@@ -62,6 +62,14 @@ if (any(duplicated(rownames(counts)))) {
 saveRDS(counts, "data/counts_qc.rds")
 saveRDS(metadata, "data/metadata_qc.rds")
 
+# Filter to samples with non-NA gender for covariate-adjusted model
+has_gender <- !is.na(metadata$gender)
+if (sum(!has_gender) > 0) {
+  message("Excluding ", sum(!has_gender), " samples with missing gender annotation")
+  metadata <- metadata[has_gender, , drop = FALSE]
+  counts <- counts[, rownames(metadata), drop = FALSE]
+}
+
 set.seed(analysis_config$sample_seed)
 n <- analysis_config$samples_per_group
 pos_samples <- rownames(metadata)[metadata$condition == "positive"]
@@ -92,7 +100,7 @@ message("Final: ", sum(metadata$condition == "negative"), " neg, ",
 saveRDS(counts, "data/counts_clean.rds")
 saveRDS(metadata, "data/metadata_clean.rds")
 
-dds <- DESeqDataSetFromMatrix(counts, metadata, design = ~ condition)
+dds <- DESeqDataSetFromMatrix(counts, metadata, design = ~ condition + gender)
 vsd <- varianceStabilizingTransformation(dds, blind = TRUE)
 saveRDS(vsd, "data/vst_data.rds")
 

scripts/03_deseq2.R

@@ -21,7 +21,7 @@ metadata <- readRDS("data/metadata_clean.rds")
 
 message("Testing ", nrow(counts), " genes across ", ncol(counts), " samples")
 
-dds <- DESeqDataSetFromMatrix(counts, metadata, design = ~ condition)
+dds <- DESeqDataSetFromMatrix(counts, metadata, design = ~ condition + gender)
 dds$condition <- relevel(dds$condition, ref = "negative")
 
 dds <- DESeq(dds)

scripts/09_sensitivity_analysis.R

@@ -29,6 +29,15 @@ if (!requireNamespace("apeglm", quietly = TRUE)) {
 
 counts_qc <- readRDS("data/counts_qc.rds")
 metadata_qc <- readRDS("data/metadata_qc.rds")
+
+# Filter to samples with non-NA gender (matching the balanced subset model)
+has_gender <- !is.na(metadata_qc$gender)
+if (sum(!has_gender) > 0) {
+  message("Excluding ", sum(!has_gender), " full-cohort samples with missing gender annotation")
+  metadata_qc <- metadata_qc[has_gender, , drop = FALSE]
+  counts_qc <- counts_qc[, rownames(metadata_qc), drop = FALSE]
+}
+
 balanced_raw <- read.csv("results/tables/deseq2_results.csv", stringsAsFactors = FALSE)
 balanced_shrunk <- read.csv("results/tables/deseq2_results_shrunken.csv", stringsAsFactors = FALSE)
 go_df <- read.csv("results/tables/go_biological_process.csv", stringsAsFactors = FALSE)
@@ -38,7 +47,7 @@ counts_balanced <- readRDS("data/counts_clean.rds")
 keep_genes_full <- rowSums(cpm(counts_qc) >= analysis_config$gene_cpm_cutoff) >= analysis_config$gene_min_samples
 counts_full <- counts_qc[keep_genes_full, , drop = FALSE]
 
-dds_full <- DESeqDataSetFromMatrix(counts_full, metadata_qc, design = ~ condition)
+dds_full <- DESeqDataSetFromMatrix(counts_full, metadata_qc, design = ~ condition + gender)
 dds_full$condition <- relevel(dds_full$condition, ref = "negative")
 dds_full <- DESeq(dds_full)