---
title: "Using ggpicrust2"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Using ggpicrust2}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)
```

## Introduction

`ggpicrust2` provides a practical workflow for PICRUSt2 downstream analysis:

- convert KO profiles to pathway-level abundance when needed
- run differential abundance analysis with multiple methods
- annotate and visualize pathway-level results
- inspect taxa-level contribution using PICRUSt2 per-sequence outputs
- perform GSEA when pathway-set analysis is more appropriate than single-feature testing

This vignette focuses on the general package workflow. For a deeper GSEA walkthrough, see the dedicated `gsea_analysis` vignette.

## Installation and example data

Install the package and the optional backends used by this tutorial:

```{r installation, eval = FALSE}
install.packages(c("ggpicrust2", "MicrobiomeStat", "BiocManager"))
BiocManager::install(c("KEGGREST", "limma"))
```

```{r example-data, eval = FALSE}
library(ggpicrust2)
library(tibble)
data("ko_abundance")
data("metadata")
alpha <- 0.05
```

Both workflows below use the same data, LinDA method, and adjusted p-value
threshold. LinDA accepts the continuous abundance estimates produced by
`ko2kegg_abundance()`, but its default count-data winsorization rounds values
inside MicrobiomeStat. In a direct `pathway_daa()` call, set
`linda_winsor = FALSE` to preserve fractional values and
`linda_adaptive = FALSE` for fixed pseudo-count handling. A fixed pseudo-count
still depends on the units of the input when zeros are present.
Count-based methods such as ALDEx2 require integer input and the package
rounds non-integer input with a warning. Choose the
method for its assumptions and your study design, not to obtain significance.
KEGG pathway annotation requires internet access and `KEGGREST`.

## One-command workflow

```{r one-command, eval = FALSE}
results <- ggpicrust2(
  data = ko_abundance,
  metadata = metadata,
  group = "Environment",
  pathway = "KO",
  daa_method = "LinDA",
  ko_to_kegg = TRUE,
  order = "pathway_class",
  p_values_bar = TRUE,
  p_values_threshold = alpha,
  x_lab = "pathway_name"
)

# A method's plot is NULL when no pathways can be plotted.
results[[1]]$plot
head(results[[1]]$results)
```

## Stepwise pathway workflow

Run the installation and example-data setup above first. This workflow uses
the same analysis settings as the one-command workflow.

### Convert KO abundance to KEGG pathway abundance

```{r ko-to-kegg, eval = FALSE}
kegg_pathway_abundance <- ko2kegg_abundance(data = ko_abundance)
head(kegg_pathway_abundance[, 1:3])
```

### Match group labels to sample identifiers

`pathway_daa()`, `pathway_heatmap()`, and `pathway_pca()` accept `metadata`
and a column name, such as `group = "Environment"`. In contrast,
`pathway_errorbar()` has no metadata argument: its capitalized `Group`
parameter requires one group label per abundance column, not a column name.

The bundled metadata and abundance table have different sample orders.
Name the group vector with sample IDs so that plotting aligns labels to the
correct samples. For your own data, replace `sample_name` and `Environment`
with your sample-ID and grouping columns. Do not pass an unnamed metadata
column unless you have already verified its order against the abundance columns.

```{r sample-groups, eval = FALSE}
stopifnot(
  !anyNA(metadata$sample_name),
  !anyDuplicated(metadata$sample_name),
  setequal(colnames(kegg_pathway_abundance), metadata$sample_name)
)
sample_groups <- setNames(metadata$Environment, metadata$sample_name)
```

### Run differential abundance analysis

```{r daa, eval = FALSE}
daa_results <- pathway_daa(
  abundance = kegg_pathway_abundance,
  metadata = metadata,
  group = "Environment",
  daa_method = "LinDA"
)

head(daa_results)
```

### Annotate pathway results

```{r annotation, eval = FALSE}
annotated_daa <- pathway_annotation(
  pathway = "KO",
  daa_results_df = daa_results,
  ko_to_kegg = TRUE,
  p_adjust_threshold = alpha
)

head(annotated_daa)
```

### Visualize pathway-level results

`ko_to_kegg = TRUE` is required here too: the rows are KEGG pathways and
use `pathway_name` annotations. This flag does not reconvert the abundance
matrix in the plotting function.

```{r errorbar, eval = FALSE}
sig_pathways <- unique(annotated_daa$feature[
  !is.na(annotated_daa$p_adjust) & annotated_daa$p_adjust < alpha
])

p <- NULL
if (length(sig_pathways) > 0) {
  p <- pathway_errorbar(
    abundance = kegg_pathway_abundance,
    daa_results_df = annotated_daa,
    Group = sample_groups,
    ko_to_kegg = TRUE,
    p_values_threshold = alpha,
    order = "pathway_class",
    x_lab = "pathway_name"
  )
} else {
  message("No pathways pass the adjusted p-value threshold; skipping the error bar plot.")
}
p
```

No significant pathways is a valid analysis outcome. Keep the results table;
do not increase the threshold or change methods just to produce a plot.
Missing KEGG annotations can also prevent plotting even when significant
results exist; check the annotation warnings separately.

```{r heatmap-pca, eval = FALSE}
if (length(sig_pathways) > 0) {
  pathway_heatmap(
    abundance = kegg_pathway_abundance[sig_pathways, , drop = FALSE],
    metadata = metadata,
    group = "Environment"
  )
}

