Streamlining feature elaboration and statistics analysis in metabolomics: the GetFeatistics R-package

2.1 Targeted analyses

The package provides straightforward calculation of absolute concentrations from signal intensities via the function get_targeted_elaboration. The function’s arguments include: (I) a table of intensities for samples, quality controls, and calibration standards used to build the curve; (II) a table reporting the nominal concentrations for quality controls and calibration levels; and (III) optionally, a table mapping each analyte to its internal standard and specifying the desired weighting scheme. Because calibration data often have higher variance at higher concentrations, applying weighted least squares (e.g., 1/Y, 1/X or 1/X²), reduces bias and improves accuracy at low concentrations [18], 19]. The function returns calculated absolute concentrations and accuracy (%) for standards and quality controls. These results can be passed to plot_calibration_curves to visualise the fitted calibration curves.

2.2 Untargeted analyses

Following feature extraction, grouping, and annotation using external tools such as patRoon or MS-DIAL, the workflow of GetFeatistics is then based to be conducted on two core data frames: “featTable”, where rows represent features and columns represent samples, containing intensity values; and “featINFO”, where rows again represent features, while columns contain the feature metadata such as m/z, retention time, chemical annotations, and other descriptors of the feature or of the annotated molecule. Dedicated functions are available to import and structure these data frames from patRoon (get_feat_info_from_patRoon, get_feat_table_from_patRoon) and MS-DIAL outputs (get_feat_info_from_MSDial, get_feat_table_from_MSDial). Confidence levels of annotated compounds can be automatically assigned according to the annotation grading proposed by Schymanski et al. [20], particularly implementing the MoNAScore cut-offs, for patRoon, and the dot product and fragment presence cut-offs, for MS-DIAL, as reported by Talavera Andújar et al. [21]. Known molecules can be searched in the feature list based on their exact mass; moreover, if the retention time is also available, i.e.: when analytical standards have been analysed under the same instrumental conditions, matching features can be identified and a confidence level 1 is then assigned (checkmolecules_in_feat_table).

In order to improve the quality of the dataset and retain only biologically meaningful features, the package includes a function to filter features based on QC criteria (QCs_proces) by filtering out features with: (I) poor reproducibility in pooled QCs (based on relative standard deviation of QC intensities); (II) insufficient presence across QC injections (based on detection frequency in QCs); (III) high blank contribution (based on blank/QC mean intensity ratio); (IIII) unreasonable response in diluted QC samples (based on a range of acceptability of QC/diluted QC mean intensity ratio). When group-specific QC pools are available, the filtering retains features meeting quality thresholds in at least one group [3].

2.3 Statistical analyses

A summary table with descriptive statistics can be generated (gentab_descr) considering the sample groups. Normality is assessed via density plots (test_normality_density_plot), Q–Q plots (test_normality_q_q_plot), and the Shapiro–Wilk test (test_normality_Shapiro_table). Data can then be processed via missing value imputation, log-transformation, centering, and scaling (transf_data). A suite of functions supports univariate statistical testing, including t-tests (gentab_P.t.test), fold change analysis (gentab_FC, gentab_FC_more_than2levels), one- and two-way ANOVA, with or without interaction terms (gentab_P.1wayANOVA_posthocTukeyHSD, gentab_P.2wayANOVA_posthocTukeyHSD). Considering the application of such statistical tests on several metabolites, p-values can be corrected using false discovery rate (FDR) procedures [22]. Results can be visualized through volcano plots (Volcano_ttest_FC) and box plots (getBoxplots). The outcome tables from statistical analyses can be enriched with corresponding feature metadata (addINFO_to_table).

For multivariate analysis, a single function (gentab_lm_long) allows the construction of regression models across the entire metabolite dataset, incorporating one or more covariates, with the option of building linear models with fixed effects, mixed-effects models with the “lmerTest” package [23], 24], or Tobit regression for censored data (specifically suited in this workflow for values below or above detection limits) via the “AER” package [25]. The output of the function is a table with complete effect estimates, with related confidence intervals and p-values (with also FDR correction), for each metabolite – covariate pair. Visual summaries of results through a Volcano plot can be created (Volcano_lm). Principal component analyses and heatmap with hierarchical clustering can be also performed with extensive customisation options (getPCA, getHeatMap).

2.4 Miscellaneous

The package also contains functions (getChemData, add_ChemData_to_featINFO) to automatically retrieve chemical properties, synonyms, and identifiers from PubChem using the Power User Gateway (PUG) [26], and chemical ontology classification via the ClassyFire taxonomy [27]. Additionally, a function (do_FELLA_enrichment_analysis) is provided for pathway enrichment analysis using the FELLA package [28].

Finally, the package includes a functionality (merge_results) to merge results from multiple chromatographic runs and/or ionization modes, prioritising overlapping compounds according to user-defined criteria (e.g., lower mass error, retention time shift, or p-value).