1 Configuration

Set the two *_OVERRIDE variables to the exact column names you saw in the inspection output.

ID_OVERRIDE   <- NULL   # e.g. "Assay"  (metadata column holding the sample/run IDs)
GENO_OVERRIDE <- NULL   # e.g. "Factor Value[genotype]"

# Regexes used to recognise the two groups. KO is tested FIRST because
# "knockout" also contains the substring "ko".
KO_RE <- "knock[ _-]?out|^ko$|mutant|^-/-$"
WT_RE <- "wild[ _-]?type|^wt$|control"

read_tab <- function(path) {
  if (grepl("\\.csv$", path, ignore.case = TRUE))
    read.csv(path,  check.names = FALSE, stringsAsFactors = FALSE)
  else
    read.delim(path, check.names = FALSE, stringsAsFactors = FALSE)
}
pick_col <- function(cands, nms, what) {
  hit <- cands[cands %in% nms]
  if (!length(hit)) stop("Could not find ", what, ". Columns are: ",
                         paste(nms, collapse = ", "))
  hit[1]
}

2 Problem 1 — Load the count matrix correctly

raw <- read_tab(params$counts_file)
names(raw) <- trimws(names(raw))

cat("Raw file dimensions (rows x cols):", dim(raw), "\n")
Raw file dimensions (rows x cols): 61860 8 
knitr::kable(head(raw[, 1:min(6, ncol(raw)), drop = FALSE]),
             caption = "First rows of the raw file (note the annotation columns)")
First rows of the raw file (note the annotation columns)
Gene ID Gene Name ERR1736465 ERR1736466 ERR1736467 ERR1736468
ENSG00000000003 TSPAN6 434 829 334 779
ENSG00000000005 TNMD 0 0 0 0
ENSG00000000419 DPM1 2202 3461 1627 4241
ENSG00000000457 SCYL3 45 72 31 119
ENSG00000000460 C1orf112 320 558 235 516
ENSG00000000938 FGR 0 10 0 0 
id_col   <- pick_col(c("Gene.ID","Gene ID","GeneID","gene_id"), names(raw), "Gene ID column")
name_col <- pick_col(c("Gene.Name","Gene Name","GeneName","gene_name","symbol"),
                     names(raw), "Gene Name column")

sample_cols <- setdiff(names(raw), c(id_col, name_col))
num_ok <- vapply(raw[, sample_cols, drop = FALSE], is.numeric, logical(1))
if (any(!num_ok))
  stop("Non-numeric columns among presumed samples: ",
       paste(sample_cols[!num_ok], collapse = ", "))

counts <- as.matrix(raw[, sample_cols, drop = FALSE])
if (anyNA(counts))     stop("NA values in the count matrix.")
if (any(counts < 0))   stop("Negative counts.")
if (any(counts != round(counts))) counts <- round(counts)   # DESeq2 needs integers
storage.mode(counts) <- "integer"

gene_id <- as.character(raw[[id_col]])
bad_id  <- is.na(gene_id) | gene_id == ""
if (any(bad_id)) gene_id[bad_id] <- paste0("noID_", seq_len(sum(bad_id)))
if (anyDuplicated(gene_id)) {
  message("Duplicate Gene IDs found -> make.unique() applied.")
  gene_id <- make.unique(gene_id)
}
rownames(counts) <- gene_id

anno <- data.frame(Gene.ID   = gene_id,
                   Gene.Name = as.character(raw[[name_col]]),
                   stringsAsFactors = FALSE)

Deliverable

cat("dim(counts):", dim(counts), "\n")
dim(counts): 61860 6 
cat("First 3 rownames (Gene.ID):", head(rownames(counts), 3), "\n")
First 3 rownames (Gene.ID): ENSG00000000003 ENSG00000000005 ENSG00000000419 
cat("First 3 colnames (samples):", head(colnames(counts), 3), "\n")
First 3 colnames (samples): ERR1736465 ERR1736466 ERR1736467 
str(counts[1:3, 1:min(3, ncol(counts))])
 int [1:3, 1:3] 434 0 2202 829 0 3461 334 0 1627
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:3] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419"
  ..$ : chr [1:3] "ERR1736465" "ERR1736466" "ERR1736467"