pathway_pca(
  abundance = kegg_pathway_abundance,
  metadata = metadata,
  group = "Environment"
)
```

### Using ALDEx2 instead

ALDEx2 is an optional Bioconductor dependency. For two groups it returns
both Welch and Wilcoxon results, so select one test before annotation and
plotting. Otherwise a pathway can occur twice with different p-values.
Specify the test in advance. ALDEx2 uses Monte Carlo sampling, so set a seed
for reproducibility; results can still differ across package versions.

```{r aldex2-alternative, eval = FALSE}
# Install once with BiocManager::install("ALDEx2").
set.seed(207)
aldex_results <- pathway_daa(
  abundance = kegg_pathway_abundance,
  metadata = metadata,
  group = "Environment",
  daa_method = "ALDEx2"
)
daa_results <- aldex_results[
  aldex_results$method == "ALDEx2_Welch's t test", , drop = FALSE
]
```

Then rerun the annotation and visualization steps with this `daa_results`.
For more than two groups or multiple contrasts, inspect `method`, `group1`,
and `group2` and select the supported test/contrast explicitly. An ALDEx2
run can have no adjusted p-values below 0.05 even when LinDA finds some;
these are different statistical procedures, not equivalent plotting modes.

## Taxa contribution workflow

PICRUSt2 contribution files attribute predicted functional abundance to taxa. `ggpicrust2` supports both gene-family-level and pathway-level contribution workflows.

### Run a synthetic contribution example

This small example illustrates input schemas and aggregation. It is not a
biological result, and its sample IDs are separate from the bundled KO dataset.

```{r contrib-example, eval = FALSE}
contrib_input <- expand.grid(
  sample = paste0("S", 1:4),
  function_id = c("K00001", "K00002"),
  taxon = c("ASV1", "ASV2"),
  stringsAsFactors = FALSE
)
contrib_input$taxon_function_abun <- seq_len(nrow(contrib_input))
contrib_data <- read_contrib_file(data = contrib_input)
contrib_metadata <- data.frame(
  sample_name = paste0("S", 1:4),
  Environment = rep(c("Control", "Treatment"), each = 2)
)
taxonomy <- data.frame(
  ASV = c("ASV1", "ASV2"),
  Genus = c("ExampleGenusA", "ExampleGenusB")
)
taxa_contrib <- aggregate_taxa_contributions(
  contrib_data, taxonomy = taxonomy, tax_level = "Genus", top_n = 2
)
head(taxa_contrib)
```

The aggregation sums the selected contribution column. By default it prefers
`norm_taxon_function_contrib` when supplied; this example supplies raw
`taxon_function_abun` only. A percentage bar subsequently normalizes within
each sample/function, so its heights describe the taxonomic composition of
that function, not a between-function abundance comparison.

```{r contrib-plots, eval = FALSE}
taxa_contribution_bar(
  contrib_agg = taxa_contrib,
  metadata = contrib_metadata,
  group = "Environment",
  facet_by = "function"
)
taxa_contribution_heatmap(contrib_agg = taxa_contrib, n_functions = 2)
```

### Read your own PICRUSt2 files

For real data, replace the synthetic input with one of the readers below and
use metadata and taxonomy for those same samples and taxa:

```r
# KO/gene-family contributions:
# contrib_data <- read_contrib_file("pred_metagenome_contrib.tsv")
# Pathway contributions (often MetaCyc):
# contrib_data <- read_pathway_contrib_file("path_abun_contrib.tsv.gz")
# Wide stratified abundance:
# contrib_data <- read_strat_file("pred_metagenome_strat.tsv")
```

When optional `daa_results_df` or `pathway_ids` filters contain KEGG pathway IDs
and the contribution table is KO-level, the function expands those pathway IDs
to KO members and retains matching KO rows. It does not turn KO contributions
into pathway contributions. The output `function_id` remains a KO identifier.
For pathway-level MetaCyc contributions, matching MetaCyc IDs are filtered
directly. Do not interpret a member-KO filter as independent evidence that a
taxon drives a reconstructed pathway's activity.

For pathway-level data, use matching pathway annotations. For example:

```{r pathway-contrib-example, eval = FALSE}
path_input <- contrib_input
path_input$function_id <- ifelse(path_input$function_id == "K00001",
                                  "GLYCOLYSIS", "PWY-5484")
path_data <- read_pathway_contrib_file(data = path_input)
path_taxa_contrib <- aggregate_taxa_contributions(
  path_data, taxonomy = taxonomy, tax_level = "Genus", top_n = 2
)
pathway_annotation_df <- pathway_annotation(
  data = data.frame(function_id = unique(path_taxa_contrib$function_id)),
  pathway = "MetaCyc"
)
```

## GSEA workflow

Use GSEA when you want pathway-set level inference from KO or EC abundance rather than testing each pathway independently.

```{r gsea, eval = FALSE}
gsea_results <- pathway_gsea(
  abundance = ko_abundance %>% column_to_rownames("#NAME"),
  metadata = metadata,
  group = "Environment",
  pathway_type = "KEGG",
  method = "camera"
)

annotated_gsea <- gsea_pathway_annotation(
  gsea_results = gsea_results,
  pathway_type = "KEGG"
)

visualize_gsea(
  gsea_results = annotated_gsea,
  plot_type = "barplot",
  n_pathways = 15
)
```

For a method-by-method GSEA explanation, covariate adjustment, and comparison with DAA, see the `gsea_analysis` vignette.

## Summary

The package is easiest to use when you choose the shortest path that matches your question:

- use `ggpicrust2()` for a fast default pathway workflow
- use the stepwise DAA functions when you need more control
- use the taxa contribution workflow when you need taxon-level attribution
- use `pathway_gsea()` when pathway-set enrichment is the primary question