ggpicrust2

🌟 **If you find ggpicrust2 helpful, please consider giving us a star on GitHub!** Your support greatly motivates us to improve and maintain this project. 🌟

ggpicrust2 is a comprehensive package designed to provide a seamless and intuitive solution for analyzing and interpreting the results of PICRUSt2 functional prediction. It offers a wide range of features, including pathway name/description annotations, advanced differential abundance (DA) methods, and visualization of DA results.

One of the newest additions to ggpicrust2 is the capability to compare the consistency and inconsistency across different DA methods applied to the same dataset. This feature allows users to assess the agreement and discrepancy between various methods when it comes to predicting and sequencing the metagenome of a particular sample. It provides valuable insights into the consistency of results obtained from different approaches and helps users evaluate the robustness of their findings.

By leveraging this functionality, researchers, data scientists, and bioinformaticians can gain a deeper understanding of the underlying biological processes and mechanisms present in their PICRUSt2 output data. This comparison of different methods enables them to make informed decisions and draw reliable conclusions based on the consistency evaluation of macrogenomic predictions or sequencing results for the same sample.

If you are interested in exploring and analyzing your PICRUSt2 output data, ggpicrust2 is a powerful tool that provides a comprehensive set of features, including the ability to assess the consistency and evaluate the performance of different methods applied to the same dataset.

News

🔬 Enhanced GSEA with Limma Camera/Fry Methods (v2.5.6)

We’ve significantly enhanced the pathway_gsea() function with limma’s competitive gene set testing methods:

• camera method (new default): Accounts for inter-gene

correlations, providing more reliable p-values than preranked methods (Wu et al., 2012)

• fry method: Fast rotation gene set test, efficient for large

gene set collections

• Covariate adjustment: Both camera and fry support adjustment for

confounding variables (age, sex, BMI, etc.)

• Contrast specification: Support for complex experimental designs

with multiple contrasts

This addresses concerns in the literature that preranked GSEA methods can produce “spectacularly wrong p-values” due to not accounting for inter-gene correlations.

📊 New Visualization Functions: Volcano Plot & Ridge Plot

We’ve added two new visualization functions for enhanced analysis and interpretation:

• pathway_volcano(): Creates publication-quality volcano plots for

differential abundance analysis, with smart label placement and color-coded significance categories.

• pathway_ridgeplot(): Creates ridge plots (joy plots) for GSEA

results interpretation, showing distribution of gene abundances within enriched pathways.

🔄 **Updated Reference Databases for Improved Pathway Annotation (v2.1.4)**

We’ve significantly enhanced the reference databases used for pathway annotation:

• EC reference data: Updated from 3,180 to 8,371 entries (163%

increase)

• KO reference data: Updated from 23,917 to 27,531 unique KO IDs

(15.4% increase)

These updates provide more comprehensive and accurate pathway annotations, especially for recently discovered enzymes and KEGG orthology entries. Users will experience improved coverage and precision in pathway analysis without needing to change any code.

🌟 **New Feature: Gene Set Enrichment Analysis (GSEA) for PICRUSt2 Data**

We’re excited to announce the addition of GSEA functionality to the ggpicrust2 package! This powerful new feature allows researchers to perform Gene Set Enrichment Analysis on PICRUSt2 predicted functional profiles, offering a more nuanced understanding of functional differences between conditions.

The new GSEA module includes:

• pathway_gsea(): Performs GSEA analysis on PICRUSt2 data
• visualize_gsea(): Creates various visualizations including

enrichment plots, dot plots, network plots, and heatmaps

• comparegseadaa(): Compares GSEA and differential abundance

analysis results

• gseapathwayannotation(): Annotates GSEA results with pathway

information

These new functions complement our existing differential abundance analysis tools, providing researchers with multiple approaches to analyze functional profiles.

🧫 **New Feature: Taxa Contribution Workflow for PICRUSt2 Per-sequence Outputs**

ggpicrust2 now supports parsing and visualizing PICRUSt2 per-sequence contribution outputs:

• readcontribfile(): Reads predmetagenomecontrib.tsv
• readstratfile(): Reads predmetagenomestrat.tsv
• aggregatetaxacontributions(): Aggregates contributions by taxon,

sample, and function

• taxacontributionbar(): Creates stacked bar plots of taxon-level

contributions

• taxacontributionheatmap(): Summarizes contribution patterns across

taxa and functions

This workflow makes it possible to move from pathway-level significance to an interpretable answer for which taxa are driving those pathway shifts.

🌟 Also Check Out: mLLMCelltype

We’re excited to introduce mLLMCelltype, our innovative framework for single-cell RNA sequencing data annotation. This iterative multi-LLM consensus framework leverages the collective intelligence of multiple large language models (including GPT-4o/4.1, Claude-3.7/3.5, Gemini-2.0, Grok-3, and others) to significantly improve cell type annotation accuracy while providing transparent uncertainty quantification.

mLLMCelltype addresses critical challenges in scRNA-seq analysis through its unique architecture:

• Multi-LLM Consensus: Overcomes single-model limitations by

harnessing diverse LLMs’ collective intelligence

• Structured Deliberation: Enables models to share reasoning,

evaluate evidence, and refine annotations through collaborative discussion

• Transparent Uncertainty Metrics: Provides quantitative measures to

identify ambiguous cell populations requiring expert review

• Hallucination Reduction: Suppresses inaccurate predictions through

cross-model critical evaluation

• No Reference Dataset Required: Performs accurate annotation

without pre-training or reference data

For researchers working with single-cell data, mLLMCelltype offers a powerful new approach to cell type annotation. Learn more about its capabilities and methodology on GitHub: mLLMCelltype Repository.

We appreciate your support and interest in our tools and look forward to seeing how they can enhance your research.

