ggpicrust2: an R package for PICRUSt2 predicted functional profile analysis and visualization

1One limitation of microbial community marker-gene sequencing is that it does not provide information about the functional composition of sample communities. In recent years, several methods have been developed to predict functions from 16S rRNA sequence data, including PICRUSt2, Tax4Fun2, MicFunPred, and PICRUSt, among others. However, the accuracy and applicability of these methods depend on the specific research question and the characteristics of the microbial community being studied. Overall, these methods have greatly enhanced our ability to understand the functional roles of microbial communities in various environments, from the human gut to soil and water ecosystems. Among the various tools available, PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) has emerged as a highly favored instrument for predicting functional profiles, as it facilitates the generation of comprehensive pathway abundances within microbial communities. By doing so, PICRUSt2 provides researchers with valuable insights into the functional roles of microbial communities.

2Nonetheless, a consensus regarding the optimal methodology for inferring and visualizing the functional abundance output generated by PICRUSt2 remains to be established within the academic community. As determining the statistically significant differences in functions and pathways between groups using differential abundance (DA) methods constitutes a critical step in the analysis, selecting an appropriate DA approach is indeed a topic of considerable importance within the scholarly discourse. The official wiki of PICRUSt2 initially recommended STatistical Analysis of taxonoMic and functional Profiles (STAMP) as the preferred software for analysis and visualization. However, STAMP has not been updated since 2015, indicating that it is unable to integrate the most recent advances in DA analysis, which are crucial for systematically making statistical inferences from PICRUSt2 output data. Furthermore, STAMP presents installation challenges on macOS platforms, making it less user-friendly and potentially hindering its adoption by researchers in the field. The performance of five DA methods supported by STAMP, including ANOVA, Kruskal–Wallis H-test, t-test (equal variance), Welch’s t-test, and White’s non-parametric test, has been shown to be relatively inferior in a recent comparison of 20 DA methods across 38 datasets. The comparison concluded that AlDEx2 and ANCOM-II produce the most consistent results across studies and agree best with the interest of results from different approaches, but still recommend that researchers should use a consensus approach based on multiple DA methods to help ensure robust biological interpretations. Despite there are several platforms or packages that support the multiple advanced DA methods such as MicrobiomeAnalyst, MicrobiomeExplorer, microbiomeMarker, they are not specifically designed for PICRUSt2 functional output data. Due to the discrepancies in format and characteristics between PICRUSt2 output data and 16S rRNA genes data, the above platforms or software intended for the analysis of 16S rRNA genes data often encounter difficulties when importing PICRUSt2 data. Although almost all DA methods can be used in R, each method creates various burdens for data import and parameter configuration which increases both the effort and time cost, and diminishes the efficiency. Additionally, these R packages often lack the ability to visualize DA results and generate publication-quality figures. Thus, developing a user-friendly R package for analyzing PICRUSt2 functional output data using various DA methods and visualizations is urgently needed to fill the gaps.

3The general workflow of the package is shown in Fig. 1. ggpicrust2 not only allows for recently developed advanced DA methods and visualization of results but also can convert PICRUSt2 output Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology (KO) abundance tables into KEGG pathway abundance tables, which cannot be performed using PICRUSt2 alone. It also provides annotation of KO, Enzyme Commission (EC), MetaCyc pathway, and KEGG pathway and enables classification of KEGG pathways. In the future, ggpicrust2 plans to incorporate a broader array of functional prediction tools, including but not limited to Tax4Fun2, in order to expand its capabilities and utility. Additionally, the package will integrate other methods that have demonstrated strong performance in simulation comparisons, ensuring continuous improvement and alignment with the latest advancements in the field.

4ggpicrust2 recommends adopting the data format of PICRUSt2 original output pred_metagenome_unstrat.tsv without reformat. But csv and txt are also acceptable. Furthermore, it is capable of supporting files that have been converted to mimic the PICRUSt2 output format, ensuring compatibility and flexibility for various data sources.

5KO is a classification system developed by the KEGG database. It uses a hierarchical structure to classify enzymes based on the reactions they catalyze. To better understand pathways’ role in different groups and classify the pathways, the KO abundance table can be converted to KEGG pathway abundance. But PICRUSt2 removes the function from PICRUSt. ko2kegg_abundance can help convert the table.

