A Robust CWM–PCA for Chemospace: Horn-based Dimensioning, Principal-angle Stability, and Dirichlet Pooling of Isomers

A Robust CWM–PCA for Chemospace: Horn-based Dimensioning, Principal-angle Stability, and Dirichlet Pooling of Isomers Xinyao Yang1, Lijuan Sun2, Jingjing Wang2 Fanbin Zeng2, Qingxia Li3 1Xi’an Jiaotong Liverpool University, Suzhou, Jiangsu, China 2Lanzhou University, Lanzhou, Gansu, China 3Fisk University, Nashville, Tennessee, USA Abstract We present a community-weighted PCA for species chemospace. Compounds are aggregated to species with two weightings: presence (detection proportion) and relative (softmax of logabundance). Isomer candidates are pooled using Dirichlet priors (A1 symmetric; A2 similarityweighted). Dimensions are selected by Horn’s parallel analysis, and robustness is assessed with leave-one-out principal angles. Of 22 descriptors, nine were near-constant or missing and were excluded. Variance concentrated in two components: presence 55.7% and 28.4%; relative 75.8% and 15.5%; both PCs exceeded the 95th percentile noise envelope. The PC1 and PC2 subspace was stable (mean cosine 0.993 to 0.994; largest deletion about 22 degrees), and A1 and A2 produced nearly identical subspaces. PC1 reflected hydrophobicity versus polarity and hydrogen bonding (higher LogP/LogD, lower H-bond acceptors/donors and polar surface area). PC2 captured adsorption and bioconcentration together with molecular flexibility and optical proxies (higher KOC and BCF, more freely rotating bonds, higher refractive index). Relative weighting is recommended as primary; presence serves as a concordant robustness check. Key Words: Horn’s parallel analysis; principal angles; leave-one-out stability; Dirichlet pooling; PCA. 1. Introduction Recent ecological research has increasingly recognized that traditional functional traits, such as root diameter, specific root length (SRL), tissue density, and nitrogen content, can only explain part of the plant’s resource acquisition and utilization strategies. In addition to these morphological and physiological traits, plants' chemical defense and longevity strategies also represent crucial aspects of their ecological function. However, traditional morphological and physiological traits often fail to reflect the diversity of metabolite compositions and chemical characteristics. As a result, integrating metabolomic data into plant functional trait frameworks has become an essential focus of recent functional ecology research. In their study published in Science Advances (Walker et al., 2023), the authors systematically revealed the independent role of leaf metabolic traits in the functional trait space across the globe for the first time. They integrated leaf metabolomic data from approximately 800 plant species from tropical and temperate regions, acquiring chemical signals through untargeted LC-MS (Liquid Chromatography-Mass Spectrometry) technology. This approach provided new insights into plant form and function that were previously unexplored through traditional morphological and physiological traits alone. These considerations motivate a statistically principled aggregation from peaks to species-level descriptors and a robust dimensioning and stability assessment. In LC-MS analysis, peak intensity refers to the signal strength of specific peaks in the chromatogram, typically quantified as peak height or area. It reflects the relative or absolute content of the target compounds (analytes) in the sample and is a core indicator for quantitative analysis. The LC-MS process separates compounds through liquid chromatography and then detects ion signals through mass spectrometry. In the chromatogram, the x-axis represents the retention time, and the y-axis represents signal intensity. Each peak corresponds to the elution signal of a particular compound. The signal intensity originates from the ion detector's response (e.g., an electron multiplier), with units usually expressed in counts or arbitrary units (AU). Two types of intensity are commonly used: peak height (maximum signal value), which is often employed for quick comparisons but is influenced by peak width, and peak area (integrated area under the peak), which more accurately represents the total amount of the compound and is frequently used for quantitative analysis. The quantitative meaning lies in the fact that peak intensity is proportional to the concentration of a compound in the sample, and within the linear range, intensity can be converted to absolute concentration through a calibration curve (e.g., ng/mL). For example, in the Extracted Ion Chromatogram (EIC), the total intensity or base peak intensity within a specific m/z window corresponds directly to the abundance of the analyte. This method can also be used to compare differences between samples, such as evaluating metabolite changes through the peak-pair intensity ratio. However, it is generally accepted that LC-MS peak intensity cannot be directly used for crosscompound quantitative comparisons. Statistically, peak values are only meaningful for “relative changes of the same compound across different samples on the same platform/method”; crosscompound comparisons regarding who is “more” or “less” do not have practical significance. The core reason that intensity lacks “quantitative meaning” is twofold: first, the response factors of different compounds vary greatly, so two compounds with the same molar amount may produce vastly different peak intensities on the chromatogram; therefore, “higher peaks” do not necessarily mean “more” of the compound; second, a high peak may simply reflect the compound's chemical/physical properties that make it more easily detected, not necessarily indicating a higher concentration or greater ecological significance. Ecological studies are concerned with how chemical properties (e.g., aromaticity, polarity, size, hydrogen bonding ability) combine to form axes like “chemical defense” or “leaf longevity,” rather than focusing on raw peak intensities of named metabolites. Due to the lack of strict quantitative comparability of peak intensities across compounds, Walker et al. (2023) chose to convert peak intensity data into a presence/absence matrix, retaining only information about whether a compound was detected in the samples. This approach avoids the systematic bias introduced by cross-compound response differences and ensures that the analysis reflects true structural differences in chemical characteristics rather than technical noise from signal intensities. In the annotation and descriptor extraction stages, the authors used the GNPS molecular network platform and LOTUS chemical database for in silico annotation of metabolites, and based on the Chemistry Development Kit (CDK), they extracted 21 chemical descriptors, including molecular weight, number of aromatic atoms, polarity, logP, hydrogen bond donors/acceptors, sp³:sp² ratio, representing key chemical dimensions such as structural complexity, reactivity, polarity, hydrophobicity, and carbon-bond saturation. Subsequently, the authors selected five representative descriptors (molecular weight, number of aromatic atoms, HBA, logP, sp³: sp² ratio) for principal component analysis (PCA), which allowed them to identify the main gradients of metabolic traits within the chemical descriptor space. By integrating these chemical characteristics, Walker et al. (2023) identified two major functional axes for leaf metabolic traits: the chemical defense spectrum, which correlates with the compound’s aromaticity, bond saturation, and polarity, representing the reactivity and detoxification potential of metabolites; and the leaf longevity axis, which is associated with molecular weight and hydrogen bond acceptors, reflecting the size, stability, and intermolecular interactions of compounds. These two chemical functional axes were found to be highly consistent across tropical and temperate plant communities, suggesting that plant metabolic chemistry follows universal ecological patterns globally. Further analysis revealed that these two axes were nearly orthogonal to the eight traditional morphological and physiological traits (e.g., specific leaf area (SLA), leaf thickness, nitrogen content), indicating that metabolic traits uncover hidden functional dimensions that traditional traits cannot explain. In other words, plant metabolic chemistry not only complements the traditional functional trait framework but also provides a new understanding of plant adaptive strategies and ecological functions. This study has methodological significance as well. First, it addresses the issue of cross-compound comparability by converting peak intensities into presence/absence information, allowing for more accurate chemical trait analysis. Second, it introduces chemical informatics indices that quantify the chemical descriptors of individual metabolites, bridging metabolomic data with ecological functions. Finally, by employing multivariate statistical analysis (PCA), the authors integrate complex chemical dimensions to construct interpretable metabolic functional gradients, providing a new methodological framework for functional trait studies. Recent work shows that metabolite-level chemical traits form axes largely orthogonal to classical functional traits (Walker et al., 2023). Moving from peak intensities to structure-informed descriptors follows established metabolomics workflows (Wang et al., 2016). Building on this concept, the current study further explores the network structure of metabolites, incorporating cooccurrence relationships and the characteristics of isolated metabolites, to investigate metabolic specialization patterns across species and habitats. This network perspective, complementing Walker et al.'s chemical trait dimensions, reveals metabolite co-variation, specialization, and ecological significance, thereby providing a more comprehensive understanding of the chemical dimensions of plant functional traits. 2. Data Preparation and Methodology 2.1 Data We analyzed root-exudate metabolomes and integrated two tabular sources from untargeted LC– MS. The metabolite table contained 219 compounds (columns include id, Peak (compound name), Similarity, R.T. (in minutes), Mass, Count) and sample columns for nine species (Tea Tree, Nephelium lappaceum tree, Syzygium samarangense tree, Phoebe puwenensis, Cornus walteri, Radermachera hainanensis, Phillyrea angustifolia, Tectona grandis tree, Pomelo Tree), each with six biological replicates. The descriptor table listed candidate chemical structures per compound id, their ChemSpiderName, and 22 physicochemical descriptors including Density, Boiling Point, Vapour Pressure, Enthalpy of Vaporization, Flash Point, Index of Refraction, Molar Refractivity, #H bond acceptors (HBA), #H bond donors (HBD), #Freely Rotating Bonds (FRB), #Rule of 5 Violations, ACD/LogP, ACD/LogD (pH 5.5), ACD/BCF (pH 5.5), ACD/KOC (pH 5.5), ACD/LogD (pH 7.4), ACD/BCF (pH 7.4), ACD/KOC (pH 7.4), Polar Surface Area (PSA), Polarizability, Surface Tension, Molar Volume (Yap, 2011). 2.2 Methodology Because a metabolite id can map to multiple candidate structures (isomers), we first integrated isomer uncertainty at the descriptor vector 𝒄𝒄𝒊𝒊𝒊𝒊 ∈ℝ𝟐𝟐𝟐𝟐 and library similarity scores 𝒔𝒔𝒊𝒊𝒊𝒊. We formed weights 𝒘𝒘𝒊𝒊𝒊𝒊 and computed an isomer-integrated descriptor 𝒄𝒄 �𝒊𝒊=�𝒘𝒘𝒊𝒊𝒊𝒊𝒄𝒄𝒊𝒊𝒊𝒊 𝑲𝑲𝒊𝒊 𝒊𝒊=𝟏𝟏 , �𝒘𝒘𝒊𝒊𝒊𝒊 𝒊𝒊 =𝟏𝟏. Two deterministic priors were evaluated. In A1(Symmetric Dirichlet), 𝒘𝒘𝒊𝒊~𝑫𝑫𝒊𝒊𝑫𝑫(𝜶𝜶, … , 𝜶𝜶) with 𝜶𝜶= 𝟎𝟎.𝟏𝟏. In A2 (similarity-informed Dirichlet), we defined 𝒔𝒔 �𝒊𝒊=𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔(𝒔𝒔𝒊𝒊/𝝉𝝉 ) with 𝝉𝝉=𝟎𝟎.𝟐𝟐 and set 𝒘𝒘𝒊𝒊∝ 𝜶𝜶𝟎𝟎𝟏𝟏+𝝀𝝀𝒔𝒔 �𝒊𝒊 with 𝜶𝜶𝟎𝟎=𝟎𝟎.𝟏𝟏, 𝝀𝝀=𝟏𝟏𝟎𝟎. We used the normalized weighted vectors (posterior means) without Monte Carlo sampling, yielding fully reproducible results. For descriptor-specific missingness among candidates, weights were renormalized over observed entries when computing 𝒄𝒄 �𝒊𝒊(𝒋𝒋). Different metabolite IDs can refer to the same chemical (ChemSpiderName) due to synonymy. To avoid double counting identical structures, we handled synonymy at the species-level weighting step by summing weights of IDs that share the same ChemSpiderName within each species before computing community-weighted means. Species weights were constructed to respect the non-comparability of absolute peak areas across compounds. We used two complementary schemes. Presence weights captured detection frequency within a species: for the six replicates 𝒍𝒍 and species 𝒔𝒔, let 𝑰𝑰𝒊𝒊,𝒔𝒔,𝒍𝒍=𝟏𝟏 if the cell is nonempty (zeros allowed) and 𝟎𝟎 otherwise; then 𝑝𝑝𝑖𝑖,𝑠𝑠=1 6�𝐼𝐼𝑖𝑖,𝑠𝑠,𝑙𝑙 ∈[0,1] 6 𝑙𝑙=1 . Relative weights captured within-compound, across-species partitioning that is invariant to crosscompound scaling. We computed median log-abundance 𝛼𝛼𝑖𝑖,𝑠𝑠=𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑛𝑛𝑙𝑙{log (𝑚𝑚𝑎𝑎𝑚𝑚𝑚𝑚𝑖𝑖,𝑠𝑠,𝑙𝑙+𝜀𝜀)} and applied a softmax across species, 𝑎𝑎𝑖𝑖,𝑠𝑠=exp�𝑚𝑚𝑖𝑖,𝑠𝑠� ∑exp (𝑚𝑚𝑖𝑖,𝑠𝑠′) 𝑠𝑠′ , so that ∑𝑎𝑎𝑖𝑖,𝑠𝑠𝑠𝑠 = 1. Combining A1/A2 with presence/relative yielded four configurations: presence_A1, presence_A2, relative_A1 and relative_A2. 2.2.1 Species-level community-weighted means (CWM) For each configuration we derived species-level community-weighted means (CWM) of descriptors. For species 𝑠𝑠 and descriptor 𝑓𝑓, with 𝑊𝑊𝑖𝑖,𝑠𝑠∈{𝑝𝑝𝑖𝑖,𝑠𝑠,𝑎𝑎𝑖𝑖,𝑠𝑠} and isomer-integrated 𝑐𝑐𝑖𝑖, we first merge synonymy by summing 𝑊𝑊𝑖𝑖,𝑠𝑠 across IDs sharing the same ChemSpiderName. We then normalized weights within species to ensure convex aggregation (see PCA conventions in Jolliffe & Cadima, 2016), 𝑊𝑊 �𝑖𝑖,𝑠𝑠=𝑊𝑊𝑖𝑖,𝑠𝑠 ∑𝑊𝑊𝑖𝑖′,𝑠𝑠𝑖𝑖′ , and computed the CWM 𝑀𝑀𝑠𝑠,𝑓𝑓=�𝑊𝑊 �𝑖𝑖,𝑠𝑠𝑐𝑐𝑖𝑖(𝑓𝑓). 𝑖𝑖 This produced a 9 × 22 species-by-descriptor matrix 𝑀𝑀 per configuration. Columns of 𝑀𝑀 were standardized (z-scores) across species to yield 𝑍𝑍, and descriptors with zero between-species variance (or all missing) were removed. Nine descriptors—Density, Boiling Point, Vapour Pressure, Enthalpy of Vaporization, Flash Point, Molar Refractivity, Polarizability, Surface Tension, and Molar Volume—were consistently dropped as near-constants at the species level; this filtering prevents numerical instability without affecting the chemical interpretation of retained axes. We applied principal component analysis (PCA; full SVD) to 𝑍𝑍 after column standardization. Let 𝑃𝑃 denote the loading matrix (descriptors × PCs) and 𝑇𝑇 the score matrix (species × PCs); the decomposition satisfies 𝑍𝑍 ≈𝑇𝑇𝑃𝑃𝑇𝑇, 𝑃𝑃𝑇𝑇𝑃𝑃=𝐼𝐼. We interpreted the leading components and exported both scores and loadings for downstream visualization (biplots, correlation circles, squared-loading bar plots). The number of components to retain was determined by Horn’s parallel analysis (Horn, 1965). For each configuration we generated 𝑛𝑛𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 =1000 null matrices with the same dimension as 𝑍𝑍, zscored columns, and obtained empirical eigenvalue percentiles 𝑞𝑞95,𝑘𝑘. All observed components with 𝜆𝜆 󰆹𝑘𝑘>𝑞𝑞95,𝑘𝑘 were retained; cumulative variance for PC1–PC2 (and PC3 when applicable) was reported alongside the decision. Robustness was assessed in two ways. First, we performed leave-one-species-out (LOO) PCA: for each species 𝑢𝑢, we recomputed PCA on 𝑍𝑍∖{𝑢𝑢} and quantified stability of the PC1–PC2 loading subspace via principal angles {𝜃𝜃1,𝜃𝜃2} between the subspaces spanned with and without 𝑢𝑢; results are summarized by cos 𝜃𝜃 and angles in degrees (Björck & Golub, 1973). Second, we evaluated sensitivity to modeling choices by computing principal angles between PC1–PC2 loading subspaces across configurations (A1 vs A2; presence vs relative) and by correlating species scores (PC1/PC2) pairwise across configurations; near-identity indicates invariance of axis definitions and species ordering to the isomer prior and weighting scheme. All computations were implemented in Python (pandas, numpy, scikit-learn, scipy, matplotlib). Parameter settings were A1: α = 0.1; A2: α₀ = 0.1, λ = 10, τ = 0.2; Horn permutations = 1,000; random state = 42. Posterior means (no Monte Carlo) were used for isomer weights to guarantee determinism. 3. Results In this study, we applied principal component analysis (PCA) to explore the chemical diversity across species based on both presence/absence and relative abundance of compounds. The results from both analytical approaches (A1 and A2) were highly similar, suggesting that the relative abundance and presence-based models provide the most robust and interpretable insights. Consequently, the following results are presented based on these two models. To assess the impact of different modeling choices, we compared the principal angles between the loading subspaces derived from different configurations (e.g., presence vs. relative and A1 vs. A2). The principal angles between the subspaces were consistently small, with a mean cosine similarity of 0.99, indicating that the chemical axes were highly consistent regardless of the method used for weighting or isomer prior. This confirms that the identified functional axes are robust to different configurations. Scree and Horn jointly support two components in all configurations; we therefore interpret PC1 and PC2 throughout. 3.1 Variance Explained by Principal Components The first principal component (PC1) consistently explained the majority of the variance in the data across all models, regardless of whether presence or relative abundance was considered. The cumulative explained variance for the first few components reinforces the importance of PC1 and PC 2 in capturing the underlying structure of the data. Table 1: Variance Explained by Principal Components (PC1, PC2, PC3) in Presence and Relative Abundance Models Model PC1 Explained Variance PC2 Explained Variance PC3 Explained Variance Presence_A1 55.71% 28.38% 11.77% Presence_A2 55.71% 28.38% 11.77% Relative_A1 75.81% 15.54% 4.89% Relative_A2 75.81% 15.54% 4.89% The data indicates that PC1 explains the greatest proportion of variance across both the presence and relative abundance models. Relative abundance-based models (e.g., Relative_A1 and Relative_A2) exhibit higher explained variance for PC1, suggesting that accounting for compound concentrations offers a more detailed understanding of the chemical variation across species. 3.2 Significance of Principal Components Based on Horn Parallel Analysis Further validation of the significance of PC1–PC2 was provided by Horn parallel analysis, which compares observed eigenvalues to those derived from random permutations. Scree profiles together with Horn’s parallel analysis supported retaining two axes: the empirical PC1–PC2 variances lay above the 95% noise envelope, so keeping PCs whose explained variance exceeds the q95 threshold is warranted. In practice we interpret PC1 and PC2 throughout. Figure 1: Scree plot and Horn’s parallel analysis for component retention. The scree plot shows the variance explained by each principal component, with the 95% noise threshold (Horn’s parallel analysis) used to determine the number of components to retain. PC1 and PC2 exceed the threshold in all configurations, justifying their retention for subsequent analysis. 3.3 Descriptors Driving Principal Components Nine chemical attributes were excluded in all four settings due to zero variance or insufficient information in species-level CWMs, and thus contributed nothing to PCA: Density, Boiling point, Vapour pressure, Enthalpy of vaporization, Flash point, Molar refractivity, Polarizability, Surface tension, and Molar volume. Their near-constancy or absence at the CWM stage explains the lack of signal. Leading contributors (top six absolute loadings, signs kept) were nearly identical between A1 and A2 and highly concordant between presence and relative weightings: Table 2: PCA Loadings for Key Chemical Descriptors Influencing PC1 in Presence_A1/A2 Model Chemical Descriptor Contribution to PC1 ACD/LogD (pH 5.5) +0.358 ACD/LogP +0.353 #H-bond acceptors –0.353 ACD/LogD (pH 7.4) +0.349 #H-bond donors –0.348 Polar Surface Area –0.348 Table 3: PCA Loadings for Key Chemical Descriptors Influencing PC2 in Presence_A1/A2 Model Chemical Descriptor Contribution to PC2 ACD/BCF (pH 7.4) +0.453 ACD/BCF (pH 5.5) +0.430 ACD/KOC (pH 5.5) +0.406 ACD/KOC (pH 7.4) +0.384 #Freely rotating bonds +0.261 Index of refraction +0.258 Table 4: PCA Loadings for Key Chemical Descriptors Influencing PC1 in Relative_A1/A2 Model Chemical Descriptor Contribution to PC1 ACD/LogP +0.315 ACD/LogD (pH 5.5) +0.311 ACD/LogD (pH 7.4) +0.296 #H-bond acceptors –0.293 Polar Surface Area –0.292 #H-bond donors –0.290 Table 5: PCA Loadings for Key Chemical Descriptors Influencing PC2 in Relative_A1/A2 Model Chemical Descriptor Contribution to PC2 ACD/KOC (pH 5.5) +0.384 #Rule of 5 violations +0.376 #Freely rotating bonds +0.364 Index of refraction +0.317 ACD/KOC (pH 7.4) +0.307 #H-bond donors +0.285 Axis meanings were stable across all four runs. PC1 captured a hydrophobicity–polarity trade-off: higher ACD/LogP and LogD (pH 5.5/7.4) contrasted with lower hydrogen-bonding capacity and polarity (#H-bond acceptors, #H-bond donors, Polar Surface Area). This is the expected hydrophobic partitioning vs. H-bond/polar surface axis. PC2 aggregated environmental adsorption/bioconcentration and molecular flexibility/optical density proxies, loading positively on KOC (pH 5.5/7.4) and BCF (pH 5.5/7.4) together with #Freely rotating bonds and Index of refraction—an axis of “organophilic/biophilic tendency + conformational flexibility.” 3.4 Stability of PCA Results: Leave-One-Out (LOO) Analysis The Leave-One-Out (LOO) stability analysis confirms the robustness of the PCA results. Cosine similarity values close to 1 and small angular differences indicate that the removal of any individual species does not significantly alter the PCA results, further validating the stability and generalizability of the identified chemical diversity patterns. Table 6: LOO Stability for Key Species in Presence and Relative Abundance Models Species Cosine Similarity (PC1) Cosine Similarity (PC2) Principal angle 1 Principal angle 2 Tea Tree 0.999934 0.999537 0.66° 1.74° Nephelium lappaceum tree 0.999775 0.997077 1.22° 4.38° Syzygium samarangense tree 0.998541 0.985718 3.10° 9.70° Phoebe puwenensis 0.999199 0.991288 2.29° 7.57° Cornus walteri 0.999402 0.994312 1.98° 6.11° Radermachera hainanensis 0.999266 0.998193 2 .19° 3.44 ° Phillyrea angustifolia 0.999762 0.995144 1.25° 5.64° Tectona grandis tree 0.996915 0.927530 4.50° 21.94° Pomelo Tree 0.999454 0.996821 1.89° 4.57° Model concordance and robustness were high. Subspace comparisons showed A1 and A2 were indistinguishable (PC-wise cosθ≈1.0), while presence vs. relative differed only by a modest rotation (PC1 cosθ≈0.988; PC2 cosθ≈0.959), indicating the same ecological chemistry with slightly different emphasis. Leave-one-out (LOO) analysis of the PC1–PC2 plane confirmed strong resilience to sample deletion. For presence, the mean cosine similarity was 0.993 (mean angle ≈ 4.68°), with a maximum angle ≈ 21.95° when omitting Tectona grandis tree; for relative, the mean cosine was 0.994 (mean angle ≈ 3.78°), maximum ≈ 22.46° (Tectona grandis tree). Other relatively “sensitive” leave-one-outs were typically Syzygium samarangense tree and Phoebe puwenensis, but with angles well below Tectona grandis tree. The interpretation is straightforward: the PC1–PC2 subspace is robust, Tectona grandis tree carries a pronounced chemical signal, and changing the isomer prior (A1 or A2) does not move the subspace. 3.5 PCA Dimension Decision The PCA dimension decision table is shown as below: Table 7: PCA Dimension Decision Based on Explained Variance Model n_vars_used Keep PC by Horn95 Cumulative PC1_PC2 Cumulative PC1_PC2_PC3 Presence_A1 13 2 0.8409 0.9586 Presence_A2 13 2 0.8409 0.9586 Relative_A1 13 2 0.9135 0.9624 Relative_A2 13 2 0.9135 0.9624 These results suggest that retaining PC1 and PC2 is sufficient to explain the majority of the variance across species, with little additional contribution from higher-order components.