Table of Contents

• Citation
• Installation
• Project Links
• Workflow
• Output
• Function Details

- ko2keggabundance() - pathwaydaa() - comparedaaresults() - pathwayannotation() - pathwayerrorbar() - pathwayheatmap() - pathwaypca() - comparemetagenomeresults() - taxa contribution workflow - pathwaygsea() - visualizegsea() - comparegseadaa() - gseapathwayannotation()

• FAQ
• Author’s Other Projects

Citation

If you use ggpicrust2 in your research, please cite the following paper:

Chen Yang and others. (2023). ggpicrust2: an R package for PICRUSt2 predicted functional profile analysis and visualization. Bioinformatics, btad470. DOI link

The package citation is also available directly in R:

r
citation("ggpicrust2")

Installation

You can install the development version of ggpicrust2 from GitHub with:

r
install.packages("devtools")
devtools::install_github("cafferychen777/ggpicrust2")

Dependent CRAN Packages

| Package | Description | |----|----| | aplot | Create interactive plots | | dplyr | A fast consistent tool for working with data frame like objects both in memory and out of memory | | ggplot2 | An implementation of the Grammar of Graphics in R | | grid | A rewrite of the graphics layout capabilities of R | | MicrobiomeStat | Statistical analysis of microbiome data | | readr | Read rectangular data (csv tsv fwf) into R | | stats | The R Stats Package | | tibble | Simple Data Frames | | tidyr | Easily tidy data with spread() and gather() functions | | ggprism | Interactive 3D plots with ‘prism’ graphics | | cowplot | Streamlined Plot Theme and Plot Annotations for ‘ggplot2’ | | ggforce | Easily add secondary axes, zooms, and image overlays to ‘ggplot2’ | | ggplotify | Convert complex plots into ‘grob’ or ‘ggplot’ objects | | magrittr | A Forward-Pipe Operator for R | | utils | The R Utils Package |

Optional Bioconductor Packages

The package works with a minimal CRAN installation, but several workflows rely on Bioconductor packages that should be installed only when you use the corresponding analysis methods or visualizations.

| Package | Description | |----|----| | phyloseq | Handling and analysis of high-throughput microbiome census data | | ALDEx2 | Differential abundance analysis of taxonomic and functional features | | SummarizedExperiment | SummarizedExperiment container for storing data and metadata together | | Biobase | Base functions for Bioconductor | | devtools | Tools to make developing R packages easier | | ComplexHeatmap | Making Complex Heatmaps in R | | BiocGenerics | S4 generic functions for Bioconductor | | BiocManager | Access the Bioconductor Project Package Repositories | | metagenomeSeq | Statistical analysis for sparse high-throughput sequencing | | Maaslin2 | Tools for microbiome analysis | | edgeR | Empirical Analysis of Digital Gene Expression Data in R | | lefser | R implementation of the LEfSE method for microbiome biomarker discovery | | limma | Linear Models for Microarray and RNA-Seq Data | | KEGGREST | R Interface to KEGG REST API | | DESeq2 | Differential gene expression analysis using RNA-seq data |

r
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
pkgs <- c("phyloseq", "ALDEx2", "SummarizedExperiment", "Biobase", "devtools",           "ComplexHeatmap", "BiocGenerics", "BiocManager", "metagenomeSeq",           "Maaslin2", "edgeR", "lefser", "limma", "KEGGREST", "DESeq2")
for (pkg in pkgs) {   if (!requireNamespace(pkg, quietly = TRUE))     BiocManager::install(pkg) }

Project Links

Use the project resources below for stable updates and support:

• Package website:
• Source repository:
• Issue tracker:
• Release notes: NEWS.md

Workflow

The easiest way to analyze the PICRUSt2 output is using ggpicrust2() function. The main pipeline can be run with ggpicrust2() function.

ggpicrust2() integrates ko abundance to kegg pathway abundance conversion, annotation of pathway, differential abundance (DA) analysis, part of DA results visualization. When you have trouble running ggpicrust2(), you can debug it by running a separate function, which will greatly increase the speed of your analysis and visualization.

ggpicrust2()

You can download the example dataset from the provided Github link and Google Drive link or use the dataset included in the package.

r
If you want to analyze the abundance of KEGG pathways instead of KO within the pathway, please set kotokegg to TRUE.
KEGG pathways typically have more descriptive explanations.
library(readr) library(ggpicrust2) library(tibble) library(tidyverse) library(ggprism) library(patchwork)
Load necessary data: abundance data and metadata
abundancefile <- "path/to/your/abundancefile.tsv"
metadata <- read_delim(
    "path/to/your/metadata.txt",
    delim = "\t",
    escape_double = FALSE,
    trim_ws = TRUE
)
Run ggpicrust2 with input file path
resultsfileinput <- ggpicrust2(file = abundance_file,
                                 metadata = metadata,
                                 group = "yourgroupcolumn", # For example dataset, group = "Environment"
                                 pathway = "KO",
                                 daa_method = "LinDA",
                                 kotokegg = TRUE,
                                 order = "pathway_class",
                                 pvaluesbar = TRUE,
                                 xlab = "pathwayname")
Run ggpicrust2 with imported data.frame
abundancedata <- readdelim(abundancefile, delim = "\t", colnames = TRUE, trim_ws = TRUE)
Run ggpicrust2 with input data
resultsdatainput <- ggpicrust2(data = abundance_data,
                                 metadata = metadata,
                                 group = "yourgroupcolumn", # For example dataset, group = "Environment"
                                 pathway = "KO",
                                 daa_method = "LinDA",
                                 kotokegg = TRUE,
                                 order = "pathway_class",
                                 pvaluesbar = TRUE,
                                 xlab = "pathwayname")