The count object is now a pure integer matrix (genes x samples); gene symbols live only in anno, so no character column can leak into the model.

3 Problem 2 — Load metadata and confirm sample matching

meta_raw <- read_tab(params$meta_file)
cat("Metadata dimensions:", dim(meta_raw), "\n")
Metadata dimensions: 6 16 
knitr::kable(head(meta_raw[, 1:min(6, ncol(meta_raw)), drop = FALSE]),
             caption = "First rows of the metadata")
First rows of the metadata
Run Sample Characteristic[cell line] Sample Characteristic Ontology Term[cell line] Sample Characteristic[cell type] Sample Characteristic Ontology Term[cell type] Sample Characteristic[disease]
ERR1736465 HS578T http://www.ebi.ac.uk/efo/EFO_0001192 epithelial cell http://purl.obolibrary.org/obo/CL_0000066 breast carcinoma
ERR1736466 HS578T http://www.ebi.ac.uk/efo/EFO_0001192 epithelial cell http://purl.obolibrary.org/obo/CL_0000066 breast carcinoma
ERR1736467 HS578T http://www.ebi.ac.uk/efo/EFO_0001192 epithelial cell http://purl.obolibrary.org/obo/CL_0000066 breast carcinoma
ERR1736468 HS578T http://www.ebi.ac.uk/efo/EFO_0001192 epithelial cell http://purl.obolibrary.org/obo/CL_0000066 breast carcinoma
ERR1736469 HS578T http://www.ebi.ac.uk/efo/EFO_0001192 epithelial cell http://purl.obolibrary.org/obo/CL_0000066 breast carcinoma
ERR1736470 HS578T http://www.ebi.ac.uk/efo/EFO_0001192 epithelial cell http://purl.obolibrary.org/obo/CL_0000066 breast carcinoma 
# Which metadata column holds the sample IDs? (highest overlap with count columns)
if (is.null(ID_OVERRIDE)) {
  overlap <- vapply(meta_raw, function(x) sum(colnames(counts) %in% as.character(x)),
                    integer(1))
  print(sort(overlap, decreasing = TRUE)[1:min(5, length(overlap))])
  if (max(overlap) == 0)
    stop("No metadata column matches colnames(counts). Set ID_OVERRIDE manually.")
  ID_COL <- names(meta_raw)[which.max(overlap)]
} else ID_COL <- ID_OVERRIDE
                                           Run 
                                             6 
              Sample Characteristic[cell line] 
                                             0 
Sample Characteristic Ontology Term[cell line] 
                                             0 
              Sample Characteristic[cell type] 
                                             0 
Sample Characteristic Ontology Term[cell type] 
                                             0 
cat("Using sample-ID column:", ID_COL, "\n")
Using sample-ID column: Run 
if (anyDuplicated(meta_raw[[ID_COL]])) stop("Duplicated sample IDs in metadata.")
rownames(meta_raw) <- as.character(meta_raw[[ID_COL]])

if (!all(colnames(counts) %in% rownames(meta_raw)))
  stop("Samples missing from metadata: ",
       paste(setdiff(colnames(counts), rownames(meta_raw)), collapse = ", "))

metadata <- meta_raw[match(colnames(counts), rownames(meta_raw)), , drop = FALSE]

Deliverable — alignment check

alignment <- data.frame(
  position   = seq_len(ncol(counts)),
  counts_col = colnames(counts),
  meta_row   = rownames(metadata),
  identical  = colnames(counts) == rownames(metadata)
)
knitr::kable(alignment, caption = "Sample alignment check")
Sample alignment check
position counts_col meta_row identical
1 ERR1736465 ERR1736465 TRUE
2 ERR1736466 ERR1736466 TRUE
3 ERR1736467 ERR1736467 TRUE
4 ERR1736468 ERR1736468 TRUE
5 ERR1736469 ERR1736469 TRUE
6 ERR1736470 ERR1736470 TRUE 
stopifnot(identical(colnames(counts), rownames(metadata)))
cat("ALL SAMPLES ALIGNED:", all(alignment$identical), "\n")
ALL SAMPLES ALIGNED: TRUE

4 Problem 3 — Condition factor with WT as reference

