ggpicrust2 2.5.19
Analysis options
-
pathway_gsea()now accepts explicitgene_setsand an opt-inlogCPMtransformation for camera/fry. The latter tests log relative abundances with variance moderation and is invariant to positive per-sample rescaling. The existing voom default is unchanged. Neither transformation guarantees calibration for predicted functional abundances or removes their upstream prediction uncertainty. -
pathway_daa()exposes LinDA’s winsorization, adaptive zero-handling flag, and fixed pseudo-count without changing the historical defaults. Disabling winsorization preserves fractional abundances that the count backend would otherwise round. LinDA output now includes its native standard errors, statistics, degrees of freedom, and nominal pointwise 95% t intervals.
Documentation and diagnostics
Fixed the existing CRAN check findings: restored generated usage sections for three plotting/comparison helpers and made backend-dependent tests respect unavailable optional packages, including ALDEx2 on macOS oldrel. Help examples also check for their optional analysis/plotting dependencies.
Restored the installed-package test entry point so
R CMD checkactually executes the bundled regression tests.Camera tutorial examples now estimate within-set correlation explicitly. A fixed correlation of 0.01 can underestimate variance for strongly correlated features; changing the abundance transformation alone does not address this assumption. The public function default is unchanged.
Restored the documentation website build and deployment from maintained sources, with generated HTML kept out of the source branch.
GSEA/DAA scatter plots now use the method-specific score label instead of labeling camera/fry signed p-value scores as normalized enrichment scores.
Audited all help examples, README workflows, and both vignettes against executable examples and underlying data transformations. Examples now use available metadata columns and self-contained inputs; README points to canonical walkthroughs instead of duplicating drifting workflows.
Clarified method-specific GSEA scores, contrasts, null hypotheses, leading-edge plots, shared tested universes, and the descriptive transformations in PCA, heatmaps, and ridge plots. Annotation adds labels and does not aggregate abundance or convert KO-level p-values to pathway-level evidence.
Corrected the contribution tutorial: KEGG filters can select member KOs, but the returned contribution table remains KO-level. Added a synthetic example with consistent taxonomy and sample IDs.
Fixed the stepwise tutorial’s
pathway_errorbar()call (#207): supply group labels named by sample ID, set KEGG pathway annotation parameters, and skip significance plots when no pathways pass the threshold. Both main workflows now use LinDA; the ALDEx2 alternative selects one test.Consolidated duplicated README workflows into the main vignette, corrected unsafe group-vector examples in plotting/table help, and removed FAQ code that rounded p-values or indexed past the available significant features.
Corrected the GSEA tutorial’s KO-versus-pathway comparison and clarified identifier matching for optional taxa-contribution filtering.
Group-length errors now report the label/sample counts and explain how to supply sample-ID-named labels instead of a metadata column name.
ggpicrust2 2.5.18
Bug Fixes
Color themes now reject unknown names instead of silently falling back to the default theme, expand large palettes without recycling identical colors, and keep accessibility mode authoritative during smart theme selection.
visualize_gsea()now validates named network/heatmap parameter bundles, applies custom scales consistently across supported plot types, rejects incompatible heatmap scale combinations, preserves group-to-color identity, and disambiguates duplicate pathway display labels with pathway IDs.pathway_ridgeplot()now requires an explicit pathway-reference ID schema, parses member lists consistently across supported table layouts, rejects ambiguous member columns, and disambiguates duplicate pathway labels.pathway_heatmap()now validates secondary groups before concatenation, rejects conflicting deprecated/new grouping arguments, avoids duplicate correlation-distance work, and keeps explicit legend/row-name visibility flags authoritative over custom themes.pathway_pca()now maps named colors by group identity and omits marginal density estimates for singleton groups instead of constructing invalid or unused density layers.pathway_volcano()now requires an available label column when labels are requested, removes blank labels, maps named colors by significance class, and preserves positive subnormal p-values on the negative-log10 scale.compare_gsea_daa()now rejects missing scatter probabilities and partial group-direction schemas, preserves positive subnormal probabilities, and returns the same explicit count plot for empty Venn and UpSet universes.Deprecated
p.adjustaliases now share one compatibility path across the top-level, DAA, GSEA, and metagenome-comparison entry points and can no longer silently override a conflicting explicitp_adjust_method.Taxa-contribution plots now preserve requested function order, reject unsupported annotation schemas instead of silently retaining raw IDs, and normalize annotation-label whitespace before conflict detection.
compare_daa_results()now rejects invalid self-comparisons wheregroup1 == group2and constructs its result table without repeated row-binding.pathway_heatmap()andtaxa_contribution_heatmap()now validate hierarchical clustering method/distance combinations. Ward linkage methods ("ward.D"and"ward.D2") now require Euclidean distance because Ward clustering has a within-cluster variance interpretation in Euclidean space; correlation or rank-based distances should be paired with linkages such as"average"or"complete".Count-like controls such as
n_pathways,top_n, andcorrelation_permutationsnow reject values larger than R’s integer range without emitting coercion warnings.visualize_gsea()also applies the same no-coercion integer check to GSEAsizevalues before plotting.taxa_contribution_bar(show_percentage = TRUE)now treats absent sample/function rows as zero-total combinations when validating percentage denominators, instead of silently dropping those sample/function bars from the plot. Explicitfunction_idsare also checked after sample alignment so partially missing requests fail with the missing IDs instead of being silently reduced to the matching subset.pathway_annotation()now normalizes logical-likeko_to_keggstrings, preventingko_to_kegg = "TRUE"from silently taking the local-reference branch instead of the requested KEGG annotation branch.pathway_gsea(method = "camera"|"fry")now rejects non-finite or rank-deficient design matrices before limma is called. Constant covariates or covariates perfectly confounded with the group variable now fail with an actionable error instead of entering an unestimable camera/fry contrast.pathway_gsea(method = "camera"|"fry")now builds limma design formulas from literal metadata column names, so non-syntactic group or covariate names such as"treatment group"and"age years"work without requiring users to rename their metadata.Fixed R scoping bug in
pathway_annotation(): error counting insidetryCatch()error handler used<-(local assignment) instead of<<-, soerror_countanderror_idswere silently never updated when KEGG API calls failed.Fixed row-specific KEGG annotation merge-back in
pathway_annotation()when the same feature appears in multiple DAA rows. Non-significant rows with the same feature ID as a significant row now keepNAannotation fields instead of inheriting annotations by feature-name matching.ko2kegg_abundance()now rejects duplicated KO identifiers after cleaning optionalko:prefixes. Duplicate KO rows previously entered pathway aggregation as repeated evidence and could distort both upper-half mean and legacy sum abundances.ko2kegg_abundance()now rejects missing or non-finite KO abundance values before pathway aggregation. The previous upper-half mean path could silently drop missing KO values during sorting, changing the number of KOs contributing to a pathway abundance.visualize_gsea()now rejects missing, empty, or duplicatepathway_idvalues in the selected rows for every plot type, preventing low-level duplicate row-name errors in heatmaps and ambiguous multi-method/contrast visualizations. Missing or empty pathway labels now fall back topathway_idrow-wise.pathway_ridgeplot()now calculates fold changes using explicit group-level semantics. Two-group inputs use factor-level order, while multi-group inputs must supplycomparison = c(group1, group2), preventing sample-column order from silently changing the log2 fold-change direction or plotting a comparison that does not match the GSEA contrast.-
Fixed dead auto-adjust logic in
create_legend_theme():missing(direction)was checked aftermatch.arg(direction, ...)which always evaluates theformal, so the auto-adjust based on legend position never executed.
Aligned
run_fgsea()internalmin_sizedefault (was 10) with the publicpathway_gsea()API default of 5.Aligned
perform_aldex2_analysis()internalinclude_effect_sizedefault (wasFALSE) with the publicpathway_daa()API default ofTRUE.Two-group ALDEx2 analyses now stop if explicitly requested
aldex.effect()output fails or is malformed, instead of warning and returning p-value-only rows without the documented effect-size columns. Effect rows are matched to test rows by feature identifier and the required finiteeffect,diff.btw,rab.all, and probability-valuedoverlapcolumns are validated before use.Paired ALDEx2 results from
compare_metagenome_results()now use the same feature-ID alignment and effect-output validation. Previously, a reorderedaldex.effect()table could attach one feature’s effect size and fold change to another feature’s p-values by row position.Fixed
run_fgsea()empty-result path: missing method argument tocreate_empty_gsea_result()producedmethod = "unknown"instead of"fgsea".Collapsed dead contrast-selection branches in
run_limma_gsea()where both if/else paths assigned the same value; the unreachable finalelsenow raises an informative error for unexpectedcontrasttypes.Fixed
pathway_gsea(method = "camera"|"fry")contrast resolution: multi-group designs now require an explicit contrast, character contrasts must exactly match a design column or non-reference group level, and numeric contrasts are validated against the design matrix. This prevents substring matching from silently testing the wrong coefficient. Named numeric contrast vectors are now matched to design column names before being passed to limma, preventing out-of-order named vectors from silently testing the wrong coefficient combination.Fixed
pathway_daa(daa_method = "ALDEx2", reference = ...): ALDEx2 condition levels are now releveled before converting them to numeric conditions, sodiff_btw/log2_fold_changedirection andgroup1/group2labels honor the requested reference level.Fixed
pathway_daa(daa_method = "Lefser", reference = ...): Lefser now relevels the grouping factor before analysis, converts counts to lefser’s expected relative-abundance assay withlefser::relativeAb(), and reports all-feature Kruskal-Wallis p-values from that same relative-abundance scale. Previously raw-count library-size differences could drive the reported p-values and the user-supplied reference was ignored. Actual Lefser backend errors now fail the call instead of returning a table with missing LDA scores.The top-level
ggpicrust2()wrapper no longer rejectsdaa_method = "Lefser"with the stale claim that Lefser does not output p-values.pathway_daa()now owns Lefser validation and supplies p-values/adjusted p-values, whilepathway_errorbar()supplies the display-only log2 fold-change fallback when no method-native log2 fold change is returned.pathway_daa(daa_method = "Lefser")now fails clearly if a per-feature Kruskal-Wallis p-value cannot be computed on the Lefser relative-abundance scale, instead of silently converting the failed test top = 1.pathway_daa(daa_method = "LinDA")now surfaces backend errors as errors instead of returning an empty DAA table, so method failures cannot be mistaken for “no differential features”.pathway_daa()now validatesreferenceonce after sample alignment andselectfiltering. Invalid, empty, or filtered-out reference levels now fail with an actionable error instead of silently falling back to the first observed group level.pathway_daa(select = ...)now requires unique, non-empty sample names and rechecks that at least four samples remain after filtering, preventing duplicated or undersized selected sample sets from reaching backend fitting.pathway_daa()now rejects missing/empty group labels and requires at least two samples per retained group after sample alignment andselectfiltering. This avoids fitting DAA models or tests on singleton groups where within-group variation cannot be assessed.pathway_daa()now rejects missing, non-finite, or negative abundance values through shared DAA input validation before backend fitting. Invalid abundance cells no longer reach different DAA backends with backend-specific low-level failures or partial statistics.pathway_daa()now preserves method-native adjusted p-values for DESeq2 (padjfromresults()), LinDA (padj), and Maaslin2 (qval) instead of discarding them and recomputing a generic wrapper-level adjustment. The requestedp_adjust_methodis forwarded to those backends where supported.pathway_daa(daa_method = "DESeq2")no longer suppresses all backend warnings. DESeq2 diagnostics such as identical values across every sample now reach the caller instead of being hidden while a result table is returned.pathway_daa(daa_method = "DESeq2"|"limma voom"|"edgeR")now validates backend feature identifiers and aligns p-values, adjusted p-values, and log fold changes by feature ID before constructing the result table. A reordered, incomplete, or malformed backend result now fails clearly instead of attaching statistics to features by row position.pathway_daa(daa_method = "metagenomeSeq")no longer silently replaces all CSS normalization-quantile failures withp = 0.5. Samples with only one positive feature now fail with an actionable error, matching metagenomeSeq’s owncumNormStatFast()requirement, while the existing degenerate-search fallback top = 0.5now emits a warning.pathway_daa(daa_method = "Maaslin2")now mirrors MaAsLin2’smake.names()feature-name sanitization when mapping results back to original feature IDs, rejects ambiguous sanitized feature IDs before model fitting, and validates MaAsLin2 result columns before returning them.pathway_daa()now validates raw and adjusted p-value columns at the DAA wrapper boundary, including method-specific adjusted p-values supplied by individual backends, so malformed probability columns fail instead of propagating into downstream summaries.pathway_daa(daa_method = "DESeq2"|"metagenomeSeq"|"LinDA"|"Maaslin2")now models against an internal syntactic metadata column, so valid user-supplied group columns with names such as"treatment group"work without forcing users to rename their metadata. Outputgroup1andgroup2labels continue to use the original group levels.pathway_daa(daa_method = "LinDA")now validates the documented MicrobiomeStat::linda output columns consumed by ggpicrust2 (pvalue,padj, andlog2FoldChange) and fails on missing, malformed, non-finite, or out-of-range values instead of filling them withNA.Wrapper-computed DAA adjusted p-values are now calculated within each method/contrast (
method,group1,group2) instead of across all pairwise contrasts at once, matching the per-contrast result semantics of DESeq2, limma, and edgeR-style outputs.Count-based DAA backends that require or assume integer counts (ALDEx2, DESeq2, edgeR, and metagenomeSeq) now warn when non-integer abundance values are rounded before fitting, replacing previously silent data mutation.
pathway_daa()now requires explicit, non-empty, unique feature identifiers in abundance row names (or a leading feature-ID column that is normalized to row names). Missing IDs fail before backend fitting, and duplicated IDs no longer propagate into ambiguous DAA result rows or downstream annotation merges.pathway_daa()now rejects additional arguments supplied through...instead of silently ignoring them. The previous documentation claimed these arguments were passed to backend DAA methods, which could mislead users into believing formula, fixed-effect, or covariate-adjustment arguments had been applied when the fitted model was actually unchanged.pathway_gsea()now validatesmin_size,max_size,nperm,seed,p_adjust_method, andinter.gene.corat the API boundary, including rejecting impossible gene-set size ranges before downstream limma/fgsea/ clusterProfiler calls.pathway_gsea()now revalidates group structure and complete design variables after sample alignment. Missing group labels or camera/fry covariates now fail before ranking/model fitting, preventingmodel.matrix()from dropping samples and preventing preranked methods from silently using a different sample set. This validation is also applied when preranked methods use an explicitcomparison = c(group1, group2), so missing aligned group labels are not silently dropped from the ranked list.pathway_gsea(method = "camera"|"fry")now passes raw non-negative counts directly tolimma::voom(). The previous pre-voom+0.5pseudocount changed library sizes and could distort voom’s mean-variance trend. A failed voom transformation now stops the analysis instead of silently switching to a log2 transform with unit weights, which discarded voom’s observation-level mean-variance weights while still reporting camera/fry results.Camera/fry GSEA output now explicitly labels its legacy
NEScompatibility column asscore_type = "signed_log10_pvalue"withscore_label = "Signed -log10(p-value)". limma camera/fry do not estimate a true normalized enrichment score, andvisualize_gsea()now uses the explicit score label for axes and legends to avoid misinterpretation.pathway_gsea(method = "fgsea"|"GSEA"|"clusterProfiler")now acceptscomparison = c(group1, group2)to define preranked GSEA ranking direction explicitly. Positive ranking statistics and positive ES/NES now map to features higher ingroup1; multi-group preranked analyses must specifycomparisoninstead of relying on implicit factor-level order.pathway_gsea()now rejects design-only arguments that cannot be honored by preranked GSEA methods. Supplyingcovariatestofgsea,GSEA, orclusterProfiler, or supplyingcontrastto a preranked method, now fails with instructions to use limma-basedcamera/fryor the prerankedcomparisonargument instead of silently running an unadjusted analysis.pathway_gsea(method = "fgsea")now honors the requestedp_adjust_methodby recomputing adjusted p-values from fgsea raw p-values. Previously the fgsea branch always returned fgsea’s built-in BH-adjustedpadjvalues even when callers requested another correction.Preranked GSEA ranking vectors are now validated before fgsea or clusterProfiler execution. Ranking statistics must be a named numeric vector with unique feature names, finite non-missing values, and at least two distinct statistics. This prevents all-tied rankings (for example all-zero group differences) from producing enrichment results driven by arbitrary input order rather than biological signal.
pathway_gsea(method = "GSEA"|"clusterProfiler")now returns a standard empty GSEA result schema when no gene sets overlap the ranked feature list after size filtering, instead of returning a partial data frame missing columns such asleading_edgeandmethod.pathway_gsea()now requires a numeric abundance matrix with explicit, non-duplicated feature row names (or a leading non-numeric feature-ID column) before gene-set matching, so malformed input fails at the API boundary instead of producing empty or low-level GSEA errors.pathway_gsea()and its internal preranked/limma GSEA helpers now reject negative, missing, or non-finite count-like abundance values instead of replacing them with zero, preventing silent changes to ranking statistics and voom’s mean-variance model.pathway_gsea(pathway_type = "MetaCyc")now rejects pathway-level MetaCyc identifiers in the abundance row names. MetaCyc GSEA requires EC-level input for EC-to-pathway gene sets; pathway-level MetaCyc abundance tables should be analyzed withpathway_daa()instead of being treated as empty enrichment evidence.pathway_gsea(),prepare_gene_sets(),run_fgsea(), andrun_limma_gsea()now validate scalar method/pathway/ranking choices and named gene-set lists before enrichment testing. Missing, empty, or duplicated gene-set names are rejected because those names becomepathway_idvalues; this prevents downstream R/limma name repair from silently changing pathway identifiers such asset1intoset1.1.Shared choice-parameter validation now requires a single non-missing supported value. This gives clear errors for invalid
plot_type,sort_by,order, method, and pathway choices instead of low-level R condition-length errors when callers pass vectors orNA.Added shared validation for
p_adjust_methodacross DAA/GSEA-facing entry points using R’sstats::p.adjust.methods, so unsupported adjustment names fail even when a backend returns method-specific adjusted p-values.Shared abundance validation now requires numeric matrix input and numeric sample columns in data frames, while treating a leading non-numeric feature ID column as metadata rather than a sample. This prevents malformed abundance tables from reaching low-level
colSums(),rowMeans(),scale(), or backend model-fitting calls with misleading errors.Abundance-consuming entry points now normalize a leading non-numeric feature ID column into row names before sample alignment and matrix conversion. This prevents feature/pathway IDs from being dropped or replaced by default row numbers in
pathway_daa(),pathway_gsea(), heatmaps, PCA, ridge plots, and error-bar summaries.ggpicrust2(data = ..., ko_to_kegg = FALSE)now preserves abundance inputs that already store feature IDs in row names. The wrapper no longer unconditionally treats the first sample column as an ID column, preventing the first sample from being dropped and abundance values from becoming feature labels.pathway_heatmap()now handles zero-variance rows or sample profiles when row/column clustering uses correlation or Spearman distance, instead of sendingNAdistances tohclust().pathway_heatmap()now accepts finite transformed inputs with negative values or zero-sum sample columns, while explicitly rejecting missing or non-finite values before z-score scaling.pathway_heatmap()now revalidates primary and secondary grouping variables after sample alignment and rejects missing aligned group labels, preventing extra metadata rows from making an otherwise invalid heatmap grouping appear valid.ko2kegg_abundance(method = "abundance")now matches PICRUSt2’s unstructured pathway upper-half indexing for even numbers of matched KOs (floor(n / 2) + 1in R). The previousceiling(n / 2)start included one lower-half KO for even pathway sizes and could underestimate pathway abundance.Clarified
ko2kegg_abundance()documentation: the default method is an offline KO-to-KEGG aggregation approximation based on PICRUSt2’s unstructured pathway rule, not the full PICRUSt2 pathway pipeline with MinPath and structured MetaCyc pathway inference.Fixed
pathway_errorbar()documentation defaults forpvalue_format(was “smart”, code uses “numeric”) andpathway_class_position(was “left”, code uses “right”).Fixed significance star/color assignment in
get_significance_stars()andget_significance_colors(): iteration order now goes from least to most significant threshold so the tightest match wins.pathway_pca()now correctly requires at least 2 groups (min_groups = 2) instead of allowing single-group input that produces a degenerate PCA.pathway_pca()now rejects missing or empty group labels after sample alignment. Previously a metadata table could retain ungrouped samples in the PCA coordinates as long as the non-missing labels still contained two groups, making color and marginal-density interpretation incomplete.pathway_pca()now draws confidence ellipses only for groups with at least four samples, matching ggplot2’sstat_ellipse()implementation. Smaller groups remain in the PCA scatter plot and trigger an explicit warning instead of relying on ggplot2’s low-level “Too few points to calculate an ellipse” message.pathway_pca()no longer rejects finite zero-sum sample columns and no longer drops samples whose values are constant across pathways. Inprcomp(t(abundance)), samples are observations and pathways are variables; only zero-variance pathways are removed before scaling.pathway_ridgeplot()now guards againstNAvalues in direction and NES columns instead of silently propagating them into ggplot aesthetics.pathway_ridgeplot()now aligns metadata to abundance columns by sample ID before calculating displayed gene/KO log2 fold changes, validates GSEA p-value/FDR/NES columns and display parameters, and rankssort_by = "NES"by absolute effect size.pathway_ridgeplot()now rejects missing, empty, or duplicatedpathway_idvalues and invalidpathway_typevalues before reference mapping. Missing or emptypathway_namelabels now fall back topathway_idinstead of causing low-level string-length errors.visualize_gsea()now validates required GSEA statistic columns before plotting, including finite NES values, probability-valued p-values/FDRs, and positive integer gene set sizes.gsea_pathway_annotation()now rejects missing or emptypathway_idvalues and invalidpathway_typevalues at the API boundary, preventing annotated GSEA tables with blank orNApathway labels.visualize_gsea(plot_type = "network"|"heatmap")now parses missing or emptyleading_edgevalues as empty sets, preventing missing leading-edge annotations from creating false pathway similarity edges.visualize_gsea(plot_type = "heatmap")now recognizes abundance input with a leading non-numeric feature-ID column, validates the aligned abundance matrix, rejects missing group annotations, and fails clearly when non-empty leading-edge genes do not match abundance row names instead of drawing an all-zero heatmap.import_MicrobiomeAnalyst_daa_results()now parses MicrobiomeAnalyst result columns by semantic names such asPvalues,FDR,Statistics, andlog2FCinstead of assuming a fixed four-column order. Feature IDs may come from a feature/name column or from row names, preventing method-specific DE outputs from being silently misread as p-values/FDR values. Imported feature IDs, p-values, FDR values, optional statistics/fold changes, method labels, and group labels are validated before returning a ggpicrust2-style DAA table.P-value annotation helpers now validate probability values, significance thresholds, star symbols, and colors.
pathway_errorbar()validates its p-value display parameters at the API boundary.pathway_volcano()now validates thatfc_thresholdis a non-negative number.pathway_volcano()andcompare_gsea_daa()now validate probability-valued p-value/FDR columns and p-value thresholds, and keep zero adjusted p-values finite in scatter/volcano log-scale displays.compare_gsea_daa(plot_type = "scatter")now rejects duplicated pathway effect-size rows before merging GSEA and DAA results, preventing many-to-many joins from inflating plotted associations. Pathway identifiers must also be non-empty.compare_gsea_daa(plot_type = "scatter")now requires explicit GSEA and DAA group-direction columns and aligns DAAlog2_fold_changevalues to the GSEA-positive NES direction before plotting. This prevents default two-group analyses from showing an artificial sign reversal between preranked GSEA (group1vsgroup2) and DAA (group2/group1) effects.compare_gsea_daa(plot_type = "venn"|"upset")now rejects direction-aware inputs that contain multiple or incompatiblegroup1/group2pairs, preventing pathways from different biological contrasts from being counted as method agreement.compare_daa_results()now compares multi-group DAA discoveries as feature/group-pair units rather than feature IDs alone, preventing agreement from being overstated when different methods flag the same feature in different pairwise contrasts. Group pairs are canonicalized as unordered for this set-level comparison, soA vs BandB vs Aare treated as the same biological comparison whileA vs BandA vs Cremain distinct. Method labels and feature/group identifiers are validated for unambiguous output.compare_metagenome_results()now drops undefined per-feature Spearman correlations before summarizing and fails with a clear error when all shared features are constant across aligned samples, instead of surfacing a low-levelwilcox.test()error.compare_metagenome_results()now preserves the paired-sample design during differential abundance analysis. The function uses paired ALDEx2 tests by default and also supports paired Wilcoxon signed-rank tests on relative abundance. Independent-group backends are rejected because treating repeated quantifications of the same samples as independent replicates produces pseudoreplication.Correlation p-values in
compare_metagenome_results()now use joint sample-label permutations of the median per-feature Spearman correlation, preserving dependence among features while breaking cross-metagenome sample correspondence. Diagonal p-values areNA, Monte Carlo p-values use the plus-one correction, and a multiplicity-adjustedp_adjust_matrixplus the contributing-feature counts are returned.compare_metagenome_results()now validates metagenome labels, feature row names, sample column names, and finite non-negative matrix values before intersecting matrices. Duplicate feature/sample IDs now fail instead of creating ambiguous by-name alignments.compare_daa_results(),pathway_errorbar(),pathway_errorbar_table(),pathway_annotation(), and the top-levelggpicrust2()wrapper now use shared validation for adjusted p-value columns and significance thresholds before filtering significant features.pathway_errorbar_table()now aligns namedGroupvectors to abundance sample columns when metadata is not provided, matchingpathway_errorbar()behavior and preventing order-dependent mean/log2 fold-change errors. Both functions now reject duplicated sample names inGroup.pathway_errorbar()andpathway_errorbar_table()now reject duplicated feature IDs within a selected DAA method/group pair before merging abundance statistics back to annotations, preventing duplicated plot/table rows from many-to-one DAA result inputs.pathway_errorbar()andpathway_errorbar_table()now reject missing, empty, or incompatibleGrouplabels before calculating group means and standard deviations. Previously samples withNAgroup labels could be silently dropped from plotted/table abundance summaries, and mismatched group labels could produce summaries unrelated to the DAA contrast.read_contrib_file(),read_strat_file(), andaggregate_taxa_contributions()now validate contribution identifier keys (sample,function_id, andtaxon) before aggregation. Missing or empty IDs, duplicated wide-format stratified sample columns, and invalid requested pathway/function IDs now fail fast instead of being silently dropped or mislabeled by downstream aggregation.read_contrib_file()andaggregate_taxa_contributions()now reject contribution tables that mix gene-family-level identifiers (such as KOs/ECs) with pathway-level identifiers in the samefunction_idcolumn. This avoids silently treating direct pathway contributions and KEGG pathway-to-KO expansions as one comparable biological unit.aggregate_taxa_contributions()now validates DAA adjusted p-values, significance thresholds,top_n, and contribution metric values. Missing, infinite, or negative contribution values now fail fast instead of being silently dropped or propagated into taxa contribution summaries.aggregate_taxa_contributions()now rankstop_ntaxa by total contribution mass instead of mean contribution over observed rows, avoiding over-selection of sparsely observed taxa with a single high contribution.taxa_contribution_bar()andtaxa_contribution_heatmap()now validate contribution values, identifier keys, requestedfunction_ids, andn_functionsbefore plotting.taxa_contribution_bar(show_percentage = TRUE)now rejects zero-total sample/function combinations instead of displaying undefined relative contributions as all-zero percentage bars.taxa_contribution_heatmap()now treats absent sparse contribution combinations as zero when computing sample means, matching PICRUSt2’s sparse contribution outputs and avoiding inflated mean heatmap intensities for taxa observed in only a subset of samples.taxa_contribution_heatmap(annotation_data = ...)now rejects conflicting labels for the same plotted function ID instead of silently choosing the first duplicate annotation label.taxa_contribution_bar()now ranks defaultn_functionsby between-sample variance in total function contribution, after summing taxa and filling absent sparse sample/function combinations with zero. This prevents taxon composition heterogeneity within a function from being mistaken for between-sample functional variation.Corrected the multi-group ALDEx2 method label from
ALDEx2_Kruskal-Wallace testtoALDEx2_Kruskal-Wallis test; package internals still recognize the legacy spelling as an alias for backward compatibility.
Internal
- Renamed
ordervariable tosort_idxinpathway_errorbar()to avoid shadowing the function parameter; removed unreachable defaultswitch()branch. - Removed redundant pre-sort in
visualize_gsea()enrichment plot (reorder()in the aesthetic handles display ordering). - Trimmed verbose per-sample/per-row log messages in
pathway_heatmap().
ggpicrust2 2.5.16
CRAN release: 2026-05-20
New Features
- Added pathway-level taxa contribution support for PICRUSt2
path_abun_contrib.tsvoutput viaread_pathway_contrib_file()and expandedread_contrib_file(type = "auto")parsing. -
aggregate_taxa_contributions()now supports contribution tables that do not containnorm_taxon_function_contrib, using available PICRUSt2 contribution metrics without requiring a differential abundance step. -
pathway_annotation()now accepts data-frame input, including rowname-basedko2kegg_abundance()output, and can annotate local KEGG pathway IDs withpathway = "KEGG".
ggpicrust2 2.5.14
CRAN release: 2026-04-29
Behavior Changes
-
pathway_daa()now defaultsinclude_effect_size = TRUE. ALDEx2 results includeeffect_size,diff_btw,log2_fold_change,rab_all, andoverlapcolumns by default, aligning ALDEx2 output with DESeq2, edgeR, limma voom, LinDA, Maaslin2, and metagenomeSeq, which all return log2 fold changes without an opt-in flag. The extraALDEx2::aldex.effect()call reuses the already-computed CLR object, so the cost is modest. Passinclude_effect_size = FALSEto restore the p-value-only output (requested in #181). - Multi-group ALDEx2 runs no longer warn that “effect size only available for two-group comparisons”; the flag is silently ignored when
aldex.effect()cannot apply. -
pathway_daa(daa_method = "DESeq2", reference = ...)now honors the user-supplied reference level and supports multi-group designs (one row-block per non-reference contrast), matching the shape returned byedgeRandlimma voom. Previously thereferenceargument was silently ignored and multi-group input yielded only a singleLevel[2]vsLevel[1]contrast.
Bug Fixes
-
ko2kegg_abundance(filter_for_prokaryotes = TRUE)now actually removes eukaryotic / human-system pathway classes from the bundled KEGG hierarchy. The previous filter matched obsolete Level 2 labels without the091xxBRITE prefixes used inko_to_kegg_reference, so it removed no rows.ko2kegg_abundance()also now excludes KEGG BRITE hierarchies and “Not Included in Pathway or Brite” buckets such asko99980before abundance calculation, because they are not pathway maps and cannot be consistently annotated as pathways. Bacterial infection and antimicrobial resistance pathways remain available in the default prokaryotic mode. -
pathway_daa()withinclude_abundance_stats = TRUEno longer produceslog2_fold_change.x/log2_fold_change.ycolumns when the DAA method already returns its ownlog2_fold_change(ALDEx2 with effect size, DESeq2, edgeR, limma voom, LinDA, Maaslin2, metagenomeSeq). The method-native log2 fold change is kept and the relative-abundance ratio is not recomputed, so model-based and ratio-based effect sizes are never conflated under the same column name. -
pathway_daa(include_abundance_stats = TRUE)now fails if the requested abundance summaries cannot be calculated for the returned feature/group pairs, instead of warning and returning a result table that omits the user-requested summary columns. -
pathway_daa(select = ...)now reorders metadata rows to match the reordered abundance columns. Previously theselectbranch reordered abundance but only filtered metadata with%in%, leaving group labels desynchronized from samples wheneverselectwas not in natural order, which silently produced wrong p-values and log fold changes. -
pathway_daa(daa_method = "limma voom")with three or more groups now labelsgroup2correctly. Previously the length-(k-1) contrast vector was recycled into a result ofn_features * (k-1)rows, producing interleavedB,C,B,C,...labels that no longer corresponded to theas.vector()-flattened p-values and coefficients. -
pathway_daa(daa_method = "Maaslin2")with three or more groups now emits one row per (feature, non-reference level) contrast instead of flattening the output viamatch()to a single row per feature. Thegroup2column is taken from Maaslin2’s ownvaluecolumn rather than being recycled from a length-(k-1) vector. -
pathway_daa(daa_method = "metagenomeSeq")no longer hardcodesmetadata$sample; sample-column autodetection fromalign_samples()is honored so metadata with a non-default sample identifier (e.g.SampleID) works end-to-end. -
pathway_daa(daa_method = "metagenomeSeq")now extracts p-values and log-fold changes from the completefitFeatureModel()result using feature identifiers. The previous code calledMRcoefs(), a top-table display helper that can return sorted or partial feature tables and currently fails onfitFeatureModelResultsin supported metagenomeSeq versions; ggpicrust2 then returned all-NAlog2_fold_changevalues with only a warning. -
pathway_daa(daa_method = "metagenomeSeq")with three or more groups now emits one row-block per (reference, non-reference) contrast –(k - 1) * n_featuresrows total – matching the shape returned by DESeq2 / edgeR / limma voom / LinDA / Maaslin2. Previously the function built a full k-column model matrix, calledfitFeatureModel()once, readcoef = 2, and hard-coded the labels asgroup1 = Level[1] / group2 = Level[2], so any contrast beyond the first non-reference level was silently dropped while the output shape looked like a two-group result.fitFeatureModel()is also metagenomeSeq’s documented two-group entry point (it tests a single coefficient and returns one p-value per feature), so each pairwise contrast is now refit on the subset of samples in the two levels of interest. Thereferenceargument governs which level is held fixed asgroup1across all contrasts. - Calling
pathway_daa(daa_method = "Maaslin2")twice in the same R session no longer fails withcannot open the connection. Stale handlers left by Maaslin2’sloggingpackage after its first-call tempdir is cleaned up are now cleared before each invocation. -
pathway_daa()now rejects abundance matrices containing negative values or duplicate sample identifiers at the validation layer, instead of letting them propagate into method-specific failures with cryptic messages. -
pathway_daa()andpathway_errorbar()now refuse sample columns with a total abundance of zero (or NA) instead of silently producing NaN inside thex / sum(x)relative-abundance step. The NaN used to be absorbed by downstreammean(..., na.rm = TRUE)aggregations, so group statistics and error-bar plots were computed from fewer samples than supplied with no warning. The error now names the offending sample(s). The sharedcompute_relative_abundance()helper replaces the duplicatedapply(., 2, function(x) x / sum(x))idiom, andvalidate_abundance()gained acheck_zero_columnsgate so every entry point shares the same contract. -
pathway_daa()now validatesdaa_methodagainst the supported set up front and suggests the canonical spelling for common typos (e.g."linDA"->"LinDA","Lefse"->"Lefser","aldex"->"ALDEx2"). Previously the method dispatch fell throughswitch()with no default branch and returnedNULLsilently, breaking downstream annotation and plotting with opaque errors. -
compare_metagenome_results()now aligns every input metagenome on the shared feature set by row name before thecbindand per-feature Spearman correlation steps. Previously both steps indexed rows by position, so two metagenomes with identical row names but different row orders were compared feature-by-position and produced meaningless (often negative) correlations. A by-name alignment now correctly returns correlation = 1 for identical inputs regardless of row order. The function also errors cleanly when a metagenome is missing row names or when the metagenomes share no features. -
compare_metagenome_results()now aligns every input metagenome on the shared sample set by column name in addition to the existing feature alignment. Previously the per-feature Spearman correlation computedstats::cor(m1[k, ], m2[k, ])by indexing columns by position, so two metagenomes with identical content but columns in different orders were compared “sample i of metagenome A” against “sample i of metagenome B” as if they were the same biological sample, producing median correlations that could be strongly negative (e.g.-0.6) on data that was really identical under by-name alignment. Mismatched column counts also used to fall through tostats::cor()and abort mid-loop with “incompatible dimensions”; the function now stops at the boundary with an actionable message when column names are missing, sample sets are disjoint, and warns when the intersection drops samples. Per-sample cross-metagenome correlation is only defined for parallel samples; this makes that contract explicit instead of enforcing it by happy-path coincidence. -
find_sample_column()now requires the Priority 1 standard-named column (e.g.sample,Sample,sample_id,sample_name) to contain unique values that match the abundance sample IDs. Previously a standard-named column with duplicate values (e.g.sample = c("S1","S1","S2","S2")) was picked purely on its name, silently misaligning every downstream function that relies onalign_samples(). -
ggpicrust2()’s internal call topathway_daa()now uses the modernp_adjust_methodargument. Previously it still forwarded the deprecatedp.adjustargument, so every normalggpicrust2()call emitted the deprecation warning that is meant to fire only when a user explicitly supplies the legacy name. The legacyp.adjustparameter remains accepted with a deprecation warning for backward compatibility. -
pathway_errorbar_table()now delegates sample-column detection and Group reordering to the sharedalign_samples()helper, removing a parallel implementation of the same logic as well as a duplicatedlength(Group) != ncol(abundance)check. Accepted metadata shapes are now identical topathway_daa()andggpicrust2(). -
DESCRIPTIONnow listsMaaslin2andmetagenomeSequnderSuggests. Both packages are documented as supported values in and are dispatched via , but neither appeared inSuggests, so the declared dependency graph, the public method list, and the runtime dispatch had drifted apart. Declaring them brings the three sources back into alignment. -
pathway_gsea(organism = ...)andprepare_gene_sets(organism = ...)now warn when a non-default value is supplied. Both functions advertised anorganismargument but the KEGG and GO branches read KO-based reference tables (ko_to_kegg_reference,ko_to_go_reference) that are organism-independent by construction, so a caller passing e.g.organism = "hsa"silently got the same gene sets as"ko". The argument is retained for signature compatibility with a deprecation warning, its documentation now records the no-op, and the parameter will be removed in a future release. Callers that rely on the default value are unaffected. -
ggpicrust2()no longer forwards itsselectargument intopathway_daa(). The wrapper’s@param selectdocuments a vector of pathway names for plot-time feature selection, but the same value was also being passed intopathway_daa()whereselectmeans sample names, so any call withselect = <pathway names>aborted immediately insidepathway_daa()with “Some selected samples not in abundance data”. The user-suppliedselectis now only forwarded topathway_errorbar()– where feature-level filtering for the figure actually happens – andpathway_daa()runs on the full sample set as documented. -
pathway_errorbar()no longer silently overwrites a method-nativelog2_fold_changecolumn with a relative-abundance mean ratio. The previous code added the column as NA only when missing, then unconditionally overwrote every row inside aforloop, so the effect size displayed in the side panel disagreed with the model output that produced the p_adjust shown next to it. The bar is now taken as-is when the DAA method supplieslog2_fold_change(DESeq2, edgeR, limma voom, LinDA, Maaslin2, metagenomeSeq, and ALDEx2 withinclude_effect_size = TRUE), so both panels of the same figure report the same model-based estimate. The mean-ratio fallback still runs when nolog2_fold_changecolumn is supplied (ALDEx2 withinclude_effect_size = FALSE, Lefser, or custom DAA frames). -
pathway_errorbar()andpathway_errorbar_table()now require every displayed significant DAA feature to be present in the abundance row names. The previous intersection-based subsetting could silently drop significant DAA rows from abundance summaries or leave plot panels describing different feature sets.pathway_errorbar()also validates displayed method-nativelog2_fold_changevalues, so non-finite effect sizes fail before plotting. -
pathway_errorbar_table()no longer derives its two group names viaunique(daa_results_filtered_sub_df$group1)[1]/unique(daa_results_filtered_sub_df$group2)[1]. The surroundingvalidate_daa_results()call already hard-rejects multi-contrast input, so theunique(...)[1]idiom was dead defensive code – but it was also shaped exactly like “silently pick the first of many”, which would have masked any future validator bypass by collapsing a(k-1) * n_featuresmulti-contrast DAA result (as produced bypathway_daa()for >=3 groups with DESeq2 / edgeR / limma voom / LinDA / Maaslin2 / metagenomeSeq) down to a single contrast without warning. Both names are now read via direct[1]indexing, keeping the fast-fail path routed through the validator where contract violations belong. -
pathway_errorbar()andpathway_errorbar_table()(viacalculate_abundance_stats()) now derive per-feature, per-group mean and standard deviation from a single shared helper,summarize_abundance_by_group(). Previouslypathway_errorbar()rolled its ownpivot_longer() %>% group_by(name, group) %>% summarise(mean(value), sd(value))path withoutna.rm = TRUE, so the same abundance matrix could produce different bar heights in the plot versus the companion table if any NA slipped through the pipeline. Unifying the aggregation removes that latent divergence and guarantees both entry points evolve together.
Internal
- Documented the intentional duplicate
align_samples()call inggpicrust2()andpathway_daa(). The wrapper must pre-align abundance/metadata before Step 4 buildsGroup_vecby positional zipping ofmetadata[[group]]withcolnames(abundance);pathway_daa()must also align independently to honor its standalone-caller contract. The two call sites are deliberately invoked with identical arguments, andalign_samples()is deterministic and idempotent, so running it twice on the same inputs is a cheap no-op and cannot drift. A regression test intest-data_utils.Rnow locks the idempotency invariant so any future change that breaks it fails loudly at test time. - Removed the
"nonsense"placeholder columns and values frompathway_errorbar()’s internal data frames. The log2-fold-change bar now sets its fill directly ongeom_bar()instead of routing a single color throughaes(fill = group_nonsense)+scale_fill_manual(), and the pathway-class / p-value side panels anchor all labels at a single x viaaes(x = "")instead of padding each data frame with a constant dummy column. A dead$group2 <- "nonsense"column that was written but never read downstream is also gone. No user-visible change. -
pathway_annotation(file = ..., ko_to_kegg = FALSE)now actually populates thedescriptioncolumn. The file-mode branch previously extracted features from sample column names (skipping columns 1 and 2 and taking the rest as IDs), so the description column was always filled withNA. Feature IDs are now read from the first column, in one unified code path shared with the DAA-results branch. -
pathway_daa(daa_method = "metagenomeSeq")no longer aborts withmissing value where TRUE/FALSE neededon small or near-uniform inputs (e.g. the minimum 4-sample / 2-group case). The normalization quantile returned bycumNormStatFast()is now checked for NA/NaN and falls back to metagenomeSeq’s documented default (p = 0.5). -
pathway_errorbar_table(sample_col = ...)now defaults to auto-detection via the same logicalign_samples()uses, so metadata usingsample,Sample,sample_id, etc. works without passingsample_colexplicitly. Previously the default was hardcoded to"sample_name", which caused the commonsampleconvention to fail withColumn 'sample_name' not found in metadata. -
pathway_daa(daa_method = "LinDA", reference = ...)now actually contrasts against the user-specified reference level. The formula passed toMicrobiomeStat::linda()did not relevel the grouping factor, so LinDA kept the factor’s natural first level as reference while thegroup1label in the result used the user’sreferenceargument. That produced rows withgroup1 == group2and anlog2_fold_changewhose sign did not reflect the requested direction. -
pathway_daa(daa_method = "Maaslin2", reference = ...)now honors the requested reference in the two-group case. Previously thereferenceargument was only forwarded to Maaslin2 for k > 2 groups and the two-group branch passedNULL, so Maaslin2 fell back to its alphabetical default while the result labeledgroup1 = <user reference>– producing the samegroup1 == group2/no-sign-flip symptom as the LinDA bug above. -
pathway_daa()now re-validates the group count after sample alignment andselectfiltering. A narrowselect =that removes every sample of a level, oralign_samples()dropping the only samples for a group, would previously let a single-group dataset reach the backends with a less actionable downstream error. -
pathway_annotation(ko_to_kegg = TRUE)now returns the full inputdaa_results_dfwith annotation columns, populatingpathway_name,pathway_description,pathway_class, andpathway_maponly for rows wherep_adjust < p_adjust_thresholdand leaving the rest asNA. Previously it returned just the significant subset, silently dropping non-significant rows that downstream code (e.g.ggpicrust2()’splot_result_list$daa_results_df) expected to still be present. -
pathway_annotation(ko_to_kegg = TRUE, organism = ...)no longer picks the first organism-specific gene linked to a KO as a “representative” and fetches its record in place of the KO entry. Gene order fromKEGGREST::keggLink()is not semantically meaningful, so isozymes or paralogs participating in different pathways could yield different annotations across KEGG builds. The function now always fetches the generic KO entry (the authoritative KO-level record) and rewrites pathway IDs from thekoprefix to the organism prefix (e.g.ko00010→hsa00010), which is KEGG’s own convention for organism-specific pathway projection. -
find_sample_column()(internal) no longer picks up categorical columns or columns with only partial overlap when auto-detecting the sample identifier column in metadata. The scan-every-column fallback now requires unique values and >= 90% overlap with the abundance sample names; standard-named columns (sample,sample_id, etc.) retain the previous lenient threshold because the column name is itself strong evidence. This preventsalign_samples()from mistakenly treating asubject_idorbatchcolumn as the sample ID when it happens to share a few strings with the sample names. -
pathway_daa(daa_method = "edgeR", reference = ...)now honors the user-supplied reference. edgeR’sexactTest()tookpair = c(1, 2)against the raw factor order, soreferencewas silently ignored and the result always labeledgroup1 = Level[1]/group2 = Level[2]. The grouping factor is now releveled so the tested contrast islog(non-ref / ref)and the labels reflect the requested direction. Multi-group edgeR runs now emit one block per (reference, non-reference) contrast rather than every pairwise combination, matching the shape returned by DESeq2 / limma voom / LinDA / Maaslin2. -
pathway_daa(daa_method = "metagenomeSeq", reference = ...)now honors the user-supplied reference. The model matrix used the raw factor order and the result labels were hardcoded toLevel[1]/Level[2], so flippingreferenceleft both labels and coefficients unchanged. The grouping factor is releveled beforefitFeatureModel()so the contrast and labels track the requested direction. -
taxa_contribution_heatmap(annotation_data = ...)now accepts the column shape actually produced bypathway_annotation(). The heatmap used to look upannotation_data$pathway/$description, butpathway_annotation()emitsfeature(ID) anddescription(non-ko_to_kegg) orpathway_name(ko_to_kegg). The lookup silently returned all-NA, leaving raw IDs on the axis. The heatmap now detectsfeature/descriptionandpathway/pathway_namecolumn pairs and relabels correctly.
ggpicrust2 2.5.13
Bug Fixes
-
ggpicrust2()now rejects the incompatible combination ofko_to_kegg = TRUEwithpathway = "EC"or"MetaCyc"up front, with an actionable error message. Previously this misuse produced a crypticNo features in abundance dataerror thrown from deep insidepathway_daa()(reported in #198). -
ko2kegg_abundance()now errors when none of the input feature IDs match the expected KO format (e.g. when EC numbers are passed in), instead of silently producing an empty result. Total-mismatch detection in the internalvalidate_feature_ids()helper was previously a blind spot. -
ko2kegg_abundance()also errors when the input contains only KO IDs that are absent from the KEGG reference, rather than returning an empty data frame that would break downstream DAA.
ggpicrust2 2.5.11
Documentation
- Unified the package website URL around the canonical
https://cafferyang.com/ggpicrust2/. - Simplified README citation guidance to use the paper DOI and
citation("ggpicrust2")instead of redundant publisher-specific links. - Replaced a moved LinDA reference URL with its DOI-based link and cleaned a KO-to-GO manual encoding warning.
ggpicrust2 2.5.10
CRAN release: 2026-02-12
Bug Fixes
- Removed
Maaslin2fromSuggeststo avoid CRAN/BioC availability failures on special check flavors where the package is not in mainstream repos. - Switched the internal MaAsLin2 call in
pathway_daa()to dynamic lookup so the method remains optional without forcing repository availability checks.
ggpicrust2 2.5.9
Bug Fixes
- Refactored
compare_daa_results()examples to use minimal in-memory DAA-like result tables instead of running external method pipelines. - Removed hard dependency on optional method packages in examples (notably
Maaslin2) so CRAN specialdonttestchecks can run without failing.
ggpicrust2 2.5.8
Bug Fixes
- Fixed CRAN pretest failure in
tests/testthat/test-pathway_daa.Rby splitting default/core method coverage from optional extended methods that require non-mainstream dependencies (e.g.,Maaslin2). - Kept extended DAA method coverage available behind
GGPICRUST2_RUN_EXTENDED_DAA_TESTS=true. - Corrected
metacyc_referencedocumentation to match actual data columns (id,description), resolvingcodocmismatch warnings.
ggpicrust2 2.5.7
Bug Fixes
- Added a regression test for
pathway_errorbar()to explicitly cover theko_to_kegg = TRUE+order = "pathway_class"path. - This guards against reintroducing the historical
tibble::column_to_rownames()/Can't find column '.'failure mode.
ggpicrust2 2.5.6
Breaking Changes
-
Removed
ggpicrust2_extended()function:- This wrapper function provided minimal value over calling
ggpicrust2()andpathway_gsea()separately - Users can achieve the same functionality by calling the individual functions directly
- This change reduces maintenance burden and improves code clarity
- This wrapper function provided minimal value over calling
Major Features
Covariate Adjustment & Improved Statistical Methods for pathway_gsea() (#193)
-
Added limma camera and fry methods as new GSEA options:
-
method = "camera"(now default): Competitive gene set test using limma’s camera function -
method = "fry": Fast rotation gene set test (self-contained) - Both methods account for inter-gene correlations, providing more reliable p-values than preranked GSEA
-
-
Added covariate adjustment support:
- New
covariatesparameter to adjust for confounding factors (age, sex, BMI, etc.) - Covariates are incorporated into the design matrix for proper statistical adjustment
- Essential for microbiome studies where host factors can confound results
- New
-
New parameters:
-
covariates: Character vector of covariate column names from metadata -
contrast: For multi-group comparisons, specify the contrast to test -
inter.gene.cor: Inter-gene correlation for camera method (default: 0.01)
-
-
Scientific background:
- Wu et al. (2012) demonstrated that preranked GSEA methods can produce “spectacularly wrong p-values” due to not accounting for inter-gene correlations
- The camera and fry methods from limma address this limitation
- Reference: Wu, D., & Smyth, G. K. (2012). Nucleic Acids Research, 40(17), e133.
-
Backward compatibility:
- Existing
fgseaandclusterProfilermethods remain available - Added informational message when using preranked methods about p-value reliability
- Existing
-
New internal functions:
-
run_limma_gsea(): Core implementation for camera/fry methods -
build_design_matrix(): Constructs design matrix with covariates
-
-
Updated documentation and vignettes:
- Method selection guide with comparison table
- Covariate adjustment examples
- Updated gsea_analysis.Rmd vignette
New Visualization Functions
-
Added
pathway_volcano()function:- Creates publication-quality volcano plots for differential abundance analysis
- Visualizes both statistical significance (-log10 p-value) and effect size (log2 fold change)
- Smart label placement using ggrepel to avoid overlapping labels
- Color-coded significance categories (Up/Down/Not Significant)
- Automatic handling of NA pathway names (won’t display “NA” labels)
- Handles infinite p-values (when p = 0) gracefully
- Customizable thresholds, colors, and appearance
-
Added
pathway_ridgeplot()function:- Creates ridge plots (joy plots) for GSEA results interpretation
- Shows distribution of gene abundances/fold changes within enriched pathways
- Color-coded by enrichment direction (Up/Down)
- Automatic pathway-KO mapping using built-in ko_to_kegg_reference data
- Supports KEGG and GO pathway types
- Helps identify whether pathways are predominantly up- or down-regulated
- Requires ggridges package (added to Suggests)
-
New dependencies added to Suggests:
-
ggridges: Required for ridge plot visualization -
ggrepel: Used for smart label placement in volcano plots
-
ggpicrust2 2.5.5
Major Features
Prokaryote-Specific Pathway Filtering (#191)
-
New
filter_for_prokaryotesparameter inko2kegg_abundance():- Defaults to TRUE, automatically filtering out eukaryote-specific pathways
- Removes biologically irrelevant pathways from bacterial/archaeal analysis:
- Cancer pathways (overview and specific types)
- Neurodegenerative diseases (Alzheimer’s, Parkinson’s, etc.)
- Substance dependence (addiction pathways)
- Cardiovascular diseases
- Endocrine and metabolic diseases (human-specific)
- Immune diseases (human-specific)
- Organismal systems (immune, nervous, endocrine, digestive, etc.)
- Retains prokaryote-relevant pathways:
- All Metabolism pathways
- Infectious disease: bacterial (Salmonella, E. coli, Tuberculosis, etc.)
- Drug resistance: antimicrobial (antibiotic resistance)
- Genetic/Environmental Information Processing
- Cellular Processes
- Set
filter_for_prokaryotes = FALSEfor eukaryotic analysis or to include all pathways - Reduces pathway count from ~370 to ~290 for typical bacterial analyses
Local KEGG Database Implementation (#113)
-
Replaced KEGG API dependency with local database:
- Implemented comprehensive local KO-to-KEGG pathway mapping (61,655 mappings)
- Covers 557 pathways and 27,127 KO IDs
- Eliminates dependency on external KEGG API
- Provides 100% data coverage with 0% missing values (vs 84% NA in previous format)
- Significantly improved performance and reliability
-
Enhanced data structure:
- Migrated from wide format (306×326 matrix) to long format (61,655×9 table)
- Added rich pathway metadata including hierarchical classification (Level1-3)
- Includes pathway names, KO descriptions, and EC numbers
- Optimized with fast lookup index for O(N) performance
-
Updated functions:
-
ko2kegg_abundance(): Now uses internal database instead of KEGG API -
pathway_gsea(): Updated to use long-format data for gene set enrichment - Both functions maintain backward compatibility
- Added comprehensive input validation and error handling
-
-
PICRUSt 2.6.2 compatibility:
- Automatic detection and cleaning of “ko:” prefix in KO IDs
- Handles both old (K##### format) and new (ko:K##### format) PICRUSt2 outputs
- Seamless migration path for existing users
-
Data quality improvements:
- Validates KO ID format (K##### pattern)
- Detects negative abundance values
- Reports missing values with detailed statistics
- Identifies all-zero KOs across samples
- Checks for duplicate column names
-
Performance:
- Processing speed: 2,145-6,452 KOs/second
- Handles large datasets (1000+ KOs, 50+ samples) efficiently
- Progress bar for long-running operations
-
Testing:
- Added comprehensive test suite with 64 tests
- Covers data structure, functionality, performance, and edge cases
- All tests passing with 100% coverage of new features
This resolves Discussion #113 and provides a robust, API-independent solution for KEGG pathway analysis.
ggpicrust2 2.5.4
Improvements
Enhanced Error Messages for pathway_annotation() (#142)
-
Significantly improved diagnostic messages when no significant pathways are found:
- Provides clear explanation when p_adjust < 0.05 filter removes all pathways
- Shows detailed statistics (total features, significant count, minimum p-value)
- Lists possible reasons (sample size, effect size, variability, method choice)
- Offers concrete recommendations (data quality checks, alternative methods)
- Suggests using ko_to_kegg = FALSE for local annotation without KEGG API
-
Enhanced error messages for KEGG API failures:
- Detailed diagnostic information when all queries fail
- Distinguishes between “not found (HTTP 404)” and “network errors”
- Lists affected KO IDs for troubleshooting
- Provides specific recommendations based on error type
- Suggests local annotation as reliable alternative
-
Updated documentation:
- Clarified p_adjust < 0.05 filtering behavior in function description
- Added note about NA columns when no significant pathways exist
- Improved parameter descriptions for ko_to_kegg parameter
- Enhanced return value documentation with filtering behavior
-
User experience improvements:
- All messages follow R package standards (warning(), stop() instead of cat())
- No emoji characters (professional text output)
- Clear formatting for readability
- Actionable information to help users resolve issues quickly
This addresses Discussion #142 where users received NA annotations without understanding why.
Bug Fixes
pathway_gsea() MetaCyc Support Fix (#174)
-
Fixed undefined organism parameter bug:
- Added
organism = "ko"parameter with default value - Prevents “object ‘organism’ not found” error
- Properly passes organism to prepare_gene_sets() function
- Added
-
Added input validation for MetaCyc data:
- Detects when MetaCyc pathway IDs are provided instead of EC numbers
- Issues clear warning directing users to use pathway_daa() for pathway-level data
- Helps users understand the difference between gene-level and pathway-level analysis
-
Enhanced documentation:
- Clarified that GSEA requires gene-level data (EC numbers for MetaCyc)
- Added cross-reference to pathway_daa() for pathway abundance analysis
- Improved parameter descriptions to prevent data type confusion
-
Scientific integrity maintained:
- Ensures correct analysis method for each data type
- Prevents misleading results from incorrect data usage
- Guides users to appropriate functions based on their data
Critical annotation_custom() Fix (#184)
-
Fixed annotation_custom() parameter type issue in pathway_errorbar():
- Removed unit object wrappers from
annotation_custom()position parameters - Now uses numeric values directly for xmin, xmax, ymin, ymax parameters
- Resolves “no applicable method for ‘rescale’ applied to an object of class ‘c(’simpleUnit’, ‘unit’, ‘unit_v2’)’” error
- Fixes pathway class background color rendering failures
- Removed unit object wrappers from
-
Root cause identified and resolved:
-
annotation_custom()expects numeric values, not unit objects for position parameters - Previous fix attempt (changing ggplot2::unit to grid::unit) was incorrect
- Both ggplot2::unit() and grid::unit() return identical objects - the issue was using unit objects at all
- Theme-related unit usage (legend.key.size, plot.margin) remains unchanged and correct
-
This fix resolves the critical rendering issue where users encountered errors when generating pathway error bar plots with pathway class backgrounds, particularly with R 4.4+ and ggplot2 4.0.0.
ggpicrust2 2.5.3
Previous Release
- Initial attempt to fix issue #184 (superseded by v2.5.4)
ggpicrust2 2.5.2
CRAN release: 2025-08-25
Bug Fixes
Critical Edge Case Resolution
-
Fixed pathway_errorbar empty data handling:
-
pathway_errorbar()now returns NULL instead of crashing when no annotation data is available - Main
ggpicrust2()function gracefully handles NULL plot objects - Users still receive complete results data even when plots cannot be generated
- Improved warning messages for better user experience
-
-
Enhanced robustness for datasets with no significant pathways:
- Prevents crashes when all pathways have p_adjust > 0.05
- Returns meaningful results with empty annotation columns for visualization
- Maintains data integrity throughout the analysis pipeline
- Works seamlessly with all PICRUSt2 versions including 2.6.2
-
Function signature consistency:
- Fixed
pathway_annotation()function calls in main function - Ensures proper parameter passing throughout the workflow
- Resolves compatibility issues between internal functions
- Fixed
These fixes ensure the package works reliably with all types of microbiome data, including edge cases where no statistically significant pathways are found.
ggpicrust2 2.5.1
Bug Fixes
PICRUSt 2.6.2 Compatibility (#174) - Complete Resolution
-
Added automatic KO ID format detection and conversion:
- Automatically detects PICRUSt 2.6.2 format with “ko:” prefixes
- Transparently removes “ko:” prefixes during data loading
- Maintains full backward compatibility with PICRUSt 2.5.2 format
- Eliminates “subscript out of bounds” errors caused by format mismatches
-
Enhanced core functions for seamless compatibility:
- Updated
ko2kegg_abundance()with automatic format detection - Updated
ggpicrust2()main function with compatibility layer - Added clear informational messages about format conversion
- Zero manual preprocessing required for users
- Updated
-
Comprehensive testing and validation:
- Tested with real PICRUSt 2.6.2 output files (5,000+ KO features)
- Verified compatibility with all major DAA methods
- Confirmed backward compatibility with existing workflows
- Performance optimized for large datasets
ggpicrust2 2.5.0
Bug Fixes
PICRUSt 2.6.2 Compatibility (#174)
-
Fixed compatibility issues with PICRUSt 2.6.2 output:
- Added comprehensive data validation for PICRUSt compatibility
- Improved zero-abundance data filtering with fallback strategies
- Enhanced error handling with specific PICRUSt version guidance
- Added graceful handling of sparse data scenarios
-
Enhanced ALDEx2 error detection:
- Better error messages for insufficient or invalid data
- Improved handling of edge cases (single features, all-zero data)
- Added PICRUSt version-specific troubleshooting guidance
-
Improved LinDA analysis robustness:
- Better filtering of zero-abundance features before analysis
- Enhanced data validation and error reporting
- Added warnings for very sparse data scenarios
-
Enhanced ko2kegg_abundance function:
- Added fallback strategies for zero-abundance KO data
- Improved compatibility warnings and error messages
- Better handling of PICRUSt format variations
-
Added comprehensive compatibility guide:
- Created detailed troubleshooting documentation
- Provided alternative analysis strategies
- Added version-specific recommendations
ggpicrust2 2.4.0
New Features
Enhanced Multi-Grouping Support for pathway_heatmap()
-
Added secondary_groups parameter (#171):
- Enables multi-level grouping in pathway heatmaps
- Supports nested faceting with multiple grouping variables
- Maintains full backward compatibility with existing code
- Automatic color generation for complex grouping structures
-
Deprecated facet_by parameter:
- Replaced with more flexible secondary_groups parameter
- Shows deprecation warning with migration guidance
- Will be removed in future major version
-
Enhanced Examples and Documentation:
- Added comprehensive examples for multi-grouping usage
- Migration guide from facet_by to secondary_groups
- Updated function documentation with new parameter descriptions
-
Comprehensive Testing:
- Added unit tests for all new functionality
- Integration tests with real data
- Backward compatibility validation
- Error handling and edge case testing
ggpicrust2 2.3.3
Major Bug Fixes and Improvements
DAA Methods Comprehensive Fixes
-
Fixed limma voom compatibility with compare_daa_results (#163):
- Added missing group1 and group2 columns for two-group comparisons
- Ensures consistent output format with other DAA methods
- Resolves “Unknown comparison type” error when using limma voom with compare_daa_results
-
Fixed Maaslin2 multiple critical issues (#164):
- Fixed sample name matching between abundance matrix and metadata
- Implemented intelligent feature name matching to handle Maaslin2’s hyphen-to-dot conversion
- Resolves “Unable to find samples in data and metadata files” error
- Eliminates NA p-values caused by feature name mismatches
- Robust handling of various feature naming conventions (hyphens, dots, underscores)
-
Enhanced DESeq2 robustness:
- Added fallback mechanism for dispersion estimation failures
- Uses gene-wise estimates when standard dispersion estimation fails
- Improved compatibility with small sample sizes and low-variance data
-
Standardized Lefser implementation:
- Added lefser package to method detection list
- Standardized output format to match other DAA methods
- Returns results for all features, not just significant ones
- Converts effect scores to p-values for consistency
Comprehensive Testing and Validation
-
100% success rate achieved for all 8 DAA methods:
- ALDEx2, DESeq2, edgeR, limma voom, metagenomeSeq, Maaslin2, LinDA, Lefser
- All methods now fully compatible with compare_daa_results function
- Consistent output format across all methods
- Extensive testing with various feature naming patterns and edge cases
ggpicrust2 2.3.2
CRAN release: 2025-07-15
Bug Fixes
- Fixed MetaCyc pathway annotation NA description issue (#154):
- Standardized MetaCyc reference data column names from ‘X1’/‘X2’ to ‘id’/‘description’
- Applied fix in all three loading paths of load_reference_data function
- Resolves column name mismatch that caused all MetaCyc annotations to return NA
- Tested with 100% success rate on sample data
- Maintains backward compatibility with KO and EC pathway types
ggpicrust2 2.3.1
Bug Fixes
- Fixed MetaCyc reference data loading issue:
- Enhanced the file search mechanism in the
load_reference_datafunction - Added multiple search paths for reference data files
- Improved error messages with more diagnostic information
- Fixed “Reference data file not found” error that some users encountered
- Enhanced the file search mechanism in the
ggpicrust2 2.3.0
Bug Fixes
- Fixed NA handling in pathway_errorbar function:
- Improved handling of NA values in the feature column
- Added robust error checking for group ordering option
- Prevents errors when processing MetaCyc pathway data with missing annotations
- Enhanced pathway_pca function to handle zero variance data:
- Automatically detects and filters out rows (pathways) with zero variance
- Keeps sample profiles as PCA observations, including samples with zero variance across pathways
- Provides informative warnings about removed pathways
- Improves error handling with clear diagnostic messages
- Fixed file extension handling in ko2kegg_abundance function:
- Resolved “the condition has length > 1” error when processing files
- Improved extension detection for different file types
- Enhanced robustness for handling various input formats
ggpicrust2 2.2.1
Bug Fixes
- Added aggregate_by_group parameter to pathway_heatmap function:
- Allows displaying representative samples (e.g., mean) for each group instead of all individual samples
- Improved handling of NA values in heatmap visualization
- Added aggregate_fun parameter for customizing the aggregation function
ggpicrust2 2.1.4
Reference Data Updates
- Updated reference databases for improved 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
ggpicrust2 2.1.1
Bug Fixes
- Improved error handling for KEGG database connections (Issue #138):
- Added robust error handling for HTTP 404 errors when connecting to the KEGG database
- Function now continues processing other KO IDs even when some IDs return HTTP 404 errors
- Added detailed logging about which KO IDs were not found
- Only throws a fatal error when all KO IDs fail to be processed
- Includes summary statistics about successful, not found, and error counts
- Fixed the LinDA analysis for multi-group comparisons (Issue #144):
- Modified the
perform_linda_analysisfunction to handle multi-group comparisons correctly - The function now creates separate result entries for each pairwise comparison
- Added the
log2FoldChangecolumn to the results for effect size information
- Modified the
ggpicrust2 2.1.0
Major Changes
-
Added Gene Set Enrichment Analysis (GSEA) functionality with the following new functions:
-
pathway_gsea(): Performs GSEA analysis, supporting KEGG, MetaCyc, and GO pathways -
visualize_gsea(): Creates visualizations of GSEA results, including enrichment plots, dotplots, barplots, network plots, and heatmaps -
compare_gsea_daa(): Compares GSEA and Differential Abundance Analysis (DAA) results -
gsea_pathway_annotation(): Adds pathway annotations to GSEA results
-
Improved network and heatmap visualization capabilities with richer parameter options and better error handling
Added preliminary support for MetaCyc and GO pathways
Fixed various bugs and optimized code structure
ggpicrust2 2.0.0
CRAN release: 2025-04-03
主要变更
Refactored the package dependencies, moving most Bioconductor packages from Imports to Suggests, reducing mandatory dependencies.
Added conditional checks to ensure that packages are only used when they are available.
Fixed the function export issue, ensuring that all public API functions are correctly exported.
Updated the example code, improving its stability and compatibility.
Updated the documentation format to comply with the latest CRAN standards.
Fixed invalid URL links in the README.
Optimized code quality, removing unused variables.
Added missing import declarations to ensure package integrity.
