Reproducible analysis of disease space via principal components using the novel R package syndRomics
Abstract
Biomedical data are usually analyzed at the univariate level, focused on a single primary outcome measure to provide insight into systems biology, complex disease states, and precision medicine opportunities. More broadly, these complex biological and disease states can be detected as common factors emerging from the relationships among measured variables using multivariate approaches. ‘Syndromics’ refers to an analytical framework for measuring disease states using principal component analysis and related multivariate statistics as primary tools for extracting underlying disease patterns. A key part of the syndromic workflow is the interpretation, the visualization, and the study of robustness of the main components that characterize the disease space. We present a new software package, syndRomics, an opensource R package with utility for component visualization, interpretation, and stability for syndromic analysis. We document the implementation of syndRomics and illustrate the use of the package in case studies of neurological trauma data.
Introduction
The goal of the burgeoning field of precision medicine is to understand complex disease states and provide opportunities for deep patient phenotyping and highly targeted therapeutics. Precision medicine requires an understanding of multidimensional disease states. Yet, the analysis of biomedical data remains largely univariate, with response variables considered individually and reports involving several distinct analyses. This analytical approach limits our interpretation of the complexity of a disease by not considering the shared information across variables and potentially contributing to irreproducibility due to statistical limitations of multiple comparison testing. Understanding the full set of interrelated disease features through multivariate statistics is the goal of the growing domain of 'syndromics' (Ferguson et al., 2011). In particular, principal component analysis (PCA) and related multivariate statistics such as nonlinear PCA or factor analysis have been proposed as tools for extracting underlying factors or patterns (principal components [PCs]) reflecting disease states (Ferguson et al., 2013; Haefeli et al., 2017a; Haefeli et al., 2017b; Kutcher et al., 2013; Nielson et al., 2014; Nielson et al., 2015; Panaretos et al., 2017; Rosenzweig et al., 2010; Rosenzweig et al., 2018; Rosenzweig et al., 2019; Zhang and Castelló, 2017). There are several other multivariate methods that could be used for multivariate pattern detection: other ordination and dimension reduction techniques, cluster analysis, discrimination analysis, or the plethora of more recent machine learning methods. The use of any of these methods has its advantages and pitfalls (Everitt and Hothorn, 2011). We focus on PCA as being one of the most widely used method for pattern detection. PCA is a multivariate statistical procedure that allows for the generation of new uncorrelated variables, called PCs, as a weighted combination of the original variables (Abdi and Williams, 2010; Hotelling, 1933; Jolliffe and Cadima, 2016). These components are ordered such that the first component explains the major source of variance in the data, the second component the second largest source of variance, etc. The extracted components reflect the interrelation between all the original variables or features, allowing for disease pattern detection, guiding in the interpretation of disease complex space and overcoming univariate analysis limitations.
Despite the extensive use of PCA in some subfields of biological research and the increasing use of PCA for disease pattern discovery, there is very limited information in the literature that can guide applied biomedical researchers about its implementation and interpretation. Here, we offer a practical guide to the application of PCA for the extraction of disease patterns that conform the disease space, with focus on reproducibility. By no means can we cover the extensive field of PCA in the present document. Rather, we aim to provide an introductory manual to extraction of reproducible disease patterns using multidimensional analytics, directed to biomedical researcher practitioners while pointing to additional relevant sources of information. We introduce a software package for the R programming language called syndRomics, implementing some of the tools described here. We will illustrate the analysis workflow and the use of the package in experimental data from case studies in neurotrauma.
The key steps in disease pattern detection by PCA are shown in Figure 1. The syndRomics package offers functionalities that aid in these steps, building on the extensive PCA framework developed by the R opensource community. The package implements a novel visualization tool, the syndromic plot first published by Ferguson et al., 2013, as well as functions to quickly generate two other publicationready visualizations (a heatmap and a barmap). In addition, the package implements resampling strategies, providing datadriven approaches to analytical decisionmaking aimed to reduce researcher subjectivity and increase reproducibility. In particular, the package offers a function to extract metrics for component and variable significance by using nonparametric permutation methods (Landgrebe et al., 2002; Linting et al., 2011; PeresNeto et al., 2003), to inform component selection and component interpretation. Finally, the package incorporates functions to study component stability toward understanding the generalizability and robustness of the analysis (Cattell and Baggaley, 1960; Cattell et al., 1969; LorenzoSeva and ten Berge, 2006).
Results
We will describe the general steps to use PCA for syndromic analysis and illustrate the use of the syndRomic package along the analytical steps with two case studies of neurotrauma data. Details of the usage and implementation of the package and functions are described in the Materials and methods section. The full code reproducing the analysis can be found in the supplementary material. Code boxes in the text provide snippets illustrating the main sections of the code. The first case study is used as a tutorial to illustrate the steps of analysis; the second case study is discussed at the end of the results section. For the first case, we used a publicly available preclinical dataset on the Open Data Commons for Spinal Cord Injury (odcsci.org) (Callahan et al., 2017; Fouad et al., 2019). We selected a subset of the dataset with accession number ODCSCI: 26 (Ferguson et al., 2018) that has been previously used for deriving the socalled spinal cord injury (SCI) syndromics (Ferguson et al., 2013). The dataset contains 159 subjects (rats) that have been studied on different motor functional outcomes across time after cervical spinal cord injury. The subset chosen for the present analysis consists of 18 outcome variables measured at 6 weeks after injury. The included variables for this analysis are shown in Table 1. For additional details of these variables, see Ferguson et al., 2013.
Step 1: Extracting PCA solution from the data
There is extensive literature on performing PCA (Abdi and Williams, 2010; Jolliffe and Cadima, 2016; Zhang and Castelló, 2017). As a consideration, biomedical data aiming to capture the multivariate disease space usually contains variables of different types (i.e. categorical, continuous, etc.) and scales, known as ‘mixedtype’ data. Moreover, missing data is a common problem in biomedicine (Hollestein and Carpenter, 2017; Kaushal, 2014; Nielson et al., 2020) that needs to be solved to be able to apply most standard PCA algorithms. Therefore, some preprocessing transformations are usually applied before performing PCA. For example, linear PCA is sensitive to the scale of variables, thus when applying a linear PCA to continuous variables of different units or scales, a common practice is to scale the data to unit variance first (i.e. equivalent to performing the PCA on the correlation matrix). The use of the package to conduct syndromics analysis from linear PCA is illustrated on the first case study. In cases of datasets with mixed data types and/or nonlinear relationships between variables, nonlinear PCA with optimal scaling transformation (Linting et al., 2007a; Mair and Leeuw, 2019) has been previously used for disease pattern analysis (Rosenzweig et al., 2018; Rosenzweig et al., 2019). We used the syndRomics package to analyze patterns from a nonlinear PCA in the second case study. In cases with missing data, strategies such as data imputation or the use of PCA algorithms allowing missing values might be needed (Dray and Josse, 2015). While missing values analysis and dealing with missingness is an extensive topic that is not covered in detail here (Rubin, 1976; Buuren, 2018), the chosen case studies do contain missing values and illustrate how the package can help to determine the stability of the PCs when imputing missing values (see component stability section).
Another consideration is selecting which variables to include in the analysis. For PCA of experimental data where there are stratifying factors (e.g. control vs. treatment), it is important to leave out variables that directly capture the variance of these factors, which would bias PCA results toward separating the experimental groups. This bias is problematic since in syndromic analysis, the goal is to find the relationship between variables describing different diseases states in an unsupervised (i.e. not guided by our design) manner. For instance, if treatment indicators are included and the variance between treatment groups is high, the PCA solution would directly capture the experimental design and confound the multivariate patterns.
The disease components can be used in subsequent analysis as multivariate outcomes or predictor indicators (Haefeli et al., 2017a; Nielson et al., 2015; Rosenzweig et al., 2018; Rosenzweig et al., 2019). PCA is used to extract the correlation structure between variables, generating new independent variables as linear combinations. Beyond the use of PCs as proxies for disease patterns, the PCs can help mitigate issues that might appear when analyzing several variables such as multicollinearity, overfitting, and multiple testing (Altman and Krzywinski, 2018; Johnson et al., 1973; Lever et al., 2017).
The reader is referred to some materials of interest on considerations and limitations when conducting PCA and related methods for biomedical research (Jiang and Eskridge, 2000; Konishi, 2015; Nguyen and Holmes, 2019; Zhang and Castelló, 2017).
Case study: In the first case study, the goal is to run a linear PCA to study the motor function components 6 weeks after cervical spinal cord injury. This will summarize all motor function variables as a small set of independent components explaining different aspects of the motor behavior after an SCI. The data contains missing values (Figure 4—figure supplement 1), and therefore we performed missing values analysis before continuing with the workflow. Typically, the first step in missing values analysis is to determine patterns of missingness and classify missing values as missing completely at random (MCAR), missing at random (MAR) or missing not at random (MNAR) (Rubin, 1976). The type of missingness will guide the decision on which is an acceptable procedure to deal with missing values. For instance, deleting all subjects that contain at least one missing observation is common practice (aka listwise deletion or completecase analysis), but it is only acceptable if missing values are MCAR. Otherwise, the robustness and proper estimation of the missing values can not be guaranteed (Schafer and Graham, 2002; Buuren, 2018). In the example data, subjects have been pooled together from different experiments. We know that the observed pattern of missingness (Figure 4—figure supplement 1) is due to a set of animals where some of the outcome measures were not studied, suggesting that missing values are MNAR. We confirmed that missing values are not MCAR using a previously described test of MCAR (Jamshidian and Jalal, 2010) implemented in the MissMech package in R (Jamshidian et al., 2014), which rejected the hypothesis of MCAR missingness in our data. Thus, excluding subjects from the analysis is not justified. Instead, we have used multiple imputation through the mice R package (Buuren and GroothuisOudshoorn, 2011) to generate 50 imputed datasets and pooled them using the mean of each observation. We will illustrate on the component stability section how the syndRomics package can be used to determine the robustness of multiple imputation for disease pattern analysis. We extracted the PCA solution of the pooled imputed data using the prcomp() function in R after centering and scaling the data to unit variance (R code box 1). Other similar functions in R or other software can be used.
R Code Box 1
pca<prcomp (pca_data, center = TRUE, scale. = TRUE).
Step 2: Component selection: how many components to keep?
The first question, after running PCA for extracting the disease components is usually to determine how many PCs are relevant. As a general consideration, the PCs with lower eigenvalues (i.e. explain less variance) have a higher chance of representing noise in the data (Jolliffe and Cadima, 2016), questioning their generality and value. The goal is to determine the minimal set of components that can be used to describe the disease space. Importantly, there is not a single, specific rule for this determination. A common method in PCA and related methods is the Scree test by Cattell, 1966, where all PCs are ordered in descending rank by their eigenvalues, and PCs above the ‘elbow’ are retained. Another criterion is the eigenvalue greater than one rule which is applied to standardized PCAs (from the correlation matrix) with the criteria of only keeping PCs with an eigenvalue (i.e. the variance of a component) above 1 (Guttman, 1954; Kaiser, 1960). A more thorough description of these and others methods can be found elsewhere (Glorfeld, 1995; Horn, 1965; Vitale et al., 2017; Zwick and Velicer, 1986). Simulations have shown these methods (specially the eigenvalue greater than one rule) to be less robust than a resampling approach for selecting the number of relevant components (Zwick and Velicer, 1986). The syndRomics package incorporates a nonparametric permutation test approximated through Monte Carlo resampling of the total ‘variance accounted for’ (VAF) of each PC to aid in the selection of relevant PCs (Buja and Eyuboglu, 1992; Glorfeld, 1995; Horn, 1965; Landgrebe et al., 2002). The permutation test can also assist in component interpretation by studying the contribution of each variable to the PCA solution (Buja and Eyuboglu, 1992; Linting et al., 2011) as we will see in the next section.
The goal of the permutation test is to determine whether the extracted PCs can be considered to be generated notatrandom. This method has been shown to outperform parametric tests for PCA in situations similar to biomedical data where sample sizes are relatively small and the data rarely comply with the assumptions of the models (Buja and Eyuboglu, 1992; Horn, 1965; Zwick and Velicer, 1986). In that regard, a hypothesis test is defined as:
H_{(null)}:PC VAF is indistinguishable from a random generation
H_{(alternative)}:PC VAF is different from random
The p values are calculated by:
where $q$ is the number of times the chosen metric is higher in the permuted distribution than in the original PCA solution and $P$ is the number of permutations (Buja and Eyuboglu, 1992). Rejecting the null hypothesis is interpreted as evidence of the tested PC being generated from true signal and not by random noise. This sets a lower bound for which PCs to consider 'important' above noise, but does not indicate the magnitude of the 'importance', which is represented by VAF. Importantly, for datasets with several directions of variance and high signaltonoise ratio, PCs with low VAF can still be statistically significant. The value of interpreting such PCs must be judged by the researcher in the context analysis in question. It is also important to consider how big $P$ needs to be when performing resampling, such as with the permutation test incorporated in the package. The reader should note that the lowest $p$ value that can be calculated is dependent on $P$. For example, if $P$ is set to a value of 10 (a relatively low value), the smallest p value that can be detected is 0.09, which occurs when $q=0$. Accordingly, $P$ should be set high enough to reach the desired minimum p value. Moreover, simulation studies have shown that $P$ under 99 have low power and a minimum of 499 permutations is recommended (Buja and Eyuboglu, 1992; Abdi and Williams, 2010; Linting, 2007). By default, we have set the number of permutations to 1000 (smallest p value approximately equal to 0.001) as this has been shown to produce good results (Landgrebe et al., 2002; Linting et al., 2011). Users of the package should keep in mind that higher numbers of permutations will increase computation time with potentially only a small gain on the approximation. Our simulations indicate that between 500 and 1000 permutations provide a good compromise between computing time and precision in estimating confidence intervals, depending on the data volume (Figure 4—figure supplement 1). The package implements a single permutation strategy for VAF, the socalled permD (permutation of the entire data set) (Buja and Eyuboglu, 1992; Linting et al., 2011) where variables are permuted independently and concomitantly (Figure 2A) opposed to permV (permutation of a single variable) (Linting et al., 2011) where variables are permuted one at the time (Figure 2B). These methods are further discussed on the component interpretation section.
R Code Box 2.
permut_pc_test (pca, pca_data, p=10000, ndim = 5, statistic = 'VAF', perm.method = 'permD').
Case study: After performing a PCA, we first determined the number of components that can be regarded as informative. Several criteria can be used as mentioned earlier. Here, we opted for the permutation test of VAF, computed using the permut_pc_test() function (R Code Box 2). We have applied this test to the data using 10,000 permutations. The results show that the three first PCs (PC1, PC2, and PC3) are significantly different from random at an alpha of 0.05 adjusting the p value (Figure 3), and therefore we will keep these three PCs for subsequent analysis. PC1 accounts for 32.9% of the variance, PC2 18.3% and PC3 9.8%.