Access the plot and results dataframe for the first DA method
exampleplot <- resultsfile_input[[1]]$plot
exampleresults <- resultsfile_input[[1]]$results
Use the example data in ggpicrust2 package
data(ko_abundance)
data(metadata)
resultsfileinput <- ggpicrust2(data = ko_abundance,
                                 metadata = metadata,
                                 group = "Environment",
                                 pathway = "KO",
                                 daa_method = "LinDA",
                                 kotokegg = TRUE,
                                 order = "pathway_class",
                                 pvaluesbar = TRUE,
                                 xlab = "pathwayname")
Analyze the EC or MetaCyc pathway
data(metacyc_abundance)
resultsfileinput <- ggpicrust2(data = metacyc_abundance,
                                 metadata = metadata,
                                 group = "Environment",
                                 pathway = "MetaCyc",
                                 daa_method = "LinDA",
                                 kotokegg = FALSE,
                                 order = "group",
                                 pvaluesbar = TRUE,
                                 x_lab = "description")
resultsfileinput[[1]]$plot
resultsfileinput[[1]]$results

If an error occurs with ggpicrust2, please use the following workflow.

r
library(readr)
library(ggpicrust2)
library(tibble)
library(tidyverse)
library(ggprism)
library(patchwork)
If you want to analyze KEGG pathway abundance instead of KO within the pathway, turn kotokegg to TRUE.
KEGG pathways typically have more explainable descriptions.
Load metadata as a tibble
data(metadata)
metadata <- readdelim("path/to/your/metadata.txt", delim = "\t", escapedouble = FALSE, trim_ws = TRUE)
Load KEGG pathway abundance
data(kegg_abundance)
keggabundance <- ko2keggabundance("path/to/your/predmetagenomeunstrat.tsv")
Perform pathway differential abundance analysis (DAA) using ALDEx2 method.
Please change group to "yourgroupcolumn" if you are not using example dataset.
From v2.5.14 ALDEx2 results include effectsize, diffbtw, log2foldchange,
rab_all, and overlap columns by default (via ALDEx2::aldex.effect()), matching
the other DAA methods that return log2 fold changes by default. Pass
includeeffectsize = FALSE to skip that step.
daaresultsdf <- pathwaydaa(abundance = keggabundance, metadata = metadata, group = "Environment", daa_method = "ALDEx2", select = NULL, reference = NULL)
Filter results for ALDEx2_Wilcoxon rank test method
Please check the unique(daaresultsdf$method) and choose one
daasubmethodresultsdf <- daaresultsdf[daaresultsdf$method == "ALDEx2_Wilcoxon rank test", ]
Ranking by |log2foldchange| is generally more biologically informative than
ranking by p-value, especially for large datasets where small effects can
reach statistical significance without being biologically meaningful.
tophits <- daasubmethodresultsdf[order(-abs(daasubmethodresultsdf$log2fold_change)), ]
Annotate pathway results using KO to KEGG conversion
daaannotatedsubmethodresultsdf <- pathwayannotation(pathway = "KO", daaresultsdf = daasubmethodresultsdf, kotokegg = TRUE)
Generate pathway error bar plot
Please change Group to metadata$yourgroupcolumn if you are not using example dataset
p <- pathwayerrorbar(abundance = keggabundance, daaresultsdf = daaannotatedsubmethodresultsdf, Group = metadata$Environment, pvaluesthreshold = 0.05, order = "pathwayclass", select = NULL, kotokegg = TRUE, pvaluebar = TRUE, colors = NULL, xlab = "pathwayname")
If you want to analyze EC, MetaCyc, and KO without conversions, turn kotokegg to FALSE.
Load metadata as a tibble
data(metadata)
metadata <- readdelim("path/to/your/metadata.txt", delim = "\t", escapedouble = FALSE, trim_ws = TRUE)
Load KO abundance as a data.frame
data(ko_abundance)
koabundance <- read.delim("path/to/your/predmetagenome_unstrat.tsv")
Perform pathway DAA using ALDEx2 method
Please change columntorownames() to the feature column if you are not using example dataset
Please change group to "yourgroupcolumn" if you are not using example dataset
ALDEx2 effect size columns (effectsize, diffbtw, log2foldchange, rab_all,
overlap) are included by default; see the section above.
daaresultsdf <- pathwaydaa(abundance = koabundance %>% columntorownames("#NAME"), metadata = metadata, group = "Environment", daa_method = "ALDEx2", select = NULL, reference = NULL)
Filter results for ALDEx2_Wilcoxon rank test method
daasubmethodresultsdf <- daaresultsdf[daaresultsdf$method == "ALDEx2_Wilcoxon rank test", ]
Annotate pathway results without KO to KEGG conversion
daaannotatedsubmethodresultsdf <- pathwayannotation(pathway = "KO", daaresultsdf = daasubmethodresultsdf, kotokegg = FALSE)
Generate pathway error bar plot
Please change columntorownames() to the feature column
Please change Group to metadata$yourgroupcolumn if you are not using example dataset
p <- pathwayerrorbar(abundance = koabundance %>% columntorownames("#NAME"), daaresultsdf = daaannotatedsubmethodresultsdf, Group = metadata$Environment, pvalues_threshold = 0.05, order = "group",
select = daaannotatedsubmethodresultsdf %>% arrange(padjust) %>% slice(1:20) %>% dplyr::select(feature) %>% pull(),
kotokegg = FALSE,
pvaluebar = TRUE,
colors = NULL,
x_lab = "description")
Workflow for MetaCyc Pathway and EC
Load MetaCyc pathway abundance and metadata
data("metacyc_abundance")
data("metadata")
Perform pathway DAA using LinDA method
Please change columntorownames() to the feature column if you are not using example dataset
Please change group to "yourgroupcolumn" if you are not using example dataset
metacycdaaresultsdf <- pathwaydaa(abundance = metacycabundance %>% columntorownames("pathway"), metadata = metadata, group = "Environment", daamethod = "LinDA")
Annotate MetaCyc pathway results without KO to KEGG conversion
metacycdaaannotatedresultsdf <- pathwayannotation(pathway = "MetaCyc", daaresultsdf = metacycdaaresultsdf, kotokegg = FALSE)
Generate pathway error bar plot
Please change columntorownames() to the feature column
Please change Group to metadata$yourgroupcolumn if you are not using example dataset
pathwayerrorbar(abundance = metacycabundance %>% columntorownames("pathway"), daaresultsdf = metacycdaaannotatedresultsdf, Group = metadata$Environment, kotokegg = FALSE, pvaluesthreshold = 0.05, order = "group", select = NULL, pvaluebar = TRUE, colors = NULL, x_lab = "description")
Generate pathway heatmap
Please change columntorownames() to the feature column if you are not using example dataset
Please change group to "yourgroupcolumn" if you are not using example dataset
featurewithp0.05 <- metacycdaaresultsdf %>% filter(p_adjust < 0.05)
pathwayheatmap(abundance = metacycabundance %>% filter(pathway %in% featurewithp0.05$feature) %>% columnto_rownames("pathway"), metadata = metadata, group = "Environment")
Generate pathway PCA plot
Please change columntorownames() to the feature column if you are not using example dataset
Please change group to "yourgroupcolumn" if you are not using example dataset
pathwaypca(abundance = metacycabundance %>% columntorownames("pathway"), metadata = metadata, group = "Environment")
Run pathway DAA for multiple methods
Please change columntorownames() to the feature column if you are not using example dataset
Please change group to "yourgroupcolumn" if you are not using example dataset
methods <- c("ALDEx2", "DESeq2", "edgeR")
daaresultslist <- lapply(methods, function(method) {
  pathwaydaa(abundance = metacycabundance %>% columntorownames("pathway"), metadata = metadata, group = "Environment", daa_method = method)
})
Compare results across different methods
comparisonresults <- comparedaaresults(daaresultslist = daaresultslist, methodnames = c("ALDEx2Welch's t test", "ALDEx2Wilcoxon rank test", "DESeq2", "edgeR"))

