🌟 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 analyzes PICRUSt2 predicted functional profiles with pathway aggregation, differential abundance analysis, gene-set testing, annotation, and visualization. It also summarizes taxon contributions and compares results across analysis methods. Method agreement is a sensitivity analysis on the same data; it does not validate measured pathway activity.
See release notes for release notes and fixes. The tutorials below cover DAA, camera/fry and preranked gene-set tests, and taxon contribution analysis.
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:
citation("ggpicrust2")You can install the development version of ggpicrust2 from GitHub with:
# install.packages("devtools")
devtools::install_github("cafferychen777/ggpicrust2")Required dependencies are installed with the package. Optional analysis backends must be installed separately for the methods you use:
| Workflow | Optional packages used in the tutorial |
|---|---|
| LinDA and KEGG annotation | MicrobiomeStat, KEGGREST |
| ALDEx2 alternative | ALDEx2 |
| camera/fry GSEA | limma |
| Preranked GSEA and its displays | fgsea, ggridges,
ComplexHeatmap, circlize,
igraph |
| GSEA/DAA Venn plot | ggVennDiagram |
The general workflow and GSEA tutorial include installation commands for their examples. The full dependency declarations are in DESCRIPTION.
Use the project resources below for stable updates and support:
Start with the bundled data. No external example files are needed.
Install MicrobiomeStat for LinDA and KEGGREST
for online KEGG annotation as shown in Using
ggpicrust2.
library(ggpicrust2)
data("ko_abundance")
data("metadata")
results <- ggpicrust2(
data = ko_abundance,
metadata = metadata,
group = "Environment",
pathway = "KO",
daa_method = "LinDA",
ko_to_kegg = TRUE,
x_lab = "pathway_name"
)
results[[1]]$plot
head(results[[1]]$results)A plot can be NULL when there are no significant,
annotated pathways. For your own PICRUSt2 files, replace
data = ko_abundance with
file = "your_abundance.tsv" and supply metadata with
matching sample IDs.
Using ggpicrust2 is the canonical complete walkthrough. It uses the same LinDA settings as the example above and covers sample alignment, annotation, plotting, ALDEx2 test selection, and contribution analysis. The GSEA tutorial covers gene-set tests, contrasts, covariates, and method-specific visualization.
Each numbered result contains a method’s plot and
annotated results table. The wrapper also returns aligned
abundance, metadata, group, and
daa_results_df for further analysis. These remain
predictions from PICRUSt2; statistical differences do not establish
measured pathway activity or causation.
data("ko_abundance")
kegg_abundance <- ko2kegg_abundance(data = ko_abundance)The default is the mean of the upper half of observed KO abundances
in each pathway/sample, not a sum or a complete reconstruction of
pathway activity. Overlapping pathways share KOs.
method = "sum" is a different aggregation and changes the
analysis scale. See ?ko2kegg_abundance for the mapping and
prokaryote-filtering rules. A retained KEGG label is not proof of that
pathway’s activity or organism specificity; shared KOs can map to
several biological contexts.
data("metacyc_abundance")
data("metadata")
metacyc_matrix <- tibble::column_to_rownames(metacyc_abundance, "pathway")
metacyc_results <- pathway_daa(
abundance = metacyc_matrix, metadata = metadata,
group = "Environment", daa_method = "LinDA"
)Features are rows and samples are columns. Choose a method for the input scale and design. Count-based methods may round non-integer predicted abundances; that is not an equivalent treatment of normalized or continuous data. For multiple groups, specify the reference and inspect the returned group pairs.
This compares discoveries at an adjusted p-value threshold. ALDEx2
returns multiple tests; split its results by method before
passing one table per test. The function treats each list element as one
discovery set. Use ?compare_daa_results for a runnable
example. Agreement is sensitivity analysis on the same data, not
independent biological validation.
annotated_metacyc <- pathway_annotation(
pathway = "MetaCyc", daa_results_df = metacyc_results,
ko_to_kegg = FALSE
)Annotation adds labels without converting the abundance or changing
the unit of a p-value. To test KEGG pathways, aggregate KO abundance
before DAA. ko_to_kegg = TRUE then selects online KEGG
annotation for those result rows; only rows below
p_adjust_threshold are annotated. Local KO, EC, MetaCyc,
and KEGG reference annotation is also available; see
?pathway_annotation.
Use the stepwise
tutorial for a complete example. Group is a vector of
labels, unlike group, which is a metadata column name.
setNames(metadata$Environment, metadata$sample_name)
preserves sample identity when the metadata and abundance table have
different orders.
The abundance panel summarizes per-sample relative abundances using
the full input feature matrix as the denominator. Use
select to limit displayed rows; pre-filtering the abundance
matrix changes the denominator. A supplied model
log2_fold_change is preserved; otherwise the plot derives a
descriptive ratio of group mean relative abundances. These are not
necessarily the same quantity. See ?pathway_errorbar and
?pathway_errorbar_table.
sig_features <- unique(metacyc_results$feature[
!is.na(metacyc_results$p_adjust) & metacyc_results$p_adjust < 0.05
])
if (length(sig_features) > 0) {
pathway_heatmap(metacyc_matrix[sig_features, , drop = FALSE],
metadata, group = "Environment")
}Keep stable feature IDs as row names; descriptions can be duplicated or missing. Colors are within-feature z-scores across samples. They do not compare absolute abundance between different pathways or establish statistical significance.
set.seed(123)
abundance_example <- matrix(rexp(30), nrow = 3, ncol = 10,
dimnames = list(c("PathwayA", "PathwayB", "PathwayC"), paste0("Sample", 1:10)))
metadata_example <- data.frame(
sample_name = colnames(abundance_example),
group = factor(rep(c("Control", "Treatment"), each = 5))
)
pathway_pca(abundance_example, metadata_example, group = "group")PCA centers and scales features. It does not normalize sample library sizes. This simulated example illustrates the display, not a biological group effect. Use preprocessing appropriate for the study; visible separation is not a test.
Use this for matrices measuring the same samples, with shared sample
IDs. It retains the pairing in DAA and reports median per-feature
Spearman correlations with joint sample-label permutation p-values. See
?compare_metagenome_results for a self-contained example
and the exact feature-intersection rules. Correlation is not proof of
unbiased agreement.
The contribution
tutorial includes synthetic inputs and examples for real PICRUSt2
files. read_contrib_file(),
read_pathway_contrib_file(), and
read_strat_file() standardize their formats;
aggregate_taxa_contributions() sums contributions at the
requested rank.
For KO-level contributions, KEGG pathway filters are expanded to their KO members. Output remains KO-level contributions; this operation does not infer pathway contributions. For pathway-level MetaCyc input, matching MetaCyc IDs are filtered directly. Use metadata and taxonomy from the same dataset.
Use gene/enzyme-level abundance, not already aggregated pathways. The GSEA tutorial explains the input assumptions:
camera tests a competitive null with an inter-gene
correlation adjustment.fry tests a self-contained null; it is not an
interchangeable faster camera.fgsea and clusterProfiler use a ranked
feature list and a different null.The camera/fry implementation uses voom on count-like abundance. A successful fit is not evidence that its uncertainty model is calibrated for every kind of predicted, normalized, or relative abundance.
See the GSEA
visualization examples. For camera/fry, the compatibility column
NES contains a signed -log10(pvalue) score,
not a normalized enrichment score or effect size.
enrichment_plot is a score-summary bar chart, not a running
enrichment curve.
Network edges and the GSEA heatmap require leading-edge genes from a preranked method. camera/fry do not supply these genes. A leading-edge heatmap displays row-standardized mean supplied abundance per leading-edge set, not one row per gene or measured gene expression.
Compare matching pathway identifiers and compatible contrasts. The GSEA tutorial shows KEGG pathway DAA versus KEGG GSEA. Gene-set and pathway-abundance tests answer different questions, so neither agreement nor disagreement alone establishes that a method is correct.
gsea_pathway_annotation() adds pathway names and
classes. It preserves the method-specific scores and p-values;
annotation cannot supply missing leading-edge genes or change the
analysis level.
Group = "Environment" passes a single string, not the
samples’ group labels. Use
Group = setNames(metadata$Environment, metadata$sample_name)
and replace those column names for your data. The sample names must
cover all abundance columns. An unnamed vector is interpreted in
abundance-column order; simply using metadata$Environment
can silently mislabel samples when the two tables have different orders,
including in the bundled example data.
The stepwise
tutorial shows the complete workflow and checks. For errors
constructing a plot, include sessionInfo() and the full
error in an issue; loading additional plotting libraries does not
resolve a group-length or sample-alignment error.
For plotting errors, report the full error and
sessionInfo(). Loading extra plotting packages does not fix
mismatched data, unsupported arguments, or incompatible dependency
versions.
For KEGG HTTP errors, inspect the reported feature ID and confirm the annotation mode matches the identifier level. A bad request is not fixed by restarting R. For certificate failures, check the system clock, certificate configuration, and network access. Local reference annotation is available when live KEGG access is unnecessary.
Use select to display a subset of significant feature
IDs, or increase max_features deliberately. For the objects
created by the stepwise tutorial:
significant_results <- annotated_daa[
!is.na(annotated_daa$p_adjust) & annotated_daa$p_adjust < alpha, , drop = FALSE
]
top_features <- head(
significant_results$feature[order(significant_results$p_adjust)], 20
)
if (length(top_features) > 0) {
pathway_errorbar(
abundance = kegg_pathway_abundance,
daa_results_df = annotated_daa,
Group = sample_groups,
select = top_features,
ko_to_kegg = TRUE,
p_values_threshold = alpha,
x_lab = "pathway_name"
)
}This only limits the display; it does not rerun testing, round
p-values, or change the significance threshold. head() also
works when fewer than 20 pathways are significant.
An empty set at the chosen adjusted p-value threshold is a valid result. Retain the results table and report the method, threshold, and sample size. Check input quality and model assumptions; do not relax thresholds, remove samples, or switch tests solely to obtain a plot. Prespecify the primary analysis and label sensitivity analyses separately.
If you’re interested in helping to test and develop MicrobiomeStat, please contact cafferychen7850@gmail.com.
We look forward to sharing more updates as these projects progress.