Figure 3—source data 1
 https://cdn.elifesciences.org/articles/61812/elife61812fig3data1v2.csv

Figure 3—source data 2
 https://cdn.elifesciences.org/articles/61812/elife61812fig3data2v2.csv

Figure 3—source data 3
 https://cdn.elifesciences.org/articles/61812/elife61812fig3data3v2.csv
Step 3: Component interpretation: what do these components mean?
A key part of the analytical workflow is the interpretation of the main components, where the most relevant PCs can be used to represent the correlation between the original variables as a proxy for multivariate disease patterns. Each component is composed of a weighted combination of all the variables. Some components might be explained by only a few variables with high importance, whereas others might have several variables with important contributions to them. There are a few metrics that can be used for interpreting the relation between the original variables and the PCs (Abdi and Williams, 2010). In the syndRomics package, we use the standardized loadings or correlation vector coefficients (Jackson and Hearne, 1973), and the communalities, which are the sum of squared loadings for each variable across selected PCs representing how much of the variance of each variable can be explained by the total number of kept components. Loadings can be interpreted as the Pearson’s r correlation coefficient between a PC and a variable, and it is used to assess the contribution of individual variables on each PC and the direction on which the variable moves along the PC (i.e. opposite or same direction as in the interpretation of a correlation). Communalities can be interpreted as the global impact of a variable in the chosen PCA solution.
In general, the strategy consists of determining a threshold for the absolute value of loadings or the communalities above which variables are considered to have important contribution in the definition of a component or across the chosen PCs. For example, if a threshold of loading > 0.2 is chosen, all variables for a given PC with a loading > 0.2 or a loading < −0.2 will be considered to contribute on the PC (aka salient variable). The matter then turns to determining an appropriate threshold. Some somewhat arbitrary rules of thumb for the loadings have been established. However, those have a strong determination in psychological studies and whether they are appropriate in biomedical research has yet to be verified. An alternative ‘quasiinferential’ method is to use permutation test as discussed above for PC VAF but testing for metrics of variable contribution such as loadings (Buja and Eyuboglu, 1992; PeresNeto et al., 2003) or communalities (Linting et al., 2011). Using resampling strategies, these permutation methods offer datadriven determination of variable importance and contribution, which might reduce subjective biases. Thus, rejecting the null hypothesis for a given metric, variable and PC, suggest that such variable has a contribution onto the construction of the component that is above what is expected by random noise. As in the case of VAF, this establishes a lower bound for loadings or communalities below which they should be considered noise. In situations with stable solutions and high signaltonoise ratio, low loadings or communalities might still be statistically significant, but the contribution of the variable should be gauged respect to other variables. In the package, we have incorporated permutation test of the loadings as in Buja and Eyuboglu, 1992; PeresNeto et al., 2003 that can serve to determine the loading threshold, where the variables are permuted independently and concomitantly (Figure 2A and C). Linting et al., designed and tested an strategy for the communalities where only one variable is permuted at the time, showing great results in determining the contribution of variables using communalities (Linting et al., 2011). This method has resulted in better determination of the significant contribution of variables on the PCA solution with higher statistical power and proper type I error, and therefore has been incorporated in the package as the default method for both the communalities and the loadings (Figure 2B and D). Following Linting et al., terminology, users can specify the permutation strategy for the loadings as one variable at the time (permV, as in [Linting et al., 2011]) or as all the variable together (permD, as in [Buja and Eyuboglu, 1992; Linting et al., 2011; PeresNeto et al., 2003]). See Materials and methods for details on the permutation algorithms. In addition to permutation strategies, the package implements bootstrapping methods for constructing confidence intervals of component loadings and communalities that can also facilitate PCs interpretation (see component stability).
The selection of number of permutations in this case follows similar rationale as described above for the VAF. It is important to note that the minimal number of permutation needed to have enough statistical power and precision will depend on the size of the dataset, both on the number of variables and samples (Buja and Eyuboglu, 1992; Figure 4—figure supplement 1). There is also the understanding that while the permD strategy is less robust than permV as suggested by Linting et al., the computational time increases considerably since variables are permuted one at the time. Moreover, adjusting p values for multiple testing might be recommended depending on the sample size. Linting et al., suggested controlling for false discovery rate (FDR) using the Benjamini and Hochberg (BH) (Benjamini and Hochberg, 1995) method. As a rule of thumb, these researchers advised to only use multiple testing correction (for FDR) when the data contains at least 20 variables and 100 observations or subjects, and to use the uncorrected pvalues otherwise (Linting et al., 2011). pValue adjustment has been incorporated in the permutation function on the package, with controlling for FDR by BH as default.
The reader should be cautioned against overinterpreting or misinterpreting the meaning of a PC. The interpretation can be subjective, and unconscious biases can be reflected on the interpretation of PCs. The tools offered by the package help mitigate potential subjective biases, although data biases will affect the results. Another consideration is that it is possible that some of these metrics seem to ‘contradict’ each other. For example, there is the possibility that a component has an important contribution to the variance of the data (high VAF) and yet all the loadings be small. Contrary, a component with a small set of high loadings could be considered to be insignificant by permuting its VAF (Buja and Eyuboglu, 1992). As in any analytical approach, domain knowledge is critical for the interpretation of disease components.
Case study: After deciding to keep three components, we studied the communalities and loadings to determine their identity. Here, we applied the permut_pc_test() function (R Code Box 3) setting the argument statistic = ‘commun’ or ‘s.loadings’ and the perm.method = ‘permV’ and using the BH method for controlling for FDR. The results of the permutation test on the communalities can be seen in Figure 3B and in Table 3. We can appreciate that all variables are significantly represented by the three chosen PCs, although there are five variables with communality less than 0.5, indicating that the retained PCs only explain 50% of the variance on these variables. In PCA, communalities can suggest which variables do or do not contribute to the extracted components altogether. Considering the loadings, the results for PC1, PC2, and PC3 are shown in Figure 3C and in Tables 4, 5 and 6, respectively. One can appreciate that the cutoff loading for significance at alpha 0.05 using the adjusted p value is approximately 0.21 for PC1, 0.25 for PC2 and 0.4 for PC3. This behavior of different thresholds for significance has been previously described (Buja and Eyuboglu, 1992) and reflects the fact that PCs accounting for less variance might contain more random noise, thus needing a higher loading for a variable to be considered as an important contributor. Loadings are indicative of both strength of association between a variable and a PC and the direction in which they interact. For reading on the interpretation of the loadings and components in this case study, see Ferguson et al., 2013.
R Code Box 3
permut_pc_test (pca, pca_data, p=1000, ndim = 3, statistic = 's.loadings', perm.method = 'permV').
permut_pc_test (pca, pca_data, p=1000, ndim = 3, statistic = 'communa', perm.method = 'permV').
Step 4: Component stability: how robust are the components?
The presence of a syndrome or disease pattern, represented by a component, should hold true regardless of variations in experiments or metrics meant to measure that same pattern. For example, two experiments with different subjects but the same collected variables should result in inferentially equivalent components if they are true features of the disease and not experimental artifacts. The sensitivity of PCs to experimental, metric, or other forms of variation is termed ‘component stability’. Components from different PCAs (from different experiments as an example) that are extremely similar are considered to be a stable, and characterizing component stability is important to determine the robustness of the initial PCA (Guadagnoli and Velicer, 1988; Linting, 2007). A robust PC would be largely unaffected by data variations (i.e. low sensitivity). The goal of the stability analysis is to determine such sensitivity.
Given that performing multiple replication experiments in biomedical research is not always possible, component stability can be approximated by resampling techniques such as bootstrapping (Babamoradi et al., 2013; Linting et al., 2007a; Timmerman et al., 2007; Zientek and Thompson, 2007). Bootstrap methods for component stability have been extensively studied, but users should be aware of the limitations and advantages of these methods and their performance for component stability depending on the use case (Babamoradi et al., 2013; Guadagnoli and Velicer, 1988; Linting et al., 2007b; Timmerman et al., 2007; Zientek and Thompson, 2007).
The package implements functionalities to help study the component stability affected by data selection variability by implementing bootstrapping methods (Babamoradi et al., 2013; Linting et al., 2007b; Timmerman et al., 2007; Zientek and Thompson, 2007) and stability metrics. The default method used in the package is the simple or ordinary bootstrap consisting of generating a new sample that has the same size (i.e. same number of subjects or observations) and same variables as the original data, but where the subjects have been randomly selected from the data with replacement (Figure 4A). This process is repeated several times (here referred as b times) to generate a sample of bootstrapped data. In the first case example, each of the b bootstrapped samples contain 159 subjects and 18 variables, but one subject might appear more than once and another subject might not show up in a specific sample.
Component stability can be studied at the whole component level or at the level of the individual variables through the loadings and communalities. The package implements component similarity indexes (aka factor matching indexes)(Cattell and Baggaley, 1960; Cattell et al., 1969; Guadagnoli and Velicer, 1991) as metrics to study the stability of PCs. These metrics can be used to determine the similarity between the different bootstrapped samples, to test the validity of the extracted component under two or more experimental conditions, to assess the multidimensional equivalence of two or more replication experiments, or to determine the impact of imputing missing values.
Case study: To understand the sensitivity of variables and components to experimental variations, we used the pc_stability() function with b = 1000 bootstrapped samples (R Code Box 4), setting the sim argument to ‘balanced’ to perform balanced bootstrapping. The function will return the average of the loadings and the specified similarity metric across all the b samples as well as the specified confidence interval. For this example, the 95% CI (accelerated and biascorrected, see Materials and methods) and the bootstrapped average can be seen in Figure 5. In general, the original loadings are close to the bootstrapped average which indicates that the results are unbiased. Moreover, the confidence regions for most higher value loadings are reasonably small, suggesting that these loadings are stable to experimental variation. In addition, the similarity metrics for the three PCs suggest component stability, meaning that the composition of the components is also stable. The accepted values for these metrics indicating stability might vary by field and the metric of interest. Some indicative values are mentioned in the respective method section, but the user should be aware of subjective biases when determining a threshold for considering stability (LorenzoSeva and ten Berge, 2006). Finally, the function will return the average and percentile CI of the communalities, which can be used to assess which variables are more stable in the selected PCA solution.
R Code Box 4.
pc_stability (pca, pca_data, B = 1000, ndim = 3, s_cut_off = 0.1, test_similarity = T, similarity_metric = 'all', sim = 'balanced', barmap_plot = T).
Assessing the impact of imputation methods on introducing noise when dealing with missing data should be considered. As described earlier, we used multiple imputation to generate the dataset for analysis. Multiple imputation generates m complete datasets where the imputed values might vary, but the observed values are the same (in the case study m = 50). We used the stability analysis described above to determine the sensitivity of the PCA solution to variations introduced by imputing missing values. We calculated the similarity metrics between all the 50 imputed datasets for the first three PCs, as well as the loadings. We observed high similarities between the PCs obtained from the imputed datasets (Table 7) and the loadings showed narrow CIs (Figure 3—figure supplement 1). We concluded that multiple imputation has produced stable solutions with acceptable impact on both the component and variables. A future version of the package might include more robust methods for pooling and testing multiple imputation in PCA context (van Ginkel and Kroonenberg, 2014). Altogether, the results suggest reliable and robust PCs extracted from the original data.
Step 5: Component visualization
Communicating the analysis is a necessary part of the workflow. Although we have included this at the end of the use case, visualization can be also used for aiding in component selection, component interpretation and component stability analysis. There are several ways a PCA solution can be visualized. Here, we describe the plots implemented in the syndRomics package.
We have coded three types of plots (syndromics plot, heatmap, and barmap) using the grammar of the graphics framework (Wilkinson, 2005) implemented in R by the ggplot2 package. This allow users to customize the plots using the rich landscape of the ggplot2 universe. The syndromic plot was first published by Ferguson et al., 2013 and represents PCs as the center of a Venn diagram (Figure 1A), consisting of (1) a middle convex triangle displaying the ‘variance accounted for’ (VAF) for a given PC and (2) radial arrows pointing to the center of the triangle for each variable with a standardized loading above a certain threshold (Figure 6A–C). The width of each arrow and the color saturation are proportional to the magnitude of the standardized loading they represent. The color of each arrow additionally differentiates between positive or negative loadings (e.g. blue represents a loading of +1, red represents a loading of −1, and white represents a loading of 0). Syndromic plots are especially useful for conveying PC identity in an easy to understand, concise way for publication. Heatmap and barmap plots are alternative visualizations of the loadings beyond the syndromic plot. The major difference between these two plots and the syndromic plot is that both the barmap (Figure 5A) and heatmap (Figure 6D) plots display all variables (or a manually selected subset) instead of only the ones with loadings above a given threshold. The absolute loadings that exceed a cutoff threshold can be noted (e.g. by a star *). Moreover, in the case of barmap plots, the cutoff is represented in the graph by vertical lines. This is particularly useful when there are too many abovethreshold variables, which would crowd the syndromic plot visualization, or when comparing loadings between PCs more easily. In addition, barmaps are useful for documenting the results of the resampling procedures since error bars can be used to represent the variation of the metrics over the resampling. The permut_pc_test() and pc_stability() functions return such plots.
Case study (R Code Box 5): We can visualize the three selected PCs using the plotting functions in syndRomics (Figure 6). In this case, we chose to represent PC1, PC2 and PC3 using the syndromics plots (Figure 6A, B and C, respectably) using a cutoff threshold of 0.45. Notice that this is higher than the threshold for significance found by the permutation analysis given the high number of variables. The full loading pattern of the three first PCs can be visualized by a heatmap (Figure 6D), where we have chosen a different cutoff for each PC (0.21, 0.25, and 0.4 for PC1, PC2, and PC3 respectively), or a Barmap (Figures 3B and 5A). Barmaps can be obtained for the loadings (barmap_loadings()) or for the communalities (barmap_commun()).
R Code Box 5
syndromic_plot (pca, pca_data, cutoff = 0.45).
heatmap_loading (pca, pca_data, ndim = 3, cutoff = c(0.21,0.25,0.4), star_values = T, text_values = F).
Case study 2
In the second case study, we used selected variables from the Transforming Research And Clinical Knowledge in Traumatic Brain Injury (TRACKTBI) pilot study (Yue et al., 2013) that were analyzed previously and made publicly available (Nielson et al., 2017). The released dataset version contains 586 deidentified human subjects who were enrolled in the TRACKTBI pilot study and the 26 selected variables previously analyzed (Nielson et al., 2017). These variables are a subset of brain imaging results, outcome metrics and genetic polymorphism (Table 2). The goal is to describe patterns of association between these three categories of variables. A noticeable difference between this dataset and the one used in the first case study is that here we are dealing with a mixed type dataset, where some variables are continuous, some nominal and some ordinal. Therefore, we performed a version of nonlinear PCA that allows for the extraction of patterns in these kinds of data. The syndRomics package has been programmed to work with the results of the princals() function from the Gifi R package. The code for this analysis is found in the supplementary material.
Missing data analysis showed an overall 21.2% missingness distributed between the outcomes and genetic polymorphism variables (Figure 7—figure supplement 1). With the exception of ‘MRI results’ that has high missingness (61.7% of the observations), all imaging variables are complete. ‘MRI results’ variable was excluded from the analysis. The subsequent test for MCAR suggest that there are 17 different patterns of missingness and that the hypothesis of MCAR can be rejected overall (pvalue<0.001). Thus, excluding subjects from the analysis is not justified (Schafer and Graham, 2002; Buuren, 2018). We instead performed 50 multiple imputations using the mice R package as in the first case study. The 50 imputed datasets where then aggregated to perform nonlinear PCA using princals() (see Materials and methods for details).
Permutation test of PC VAF suggests that the first 6 PCs contain information that can be regarded as significant above random chance. Although a deep analysis of these six PCs might be of interest, the first three PCs explain the major variance (25.6%, 10.6%, and 9.8%, respectively). Therefore, we focused on interpreting these for illustration purposes (Figure 7 and Tables 8–10). The first PC significantly loaded highly on two genetic variants in opposite directions (SNP_DRD2 loading = −0.677, SNP_ANKK1_Gly318AR loading = 0.661) as well as outcomes of neuropsychological function at 6 months after TBI (CVLT_long loading = −0.614, CVLT_short loading = −0.542)(Figure 7A–C). All other variables also significantly loaded on to PC1, but with loadings ~ 0.3 (Tables 8–10) suggesting that their contribution in PC1 identity is less important. Lower values in CVLT (California Verbal Learning Test) suggest learning and memory impairments, which are well known after TBI. Given the negative loadings for the included CVLT measures (short and long recall), negative values in PC1 might reflect better CVLT outcomes at 6 months after TBI. The stability of the PC1 pattern to multiple imputation is relatively low, with higher loadings showing high variation (Figure 7—figure supplement 1, Table 11), emphasizing the importance of studying stability of components to multiple imputation. Nonetheless, the bootstrapped loadings were stable (Figure 7B), and decay in CVLT performance after TBI has been previously associated to polymorphisms in DRD2 and ANKK1 genes (Failla et al., 2015; McAllister et al., 2008; Nielson et al., 2017; Yue et al., 2017), providing literature validation of PC1. The variables with higher positive loadings in PC2 were related to the imaging findings and negative loadings with global function outcomes at 3 and 6 months after TBI (GOSE score)(Figure 7. A, D). Lower scores in GOSE are indicative of lower global function and positive values in imaging findings are suggestive of a bigger or more noticeable brain damage. PC2 presented the higher stability to both resampling and to multiple imputation (Figure 7—figure supplement 1, Table 11). Altogether, PC2 might be interpreted as a surrogate for ‘TBI severity’, where higher positive values would indicate higher brain damage with less function at 3 and 6 months after injury, a signature described in the previous analysis of this data (Nielson et al., 2017). Finally, given the instability of PC3, with most loadings being considered nonsignificant by the permutation test and the high variance to multiple imputation (Figure 7—figure supplement 1, Table 11), PC3 can not be interpreted with certainty, and we should not attempt its explanation.
Discussion
Biomedical research needs more multivariate analytics to help realize the potential of precision medicine. While multiple variables are collected in typical preclinical experiments and clinical trials, univariate statistics continue to be the major analytical and decisionmaking approaches across the different biomedical fields, narrowing our understanding of the complexity of any disease. With the advent of ‘omics’, analytical approaches for highdimensional data have started to become more prevalent for the analysis of biological data. Yet, outside the realm of medical bioinformatics, biomedical research continues to be, for most part univariate. The lack of multivariate approaches in analyzing biomedical data can cause biases and constraints to the interpretation of the results and contribute to the lack of reproducibility and benchtobedside translation (Ferguson et al., 2011; Huie et al., 2018).
The extraction of disease space, through the use of multivariate methods, can increase our understanding of complex relationships commonly present in biomedical data while preventing some of the issues associated with an excessive use of univariate analytics such as multiple comparison testing and associated false discoveries by chance (Benjamini and Hochberg, 1995; Ferguson et al., 2013; Krzywinski and Altman, 2014). For example, it is common in biomedicine to measure several behavioral and histopathological outcomes that are analyzed independently at the univariate level. This approach increases the chance of falsepositive results due to the accumulation of type I testing errors (Benjamini and Hochberg, 1995; Dunn, 1961; Krzywinski and Altman, 2014). Although there are methods to correct for errors when running numerous tests such as multipletesting correction, their use in biomedicine outside of bioinformatic analysis is scarce. Even when correcting for multiple testing, performing several univariate analyses limits our understanding since univariate analysis does not allow us to study and infer the relationship between measures that might capture different aspects of the matter of study. In our first example case study, several functional tests can be used to study the recovery of forelimb motor function after cervical spinal cord injury in animal models. Each test further contains multiple measures about particular aspects of recovery. Knowing the relationship between these measures through multivariate approaches can increase our understanding of the matter of study while reducing the burden of multiple testing (Ferguson et al., 2013). Importantly, it is also possible that a single univariate test that does not produce significant results misses true biological effects, while a multivariate analysis including the same variables can find patterns and relationship between variables that are significant. Syndromic analysis is, therefore, a framework that uses multivariate analysis of biomedical data in a holistic way, aiming to reveal interactions within complex (patho)physiological niches, that would be otherwise challenging to discern. Applying syndromic analysis to biomedical data will help uncover the complex relationships of variables and features that constitute different disease and biological states and ultimately accelerate research toward precision medicine.
The syndRomics package implements several functionalities for the visualization, the interpretation, and the analysis of the stability of principal components to facilitate the extraction and analysis of disease patterns. We have demonstrated its usage, showing the potential of the package to support PCAbased analysis in understanding disease complexity. Although the core functionalities of the package are included, future versions might also implement outputs from other PCA functions as inputs, such as those from the PCA functions in the FactoMineR package (Lê et al., 2008) or the psych package (Revelle, 2017), allowing for better integration to the PCA landscape in R. In addition, other algorithms of bootstrapping and permutation methods for PCA solutions could be incorporated to increase the options and better adapt to the specifics needs of the user (Hong et al., 2006; Linting et al., 2011; Vitale et al., 2017; Zientek and Thompson, 2007).
Here, we emphasize guidance and tools for robust determination of PCAbased disease patterns. We have incorporated resampling methods aiming to reduce subjective biases and to study the stability and generality of the analysis. Although we have shown the use of these functions in different contexts along the process, much more work can be done to extend syndRomics. For example, we demonstrated the stability of our analysis under multiple imputation, and future research could investigate number of multiple imputations or missing conditions necessary for stable disease pattern detection. In addition, visualization features of syndRomics may be extended to help interpret disease patterns resolved by other multivariate or machine learning tools involving structure coefficients or feature impact scores. The syndRomics resampling methods could also be used to estimate the sample size required for stable PCs in the context of syndromic analysis, allowing for sample planning. The implementations in the package are thus positioned to empower both biological and statistical research toward understanding complex biology and diseases.
Materials and methods
Availability and requirements
Request a detailed protocolThe code to reproduce this analysis can be found in the supplementary material. The data for the first use case comes from the ODCSCI (Open Data Commons for Spinal Cord Injury, RRID:SCR_016673, http://odcsci.org), ODCSCI:26 dataset (https://scicrunch.org/odcsci/about/odcsci_26). The data for the second use case comes from TRACKSCI and can be downloaded from 10.1371/journal.pone.0169490. The package can be installed from GitHub (https://github.com/ucsffergusonlab/syndRomics) where installation instructions, package manual and examples of usage are provided. Descriptions of the arguments and function usage can be found in the internal package documentation once installed or in the package manual. The package has been implemented in R (R Development Core Team, 2019) through RStudio (Team RS, 2018) using a few other packages beyond the ones bundled in R as dependencies: dplyr (Wickham et al., 2018), ggplot2 (Wickham, 2016), stringr (Wickham, 2019), tidyr (Wickham and Henry, 2020), ggrepel (Slowikowski, 2019), ggnewscale (Elio Campitelli, 2020), pracma (Borchers, 2019), png (Urbanek, 2013), boot (Canty and Ripley, 2019; Davison and Hinkley, 1997), rlang (Henry and Wickham, 2020), and Gifi (Mair and Leeuw, 2019).
Package implementation
Request a detailed protocolThe syndRomics package offers two major functionalities for the purpose of aiding in the process of syndromics analysis: (1) visualization functions and (2) functions incorporating resampling methods to determine stability and inference of PCs.
Visualization functions
The visualization functions are: syndromic_plot(), heatmap_loadings(), barmap_loadings(), barmap_commun() and VAF_plot(). For the visualization functions, the user can pass an R data.frame object with the standardized loadings (or other metrics) obtained by running PCA and related multivariate methods in their preferred software. We opted for this approach to avoid requiring specific implementations of PCA. Loadings obtained from any PCA solution can be easily formatted for usage with the syndRomics visualization functions. All functions in the package that takes a data frame as argument use the same format (Table 12): variables are organized as rows, and the first column is called ‘Variables’ and contains the names of the respective variables. The other columns contain the PC loadings and are named ‘PC1’, ‘PC2’, etc. Alternatively, the visualizations can also be obtained from the output of the prcomp() function in the stats package in R (linear PCA) or from the output of the princals() function in the Gifi package in R (nonlinear PCA by categorical PCA). Finally, the results from pc_stability() and permut_pc_test() can be passed to the plot() generic function in R as the package incorporate the corresponding S3 method for ‘syndromics’ class object.
syndromic_plot ()
Request a detailed protocolThe list of arguments for the syndromic_plot() function are presented in the package manual. The syndromic_plot() function will internally call extract_syndromic_plot() function (see utility functions) and return a list of ggplot2 objects containing the syndromic plot for the first ndim PCs. For example, if ndim = 5, a syndromic plot for PCs 1 to 5 will be generated. Another important argument is the cut_off, which determines the threshold of absolute standardized loadings to consider for plotting. This argument is chosen by the user and is required (with no default). Another required argument is VAF in case the syndromic_plot() function is called using a data.frame input. If the output of the prcomp() or princals() functions is used, the syndromic_plot() function extracts VAF internally and the userdefined VAF will be ignored. When required, VAF is a character vector of the form ‘XX%”,”XX%”, etc., where XX is the VAF for each PC to plot, starting with the first PC, followed by the second, etc. (e.g. c(‘60.1%”,”25.3%”) for PC1 and PC2, respectively). An issue we found during the implementation is that the arrow visualization does not display correctly in the R graphical device on Windows machines. Rendering the plot into *.pdf format, for instance using the ggsave() function from the ggplot2 package, solves the problem.
heatmap_loadings(), barmap_loadings() and barmap_commun()
Request a detailed protocolMost of the functionalities described for the syndromic_plot() function also apply for the heatmap_loading(), the barmap_loading(), and the barmap_commun() functions. A noticeable difference in barmap_loading() is that the function will plot the PCs specified in ndim instead of the first ndim components. For example, if ndim = c(3,4,5), components 3, 4, and 5 will be plotted. This allows for more flexibility on which components to plot, such as isolating a single component (e.g. ndim = 3 will only plot component 3).
VAF_plot()
Request a detailed protocolThis function can be used to plot a VAF plot from a prcomp() or princals() output. There are two style options, ‘line’ or ‘reduced’.
Resampling functions
There are two major functions using resampling methods, the permut_pc_test() function that implements nonparametric permutation test for either PC VAF for aiding in component selection or PC loadings and communalities for aiding in component interpretation, and the pc_stability() function that implements bootstrapping of PC loadings for stability analysis. These functions take as input the output of the prcomp() or the princals() functions in R as well as the original dataset used on these functions as inputs. The specific call of prcomp() or the princals() used to obtain the original PCA solution is passed down to the resampling functions in the syndRomics package, ensuring that the same arguments are used for resampling (with the exception of the data argument on the original prcomp() or the princals() call, that will be internally changed for each resampling iteration).
permut_pc_test()
Request a detailed protocolIn the syndRomics package, the null distribution for the permutation test is generated by permuting the values of each variable independently and concomitantly several times (permD) or permuting one variable at the time (permV) and rerunning the PCA on each permuted sample (Figure 2; Buja and Eyuboglu, 1992; Glorfeld, 1995; Linting et al., 2011). When permV method is selected to measure the impact of permuting on loadings, a step of Procrustes rotation of each loading matrix toward the original loading matrix is added to resolve sign reflection, rotation indeterminacy and component translocation (Figure 2D, see pc_stability for detailed explanation). This step is not performed when the analysis is performed on the communalities since are invariant to such PCA resampling issues (Linting et al., 2011). Confidence intervals of the permuted distribution (null distribution) are calculated using the (1α)x100% (percentile) of the distribution (Buja and Eyuboglu, 1992).
The function calls the permut_pca_D() or permut_pca_V() utility generic function internally to generate the permuted distribution of the selected metric (either “VAF”,“s.loadings” or “comuna”) using either the prcomp() function for linear PCA or the princals() function for nonlinear PCA implemented as S3 R method for the class “prcomp” or “princals”. If “VAF” is specified the permD permutation will be used, ignoring the input of the user on the perm.method argument, returning a matrix containing the VAF for the original PCs, as well as the average and the CI of the permuted VAF distribution. In case “s. loadings” or “communa” are specified, the specified permutation method will be considered (i.e. permD or permV) and the function will return a list of matrices, one for each selected PC, with the original loadings, and the average and CI of the permuted loadings distribution. In both cases, p values are calculated as described in the main text and returned. Adjusted p values using the specified method in the adjust.method argument are also returned.
pc_stability()
Request a detailed protocolComponent stability can be studied at the whole component level, known as factor invariance, or at the level of the individual loadings. We have implemented both options in the package. By default, the pc_stability() function returns the average and the accelerated and biascorrected 95% confidence intervals (CI) of the loadings of the bootstrap distribution (Efron, 1987). Depending on the sample size and the number of chosen resamples, the biascorrected CI will fail and the percentile (1α)x100% CI will be returned (with corresponding notification). In addition, component similarity or factor matching metrics can be computed by setting the test_similarity=TRUE, which will call the component_similarity() function. For each of the specified similarity metrics, this function returns the average of the metric and its confidence interval (95% CI by default) by the percentiles of the bootstrap distribution. The confidence level and the CI method for the loadings can be changed by changing the conf and ci_type arguments. The function uses the boot() function for generating the bootstrapped samples and the boot.ci() function for extracting the confidence intervals of the loadings. Both boot() and boot.ci() are from the boot package in R. This allows the use of different bootstrapping strategies such as simple or ordinary bootstrapping (by default) or balanced bootstrapping. The reader is referred to the boot package documentation for more details on the different sim methods.
A major problem of performing resampling procedures in PCA is what is known as indeterminacies that can invalidate comparing between bootstrapped samples (Babamoradi et al., 2013; Chan et al., 1999; Linting, 2007; Timmerman et al., 2007; Zabala and Pascual, 2016). Sign reflection refers to the change of sign on the component loadings in a PC given slight variation of the data. In addition, slight data variation can also cause component/factor translocation, the change in the position of a component in the PCA solution (e.g. PC1 shifts to the position of PC2), especially when two components have similar VAF. Another problem on performing PCAs with variations in the data is the possibility of rotation indeterminacy when the PCA solution of a resampled data presents with a different rotation of the original PCA solution. These issues generate artificially biased bootstrapped distributions, potentially invalidating the procedure (Timmerman et al., 2007; Zientek and Thompson, 2007). We have implemented a step of procrustes rotation between the original loadings (target) and the bootstrapped sample, as has been previously demonstrated to be a reasonable method to deal with such issues (Timmerman et al., 2007; Zientek and Thompson, 2007). The Procrustes rotation is obtained by the procrustes() function from the pracma package. The algorithm for bootstrapping the PCA solutions is represented in Figure 4 and implemented in the utility function boot_pca_sample(). The number of bootstrap samples is set to 1000 by default. The user must be careful on setting the number too low, reducing the performance of the approximation (Efron, 1987). However, setting the number of bootstrap samples too high might unnecessarily increase computing time with little gain (Figure 4—figure supplement 1).
Indexes of component similarity
We have included several component similarity indexes for determining component/factor invariance in syndRomics. The function component_similiarity () returns the specified similarity metrics as well as their summary statistics (average and standard deviation, if applicable) from a list of loading matrices (load.list). The argument s_cut_off is used in the calculation of the Cattell’s s index (see below) and ndim is used to limit the number of components from which to compute the indexes from. Each index has been programmed in a separate utility function for convenience. Although they are not meant to be manually called, users can call them to calculate any of these metrics for a given set of two component loadings. The similarity_metric argument takes a single character or a vector of characters to specify which metrics to compute. These can be: ‘cc_index’, ‘r_correlation’, ‘rmse’ and/or ‘s_index’. The user can also specify ‘all’ to get all metrics. Their definitions are documented below.
Congruence coefficient (CC, ‘cc_index’)
Request a detailed protocolFirst suggested by Burt, 1948, it was popularized by Tucker, 1951 and therefore is also known as Tucker’s congruence coefficient. It is calculated as (2):
where ${x}_{i}$ and ${y}_{i}$ are the loadings of the variable $i$ on the component or factor $x$ and $y$ respectively. $\varnothing \left(x,y\right)$ is equivalent to the cosine of the angle between two vectors and is also referred to as the cosine similarity metric. CC is a measure of proportional similarity between two components, and technically the index has a range from 1 (perfect negative congruence) to 1 (perfect positive congruence). In practice, because the all the loadings of a PC can be multiplied by 1 without changing the meaning of the PC, the absolute value of CC is considered, which correspondingly ranges from 0 to 1. The closer to 1, the more similar the two components are. Chan et al. discussed the 0.9 rule of thumb as an indicator of good matching between PCs (Chan et al., 1999). The application of CC as a similarity metric for factor invariance has been extensively studied (Chan et al., 1999; LorenzoSeva and ten Berge, 2006).
Pearson’s correlation coefficient (r, ‘r_correlation’)
Request a detailed protocolThe calculation of r between two vectors of component loadings has also been used as a pattern matching metric (Guadagnoli and Velicer, 1991). It is computed as (3):
In the syndRomics package, the Pearson’s correlation coefficient is calculated using the cor() function of the stats package.
Root mean square error (RMSE, ‘rmse’)
Request a detailed protocolRMSE has also been used as a metric for factor matching (Guadagnoli and Velicer, 1991). It is calculated as the square root of the average squared difference of the loadings of the variables as (4):
Component 2  

Component 1  PS  H  NS  
PS  ${f}_{11}$  ${f}_{12}$  ${f}_{13}$  
H  ${f}_{21}$  ${f}_{22}$  ${f}_{23}$  
NS  ${f}_{31}$  ${f}_{32}$  ${f}_{33}$ 
where $n$ is the number of variables in both components $x$ and $y$. A RMSE of 0 determines a perfect matching, and therefore the smaller the RMSE is, the more equivalent the two components $x$ and $y$ are.
Cattell’s s index (‘s_index’)
Request a detailed protocolThe s index was first suggested by Cattell and Baggaley, 1960; Cattell et al., 1969. It is based on the factor mandate matrix (Cattell and Baggaley, 1960) where loadings are either one if a component is considered to act on a variable, called a salient variable, or 0 if not (forming the hyperplane space). Cattell’s suggested an arbitrary ± 0.1 cutoff where variables with loadings outside the cutoff range are removed from the hyperplane and considered to be salient variables. In practice, one might want to alter the threshold depending on the experimental conditions. Any loading inside the cutoff range is then interpreted as having been produced by chance. The s index is calculated from the crossclassification of the common variables of two components/factors:
where PS = positive salient variable; H = hyperplane variable; NS = negative salient variable; ${f}_{ij}$ is the joint frequency. Positive and negative salient variables are variables outside the cutoff range with positive or negative loadings respectively.
Pattern matching is determined by comparing the cell frequencies in the crossclassification table. Here we implement the simplified form of calculating s (5):
The reader is referred to Cattell and Baggaley, 1960; Cattell et al., 1969; Guadagnoli and Velicer, 1991 for details on reasoning and calculations. s ranges from 1 (perfect similarity) to −1 (perfect dissimilarity) centered at 0 (pattern due to chance). Similar to CC, the absolute value of s is considered.
Internal functions
There are internal functions used by the package that the user might never have to call directly, although they are accessible in case the user needs them. Here, we provided a general description of those, leaving the details to the package documentation. All the internal functions to extract similarity metrics are: extract_cc(), extract_s() and extract_rmse(). They all take two numeric vectors and return the corresponding similarity metric between them.
new_syndromics()
Request a detailed protocolHelper function to construct the ‘syndromics’ class object that will be use in the S3 generic and method functions. It returns an object of class ‘syndromics’ of the type list.
stand_loadings()
Request a detailed protocolThis function extracts the standardized loadings from the output of the prcomp() or the princals() functions. In the case of the prcomp() solution, the standardized loadings are calculated as: $s.loadings=eigenvectors\times \sqrt{eigenvalues}$ if the PCA was performed on the standardized (scaled to unit variance) data or $s.loadings=(eigenvector\times \sqrt{eigenvalues})/S$ where $S$ is the vector of the variables standard deviation. In the case of princals(), standardized loadings are returned directly in its output and therefore stand_loadings() returns those. The function returns a data frame with the standardized loadings in the form of variables as rows and PCs as columns.
extract_loadings()
Request a detailed protocolThis is a wrapper function for stand_loadings() with added functionalities such as error breakers that is used by most functions in the package.
extract_syndromic_plot()
Request a detailed protocolThis function is internally called by the syndromic_plot() function and returns a ggplot2 object with the syndromic plot for the specified PC. The only argument that is not present in the syndromic_plot() function is the pc argument that specifies which PC to plot. Users should always use syndromic_plot() function instead of extract_syndromic_plot() since syndromic_plot() automatically incorporates other functionalities.
component_similarity()
Request a detailed protocolThis function is called by the pc_stability() function to calculate the specified similarity metric (see above) between the given list of data frames of loadings. While pc_similarity() uses this function to calculate similarity between the original (parent) loadings and a B sample loadings, the passed list of loadings can be n > 2. Then, the similarity metrics will be calculated between all combinations of n. It returns a list of objects containing a list of the comparisons, a data frame with the averaged metric and the bounds of confidence interval for each specified metric and PC.
boot_pca_sample()
Request a detailed protocolThis generic function is passed to the statistic argument of the boot() function internally called by the pc_stability() function. It implements the bootstrapping algorithm described above (Figure 2A). Then the boot() function will call boot_pca_sample() B times from the specified data and the pca output of the prcomp() (through the method boot_pca_sample.prcomp()) or princals() (through the method boot_pca_sample.princals()) function, returning a list of B data frames of loadings. The bootstrapping method can be specified using the sim argument.
permut_pca_D() or permut_pca_V()
Request a detailed protocolThis is a generic function internally called by permut_pca_test() to produce P permutations of the given output of the prcomp() or the princals() functions using permD or permV method. Four S3 R function methods are implemented: permut_pca_D.prcomp(), permut_pca_D.princals(), permut_pca_V.prcomp(), permut_pca_V.princals(). It returns a list of the results of permuting the data, conducting a PCA and extracting either the VAF or the standardized loadings for each P as in Figure 2.
Plot.syndromics()
Request a detailed protocolThis function implement the S3 method for plotting ‘syndromics’ class objects generated by pc_stability() and permut_pc_test() functions using the R generic plot(). It returns specific plots calling the visualization functions implemented in the package.
Nonlinear PCA
Request a detailed protocolNonlinear PCA by optimal scaling and alternating least square was obtained using the princals() function from the {Gifi} package in R. We specified to analyze all variables with nominal restriction scaling, allowing for nonmonotonic transformations, and set a restriction of 3 degrees in polynomial transformations for nonlinearity. The corresponding instruction was: princals(nlpca_data, ndim = ncol(nlpca_data), ordinal = FALSE, degrees = 3, knots = knotsGifi(nlpca_data, type=‘E’)), where nlpca_data is the imputed dataset for case study 2 (see supplementary code for more details).
Missing data analysis
Request a detailed protocolDetails on the code are available as supplementary material. Data wrangling for the two case studies was performed using R packages included in the Tidyverse package. Missing pattern visualization were obtained using the naniar (Tierney et al., 2020) R packages. Test for MCAR was performed using the TestMCARNormality() function from the MissMech package (Jamshidian et al., 2014). Multiple imputation was performed using predicting mean matching method available in the mice (Buuren and GroothuisOudshoorn, 2011) R package, setting the number of imputations to m = 50. A list of 50 complete datasets were then obtained and processed by PCA as specified in the main text. For each m dataset, the loadings where extracted and rotated using Procrustes rotation (pracma package) toward the average of the imputed datasets. The distributions of loadings and component similarities for the first three PCs where calculated using the syndRomics package as described above.
Data availability
This work used already available data at the Open Data Commons for Spinal Cord Injury (http://odcsci.org/) and Plos One (https://doi.org/10.1371/journal.pone.0169490).

Open Data Common for Spinal Cord InjuryID (ODCSCI:26). Cervical (C5), unilateral spinal cord injury with diverse 740 injury modalities, multiple behavioral outcomes, and histopathology.

journalUncovering precision phenotypebiomarker associations in traumatic brain injury using topological data analysis.https://doi.org/10.1371/journal.pone.0169490
References

Principal component analysisWiley Interdisciplinary Reviews: Computational Statistics 2:433–459.https://doi.org/10.1002/wics.101

Bootstrap based confidence limits in principal component analysis — A case studyChemometrics and Intelligent Laboratory Systems 120:97–105.https://doi.org/10.1016/j.chemolab.2012.10.007

Controlling the false discovery rate: a practical and powerful approach to multiple testingJournal of the Royal Statistical Society 57:289–300.https://doi.org/10.1111/j.25176161.1995.tb02031.x

SoftwarePracma: Practical Numerical Math Functions, version 2.2.9Pracma: Practical Numerical Math Functions.

Remarks on parallel analysisMultivariate Behavioral Research 27:509–540.https://doi.org/10.1207/s15327906mbr2704_2

The factorial study of temperamental traitsBritish Journal of Statistical Psychology 1:178–203.https://doi.org/10.1111/j.20448317.1948.tb00236.x

Mice : Multivariate Imputation by Chained Equations in RJournal of Statistical Software 45:1–67.https://doi.org/10.18637/jss.v045.i03

Developing a data sharing community for spinal cord injury researchExperimental Neurology 295:135–143.https://doi.org/10.1016/j.expneurol.2017.05.012

The scree test for the number of factorsMultivariate Behavioral Research 1:245–276.https://doi.org/10.1207/s15327906mbr0102_10

Factor matching procedures: an improvement of the s index; with tablesEducational and Psychological Measurement 29:781–792.https://doi.org/10.1177/001316446902900405

The salient variable similarity index for factor matchingBritish Journal of Statistical Psychology 13:33–46.https://doi.org/10.1111/j.20448317.1960.tb00037.x

BookBootstrap Methods and Their ApplicationCambridge, United KIngdom: Cambridge University Press.

Multiple comparisons among meansJournal of the American Statistical Association 56:52–64.https://doi.org/10.1080/01621459.1961.10482090

Better bootstrap confidence intervalsJournal of the American Statistical Association 82:171–185.https://doi.org/10.1080/01621459.1987.10478410

SoftwareGgnewscale: Multiple Fill and Colour Scales in “Ggplot2”, version 0.4.1Ggnewscale: Multiple Fill and Colour Scales in “Ggplot2”.

BookAn Introduction to Applied Multivariate Analysis with RBerlin, Germany: SpringerVerlag.

Posttraumatic brain injury cognitive performance is moderated by variation within ANKK1 and DRD2 genesJournal of Head Trauma Rehabilitation 30:E54–E66.https://doi.org/10.1097/HTR.0000000000000118

Syndromics: a bioinformatics approach for neurotrauma researchTranslational Stroke Research 2:438–454.https://doi.org/10.1007/s1297501101211

DataCervical (C5), unilateral spinal cord injury with diverse injury modalities, multiple behavioral outcomes, and histopathologyOpen Data Common for Spinal Cord Injury.https://doi.org/10.7295/W9T72FMZ

An improvement on Horn's Parallel Analysis Methodology for Selecting the Correct Number of Factors to RetainEducational and Psychological Measurement 55:377–393.https://doi.org/10.1177/0013164495055003002

Relation of sample size to the stability of component patternsPsychological Bulletin 103:265–275.https://doi.org/10.1037/00332909.103.2.265

A comparison of pattern matching indicesMultivariate Behavioral Research 26:323–343.https://doi.org/10.1207/s15327906mbr2602_7

Some necessary conditions for commonfactor analysisPsychometrika 19:149–161.https://doi.org/10.1007/BF02289162

Multivariate analysis of MRI biomarkers for predicting neurologic impairment in cervical spinal cord injuryAmerican Journal of Neuroradiology 38:648–655.https://doi.org/10.3174/ajnr.A5021

SoftwareRlang: Functions for Base Types and Core R and “Tidyverse” Features, version 0.4.4Rlang: Functions for Base Types and Core R and “Tidyverse” Features.

Missing data in clinical research: an integrated approachBritish Journal of Dermatology 177:1463–1465.https://doi.org/10.1111/bjd.16010

Bootstrap scree tests: a monte carlo simulation and applications to published dataBritish Journal of Mathematical and Statistical Psychology 59:35–57.https://doi.org/10.1348/000711005X66770

Analysis of a complex of statistical variables into principal componentsJournal of Educational Psychology 24:417–441.https://doi.org/10.1037/h0071325

Neurotrauma as a bigdata problemCurrent Opinion in Neurology 31:702–708.https://doi.org/10.1097/WCO.0000000000000614

MissMech : AnR package for testing homoscedasticity, multivariate Normality, and missing completely at random (MCAR)Journal of Statistical Software 56:1–31.https://doi.org/10.18637/jss.v056.i06

ConferenceBias in principal components analysis due to correlated observationsConference on Applied Statistics in Agriculture.https://doi.org/10.4148/24757772.1247

Principal component analysis: a review and recent developmentsPhilosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 374:20150202.https://doi.org/10.1098/rsta.2015.0202

The Application of Electronic Computers to Factor AnalysisEducational and Psychological Measurement 20:141–151.https://doi.org/10.1177/001316446002000116

Missing data in clinical trials: pitfalls and remediesInternational Journal of Applied & Basic Medical Research 4:S6–S7.

Principal component analysis for designed experimentsBMC Bioinformatics 16 Suppl 18:S7.https://doi.org/10.1186/1471210516S18S7

A principal component analysis of coagulation after traumaJournal of Trauma and Acute Care Surgery 74:1223–1230.https://doi.org/10.1097/TA.0b013e31828b7fa1

Permutationvalidated principal components analysis of microarray dataGenome Biology 3:research0019.1.https://doi.org/10.1186/gb200234research0019

FactoMineR : An R package for multivariate analysisJournal of Statistical Software 25:1–18.https://doi.org/10.18637/jss.v025.i01

SoftwareDoctoral Thesis: Nonparametric Inference in Nonlinear Principal Components Analysis: Exploration and BeyondDoctoral Thesis: Nonparametric Inference in Nonlinear Principal Components Analysis: Exploration and Beyond.

Nonlinear principal components analysis: introduction and applicationPsychological Methods 12:336–358.https://doi.org/10.1037/1082989X.12.3.336

SoftwareGifi: Multivariate Analysis with Optimal Scaling, version 0.39Gifi: Multivariate Analysis with Optimal Scaling.

Ten quick tips for effective dimensionality reductionPLOS Computational Biology 15:e1006907.https://doi.org/10.1371/journal.pcbi.1006907

Development of a database for translational spinal cord injury researchJournal of Neurotrauma 31:1789–1799.https://doi.org/10.1089/neu.2014.3399

SoftwareFactor Analysis as a Tool for Pattern Recognition in Biomedical Research; a Review with Application in R SoftwareFactor Analysis as a Tool for Pattern Recognition in Biomedical Research; a Review with Application in R Software.

SoftwareR: A Language and Environment for Statistical ComputingR Foundation for Statistical Computing, Vienna, Austria.

SoftwarePsych: Procedures for Personality and Psychological ResearchPsych: Procedures for Personality and Psychological Research.

Extensive spontaneous plasticity of corticospinal projections after primate spinal cord injuryNature Neuroscience 13:1505–1510.https://doi.org/10.1038/nn.2691

Missing data: Our view of the state of the artPsychological Methods 7:147–177.https://doi.org/10.1037/1082989X.7.2.147

SoftwareGgrepel: Automatically Position NonOverlapping Text Labels with Ggplot2Ggrepel: Automatically Position NonOverlapping Text Labels with Ggplot2.

SoftwareRStudio: Integrated Development for R. RStudio, IncRStudio: Integrated Development for R. RStudio, Inc.

SoftwareNaniar: Data Structures, Summaries, and Visualisations for Missing Data (0.5.0)Naniar: Data Structures, Summaries, and Visualisations for Missing Data (0.5.0).

Estimating confidence intervals for principal component loadings: A comparison between the bootstrap and asymptotic resultsBritish Journal of Mathematical and Statistical Psychology 60:295–314.https://doi.org/10.1348/000711006X109636

Softwarepng: Read and write PNG images, version R package version 0.17png: Read and write PNG images.

Using Generalized Procrustes Analysis for Multiple Imputation in Principal Component AnalysisJournal of Classification 31:242–269.https://doi.org/10.1007/s003570149154y

Softwaredplyr: A Grammar of Data Manipulation, version R package version 0.8.4dplyr: A Grammar of Data Manipulation.

Softwarestringr: Simple, Consistent Wrappers for Common String Operationsstringr: Simple, Consistent Wrappers for Common String Operations.

Principal components analysis in clinical studiesAnnals of Translational Medicine 5:351.https://doi.org/10.21037/atm.2017.07.12

Applying the bootstrap to the multivariate case: Bootstrap component/factor analysisBehavior Research Methods 39:318–325.https://doi.org/10.3758/BF03193163

Comparison of five rules for determining the number of components to retainPsychological Bulletin 99:432–442.https://doi.org/10.1037/00332909.99.3.432
Article and author information
Author details
Funding
National Institutes of Health (NS106899)
 Adam R Ferguson
National Institutes of Health (NS088475)
 Adam R Ferguson
Department of Veterans Affairs (I01RX02245)
 Adam R Ferguson
Department of Veterans Affairs (I01RX002787)
 Adam R Ferguson
Craig H. Neilsen Foundation (Special Project)
 Adam R Ferguson
Wings for Life (Special Project)
 Adam R Ferguson
Wings for Life (Individual Grant)
 Abel TorresEspín
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This work is supported by NIH grants NS106899 (ARF), NS088475 (ARF); VA Grants 1I01R × 002245 (ARF) and I01R × 002787 (ARF)
Version history
 Received: August 5, 2020
 Accepted: January 13, 2021
 Accepted Manuscript published: January 14, 2021 (version 1)
 Version of Record published: February 2, 2021 (version 2)
Copyright
© 2021, TorresEspín et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,966
 views

 233
 downloads

 21
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Cell Biology
 Computational and Systems Biology
Bats have unique characteristics compared to other mammals, including increased longevity and higher resistance to cancer and infectious disease. While previous studies have analyzed the metabolic requirements for flight, it is still unclear how bat metabolism supports these unique features, and no study has integrated metabolomics, transcriptomics, and proteomics to characterize bat metabolism. In this work, we performed a multiomics data analysis using a computational model of metabolic fluxes to identify fundamental differences in central metabolism between primary lung fibroblast cell lines from the black flying fox fruit bat (Pteropus alecto) and human. Bat cells showed higher expression levels of Complex I components of electron transport chain (ETC), but, remarkably, a lower rate of oxygen consumption. Computational modeling interpreted these results as indicating that Complex II activity may be low or reversed, similar to an ischemic state. An ischemiclike state of bats was also supported by decreased levels of central metabolites and increased ratios of succinate to fumarate in bat cells. Ischemic states tend to produce reactive oxygen species (ROS), which would be incompatible with the longevity of bats. However, bat cells had higher antioxidant reservoirs (higher total glutathione and higher ratio of NADPH to NADP) despite higher mitochondrial ROS levels. In addition, bat cells were more resistant to glucose deprivation and had increased resistance to ferroptosis, one of the characteristics of which is oxidative stress. Thus, our studies revealed distinct differences in the ETC regulation and metabolic stress responses between human and bat cells.

 Computational and Systems Biology
 Neuroscience
Normal aging leads to myelin alterations in the rhesus monkey dorsolateral prefrontal cortex (dlPFC), which are positively correlated with degree of cognitive impairment. It is hypothesized that remyelination with shorter and thinner myelin sheaths partially compensates for myelin degradation, but computational modeling has not yet explored these two phenomena together systematically. Here, we used a twopronged modeling approach to determine how agerelated myelin changes affect a core cognitive function: spatial working memory. First, we built a multicompartment pyramidal neuron model fit to monkey dlPFC empirical data, with an axon including myelinated segments having paranodes, juxtaparanodes, internodes, and tight junctions. This model was used to quantify conduction velocity (CV) changes and action potential (AP) failures after demyelination and subsequent remyelination. Next, we incorporated the single neuron results into a spiking neural network model of working memory. While complete remyelination nearly recovered axonal transmission and network function to unperturbed levels, our models predict that biologically plausible levels of myelin dystrophy, if uncompensated by other factors, can account for substantial working memory impairment with aging. The present computational study unites empirical data from ultrastructure up to behavior during normal aging, and has broader implications for many demyelinating conditions, such as multiple sclerosis or schizophrenia.