Output

The typical output of the ggpicrust2 is like this.

function details

ko2kegg_abundance()

KEGG Orthology(KO) is a classification system developed by the Kyoto Encyclopedia of Genes and Genomes (KEGG) data-base(Kanehisa et al., 2022). It uses a hierarchical structure to classify enzymes based on the reactions they catalyze. To better understand pathways’ role in different groups and classify the pathways, the KO abundance table needs to be converted to KEGG pathway abundance. But PICRUSt2 removes the function from PICRUSt. ko2kegg_abundance() can help convert the table.

r
Sample usage of the ko2kegg_abundance function
devtools::install_github('cafferychen777/ggpicrust2')
library(ggpicrust2)
Assume that the KO abundance table is stored in a file named "ko_abundance.tsv"
koabundancefile <- "ko_abundance.tsv"
Convert KO abundance to KEGG pathway abundance
keggabundance <- ko2keggabundance(file = koabundancefile)
Alternatively, if the KO abundance data is already loaded as a data frame named "ko_abundance"
data("ko_abundance")
keggabundance <- ko2keggabundance(data = ko_abundance)
The resulting kegg_abundance data frame can now be used for further analysis and visualization.

pathway_daa()

Differential abundance (DA) analysis plays a major role in PICRUSt2 downstream analysis. pathway_daa() integrates the main DA methods used for predicted functional profiles, excluding ANCOM and ANCOMBC. It includes ALDEx2 (Fernandes et al., 2013), DESeq2 (Love et al., 2014), Maaslin2 (Mallick et al., 2021), LinDA (Zhou et al., 2022), edgeR (Robinson et al., 2010), limma voom (Ritchie et al., 2015), metagenomeSeq (Paulson et al., 2013), and Lefser (Segata et al., 2011).

r
The abundance table is recommended to be a data.frame rather than a tibble.
The abundance table should have feature names or pathway names as row names, and sample names as column names.
You can use the output of ko2kegg_abundance
koabundancefile <- "path/to/your/predmetagenomeunstrat.tsv"
keggabundance <- ko2keggabundance(koabundancefile) # Or use data(kegg_abundance)
metadata <- readdelim("path/to/your/metadata.txt", delim = "\t", escapedouble = FALSE, trim_ws = TRUE)
The default DAA method is "ALDEx2"
Please change group to "yourgroupcolumn" if you are not using example dataset
daaresultsdf <- pathwaydaa(abundance = keggabundance, metadata = metadata, group = "Environment", daamethod = "LinDA", select = NULL, padjust_method = "BH", reference = NULL)
If you have more than 3 group levels and want to use the LinDA, limma voom, or Maaslin2 methods, you should provide a reference.
metadata <- readdelim("path/to/your/metadata.txt", delim = "\t", escapedouble = FALSE, trim_ws = TRUE)
Please change group to "yourgroupcolumn" if you are not using example dataset
daaresultsdf <- pathwaydaa(abundance = keggabundance, metadata = metadata, group = "Group", daamethod = "LinDA", select = NULL, padjust_method = "BH", reference = "Harvard BRI")
Other example
data("metacyc_abundance")
data("metadata")
metacycdaaresultsdf <- pathwaydaa(abundance = metacycabundance %>% columntorownames("pathway"), metadata = metadata, group = "Environment", daamethod = "LinDA", select = NULL, padjustmethod = "BH", reference = NULL)