6DA analysis plays a major role in PICRUSt2 downstream analysis. pathway_daa integrates almost all DA methods applicable to the predicted functional profile which there excludes ANCOM and ANCOMBC. It includes ALDEx2, DEseq2, Maaslin2, LinDA, edgeR, limma voom, metagenomeSeq, and lefser which have demonstrated varying degrees of success in distinct benchmarking assessments.pathway_daa provides a convenient way to run these methods and compare the results. compare_daa_results can be used to compare the consistency and inconsistency of statistically significant features obtained using different methods. It creates a report showing the number of common and different features identified by each method, and the features themselves. This helps researchers choose the method most suitable for their dataset and research question.

7For guidance on selecting the most appropriate DA method for your study, Supplementary Table S1 provides a brief comparison and description of various DA methods, assisting researchers in selecting the best method based on the characteristics of their dataset and research question.

8pathway_annotation can annotate the KO, EC, MetaCyc pathways’ description from annotation table. The KO database is a database of molecular functions represented in terms of functional orthologs. The EC number is a numerical classification scheme for enzymes, based on the chemical reactions they catalyze. MetaCyc pathways describe the biochemical transformations of small molecules within an organism. And it can pull requests to online KEGG database to annotate KEGG pathways’ pathway_name, pathway_description, pathway_class, and pathway_map. The function can be used to annotate the output file of PICRUSt2 or the output table of pathway_daa.

9The mainstream visualization of PICRUSt2 is bar_plot, error_bar_plot, pca_plot, heatmap_plot. pathway_errorbar can show the relative abundance difference between groups and log2 fold change and P-values derived from DA results. pathway_pca can show the difference after dimensional reduction via principal component analysis. pathway_heatmap can visualize the patterns in PICRUSt2 output data which can be useful for identifying trends or highlighting areas of interests.

10ggpicrust is the integration function of pathway_daa, pathway_annotation, pathway_errorbar, and ko2kegg_abundance. This tool is designed to facilitate the entire data analysis process for those who are new to the field. However, it is also capable of being used by professional analysts in a modular fashion, allowing for increased customization and control. To further support users and promote the understanding of the package’s capabilities, we have developed a detailed user manual, which is provided as Supplementary Materials. This document includes step-by-step installation instructions, explanations of the main features, and guidance on how to effectively leverage the ggpicrust2 package within academic research. Our aim is to ensure that both novice and experienced researchers can easily access and benefit from the package’s advanced functionalities.

11After employing PICRUSt2 to perform functional profile prediction on microbiome data from C9orf72 loss of function mice, including those that underwent fecal transplantation and those that did not, our subsequent data analysis using ggpicrust2 entailed the implementation of LinDA. This approach led to the identification of KEGG pathways that demonstrated statistically significant differences between the prosurvival and proinflammatory environments across both groups of mice. Of particular interest were the pathways ko05016, which is primarily involved in the pathogenesis of Huntington’s disease, and ko05012, known for its association with Parkinson’s disease. Both pathways are linked to human diseases and neurodegenerative disorders. The DA results were meticulously annotated and the output was visualized for subsequent analysis. The visual representation of the results, which provides insights into the involvement of these pathways in the studied conditions, is depicted in Fig. 1.

12compare_metagenome_results analyzes and compares the functional predictions from different methods and sequencing metagenomes. It accepts a list of metagenome count matrices, their corresponding sample names, a method for DA, a method for P-value adjustment, and a reference group level for DA. The function concatenates all metagenome count matrices, creates new sample metadata, performs DA, and calculates Spearman correlation coefficients and corresponding P-values between each pair of metagenomes.

13ggpicrust2, available at CRAN and https://github.com/cafferychen777/ggpicrust2, is an R package developed explicitly for PICRUSt2 predicted functional profile to do advanced DA analysis and visualization of the DA results. This package effectively addresses the limitations of existing tools in terms of methods and visualization, and its integrated and distributed design caters to both professionals and beginners by meeting the needs of both groups. By providing a seamless experience for analyzing and visualizing DA results, ggpicrust2 has the potential to significantly enhance the quality and efficiency of research involving functional profile predictions. ggpicrust2 has already been incorporated into the PICRUSt2 wiki documentation, reflecting its growing recognition and adoption within the research community.