Introduction
Compositional data comprise vectors of quantitative measures that describe the constituent parts of a whole. In compositional data, there is relevant information found in the ratios between the different parts of the vector [ 1 ]. Data arising from various -omics platforms are compositional in nature—for instance, RNA sequencing (RNA-seq), 16S amplicon sequencing, and single-cell RNA-seq are all examples of platforms that yield compositional data. The primary reason for the compositional classification of these platforms is that the number of sequences varies across samples due to technical artifacts in the sequencing process, meaning that within-sample counts must either be normalized to a unit sum or log-ratio transformed in order to effect any kind of meaningful interpretation [ 2 ]. Some microbiome studies have relied onrarefying such data, which involves resampling the counts to a constant read depth across all samples. However, it has been recognized that such a procedure essentially amounts to discarding valid data and should therefore be discouraged [ 3 ].
This interpretability problem is particularly worrisome when considering associations between features—the prime example being Pearson’s correlation. Correlations on raw counts have no meaningful interpretation. Similarly, when normalizing the counts to proportions, a negative bias is introduced in the correlations as a result of the unit-sum constraint.
A vector containing non-negative components is a D-part composition if there is relevant information found in the ratios between its non-zero elements [ 4 ]. Two vectors and are compositionally equivalent if for some constant k> 0 . This property thus defines a set of equivalence classes for compositional vectors. As such, we proceed to describe compositions that have been subject to the closure operation: . We refer to compositional vectors having undergone the closure operation asproportions.
An appropriate class of measures based on proportionality has been established to address the issue of quantifying association in compositional data. The most basic such measure is the variance of the log-ratio of two proportions. Suppose is a random vector for individual i in an independent, identically distributed random sample for i=1,… ,n containing proportions such that and . Though compositional data theory in general includes the possibility of proportions equal to zero or one, in this article we exclude these two cases as the metrics of proportionality under consideration are only defined for proportions strictly between zero and one. Note that this does not preclude zero counts appearing in the data (i.e. sampling zeros); it only rules out zeros in the true, unobserved compositions. The probability distribution on is
where is some distribution on the simplex (e.g. multivariate logit-normal, Dirichlet, etc.) defined by a parameter vector . The variance of the log-ratio between proportions j and k is given by
When is exactly proportional to , we have that . Hence, the smaller the value of , the stronger the association between the log-ratio abundances of features j and k . Note that v is non-negative, and thus this metric does not give information about the direction of the associations of the pairs of features.
Note that we calculate the remaining metrics of proportionality on the log-ratio transformed proportions; e.g. we could use the additive log-ratio (ALR) transformation with reference D , given by
Alternatively, one could use the centered log-ratio (CLR) transformation, given by
where, is the geometric mean of the proportions for sample i . The authors in [ 5 , 6 ] emphasized using log-ratio transformed proportions in proportionality metrics in order to properly define an accordance between the relative and underlying absolute abundances.
One important issue with is that it lacks an interpretable scale. To address this issue, [ 5 , 6 ] developed alternative measures of proportionality.
The first is called ϕ , and is simply scaled by the variance of the logarithm of its first argument. In this case, we must keep in mind the choice of log-ratio transformation. For the ALR transformation with reference D , we denote ϕ as:
Similarly, for the CLR transformation, we have
We have that and . The interpretation is the same as for v , only now it is scaled based on its first argument. The main downside of this metric is that it is not symmetric with regards to its arguments. A final proportionality metric, ρ is defined, assuming ALR transformation with reference D as:
For the CLR transformation, we have
with and . If , then the log-ratio transformed and are strongly proportional, and if , then they are strongly inversely proportional. There exists another metric of proportionality commonly referred to as , which is a symmetrized version of ϕ . There is, however, a one-to-one correspondence between ρ and ; as such, since ρ admits a more interpretable scale, we will not consider in this article.
Metrics of proportionality recently gained more recognition when [ 7 ] suggested their use to detect gene–gene associations in single-cell sequencing experiments. This publication compared numerous measures of association across multiple datasets by evaluating their ability to create functionally coherent single-cell gene co-expression networks. There was a clear advantage for the proportionality metrics across the datasets. Additional work has been done in recent years on proportionality metrics in subcompositions [ 8 ], in regularization of compositional covariance [ 9 ], and in variable selection in log-ratios [ 10 ].
The one common property of these kinds of metrics is that they arescale invariant, which is the main property allowing their use in compositional data. However, classical techniques for the analysis of compositional data were developed with continuous measures in mind; e.g. the constituent parts of a soil sample. Conversely, sequencing platforms provide data in the form of counts. A number of recent studies have shed light on problems that arise when applying traditional methods for compositional data analysis on count-based compositional data. It was shown by [ 11 ] that proportionality cannot be exactly represented in count-based compositional data, which they referred to as lattice compositional data. [ 12 ] showed a similar phenomenon in logarithmic transformations in single-cell sequencing; varying read depths among cells led to systematic errors and spurious differences in expression. Relatedly, [ 13 ] presented novel modeling techniques in light of these issues in count data.
In addition to the representation problems outlined in [ 11 ], we further claim that there is in fact a bias that arises when applying proportionality metrics on count-based data, which is a result of added variability from the randomness of the sequencing process. In this paper, we quantify this bias and provide an alternative model-based means of estimating proportionality metrics in count-based compositional data.
Materials and methods
Bias
In this section, we outline the inherent bias in estimating metrics of proportionality using observed counts. Suppose that we have a matrix containing counts of D features (e.g. genes, species, etc.) from n samples. Let the entries of be , representing the observed count of feature j in sample i . The total number of reads for sample i is denoted by . We seek to compare the result of plugging counts (or equivalently observed proportions ) into v , ϕ , and ρ against the true (unknown) values which instead takeunobserved proportions as arguments to these metrics. In each case, there is no exact closed form for the bias, and we instead use an approximation based on a second-order Taylor series expansion. We summarize the results here, but full derivations can be found in Section S2. In each case, it is clear that the bias is more pronounced when the read depths are small, on average. However, it is also evident that the amount of variation in the read depths itself plays into the bias, with more variability leading to a larger discrepancy between the empirical and true values of the metrics. The underlying distribution of the unobserved proportion vector is also relevant. We also note that the bias takes on a different form depending on the choice of log-ratio transformation.
First, we consider . This particular metric is invariant with respect to the choice of the initial log-ratio transformations we consider. This metric applied to observed counts is approximately related to the same metric applied to the unobserved true proportions as follows:
where,
As expected, the expectation of is an important factor; however, the expectations of and are also present, suggesting that the discrepancy will be more pronounced in less abundant features.
Unlike v , the metric ϕ is sensitive to the choice of log-ratio transformation, and thus the bias results are shown separately for the ALR and CLR cases. First, if we consider the ALR transformation with reference D , the bias can be expressed as:
Likewise, for CLR the bias is given by
where
Finally, we consider the analogous result for :
The bias for CLR is given by
where,
We can once again see that the bias is less pronounced when the mean of the read depths is large. However, the variation of , as well as the distribution of the proportions is also important. This suggests that there could still be bias present for estimates for certain feature pairs even when the read depths are large.
To numerically investigate these biases, we run a small simulation from a multinomial logit-normal distribution (see the next section for the definition of this distribution). We simulate D=7 features on n=2000 samples, with 1000 simulation replications. In each replication, the estimated proportionality metrics on the raw counts , , and are calculated (after imputing zeros using Bayesian-multiplicative replacement [ 14 ], if necessary) and compared to their respective true values , , and . Additionally, a bias-corrected version of the estimator is calculated by adjusting the empirical estimates using equations 5 , 7 , 8 , 10 , and 11 .
Figure 1 shows the difference between the estimated values and the true values for a single element of the matrix. It is evident that all of v , ϕ , and ρ exhibit strong biases when calculated on the raw counts. In all cases, the sampling distribution of the corrected estimator is approximately centered around zero. Results for the full matrix can be seen in Supplementary Figs S1 – S5 .
It should be noted that this bias-correction approach is not feasible in practice, as it would require knowledge of the true population parameters, which happen to be known in this simulation. As an alternative, we propose directly modeling the distribution of the proportions . The parameters of this distribution then define the values of v , ϕ , and ρ , allowing the proportionality metrics to be estimated free from the variation induced by the sequencing process itself.
Model
We proceed in the context of the multinomial logit-normal model, posited by [ 15 ]. This model was revisited by [ 16 ] to impose a penalty term to the covariance matrix on the ALR scale; here we extend the model to penalize either on the ALR or CLR scale. This, in turn, allows penalized estimation of metrics of proportionality on either scale. The logit-normal multinomial model has been extremely popular in statistical modeling of the microbiome [ 17 – 19 ]. In this formulation, we assume that the count vector for sample i , denoted by , is distributed as
where D is the total number of features observed among all samples, and is the read depth. The proportion vector for individual i is assumed to follow the logit-normal distribution, which is characterized by the inverse ALR transformation
where is the closure operation. Each is an unobserved, latent vector assumed to be multivariate-normal:
We set the feature as the reference feature in the ALR transformation (without loss of generality; any feature could be moved to that column). This implicitly prescribes a multivariate logit-normal distribution on the proportion vectors . The read depths themselves are random variables on which we assume a log-normal distribution
The log-normal distribution is appropriate given that the read depth distribution in a dataset can span multiple orders of magnitude. This assumption has been used in previous literature [ 20 ].
To handle the possibility of a large number of features D relative to the sample size n , we apply the Graphical Lasso (GLasso) penalty [ 21 ] to the multivariate normal log-likelihood for the vectors:
where denotes the -norm, and is the penalty parameter. The penalty in this case is applied to the precision matrix on the ALR scale. This penalization is very important in guarding against spurious associations when the number of features is large. Alternatively, if we are interested in estimating proportionality metrics on the CLR scale, we can instead include the penalty term , where , denotes the D-1 identity matrix, denotes a D-1 vector of ones, and .
Estimating proportionality metrics
In the logit-normal framework, the elements have a convenient form. That is, if the diagonal elements of are for j=1,… ,D and the off-diagonal elements are for j≠ k , then the v elements can be written as
This formula will apply to both ALR- and CLR-transformed proportions. For ϕ and ρ in the ALR case, we have
For the CLR case we have
To fit the logit-normal multinomial model with Graphical Lasso penalty applied to either the ALR- or CLR-scale precision matrix, we perform maximum likelihood estimation through a modification of the EM algorithm posited in [ 16 , 22 ]. Zero counts can be input directly into the model without imputation; the corresponding elements of each are estimated in the E-step of the algorithm. Model selection (through the choice of tuning parameter λ ) is performed using the Extended Bayesian Information Criterion [ 23 ]. The steps of the algorithm can be found in Section S4, equations 14 – 17 .
Simulation study
We performed an extensive simulation study to compare the performance of model-based versus empirical estimates. We focus on ρ in the simulations since it is the most interpretable metric, and we can set it to zero for certain feature pairs, allowing precision-recall analysis (whereas the null for v and ϕ is infinity).
Our simulation parameters are based on an initial multinomial logit-normal model fitted to single-cell RNA-seq data from a mouse embryonic stem-cell (mESC) dataset from Buettner et al . [ 24 ]. The model parameters , , , and , are estimated from the data and are then used as the ground truth in the simulations in order to ensure the simulated data are realistic. The read depths are simulated from . We run two separate simulation scenarios; one assuming is sparse and another assuming is sparse. See Section S5 for more details about the simulation procedures.
For both ALR and CLR simulation studies, we consider multiple scenarios. We simulate datasets under every combination of sample sizes n=100,200,500 , feature sizes D=50,200,200 , and sparsity proportions in ρ (on either ALR or CLR scale) as 0.15, 0.5, and 0.85. We ran 50 replications of each simulation scenario.
Once we have obtained working values of all the parameters, we modify the true covariance matrices on either the ALR or CLR scale to be sparse. This will, in turn, allow some of the true and to be zero. For the ALR transformation, we let be a sparse matrix. Likewise, for the CLR transformation, we set values in the CLR covariance matrix to zero and then transform back to the ALR scale to simulate the data. This lets us ascertain the ability of the model-based and empirical estimates to find the non-zero elements of ρ . Since the empirical estimator will not set values exactly to zero, we declare each feature pair association as significant via the test for Pearson’s correlation on the appropriate log-ratio scale. A zero correlation on that scale is equivalent to a zero value of ρ calculated from the same transformation. This procedure has precedent in the literature [ 25 ]. The performance measures for ρ we consider are (i) F-score, which combines precision and recall into a single measure, and (ii) mean-square error (MSE). A higher value of the F -score and a lower value of MSE correspond to better performance for ρ estimation.
Data application
To demonstrate the applicability of the proportionality metric estimation technique introduced in this paper, we apply the methods to two datasets—one is a single-cell RNA-seq dataset, and the other is a metagenomic dataset. In both cases, we focus on estimating ρ on the ALR scale.
The first dataset is a mESC single-cell RNA-seq dataset from [ 24 ]. Flow cytometry was used to sort by cell-type, with cells sorted into G1, S, and G2M cell cycle stages. The original dataset contained 96 samples and 38 390 genes. The goal of our analysis is to focus on genes that were previously shown to be differentially expressed among cell types. In this analysis, we consider 570 genes associated with cell cycle based on GO annotations; this list was provided in the original publication. We took a further subset of genes by filtering out genes having >20% dropouts in at least one of the cell-cycle stages. There were 303 genes remaining after filtering. This analysis allows investigation of whether gene–gene associations differ between cell-cycle stages in a subset of genes that were previously known to be related to cell cycle.
The second dataset comprises shotgun metagenomic sequencing information derived from stool samples procured from participants with pediatric-onset (symptom onset < 18 years of age) multiple sclerosis (MS) and unaffected controls [ 26 – 28 ]. The goal of the study was to determine how suppression of dietary fibre fermentation can induce inflammation in MS. Data were collected through a larger network called the Canadian Pediatric Demyelinating Disease Network; in the data specifically collected for the dietary fibre study, all participants were under 22 years of age at the time of stool sample procurement. There were 30 MS participants (24 female) and 52 unaffected controls (36 female) who provided a stool sample. For our analyses, we proceed at the genus level; initially there were 2101 genera available in the dataset. We filtered out genera with >10% zeros, leaving 204 genera. This facilitated more stable estimates of the model parameters in the small sample size case. The analysis focuses on which genus pairs have differing associations between MS participants and unaffected controls.
Results
Simulation results
We now show the results from the simulations considering sparsity among the ρ values. ALR results can be seen in Fig. 2 , where we compare the distributions of the F -score and MSE for the different estimates. In the majority of scenarios, the model-based estimates outperform the empirical estimates in terms of F -score. In a small number of cases, the empirical estimates are comparable or marginally better; this may be due to over-penalization in the Graphical LASSO step. In all cases, the MSEs for model-based estimates are smaller than those of the empirical estimates. In the CLR case, shown in Fig. 3 , we again see an improvement in F -score for the majority of scenarios, though there may again be over-penalization for the model-based estimates in the high sparsity case when D=50 . The MSE is smaller for model-based in the vast majority of scenarios, with a few exceptions in cases with a smaller number of features.
This simulation study shows strong evidence of the superiority of the model-based estimates over the empirical estimates for ρ elements. Thus, the model-based estimates can be used to guard against spurious correlations when considering the proportionality metric ρ .
Data application results
In this section, we outline the results from the two data analyses described in the “Materials and methods” section. First we discuss the results from the murine single-cell RNA-seq dataset. All gene descriptions were obtained from GeneCards [ 29 ].
We consider the estimated ρ values for each of the pairs of the 303 genes retained after the aforementioned filtering steps. Estimates were obtained using both the empirical estimator and the model-based estimator. We compare the top 100 gene–gene associations (in absolute value) detected using the empirical and model-based estimators. In the G1 phase, there were 56 gene pairs in common, in G2M phase, there were 38 gene pairs in common, and in S phase, there were 40 gene pairs in common ( Supplementary Fig. S6 ). Relatedly, the empirical and model-based estimates are compared in Fig. 4 . This shows the attenuation toward the null in all three subgroups for the model-based estimates, which is an important guard against spurious correlations. These comparisons between the empirical and model-based estimates show how different the results can be using the empirical estimator for ρ . In the S stage analysis, there is more shrinkage of ρ for certain gene pairs. This is likely due to a higher proportion of zeros in the S counts (11.5%) compared to the G1 and G2M stages (4.7% and 2.7%, respectively).
The top 100 gene pairs within each cell stage are shown in Supplementary Tables S4 – S6 . Additionally, Table 1 summarizes the top differences (in absolute value) in ρ values between phases G1 and G2M; more extensive tables showing the top 100 gene pairs for each of the cell phase comparisons are shown in Supplementary Tables S1 – S3 . In the G1 versus G2M comparison, the top 25 ρ differences all involve CINP (cyclin dependent kinase 2 interacting protein), which is part of the DNA replication complex and binds chromatin in the G1 phase. The largest ρ difference is seen between genes KIF20B and CINP, with . KIF20B is involved in cellular proliferation and regulation of cytokinesis. Silencing KFI20B has also been shown to increase cell arrest in G1 phase[ 30 ].
Table 1.
Top 25 gene pairs whose ρ value differs between G1 and G2M stages
| Gene 1 | Gene 2 | |||
| KIF20B | CINP | 0.4516 | −0.0832 | 0.5349 |
| TIPIN | CINP | 0.3992 | −0.1094 | 0.5085 |
| BUB3 | CINP | 0.3991 | −0.1002 | 0.4992 |
| SON | CINP | 0.3793 | −0.1186 | 0.4979 |
| CKS2 | CINP | 0.3827 | −0.1131 | 0.4958 |
| NASP | CINP | 0.3767 | −0.1174 | 0.4941 |
| RUVBL1 | CINP | 0.3439 | −0.1485 | 0.4924 |
| CINP | SPC25 | 0.3673 | −0.1247 | 0.492 |
| CINP | CCNE1 | 0.3996 | −0.0912 | 0.4908 |
| CINP | SMARCB1 | 0.4548 | −0.0338 | 0.4885 |
| ESCO2 | CINP | 0.3814 | −0.1062 | 0.4876 |
| SPIN1 | CINP | 0.3908 | −0.096 | 0.4869 |
| USP9X | CINP | 0.4047 | −0.0798 | 0.4845 |
| RCC2 | CINP | 0.389 | −0.0952 | 0.4842 |
| UPF1 | CINP | 0.353 | −0.1311 | 0.4841 |
| CINP | DBF4 | 0.3502 | −0.1331 | 0.4833 |
| CINP | CDK1 | 0.4059 | −0.0759 | 0.4818 |
| WTAP | CINP | 0.3834 | −0.0931 | 0.4765 |
| LIG1 | CINP | 0.3865 | −0.0864 | 0.4729 |
| MCM4 | CINP | 0.3894 | −0.0826 | 0.472 |
| CINP | CCAR1 | 0.3655 | −0.1052 | 0.4707 |
| CINP | UHRF1 | 0.3891 | −0.0809 | 0.47 |
| CENPW | CINP | 0.3579 | −0.1119 | 0.4698 |
| SMC2 | CINP | 0.3774 | −0.0921 | 0.4695 |
| CINP | FBXO5 | 0.363 | −0.1049 | 0.468 |
In the G1 versus S comparison ( Supplementary Table S2 ), almost all of the top ρ differences involve the gene CDC45, which is essential for initiation of DNA replication. The top difference is between CDC45 and TIPIN, whose encoded protein is part of the replisome complex. Another top difference is between CDC45 and UHRF1, which plays an important role in the transition from G1 to S. One notable difference is between CDC45 and MCM7, which have been shown to interact with one another. Several top pairs include the ARHGEF2 gene, which is implicated in Rho-GTPase activation.
In the G2M versus S comparison ( Supplementary Table S3 ), most of the top differences involve CINP, whose encoded protein is a component of the DNA replication complex. The top difference is between CINP and NUP37, which is related to the transportation of macromolecules between the cytoplasm and nucleus. CINP and MCM5 is another important difference, where the associated protein for MCM5 is upregulated during transition from G0 to G1 and S cycles. The difference for CINP and PPP1CC was also notable; the protein associated with PPP1CC is involved in the regulation of cell division.
Next, we shift focus to the pediatric-onset MS data. We compare ρ values between MS participants and unaffected controls at the genus level. Supplementary table S7 shows the highest ρ values (in absolute value) in both MS and unaffected controls. The top 25 differences between MS and controls are shown in Table 2 ; a larger table with more genus pairs is shown in Supplementary Table S8 . Genera appearing among the top 25 differences that were previously shown to be associated with MS include Dialister [ 31 ]; Adlercreutzia [ 32 ], Alistipes [ 33 ], Slackia [ 34 ], Barnesiella [ 35 ], Odoribacter [ 36 ], and Collinsella [ 37 ]. Escherichia appears in many of the top genus–genus associations, though there is limited existing evidence of its link with MS other than in a small experiment where Experimental Autoimmune Encephalomyelitis, which is an animal model for MS, was modulated through a colonization of Escherichia coli [ 38 ].
Table 2.
Top 25 genus pairs whose ρ value differs between MS ( ) and control samples ( )
| Genus 1 | Genus 2 | |||
| Gallibacterium | Escherichia | 0.7002 | 0.3253 | 0.3749 |
| Dialister | Phascolarctobacterium | 0.5322 | 0.1975 | 0.3347 |
| Escherichia | Adlercreutzia | 0.5599 | 0.231 | 0.329 |
| Serratia | Escherichia | 0.6953 | 0.3716 | 0.3237 |
| Escherichia | Leuconostoc | 0.676 | 0.3593 | 0.3167 |
| Escherichia | Sphingomonas | 0.6159 | 0.2993 | 0.3166 |
| Escherichia | Vagococcus | 0.7101 | 0.4022 | 0.3079 |
| Klebsiella | Escherichia | 0.7216 | 0.4202 | 0.3013 |
| Escherichia | Rhizobium | 0.6518 | 0.352 | 0.2997 |
| Escherichia | Synechococcus | 0.6531 | 0.3569 | 0.2962 |
| Dialister | Slackia | 0.7149 | 0.4192 | 0.2957 |
| Escherichia | Cupriavidus | 0.6319 | 0.3373 | 0.2946 |
| Segatella | Barnesiella | 0.6308 | 0.3364 | 0.2944 |
| Escherichia | Alistipes | 0.5595 | 0.2653 | 0.2942 |
| Escherichia | Ruminiclostridium | 0.6133 | 0.3196 | 0.2938 |
| Citrobacter | Odoribacter | 0.6841 | 0.3911 | 0.293 |
| Dialister | Collinsella | 0.5686 | 0.278 | 0.2906 |
| Escherichia | Methylobacterium | 0.6292 | 0.3391 | 0.2901 |
| Stenotrophomonas | Escherichia | 0.6012 | 0.3172 | 0.284 |
| Escherichia | Brevibacillus | 0.6462 | 0.3637 | 0.2825 |
| Escherichia | Staphylococcus | 0.6201 | 0.3421 | 0.278 |
| Dialister | Parafannyhessea | 0.7019 | 0.424 | 0.2779 |
| Escherichia | Aristaeella | 0.5216 | 0.2437 | 0.2779 |
| Lacticaseibacillus | Raoultibacter | 0.3834 | 0.6736 | −0.2903 |
| Turicibacter | Raoultibacter | 0.4454 | 0.7833 | −0.338 |
In 23 of the top 25 genus pairs, the estimated ρ -values are greater in MS compared to controls. This trend is present in the top 100 genus pairs as well, with a few exceptions. In Supplementary Fig. S7 , we see that there are only 30 genus pairs common in the top 100 pairs for MS and unaffected controls. This analysis has uncovered a change in community-level dynamics in the MS participants among a number of aforementioned genera whose abundances were already known to be related to MS. The vast majority of the top ρ differences are stronger for MS. This appears to be consistent with previous work in MS. One study showed enrichment of microbial co-occurrences in MS compared to controls from the same household [ 39 ]. Networks in the MS samples showed more connections overall compared to the controls. In a recent meta-analysis, there were several positive network connections observed in MS that were not present in controls [ 40 ].
To investigate the effect of the choice of zero-handling procedures, we reran the MS analysis after imputing zeros using Bayesian-Multiplicative replacement [ 14 ], as implemented in the zCompositions package in R. The top 25 differences between MS and unaffected controls can be seen in Supplementary Table S9 . Of the genus pairs appearing in this new table, 13 were in the original top 25 table. The remaining 12 pairs appear in the top 100 table from the original data analysis. We note that the values in the top 100 differences table have quite a small range (0.2327–0.3749), and so some reordering in the top results is expected after imputation.
Discussion
We have demonstrated that empirical estimates of metrics of proportionality in count-based platforms can lead to bias. Though the bias is mitigated when the mean read depth is large, there may still exist a bias in less-abundant features regardless. We therefore designed a multinomial logit-normal model to calculate model-based estimates of the proportionality metrics. The simulation shows that with high sparsity level and a large number of features, the model-based estimates are much more precise and accurate, as shown by the MSEs. Importantly, empirical versus model-based estimates differed in a way that could drastically change results regarding which features are most strongly correlated with one another.
We performed a timing and memory analysis using data from the simulation study. We ran our model on a candidate dataset from each of n=100,200,500 , and p=50,100,200 . We also added one scenario where n=2000 and p=500 . Results are shown in Supplementary Tables S10 and S11. For the original simulated datasets, the timing and memory usage are quite reasonable, with the n=500 and p=200 scenario requiring 36 min and 2.56 GB of memory. The n=2000 and p=500 scenario took 120 min and required 44.31 GB of memory. The high memory usage in this case was due to running multiple processors simultaneously. Memory usage could be decreased by using fewer processors, albeit at the expense of running time. These timing and memory analyses were run on Digital Research Alliance of Canada’s “Nibi” cluster.
We performed analyses on two datasets using the novel estimation technique: (i) a mESC single-cell RNA dataset, and (ii) a pediatric-onset MS metagenomic dataset. In the mESC dataset, we saw differences in gene–gene associations between the different cell stages that were biologically meaningful; in particular, top differences involved genes that were important in DNA replication. In the MS dataset, we found differences in genus–genus associations involving genera that were previously shown to be associated with MS status. These results show the potential for uncovering new biological insights in count-based datasets when using more appropriate means of calculating metrics of proportionality than the existing empirical calculations.
In both datasets, we have demonstrated that conclusions from our novel proportionality metrics estimation framework can differ from naïvely estimated metrics. In particular, we see attenuation toward the null, on average, due to the shrinkage penalty included in the model. This is an important guard against spurious correlations, given there are a very high number of pairwise measures being estimated in datasets with a large number of features.
We note that, in the data applications in this work, we have applied relatively strict filtering criteria for zeros. There are several reasons for this. First, in the MS dataset, we note that 21% of genera have >90% zeros, some of which were only observed in a single individual. The median zero proportion over all genera is 78%. To put this into perspective, the cmultRepl() zero imputation function in the zCompositions package issues a warning when the proportion of zeros exceeds 80%. Additionally, pairwise measures (such as log-ratio correlations) suffer a high degree of error after imputation on data that have a high proportion of zeros [ 41 ]. In particular, the error increases rapidly after the proportion of zeros is >0.5 in the data.
One limitation of our approach is that the model does not differentiate between structural and sampling zeros, which further underscores the need for filtering out features with many zero counts. We note that the EM algorithm takes longer for convergence when the number of zero proportions is high, and performance is generally worse. Though this could be remedied in the model by introducing zero-inflated parameters as done in [ 20 , 42 ], this introduces a new problem regarding the interpretation of the proportionality metrics, which are not defined for zero-valued arguments. One could conceivably define these metrics conditional on non-zero values, but this would ignore valid information contained in the probability of a zero proportion. To address this, zero-inflated versions of Kendall’s τ and Spearman’s ρ could be used in the context of a zero-inflated model. This is left for future work.
Supplementary Material
Acknowledgements
We are grateful for all the participants’ involvement, especially children and teenagers with MS and their parents. We are also grateful to all the investigators and study teams at each site involved in the Canadian Paediatric Demyelinating Disease Network study, without whom this study would not have been possible. We acknowledge the important contributions of the H.A. team and H.T. team; Thomas Duggan in facilitating study set-up, coordination, and data collection; Bonnie Leung for study coordination; Michael Sargent (Department of Internal Medicine, and the University of Manitoba IBD Clinical and Research Centre laboratory) for managing the stool biobank; and Jessica D. Forbes (University of Toronto) for assisting with the original grant.
Author contributions : Kevin McGregor (Conceptualization [lead], Formal analysis [lead], Methodology [lead], Software [lead], Writing—original draft [lead]), Nneka Okaeme (Formal analysis [supporting], Writing—review & editing [supporting]), Reihane Khorasaniha (Data curation [supporting]), Simona Veniamin (Data curation [supporting]), Juan Jovel (Data curation [supporting]), Richard Miller (Data curation [supporting]), Ramsha Mahmood (Data curation [supporting]), Morag Graham (Data curation [supporting]), Christine Bonner (Data curation [supporting]), Charles N. Bernstein (Data curation [supporting]), Douglas L. Arnold (Data curation [supporting]), Amit Bar-Or (Data curation [supporting]), Ruth Ann Marrie (Data curation [supporting]), Julia O’Mahony (Data curation [supporting]), E. Ann Yeh (Data curation [supporting]), Yinshan Zhao (Data curation [supporting]), Brenda Banwell (Data curation [supporting]), Emmanuelle Waubant (Data curation [supporting]), Gary Van Domselaar (Data curation [supporting]), Feng Zhu (Data curation [supporting]), Ali I. Mirza (Data curation [supporting]), Helen Tremlett (Data curation [lead]), and Heather Armstrong (Data curation [lead])
Contributor Information
Kevin McGregor, Department of Statistics, University of Manitoba, Manitoba, R3T 2M5, Canada.
Nneka Okaeme, Department of Statistics and Actuarial Science, University of Waterloo, Ontario, N2L 3G1, Canada.
Reihane Khorasaniha, Department of Food and Human Nutritional Sciences, University of Manitoba, Manitoba, R3T 2N2, Canada.
Simona Veniamin, Department of Medicine, Division of Gastroenterology, University of Alberta, Alberta, T6G 2B7, Canada.
Juan Jovel, Faculty of Veterinary Medicine, University of Calgary, Alberta, T2N 1N4, Canada.
Richard Miller, Department of Internal Medicine, University of Manitoba, Manitoba, R3A 1R9, Canada.
Ramsha Mahmood, Department of Internal Medicine, University of Manitoba, Manitoba, R3A 1R9, Canada.
Morag Graham, National Microbiology Laboratory, Public Health Agency of Canada, Manitoba, R3E 3M4, Canada.
Christine Bonner, National Microbiology Laboratory, Public Health Agency of Canada, Manitoba, R3E 3M4, Canada.
Charles N Bernstein, Department of Internal Medicine, University of Manitoba, Manitoba, R3A 1R9, Canada.
Douglas L Arnold, Department of Neurology and Neurosurgery, McGill University, Québec, H3A 1A1, Canada.
Amit Bar-Or, Center for NeuroInflammation and Experimental Therapeutics and Children’s Hospital of Philadelphia, University of Pennsylvania, Pennsylvania, 19104, United States.
Ruth Ann Marrie, Departments of Medicine and Community Health and Epidemiology, Dalhousie University, Nova Scotia, B3H 1V7, Canada; MS Clinic at Health Sciences Centre, University of Manitoba, Manitoba R3E 1B2, Canada.
Julia O’Mahony, Mellen Center for Multiple Sclerosis, Cleveland Clinic, Ohio, 44195, United States.
Eluen Ann Yeh, Division of Neurology, Department of Pediatrics, Hospital for Sick Children, University of Toronto, Ontario, M5G 1X8, Canada.
Yinshan Zhao, Department of Medicine, University of British Columbia, British Columbia, V6T 2B5, Canada.
Brenda Banwell, Department of Pediatrics, Johns Hopkins University, Maryland, 21287, United States.
Emmanuelle Waubant, Weill Institute for Neurosciences, University of California San Francisco, California, 94158, United States.
Natalie Knox, National Microbiology Laboratory, Public Health Agency of Canada, Manitoba, R3E 3M4, Canada.
Gary Van Domselaar, National Microbiology Laboratory, Public Health Agency of Canada, Manitoba, R3E 3M4, Canada.
Feng Zhu, Department of Medicine, University of British Columbia, British Columbia, V6T 2B5, Canada.
Ali I Mirza, Department of Medicine, University of British Columbia, British Columbia, V6T 2B5, Canada.
Helen Tremlett, Department of Medicine, University of British Columbia, British Columbia, V6T 2B5, Canada.
Heather Armstrong, Department of Medicine, Division of Gastroenterology, University of Alberta, Alberta, T6G 2B7, Canada; Department of Internal Medicine, University of Manitoba, Manitoba, R3A 1R9, Canada; Centre of Excellence for Gastrointestinal Inflammation and Immunity Research, University of Alberta, Alberta, T6G 2R3, Canada.
Supplementary data
Supplementary data is available at NAR Genomics & Bioinformatics online.
Conflict of interest
H.T. has, in the last 5 years, received research support from the Canada Research Chair Program, the National Multiple Sclerosis Society, the Canadian Institutes of Health Research, the Multiple Sclerosis Society of Canada, the Multiple Sclerosis Scientific Research Foundation, and the EDMUS Foundation (Fondation EDMUS contre la sclérose en plaques). In addition, in the last 5 years, H.T. has had travel expenses or registration fees prepaid or reimbursed to present at CME conferences from the Consortium of MS Centres (2018, 2023), the Canadian Neurological Sciences Federation (2023), National MS Society (2018, 2022), ECTRIMS/ACTRIMS (2017–2023), and the American Academy of Neurology (2019). Speaker honoraria are either declined or donated to an MS charity or to an unrestricted grant for use by H.T.’s research group.
Funding
K.M. acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC Discovery Grant: RGPIN-2021-03634) as well as an NSERC Undergraduate Student Research Award for N.O. The MS study was supported by funding from Canada Research Chair (PI: H.A.), the Multiple Sclerosis Society (PI: H.A.), and the Multiple Sclerosis Scientific and Research Foundation (#EGID: 2636; PI: H.T.). These latter funding sources were not involved in the study design, the collection, analysis, and interpretation of the data, or in the decision to submit this article for publication.
Data availability
The mouse embryonic stem cell dataset is publicly available in BioStudies ( https://www.ebi.ac.uk/biostudies/ ) under accession number E-MTAB-2512. The pediatric-onset MS dataset was submitted to the SRA archive of NCBI ( https://www.ncbi.nlm.nih.gov/sra ) and is publicly available under accession number PRJNA1022308. An R package implementation is available at https://cran.r-project.org/package=countprop . For the analyses in this paper, data are available at https://doi.org/10.5281/zenodo.17782904 . Code is available at https://github.com/kevinmcgregor/proportionality_code .
References
1. Aitchison J. The Statistical Analysis of Compositional Data. Chapman and Hall, 1986.
2. Gloor GB, Macklaim JM, Pawlowsky-Glahn V et al. Microbiome datasets are compositional: and this is not optional. Front Microbiol. 2017;8:2224. 10.3389/fmicb.2017.02224
3. McMurdie PJ, Holmes S. Waste not, want not: why rarefying microbiome data is inadmissible. PLoS Comput Biol. 2014;10:e1003531. 10.1371/journal.pcbi.1003531
4. Egozcue JJ, Pawlowsky-Glahn V. Compositional data: the sample space and its structure. TEST. 2019;28:599–638. 10.1007/s11749-019-00670-6
5. Lovell D, Pawlowsky-Glahn V, Egozcue JJ et al. Proportionality: a valid alternative to correlation for relative data. PLoS Comput Biol. 2015;11:e1004075. 10.1371/journal.pcbi.1004075
6. Erb I, Notredame C. How should we measure proportionality on relative gene expression data?. Theory Biosci. 2016;135:21–36. 10.1007/s12064-015-0220-8
7. Skinnider MA, Squair JW, Foster LJ. Evaluating measures of association for single-cell transcriptomics. Nat Methods. 2019;16:381–6. 10.1038/s41592-019-0372-4
8. Egozcue JJ, Pawlowsky-Glahn V. Subcompositional coherence and a novel proportionality index of parts. SORT. 2023;47:229–44.
9. Jin S, Notredame C, Erb I. Compositional covariance shrinkage and regularised partial correlations. SORT. 2023;47:245–68. 10.57645/20.8080.02.8
10. Ma J, Xie P, Pantoja K et al. Variable selection in balance regression with applications to microbiome compositional data. arXiv, https://arxiv.org/abs/2304.00143 , 31 March 2023, preprint: not peer reviewed.
11. Lovell DR, Chua XY, McGrath A. Counts: an outstanding challenge for log-ratio analysis of compositional data in the molecular biosciences. NAR Genom Bioinform. 2020;2:lqaa040. 10.1093/nargab/lqaa040
12. Lun A. Overcoming systematic errors caused by log-transformation of normalized single-cell RNA sequencing data. bioRxiv, 10.1101/404962, 31 August 2018, preprint: not peer reviewed.
13. Egozcue JJ, Graffelman J, Ortego MI et al. Some thoughts on counts in sequencing studies. NAR Genom Bioinform. 2020;2:lqaa094.
14. Palarea-Albaladejo J, Martín-Fernández JA. zCompositions—R package for multivariate imputation of left-censored data under a compositional approach. Chemomet Intell Lab Syst. 2015;143:85–96.
15. Xia F, Chen J, Fung WK et al. A logistic normal multinomial regression model for microbiome compositional data analysis. Biometrics. 2013;69:1053–63.
16. Tian C. Microbial network recovery by compositional graphical lasso under additive log-ratio transformation. PhD Thesis, Oregon State University, 2020.
17. Willis AD, Martin BD. Divnet: estimating diversity in networked communities. Biostatistics. 2020;21:586–603.
18. Silverman JD, Roche K, Holmes ZC et al. Bayesian multinomial logistic normal models through marginally latent matrix-T processes. J Mach Learn Res. 2022;23:7–1.
19. McGregor K, Labbe A, Greenwood CMT. MDiNE: a model to estimate differential co-occurrence networks in microbiome studies. Bioinformatics. 2020;36:1840–7.
20. Ma S, Ren B, Mallick H et al. A statistical model for describing and simulating microbial community profiles. PLOS Comput Biol. 2021;17:e1008913.
21. Friedman J, Hastie T, Tibshirani R. Sparse inverse covariance estimation with the graphical lasso. Biostatistics. 2008;9:432–41. 10.1093/biostatistics/kxm045
22. Hoff PD. Nonparametric modeling of hierarchically exchangeable data. Technical Report 421, University of Washington, Statistics Department, 2003.
23. Chen J, Chen Z. Extended BIC for small-n-large-p sparse GLM. Stat Sinica. 2012;22:555–74.
24. Buettner F, Natarajan KN, Casale FP et al. Computational analysis of cell-to-cell heterogeneity in single-cell RNA-sequencing data reveals hidden subpopulations of cells. Nat Biotechnol. 2015;33:155–60. 10.1038/nbt.3102
25. Kurtz ZD, Müller CL, Miraldi ER et al. Sparse and compositionally robust inference of microbial ecological networks. PLoS Comput Biol. 2015;11:e1004226. 10.1371/journal.pcbi.1004226
26. Tremlett H, Zhu F, Arnold D et al. The gut microbiota in pediatric multiple sclerosis and demyelinating syndromes. Ann Clin Transl Neur. 2021;8:2252–69. 10.1002/acn3.51476
27. Mirza AI, Zhu F, Knox N et al. Metagenomic analysis of the pediatric-onset multiple sclerosis gut microbiome. Neurology. 2022;98:e1050–63. 10.1212/WNL.0000000000013245
28. Mirza AI, Zhu F, Knox N et al. The metabolic potential of the paediatric-onset multiple sclerosis gut microbiome. Mult Scler Relat Dis. 2022;63:103829. 10.1016/j.msard.2022.103829
29. Stelzer G, Rosen N, Plaschkes I et al. The GeneCards suite: from gene data mining to disease genome sequence analyses. Curr Protoc Bioinform. 2016;54:1–30. 10.1002/cpbi.5
30. Li X, Huang W, Huang W et al. Kinesin family members KIF2C/4A/10/11/14/18B/20A/23 predict poor prognosis and promote cell proliferation in hepatocellular carcinoma. Am J Trans Res. 2020;12:1614.
31. Bonnechère B, Amin N, van Duijn C. What are the key gut microbiota involved in neurological diseases? A systematic review. Int J Mol Sci. 2022;23:13665. 10.3390/ijms232213665
32. Schepici G, Silvestro S, Bramanti P et al. The gut microbiota in multiple sclerosis: an overview of clinical trials. Cell Transplant. 2019;28:1507–27. 10.1177/0963689719873890
33. Schwerdtfeger LA, Montini F, Lanser TB et al. Gut microbiota and metabolites are linked to disease progression in multiple sclerosis. Cell Rep Med. 2025;6:102055.
34. Duarte-Silva E, Meuth SG, Peixoto CA. Microbial metabolites in multiple sclerosis: implications for pathogenesis and treatment. Front Neurosci. 2022;16:885031. 10.3389/fnins.2022.885031
35. Nitzan Z, Staun-Ram E, Volkowich A et al. Multiple sclerosis-associated gut microbiome in the Israeli diverse populations. Int J Mol Sci. 2023;24:15024. 10.3390/ijms241915024
36. Horton MK, McCauley K, Fadrosh D et al. Gut microbiome is associated with multiple sclerosis activity in children. Ann Clin Transl Neur. 2021;8:1867–83. 10.1002/acn3.51441
37. Jangi S, Gandhi R, Cox LM et al. Alterations of the human gut microbiome in multiple sclerosis. Nat Commun. 2016;7:12015. 10.1038/ncomms12015
38. Libbey JE, Sánchez JMS, Fleming BA et al. Modulation of experimental autoimmune encephalomyelitis through colonisation of the gut with Escherichia coli . Benef Microbes. 2020;11:669–84. 10.3920/BM2020.0083
39. Zhou X, Baumann R, Gao X et al. Gut microbiome of multiple sclerosis patients and paired household healthy controls reveal associations with disease risk and course. Cell. 2022;185:3467–86. 10.1016/j.cell.2022.08.021
40. Lin Q, Dorsett Y, Mirza A et al. Meta-analysis identifies common gut microbiota associated with multiple sclerosis. Genome Med. 2024;16:94. 10.1186/s13073-024-01364-x
41. Lubbe S, Filzmoser P, Templ M. Comparison of zero replacement strategies for compositional data with large numbers of zeros. Chemomet Intell Lab Syst. 2021;210:104248. 10.1016/j.chemolab.2021.104248
42. Zeng Y, Pang D, Zhao H et al. A zero-inflated logistic normal multinomial model for extracting microbial compositions. J Am Stat Assoc. 2022;118:2356–69.
Associated Data
Supplementary Materials
Data Availability Statement
The mouse embryonic stem cell dataset is publicly available in BioStudies ( https://www.ebi.ac.uk/biostudies/ ) under accession number E-MTAB-2512. The pediatric-onset MS dataset was submitted to the SRA archive of NCBI ( https://www.ncbi.nlm.nih.gov/sra ) and is publicly available under accession number PRJNA1022308. An R package implementation is available at https://cran.r-project.org/package=countprop . For the analyses in this paper, data are available at https://doi.org/10.5281/zenodo.17782904 . Code is available at https://github.com/kevinmcgregor/proportionality_code .