if (is.null(GENO_OVERRIDE)) {
  cand <- grep("genotype|condition|group|treatment|phenotype|factor value",
               names(metadata), ignore.case = TRUE, value = TRUE)
  if (!length(cand)) stop("No condition-like column. Set GENO_OVERRIDE manually.")
  score <- vapply(cand, function(cn) {
    v <- tolower(trimws(as.character(metadata[[cn]])))
    sum(grepl(KO_RE, v) | grepl(WT_RE, v))
  }, integer(1))
  print(score)
  GENO_COL <- cand[which.max(score)]
} else GENO_COL <- GENO_OVERRIDE
              Sample Characteristic[genotype] 
                                            6 
Sample Characteristic Ontology Term[genotype] 
                                            0 
                       Factor Value[genotype] 
                                            6 
         Factor Value Ontology Term[genotype] 
                                            0 
cat("Using genotype column:", GENO_COL, "\n")
Using genotype column: Sample Characteristic[genotype] 
print(table(metadata[[GENO_COL]], useNA = "ifany"))

    Snai1 knockout wild type genotype 
                 3                  3 
v <- tolower(trimws(as.character(metadata[[GENO_COL]])))
cond <- ifelse(grepl(KO_RE, v), "KO",
        ifelse(grepl(WT_RE, v), "WT", NA_character_))   # KO tested first!
if (anyNA(cond))
  stop("Unmapped genotype labels: ", paste(unique(v[is.na(cond)]), collapse = ", "),
       " -> hard-code the mapping here.")

metadata$Condition <- factor(cond, levels = c("WT", "KO"))

Deliverable

print(levels(metadata$Condition))
[1] "WT" "KO"
print(table(metadata$Condition))

WT KO 
 3  3

Interpretation. WT is the reference (baseline / intercept) level, so the model coefficient Condition_KO_vs_WT is the log2 fold change of KO relative to WT. A positive log2FoldChange therefore means higher expression in the knockout.

5 Problem 4 — Build the DESeqDataSet

stopifnot(!anyNA(metadata$Condition), nlevels(droplevels(metadata$Condition)) == 2)
dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData   = metadata,
                              design    = ~ Condition)
dds
class: DESeqDataSet 
dim: 61860 6 
metadata(1): version
assays(1): counts
rownames(61860): ENSG00000000003 ENSG00000000005 ... ENSG00000290165
  ENSG00000290166
rowData names(0):
colnames(6): ERR1736465 ERR1736466 ... ERR1736469 ERR1736470
colData names(17): Run Sample Characteristic[cell line] ... Analysed
  Condition

6 Problem 5 — Filter lowly expressed genes

n_before <- nrow(dds)
keep     <- rowSums(counts(dds)) > params$filter_main
dds      <- dds[keep, ]
n_after  <- nrow(dds)

knitr::kable(data.frame(
  Genes_before = n_before, Threshold = paste0("rowSums > ", params$filter_main),
  Genes_after  = n_after,  Removed   = n_before - n_after,
  Percent_kept = round(100 * n_after / n_before, 1)),
  caption = "Filtering summary")
Filtering summary
Genes_before Threshold Genes_after Removed Percent_kept
61860 rowSums > 10 14548 47312 23.5

Why filter. Genes with almost no reads carry essentially no evidence for a difference; they cannot reach significance but still enter the Benjamini–Hochberg correction, inflating the multiple-testing burden and reducing power for real hits. They also destabilise dispersion estimation and shrinkage. Removing them is independent of the test statistic (it uses only overall counts, not the group labels), so it does not bias the p-values.

Stricter alternative worth mentioning: keep genes with rowSums(counts >= 10) >= min(table(metadata$Condition)), i.e. require a usable count in at least one whole group rather than one read summed over all samples.

7 Problem 6 — Run DESeq2 and extract results

dds <- DESeq(dds)
resultsNames(dds)
[1] "Intercept"          "Condition_KO_vs_WT"
saveRDS(dds, "results/rds/dds_filterB.rds")
res_05  <- results(dds, contrast = c("Condition","KO","WT"), alpha = 0.05)
res_1e5 <- results(dds, contrast = c("Condition","KO","WT"), alpha = 1e-5)

cat("=========== summary(res), alpha = 0.05 ===========\n"); summary(res_05)
=========== summary(res), alpha = 0.05 ===========

