Error in pathway_heatmap(): unused arguments (show_rownames = TRUE, show_colnames = TRUE) and missing "group" argument

[ShahariarArif]

❓ Question / Help Request

📋 What do you need help with?

• [+] Understanding how to use a specific function
• Choosing the right analysis method
• Interpreting results
• Data preparation and formatting
• Installation or setup issues
• Best practices for my research question
• Troubleshooting unexpected results
• Other: ___________

🔬 Research Context

Briefly describe your research question and study design:

Study type:

• Case-control study
• Time series / longitudinal
• Cross-sectional
• [+] Experimental design
• Other: ___________

📊 Your Data

Data source:

• [ +] PICRUSt2 output (pred_metagenome_unstrat.tsv)
• Processed abundance data
• Example data from package
• Other: ___________

Data characteristics:

• Number of samples: 9
• Number of groups/conditions: 3
• Sample size per group: 3
• Type of microbiome data: [16S]

💻 What you've tried

Current approach:

# Please share your current code (even if it's not working)
library(ggpicrust2)

# Your code here...

library(ggpicrust2)
library(tidyverse)
library(patchwork)
library(RColorBrewer)

Paths

path_abun <- "/home/arif/16S_QIIME2/picrust2_input/picrust2_out/pathways_out/path_abun_unstrat.tsv"
path_meta <- "/home/arif/16S_QIIME2/metadata.tsv"

Load abundance

abundance <- read_tsv(path_abun, show_col_types = FALSE)

Load metadata – force "group" column

metadata <- read_tsv(path_meta, show_col_types = FALSE) %>%
column_to_rownames("sample-id") %>%
select(group) # ← exactly "group"

Color palette

my_colors <- colorRampPalette(rev(brewer.pal(11, "RdYlBu")))(100)

Level 1

p1 <- pathway_heatmap(
abundance = abundance,
metadata = metadata,
cluster_rows = TRUE,
cluster_cols = TRUE
) + ggtitle("KEGG Level 1 Pathways")

Level 2

p2 <- pathway_heatmap(
abundance = abundance,
metadata = metadata,
cluster_rows = TRUE,
cluster_cols = TRUE
) + ggtitle("KEGG Level 2 Pathways")

Level 3 (top 50)

top50 <- abundance %>%
column_to_rownames("pathway") %>%
rowSums() %>%
sort(decreasing = TRUE) %>%
names() %>%
head(50)

p3 <- pathway_heatmap(
abundance = abundance %>% filter(pathway %in% top50),
metadata = metadata,
cluster_rows = TRUE,
cluster_cols = TRUE
) + ggtitle("KEGG Level 3 (Top 50)")

Combine

combined <- p1 + p2 + p3 + plot_layout(ncol = 1)

Save

dir.create("kegg_heatmaps", showWarnings = FALSE)
ggsave("kegg_heatmaps/Level1.png", p1, width = 12, height = 10, dpi = 300)
ggsave("kegg_heatmaps/Level2.png", p2, width = 16, height = 14, dpi = 300)
ggsave("kegg_heatmaps/Level3_Top50.png", p3, width = 20, height = 18, dpi = 300)
ggsave("kegg_heatmaps/All_Levels.png", combined, width = 20, height = 30, dpi = 300)

cat("Heatmaps saved! Check folder: kegg_heatmaps/\n")

Current results or issues:
Describe what you're getting vs. what you expected.

🎯 Specific Questions

1. 🎯 Specific Questions

2. Why does pathway_heatmap() throw "unused arguments (show_rownames = TRUE, show_colnames = TRUE)" even though these are listed in some documentation/examples?
3. When I remove those arguments, I get "argument 'group' is missing, with no default" — but my metadata has a column exactly named "group".
4. Is there a different function or updated syntax/script

for creating level-wise (Level 1/2/3) KEGG pathway heatmaps in the current version?
4. What is the recommended way to create clustered heatmaps grouped by sample metadata (e.g., Control vs Treatment groups)?

📚 Background Research

• I've read the package documentation
• I've reviewed the vignettes
• [ +] I've searched existing issues
• I've consulted relevant literature

Helpful resources you've already checked:

🖥️ Environment (if relevant)

• **Operating System: Ubuntu 24.04.3 LTS
• **R Version: 4.5.2
• **ggpicrust2 Version: 2.5.2

🔄 Analysis Goals

What do you ultimately want to achieve?

• Identify differentially abundant pathways
• [ +] Create publication-quality figures
• Compare different methods
• Perform statistical analysis
• Data exploration and visualization
• Other: ___________

⏰ Timeline

• This is urgent for a deadline
• [+] I have some flexibility with timing
• This is for future planning

📖 Additional Context

Any other information that might be helpful for understanding your question or providing better guidance.

metadata.tsv

[cafferychen777]

Hi @ShahariarArif,

Thank you for reporting this issue. I've investigated and this appears to be a usage issue rather than a bug in the package.

Issue 1: "unused arguments (show_rownames = TRUE, show_colnames = TRUE)"

The pathway_heatmap() function uses show_row_names (with underscore), not show_rownames. Also, there is no show_colnames parameter - column labels are displayed through facet grouping.

You may have confused this with ComplexHeatmap package syntax, which uses different parameter names.

Issue 2: "argument 'group' is missing, with no default"

The group parameter is required and has no default value. You must specify which column in your metadata contains the grouping information.

Looking at your code:

p1 <- pathway_heatmap(
  abundance = abundance,
  metadata = metadata,
  cluster_rows = TRUE,
  cluster_cols = TRUE
)  # ← Missing group parameter!

Correct Usage

p1 <- pathway_heatmap(
  abundance = abundance,
  metadata = metadata,
  group = "group",           # Required - specify your grouping column name
  cluster_rows = TRUE,
  cluster_cols = TRUE,
  show_row_names = TRUE      # Note: use underscore, not camelCase
)

