Multivariate Independent Component Analysis Identifies Patients in Newborn Screening Equally to Adjusted Reference Ranges

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Multivariate Independent Component Analysis Identifies Patients in Newborn Screening Equally to Adjusted Reference Ranges © 2023 the Authors Published version Kouři,l Štěpán; de Sousa, Julie; Fačevicová Kamila; Gardlo, Alžběta; Muehlmann, Christoph; Nordhausen, Klaus; Friedecký, David; Adam, Tomáš Kouři,l Štěpán, de Sousa, Julie, Fačevicová Kamila, Gardlo, Alžběta, Muehlmann, Christoph, Nordhausen, Klaus, Friedecký, David, Adam, Tomáš. (2023). Multivariate Independent Component Analysis Identifies Patients in Newborn Screening Equally to Adjusted Reference Ranges. International Journal of Neonatal Screening, 9(4), Article 60. https://doi.org/10.3390/ijns9040060 2023 Citation: Kouˇril, Š.; de Sousa, J.; Faˇcevicová, K.; Gardlo, A.; Muehlmann, C.; Nordhausen, K.; Friedecký, D.; Adam, T. Multivariate Independent Component Analysis Identifies Patients in Newborn Screening Equally to Adjusted Reference Ranges. Int. J. Neonatal Screen. 2023,9, 60. https://doi.org/ 10.3390/ijns9040060 Academic Editor: Ralph Fingerhut Received: 22 August 2023 Revised: 22 September 2023 Accepted: 17 October 2023 Published: 20 October 2023 Copyright: © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). International Journal of Neonatal Screening Article Multivariate Independent Component Analysis Identifies Patients in Newborn Screening Equally to Adjusted Reference Ranges Štˇepán Kouˇril 1, Julie de Sousa 2,3, Kamila Faˇcevicová3, Alžbˇeta Gardlo 1,2, Christoph Muehlmann 4, Klaus Nordhausen 5, David Friedecký1and Tomáš Adam 1,2,6,* 1Department of Clinical Biochemistry, University Hospital Olomouc, 779 00 Olomouc, Czech Republic; [email protected] (D.F.) 2Laboratory of Metabolomics, Institute of Molecular and Translational Medicine, PalackýUniversity Olomouc, 779 00 Olomouc, Czech Republic 3Department of Mathematical Analysis and Applications of Mathematics, Faculty of Science, Palacký University Olomouc, 779 00 Olomouc, Czech Republic; [email protected] 4Institute of Statistics & Mathematical Methods in Economics, Vienna University of Technology, 1040 Vienna, Austria 5Department of Mathematics and Statistics, University of Jyväskylä, 40014 Jyväskylä, Finland 6Faculty of Health Care, The Slovak Medical University in Bratislava, 974 05 BanskáBystrica, Slovakia *Correspondence: [email protected]; Tel.: +42-0736-737-371 Abstract: Newborn screening (NBS) of inborn errors of metabolism (IEMs) is based on the reference ranges established on a healthy newborn population using quantile statistics of molar concentrations of biomarkers and their ratios. The aim of this paper is to investigate whether multivariate independent component analysis (ICA) is a useful tool for the analysis of NBS data, and also to address the structure of the calculated ICA scores. NBS data were obtained from a routine NBS program performed between 2013 and 2022. ICA was tested on 10,213/150 free-diseased controls and 77/20 patients (9/3 different IEMs) in the discovery/validation phases, respectively. The same model computed during the discovery phase was used in the validation phase to confirm its validity. The plots of ICA scores were constructed, and the results were evaluated based on 5sd levels. Patient samples from 7/3 different diseases were clearly identified as 5sd-outlying from control groups in both phases of the study. Two IEMs containing only one patient each were separated at the 3sd level in the discovery phase. Moreover, in one latent variable, the effect of neonatal birth weight was evident. The results strongly suggest that ICA, together with an interpretation derived from values of the “average member of the score structure”, is generally applicable and has the potential to be included in the decision process in the NBS program. Keywords: newborn screening; independent component analysis; mass spectrometry; multivariate statistical analysis; inborn errors of metabolism; compositional data analysis 1. Introduction Inborn errors of metabolism (IEMs) are caused by an enzyme deficiency that usually leads to the accumulation of substrates of the defective enzyme. Nearly 1900 diseases are currently classified in this group. Newborn screening (NBS) is an active and widespread search for diseases in their early, preclinical stages so that these diseases can be diagnosed and treated before they become manifest and cause irreversible damage to a child’s health. In the screening of IEMs, flow injection–tandem mass spectrometry analysis is used to quantify specific analytes (amino acids, acylcarnitines) which are diagnostically relevant to the diseases in question [ 1 ]. The data evaluation is performed using reference ranges established on a healthy newborn population using quantile statistics of untransformed data. The key diagnostic biomarkers of IEMs are elevated levels of substrates of defective Int. J. Neonatal Screen. 2023,9, 60. https://doi.org/10.3390/ijns9040060 https://www.mdpi.com/journal/ijns Int. J. Neonatal Screen. 2023,9, 60 2 of 14 enzymes. The first introduced screening program detected phenylketonuria (PKU) patients using a bacterial inhibition assay for (semi)quantitative analysis of phenylalanine (Phe) in dried blood spots [ 2 ]. Besides reference ranges for individual diagnostic metabolites, their mutual ratios [ 3 ] are applied as a post-analytical tool. The cutoff target ranges of analytes and ratios are then defined as selected percentiles of the control and diseased populations. When overlaps occur, adjustments are made to optimize sensitivity and specificity considering Wilson and Jungner criteria [4]. The first application of amino acid relations used in diagnostics was the phenylalanine/tyrosine (Phe/Tyr) ratio (substrate and product of the defective enzymatic conversion, respectively) for the detection of PKU heterozygous carriers [ 5 ]. This approach was successfully applied to NBS, lowering substantially (to 25%) false positives [ 6 , 7 ] compared to evaluation based on the absolute concentration of Phe alone. Since its introduction, several ratios were suggested for screening of medium-chain acyl-CoA dehydrogenase deficiency (MCAD) [ 8 ] and/or carnitine palmitoyltransferase II deficiency (CPTII) [ 9 , 10 ]. Similarly, ratios for screening pyruvate dehydrogenase complex deficiencies and other mitochondrial disorders associated with lactic acidosis and variably elevated alanine (Ala) and proline (Pro) levels have recently been proposed [ 11 ]. The authors used strictly ketogenic amino acid leucine (Leu) and lysine (Lys), not involved in the glycolytic pathway generating pyruvate, for “normalizing metabolites in quantitative analysis of Ala and Pro” [ 11 ] and forming diagnostically effective Ala/Leu and Pro/Leu ratios. In the field of IEMs, the developed diagnostic ratios generally correct for sampling variability (including the hematocrit effect) and elucidate substrate/product relations (e.g., proximal Phe/Tyr ratio; distant in case of linear enzymatic pathways, e.g., octanoylcarnitine (C8) and acetylcarnitine (C2) acylcarnitines). In the Czech Republic, the NBS program targets 15 IEMs listed in Table 1. Table 1. List of diseases measured using mass spectrometry included in newborn screening in the Czech Republic. Disease * Primary Biomarkers Secondary Biomarkers Number of Patients (Discovery/Validation Set) PKU Phe, Phe/Tyr - 54/16 MSUD Xle, Xle/Ala, (Xle + Val)/Pro + Tyr) Val 7/0 MCAD C8, C8/C2 C10, C10:1, C6, C8/C10 7/3 LCHAD C16-OH, C18:1-OH C18-OH 4/0 VLCAD C14:1, C14:1/C16 C14 1/0 CPT I C0, C0/(C16 + C18) C18, C18:1, C16 - CPT II/CACT ** C16, (C16 + C18:1)/C2 C18, C18:1, C0 - GA I C5DC, C5DC/C16 C5DC/C8 1/0 IVA C5, C5/C8 C5/C2 1/0 HCY(CBS) Met, Met/Phe - - HCY(MTHFR) Met, Met/Phe, - 1/0 ARG Arg, Arg/Orn, Arg/Phe - - CIT/ASA ** Cit, Cit/Phe, Orn/Cit, ArgSucc - 1/1 * For a complete list of OMIM numbers, names of diseases, and their incidence see Table S1. Full names of metabolites: Phenylalanine (Phe), Tyrosine (Tyr), Leucine/Isoleucine (Xle), Alanine (Ala), Valine (Val), Proline (Pro), Octanoylcarnitine (C8), Acetylcarnitine (C2), Decanoylcarnitine (C10), Decenoylcarnitine (C10:1), Hexanoylcarnitine (C6), 3Hydroxypalmitoylcarnitine (C16-OH), 3-Hydroxyoleoylcarnitine (C18:1-OH), 3-Hydroxystearoylcarnitine (C18-OH), Tetradecenoylcarnitine (C14:1), Palmitoylcarnitine (C16), Tetradecanoylcarnitine (C14), Carnitine free (C0), Stearoylcarnitine (C18), Oleoylcarnitine (C18:1), Glutarylcarnitine (C5DC), Isovalerylcarnitine/Methylbutyrylcarnitine (C5), Methionine (Met), Arginine (Arg), Ornithine (Orn), Citrulline (Cit), Argininosuccinate (ArgSucc). ** two distinctive metabolic diseases indistinguishable by screening biomarkers. In general, it is expected that the metabolic findings in patients are atypical compared to the healthy population. For this reason, the data from patients can be understood as outlying observations. The above-described routine approach to NBS, based on reference ranges of individual metabolites or their mutual ratios, uses only marginal information from the metabolic profile, which technically represents a multivariate observation. Therefore, it essentially provides an incomplete picture, and univariate outlier detection methods are missing (for example, outliers where the interplay between variables is atypical). Classical Int. J. Neonatal Screen. 2023,9, 60 3 of 14 multivariate outlier detection tools are usually based on the Mahalanobis distance which, however, scales badly with increasing dimension [ 12 ]. Hence, given the dimensionality of the NBS data, dimension reduction prior to outlier detection seems a natural way to proceed. On the other hand, traditionally used unsupervised dimension reduction methods like principal component analysis (PCA) might not always be suitable as, for example, outliers do not necessarily need to be found in the direction of high variations, representing the main aim of the method. The aim of this study is the application of independent component analysis (ICA) on metabolic NBS data, as a multivariate alternative to the traditional methods based on univariate reference ranges established on prior knowledge of the disease biomarker behavior and percentile statistics. The structure of the ICA scores with the aim of investigating their potential to elucidate metabolic relations of the probands is addressed and a comparison with the traditional methods is provided. We applied ICA on data obtained in a routine newborn screening program performed between 2013 and 2022 in the Czech metabolic screening center in Olomouc. The reason for choosing ICA is that it maximizes non-Gaussianity in the data while searching for independent components, which is beneficial in our context where atypical observations are usually acknowledged regarding the Gaussian distribution. Therefore, we suggest considering the mean ± 5sd rule for the computed independent components as an alternative to the univariate reference ranges of the original variables (i.e., metabolites). Moreover, ICA could potentially reveal other unknown important metabolite ratios theoretically usable in classical newborn screening. 2. Materials and Methods 2.1. Patients and Samples The anonymized samples were processed as separate discovery (years 2011–2020, 10,213 disease-free controls and 77 patients suffering from nine different diseases) and validation (years 2021–2022, 150 controls and 20 patients; see Table 1) studies. Blood samples were collected from newborns to the screening cards (Whatman 903) 48–72 h after birth and transported to the laboratory. Samples were prepared according to the CE-IVD kit from Chromsystems (order no. 57000). Discs (3 mm diameter) were punched out from the cards and placed in the 96-well microtiter plate. The extraction buffer (100 µ L) with internal standard was added to each sample. The plates were covered with foil, shaken (600 rpm) for 20 min at laboratory temperature, and then centrifuged (10 min at 2000 rpm). The supernatant was used for the analysis. 2.2. DBS Analysis by Mass Spectrometry Mass spectrometric analysis was performed using LC-MS API 4000 Ultimate 3000 RS (Sciex). Amino acids and acylcarnitines were determined in dry blood spots using the above-mentioned CE-IVD kit. The analytical method is routinely used in the laboratory for screening more than 30,000 newborns per year and is accredited according to ISO 15,189 and participates in ERNDIM and Newborn Screening Quality Assurance Program (NSQAP) quality control schemes. Concentrations of 26 amino acids and acylcarnitines used for screening purposes in the Czech Republic were determined and used for further analyses (Table S2). Data were centered and normalized by log-ratio (clr) transformation without any additional pre-processing steps. 2.3. Data Analysis—ICA The statistical processing of the data is based on ICA and was done with the help of the software R [ 13 ] using fICA [ 14 ] and robCompositions [ 15 ] packages. ICA in general looks for a set of latent variables zi , which have the form of linear combinations of the measured variables (metabolites) xk . The main property of the latent vector z is that its components are standardized, independent, and likely to follow a non-Gaussian distribution looking as non-Gaussian as possible. Therefore, the method can typically reveal complex sources Int. J. Neonatal Screen. 2023,9, 60 4 of 14 of outlyingness or groupings within the data [ 16 ]. Formally written, the independent component model assumes that the observable p-variate vector x is linked with the platent components as Az+b , where A is a nonsingular matrix (i.e., square matrix with non-zero determinant) and b is a location vector, which can be understood as a vector of means of x . The method, in consequence, searches for a matrix W, called unmixing matrix, which would identify the latent structures as W(x−b) . The resulting vector z is then given uniquely up to the sign and order of its components. There are several strategies on how to estimate the unmixing matrix W ; some of them are listed in [ 17 ] together with a more detailed description of the method. Our case study is based on the adaptive deflation-based FastICA approach [ 18 ] which maximizes the non-Gaussianity of the latent components using the most suitable non-Gaussianity measure for each component separately. The first step of almost all ICA methods is whitening where the data are standardized and uncorrelated. For that purpose, the covariance matrix of the data needs to be full rank, i.e., nonsingular. Similar to the approach used in [ 17 ], we treat the metabolomic data as relative-valued and the final algorithm thus combines ICA with the principles of compositional data analysis [ 19 ]. So-called centered log-ratio (clr) coefficients give the natural representation of relative-valued (compositional) data. Within this representation, each of the original variables xkcorresponds to one coefficient clr(xk)of the form ln(xk/g(x)), with g(x) denoting the geometric mean of the whole vector x . This results in a favorable interpretation in the sense of relative dominance of the given part xk over the whole composition (metabolic profile). On the other hand, the construction of the clr coefficients also implies its zero-sum property, which prevents its direct use in the ICA algorithm as it results in a singular data matrix. An alternative representation is given by the family of isometric log-ratio (ilr) coordinates. This representation characterizes the compositional vector x by a system of p− 1 orthonormal (i.e., orthogonal, with a unit norm) real-valued coordinates, given as ilr(x)=VTln(x) , with V being a p×(p−1) matrix, called contrast matrix. This coordinate system overcomes the problem of singularity and is, therefore, popularly used in a wide range of multivariate statistical methods, including those relying on the full rank assumption. A detailed description of the construction and properties of this coordinate system is given (e.g., in [ 19 ]); however, let us emphasize that the clr and ilr representations are mutually transferable through the contrast matrix V as ilr(x)=VTclr(x) . Accordingly, the results obtained from the ilr representation can be transformed into the clr space and interpreted there. Note also here that the clr results are invariant to the chosen ilr basis. Following the strategy described in [ 17 ], first, the unmixing matrix Wilr is estimated for the whitened ilr representation of the data (we used the system of pivot coordinates, see [ 20 ] for details) and consequently rotated as Wclr =WilrVT . The rows of Wclr can be then understood as loading vectors, specifying the contribution of the clr coefficients (i.e., the relative dominance of the respective compositional parts) to the overall values of the latent components called scores. A positive loading determines that the increase in the relative dominance of the respective compositional part results in an increase in the score. The increase in the relative dominance of parts with a negative loading corresponds to the decrease in the score. More specifically, the (p−1) -dimensional vector of scores z is for a patient/control sample with measurement x equal to Wclrclr(x)−bclr , with bclr standing for the mean clr vector. For formal decision making, the component-wise rule mean(zi)±κ·sd(zi) , i= 1, . . . , p− 1 is employed to decide if an observation is in the direction of the i-th latent variable (independent component, ICi) appearing as an outlier. We follow here [ 21 ] and use median and median absolute deviation (mad) instead of classical estimates of mean and sd as this avoids the masking effects of outliers. Consequently, the estimates of sd(zi) are derived from mads as 1.4826 ·mad(zi) [ 22 ]. The tuning parameter value κ is used to decide how extreme observations are considered outlying. Based on the training data, this value was set to be 5, as the 5sd rule reflects the low proportion of patients within the sample and, in comparison to the 3 or 4sd rule, yields the best separation of patients from controls. Int. J. Neonatal Screen. 2023,9, 60 5 of 14 3. Results and Discussion It is to be expected that ICA, by its nature, projects patients with individual diseases in separate ICs, as each of the diseases affects the metabolic profile differently. In order to gain better insight into which roles the individual metabolites play in the calculations and to be able to relate them to the routine criteria based on reference ranges of biomarkers or their ratios, we propose looking at the structure of an “average score”. Average score zj i of the i-th IC, computed for the j-th group of patients with the individual disease or the group of controls can be written as: zj i= p ∑ k=1 Wclr [i,k]clr(xk)j−bclr k= p ∑ k=1kaj ICi , where i= 1, . . . , p− 1 is an index of the IC of the interest and jranges over the groups of patients/controls. Wclr [i,k] stands for the entry of the clr unmixing matrix at the position [i,k] ; therefore, it gives the value of loading respective to the i-th IC and the k-th metabolite, k= 1, . . . , p . Finally, clr(xk)j equals the mean clr value of the k-th metabolite computed within the j-th group, while bclr k denotes its overall mean (the k-th entry of the clr mean vector bclr). Even though it is not a general feature of ICA, which is within its estimation phase completely unaware of the grouping (i.e., unsupervised), based on the above formula, we can understand the loadings as weights allowing for emphasizing the differences among the individual studied groups within the given IC. Let us further introduce the concept of “an average member of the score structure” (AMSS). To be able to better describe the sources of (possible) outlyingness of the j-th group of patients/controls in the direction of the i-th latent variable, a system of AMSS values kaj ICi is computed as Wclr [i,k]clr(xk)j−bclr k , k= 1, . . . , p . Each of these values quantifies the contribution of the k-th metabolite to the respective mean score zj i. 3.1. Interpretation of Components In Table S2, there are values of AMSS (including classical sd where applicable; calculated as variance of kaj ICi ) computed for all the studied groups in the first 16 ICs (out of 25), which leads to the best separation of patients/controls. In the following paragraphs, AMSS values for those groups which create separated clusters in the given IC (i.e., a control group vs. one of the patient groups) are sorted and compared. This provides information on which metabolites contribute the most to the separation in terms of the clr coefficients of the measurements. Furthermore, the investigation of AMSS across all the studied groups presents an opportunity for a direct comparison of the effects revealed by ICA to the routinely used approach based on NBS biomarkers. 3.1.1. IC1—MCADD ×VLCAD/LCHAD/IVA/GAI In NBS, patients with MCAD are identified by elevated levels of octanoylcarnitine (C8) and its ratio to acetylcarnitine (C8/C2) as primary markers. Other secondary markers C8/decanoylcarnitine (C10), hexanoylcarnitine (C6), and decenoylcarnitine (C10:1) were suggested; however, these biomarkers are only supportive in the decision-making process [8]. The first latent variable IC1 clearly separates the group of MCAD patients, characterized by its highly negative value (Figure 1A). According to the AMSS values, collected in Table S2, this separation is mainly given by significantly lower C8aMCAD IC1 in comparison to the same value computed for the control group, C8acontrol IC1 . The second primary applied diagnostic biomarker of MCAD is C8/C2 ratio. In our data, values of AMSS for C2 are similar between the patients and healthy controls, C2aMCAD IC1 of − 0.028 vs. C2acontrol IC1 of 0.000, and the sd is very low for both groups (0.012 and 0.010, respectively). In this ratio, C2 Int. J. Neonatal Screen. 2023,9, 60 6 of 14 performs as a “reference/anchoring metabolite” (unaffected by the disease and reflecting the general status of the organism) whose level is stable between controls and patients. Other metabolites and their ratios are used as secondary diagnostic markers and serve mainly to confirm the diagnosis, such as C6, C10, C10:1, and C8/C10 ratio. The differences in the AMSS values of these metabolites are not as high as for the C8 primary marker (see Table S2), but they also affect the separation of the two groups. Contrary to C6 and C8, C10aMCAD IC1 alone shifts the MCAD group more into positive values, thus worsening the separation of the two groups, although patients show elevated concentration levels compared to healthy controls. Int. J. Neonatal Screen. 2023, 9, x FOR PEER REVIEW 6 of 14 C8/decanoylcarnitine (C10), hexanoylcarnitine (C6), and decenoylcarnitine (C10:1) were suggested; however, these biomarkers are only supportive in the decision-making process [8]. The first latent variable IC1 clearly separates the group of MCAD patients, characterized by its highly negative value (Figure 1A). According to the AMSS values, collected in Table S2, this separation is mainly given by significantly lower 𝑎   in comparison to the same value computed for the control group, 𝑎   . The second primary applied diagnostic biomarker of MCAD is C8/C2 ratio. In our data, values of AMSS for C2 are similar between the patients and healthy controls, 𝑎   of −0.028 vs. 𝑎   of 0.000, and the sd is very low for both groups (0.012 and 0.010, respectively). In this ratio, C2 performs as a “reference/anchoring metabolite” (unaffected by the disease and reflecting the general status of the organism) whose level is stable between controls and patients. Other metabolites and their ratios are used as secondary diagnostic markers and serve mainly to confirm the diagnosis, such as C6, C10, C10:1, and C8/C10 ratio. The differences in the AMSS values of these metabolites are not as high as for the C8 primary marker (see Table S2), but they also affect the separation of the two groups. Contrary to C6 and C8, 𝑎   alone shifts the MCAD group more into positive values, thus worsening the separation of the two groups, although patients show elevated concentration levels compared to healthy controls. In the loadings table (included in Table S2), some metabolites with no obvious pathobiochemical connection to MCAD, e.g., tetradecanoylcarnitine (C14) and tetradecenoylcarnitine (C14:1), are also enhanced with a non-zero loading value. It can be observed that the increased loadings correspond to different AMSS values for patients suffering from Very long-chain acyl-CoA dehydrogenase deficiency (VLCAD), i.e., slightly higher 𝑎   and 𝑎  : . The result of this is a minor separation of the VLCAD group in IC1 (less than 5sd; similarly with 𝑎  : for Long-chain 3hydroxyacyl-coenzyme A dehydrogenase deficiency (LCHAD), 𝑎   for Isovaleric acidemia (IVA), and 𝑎   for Glutaric aciduria I (GAI) patients), since loadings perform as “general weights” of the metabolites. Hence, it should be noted that by looking only at loadings, as is typically done, for example, for PCA, one may not be able to reveal the sources for the separation of a specific group of patients. In such a case, it is the table of AMSS values, which are computed for each group separately, that helps with a better understanding of the general picture. Int. J. Neonatal Screen. 2023, 9, x FOR PEER REVIEW 7 of 14 Figure 1. ICA score plots for components IC1 to IC8 of the discovery study: (A) Components IC1 and IC2. (B) Components IC3 and IC4. (C) Components IC5 and IC6. (D) Components IC7 and IC8. The dashed lines represent 3 (red),4 (green), and 5sd (blue) rule thresholds. 3.1.2. IC2—PKU × MSUD In IC2, there is a clear separation of PKU patients and controls. The routine diagnostic biomarkers based on the logic of the enzyme defects (Phenylalanine hydroxylase (PAH) in classical PKU and other tetrahydrobiopterin-turnover-related hyperphenylalaninemias (Table S1)) are Phe, and Phe to Tyr ratio. The concentration levels of Phe are increased since Phe is not converted to Tyr by PAH. Therefore, Tyr levels may be disproportional to increased Phe levels in PKU patients (statistically significant decrease in Tyr observed in PKU patients in our data; Welch Two Sample t-test, p = 1.098 × 10−7), which is reflected in Phe/Tyr ratio. In Table S2, the biggest difference in AMSS values can be observed between 𝑎   equal to −9.149 and 𝑎   equal to 0.043. Thus, PKU patients are shifted to negative values in the ICA score plot (Figure 1A). Values of 𝑎   and 𝑎   are almost similar for the two groups (−0.201 for patients and 0.001 for healthy controls, respectively) pointing to Tyr generally working in NBS as a reference metabolite in the Phe/Tyr ratio. One healthy control sample from NBS is located in the area exceeding the 5sd threshold near the PKU group. This patient has a Phe concentration of 119.99 µmol L−1 and Phe/Tyr ratio of 2.04. Since the cutoff values in NBS routine procedure for Phe and Phe/Tyr ratio are, respectively, 120 µmol L−1 and 2, the patient did not exceed both parameters and was therefore flagged as “negative” in the screening despite his values being borderline. A secondary effect in IC2 is a minor separation of Maple syrup urine disease (MSUD) patients, who are further better clustered in IC4 and described under that section. The separation here is due to increased loadings of clr coefficients of metabolites leucine/isoleucine (Xle) and valine (Val) (biomarkers of MSUD). The values of AMSS in MSUD patients compared to all other groups (Table S2), namely 𝑎   and 𝑎   , corroborate the similar nature of this effect as was the case with 𝑎   and 𝑎  : . Figure 1. ICA score plots for components IC1 to IC8 of the discovery study: ( A ) Components IC1 and IC2. ( B ) Components IC3 and IC4. ( C ) Components IC5 and IC6. ( D ) Components IC7 and IC8. The dashed lines represent 3 (red),4 (green), and 5sd (blue) rule thresholds. Int. J. Neonatal Screen. 2023,9, 60 7 of 14 In the loadings table (included in Table S2), some metabolites with no obvious pathobiochemical connection to MCAD, e.g., tetradecanoylcarnitine (C14) and tetradecenoylcarnitine (C14:1), are also enhanced with a non-zero loading value. It can be observed that the increased loadings correspond to different AMSS values for patients suffering from Very long-chain acyl-CoA dehydrogenase deficiency (VLCAD), i.e., slightly higher C14aVLCAD IC1 and C14:1aVLCAD IC1 . The result of this is a minor separation of the VLCAD group in IC1 (less than 5sd; similarly with C14:1aLCHAD IC1 for Long-chain 3-hydroxyacyl-coenzyme A dehydrogenase deficiency (LCHAD), C5aIVA IC1 for Isovaleric acidemia (IVA), and C5DCaGAI IC1 for Glutaric aciduria I (GAI) patients), since loadings perform as “general weights” of the metabolites. Hence, it should be noted that by looking only at loadings, as is typically done, for example, for PCA, one may not be able to reveal the sources for the separation of a specific group of patients. In such a case, it is the table of AMSS values, which are computed for each group separately, that helps with a better understanding of the general picture. 3.1.2. IC2—PKU ×MSUD In IC2, there is a clear separation of PKU patients and controls. The routine diagnostic biomarkers based on the logic of the enzyme defects (Phenylalanine hydroxylase (PAH) in classical PKU and other tetrahydrobiopterin-turnover-related hyperphenylalaninemias (Table S1)) are Phe, and Phe to Tyr ratio. The concentration levels of Phe are increased since Phe is not converted to Tyr by PAH. Therefore, Tyr levels may be disproportional to increased Phe levels in PKU patients (statistically significant decrease in Tyr observed in PKU patients in our data; Welch Two Sample t-test, p= 1.098 × 10 −7 ), which is reflected in Phe/Tyr ratio. In Table S2, the biggest difference in AMSS values can be observed between PheaPKU IC2 equal to − 9.149 and Pheacontrol IC2 equal to 0.043. Thus, PKU patients are shifted to negative values in the ICA score plot (Figure 1A). Values of TyraPKU IC2 and Tyracontrol IC2 are almost similar for the two groups ( − 0.201 for patients and 0.001 for healthy controls, respectively) pointing to Tyr generally working in NBS as a reference metabolite in the Phe/Tyr ratio. One healthy control sample from NBS is located in the area exceeding the 5sd threshold near the PKU group. This patient has a Phe concentration of 119.99 µ mol L −1 and Phe/Tyr ratio of 2.04. Since the cutoff values in NBS routine procedure for Phe and Phe/Tyr ratio are, respectively, 120 µ mol L −1 and 2, the patient did not exceed both parameters and was therefore flagged as “negative” in the screening despite his values being borderline. A secondary effect in IC2 is a minor separation of Maple syrup urine disease (MSUD) patients, who are further better clustered in IC4 and described under that section. The separation here is due to increased loadings of clr coefficients of metabolites leucine/isoleucine (Xle) and valine (Val) (biomarkers of MSUD). The values of AMSS in MSUD patients compared to all other groups (Table S2), namely XleaMSUD IC2 and ValaMSUD IC2 , corroborate the similar nature of this effect as was the case with C14aVLCAD IC1and C14:1aVLCAD IC1. 3.1.3. IC3—LCHAD ×GAI Elevated levels of C16-OH, C18-OH, and C18:1-OH are hallmarks of LCHAD deficiency patients in screening programs. These biomarkers are reflected in the reduced values of C16−0HaLCHAD IC3 , C18−0HaLCHAD IC3 , and C18:1−0HaLCHAD IC3 compared to all other groups (Table S2), which shift the patient’s group into negative values in the ICA score plot (Figure 1B). GAI patients are separated in the opposite direction compared to LCHAD patients in IC3. Due to the increased concentration levels of GAI biomarker glutarylcarnitine (C5DC) and its positive C5DCaGAI IC3 value (Table S2), the GAI patient is well separated from all other patients and the control group (Figure 1B) in the positive scores. 3.1.4. IC4—MSUD ×IVA In IC4, the separation of MSUD and IVA patients can be observed. MSUD patients are characterized by increased Xle and Val values. The positive AMSS values of the biomarkers Int. J. Neonatal Screen. 2023,9, 60 8 of 14 in these patients, XleaMSUD IC4 and ValaMSUD IC4 (Table S2), compared to the other groups result in their separation from the mass by shifting them further to the positive values in ICA score plot (Figure 1B). Among the highest AMSS values, there is also an increased C5aMSUD IC4 which corresponds to a reduction in C5 concentrations in MSUD patients (statistically significant decrease; Welch Two Sample t-test, p= 2.311 × 10 −5 ). With respect to C5 having the highest negative loading value in IC4 (see Table S2, note the opposite sign of C5 and Xle and Val loadings), it significantly contributes to the resulting high scores of MSUD patients compared to the rest of the data. Moreover, the C5/Xle ratio has been previously described [ 23 ] as a ratio detected by the CLIR Productivity Tools [ 24 ] in the data from NBS for MSUD in the Netherlands. The authors concluded that “the C5/Xle ratio is predominantly determined by the Xle concentration, and is not of added value”. This statement is slightly in disagreement with our results, where the effect of a reduction in C5 concentration itself can be observed. Thus, the C5/Xle ratio detected by CLIR may be a suitable parameter for the screening. IVA patients show elevated C5 values. Strongly negative loading of C5 shifts the respective C5aIVA IC4 value as well as the IC4 scores of the IVA patient to the negative values. IVA is described in more detail under section IC8 where the patient is separated without the loadings being highly influenced by other diseases. 3.1.5. IC5—Weights The clustering visible in IC5 (Figure 1C) was described in a previous publication [ 17 ], where loadings of C16, Val, C18:1, C18OH, and C0 discriminated patients with low birth weight. This finding is not at all straightforward and further research is necessary. 3.1.6. IC6—GAI The highly negative value of C5DCaGAI IC6 for the analyte C5DC (primary routine NBS biomarker increased in GAI patients) causes a clear separation of the GAI patient in the ICA score plot from all other groups (Figure 1C, Table S2). The second screening measure routinely used in NBS to diagnose GAI patients is C5DC/C16 ratio (Table 1). The C16aj IC6 (where jranges over all studied groups) is very similar among all patients and controls, showing C16 together with low sd as one of the viable “reference metabolites” in this component (Table S2). In the opposite direction to the GAI patient, one PKU patient exceeds 5sd due to low C5DC concentration. 3.1.7. IC7—ASA In IC7, a patient with Argininosuccinic aciduria (ASA; detectable in NBS based on increased levels of Citrulline (Cit) and Argininosuccinate (ArgSucc)) is separated. AMSS values of the two associated analytes, CitaASA IC7 and ArgSuccaASA IC7 , are higher for this patient compared to all other groups, assigning the ASA patient a positive score in ICA score plot (Figure 1D, Table S2). 3.1.8. IC8—IVA In NBS, IVA patients are detected using increased levels of C5 concentration and C5/C8 and C5/C2 ratios. The only distinctly separated patient on the y-axis of Figure 1D is an IVA case. The negative C5aIVA IC8 value compared to other groups and similar values for “reference metabolites” C8 (except for the group of MCAD patients, where it is a substrate of defective enzymatic reaction; see Section 3.1.1) and C2 across the groups, C8aj IC8 and C2aj IC8 , support the use of these analytes in denominators of the biomarker ratios in NBS (Table S2).