out of 14548 with nonzero total read count
adjusted p-value < 0.05
LFC > 0 (up)       : 2155, 15%
LFC < 0 (down)     : 2068, 14%
outliers [1]       : 32, 0.22%
low counts [2]     : 1411, 9.7%
(mean count < 3)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
cat("=========== summary(res), alpha = 1e-5 ===========\n"); summary(res_1e5)
=========== summary(res), alpha = 1e-5 ===========

out of 14548 with nonzero total read count
adjusted p-value < 1e-05
LFC > 0 (up)       : 769, 5.3%
LFC < 0 (down)     : 538, 3.7%
outliers [1]       : 32, 0.22%
low counts [2]     : 2257, 16%
(mean count < 5)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
knitr::kable(data.frame(
  alpha        = c(0.05, 1e-5),
  non_NA_padj  = c(sum(!is.na(res_05$padj)), sum(!is.na(res_1e5$padj))),
  sig_at_alpha = c(sum(res_05$padj  < 0.05,  na.rm = TRUE),
                   sum(res_1e5$padj < 1e-5,  na.rm = TRUE))),
  caption = "Effect of alpha")
Effect of alpha
alpha non_NA_padj sig_at_alpha
0.05 13105 4223
0.00 12259 1307

What changes and why. The raw p-values are identical in both calls — alpha does not touch the Wald test. alpha is the target FDR that DESeq2 uses to optimise the independent filtering threshold (it picks the mean-count cut-off maximising rejections at that FDR) and it sets the cut-off reported by summary(). At alpha = 1e-5 the optimised filter is more aggressive, so a different set of genes receives a non-NA padj, and BH is applied over that different set — hence slightly different adjusted p-values and far fewer “significant” calls. 0.05 is the conventional, reportable choice; 1e-5 is an ultra-conservative screen.

8 Problem 7 — Results data frame + gene names

res_df <- as.data.frame(res_05)
res_df$Gene.ID <- rownames(res_df)
res_df <- merge(res_df, anno, by = "Gene.ID", all.x = TRUE, sort = FALSE)
res_df <- res_df[order(res_df$padj, na.last = TRUE), ]
stopifnot("Gene.Name" %in% names(res_df), nrow(res_df) == nrow(res_05))

Deliverable

knitr::kable(head(res_df[, c("Gene.ID","Gene.Name","baseMean",
                             "log2FoldChange","pvalue","padj")]),
             digits = 4, caption = "head() of the merged results table")
head() of the merged results table
Gene.ID Gene.Name baseMean log2FoldChange pvalue padj
6885 ENSG00000151632 AKR1C2 1342.2 -6.164 0 0
2676 ENSG00000108691 CCL2 2157.8 8.493 0 0
1561 ENSG00000096696 DSP 798.2 -4.300 0 0
4076 ENSG00000124212 PTGIS 467.7 -5.393 0 0
12769 ENSG00000225383 SFTA1P 773.0 4.750 0 0
7848 ENSG00000163017 ACTG2 1886.6 3.004 0 0 
cat("Gene.Name present:", "Gene.Name" %in% names(res_df), "\n")
Gene.Name present: TRUE

9 Problem 8 — Top 10 most significant genes

sig   <- res_df[!is.na(res_df$padj), ]
sig   <- sig[order(sig$padj, -abs(sig$log2FoldChange)), ]
top10 <- head(sig[, c("Gene.Name","log2FoldChange","padj")], 10)
knitr::kable(top10, digits = 4, row.names = FALSE,
             caption = "Top 10 genes by adjusted p-value (ties broken by |log2FC|)")
Top 10 genes by adjusted p-value (ties broken by |log2FC|)
Gene.Name log2FoldChange padj
AKR1C2 -6.164 0
CCL2 8.493 0
DSP -4.300 0
PTGIS -5.393 0
SFTA1P 4.750 0
ACTG2 3.004 0
NPR3 5.662 0
AKR1C1 -2.929 0
OPHN1 -3.563 0
RPL22L1 -2.468 0

10 Problem 9 — Genes of interest

genes_to_check <- c("THY1","SFMBT2","PASD1","SNAI1")
goi <- res_df[res_df$Gene.Name %in% genes_to_check,
              c("Gene.Name","baseMean","log2FoldChange","pvalue","padj")]