For Your Specific Case

Since your metadata has a column named "group", use:

library(ggpicrust2)
library(tidyverse)

# Load data
abundance <- read_tsv(path_abun, show_col_types = FALSE) %>%
  column_to_rownames("pathway")

metadata <- read_tsv(path_meta, show_col_types = FALSE)

# Create heatmap - note the required 'group' parameter
p1 <- pathway_heatmap(
  abundance = abundance,
  metadata = metadata,
  group = "group",
  cluster_rows = TRUE,
  cluster_cols = TRUE
)

The documentation in the README and package help pages (?pathway_heatmap) shows the correct usage with all required parameters.

Best regards

[ShahariarArif]

Hi cafferychen777,
Thanks for your reply. Mainly I am a self learner and very new in metagenomic data analysis... i would like to know one more thing. As my path_abun_unstrat.tsv only contain pathway name how will be they assigned to particular hierarchial level (1, 2, 3)? does ggpicrust2 have any innate algorithom to assigne them to particular level? i am asking this as i am still facing error for plotting the graph in level 1, level 2 and level 3. I hope i will get more guide line from you. Thanks again.

[cafferychen777]

Hi @ShahariarArif,

Thanks for the screenshots - now I can see the issue clearly!

The Problem

Looking at your path_abun_unstrat.tsv file (Image 2), your data contains MetaCyc pathway IDs such as:

• 12DICHLORETHDEG-PWY
• AEROBACTINSYN-PWY
• BRANCHED-CHAIN-AA-SYN-PWY

However, your code (Image 1) is trying to create "KEGG Level 1 Pathways" heatmaps.

This won't work because MetaCyc and KEGG are completely different databases.

Database	Example ID	Has Level 1/2/3?
MetaCyc	AEROBACTINSYN-PWY	❌ No
KEGG	ko00010	✅ Yes

The error There are only 0 columns occurs because the function cannot find any KEGG hierarchy information in your MetaCyc data.

Solution

You have two options:

Option 1: Analyze MetaCyc pathways directly (without hierarchy)

library(ggpicrust2)
library(tidyverse)

# Load your MetaCyc data
abundance <- read_tsv("path_abun_unstrat.tsv") %>%
  column_to_rownames("pathway")

metadata <- read_tsv("metadata.tsv")

# Run DAA analysis
daa_results <- pathway_daa(
  abundance = abundance,
  metadata = metadata,
  group = "group"
)

# Annotate with MetaCyc descriptions
annotated_results <- pathway_annotation(
  pathway = "MetaCyc",
  daa_results_df = daa_results,
  ko_to_kegg = FALSE
)

# Create heatmap (no hierarchy levels for MetaCyc)
p <- pathway_heatmap(
  abundance = abundance,
  metadata = metadata,
  group = "group"
)

Option 2: Use KO data for KEGG hierarchy analysis

If you want Level 1/2/3 hierarchy, you need KO (KEGG Orthology) abundance from PICRUSt2, not MetaCyc.

1. Check if you have a file like pred_metagenome_unstrat.tsv (KO abundance) from your PICRUSt2 output
2. Use ko2kegg_abundance() to convert:

# Convert KO to KEGG pathways (this gives you level1, level2, level3)
kegg_abundance <- ko2kegg_abundance(
  file = "pred_metagenome_unstrat.tsv"
)

Summary

• Your current data: MetaCyc pathways → Use Option 1
• For KEGG Level 1/2/3: Need KO abundance data → Use Option 2

Let me know which approach works better for your research goals!

Best regards

[ShahariarArif]

Thank you again... I would like to follow the option 2. Can i get the full pipeline for option 2?

[cafferychen777]

Hi @ShahariarArif,

Great choice! Here's the complete pipeline for Option 2 (KEGG pathway analysis with hierarchy levels).

Prerequisites

First, make sure you have the KO abundance file from PICRUSt2. This is typically named pred_metagenome_unstrat.tsv and is located in your PICRUSt2 output folder (usually in the KO_metagenome_out directory).

The file should look like this:

function	Sample1	Sample2	Sample3	...
K00001	100	150	200	...
K00002	50	75	100	...

Complete Pipeline

# Load required libraries
library(ggpicrust2)
library(tidyverse)
library(ggprism)

# ============================================
# STEP 1: Load your data
# ============================================

# Path to your PICRUSt2 KO abundance file
ko_file <- "/home/arif/16S_QIIME2/picrust2_input/picrust2_out/KO_metagenome_out/pred_metagenome_unstrat.tsv"

# Load metadata
metadata <- read_tsv("path/to/your/metadata.tsv", show_col_types = FALSE)

# ============================================
# STEP 2: Convert KO abundance to KEGG pathway abundance
# ============================================

# This function converts KO IDs to KEGG pathways
# The output includes pathway hierarchy information (level1, level2, level3)
kegg_abundance <- ko2kegg_abundance(file = ko_file)

# View the result
head(kegg_abundance)
# Rows are KEGG pathway IDs (e.g., ko00010, ko00020)
# Columns are your samples

# ============================================
# STEP 3: Perform Differential Abundance Analysis (DAA)
# ============================================

daa_results <- pathway_daa(
  abundance = kegg_abundance,
  metadata = metadata,
  group = "group",        # Change to your grouping column name
  daa_method = "LinDA"    # You can also use "ALDEx2", "DESeq2", etc.
)

# Filter for significant results
significant_results <- daa_results %>% 
  filter(p_adjust < 0.05)

print(paste("Found", nrow(significant_results), "significant pathways"))

# ============================================
# STEP 4: Annotate pathways with KEGG information
# ============================================