comparedaaresults()

r
library(ggpicrust2)
library(tidyverse)
data("metacyc_abundance")
data("metadata")
Run pathway_daa function for multiple methods
Please change columntorownames() to the feature column if you are not using example dataset
Please change group to "yourgroupcolumn" if you are not using example dataset
methods <- c("ALDEx2", "DESeq2", "edgeR")
daaresultslist <- lapply(methods, function(method) {
  pathwaydaa(abundance = metacycabundance %>% columntorownames("pathway"), metadata = metadata, group = "Environment", daa_method = method)
})
method_names <- c("ALDEx2","DESeq2", "edgeR") 
Compare results across different methods
 comparisonresults <- comparedaaresults(daaresultslist = daaresultslist, methodnames = method_names)

pathway_annotation()

**If you are in China and you are using kegg pathway annotation, Please make sure your internet can break through the firewall.**

New Feature (v2.1.4): The pathway_annotation() function now supports species-specific KEGG pathway annotation through the new organism parameter. You can specify KEGG organism codes (e.g., “hsa” for human, “eco” for E. coli) to get species-specific pathway information. If no organism is specified (default), the function retrieves generic KO information not specific to any organism.

Note: When kotokegg = TRUE, only pathways with padjust < padjust_threshold are sent to the KEGG API for annotation. The padjustthreshold parameter defaults to 0.05 and can be customized. When called from ggpicrust2(), this threshold is automatically set to match the pvaluesthreshold parameter for consistency.

r
Make sure to check if the features in daaresultsdf correspond to the selected pathway
Annotate KEGG Pathway
data("kegg_abundance")
data("metadata")
Please change group to "yourgroupcolumn" if you are not using example dataset
daaresultsdf <- pathwaydaa(abundance = keggabundance, metadata = metadata, group = "Environment", daa_method = "LinDA")
Generic KO to KEGG pathway annotation (not specific to any organism)
daaannotatedresultsdf <- pathwayannotation(pathway = "KO", daaresultsdf = daaresultsdf, kotokegg = TRUE)
Species-specific KEGG pathway annotation (e.g., for human)
humanannotatedresultsdf <- pathwayannotation(pathway = "KO", daaresultsdf = daaresultsdf, kotokegg = TRUE, organism = "hsa")
Species-specific KEGG pathway annotation (e.g., for E. coli)
ecoliannotatedresultsdf <- pathwayannotation(pathway = "KO", daaresultsdf = daaresultsdf, kotokegg = TRUE, organism = "eco")
Annotate KO
data("ko_abundance")
data("metadata")
Please change columntorownames() to the feature column if you are not using example dataset
Please change group to "yourgroupcolumn" if you are not using example dataset
daaresultsdf <- pathwaydaa(abundance = koabundance %>% columntorownames("#NAME"), metadata = metadata, group = "Environment", daa_method = "LinDA")
daaannotatedresultsdf <- pathwayannotation(pathway = "KO", daaresultsdf = daaresultsdf, kotokegg = FALSE)
Annotate KEGG
daaannotatedresultsdf <- pathwayannotation(pathway = "EC", daaresultsdf = daaresultsdf, kotokegg = FALSE)
Annotate MetaCyc Pathway
data("metacyc_abundance")
data("metadata")
Please change columntorownames() to the feature column if you are not using example dataset
Please change group to "yourgroupcolumn" if you are not using example dataset
metacycdaaresultsdf <- pathwaydaa(abundance = metacycabundance %>% columntorownames("pathway"), metadata = metadata, group = "Environment", daamethod = "LinDA")
metacycdaaannotatedresultsdf <- pathwayannotation(pathway = "MetaCyc", daaresultsdf = metacycdaaresultsdf, kotokegg = FALSE)

pathway_errorbar()

r
data("ko_abundance")
data("metadata")
keggabundance <- ko2keggabundance(data = koabundance) # Or use data(keggabundance)
Please change group to "yourgroupcolumn" if you are not using example dataset
daaresultsdf <- pathwaydaa(keggabundance, metadata = metadata, group = "Environment", daa_method = "LinDA")
daaannotatedresultsdf <- pathwayannotation(pathway = "KO", daaresultsdf = daaresultsdf, kotokegg = TRUE)
Please change Group to metadata$yourgroupcolumn if you are not using example dataset
p <- pathwayerrorbar(abundance = keggabundance,
           daaresultsdf = daaannotatedresults_df,
           Group = metadata$Environment,
           kotokegg = TRUE,
           pvaluesthreshold = 0.05,
           order = "pathway_class",
           select = NULL,
           pvaluebar = TRUE,
           colors = NULL,
           xlab = "pathwayname")
If you want to analysis the EC. MetaCyc. KO without conversions.
data("metacyc_abundance")
data("metadata")
metacycdaaresultsdf <- pathwaydaa(abundance = metacycabundance %>% columntorownames("pathway"), metadata = metadata, group = "Environment", daamethod = "LinDA")
metacycdaaannotatedresultsdf <- pathwayannotation(pathway = "MetaCyc", daaresultsdf = metacycdaaresultsdf, kotokegg = FALSE)
p <- pathwayerrorbar(abundance = metacycabundance %>% columntorownames("pathway"),
           daaresultsdf = metacycdaaannotatedresultsdf,
           Group = metadata$Environment,
           kotokegg = FALSE,
           pvaluesthreshold = 0.05,
           order = "group",
           select = NULL,
           pvaluebar = TRUE,
           colors = NULL,
           x_lab = "description")

pathway_heatmap()

In this section, we will demonstrate how to create a pathway heatmap using the pathway_heatmap function in the ggpicrust2 package. This function visualizes the relative abundance of pathways in different samples.

Use the fake dataset