goi$Call <- ifelse(is.na(goi$padj), "not tested (filtered)",
            ifelse(goi$padj < 0.05 & goi$log2FoldChange > 0, "Up in KO",
            ifelse(goi$padj < 0.05 & goi$log2FoldChange < 0, "Down in KO", "not significant")))
knitr::kable(goi, digits = 4, row.names = FALSE, caption = "Candidate genes")
Candidate genes
Gene.Name baseMean log2FoldChange pvalue padj Call
THY1 3146.01 12.372 0 0 Up in KO
SFMBT2 260.22 -11.807 0 0 Down in KO
PASD1 251.23 11.076 0 0 Up in KO
SNAI1 53.89 -4.155 0 0 Down in KO 
missing <- setdiff(genes_to_check, res_df$Gene.Name)
if (length(missing)) {
  cat("Not in the results table:", paste(missing, collapse = ", "), "\n")
  ids <- anno$Gene.ID[anno$Gene.Name %in% missing]
  if (length(ids))
    print(data.frame(Gene.Name = anno$Gene.Name[match(ids, anno$Gene.ID)],
                     total_counts = rowSums(counts[ids, , drop = FALSE])))
}

Interpretation (edit to match your numbers). Read the sign of log2FoldChange: positive = up in KO, negative = down in KO — but only call a gene differentially expressed when padj < 0.05. A large fold change with padj near 1 (or a small baseMean) is noise, not a result. Genes absent from the table were removed by the count filter, i.e. they are essentially not expressed in this dataset.

11 Problem 10 — Volcano plot

vol <- EnhancedVolcano(res_df,
  lab = res_df$Gene.Name, x = "log2FoldChange", y = "padj",
  pCutoff = 0.05, FCcutoff = 1,
  pointSize = 2, labSize = 3.5, drawConnectors = TRUE, max.overlaps = 20,
  title = "KO vs WT", subtitle = "DESeq2 | FDR 0.05, |log2FC| > 1",
  caption = paste0(sum(!is.na(res_df$padj)), " genes tested"))
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
ℹ The deprecated feature was likely used in the EnhancedVolcano package.
  Please report the issue to the authors.
This warning is displayed once per session.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.
Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
ℹ Please use the `linewidth` argument instead.
ℹ The deprecated feature was likely used in the EnhancedVolcano package.
  Please report the issue to the authors.
This warning is displayed once per session.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.
print(vol)

ggplot2::ggsave("figures/volcano_KO_vs_WT.png", vol, width = 9, height = 7, dpi = 300)
systemfonts and textshaping have been compiled with different versions of Freetype. Because of this, textshaping will not use the font cache provided by systemfonts
EnhancedVolcano(res_df,
  lab = res_df$Gene.Name, x = "log2FoldChange", y = "pvalue",
  pCutoff = 1e-5, FCcutoff = 2,
  pointSize = 2, labSize = 3.5,
  title = "KO vs WT (strict: raw p < 1e-5, |log2FC| > 2)")

Note: with y = "padj" the pCutoff line is on the adjusted-p scale; genes with NA padj are not drawn. Tightening pCutoff/FCcutoff reduces label clutter.

12 Reusable analysis function (Problems 11–12)

run_deseq <- function(counts, metadata, filter_min = 10, alpha = 0.05, anno, quiet = TRUE) {
  keep <- rowSums(counts) > filter_min
  d <- DESeqDataSetFromMatrix(counts[keep, , drop = FALSE],
                              colData = metadata, design = ~ Condition)
  d <- tryCatch(DESeq(d, quiet = quiet), error = function(e) {
    message("Default size factors failed (", conditionMessage(e), ") -> poscounts.")
    DESeq(estimateSizeFactors(d, type = "poscounts"), quiet = quiet)
  })
  r <- as.data.frame(results(d, contrast = c("Condition","KO","WT"), alpha = alpha))
  r$Gene.ID <- rownames(r)
  r <- merge(r, anno, by = "Gene.ID", all.x = TRUE, sort = FALSE)
  list(dds = d, df = r,
       n_tested = nrow(d),
       n_padj   = sum(!is.na(r$padj)),
       n_sig    = sum(r$padj < 0.05, na.rm = TRUE))
}
top_names <- function(x, n = 10) head(x$df$Gene.Name[order(x$df$padj, na.last = TRUE)], n)

