Skip to content

Simplify dream contrast handling #11584

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
grst wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Simplify dream contrast handling #11584

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
2486e25
Simplify dream contrast handling
grst May 11, 2026
d995b36
fix pre-commit
grst May 11, 2026
f30d573
Use lme4 as source for nobars
grst May 11, 2026
3e79fb3
Apply suggestion from @delfiterradas
grst May 11, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
84 changes: 36 additions & 48 deletions modules/nf-core/variancepartition/dream/templates/dream.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -272,8 +272,39 @@ if (as.logical(opt\$apply_voom)) {
vobjDream <- countMatrix
}

# Fit the DREAM model with ddf and reml options
fitmm <- dream(vobjDream, form, metadata, ddf = opt\$ddf, reml = opt\$reml, BPPARAM = bp)
# If contrast_string is provided, use that for makeContrast
if (!is.null(opt\$contrast_string)) {
contrast_string <- opt\$contrast_string
} else {
# Construct the contrast_string
treatment_target <- make.names(paste0(opt\$contrast_variable, opt\$contrast_target))
treatment_reference <- make.names(paste0(opt\$contrast_variable, opt\$contrast_reference))
contrast_string <- paste0(treatment_target, "-", treatment_reference)
}

# define contrasts outside DREAM
design = model.matrix(lme4::nobars(form), metadata)

# print diagnostic output of design matrix
head(design, 3)
cat("Raw coefficient names:\n")
print(colnames(design))

# Use "normalized" coefficient names to compute contrast vector
normalized_coef_names = make.names(colnames(design))
cat("Coefficient names used for contrasts:", paste(normalized_coef_names, collapse = ", "), "\n")

contrast_mat = makeContrasts(contrast = contrast_string, levels = normalized_coef_names)
# however, reset the coefficient names in the contrast matrix to use the "original" contrast names used by DREAM
# DREAM will check if the coefficient names match its internal design matrix, so there is no danger that
# design() and dream() will do things differently without noticing.
rownames(contrast_mat) = colnames(design)
stopifnot("exactly one contrast defined"=ncol(contrast_mat) == 1)

# Fit the DREAM model with ddf and reml options.
# Specifying contrast matrix directly here saves time refitting the model and avoids downstream issues
# with repeated calling of eBayes
fitmm <- dream(vobjDream, form, metadata, ddf = opt\$ddf, reml = opt\$reml, BPPARAM = bp, L=contrast_mat)

# Parse stdev_coef_lim and winsor_tail_p into numeric vectors
stdev_coef_lim_vals <- as.numeric(unlist(strsplit(opt\$stdev_coef_lim, ",")))
Expand All @@ -285,52 +316,9 @@ fitmm <- eBayes(fitmm, proportion = opt\$proportion,
trend = opt\$trend, robust = opt\$robust,
winsor.tail.p = winsor_tail_p_vals)

# Display design matrix
head(fitmm\$design, 3)
cat("Raw coefficient names from dream():\n")
print(colnames(fitmm\$design))

# Normalize coefficient names consistently before building contrasts
normalized_coef_names <- make.names(colnames(fitmm\$design))
colnames(fitmm\$design) <- normalized_coef_names

if (!is.null(colnames(fitmm\$coefficients))) {
colnames(fitmm\$coefficients) <- normalized_coef_names
}
if (!is.null(colnames(fitmm\$stdev.unscaled))) {
colnames(fitmm\$stdev.unscaled) <- normalized_coef_names
}
if (!is.null(fitmm\$cov.coefficients)) {
rownames(fitmm\$cov.coefficients) <- normalized_coef_names
colnames(fitmm\$cov.coefficients) <- normalized_coef_names
}

cat("Coefficient names used for contrasts:", paste(normalized_coef_names, collapse = ", "), "\n")

# If contrast_string is provided, use that for makeContrast
if (!is.null(opt\$contrast_string)) {
contrast_string <- opt\$contrast_string
} else {
# Construct the contrast_string
treatment_target <- make.names(paste0(opt\$contrast_variable, opt\$contrast_target))
treatment_reference <- make.names(paste0(opt\$contrast_variable, opt\$contrast_reference))
contrast_string <- paste0(treatment_target, "-", treatment_reference)
}

# Use makeContrasts if contrast_string exists
if (is_valid_string(contrast_string)) {
Copy link
Copy Markdown
Member Author

@grst grst May 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@delfiterradas would I need this check, too? I.e. is there a valid usecase to call this module without contrast, e.g. for normalization only?

Is there any way this could have been false, als both branches of the previous if set a contrast_string?

Copy link
Copy Markdown
Contributor

@delfiterradas delfiterradas May 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are right, that check was redundant!

cat("Using contrast string:", contrast_string, "\n")

contrast_matrix <- makeContrasts(contrast = contrast_string, levels = normalized_coef_names)
fit2 <- contrasts.fit(fitmm, contrast_matrix)
fit2 <- eBayes(fit2, proportion = opt\$proportion,
stdev.coef.lim = stdev_coef_lim_vals,
trend = opt\$trend, robust = opt\$robust,
winsor.tail.p = winsor_tail_p_vals)
results <- topTable(fit2, number = Inf,
adjust.method = opt\$adjust.method,
p.value = opt\$p.value, lfc = opt\$lfc, confint = opt\$confint)
}
results <- topTable(fitmm, coef=colnames(contrast_mat), number = Inf,
adjust.method = opt\$adjust.method,
p.value = opt\$p.value, lfc = opt\$lfc, confint = opt\$confint)

results\$gene_id <- rownames(results)
results <- results[, c("gene_id", setdiff(names(results), "gene_id"))]
Expand Down
Loading