r
Create example functional pathway abundance data
abundance_example <- matrix(rnorm(30), nrow = 3, ncol = 10)
colnames(abundance_example) <- paste0("Sample", 1:10)
rownames(abundance_example) <- c("PathwayA", "PathwayB", "PathwayC")
Create example metadata
Please change your sample id's column name to sample_name
metadataexample <- data.frame(samplename = colnames(abundance_example),
                               group = factor(rep(c("Control", "Treatment"), each = 5)))
Create a heatmap
pathwayheatmap(abundanceexample, metadata_example, "group")

Use the real dataset

r
library(tidyverse)
library(ggh4x)
library(ggpicrust2)
Load the data
data("metacyc_abundance")
Load the metadata
data("metadata")
Perform differential abundance analysis
metacycdaaresultsdf <- pathwaydaa(
  abundance = metacycabundance %>% columnto_rownames("pathway"),
  metadata = metadata,
  group = "Environment",
  daa_method = "LinDA"
)
Annotate the results
annotatedmetacycdaaresultsdf <- pathway_annotation(
  pathway = "MetaCyc",
  daaresultsdf = metacycdaaresults_df,
  kotokegg = FALSE
)
Filter features with p < 0.05
featurewithp0.05 <- metacycdaaresultsdf %>%
  filter(p_adjust < 0.05)
Create the heatmap
pathway_heatmap(
  abundance = metacyc_abundance %>%
    right_join(
      annotatedmetacycdaaresultsdf %>% select(all_of(c("feature","description"))),
      by = c("pathway" = "feature")
    ) %>%
    filter(pathway %in% featurewithp_0.05$feature) %>%
    select(-"pathway") %>%
    columntorownames("description"),
  metadata = metadata,
  group = "Environment"
)

pathway_pca()

In this section, we will demonstrate how to perform Principal Component Analysis (PCA) on functional pathway abundance data and create visualizations of the PCA results using the pathway_pca function in the ggpicrust2 package.

Use the fake dataset

r
Create example functional pathway abundance data
abundance_example <- matrix(rnorm(30), nrow = 3, ncol = 10)
colnames(keggabundanceexample) <- paste0("Sample", 1:10)
rownames(keggabundanceexample) <- c("PathwayA", "PathwayB", "PathwayC")
Create example metadata
metadataexample <- data.frame(samplename = colnames(keggabundanceexample),
                                group = factor(rep(c("Control", "Treatment"), each = 5)))
Perform PCA and create visualizations
pathwaypca(abundance = abundanceexample, metadata = metadata_example, "group")

Use the real dataset

r
Create example functional pathway abundance data
data("metacyc_abundance")
data("metadata")
pathwaypca(abundance = metacycabundance %>% columntorownames("pathway"), metadata = metadata, group = "Environment")

comparemetagenomeresults()

r
library(ComplexHeatmap)
set.seed(123)
First metagenome
metagenome1 <- abs(matrix(rnorm(1000), nrow = 100, ncol = 10))
rownames(metagenome1) <- paste0("KO", 1:100)
colnames(metagenome1) <- paste0("sample", 1:10)
Second metagenome
metagenome2 <- abs(matrix(rnorm(1000), nrow = 100, ncol = 10))
rownames(metagenome2) <- paste0("KO", 1:100)
colnames(metagenome2) <- paste0("sample", 1:10)
Put the metagenomes into a list
metagenomes <- list(metagenome1, metagenome2)
Define names
names <- c("metagenome1", "metagenome2")
Compare the same biological samples with paired inference
results <- comparemetagenomeresults(
  metagenomes,
  names,
  daa_method = "paired Wilcoxon",
  correlation_permutations = 999
)
Print the correlation matrix
print(results$correlation$cor_matrix)
Print raw and multiplicity-adjusted permutation p-values
print(results$correlation$p_matrix)
print(results$correlation$padjustmatrix)

taxa contribution workflow

The taxa contribution workflow connects PICRUSt2 contribution outputs with downstream pathway interpretation. Use readcontribfile() for predmetagenomecontrib.tsv, readpathwaycontrib_file() for pathway-level pathabuncontrib.tsv, or readstratfile() for predmetagenomestrat.tsv, then aggregate to the taxonomic level you want to visualize.

r
library(ggpicrust2)
Parse PICRUSt2 per-sequence contribution output
contribdata <- readcontribfile("predmetagenome_contrib.tsv")
Or parse pathway-level contribution output without running DA analysis
PICRUSt2 pathway contribution files are often gzipped and use MetaCyc IDs.
pathcontribdata <- readpathwaycontribfile("pathabun_contrib.tsv.gz")
Optional: use pathway-level DAA results to keep only significant pathways
data("kegg_abundance")
data("metadata")
daaresults <- pathwaydaa(   abundance = kegg_abundance,   metadata = metadata,   group = "Environment",   daa_method = "ALDEx2" )
Aggregate contributions to genus level and keep top taxa
taxacontrib <- aggregatetaxa_contributions(
  contribdata = contribdata,
  taxonomy = yourtaxonomytable,
  tax_level = "Genus",
  top_n = 10,
  daaresultsdf = daa_results
)
Visualize per-sample contributions
taxacontributionbar(
  contribagg = taxacontrib,
  metadata = metadata,
  group = "Environment",
  facet_by = "function"
)
Summarize mean contribution patterns across taxa and functions
taxacontributionheatmap(
  contribagg = taxacontrib,
  n_functions = 20
)
Pathway-level contribution workflow without DA filtering
pathtaxacontrib <- aggregatetaxacontributions(
  contribdata = pathcontrib_data,
  taxonomy = yourtaxonomytable,
  tax_level = "Genus",
  top_n = 10
)
pathwayannotationdf <- pathway_annotation(   data = data.frame(functionid = unique(pathtaxacontrib$functionid)),   pathway = "MetaCyc" )
taxacontributionheatmap(   contribagg = pathtaxa_contrib,   annotationdata = pathwayannotation_df,   n_functions = 20 )