13 Problem 11 — Diagnose a sample mismatch bug

ok_run <- run_deseq(counts, metadata, params$filter_main, params$alpha_main, anno)

set.seed(42)
meta_bad <- metadata[sample(nrow(metadata)), , drop = FALSE]
# Simulate the *silent* bug: labels look correct, the content is scrambled.
# (Some DESeq2 versions refuse mismatched rownames outright -- that is the safety
#  check working. Relabelling reproduces the real-world mistake of loading
#  metadata without IDs and assuming the order.)
rownames(meta_bad) <- colnames(counts)
bad_run <- run_deseq(counts, meta_bad, params$filter_main, params$alpha_main, anno)

knitr::kable(data.frame(rank = 1:10,
                        correct = top_names(ok_run),
                        shuffled = top_names(bad_run)),
             caption = "Top 10 genes: correct vs shuffled metadata")
Top 10 genes: correct vs shuffled metadata
rank correct shuffled
1 AKR1C2 SERP1
2 CCL2 GADD45GIP1
3 DSP KNTC1
4 PTGIS MT2A
5 SFTA1P SQSTM1
6 ACTG2 RPL36
7 NPR3 PDLIM7
8 AKR1C1 AIFM1
9 OPHN1 RPL26L1
10 RPL22L1 DNAJB12 
knitr::kable(data.frame(
  run = c("correct","shuffled"),
  genes_tested = c(ok_run$n_tested, bad_run$n_tested),
  significant_padj_0.05 = c(ok_run$n_sig, bad_run$n_sig),
  top10_overlap = c(length(intersect(top_names(ok_run), top_names(bad_run))),
                    length(intersect(top_names(ok_run), top_names(bad_run)))),
  median_dispersion = c(median(mcols(ok_run$dds)$dispersion, na.rm = TRUE),
                        median(mcols(bad_run$dds)$dispersion, na.rm = TRUE))),
  caption = "Consequences of misalignment")
Consequences of misalignment
run genes_tested significant_padj_0.05 top10_overlap median_dispersion
correct 14548 4223 0 0.0340
shuffled 14548 0 0 0.0725 
p_ok  <- plotPCA(vst(ok_run$dds,  blind = FALSE), intgroup = "Condition") +
         ggplot2::ggtitle("Correct alignment")
using ntop=500 top features by variance
Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
ℹ Please use tidy evaluation idioms with `aes()`.
ℹ See also `vignette("ggplot2-in-packages")` for more information.
ℹ The deprecated feature was likely used in the DESeq2 package.
  Please report the issue to the authors.
This warning is displayed once per session.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.
p_bad <- plotPCA(vst(bad_run$dds, blind = FALSE), intgroup = "Condition") +
         ggplot2::ggtitle("Shuffled labels")
using ntop=500 top features by variance
print(p_ok); print(p_bad)

Explanation (edit to match your numbers). DESeq2 does not error, it silently fits the wrong model: each column of counts is now paired with another sample’s genotype. The “WT” and “KO” groups become mixtures of the same biological populations, so true between-group differences are absorbed into within-group variance. Dispersion estimates rise, Wald statistics shrink, and the number of significant genes collapses; the few remaining “top hits” are driven by individual outlier samples rather than a consistent group effect, so the top-10 lists barely overlap. The PCA makes it visible: samples separate by genotype when aligned and are intermingled when shuffled. Alignment must be enforced by name (match() + stopifnot(identical(...))), never assumed from row order.

14 Problem 12 — Compare filtering strategies

filters <- c(A = 0, B = 10, C = 50)
runs <- lapply(filters, function(f)
  run_deseq(counts, metadata, filter_min = f, alpha = params$alpha_main, anno = anno))
saveRDS(runs, "results/rds/filter_runs.rds")

summary_tbl <- data.frame(
  filter          = names(filters),
  threshold       = paste0("rowSums > ", unname(filters)),
  genes_tested    = vapply(runs, `[[`, integer(1), "n_tested"),
  genes_with_padj = vapply(runs, `[[`, integer(1), "n_padj"),
  sig_padj_0.05   = vapply(runs, `[[`, integer(1), "n_sig"))