# This adds pathway names, descriptions, and hierarchy levels
annotated_results <- pathway_annotation(
  pathway = "KO",
  daa_results_df = daa_results,
  ko_to_kegg = TRUE
)

# View annotated results
head(annotated_results)

# ============================================
# STEP 5: Visualize results
# ============================================

# --- Option A: Error bar plot (grouped by pathway class) ---
p_errorbar <- pathway_errorbar(
  abundance = kegg_abundance,
  daa_results_df = annotated_results,
  Group = metadata$group,      # Change to your grouping column
  p_values_threshold = 0.05,
  order = "pathway_class",     # Order by KEGG hierarchy
  ko_to_kegg = TRUE,
  p_value_bar = TRUE,
  x_lab = "pathway_name"
)
print(p_errorbar)

# Save the plot
ggsave("kegg_errorbar.pdf", p_errorbar, width = 12, height = 10)

# --- Option B: Heatmap ---
# Filter abundance for significant pathways only
significant_pathways <- significant_results$feature

p_heatmap <- pathway_heatmap(
  abundance = kegg_abundance[rownames(kegg_abundance) %in% significant_pathways, ],
  metadata = metadata,
  group = "group",             # Change to your grouping column
  cluster_rows = TRUE,
  cluster_cols = TRUE
)
print(p_heatmap)

# Save the heatmap
ggsave("kegg_heatmap.pdf", p_heatmap, width = 10, height = 8)

# --- Option C: PCA plot ---
p_pca <- pathway_pca(
  abundance = kegg_abundance,
  metadata = metadata,
  group = "group"              # Change to your grouping column
)
print(p_pca)

# ============================================
# STEP 6: Analyze by KEGG hierarchy levels (OPTIONAL)
# ============================================

# Load the reference data to see hierarchy information
data("ko_to_kegg_reference")

# View available hierarchy levels
unique(ko_to_kegg_reference$level1)
# Example output:
# "09100 Metabolism"
# "09120 Genetic Information Processing"
# "09130 Environmental Information Processing"
# etc.

# Filter for specific Level 1 category (e.g., Metabolism only)
metabolism_pathways <- ko_to_kegg_reference %>%
  filter(level1 == "09100 Metabolism") %>%
  pull(pathway_id) %>%
  unique()

# Create heatmap for Metabolism pathways only
metabolism_abundance <- kegg_abundance[rownames(kegg_abundance) %in% metabolism_pathways, ]

p_metabolism <- pathway_heatmap(
  abundance = metabolism_abundance,
  metadata = metadata,
  group = "group",
  cluster_rows = TRUE
) + ggtitle("KEGG Level 1: Metabolism Pathways")
print(p_metabolism)

Understanding KEGG Hierarchy Levels

Level	Description	Example
Level 1	Broad categories	"Metabolism", "Genetic Information Processing"
Level 2	Sub-categories	"Carbohydrate metabolism", "Amino acid metabolism"
Level 3	Specific pathways	"00010 Glycolysis", "00020 Citrate cycle"

Troubleshooting Tips

1. File not found: Make sure your KO abundance file path is correct
2. No significant results: Try different DAA methods or adjust your p-value threshold
3. Sample name mismatch: Ensure sample names in abundance data match those in metadata

Summary of Key Functions

Step	Function	Purpose
1	ko2kegg_abundance()	Convert KO → KEGG pathways
2	pathway_daa()	Differential abundance analysis
3	pathway_annotation()	Add pathway descriptions
4	pathway_errorbar()	Visualize DA results
5	pathway_heatmap()	Create heatmap
6	pathway_pca()	PCA visualization

Let me know if you have any questions or run into any issues!

Best regards

[ShahariarArif]

Still facing problem may be due to no significant difference... However without filtering for significance and other analysis is it possible to create hierarchial heat map? With the same data an analist created a heatmap, bar graph, flowerchart

Load required libraries

library(ggpicrust2)
library(tidyverse)
library(ggprism)
library(MicrobiomeStat)

============================================

STEP 1: Load your data

============================================

Path to your PICRUSt2 KO abundance file

ko_file <- "/home/arif/16S_QIIME2/picrust2_input/picrust2_out/KO_metagenome_out/pred_metagenome_unstrat.tsv"

Load metadata

metadata <- read_tsv("/home/arif/16S_QIIME2/picrust2_input/metadata.tsv", show_col_types = FALSE)

============================================

STEP 2: Convert KO abundance to KEGG pathway abundance

============================================

This function converts KO IDs to KEGG pathways

The output includes pathway hierarchy information (level1, level2, level3)

kegg_abundance <- ko2kegg_abundance(file = ko_file)
Loading data from file...
Rows: 8070 Columns: 10
── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (1): function
dbl (9): CR1, CR2, CR3, T3R1, T3R2, T3R3, T6R1, T6R2, T6R3

ℹ Use spec() to retrieve the full column specification for this data.
ℹ Specify the column types or set show_col_types = FALSE to quiet this message.
Checking KO ID format compatibility...
Detected PICRUSt 2.6.2 format with 'ko:' prefix. Applying compatibility fix...
Cleaned 8070 KO IDs by removing 'ko:' prefix
Loading KEGG reference data. This might take a while...
Processing 306 KEGG pathways...
|========================================================================================================================| 100%
Total KO matches found: 4693
Number of non-zero pathways before filtering: 249
Removing KEGG pathways with zero abundance across all samples...
Final number of KEGG pathways: 249

View the result