pathway_gsea()

The pathway_gsea() function performs Gene Set Enrichment Analysis (GSEA) on PICRUSt2 predicted functional profiles. GSEA is a powerful method for identifying enriched pathways between different conditions, offering a more nuanced understanding of functional differences compared to traditional differential abundance analysis.

New in v2.5.6: The function now supports limma’s camera and fry methods, which provide more reliable p-values by accounting for inter-gene correlations (Wu et al., 2012). The camera method is now the default and recommended approach.

| Method | Type | Covariate Support | Description | |----|----|----|----| | camera (default) | Competitive | ✅ Yes | Recommended. Accounts for inter-gene correlations | | fry | Self-contained | ✅ Yes | Fast rotation test, efficient for large gene set collections | | fgsea | Preranked | ❌ No | Fast preranked GSEA (legacy) | | clusterProfiler | Preranked | ❌ No | Traditional GSEA implementation (legacy) |

r
library(ggpicrust2)
library(tidyverse)
Load example data
data("ko_abundance")
data("metadata")
Perform GSEA analysis with camera (recommended)
gsearesults <- pathwaygsea(
  abundance = koabundance %>% columnto_rownames("#NAME"),
  metadata = metadata,
  group = "Environment",
  method = "camera",   # Recommended: accounts for inter-gene correlations
  pathway_type = "KEGG",
  min_size = 5,
  max_size = 500
)
With covariate adjustment (powerful feature of camera/fry)
gsearesultsadjusted <- pathway_gsea(
  abundance = koabundance %>% columnto_rownames("#NAME"),
  metadata = metadata,
  group = "Disease",
  covariates = c("age", "sex"),  # Adjust for confounders
  method = "camera",
  pathway_type = "KEGG"
)
View the results
head(gsea_results)

visualize_gsea()

The visualize_gsea() function creates various visualizations for GSEA results, including enrichment plots, dot plots, network plots, and heatmaps.

r
library(ggpicrust2)
library(tidyverse)
Load example data and perform GSEA
data("ko_abundance")
data("metadata")
gsearesults <- pathwaygsea(   abundance = koabundance %>% columnto_rownames("#NAME"),   metadata = metadata,   group = "Environment" )
Create an enrichment plot for a specific pathway
Annotate results so plots can use readable pathway names
annotatedresults <- gseapathway_annotation(
  gsearesults = gsearesults,
  pathway_type = "KEGG"
)
Create an enrichment-style summary plot
enrichmentplot <- visualizegsea(
  gsearesults = annotatedresults,
  plottype = "enrichmentplot",
  n_pathways = 10
)
Create a dot plot showing top enriched pathways
dotplot <- visualizegsea(
  gsearesults = annotatedresults,
  plot_type = "dotplot",
  n_pathways = 20,  # Show top 20 pathways
  sort_by = "NES"   # Sort by Normalized Enrichment Score
)
Create a network plot showing pathway relationships
networkplot <- visualizegsea(
  gsearesults = gsearesults,
  plot_type = "network",
  n_pathways = 15,
  network_params = list(
    similarity_measure = "jaccard",
    similarity_cutoff = 0.2,
    layout = "fruchterman",
    nodecolorby = "NES"
  )
)
Create a heatmap showing pathway gene expression
heatmapplot <- visualizegsea(
  gsearesults = gsearesults,
  plot_type = "heatmap",
  abundance = koabundance %>% columnto_rownames("#NAME"),
  metadata = metadata,
  group = "Environment",
  n_pathways = 10,
  heatmap_params = list(
    cluster_rows = TRUE,
    cluster_columns = TRUE,
    showcolumnnames = TRUE,
    showrownames = FALSE
  )
)

comparegseadaa()

The comparegseadaa() function compares results from GSEA and differential abundance analysis (DAA) to identify pathways that are consistently identified by both methods or uniquely identified by each method.

r
library(ggpicrust2)
library(tidyverse)
Load example data
data("ko_abundance")
data("metadata")
Prepare pathway-level abundance for DAA so identifiers match GSEA pathway IDs
keggpathwayabundance <- ko2keggabundance(data = koabundance)
Perform GSEA analysis
gsearesults <- pathwaygsea(
  abundance = koabundance %>% columnto_rownames("#NAME"),
  metadata = metadata,
  group = "Environment"
)
Perform DAA analysis
daaresults <- pathwaydaa(
  abundance = keggpathwayabundance,
  metadata = metadata,
  group = "Environment",
  daa_method = "ALDEx2"
)
Compare GSEA and DAA results
comparison <- comparegseadaa(
  gsearesults = gsearesults,
  daaresults = daaresults,
  p_threshold = 0.05,
  plot_type = "venn"  # Can be "venn", "upset", or "scatter"
)
View the comparison plot
comparison$plot
View the overlapping pathways
head(comparison$results$overlap)

gseapathwayannotation()

The gseapathwayannotation() function annotates GSEA results with pathway information, including pathway names, descriptions, and classifications.

r
library(ggpicrust2)
library(tidyverse)
Load example data and perform GSEA
data("ko_abundance")
data("metadata")
gsearesults <- pathwaygsea(   abundance = koabundance %>% columnto_rownames("#NAME"),   metadata = metadata,   group = "Environment" )
Annotate GSEA results
annotatedresults <- gseapathway_annotation(
  gsearesults = gsearesults,
  pathway_type = "KEGG"
)
View the annotated results
head(annotated_results)

FAQ

Issue 1: pathway_errorbar error