knitr::kable(summary_tbl, row.names = FALSE, caption = "Filtering comparison")
Filtering comparison
filter threshold genes_tested genes_with_padj sig_padj_0.05
A rowSums > 0 19580 13489 4129
B rowSums > 10 14548 13105 4223
C rowSums > 50 11698 11666 4238 
sigsets <- lapply(runs, function(r) r$df$Gene.ID[which(r$df$padj < 0.05)])
ex <- setdiff(sigsets$A, sigsets$C)
if (!length(ex)) ex <- setdiff(sigsets$C, sigsets$A)

if (length(ex)) {
  ex <- ex[1]
  ex_tbl <- do.call(rbind, lapply(names(runs), function(nm) {
    i <- match(ex, runs[[nm]]$df$Gene.ID)
    if (is.na(i))
      data.frame(filter = nm, Gene.Name = anno$Gene.Name[match(ex, anno$Gene.ID)],
                 baseMean = NA, log2FoldChange = NA, pvalue = NA, padj = NA)
    else
      data.frame(filter = nm, Gene.Name = runs[[nm]]$df$Gene.Name[i],
                 baseMean = runs[[nm]]$df$baseMean[i],
                 log2FoldChange = runs[[nm]]$df$log2FoldChange[i],
                 pvalue = runs[[nm]]$df$pvalue[i], padj = runs[[nm]]$df$padj[i])
  }))
  knitr::kable(ex_tbl, digits = 5, row.names = FALSE,
               caption = "Example gene whose significance depends on the filter")
} else cat("No gene changed significance between filters A and C.")
Example gene whose significance depends on the filter
filter Gene.Name baseMean log2FoldChange pvalue padj
A NADK 5.847 4.596 0.00349 0.01474
B NADK 5.847 4.597 0.00295 0.01213
C NADK NA NA NA NA

Interpretation (edit to match your numbers). Loosening the filter (A) admits thousands of low-count genes; they add tests to the BH correction without adding discoveries, so borderline genes lose significance. Tightening it (C) removes noisy genes and can rescue borderline genes by shrinking the test family, but it also discards genuinely low-expressed genes — some real biology is lost. Note the two distinct numbers: genes tested (rows entering DESeq2) versus genes with a padj (those surviving DESeq2’s own independent filtering / outlier removal). The differences between A/B/C are a multiple-testing effect, not new evidence about individual genes; rowSums > 10 is a reasonable, pre-specified compromise.

15 Problem 13 — Collaborator results table

final <- ok_run$df[, c("Gene.ID","Gene.Name","baseMean","log2FoldChange",
                       "lfcSE","stat","pvalue","padj")]
sg <- !is.na(final$padj) & final$padj < 0.05          # NA guard is essential
final$Direction <- ifelse(sg & final$log2FoldChange > 0, "Up_in_KO",
                   ifelse(sg & final$log2FoldChange < 0, "Down_in_KO", "Not_Significant"))
final <- final[order(final$padj, na.last = TRUE), ]

write.csv(final, "results/DESeq2_WT_vs_KO_results.csv", row.names = FALSE)

shortlist <- head(final[final$Direction != "Not_Significant" &
                        abs(final$log2FoldChange) > 1, ], 20)
write.csv(shortlist, "results/DESeq2_top_hits_shortlist.csv", row.names = FALSE)

knitr::kable(head(final, 10), digits = 4, row.names = FALSE,
             caption = "First 10 rows of the exported CSV")
First 10 rows of the exported CSV
Gene.ID Gene.Name baseMean log2FoldChange lfcSE stat pvalue padj Direction
ENSG00000151632 AKR1C2 1342.2 -6.164 0.1991 -30.96 0 0 Down_in_KO
ENSG00000108691 CCL2 2157.8 8.493 0.2830 30.01 0 0 Up_in_KO
ENSG00000096696 DSP 798.2 -4.300 0.1601 -26.86 0 0 Down_in_KO
ENSG00000124212 PTGIS 467.7 -5.393 0.2013 -26.79 0 0 Down_in_KO
ENSG00000225383 SFTA1P 773.0 4.750 0.1811 26.23 0 0 Up_in_KO
ENSG00000163017 ACTG2 1886.6 3.004 0.1299 23.13 0 0 Up_in_KO
ENSG00000113389 NPR3 548.9 5.662 0.2479 22.84 0 0 Up_in_KO
ENSG00000187134 AKR1C1 2255.7 -2.929 0.1291 -22.69 0 0 Down_in_KO
ENSG00000079482 OPHN1 4581.3 -3.563 0.1580 -22.56 0 0 Down_in_KO
ENSG00000163584 RPL22L1 1564.2 -2.468 0.1111 -22.21 0 0 Down_in_KO