head(kegg_abundance)
CR1 CR2 CR3 T3R1 T3R2 T3R3 T6R1 T6R2 T6R3
ko05340 266374.3 254832.5 313770.4 86219.24 53283.76 62950.15 81663.56 38664.84 117002.6
ko00564 2975573.6 2765268.0 3541623.2 746825.72 652427.96 821192.97 859448.66 362074.06 1520087.1
ko00680 4729766.2 4505052.0 5638723.8 1318328.22 895432.43 1142964.70 1456312.98 696535.91 2125552.1
ko00562 523864.4 624828.8 615134.2 190885.14 72648.31 212291.87 222612.20 139713.14 254812.4
ko03030 1839376.6 1749002.6 2155985.6 932083.92 438907.96 1132638.81 551243.93 244404.67 1236185.6
ko00561 2527919.3 2282609.9 3013212.3 514461.89 321511.16 629850.61 704914.26 283705.23 882751.3

Rows are KEGG pathway IDs (e.g., ko00010, ko00020)

Columns are your samples

============================================

STEP 3: Perform Differential Abundance Analysis (DAA)

============================================

daa_results <- pathway_daa(

• abundance = kegg_abundance,
• metadata = metadata,
• group = "group", # Change to your grouping column name
• daa_method = "LinDA" # You can also use "ALDEx2", "DESeq2", etc.
• )
  Using column 'sample-id' as sample identifier
  Running LinDA analysis...
  Group variable: group with levels: Control, T3, T6
  Reference level: Control
  Number of features: 249, Number of samples: 9
  0 features are filtered!
  The filtered data has 9 samples and 249 features will be tested!
  Imputation approach is used.
  Fit linear models ...
  Completed.
  Warning message:
  In MicrobiomeStat::linda(feature.dat = feature.dat, meta.dat = meta.dat, :
  Some features have less than 3 nonzero values!
  They have virtually no statistical power. You may consider filtering them in the analysis!

Filter for significant results

significant_results <- daa_results %>%

• filter(p_adjust < 0.05)

print(paste("Found", nrow(significant_results), "significant pathways"))
[1] "Found 0 significant pathways"

============================================

STEP 4: Annotate pathways with KEGG information

============================================

This adds pathway names, descriptions, and hierarchy levels

annotated_results <- pathway_annotation(

• pathway = "KO",
• daa_results_df = daa_results,
• ko_to_kegg = TRUE
• )
  KO to KEGG is set to TRUE. Proceeding with KEGG pathway annotations...
  No organism specified. Using generic KO information across all species.
  We are connecting to the KEGG database to get the latest results, please wait patiently.
  Processing pathways in chunks...
  [2026-01-16 21:23:15] WARN: No significant pathways found. Returning data with empty annotation columns.
  Warning message:
  No statistically significant biomarkers found (p_adjust < 0.05).
  Returning results with empty annotation columns for visualization purposes.
  Consider using a less stringent threshold if significant results are expected.

View annotated results

head(annotated_results)
feature method group1 group2 p_values log2FoldChange p_adjust adj_method pathway_name pathway_description pathway_class
1 ko05340 LinDA Control T3 0.9256188 -0.07217627 0.9992453 BH
2 ko00564 LinDA Control T3 0.9240828 -0.07228183 0.9992453 BH
3 ko00680 LinDA Control T3 0.8356823 -0.15966668 0.9992453 BH
4 ko00562 LinDA Control T3 0.9630662 -0.03309169 0.9992453 BH
5 ko03030 LinDA Control T3 0.2668518 0.66276386 0.6578822 BH
6 ko00561 LinDA Control T3 0.5539739 -0.46176674 0.9799257 BH
pathway_map
1
2
3
4
5
6

============================================

STEP 5: Visualize results

============================================

--- Option A: Error bar plot (grouped by pathway class) ---

p_errorbar <- pathway_errorbar(

• abundance = kegg_abundance,
• daa_results_df = annotated_results,
• Group = metadata$group, # Change to your grouping column
• p_values_threshold = 0.05,
• order = "pathway_class", # Order by KEGG hierarchy
• ko_to_kegg = TRUE,
• p_value_bar = TRUE,
• x_lab = "pathway_name"
• )
  The following pathways are missing annotations and have been excluded: ko05340, ko00564, ko00680, ko00562, ko03030, ko00561, ko00440, ko00250, ko00740, ko04664, ko04940, ko00010, ko00195, ko03450, ko00760, ko00920, ko00311, ko00310, ko04146, ko00600, ko04141, ko04142, ko00604, ko04260, ko03040, ko05142, ko04540, ko04712, ko00909, ko00513, ko04973, ko00510, ko05110, ko04974, ko04976, ko00450, ko01051, ko00565, ko00524, ko00300, ko05222, ko00514, ko05416, ko00260, ko00190, ko03440, ko00750, ko00950, ko00592, ko00591, ko00590, ko00062, ko00061, ko03070, ko00253, ko03060, ko04370, ko04730, ko04740, ko00380, ko00500, ko05120, ko04666, ko05322, ko04964, ko04962, ko04960, ko00625, ko00624, ko00627, ko00626, ko00621, ko00620, ko00623, ko00622, ko00270, ko04380, ko00940, ko00941, ko01053, ko00943, ko00100, ko00945, ko01057, ko01056, ko05016, ko04145, ko00071, ko00072, ko04360, ko05219, ko05215, ko05214, ko05213, ko05211, ko05210, ko01055, ko00902, ko00534, ko04910, ko00531, ko04916, ko00532, ko00360, ko00633, ko00630, ko00363, ko00364, ko05130, ko04115, ko04114, ko00121, ko04914, ko00430, ko02010, ko00130, ko04330, ko03050, ko00361, ko00040, ko00730, ko00362, ko01040, ko00603, ko03018, ko04270, ko00281, ko00280, ko03013, ko04626, ko05200, ko03015, ko00312, ko05020, ko05143, ko00523, ko00520, ko00521, ko05146, ko05145, ko00052, ko00053, ko00051, ko00401, ko00400, ko04020, ko00350, ko00480, ko00643, ko00642, ko00640, ko00720, ko00960, ko02020, ko00120, ko00965, ko03008, ko05414, ko00290, ko04614, ko03410, ko05010, ko05012, ko00980, ko00410, ko00983, ko05150, ko00791, ko00790, ko05131, ko00020, ko00710, ko00196, ko02060, ko00340, ko00785, ko00550, ko00650, ko05220, ko02030, ko03320, ko04744, ko04745, ko04530, ko00522, ko04612, ko04621, ko04622, ko00910, ko04971, ko00460, ko04970, ko00830, ko00780, ko00511, ko00970, ko00030, ko04972, ko00232, ko00230, ko04122, ko00540, ko05014, ko00660, ko03020, ko04512, ko00982, ko03010, ko00140, ko05100, ko00860, ko05410, ko00331, ko00330, ko0 ...
  You can use the 'pathway_annotation' function to add annotations for these pathways.
  Warning message:
  All 498 rows were excluded due to missing 'pathway_name' annotations. This typically occurs when no statistically significant pathways were found. Cannot create error bar plot with empty data. Returning NULL.

print(p_errorbar)
NULL

Save the plot

ggsave("kegg_errorbar.pdf", p_errorbar, width = 12, height = 10)

--- Option B: Heatmap ---

Filter abundance for significant pathways only

significant_pathways <- significant_results$feature

p_heatmap <- pathway_heatmap(

• abundance = kegg_abundance[rownames(kegg_abundance) %in% significant_pathways, ],
• metadata = metadata,
• group = "group", # Change to your grouping column
• cluster_rows = TRUE,
• cluster_cols = TRUE
• )
  Error in pathway_heatmap(abundance = kegg_abundance[rownames(kegg_abundance) %in% :
  At least one pathway is required

[cafferychen777]

Hi @ShahariarArif,

I can see the issue now. Your analysis found 0 significant pathways (p_adjust < 0.05), which is causing all downstream visualizations to fail.

Why This Happened

Looking at your output:

• Sample size: Only 9 samples total (3 groups × 3 samples each)
• Warning: Some features have less than 3 nonzero values!
• Result: Found 0 significant pathways

With only 3 samples per group, statistical power is very limited. After multiple testing correction (BH method), it's common to find no significant results.

Solutions

Option 1: Explore with relaxed threshold (for visualization purposes)

# Use a more relaxed threshold for exploration
# Note: This is for exploratory analysis only, not for publication claims

# Filter with relaxed threshold
exploratory_results <- daa_results %>% 
  filter(p_values < 0.1)  # Use raw p-value, not adjusted

print(paste("Found", nrow(exploratory_results), "pathways with p < 0.1"))

# If you found some pathways, continue with visualization
if (nrow(exploratory_results) > 0) {
  # Create heatmap with top pathways (by raw p-value)
  top_pathways <- daa_results %>%
    arrange(p_values) %>%
    slice(1:20) %>%  # Top 20 pathways
    pull(feature)
  
  p_heatmap <- pathway_heatmap(
    abundance = kegg_abundance[rownames(kegg_abundance) %in% top_pathways, ],
    metadata = metadata,
    group = "group",
    cluster_rows = TRUE,
    cluster_cols = TRUE
  )
  print(p_heatmap)
}

Option 2: Visualize all pathways without significance filtering

# Create heatmap of ALL pathways (no significance filter)
# Useful for exploring overall patterns

# Select top 30 most variable pathways
pathway_variance <- apply(kegg_abundance, 1, var)
top_variable_pathways <- names(sort(pathway_variance, decreasing = TRUE))[1:30]

p_heatmap_all <- pathway_heatmap(
  abundance = kegg_abundance[top_variable_pathways, ],
  metadata = metadata,
  group = "group",
  cluster_rows = TRUE,
  cluster_cols = TRUE
)
print(p_heatmap_all)

# PCA plot (doesn't require significant pathways)
p_pca <- pathway_pca(
  abundance = kegg_abundance,
  metadata = metadata,
  group = "group"
)
print(p_pca)

Option 3: Try different DAA methods

# Some methods may be more suitable for small sample sizes
# Try Wilcoxon test (non-parametric, better for small samples)

daa_results_wilcox <- pathway_daa(
  abundance = kegg_abundance,
  metadata = metadata,
  group = "group",
  daa_method = "ALDEx2"  # ALDEx2 includes Wilcoxon test
)

# Check results from Wilcoxon test
wilcox_results <- daa_results_wilcox %>%
  filter(method == "ALDEx2_Wilcoxon rank test") %>%
  filter(p_adjust < 0.05)

print(paste("ALDEx2 Wilcoxon found", nrow(wilcox_results), "significant pathways"))

Important Notes

1. Small sample size limitation: With only 3 samples per group, even real biological differences may not reach statistical significance. This is a fundamental limitation of small studies.

2. Interpretation: If no pathways are significant, it could mean:

   • No real differences between groups
   • Differences exist but sample size is too small to detect them
   • High biological variability within groups

3. Recommendations for future studies:

   • Aim for at least 5-10 samples per group
   • More samples = more statistical power

Quick Summary

Approach	Purpose	Code
Relaxed threshold	Exploratory analysis	filter(p_values < 0.1)
Top variable pathways	Pattern exploration	top 30 by variance
Different DAA method	Alternative statistics	ALDEx2 with Wilcoxon
PCA plot	Overall patterns	pathway_pca()

Try the options above and let me know what you find!

Best regards

[ShahariarArif]

for option 2 can i get the full code? in option 1 there remain a error.:

Load required libraries

library(ggpicrust2)
library(tidyverse)
library(ggprism)
library(MicrobiomeStat)
library(ggdendro)

Path to your PICRUSt2 KO abundance file

ko_file <- "/home/arif/16S_QIIME2/picrust2_input/picrust2_out/KO_metagenome_out/pred_metagenome_unstrat.tsv"

Load metadata

metadata <- read_tsv("/home/arif/16S_QIIME2/picrust2_input/metadata.tsv", show_col_types = FALSE)

This function converts KO IDs to KEGG pathways

The output includes pathway hierarchy information (level1, level2, level3)

kegg_abundance <- ko2kegg_abundance(file = ko_file)
Loading data from file...
Rows: 8070 Columns: 10
── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (1): function
dbl (9): CR1, CR2, CR3, T3R1, T3R2, T3R3, T6R1, T6R2, T6R3

ℹ Use spec() to retrieve the full column specification for this data.
ℹ Specify the column types or set show_col_types = FALSE to quiet this message.
Checking KO ID format compatibility...
Detected PICRUSt 2.6.2 format with 'ko:' prefix. Applying compatibility fix...
Cleaned 8070 KO IDs by removing 'ko:' prefix
Loading KEGG reference data. This might take a while...
Processing 306 KEGG pathways...
|======================================================================================================================| 100%
Total KO matches found: 4693
Number of non-zero pathways before filtering: 249
Removing KEGG pathways with zero abundance across all samples...
Final number of KEGG pathways: 249

View the result

head(kegg_abundance)
CR1 CR2 CR3 T3R1 T3R2 T3R3 T6R1 T6R2 T6R3
ko05340 266374.3 254832.5 313770.4 86219.24 53283.76 62950.15 81663.56 38664.84 117002.6
ko00564 2975573.6 2765268.0 3541623.2 746825.72 652427.96 821192.97 859448.66 362074.06 1520087.1
ko00680 4729766.2 4505052.0 5638723.8 1318328.22 895432.43 1142964.70 1456312.98 696535.91 2125552.1
ko00562 523864.4 624828.8 615134.2 190885.14 72648.31 212291.87 222612.20 139713.14 254812.4
ko03030 1839376.6 1749002.6 2155985.6 932083.92 438907.96 1132638.81 551243.93 244404.67 1236185.6
ko00561 2527919.3 2282609.9 3013212.3 514461.89 321511.16 629850.61 704914.26 283705.23 882751.3
daa_results <- pathway_daa(

• abundance = kegg_abundance,
• metadata = metadata,
• group = "group", # Change to your grouping column name
• daa_method = "LinDA" # You can also use "ALDEx2", "DESeq2", etc.
• )
  Using column 'sample-id' as sample identifier
  Running LinDA analysis...
  Group variable: group with levels: Control, T3, T6
  Reference level: Control
  Number of features: 249, Number of samples: 9
  0 features are filtered!
  The filtered data has 9 samples and 249 features will be tested!
  Imputation approach is used.
  Fit linear models ...
  Completed.
  Warning message:
  In MicrobiomeStat::linda(feature.dat = feature.dat, meta.dat = meta.dat, :
  Some features have less than 3 nonzero values!
  They have virtually no statistical power. You may consider filtering them in the analysis!

Filter with relaxed threshold

exploratory_results <- daa_results %>%

• filter(p_values < 0.1) # Use raw p-value, not adjusted

print(paste("Found", nrow(exploratory_results), "pathways with p < 0.1"))
[1] "Found 132 pathways with p < 0.1"

If you found some pathways, continue with visualization

if (nrow(exploratory_results) > 0) {

• Create heatmap with top pathways (by raw p-value)

• top_pathways <- daa_results %>%

• arrange(p_values) %>%
  

• slice(1:20) %>%  # Top 20 pathways
  

• pull(feature)
  

• p_heatmap <- pathway_heatmap(

• abundance = kegg_abundance[rownames(kegg_abundance) %in% top_pathways, ],
  

• metadata = metadata,
  

• group = "group",
  

• cluster_rows = TRUE,
  

• cluster_cols = TRUE
  

• )
• print(p_heatmap)
• }
  Samples ordered by hierarchical clustering ( complete method, euclidean distance)
  Column order: T6R3, T3R1, T3R3, CR2, CR1, CR3, T3R2, T6R1, T6R2
  Pathways ordered by hierarchical clustering ( complete method, euclidean distance)
  Row order: ko00513, ko04916, ko04722, ko04720, ko04745, ko04744, ko04114, ko04740, ko05214, ko04972, ko04080, ko00902, ko00331, ko04920, ko04310, ko04330, ko05220, ko04142, ko00604, ko00531
  Warning message:
  In ggplot2::guide_colorbar(direction = if (colorbar_position %in% :
  Arguments in ... must be used.
  ✖ Problematic argument:
  • breaks = colorbar_breaks
  ℹ Did you misspell an argument name?

This adds pathway names, descriptions, and hierarchy levels

annotated_results <- pathway_annotation(

• pathway = "KO",
• daa_results_df = daa_results,
• ko_to_kegg = TRUE,
• )
  KO to KEGG is set to TRUE. Proceeding with KEGG pathway annotations...
  No organism specified. Using generic KO information across all species.
  We are connecting to the KEGG database to get the latest results, please wait patiently.
  Processing pathways in chunks...
  [2026-01-19 01:05:55] WARN: No significant pathways found. Returning data with empty annotation columns.
  Warning message:
  No statistically significant biomarkers found (p_adjust < 0.05).
  Returning results with empty annotation columns for visualization purposes.
  Consider using a less stringent threshold if significant results are expected.

View annotated results

head(annotated_results)
feature method group1 group2 p_values log2FoldChange p_adjust adj_method pathway_name pathway_description pathway_class
1 ko05340 LinDA Control T3 0.9256188 -0.07217627 0.9992453 BH
2 ko00564 LinDA Control T3 0.9240828 -0.07228183 0.9992453 BH
3 ko00680 LinDA Control T3 0.8356823 -0.15966668 0.9992453 BH
4 ko00562 LinDA Control T3 0.9630662 -0.03309169 0.9992453 BH
5 ko03030 LinDA Control T3 0.2668518 0.66276386 0.6578822 BH
6 ko00561 LinDA Control T3 0.5539739 -0.46176674 0.9799257 BH
pathway_map
1
2
3
4
5
6

; though p<0.1 Warning message:
No statistically significant biomarkers found (p_adjust < 0.05).
Returning results with empty annotation columns for visualization purposes.
Consider using a less stringent threshold if significant results are expected.

[cafferychen777]

Hi @ShahariarArif,

Great progress! Your Option 1 heatmap was successfully created (the warning about guide_colorbar is harmless and doesn't affect the plot).

Here's the complete code for Option 2 and a solution to get pathway annotations without the p_adjust < 0.05 filter.

Option 2: Complete Code (Visualize by Variance, No Significance Filter)

# Load required libraries
library(ggpicrust2)
library(tidyverse)
library(ggprism)

# ============================================
# STEP 1: Load data (you already did this)
# ============================================
ko_file <- "/home/arif/16S_QIIME2/picrust2_input/picrust2_out/KO_metagenome_out/pred_metagenome_unstrat.tsv"
metadata <- read_tsv("/home/arif/16S_QIIME2/picrust2_input/metadata.tsv", show_col_types = FALSE)
kegg_abundance <- ko2kegg_abundance(file = ko_file)

# ============================================
# STEP 2: Select top 30 most variable pathways
# ============================================
# Calculate variance for each pathway across all samples
pathway_variance <- apply(kegg_abundance, 1, var)

# Sort by variance and get top 30
top_variable_pathways <- names(sort(pathway_variance, decreasing = TRUE))[1:30]

# View selected pathways
print("Top 30 most variable pathways:")
print(top_variable_pathways)

# ============================================
# STEP 3: Get pathway annotations manually
# ============================================
# Load the reference data
data("ko_to_kegg_reference")

# Create a lookup table for pathway information
pathway_info <- ko_to_kegg_reference %>%
  select(pathway_id, pathway_name, level1, level2, level3) %>%
  distinct()

# Get annotations for your top pathways
top_pathway_annotations <- pathway_info %>%
  filter(pathway_id %in% top_variable_pathways)

# View the annotations
print("Pathway annotations:")
print(top_pathway_annotations %>% select(pathway_id, pathway_name, level1))

# ============================================
# STEP 4: Create heatmap with top variable pathways
# ============================================
p_heatmap_variance <- pathway_heatmap(
  abundance = kegg_abundance[top_variable_pathways, ],
  metadata = metadata,
  group = "group",
  cluster_rows = TRUE,
  cluster_cols = TRUE,
  show_row_names = TRUE
)
print(p_heatmap_variance)

# Save the heatmap
ggsave("kegg_heatmap_top30_variance.pdf", p_heatmap_variance, width = 12, height = 10)

# ============================================
# STEP 5: Create PCA plot (all pathways)
# ============================================
p_pca <- pathway_pca(
  abundance = kegg_abundance,
  metadata = metadata,
  group = "group"
)
print(p_pca)

# Save PCA plot
ggsave("kegg_pca.pdf", p_pca, width = 8, height = 6)

# ============================================
# STEP 6: Create heatmap with pathway names (BONUS)
# ============================================
# Create a named vector for better heatmap labels
pathway_labels <- pathway_info %>%
  filter(pathway_id %in% top_variable_pathways) %>%
  select(pathway_id, pathway_name) %>%
  distinct() %>%
  deframe()  # Convert to named vector

# Rename rows in abundance matrix
abundance_with_names <- kegg_abundance[top_variable_pathways, ]
rownames(abundance_with_names) <- pathway_labels[rownames(abundance_with_names)]

# Remove NA row names if any
abundance_with_names <- abundance_with_names[!is.na(rownames(abundance_with_names)), ]

# Create heatmap with pathway names
p_heatmap_named <- pathway_heatmap(
  abundance = abundance_with_names,
  metadata = metadata,
  group = "group",
  cluster_rows = TRUE,
  cluster_cols = TRUE,
  show_row_names = TRUE
)
print(p_heatmap_named)

# Save
ggsave("kegg_heatmap_named.pdf", p_heatmap_named, width = 14, height = 10)

# ============================================
# STEP 7: Group by KEGG Level 1 (BONUS)
# ============================================
# See what Level 1 categories your pathways belong to
level1_summary <- pathway_info %>%
  filter(pathway_id %in% rownames(kegg_abundance)) %>%
  count(level1, sort = TRUE)

print("KEGG Level 1 distribution in your data:")
print(level1_summary)

# Filter for specific Level 1 category (e.g., Metabolism)
metabolism_pathways <- pathway_info %>%
  filter(level1 == "09100 Metabolism") %>%
  filter(pathway_id %in% rownames(kegg_abundance)) %>%
  pull(pathway_id) %>%
  unique()

# Get top 20 most variable metabolism pathways
metabolism_variance <- pathway_variance[names(pathway_variance) %in% metabolism_pathways]
top_metabolism <- names(sort(metabolism_variance, decreasing = TRUE))[1:20]

# Create metabolism-only heatmap
p_metabolism <- pathway_heatmap(
  abundance = kegg_abundance[top_metabolism, ],
  metadata = metadata,
  group = "group",
  cluster_rows = TRUE
)
print(p_metabolism)

ggsave("kegg_metabolism_heatmap.pdf", p_metabolism, width = 12, height = 8)

Why pathway_annotation() Returns NA

The function pathway_annotation(..., ko_to_kegg = TRUE) only annotates pathways with p_adjust < 0.05 by default. Since your data has no significant results after multiple testing correction, all annotations are NA.

Solution: Use the manual annotation approach in Step 3 above, which gives you full control over which pathways to annotate.

Summary of Outputs

File	Description
kegg_heatmap_top30_variance.pdf	Top 30 most variable pathways
kegg_heatmap_named.pdf	Same, but with pathway names instead of IDs
kegg_pca.pdf	PCA plot showing sample clustering
kegg_metabolism_heatmap.pdf	Top 20 metabolism pathways only

About the Warning

The warning Arguments in ... must be used for guide_colorbar is a minor ggplot2 compatibility issue and does not affect your plot. Your heatmap was created successfully.

Let me know if you need any clarification!

Best regards

[ShahariarArif]

But how will i be able to generate hierarchy level heatmap or even the data file? Does the ggpicrust2 package can map the KEGG pathway?

[cafferychen777]

Thanks for the detailed report. Here is what’s going on in ggpicrust2 v2.5.x and how to fix it.

1. show_rownames / show_colnames errors
   In the current version, the argument name is show_row_names (with row), and there is no show_colnames parameter. That’s why you get ‘unused arguments’. If you want to hide column names, use a custom theme:

p <- pathway_heatmap(
  abundance = abun_mat,
  metadata = metadata,
  group = "group",
  show_row_names = TRUE,
  custom_theme = ggplot2::theme(axis.text.x = ggplot2::element_blank())
)

2. group is required (must be a column name string)
   pathway_heatmap() requires group = "group" explicitly. Also, avoid dropping sample IDs when you select:

metadata <- readr::read_tsv(path_meta, show_col_types = FALSE) %>%
  dplyr::rename(sample_id = `sample-id`) %>%
  dplyr::select(sample_id, group)

abun_mat <- abundance %>%
  dplyr::rename(pathway_id = pathway) %>%
  tibble::column_to_rownames("pathway_id") %>%
  as.matrix()

p <- pathway_heatmap(
  abundance = abun_mat,
  metadata = metadata,
  group = "group",
  cluster_rows = TRUE,
  cluster_cols = TRUE
)

3. How to get KEGG hierarchy (Level 1/2/3)
   path_abun_unstrat.tsv does not contain hierarchy columns. ggpicrust2 can map KEGG IDs to hierarchy using the built-in reference table ko_to_kegg_reference. Example: add levels, then aggregate by level for heatmaps:

data("ko_to_kegg_reference", package = "ggpicrust2")

kegg_levels <- ko_to_kegg_reference %>%
  dplyr::distinct(pathway_id, level1, level2, level3)

abun_annot <- abundance %>%
  dplyr::rename(pathway_id = pathway) %>%
  dplyr::left_join(kegg_levels, by = "pathway_id")

# Level 1 heatmap
level1_mat <- abun_annot %>%
  dplyr::group_by(level1) %>%
  dplyr::summarise(dplyr::across(where(is.numeric), sum, na.rm = TRUE), .groups = "drop") %>%
  tibble::column_to_rownames("level1") %>%
  as.matrix()

p1 <- pathway_heatmap(level1_mat, metadata, group = "group", cluster_rows = TRUE, cluster_cols = TRUE)

# Level 2 heatmap
level2_mat <- abun_annot %>%
  dplyr::group_by(level2) %>%
  dplyr::summarise(dplyr::across(where(is.numeric), sum, na.rm = TRUE), .groups = "drop") %>%
  tibble::column_to_rownames("level2") %>%
  as.matrix()

p2 <- pathway_heatmap(level2_mat, metadata, group = "group", cluster_rows = TRUE, cluster_cols = TRUE)

# Level 3 (pathway IDs)
level3_mat <- abun_annot %>%
  dplyr::filter(!is.na(level3)) %>%
  dplyr::select(-level1, -level2, -level3) %>%
  tibble::column_to_rownames("pathway_id") %>%
  as.matrix()

p3 <- pathway_heatmap(level3_mat, metadata, group = "group", cluster_rows = TRUE, cluster_cols = TRUE)

Short answer to your last question: yes, ggpicrust2 can map KEGG pathways to hierarchy levels using the internal ko_to_kegg_reference table. You just need to join it and aggregate by level1/level2/level3 before calling pathway_heatmap().

[cafferychen777]

TL;DR minimal working example

library(ggpicrust2)
library(tidyverse)

# 1) Load inputs
abundance <- read_tsv("/home/arif/16S_QIIME2/picrust2_input/picrust2_out/pathways_out/path_abun_unstrat.tsv", show_col_types = FALSE)
metadata <- read_tsv("/home/arif/16S_QIIME2/picrust2_input/metadata.tsv", show_col_types = FALSE) %>%
  rename(sample_id = `sample-id`) %>%
  select(sample_id, group)

# 2) Map KEGG levels (built-in reference)
data("ko_to_kegg_reference", package = "ggpicrust2")
levels <- ko_to_kegg_reference %>% distinct(pathway_id, level1, level2, level3)

abun_annot <- abundance %>%
  rename(pathway_id = pathway) %>%
  left_join(levels, by = "pathway_id")

# 3) Level 1 heatmap (aggregate by level1)
level1_mat <- abun_annot %>%
  group_by(level1) %>%
  summarise(across(where(is.numeric), sum, na.rm = TRUE), .groups = "drop") %>%
  column_to_rownames("level1") %>%
  as.matrix()

p1 <- pathway_heatmap(level1_mat, metadata, group = "group", cluster_rows = TRUE, cluster_cols = TRUE)
print(p1)

Notes:

• Use show_row_names (not show_rownames). show_colnames is not a supported argument.
• group is required and must be a metadata column name (string).
• If you want Level 2/3, aggregate by level2/level3 exactly the same way as level1.