16_mock_3sps: sample QC + analysis of the within-batch mock contrast (yH001, yH149, yH714)
Hubert Kicinski
2026-08-18
Purpose
The 2026 run adds a within-batch t0 (mock.0m) for the
three comparative species, so the mock-vs-t0 handling contrast is
finally prep-matched. Before using the data, I ask two things: how large
is the 2026-vs-2023 prep batch effect (can we reasonably compare), and
are the new libraries technically sound (MultiQC results indicate:
yes).
library(DESeq2); library(edgeR); library(dplyr); library(tidyr); library(tibble)
library(stringr); library(readr); library(purrr); library(ggplot2); library(ggprism)
library(pheatmap); library(RColorBrewer); library(matrixStats); library(here)
base_dir <- here("16_mock_3sps")
new_dir <- file.path(base_dir, "inputs", "counts_new")
og_dir <- file.path(base_dir, "inputs", "counts_OG")
out_dir <- file.path(base_dir, "outputs")
fig_dir <- file.path(base_dir, "figures")
species <- list(
yH001 = list(new = "yH001_mock_gene_counts.tsv", og = "CGLA-salmon.merged.gene_counts.tsv", sp = "C. glabrata"),
yH149 = list(new = "yH149_mock_gene_counts.tsv", og = "KLAC-salmon.merged.gene_counts.tsv", sp = "K. lactis"),
yH714 = list(new = "yH714_mock_gene_counts.tsv", og = "CALB-salmon.merged.gene_counts.tsv", sp = "C. albicans")
)
batch_cols <- c("2023" = "#377EB8", "2026" = "#E41A1C")
cond_cols <- c(t0 = "#4D4D4D", mock = "#E41A1C", noPi = "#377EB8")
dir_cols <- c(down = "#4A90D9", up = "#C0392B")
save_fig <- function(p, fname, width = 9, height = 6) {
ggsave(file.path(fig_dir, fname), plot = p, width = width, height = height, dpi = 300, bg = "white")
}
collapse_lanes <- function(m) {
if (!any(grepl("\\.L00\\d$", colnames(m)))) return(m)
base <- sub("\\.L00\\d$", "", colnames(m))
out <- sapply(unique(base), function(s) rowSums(m[, base == s, drop = FALSE]))
rownames(out) <- rownames(m)
out
}
read_counts <- function(path) {
raw <- read.delim(path, row.names = 1, check.names = FALSE)
raw <- raw[, !colnames(raw) %in% "gene_name", drop = FALSE]
m <- collapse_lanes(as.matrix(raw))
storage.mode(m) <- "double"
m
}
build_species <- function(code) {
cfg <- species[[code]]
new_m <- read_counts(file.path(new_dir, cfg$new))
new_meta <- tibble(sample = colnames(new_m)) |>
mutate(tp = as.integer(str_extract(sample, "(?<=\\.)\\d+(?=m\\.rep)")),
rep = paste0("rep", as.integer(str_extract(sample, "(?<=rep)\\d+"))),
rep_orig = NA_character_,
batch = "2026",
cond = if_else(tp == 0, "t0", "mock"))
stopifnot(!any(is.na(new_meta$tp)))
og_m <- read_counts(file.path(og_dir, cfg$og))
og_meta <- tibble(sample = colnames(og_m)) |>
mutate(tp = as.integer(str_extract(sample, "(?<=\\.t)\\d{4}")),
rep_orig = str_extract(sample, "(?<=\\.)b\\d+"),
batch = "2023",
cond = if_else(tp == 0, "t0", "noPi")) |>
mutate(rep = paste0("rep", dense_rank(rep_orig)), .by = tp)
stopifnot(!any(is.na(og_meta$tp)))
common <- intersect(rownames(new_m), rownames(og_m))
stopifnot(length(common) > 1000)
colnames(new_m) <- paste0("new__", colnames(new_m))
colnames(og_m) <- paste0("og__", colnames(og_m))
new_meta$sample <- paste0("new__", new_meta$sample)
og_meta$sample <- paste0("og__", og_meta$sample)
counts <- round(cbind(new_m[common, , drop = FALSE], og_m[common, , drop = FALSE]))
meta <- bind_rows(new_meta, og_meta) |>
mutate(label = gsub("^(new__|og__)", "", sample),
group = paste(batch, cond, tp, sep = "_"),
batch = factor(batch, levels = c("2023", "2026")),
cond = factor(cond, levels = c("t0", "mock", "noPi"))) |>
arrange(batch, cond, tp, rep) |>
mutate(label = factor(label, levels = unique(label)))
keep <- rowSums(counts >= 10) >= 2
dge <- calcNormFactors(DGEList(counts = counts[keep, ]), method = "TMM")
list(code = code, sp = cfg$sp,
counts = counts, counts_filt = counts[keep, ],
lcpm = cpm(dge, log = TRUE, prior.count = 2),
meta = meta)
}
dat <- map(names(species), build_species) |> set_names(names(species))
design_tbl <- map_dfr(dat, ~ count(.x$meta, batch, cond, tp) |> mutate(code = .x$code, sp = .x$sp)) |>
select(code, sp, batch, cond, tp, n)
write_csv(design_tbl, file.path(out_dir, "00_design_summary.csv"))
#Note to self: Took these from my HE26069 analysis
#adjusted a few functions to fit calb, klac, and cgla
t0_means <- map_dfr(dat, function(d) {
s_new <- as.character(d$meta$sample[d$meta$batch == "2026" & d$meta$cond == "t0"])
s_og <- as.character(d$meta$sample[d$meta$batch == "2023" & d$meta$cond == "t0"])
if (!length(s_new) || !length(s_og)) return(NULL)
tibble(code = d$code, sp = d$sp, gene_id = rownames(d$lcpm),
t0_2026 = rowMeans(d$lcpm[, s_new, drop = FALSE]),
t0_2023 = rowMeans(d$lcpm[, s_og, drop = FALSE]))
})
t0_cor <- t0_means |>
summarise(.by = c(code, sp),
n_genes = n(),
pearson = cor(t0_2023, t0_2026),
spearman = cor(t0_2023, t0_2026, method = "spearman"),
median_abs_diff = median(abs(t0_2026 - t0_2023)))
write_csv(t0_cor, file.path(out_dir, "01_t0_batch_correlation.csv"))
pca_dat <- map(dat, function(d) {
pca <- prcomp(t(d$lcpm), center = TRUE, scale. = FALSE)
ve <- (pca$sdev^2) / sum(pca$sdev^2) * 100
df <- as.data.frame(pca$x[, 1:2]) |>
rownames_to_column("sample") |>
left_join(d$meta |> select(sample, label, batch, cond, tp), by = "sample") |>
mutate(sp = d$sp, code = d$code, tp_f = factor(tp, levels = sort(unique(tp))))
list(df = df, ve = ve)
})
pca_df <- map_dfr(pca_dat, "df")
write_csv(pca_df |> select(code, sp, sample, label, batch, cond, tp, PC1, PC2),
file.path(out_dir, "01_joint_pca_coordinates.csv"))
t0_pca_dist <- map_dfr(pca_dat, function(x) {
d <- x$df
a <- d |> filter(batch == "2026", cond == "t0") |> summarise(PC1 = mean(PC1), PC2 = mean(PC2))
b <- d |> filter(batch == "2023", cond == "t0") |> summarise(PC1 = mean(PC1), PC2 = mean(PC2))
ref <- d |> filter(batch == "2023") |> summarise(spread = max(dist(cbind(PC1, PC2))))
tibble(sp = d$sp[1], code = d$code[1],
t0_gap = sqrt((a$PC1 - b$PC1)^2 + (a$PC2 - b$PC2)^2),
timecourse_spread = ref$spread,
gap_frac = t0_gap / ref$spread)
})
write_csv(t0_pca_dist, file.path(out_dir, "01_t0_pca_distance.csv"))
mock_lfc <- map_dfr(dat, function(d) {
keep <- d$meta$batch == "2026"
cd <- d$meta[keep, ] |> as.data.frame()
rownames(cd) <- cd$sample
cd$group <- factor(paste0(cd$cond, cd$tp), levels = c("t00", paste0("mock", sort(unique(cd$tp[cd$tp > 0])))))
mm <- d$counts[, as.character(cd$sample), drop = FALSE]
mm <- mm[rowSums(mm >= 10) >= 2, ]
dds <- DESeq(DESeqDataSetFromMatrix(round(mm), cd, ~ group), quiet = TRUE)
map_dfr(setdiff(levels(cd$group), "t00"), function(g) {
r <- results(dds, contrast = c("group", g, "t00"),
independentFiltering = FALSE, cooksCutoff = FALSE)
as_tibble(as.data.frame(r), rownames = "gene_id") |>
transmute(code = d$code, sp = d$sp, contrast = g,
gene_id, lfc = log2FoldChange, padj)
})
})
write_csv(mock_lfc, file.path(out_dir, "01_mock_vs_t0_lfc.csv"))
mock_summary <- mock_lfc |>
summarise(.by = c(code, sp, contrast),
n_tested = n(),
up = sum(coalesce(padj < 0.05, FALSE) & lfc > 1),
down = sum(coalesce(padj < 0.05, FALSE) & lfc < -1))
write_csv(mock_summary, file.path(out_dir, "01_mock_vs_t0_deg_counts.csv"))
libsize <- map_dfr(dat, function(d) {
tibble(code = d$code, sp = d$sp, sample = colnames(d$counts),
total = colSums(d$counts)) |>
left_join(d$meta |> select(sample, label, batch, cond, tp), by = "sample")
})
write_csv(libsize, file.path(out_dir, "01_library_size.csv"))
lcpm_long <- map_dfr(dat, function(d) {
as.data.frame(d$lcpm) |>
rownames_to_column("gene_id") |>
pivot_longer(-gene_id, names_to = "sample", values_to = "lcpm") |>
left_join(d$meta |> select(sample, label, batch, cond, tp), by = "sample") |>
mutate(code = d$code, sp = d$sp)
})
rle_long <- map_dfr(dat, function(d) {
m <- d$lcpm
for (g in unique(d$meta$group)) {
cols <- as.character(d$meta$sample[d$meta$group == g])
if (length(cols) < 2) next
m[, cols] <- d$lcpm[, cols, drop = FALSE] - rowMedians(d$lcpm[, cols, drop = FALSE])
}
as.data.frame(m) |>
rownames_to_column("gene_id") |>
pivot_longer(-gene_id, names_to = "sample", values_to = "rle") |>
left_join(d$meta |> select(sample, label, batch, cond, tp, group), by = "sample") |>
mutate(code = d$code, sp = d$sp)
})
rep_sd <- map_dfr(dat, function(d) {
map_dfr(unique(d$meta$group), function(g) {
cols <- as.character(d$meta$sample[d$meta$group == g])
if (length(cols) < 2) return(NULL)
tibble(code = d$code, sp = d$sp, group = g,
batch = d$meta$batch[match(g, d$meta$group)],
sd = rowSds(d$lcpm[, cols, drop = FALSE]))
})
})
write_csv(rep_sd |> summarise(.by = c(code, sp, group, batch), median_sd = median(sd)),
file.path(out_dir, "01_replicate_sd_median.csv"))
print(design_tbl, n = Inf)## # A tibble: 46 × 6
## code sp batch cond tp n
## <chr> <chr> <fct> <fct> <int> <int>
## 1 yH001 C. glabrata 2023 t0 0 2
## 2 yH001 C. glabrata 2023 noPi 15 2
## 3 yH001 C. glabrata 2023 noPi 30 2
## 4 yH001 C. glabrata 2023 noPi 45 2
## 5 yH001 C. glabrata 2023 noPi 60 2
## 6 yH001 C. glabrata 2023 noPi 90 2
## 7 yH001 C. glabrata 2023 noPi 120 2
## 8 yH001 C. glabrata 2023 noPi 180 2
## 9 yH001 C. glabrata 2023 noPi 240 2
## 10 yH001 C. glabrata 2023 noPi 360 2
## 11 yH001 C. glabrata 2023 noPi 480 2
## 12 yH001 C. glabrata 2026 t0 0 2
## 13 yH001 C. glabrata 2026 mock 30 2
## 14 yH001 C. glabrata 2026 mock 60 2
## 15 yH149 K. lactis 2023 t0 0 2
## 16 yH149 K. lactis 2023 noPi 15 2
## 17 yH149 K. lactis 2023 noPi 30 2
## 18 yH149 K. lactis 2023 noPi 45 2
## 19 yH149 K. lactis 2023 noPi 60 2
## 20 yH149 K. lactis 2023 noPi 90 2
## 21 yH149 K. lactis 2023 noPi 120 1
## 22 yH149 K. lactis 2023 noPi 150 2
## 23 yH149 K. lactis 2023 noPi 180 2
## 24 yH149 K. lactis 2023 noPi 210 2
## 25 yH149 K. lactis 2023 noPi 240 2
## 26 yH149 K. lactis 2023 noPi 360 2
## 27 yH149 K. lactis 2023 noPi 480 2
## 28 yH149 K. lactis 2026 t0 0 2
## 29 yH149 K. lactis 2026 mock 30 2
## 30 yH149 K. lactis 2026 mock 60 2
## 31 yH714 C. albicans 2023 t0 0 2
## 32 yH714 C. albicans 2023 noPi 15 2
## 33 yH714 C. albicans 2023 noPi 30 2
## 34 yH714 C. albicans 2023 noPi 45 2
## 35 yH714 C. albicans 2023 noPi 60 2
## 36 yH714 C. albicans 2023 noPi 90 2
## 37 yH714 C. albicans 2023 noPi 120 2
## 38 yH714 C. albicans 2023 noPi 150 2
## 39 yH714 C. albicans 2023 noPi 180 2
## 40 yH714 C. albicans 2023 noPi 210 2
## 41 yH714 C. albicans 2023 noPi 240 2
## 42 yH714 C. albicans 2023 noPi 360 2
## 43 yH714 C. albicans 2023 noPi 480 2
## 44 yH714 C. albicans 2026 t0 0 2
## 45 yH714 C. albicans 2026 mock 30 2
## 46 yH714 C. albicans 2026 mock 60 2> helpers processed data properly. I have a long format design tibble that stores the data for me by key. This will make downstream analysis easier for me.
01. Prep-batch magnitude: 2026 t0 vs 2023 t0
Now, our t0s should IDEALLY be the same. However, due to different extraction and library-prep batch, there should be reasonable expectation for batch effects. Nevertheless, for this analysis to work as intended, these batch effects should be neglible such that noPi 30 min and noPi 60 min can reasonably be compared to mock 60 min and 30 min in the 3 species. For this, we test how large the 2026-vs-2023 offset is.
print(t0_cor)## # A tibble: 3 × 6
## code sp n_genes pearson spearman median_abs_diff
## <chr> <chr> <int> <dbl> <dbl> <dbl>
## 1 yH001 C. glabrata 5231 0.929 0.888 0.482
## 2 yH149 K. lactis 5105 0.914 0.879 0.508
## 3 yH714 C. albicans 6096 0.891 0.856 0.695p_t0_scatter <- ggplot(t0_means, aes(t0_2023, t0_2026)) +
geom_point(size = 0.4, alpha = 0.2, colour = "grey30") +
geom_abline(slope = 1, intercept = 0, colour = "#E41A1C", linewidth = 0.6) +
geom_text(data = t0_cor, aes(-Inf, Inf, label = sprintf("r = %.3f", pearson)),
hjust = -0.2, vjust = 1.6, size = 3.4, inherit.aes = FALSE) +
facet_wrap(~ sp, nrow = 1) +
labs(x = "2023 t0 mean log2CPM", y = "2026 t0 mean log2CPM",
title = "Prep-batch exp. scatter") +
theme_prism(base_size = 11) +
theme(strip.text = element_text(face = "italic"))
save_fig(p_t0_scatter, "01_T0_BATCH_MAGNITUDE_SCATTER.png", width = 12, height = 4.5)
p_t0_scatter
Notes: > All sps show high correlation in normalized counts between 2023 and 2026. of greater than .89 r
01. Joint PCA: where does the 2026 t0 sit
Correlation is dominated by unchanged genes, so it can read clean
while a systematic offset hides underneath. The PCA asks whether the two
t0s occupy the same position in the joint space. gap_frac
expresses the distance between the two t0 centroids as a fraction of the
2023 timecourse’s own spread.
print(t0_pca_dist)## # A tibble: 3 × 5
## sp code t0_gap timecourse_spread gap_frac
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 C. glabrata yH001 37.2 130. 0.286
## 2 K. lactis yH149 36.9 136. 0.271
## 3 C. albicans yH714 52.3 114. 0.459p_pca <- map(names(pca_dat), function(cd) {
x <- pca_dat[[cd]]
ggplot(x$df, aes(PC1, PC2, fill = tp_f, shape = batch)) +
geom_point(size = 3.4, stroke = 0.7, colour = "grey20") +
geom_point(data = filter(x$df, cond == "t0"),
shape = 21, size = 5.4, stroke = 1.1, fill = NA, colour = "black") +
scale_shape_manual(values = c("2023" = 21, "2026" = 24), name = NULL) +
scale_fill_viridis_d(option = "turbo", name = "timepoint (min)") +
guides(fill = guide_legend(override.aes = list(shape = 21, colour = "grey20"), ncol = 2)) +
labs(x = sprintf("PC1 (%.1f%%)", x$ve[1]), y = sprintf("PC2 (%.1f%%)", x$ve[2]),
title = sprintf("%s (%s): joint PCA", x$df$sp[1], cd)) +
theme_prism(base_size = 11) +
theme(legend.position = "right",
legend.key.size = unit(0.4, "cm"),
legend.text = element_text(size = 8),
plot.margin = margin(5, 15, 5, 5))
}) |> set_names(names(pca_dat))
walk2(p_pca, names(p_pca), function(p, code) save_fig(p, sprintf("01_JOINT_PCA_%s.png", code), width = 12, height = 7))
walk(p_pca, print)
Notes: Now, what we can see is that there does exist a batch effect as shown in the principle component space. In the ideal case the 2023 and 2026 t0s would overlap, since both are unperturbed cells harvested before any treatment, so any separation between them is technical rather than biological. They do not overlap (t0 gap = 37, 37, 52 PC units for Cgla, Klac, Calb), which gives a direct estimate of the prep-batch offset.
One technical explanation could be the sequencing platform of AVITI vs. Illumina previously. The tricky thing is that I have no ideas on how to bioinformatically probe platform artifacts from expression data. Even if I did, platform is only one arm of the “batch effect”: the 2026 t0 differs from 2023 in extraction batch, library prep, and platform all at once, so any signature I found could belong to any of the three. The most I can say is that the offset exists and is technical, since the two t0s are the same biological state by construction.
The offset runs along roughly the same principal component as the starvation response, so a batch shift is indistinguishable from a shift in time-in-starvation. As such, comparing 2026 mock against 2023 noPi means any handling estimate carries the batch offset with it.
The one anchor available is the shared t0. Both batches contain the same biological state, so the t0-to-t0 difference is a direct estimate of the technical offset and could in principle be used to correct the rest. That correction rests on two assumptions I cannot test with this design: that the offset is constant across conditions (estimated only at t0), and that removing it does not also remove real signal (it lies along the response axis). Worth attempting as a sensitivity analysis, uncorrected vs corrected. If the handling fraction is stable across both, the result stands; if it moves, the design cannot support the claim.
01. Sample-to-sample distance
Do the 2026 t0 replicates pair with each other, does mock separate from t0, and where do the 2026 samples sit relative to the 2023 timecourse?
hm_dist <- map(dat, function(d) {
ord <- as.character(d$meta$sample)
m <- d$lcpm[, ord, drop = FALSE]
dm <- as.matrix(dist(t(m)))
lab <- as.character(d$meta$label)
dimnames(dm) <- list(lab, lab)
ann <- data.frame(batch = d$meta$batch, cond = d$meta$cond, row.names = lab)
pheatmap(dm,
annotation_row = ann, annotation_col = ann,
annotation_colors = list(batch = batch_cols, cond = cond_cols),
cluster_rows = TRUE, cluster_cols = TRUE,
color = colorRampPalette(rev(brewer.pal(9, "Blues")))(255),
fontsize_row = 6.5, fontsize_col = 6.5, silent = TRUE,
main = sprintf("%s (%s): sample-to-sample distance", d$sp, d$code))
})
walk2(hm_dist, names(hm_dist), function(ph, code) {
png(file.path(fig_dir, sprintf("01_SAMPLE_DISTANCE_HEATMAP_%s.png", code)),
width = 11, height = 11, units = "in", res = 300)
grid::grid.newpage(); grid::grid.draw(ph$gtable); dev.off()
})
walk(hm_dist, function(ph) { grid::grid.newpage(); grid::grid.draw(ph$gtable) })
Notes: Sample-to-sample coherence is largely as expected. replicates pair as mutual nearest neighbors, the 2023 series orders by time (early w/ early, late TPs w/ late), and the mock samples cluster together, subclustering with the early 2023 timepoints.
Two replicate-level exceptions: > Cgla mock.60m.rep2 does not pair with rep1 and sits outside the 2026 clade
Calb mock.0m.rep2 does not pair with rep1; rep1 clusters with the handled 30/60 mocks, which is the suspicious pattern for a t0. The Calb within-batch t0 is therefore internally discordant.
Klac is clean throughout, and its 2026 t0 pair clusters adjacent to the 2023 t0 pair, the expected cross-batch topology.
01. Does mock vs t0 recover a handling response
Within-batch contrast. If the new pairing is sound this should produce a response, not a flat result.
print(mock_summary)## # A tibble: 6 × 6
## code sp contrast n_tested up down
## <chr> <chr> <chr> <int> <int> <int>
## 1 yH001 C. glabrata mock30 5190 179 202
## 2 yH001 C. glabrata mock60 5190 173 64
## 3 yH149 K. lactis mock30 5069 266 328
## 4 yH149 K. lactis mock60 5069 256 206
## 5 yH714 C. albicans mock30 5949 469 520
## 6 yH714 C. albicans mock60 5949 429 369p_mock_deg <- mock_summary |>
pivot_longer(c(up, down), names_to = "dir", values_to = "n") |>
mutate(dir = factor(dir, levels = c("down", "up")),
signed = if_else(dir == "down", -n, n)) |>
ggplot(aes(contrast, signed, fill = dir)) +
geom_col(width = 0.6) +
geom_hline(yintercept = 0, colour = "grey30", linewidth = 0.4) +
geom_text(aes(label = n, vjust = if_else(signed >= 0, -0.4, 1.3)), size = 3.2) +
scale_fill_manual(values = dir_cols, name = NULL) +
facet_wrap(~ sp, nrow = 1) +
labs(x = NULL, y = "DEG count (padj < 0.05, |LFC| > 1)",
title = "Handling response: mock vs within-batch t0") +
theme_prism(base_size = 11) +
theme(strip.text = element_text(face = "italic"), legend.position = "bottom")
save_fig(p_mock_deg, "01_MOCK_VS_T0_DEG_COUNTS.png", width = 10, height = 4.5)
p_mock_deg
Notes: Using our standard DEG thresholds we have been using for the 4 sps E009, we recover differential expression of genes and therefore response. Not surprising, but I had to check. These DEG counts are ~ = to the first analysis of the HE26069 set we did earlier in Spring.
01. Library size
libsize_2026 <- libsize |>
filter(batch == "2026") |>
mutate(short = sub("^yH\\d+\\.[A-Za-z]+\\.", "", label))
p_libsize <- ggplot(libsize_2026, aes(short, total / 1e6, fill = cond)) +
geom_col() +
scale_fill_manual(values = cond_cols, name = NULL) +
facet_wrap(~ sp, scales = "free_x", nrow = 1) +
labs(x = NULL, y = "assigned reads (millions)", title = "Library size, 2026 mock arm") +
theme_prism(base_size = 10) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8),
strip.text = element_text(face = "italic"),
legend.position = "top")
save_fig(p_libsize, "01_LIBRARY_SIZE_2026_BARPLOT.png", width = 11, height = 5)
p_libsize
Notes: Our low read counts of ~5x10^6 reads are fixed from the earlier HE26069 run. All reads are above 10 million!
01. Count distribution after TMM
p_count_dist <- ggplot(lcpm_long, aes(label, lcpm, fill = batch)) +
geom_boxplot(outlier.size = 0.15, outlier.alpha = 0.1, linewidth = 0.3) +
scale_fill_manual(values = batch_cols, name = NULL) +
facet_wrap(~ sp, scales = "free_x", nrow = 3) +
labs(x = NULL, y = "log2CPM", title = "Per-library count distribution after TMM") +
theme_prism(base_size = 10) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 5.5),
strip.text = element_text(face = "italic"), legend.position = "bottom")
save_fig(p_count_dist, "01_COUNT_DISTRIBUTION_BOXPLOT.png", width = 14, height = 9)
p_count_dist
01. Per-library RLE
p_rle <- ggplot(rle_long, aes(label, rle, fill = batch)) +
geom_hline(yintercept = 0, colour = "grey60", linewidth = 0.3) +
geom_boxplot(outlier.size = 0.15, outlier.alpha = 0.1, linewidth = 0.3) +
scale_fill_manual(values = batch_cols, name = NULL) +
coord_cartesian(ylim = c(-0.6, 0.6)) +
facet_wrap(~ sp, scales = "free_x", nrow = 3) +
labs(x = NULL, y = "RLE (log2CPM vs within-group median)", title = "Per-library RLE") +
theme_prism(base_size = 10) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 5.5),
strip.text = element_text(face = "italic"), legend.position = "bottom")
save_fig(p_rle, "01_RLE_BOXPLOT.png", width = 14, height = 9)
p_rle
01. Replicate spread
p_rep_sd <- ggplot(rep_sd, aes(group, sd, fill = batch)) +
geom_boxplot(outlier.size = 0.15, outlier.alpha = 0.1, linewidth = 0.3) +
scale_fill_manual(values = batch_cols, name = NULL) +
coord_cartesian(ylim = c(0, quantile(rep_sd$sd, 0.99, na.rm = TRUE))) +
facet_wrap(~ sp, scales = "free_x", nrow = 3) +
labs(x = NULL, y = "within-replicate SD (log2CPM)", title = "Replicate spread per group") +
theme_prism(base_size = 10) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 6),
strip.text = element_text(face = "italic"), legend.position = "bottom")
save_fig(p_rep_sd, "01_REPLICATE_SD_BOXPLOT.png", width = 12, height = 9)
p_rep_sd
alpha <- 0.05
lfc_gate <- 1.0
ycap <- 30
handling_3sps <- mock_lfc |>
mutate(tp = as.integer(sub("mock", "", contrast)),
handling_gene = coalesce(padj < alpha, FALSE) & abs(lfc) > lfc_gate,
dir = case_when(handling_gene & lfc > 0 ~ "induced",
handling_gene & lfc < 0 ~ "repressed",
TRUE ~ "n.s."))
write_csv(handling_3sps |> filter(handling_gene) |>
select(code, sp, tp, gene_id, lfc, padj, dir),
file.path(out_dir, "02_handling_gene_sets_3sps.csv"))
handling_counts <- handling_3sps |> filter(handling_gene) |> count(sp, tp, dir) |>
pivot_wider(names_from = dir, values_from = n, values_fill = 0)
write_csv(handling_counts, file.path(out_dir, "02_handling_deg_counts_3sps.csv"))
vol_02 <- handling_3sps |>
filter(!is.na(padj)) |>
mutate(neglp = -log10(padj), capped = neglp > ycap, y = pmin(neglp, ycap),
dir = factor(dir, levels = c("n.s.", "repressed", "induced")))
print(handling_counts)## # A tibble: 6 × 4
## sp tp induced repressed
## <chr> <int> <int> <int>
## 1 C. albicans 30 469 520
## 2 C. albicans 60 429 369
## 3 C. glabrata 30 179 202
## 4 C. glabrata 60 173 64
## 5 K. lactis 30 266 328
## 6 K. lactis 60 256 206p_hvol <- ggplot(vol_02, aes(lfc, y)) +
geom_vline(xintercept = c(-lfc_gate, lfc_gate), linetype = 2, colour = "grey65", linewidth = 0.3) +
geom_hline(yintercept = -log10(alpha), linetype = 2, colour = "grey65", linewidth = 0.3) +
geom_point(data = filter(vol_02, dir == "n.s."), colour = "grey82", size = 0.35, alpha = 0.3) +
geom_point(data = filter(vol_02, dir != "n.s."), aes(colour = dir, shape = capped), size = 0.85, alpha = 0.8) +
scale_colour_manual(values = c(induced = "#08519C", repressed = "#6BAED6"), name = NULL,
breaks = c("induced", "repressed")) +
scale_shape_manual(values = c(`FALSE` = 16, `TRUE` = 17), guide = "none") +
facet_grid(tp ~ sp) +
labs(x = "log2FC mock vs t0", y = expression(-log[10]~adjusted~italic(p)),
title = "Handling effect per species (within-2026 contrast)",
caption = sprintf("handling gene = padj < %.2f and |log2FC| > %.0f; y capped at %d", alpha, lfc_gate, ycap)) +
theme_prism(base_size = 11) +
theme(strip.text.x = element_text(face = "italic"), legend.position = "bottom")
save_fig(p_hvol, "02_HANDLING_VOLCANO_3SPS.png", width = 12, height = 8)
p_hvol
Notes: Interestingly, C.g seems to have the weakest
response to handling. C.a and K.l have the same bidirectional spread of
DEGs as we saw in S.cer.
nopi_lfc <- map_dfr(dat, function(d) {
cd <- d$meta |> filter(batch == "2023", tp %in% c(0, 30, 60)) |> as.data.frame()
rownames(cd) <- cd$sample
cd$group <- factor(paste0("t", cd$tp), levels = c("t0", "t30", "t60"))
mm <- d$counts[, as.character(cd$sample), drop = FALSE]
mm <- mm[rowSums(mm >= 10) >= 2, ]
dds <- DESeq(DESeqDataSetFromMatrix(round(mm), cd, ~ group), quiet = TRUE)
map_dfr(c(30, 60), function(t) {
r <- results(dds, contrast = c("group", paste0("t", t), "t0"),
independentFiltering = FALSE, cooksCutoff = FALSE)
tibble(code = d$code, sp = d$sp, tp = t,
gene_id = rownames(r),
lfc_nopi = r$log2FoldChange,
padj_nopi = r$padj)
})
})
write_csv(nopi_lfc, file.path(out_dir, "03_nopi_vs_t0_2023_lfc.csv"))
overlap_tbl <- nopi_lfc |>
inner_join(handling_3sps |> select(code, tp, gene_id, lfc_mock = lfc,
padj_mock = padj, handling_gene),
by = c("code", "tp", "gene_id")) |>
mutate(nopi_deg = coalesce(padj_nopi < alpha, FALSE) & abs(lfc_nopi) > lfc_gate,
same_dir = sign(lfc_nopi) == sign(lfc_mock),
deg_class = case_when(
nopi_deg & handling_gene & same_dir ~ "handling-shared",
nopi_deg ~ "noPi-specific",
TRUE ~ "other"))
write_csv(overlap_tbl |> filter(deg_class != "other") |>
select(code, sp, tp, gene_id, lfc_nopi, padj_nopi, lfc_mock, padj_mock, deg_class),
file.path(out_dir, "03_deg_overlap_classes_3sps.csv"))
decomp_summary <- overlap_tbl |>
filter(nopi_deg) |>
summarise(.by = c(sp, tp),
nopi_degs = n(),
handling_shared = sum(deg_class == "handling-shared"),
nopi_specific = sum(deg_class == "noPi-specific"),
frac_handling = round(handling_shared / nopi_degs, 3))
write_csv(decomp_summary, file.path(out_dir, "03_decomposition_summary_3sps.csv"))
venn_03 <- overlap_tbl |>
mutate(mock_deg = handling_gene) |>
summarise(.by = c(sp, tp),
nopi_only = sum(nopi_deg & !(mock_deg & same_dir)),
shared = sum(nopi_deg & mock_deg & same_dir),
handling_only = sum(mock_deg & !(nopi_deg & same_dir)),
nopi_total = sum(nopi_deg),
handling_total = sum(mock_deg)) |>
mutate(panel_lab = sprintf("%s | %d min\nnoPi: %d handling: %d", sp, tp, nopi_total, handling_total))
write_csv(venn_03, file.path(out_dir, "03_venn_counts_3sps.csv"))
bar_03 <- overlap_tbl |>
filter(nopi_deg) |>
mutate(direction = if_else(lfc_nopi > 0, "induced", "repressed")) |>
count(sp, tp, deg_class, direction) |>
mutate(signed = if_else(direction == "repressed", -n, n),
key = paste(deg_class, direction, sep = "."))
ho_hist <- bind_rows(
overlap_tbl |> filter(handling_gene, !(nopi_deg & same_dir)) |> mutate(cls = "handling_only"),
overlap_tbl |> filter(nopi_deg, handling_gene, same_dir) |> mutate(cls = "shared")
)
print(decomp_summary)## # A tibble: 6 × 6
## sp tp nopi_degs handling_shared nopi_specific frac_handling
## <chr> <dbl> <int> <int> <int> <dbl>
## 1 C. glabrata 30 640 47 593 0.073
## 2 C. glabrata 60 2007 130 1877 0.065
## 3 K. lactis 30 1835 520 1315 0.283
## 4 K. lactis 60 792 88 704 0.111
## 5 C. albicans 30 1916 700 1216 0.365
## 6 C. albicans 60 1681 551 1130 0.328fill_03 <- c("handling-shared.induced" = "#2166AC", "handling-shared.repressed" = "#8EC4DE",
"noPi-specific.induced" = "#B2182B", "noPi-specific.repressed" = "#E8A0A8")
p_decomp <- ggplot(bar_03, aes(deg_class, signed, fill = key)) +
geom_col(width = 0.6) +
geom_hline(yintercept = 0, colour = "grey30", linewidth = 0.4) +
geom_text(aes(label = n, vjust = if_else(signed >= 0, -0.4, 1.3)), size = 3) +
scale_fill_manual(values = fill_03, guide = "none") +
facet_grid(tp ~ sp) +
labs(x = NULL, y = "noPi DEG count (induced \u2191 / repressed \u2193)",
title = "noPi DEG partition: handling-shared vs noPi-specific") +
theme_prism(base_size = 11) +
theme(strip.text.x = element_text(face = "italic"),
axis.text.x = element_text(angle = 15, hjust = 1))
save_fig(p_decomp, "03_DEG_OVERLAP_PARTITION_3SPS.png", width = 12, height = 8)
p_decomp
venn_03 <- overlap_tbl |>
mutate(mock_deg = handling_gene) |>
summarise(.by = c(sp, tp),
nopi_only = sum(nopi_deg & !(mock_deg & same_dir)),
shared = sum(nopi_deg & mock_deg & same_dir),
handling_only = sum(mock_deg & !(nopi_deg & same_dir)),
nopi_total = sum(nopi_deg),
handling_total = sum(mock_deg)) |>
mutate(panel_lab = sprintf("%s | %d min\nnoPi: %d handling: %d", sp, tp, nopi_total, handling_total))
write_csv(venn_03, file.path(out_dir, "03_venn_counts_3sps.csv"))
circ <- function(cx, r = 1.15, n = 200) {
th <- seq(0, 2 * pi, length.out = n)
tibble(x = cx + r * cos(th), y = r * sin(th))
}
c1 <- circ(-0.62); c2 <- circ(0.62)
venn_long <- venn_03 |>
pivot_longer(c(nopi_only, shared, handling_only), names_to = "region", values_to = "n") |>
mutate(x = case_when(region == "nopi_only" ~ -1.15,
region == "shared" ~ 0,
region == "handling_only" ~ 1.15),
y = 0)
p_venn <- ggplot(venn_long) +
geom_path(data = c1, aes(x, y), colour = "#B2182B", linewidth = 0.9) +
geom_path(data = c2, aes(x, y), colour = "#2166AC", linewidth = 0.9) +
geom_text(aes(x, y, label = n), size = 4.2, fontface = "bold") +
annotate("text", x = -1.15, y = 1.45, label = "noPi", colour = "#B2182B", size = 3.4, fontface = "bold") +
annotate("text", x = 1.15, y = 1.45, label = "handling", colour = "#2166AC", size = 3.4, fontface = "bold") +
facet_wrap(~ panel_lab, nrow = 3, dir = "v") +
coord_fixed(xlim = c(-2.4, 2.4), ylim = c(-1.7, 1.9)) +
labs(title = "noPi vs handling DEG overlap (same-direction)",
caption = sprintf("DEG = padj < %.2f, |log2FC| > %.0f, each arm within its own batch; shared = significant in both, same direction", alpha, lfc_gate)) +
theme_prism(base_size = 10) +
theme(axis.line = element_blank(), axis.text = element_blank(),
axis.ticks = element_blank(), axis.title = element_blank(),
strip.text = element_text(size = 8.5))
save_fig(p_venn, "03_DEG_OVERLAP_VENN_3SPS.png", width = 9, height = 10)
p_venn