Deliverable

table(final$Direction)

     Down_in_KO Not_Significant        Up_in_KO 
           2068           10325            2155 
cat("CSV written to:", normalizePath("results/DESeq2_WT_vs_KO_results.csv"), "\n")
CSV written to: /home/neolith/Data/Course/R bulk scRNAseq/rnaseq-wt-ko/results/DESeq2_WT_vs_KO_results.csv

For collaborators: keep the Gene.ID column and import the CSV via Data → From Text/CSV (set Gene columns to “Text”), otherwise Excel converts symbols such as SEPT2 / MARCH1 into dates. NA in padj = gene not tested (filtered or outlier).

17 Reproducibility

sessionInfo()
R version 4.4.3 (2025-02-28)
Platform: x86_64-conda-linux-gnu
Running under: Ubuntu 24.04.5 LTS

Matrix products: default
BLAS/LAPACK: /home/neolith/anaconda3/envs/rstudio_env/lib/libopenblasp-r0.3.33.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Asia/Dhaka
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] EnhancedVolcano_1.24.0      ggrepel_0.9.8              
 [3] ggplot2_4.0.4               DESeq2_1.46.0              
 [5] SummarizedExperiment_1.36.0 Biobase_2.66.0             
 [7] MatrixGenerics_1.18.1       matrixStats_1.5.0          
 [9] GenomicRanges_1.58.0        GenomeInfoDb_1.42.3        
[11] IRanges_2.40.1              S4Vectors_0.44.0           
[13] BiocGenerics_0.52.0        

loaded via a namespace (and not attached):
 [1] gtable_0.3.6            xfun_0.59               bslib_0.11.0           
 [4] lattice_0.23-1          numDeriv_2016.8-1.1     vctrs_0.7.3            
 [7] tools_4.4.3             generics_0.1.4          parallel_4.4.3         
[10] tibble_3.3.1            pkgconfig_2.0.3         pheatmap_1.0.13        
[13] Matrix_1.7-6            RColorBrewer_1.1-3      S7_0.2.2               
[16] lifecycle_1.0.5         GenomeInfoDbData_1.2.13 compiler_4.4.3         
[19] farver_2.1.2            textshaping_1.0.4       codetools_0.2-20       
[22] htmltools_0.5.9         sass_0.4.10             yaml_2.3.12            
[25] pillar_1.11.1           crayon_1.5.3            jquerylib_0.1.4        
[28] MASS_7.3-66             BiocParallel_1.40.2     DelayedArray_0.32.0    
[31] cachem_1.1.0            emdbook_1.3.14          abind_1.4-8            
[34] bdsmatrix_1.3-7         tidyselect_1.2.1        locfit_1.5-9.12        
[37] digest_0.6.39           mvtnorm_1.4-2           dplyr_1.2.1            
[40] labeling_0.4.3          fastmap_1.2.0           grid_4.4.3             
[43] cli_3.6.6               SparseArray_1.6.2       magrittr_2.0.5         
[46] S4Arrays_1.6.0          withr_3.0.3             scales_1.4.0           
[49] UCSC.utils_1.2.0        rmarkdown_2.32          XVector_0.46.0         
[52] httr_1.4.9              otel_0.3.0              ragg_1.5.2             
[55] coda_0.19-4.1           evaluate_1.0.5          knitr_1.52             
[58] bbmle_1.0.26            rlang_1.3.0             Rcpp_1.1.2             
[61] glue_1.8.1              apeglm_1.28.0           rstudioapi_0.19.0      
[64] jsonlite_2.0.0          plyr_1.8.9              R6_2.6.1               
[67] systemfonts_1.3.2       zlibbioc_1.52.0