---
title: "Differential Expression Analysis (Tabbed Report)"
format:
html:
toc-location: left
grid:
sidebar-width: 300px
include-after-body: fgcz-plot-finder.html
toc: true
params:
se_file: null
prolfquapp_source_path: null
fdr_threshold: 0.05
diff_threshold: 1
vignette: >
%\VignetteIndexEntry{Differential Expression Analysis (Tabbed Report)}
%\VignetteEngine{quarto::format}
%\VignetteEncoding{UTF-8}
---
```{r setup}
#| include: false
source_path <- params$prolfquapp_source_path
if (!is.null(source_path) && length(source_path) == 1 && !is.na(source_path) && nzchar(source_path)) {
r_dir <- file.path(source_path, "R")
description <- file.path(source_path, "DESCRIPTION")
package_name <- if (file.exists(description)) {
unname(read.dcf(description)[1, "Package"])
} else {
NA_character_
}
is_source <- file.exists(file.path(source_path, "DESCRIPTION")) &&
identical(package_name, "prolfquapp") &&
dir.exists(r_dir) &&
length(list.files(r_dir, pattern = "[.]R$")) > 0
if (!is_source) {
stop("`prolfquapp_source_path` does not point to a prolfquapp source tree.", call. = FALSE)
}
if (!requireNamespace("devtools", quietly = TRUE)) {
stop("Package 'devtools' is required when `prolfquapp_source_path` is set.", call. = FALSE)
}
suppressPackageStartupMessages(devtools::load_all(source_path, quiet = TRUE))
}
required_packages <- c(
"DT",
"dplyr",
"ggplot2",
"grid",
"gridExtra",
"plotly",
"prolfqua",
"prolfquapp",
"tibble",
"UpSetR"
)
missing_packages <- required_packages[
!vapply(required_packages, requireNamespace, quietly = TRUE, FUN.VALUE = logical(1))
]
if (length(missing_packages) > 0) {
stop(
"Missing package(s) required to render this report: ",
paste(missing_packages, collapse = ", "),
call. = FALSE
)
}
knitr::opts_chunk$set(
out.width = "600px",
fig.align = "center"
)
se_file <- params$se_file
if (is.null(se_file) || !nzchar(se_file)) {
se_file <- system.file("extdata", "3106962.rds", package = "prolfquapp")
}
report <- prolfquapp:::se_report_lfqdata(se_file)
se <- report$se
meta <- report$metadata
lfq_raw <- report$lfq_raw
lfq_transformed <- report$lfq_transformed
feature_annotation <- report$feature_annotation
contrast_table <- report$contrast_table
contrast_object <- report$contrast_object
is_saint_report <- identical(report$contrast_model, "saint") ||
all(c("Bait", "BFDR", "log2_EFCs") %in% colnames(contrast_table))
dir.create("plots", showWarnings = FALSE, recursive = TRUE)
empty_report_plot <- function(label) {
ggplot2::ggplot() +
ggplot2::annotate("text", x = 0, y = 0, label = label) +
ggplot2::theme_void()
}
safe_plot <- function(expr, label) {
tryCatch(expr, error = function(e) {
empty_report_plot(paste(label, conditionMessage(e), sep = "\n"))
})
}
sized_widget <- function(widget, width = 600) {
htmltools::tags$div(
style = sprintf("max-width: %dpx; margin: 0 auto;", width),
widget
)
}
draw_or_empty <- function(expr, label) {
tryCatch(
expr,
error = function(e) {
show_plot(empty_report_plot(paste(label, conditionMessage(e), sep = "\n")))
}
)
}
show_plot <- function(plot) {
if (is.null(plot)) {
plot <- empty_report_plot("No plot is available for this data set.")
}
print(plot)
invisible(plot)
}
save_if_ggplot <- function(plot, name, width = 8, height = 5) {
if (inherits(plot, "ggplot")) {
tryCatch(
{
ggplot2::ggsave(file.path("plots", paste0(name, ".png")), plot, width = width, height = height, dpi = 300)
ggplot2::ggsave(file.path("plots", paste0(name, ".pdf")), plot, width = width, height = height)
},
error = function(e) NULL
)
}
invisible(plot)
}
plot_experiment_supported_feature_counts <- function(lfq, annotation, min_children = 2) {
if (is.null(annotation) || nrow(annotation) == 0 || !"nrPeptides" %in% colnames(annotation)) {
return(NULL)
}
hierarchy_col <- lfq$relevant_hierarchy_keys()[[1]]
if (!hierarchy_col %in% colnames(annotation)) {
return(NULL)
}
supported_features <- annotation |>
dplyr::filter(!is.na(.data$nrPeptides), .data$nrPeptides >= min_children)
supported_features <- unique(supported_features[[hierarchy_col]])
if (length(supported_features) == 0) {
return(NULL)
}
filtered_data <- lfq$data_long() |>
dplyr::filter(.data[[hierarchy_col]] %in% supported_features)
if (nrow(filtered_data) == 0) {
return(NULL)
}
prolfqua::LFQData$new(
tibble::as_tibble(filtered_data),
lfq$get_config()
)$get_Summariser()$plot_hierarchy_counts_sample(nr_children = 1)
}
html_escape <- function(x) {
x <- gsub("&", "&", x, fixed = TRUE)
x <- gsub("<", "<", x, fixed = TRUE)
x <- gsub(">", ">", x, fixed = TRUE)
x <- gsub('"', """, x, fixed = TRUE)
x <- gsub("'", "'", x, fixed = TRUE)
x
}
html_link <- function(url) {
sprintf(
'<a href="%s" target="_blank" rel="noopener noreferrer">%s</a>',
html_escape(url),
html_escape(url)
)
}
dt_table <- function(df, caption = NULL, page_length = 10, escape = TRUE) {
if (is.null(df) || nrow(df) == 0) {
df <- data.frame(message = "No rows available.")
}
DT::datatable(
df,
caption = caption,
filter = "bottom",
extensions = "Buttons",
rownames = FALSE,
class = "compact stripe",
escape = escape,
options = list(pageLength = page_length, scrollX = TRUE, dom = "Blfrtip", buttons = c("csv", "excel"))
)
}
contrast_definition_table <- function(contrasts) {
if (is.null(contrasts) || length(contrasts) == 0) {
return(tibble::tibble())
}
if (is.data.frame(contrasts)) {
return(tibble::as_tibble(contrasts))
}
values <- unlist(contrasts, recursive = FALSE, use.names = TRUE)
if (length(values) == 0) {
return(tibble::tibble())
}
tibble::tibble(
contrast = names(values),
definition = as.character(values)
)
}
feature_id_column <- function(df) {
candidates <- c("feature_id", "protein_Id", "peptide_Id", "metabolite_Id", "compound_Id", "ID")
hit <- candidates[candidates %in% colnames(df)]
if (length(hit) == 0) {
return(NULL)
}
hit[[1]]
}
significant_sets_by_direction <- function(significant_contrasts, direction = c("all", "increased", "decreased")) {
direction <- match.arg(direction)
feature_col <- feature_id_column(significant_contrasts)
if (is.null(feature_col) || nrow(significant_contrasts) == 0) {
return(list())
}
contrast_col <- if ("Bait" %in% colnames(significant_contrasts)) "Bait" else "contrast"
effect_col <- if ("log2_EFCs" %in% colnames(significant_contrasts)) "log2_EFCs" else "diff"
if (!all(c(contrast_col, effect_col) %in% colnames(significant_contrasts))) {
return(list())
}
data <- significant_contrasts
if (identical(direction, "increased")) {
data <- dplyr::filter(data, .data[[effect_col]] > 0)
} else if (identical(direction, "decreased")) {
data <- dplyr::filter(data, .data[[effect_col]] < 0)
}
sets <- split(data[[feature_col]], data[[contrast_col]])
sets <- lapply(sets, unique)
sets[lengths(sets) > 0]
}
draw_upset_or_empty <- function(sets, empty_label) {
if (length(sets) < 2) {
return(show_plot(empty_report_plot(empty_label)))
}
draw_or_empty(
UpSetR::upset(UpSetR::fromList(sets), order.by = "freq"),
"Could not draw contrast agreement plot."
)
}
contrast_plotter <- NULL
if (!is.null(contrast_object)) {
contrast_plotter <- contrast_object$get_Plotter(
fc_threshold = params$diff_threshold,
fdr_threshold = params$fdr_threshold
)
} else if (is_saint_report && !is.null(feature_id_column(contrast_table))) {
contrast_plotter <- prolfqua::ContrastsPlotter$new(
contrast_table,
subject_id = feature_id_column(contrast_table),
fcthresh = params$diff_threshold,
volcano = list(list(score = "BFDR", name = "FDR", thresh = params$fdr_threshold)),
histogram = list(
list(score = "BFDR", xlim = c(0, 1, 0.05)),
list(score = "SaintScore", xlim = c(0, 1, 0.05))
),
score = list(list(score = "SaintScore", thresh = 0.75)),
modelName = "modelName",
diff = "log2_EFCs",
contrast = "Bait"
)
}
significant_contrasts <- if (
is_saint_report &&
nrow(contrast_table) > 0 &&
all(c("BFDR", "log2_EFCs") %in% colnames(contrast_table))
) {
contrast_table |>
dplyr::filter(.data$BFDR < params$fdr_threshold, .data$log2_EFCs > params$diff_threshold)
} else if (nrow(contrast_table) > 0 && all(c("FDR", "diff") %in% colnames(contrast_table))) {
contrast_table |>
dplyr::filter(.data$FDR < params$fdr_threshold, abs(.data$diff) > params$diff_threshold)
} else {
tibble::tibble()
}
feature_col <- feature_id_column(contrast_table)
report_provenance <- prolfquapp:::.report_provenance(
project_spec = meta$report_provenance
)
```
::: {.panel-tabset}
# Overview
```{r overview-data-summary}
#| results: asis
cat(prolfquapp:::.report_overview_cards(lfq_raw))
```
<img src="`r knitr::image_uri(system.file("report-assets", "differential-expression-tabset.png", package = "prolfquapp"))`" alt="Visual abstract: quantified protein abundances and sample annotation define an analysis design and contrasts, which produce quality checks and differential-abundance results." style="width: 100%; max-height: 22rem; object-fit: contain;">
This report brings together experimental design, feature detection, quality control, and differential-abundance results. The tabs make the analysis path explicit: define the design and contrasts, inspect data quality, then review the result table.
Input: quantified feature abundances from `r nrow(lfq_raw$factors())` samples. The detailed source-data reference is recorded in **Session Info**.
# Settings
::: {.panel-tabset}
## Analysis
```{r}
#| label: tbl-analysis-settings
settings <- data.frame(
parameter = c("FDR threshold", "Difference threshold", "Model formula"),
value = c(
params$fdr_threshold,
params$diff_threshold,
if (!is.null(meta$formula$formula)) meta$formula$formula[[1]] else NA
)
)
dt_table(settings, caption = "Analysis settings", page_length = 5)
```
## Design
```{r}
#| label: tbl-sample-counts
sample_counts <- as.data.frame(table(lfq_raw$factors()[lfq_raw$relevant_factor_keys()]))
colnames(sample_counts)[ncol(sample_counts)] <- "samples"
dt_table(sample_counts, caption = "Number of samples per condition")
```
```{r}
#| label: tbl-sample-annotation
sample_annotation <- lfq_raw$factors()
dt_table(sample_annotation, caption = "Sample annotation and experimental factors")
```
## Contrasts
```{r}
#| label: tbl-contrast-definitions
contrast_definitions <- contrast_definition_table(meta$contrasts)
dt_table(contrast_definitions, caption = "Contrast definitions", page_length = 10)
```
:::
# Feature Detection
::: {.panel-tabset}
## Counts
```{r}
#| label: fig-protein-counts
#| fig-cap: "Number of quantified features per sample."
#| fig-width: 10
#| fig-height: 6
protein_counts <- safe_plot(
lfq_raw$get_Summariser()$plot_hierarchy_counts_sample(nr_children = 1),
"Could not draw feature counts per sample."
)
save_if_ggplot(protein_counts, "protein_counts_per_sample", width = 10, height = 6)
show_plot(protein_counts)
```
```{r}
#| label: fig-protein-counts-two-peptides
#| fig-cap: "Number of quantified features per sample among features with at least two peptides in the experiment."
#| fig-width: 10
#| fig-height: 6
protein_counts_two_peptides <- safe_plot(
plot_experiment_supported_feature_counts(
lfq_raw,
feature_annotation,
min_children = 2
),
"Could not draw feature counts with at least two peptides in the experiment."
)
if (is.null(protein_counts_two_peptides)) {
protein_counts_two_peptides <- empty_report_plot(
"No features with at least two peptides in the experiment are available for this data set."
)
}
save_if_ggplot(protein_counts_two_peptides, "protein_counts_two_peptides_experiment", width = 10, height = 6)
show_plot(protein_counts_two_peptides)
```
## Overlap
```{r}
#| label: fig-protein-detection-overlap
#| fig-cap: "UpSet plot of feature-detection overlap between experimental groups: each vertical bar counts the features detected (at least one non-missing measurement) in one specific combination of groups, the combination is shown by the filled dots below the bar, and the horizontal set-size bars on the left give each group's total number of detected features."
#| fig-width: 9
#| fig-height: 6
draw_or_empty(
lfq_raw$get_Summariser()$upset_interaction_missing_stats(tr = 1),
"Could not draw feature detection overlap plot."
)
```
:::
# Quality Control
::: {.panel-tabset}
## Missing Values
```{r}
#| label: fig-missing-heatmap
#| fig-cap: "Missing-value heatmap of raw feature abundances: rows are the features with at least one missing value, columns are samples, and each cell marks whether the measurement is missing (black) or observed (white); columns are clustered by binary distance."
#| fig-width: 7
#| fig-height: 6
show_plot(safe_plot(
lfq_raw$get_Plotter()$na_heatmap(),
"Could not draw missing value heatmap."
))
```
```{r}
#| label: fig-missingness-per-group
#| fig-cap: "Missingness of raw feature abundances per experimental group (three panels, faceted by group). Left: number of features versus number of missing values in the group (bars); middle: the same counts as a cumulative sum over increasing number of missing values; right: kernel density of mean raw abundance (log10 x-axis) coloured by the number of missing values."
#| fig-width: 15
#| fig-height: 4.5
missing_per_group <- safe_plot(
lfq_raw$get_Summariser()$plot_missingness_per_group(),
"Could not draw missingness per group."
)
missing_cumsum <- safe_plot(
lfq_raw$get_Summariser()$plot_missingness_per_group_cumsum(),
"Could not draw cumulative missingness per group."
)
missing_hist <- safe_plot(
lfq_raw$get_Plotter()$missigness_histogram(),
"Could not draw missingness intensity histogram."
)
combined_missing <- gridExtra::arrangeGrob(missing_per_group, missing_cumsum, missing_hist, ncol = 3)
save_if_ggplot(missing_per_group, "missingness_per_group", width = 9, height = 4)
save_if_ggplot(missing_cumsum, "missingness_per_group_cumsum", width = 9, height = 4)
save_if_ggplot(missing_hist, "missingness_histogram", width = 9, height = 4)
grid::grid.draw(combined_missing)
```
## Abundance Distributions
```{r}
#| label: fig-abundance-density
#| fig-cap: "Raw and transformed abundance distributions. Each curve represents one sample; the x-axes show raw abundance and transformed intensity, and the y-axes show kernel density."
#| fig-width: 12
#| fig-height: 5
#| out-width: "600px"
raw_density <- safe_plot(
lfq_raw$get_Plotter()$intensity_distribution_density(),
"Could not draw raw abundance density."
) +
ggplot2::labs(tag = "A")
trans_density <- safe_plot(
lfq_transformed$get_Plotter()$intensity_distribution_density(),
"Could not draw transformed abundance density."
) +
ggplot2::labs(tag = "B")
save_if_ggplot(raw_density, "raw_abundance_density")
save_if_ggplot(trans_density, "transformed_abundance_density")
sized_widget(
prolfquapp::plotly_ggplot_subplot(
raw_density,
trans_density,
width = 600,
height = 520
)
)
```
## Variance
```{r}
#| label: fig-variance-violin
#| fig-cap: "Raw feature abundance CV, log2 raw feature abundance SD, and transformed feature abundance SD distributions."
#| fig-width: 15
#| fig-height: 5
raw_violin <- safe_plot(
lfq_raw$get_Stats()$violin(),
"Could not draw raw CV violin plot."
) +
ggplot2::labs(y = "CV")
log2_raw <- lfq_raw$get_Transformer()$log2()$lfq
log2_raw_violin <- safe_plot(
log2_raw$get_Stats()$violin(),
"Could not draw log2 raw SD violin plot."
) +
ggplot2::labs(y = "SD (log2 raw)")
trans_violin <- safe_plot(
lfq_transformed$get_Stats()$violin(),
"Could not draw transformed SD violin plot."
) +
ggplot2::labs(y = "SD (transformed)")
save_if_ggplot(raw_violin, "raw_cv_violin")
save_if_ggplot(log2_raw_violin, "log2_raw_sd_violin")
save_if_ggplot(trans_violin, "transformed_sd_violin")
gridExtra::grid.arrange(raw_violin, log2_raw_violin, trans_violin, nrow = 1)
```
```{r}
#| label: tbl-variance-summary
variance_summary <- dplyr::bind_rows(
raw_cv = lfq_raw$get_Stats()$stats_quantiles()$wide,
log2_sd = log2_raw$get_Stats()$stats_quantiles()$wide,
transformed_sd = lfq_transformed$get_Stats()$stats_quantiles()$wide,
.id = "source"
) |>
dplyr::mutate(dplyr::across(where(is.numeric), \(x) round(x, 3)))
dt_table(
variance_summary,
caption = paste(
"Quantiles of per-feature dispersion across samples for three data layers:",
"raw-abundance CV, log2 raw-abundance SD, and transformed-abundance SD",
"(one row per source)."
),
page_length = 15
)
```
## Sample Structure
::: {.panel-tabset}
### PCA
```{r}
#| label: fig-pca
#| fig-cap: "PCA of transformed feature abundances; each point is a sample plotted on PC1 versus PC2 (percent variance shown on the axes) and coloured by experimental group (features with missing values are removed)."
#| fig-width: 7
#| fig-height: 6
pca_plot <- safe_plot(
lfq_transformed$get_Plotter()$pca(),
"Could not draw PCA."
)
save_if_ggplot(pca_plot, "pca", width = 7, height = 6)
show_plot(pca_plot)
```
### Correlation
```{r}
#| label: fig-correlation
#| fig-cap: "Sample correlation heatmap for transformed feature abundances."
#| fig-width: 7
#| fig-height: 6
show_plot(safe_plot(
lfq_transformed$get_Plotter()$heatmap_cor(),
"Could not draw sample correlation heatmap."
))
```
### Abundance Heatmap
```{r}
#| label: fig-abundance-heatmap
#| fig-cap: "Heatmap of transformed feature abundances. Cell color encodes the row z-score on a green-black-red scale (green: low, black: average, red: high); missing values are shown in light gray."
#| fig-width: 7
#| fig-height: 6
show_plot(safe_plot(
lfq_transformed$get_Plotter()$heatmap(),
"Could not draw abundance heatmap."
))
```
:::
:::
# Differential Abundance
::: {.panel-tabset}
## Summary
```{r}
#| label: tbl-significant-summary
if (nrow(significant_contrasts) > 0) {
contrast_col <- if ("Bait" %in% colnames(significant_contrasts)) "Bait" else "contrast"
effect_col <- if ("log2_EFCs" %in% colnames(significant_contrasts)) "log2_EFCs" else "diff"
signif_summary <- significant_contrasts |>
dplyr::mutate(direction = ifelse(.data[[effect_col]] > 0, "increased", "decreased")) |>
dplyr::count(.data[[contrast_col]], .data$direction, name = "features")
} else {
signif_summary <- tibble::tibble()
}
dt_table(signif_summary, caption = "Significant features by contrast and direction")
```
```{r}
#| label: tbl-differential-abundance-summary
if (is_saint_report && nrow(contrast_table) > 0 && all(c("Bait", "BFDR", "log2_EFCs") %in% colnames(contrast_table))) {
abundance_summary <- contrast_table |>
dplyr::mutate(
significant = .data$BFDR < params$fdr_threshold & .data$log2_EFCs > params$diff_threshold
) |>
dplyr::group_by(.data$Bait) |>
dplyr::summarise(
tested = dplyr::n(),
significant = sum(.data$significant, na.rm = TRUE),
not_significant = tested - significant,
.groups = "drop"
)
} else if (nrow(contrast_table) > 0 && all(c("contrast", "FDR", "diff") %in% colnames(contrast_table))) {
abundance_summary <- contrast_table |>
dplyr::mutate(
significant = .data$FDR < params$fdr_threshold & abs(.data$diff) > params$diff_threshold
) |>
dplyr::group_by(.data$contrast) |>
dplyr::summarise(
tested = dplyr::n(),
significant = sum(.data$significant, na.rm = TRUE),
not_significant = tested - significant,
.groups = "drop"
)
} else {
abundance_summary <- tibble::tibble()
}
dt_table(abundance_summary, caption = "Differential abundance calls by contrast")
```
## Contrast Agreement
::: {.panel-tabset}
### All
```{r}
#| label: fig-significant-feature-overlap-all
#| fig-cap: "UpSet plot of significant features (both FDR and fold-change thresholds met, either direction) shared between contrasts: each vertical bar counts the features called significant in one specific combination of contrasts, shown by the filled dots below the bar, and the horizontal set-size bars on the left give each contrast's total number of significant features."
#| fig-width: 9
#| fig-height: 6
draw_upset_or_empty(
significant_sets_by_direction(significant_contrasts, "all"),
"Significant features are available for fewer than two contrasts."
)
```
### Increased
```{r}
#| label: fig-significant-feature-overlap-increased
#| fig-cap: "UpSet plot of significantly increased-abundance features (positive fold-change passing both thresholds) shared between contrasts: each vertical bar counts the up-regulated features in one specific combination of contrasts, shown by the filled dots below the bar, and the horizontal set-size bars on the left give each contrast's total number of up-regulated features."
#| fig-width: 9
#| fig-height: 6
draw_upset_or_empty(
significant_sets_by_direction(significant_contrasts, "increased"),
"Increased-abundance significant features are available for fewer than two contrasts."
)
```
### Decreased
```{r}
#| label: fig-significant-feature-overlap-decreased
#| fig-cap: "UpSet plot of significantly decreased-abundance features (negative fold-change passing both thresholds) shared between contrasts: each vertical bar counts the down-regulated features in one specific combination of contrasts, shown by the filled dots below the bar, and the horizontal set-size bars on the left give each contrast's total number of down-regulated features."
#| fig-width: 9
#| fig-height: 6
draw_upset_or_empty(
significant_sets_by_direction(significant_contrasts, "decreased"),
"Decreased-abundance significant features are available for fewer than two contrasts."
)
```
:::
## Volcano
```{r}
#| label: fig-volcano
#| fig-cap: "Volcano plot of differential abundance: log2 fold-change (x) versus -log10(FDR) (y), one panel per contrast; each point is a feature, with dashed vertical lines at the fold-change threshold and a dashed horizontal line at the FDR threshold."
#| fig-width: 9
#| fig-height: 7
if (!is.null(contrast_plotter)) {
volcano_plot <- safe_plot(contrast_plotter$volcano()$FDR, "Could not draw volcano plot.")
save_if_ggplot(volcano_plot, "volcano", width = 9, height = 7)
show_plot(volcano_plot)
} else {
show_plot(empty_report_plot("No contrast results are available."))
}
```
## MA Plot
```{r}
#| label: fig-ma
#| fig-cap: "MA plots of differential abundance (two stacked panels, one row of panels per contrast): each point is a feature with its log2 fold-change on the y-axis; the top panel plots this against mean transformed abundance and the bottom panel against abundance rank, with a horizontal line at zero fold-change and dashed lines at the fold-change threshold."
#| fig-width: 8
#| fig-height: 9
if (!is.null(contrast_plotter)) {
ma_plot <- safe_plot(contrast_plotter$ma_plot(rank = FALSE), "Could not draw MA plot.")
ma_rank_plot <- safe_plot(
contrast_plotter$ma_plot(fc = params$diff_threshold, rank = TRUE),
"Could not draw ranked MA plot."
)
save_if_ggplot(ma_plot, "ma_plot", width = 8, height = 5)
save_if_ggplot(ma_rank_plot, "ma_rank_plot", width = 8, height = 5)
gridExtra::grid.arrange(ma_plot, ma_rank_plot, ncol = 1)
} else {
show_plot(empty_report_plot("No contrast results are available."))
}
```
## Significant Features
```{r}
#| label: fig-significant-protein-heatmap
#| fig-cap: "Heatmap of the significant differential-abundance features: rows are features called significant in at least one contrast, columns are samples, and each cell encodes the per-feature (row) z-score of the transformed abundance (row labels are shown only when at most 30 features are significant)."
#| fig-width: 9
#| fig-height: 8
if (nrow(significant_contrasts) > 0 && !is.null(feature_col)) {
significant_lfq <- lfq_transformed$get_copy()
significant_lfq <- significant_lfq$get_subset(significant_contrasts)
row_names <- tryCatch(significant_lfq$hierarchy_counts()[[2]] <= 30, error = function(e) FALSE)
significant_heatmap <- safe_plot(
significant_lfq$get_Plotter()$heatmap(rownames = row_names),
"Could not draw significant feature heatmap."
)
show_plot(significant_heatmap)
} else {
show_plot(empty_report_plot("No significant differential-abundance results are available."))
}
```
```{r}
#| label: tbl-significant-features
dt_table(significant_contrasts, caption = "Significant differential abundance results", page_length = 25)
```
## Distributions
::: {.callout-note collapse="true"}
## Why look at the fold-change and p-value distributions?
These two histograms are model diagnostics read across all features at once, not per feature. The **fold-change** distribution should be centred near zero, because most features are not differentially abundant; a systematic shift away from zero points to a normalization or experimental-design problem rather than real biology. The **p-value** (or FDR/BFDR) distribution should be roughly uniform under the null hypothesis, with a peak near zero when true effects are present; a peak near one, or a strongly non-uniform shape, indicates that the statistical model does not describe the data well (for example unmodelled variance or outliers), so the significance calls should be interpreted with caution.
:::
```{r}
#| label: fig-fc-pvalue
#| fig-cap: >-
#| Distributions of the differential-abundance statistics across all features.
#| Left: histogram of the estimated log2 fold-changes (x-axis log2
#| fold-change, y-axis feature count), expected to centre near zero because
#| most features are not differentially abundant. Right: histogram of the
#| significance score (p-value, or FDR/BFDR depending on the model; x-axis
#| 0-1, y-axis feature count), expected to be roughly uniform with a peak near
#| zero when true effects are present.
#| fig-width: 10
#| fig-height: 5
if (!is.null(contrast_plotter)) {
fc_hist <- safe_plot(contrast_plotter$histogram_diff(), "Could not draw fold-change histogram.")
histograms <- contrast_plotter$histogram()
score_col <- dplyr::case_when(
"p.value" %in% colnames(contrast_table) ~ "p.value",
"FDR" %in% colnames(contrast_table) ~ "FDR",
"BFDR" %in% colnames(contrast_table) ~ "BFDR",
TRUE ~ NA_character_
)
score_hist <- if (!is.na(score_col) && score_col %in% names(histograms)) {
safe_plot(histograms[[score_col]], paste("Could not draw", score_col, "histogram."))
} else {
empty_report_plot("No p-value or FDR histogram is available.")
}
score_name <- if (is.na(score_col)) "score" else score_col
save_if_ggplot(fc_hist, "fold_change_histogram")
save_if_ggplot(score_hist, paste0(score_name, "_histogram"))
gridExtra::grid.arrange(fc_hist, score_hist, nrow = 1)
} else {
show_plot(empty_report_plot("No contrast results are available."))
}
```
:::
# Result Table
```{r}
#| label: tbl-contrast-results
dt_table(contrast_table, caption = "Differential abundance results", page_length = 25)
```
# Session Info
::: {.panel-tabset}
## Report provenance
```{r}
#| label: tbl-report-provenance
#| tbl-cap: "Compact report provenance, including the source input-data reference."
knitr::kable(prolfquapp:::.report_provenance_table(report_provenance))
```
```{r}
#| label: tbl-bfabric
urls <- meta$bfabric_urls
url_df <- data.frame(field = character(), url = character())
if (is.list(urls) && length(urls) > 0) {
url_values <- vapply(urls, function(x) {
if (is.null(x) || length(x) == 0 || all(is.na(x))) {
return(NA_character_)
}
as.character(x[[1]])
}, character(1))
url_df <- data.frame(field = names(url_values), url = unname(url_values), stringsAsFactors = FALSE)
}
url_df <- url_df[!is.na(url_df$url) & nzchar(url_df$url), , drop = FALSE]
if (nrow(url_df) > 0) {
url_df$field <- html_escape(url_df$field)
url_df$url <- html_link(url_df$url)
}
dt_table(url_df, caption = "B-Fabric links", page_length = 5, escape = FALSE)
```
## R session info
```{r}
sessionInfo()
```
:::
:::