When using pathway_errorbar with the following parameters:

r
pathway_errorbar(abundance = abundance,
                 daaresultsdf = daaresultsdf,
                 Group = metadata$Environment,
                 kotokegg = TRUE,
                 pvaluesthreshold = 0.05,
                 order = "pathway_class",
                 select = NULL,
                 pvaluebar = TRUE,
                 colors = NULL,
                 xlab = "pathwayname")

You may encounter an error:

Error in ggplot_add(): ! Can't add e2 to a object. Run rlang::last_trace() to see where the error occurred.

Make sure you have the patchwork package loaded:

r
library(patchwork)

Issue 2: guidetrain.prismoffset_minor error

You may encounter an error with guidetrain.prismoffset_minor:

Error in guidetrain.prismoffsetminor(guide, panelparams[[aesthetic]]) : No minor breaks exist, guideprismoffset_minor needs minor breaks to work

Error in get(as.character(FUN)，mode = "function"object envir = envir) guideprismoffset_minor' of mode'function' was not found

Ensure that the ggprism package is loaded:

r
library(ggprism)

Issue 3: SSL certificate problem

When encountering the following error:

SSL peer certificate or SSH remote key was not OK: [rest.kegg.jp] SSL certificate problem: certificate has expired

If you are in China, make sure your computer network can bypass the firewall.

Issue 4: Bad Request (HTTP 400)

When encountering the following error:

Error in .getUrl(url, .flatFileParser) : Bad Request (HTTP 400).

Please restart R session.

Issue 5: Error in grid.Call(C_textBounds, as.graphicsAnnot(xlabel),x$x, x$y, :

When encountering the following error:

Error in grid.Call(C_textBounds, as.graphicsAnnot(xlabel),x$x, x$y, :

Please having some required fonts installed. You can refer to this thread.

Issue 6: Visualization becomes cluttered when there are more than 30 features of statistical significance.

When faced with this issue, consider the following solutions:

Solution 1: Utilize the ‘select’ parameter

The ‘select’ parameter allows you to specify which features you wish to visualize. Here’s an example of how you can apply this in your code:

ggpicrust2::pathway_errorbar( abundance = kegg_abundance, daaresultsdf = daaresultsdf_annotated, Group = metadata$Day, pvaluesthreshold = 0.05, order = "pathway_class", select = c("ko05340", "ko00564", "ko00680", "ko00562", "ko03030", "ko00561", "ko00440", "ko00250", "ko00740", "ko04940", "ko00010", "ko00195", "ko00760", "ko00920", "ko00311", "ko00310", "ko04146", "ko00600", "ko04141", "ko04142", "ko00604", "ko04260", "ko00909", "ko04973", "ko00510", "ko04974"), kotokegg = TRUE, pvaluebar = FALSE, colors = NULL, xlab = "pathwayname" )

Solution 2: Limit to the Top 20 features

If there are too many significant features to visualize effectively, you might consider limiting your visualization to the top 20 features with the smallest adjusted p-values:

daaresultsdfannotated <- daaresultsdfannotated[!is.na(daaresultsdfannotated$pathwayname),]

daaresultsdfannotated$padjust <- round(daaresultsdfannotated$padjust,5)

lowpfeature <- daaresultsdfannotated[order(daaresultsdfannotated$p_adjust), ]$feature[1:20]

p <- ggpicrust2::pathway_errorbar( abundance = kegg_abundance, daaresultsdf = daaresultsdf_annotated, Group = metadata$Day, pvaluesthreshold = 0.05, order = "pathway_class", select = lowpfeature, kotokegg = TRUE, pvaluebar = FALSE, colors = NULL, xlab = "pathwayname")

Issue 7: There are no statistically significant biomarkers

If you are not finding any statistically significant biomarkers in your analysis, there could be several reasons for this:

• **The true difference between your groups is small or

non-existent.** If the microbial communities or pathways you’re comparing are truly similar, then it’s correct and expected that you won’t find significant differences.

• Your sample size might be too small to detect the differences.

Statistical power, the ability to detect differences if they exist, increases with sample size.

• The variation within your groups might be too large. If there’s

a lot of variation in microbial communities within a single group, it can be hard to detect differences between groups.

Here are a few suggestions:

• Increase your sample size: If possible, adding more samples to

your analysis can increase your statistical power, making it easier to detect significant differences.

• Decrease intra-group variation: If there’s a lot of variation

within your groups, consider whether there are outliers or subgroups that are driving this variation. You might need to clean your data, or to stratify your analysis to account for these subgroups.

• Change your statistical method or adjust parameters: Depending

on the nature of your data and your specific question, different statistical methods might be more or less powerful. If you’re currently using a parametric test, consider using a non-parametric test, or vice versa. Also, consider whether adjusting the parameters of your current test might help.

Remember, not finding significant results is also a result and can be informative, as it might indicate that there are no substantial differences between the groups you’re studying. It’s important to interpret your results in the context of your specific study and not to force statistical significance where there isn’t any.

With these strategies, you should be able to create a more readable and informative visualization, even when dealing with a large number of significant features.

Author’s Other Projects

• MicrobiomeStat: The

MicrobiomeStat package is a dedicated R tool for exploring longitudinal microbiome data. It also accommodates multi-omics data and cross-sectional studies, valuing the collective efforts within the community. This tool aims to support researchers through their extensive biological inquiries over time, with a spirit of gratitude towards the community’s existing resources and a collaborative ethos for furthering microbiome research.

If you’re interested in helping to test and develop MicrobiomeStat, please contact .

• MicrobiomeGallery:

This is a web-based platform currently under development, which aims to provide a space for sharing microbiome data visualization code and datasets.

We look forward to sharing more updates as these projects progress.