A heavy-tailed model for analyzing miRNA-seq raw read counts

2.1 Discrete stable distributions

Discrete stable distributions form a class of probability distributions defined on nonnegative integers and allowing heavy tails. These features make discrete stable distributions attractive for fitting miRNA-seq raw read counts. Discrete stable distributions are characterized by two parameters, a tail index parameter α ∈ (0, 1] and a parameter λ > 0, in this paper, denoted by DStable(α, λ). The Poisson distribution forms a special case of discrete stable distributions with α = 1 and corresponds to the lightest tail of the discrete stable distributions. Discrete stable distributions with α < 1 exhibit heavy tails (e.g., Soltani et al. 2009), also referred to as Paretian tails (e.g., Rémillard and Theodorescu 2000, Remark 2.7). The lower the tail index α ∈ (0, 1), the slower the decay, the heavier the tail and the more prone to extreme values is the distribution.

2.1.2 Probability calculations

Except for the Poisson distribution, probability mass functions of discrete stable distributions cannot be explicitly expressed in terms of elementary functions. In particular, when calculating Pr _θ (X = k) for a discrete stable distribution with parameters θ = (α, λ) and X ∼ DStable(α, λ), no explicit analytic form is available. In (Christoph and Schreiber 1998, Formula (8), p. 245) a recursion formula was derived (see also Steutel and van Harn 2003, Formula (5.18)):

(2) Pr θ ( X = k + 1 ) = λ ∑ m = 0 k ( − 1 ) m m + 1 k + 1 Pr θ ( X = k − m ) α m + 1 ,

where k = 1, 2, 3, … and Pr _θ (X = 0) = e^−λ . Explicit formulas follow: Pr _θ (X = 1) = αλe^−λ , Pr θ ( X = 2 ) = α λ 2 e − λ ( α ( λ − 1 ) + 1 ) and Pr θ ( X = 3 ) = α λ 3 e − λ [ α λ 2 ( α ( λ − 1 ) + 1 ) − α ( α − 1 ) λ + 1 2 ( α − 1 ) ( α − 2 ) ] .

The impact of changing λ for a fixed tail index α is illustrated in Figure 1. Here, the higher values of λ not only shift the distribution rightward, but also disperse the primary distribution mass (Table 1).

Figure 1:

For tail indices α = 0.1, 0.5, 0.9 and parameters λ = 2, 4, 5, the values of Pr _θ (X = k) of X ∼ DStable(α, λ) are calculated using the recursive formula (2).

Table 1:

For tail indices α = 1, 0.9, 0.8, 0.4, 0.2 and parameter λ = 1, the values of Pr _θ (X = k) of X ∼ DStable(α, λ). Calculations up to k ≤ 10⁴ are made using the recursive formula (2), and for values k ≥ 10⁴, the asymptotic tail formula (3) is used.

k	α = 1 ≡ Pois(1)	α = 0.9	α = 0.8	α = 0.4	α = 0.2
0	0.368	0.368	0.368	0.368	0.368
10	1.01 × 10−7	1.71 × 10−3	3.76 × 10−3	9.01 × 10−3	6.54 × 10−3
102	3.94 × 10−159	1.55 × 10−5	4.55 × 10−5	2.82 × 10−4	4.98 × 10−4
103	0	1.89 × 10−7	6.97 × 10−7	8.92 × 10−6	3.47 × 10−5
104	0	2.38 × 10−9	1.10 × 10−8	2.82 × 10−7	2.39 × 10−6
105	0	2.99 × 10−11	1.74 × 10−10	8.92 × 10−9	1.58 × 10−7
106	0	3.77 × 10−13	2.76 × 10−12	2.82 × 10−10	1.03 × 10−8

The recursive formula (2) is straightforward, but the calculation time grows rapidly as k increases, and it is computationally time consuming, especially when handling extremes. For calculating the probabilities of large values, (Christoph and Schreiber 1998, Formulas (11) and (12)) proved a tail asymptotic formula for k → ∞:

(3) Pr θ ( X = k ) = 1 π ∑ j = 1 [ ( α + 1 ) / α ] ( − 1 ) j + 1 j ! λ j ⁡ sin ( α j π ) Γ ( α j + 1 ) k − α j − 1 + O ( k − α − 2 ) .

Note that (Doray et al. 2009, pp. 2008) derived an explicit formula for calculating Pr _θ (X = k), which for k = 1, 2, 3, … can be expressed via Gamma functions as follows

(4) Pr θ ( X = k ) = ( − 1 ) k e − λ ∑ m = 0 k ∑ j = 0 m m j α j k ( − 1 ) j λ m m !

(5) = ( − 1 ) k α e − λ k ! ∑ m = 1 k λ m ∑ j = 1 m ( − 1 ) j Γ ( α j ) ( j − 1 ) ! ( m − j ) ! 1 Γ ( α j + 1 − k ) ,

where α ∈ (0, 1], λ > 0. For computations in R, the reciprocal Gamma function 1/Γ(αj − k + 1) can be calculated as follows (Prodanov 2019) 1 / Γ ( − z ) = − sin ⁡ π z π Γ ( z + 1 ) . However, for further calculations, in R, there is a limitation when using ‘factorial’ and ‘gamma’ functions for larger values (see R Core Team 2020).

Hence, for further analysis, we adopt a hybrid formula: combining the recursive formula (2) with the tail asymptotic formula (3), with a change point for values bigger than pre-defined K > 0 between the formulas. In this paper we set the change point K = 1000.

2.2 Estimating the parameters of discrete stable distributions

Modeling via discrete stable distribution means estimating its parameters α and λ. Because the probability mass function does not have a closed analytical form and not all moments exist, the popular maximum likelihood or moment-based techniques are not directly applicable for discrete stable distributions. A popular approach is to use the empirical probability generating function (EPGF): let X = {X ₁, …, X _n } be a random sample from a probability distribution, then the empirical probability generating function G _n :[−1, 1] → [0, 1], is defined as

To illustrate the behavior of the empirical probability generating function G _n (z) compared to the probability generating function of discrete stable laws G(z), as given by Equation (1), individual samples of discrete stable distributions are generated by simulation. For each sample, the absolute differences in the values of the theoretical and empirical functions, | G ̂ n ( z ) − G ( z ) | for z ∈ [−1, 1], are presented in Figure 2.

Figure 2:

The absolute error between the empirical and theoretical probability generating function, | G ̂ n ( z ) − G ( z ) | , z ∈ [ − 1,1 ] , where G ̂ n is based on a sample from DStable(α, 1) with sizes n = 10, 100, 1000.

By the strong law of large numbers, G _n (z) → G(z) almost surely as n → ∞ for each z ∈ [−1, 1]. Estimation based on empirical probability generating functions has been suggested in Doray et al. (2009) and Slámová and Klebanov (2014). However, the former has shortcomings in terms of search procedures and the latter requires the estimation of additional parameters. In this paper, computationally straightforward closed-form estimators are used, proposed in Marcheselli et al. (2008). Fix two points z ₁, z ₂ ∈ (−1, 1), z ₁ ≠ z ₂ and solve the corresponding system of equations of (1) for α and λ. Substituting G(z) with G _n (z) leads to point estimators for α and λ:

The preference for estimators (7) is rooted in their advantageous characteristics: they offer closed-form simplicity and computational efficiency. Generally, these estimators (7) are versatile and applicable for evaluation in any z ₁, z ₂ ∈ (−1, 1) where z ₁ ≠ z ₂. Although the approach is attractive from a computational perspective, determining the choice of z ₁, z ₂ to achieve the best estimates remains an open question, analogous to related studies in continuous stable distributions (Krutto 2018; Lember and Krutto 2022). However, note that G n ( − 1 ) = 1 n ∑ i = 1 n ( − 1 ) X i , G n ( 0 ) = 1 n ∑ i = 1 n I { X i = 0 } , and G _n (1) = 1 while referring to Equation (1), for discrete stable laws the probability generating function gives: G ( − 1 ) = e − λ 2 α , G(0) = e^−λ , and G(1) = 1. Consequently, at z = 1, we have G _n (1) = G(1) = 1. This can be interpreted as the values of G _n (z) and G(z) being closer as z approaches 1 (compared to other values of z in [−1,1]), a trend also observed in Figure 2. Hence, it is suggested (e.g., Marcheselli et al. 2008) to consider arguments z ₁ ≠ z ₂ close to 1.

To identify the optimal choices of z ₁, z ₂, maximum likelihood calculations can be used. Consider a realization k = (k ₁, …, k _n ) of a discrete random vector (e.g., miRNA-seq read counts) X = (X ₁, …, X _n ), where X _i represents the random variable associated with observation i. In this context, the vector k represents the observed values, while X represents the corresponding random variables.^[4] The likelihood function L( θ ) is then given by L ( θ ∣ k ) = ∏ j = 1 n p X ( k j ∣ θ ) , where p _X ( k ∣ θ ) = Pr _θ (X = k) and θ = (α, λ) represent the discrete stable distribution parameters. The maximum likelihood estimate θ ̂ n * z 1 * , z 2 * , k is the maximizer of L ( θ ̂ n ( z 1 , z 2 ) ∣ k ) , where θ ̂ n ( z 1 , z 2 , k ) = ( α ̂ n ( z 1 , z 2 , k ) , λ ̂ n ( z 1 , z 2 , k ) ) is computed via formulas (7) at z ₁, z ₂ from a (sub)grid of [0, 1) × [−1, 1). Calculating Pr _θ (X = k), where X ∼ DStable(α, λ), α ∈ (0, 1], and λ > 0, involves combining the recursive formula (2) with the tail asymptotic formula (3). For details, see Algorithm 1.

2.3 The goodness-of-fit of the Poisson, the negative binomial and the discrete stable distributions

The goodness-of-fit of the Poisson, the negative binomial and the discrete stable distributions is compared using Akaike Information Criteria (AIC). Recall that AIC = −2l( θ ∣k) + 2s, where s is the number of parameters and l( θ ∣ k ) is the log-likelihood function. Given k = (k ₁, …, k _n ), the realization of a discrete random vector X = (X ₁, …, X _n ), the log-likelihood function is given by l ( θ ∣ k ) = ∑ j = 1 n ⁡ log p X ( k j ∣ θ ) , where p(k _j ∣ θ ) = Pr _θ (X _j = k _j ) and θ is the vector of distribution parameters. For the discrete stable distribution, the values of Pr _θ (X _j = k _j ) are calculated using the formulas (2) and (3). For the Poisson and negative binomial distributions, the log-likelihood functions are computed using the R package stats. When fitting the data, both the Poisson and the negative binomial distributions are fitted using the Moment Matching Estimation (MME), available in the R package fitdistrplus (Delignette-Muller and Dutang 2015). The moment method is chosen because of computational challenges encountered when running the Maximum Likelihood Estimation (MLE) method for the negative binomial distribution. However, note that for the Poisson distribution, the MLE and moment estimators yield the same results.

To further assess goodness-of-fit, we present QQ plots and estimate Cramer von Mises statistic ω ² and Anderson–Darling statistics A (e.g., Stephens et al. 1986, pp. 100) for a selection of miRNAs. These are based on calculations performed on the cumulative distribution functions, for which the discrete stable distribution lacks a closed analytical form. Nevertheless, we can estimate the cumulative distribution function for the discrete stable distribution based on the probability mass function, calculated using the hybrid formula described in Section 1.1. We mention that for the Poisson and the negative binomial distributions, advanced forms of Cramér-von Mises and Anderson–Darling tests are implemented in the R package goftest (Faraway et al. 2021). However, to have comparable results for all, we calculate same basic statistics, which can be used as a measure of distance to assess fit. For the Cramér-von Mises w ², we use:

and for the Anderson–Darling A, we use:

where X ₍₁₎, …, X _(n) represent ordered data, and F is the estimated distribution function of the Poisson, negative binomial and discrete stable, respectively.

Caution is advised when interpreting the QQ plots and the Cramer-von Mises and Anderson–Darling statistics for the discrete stable distribution due to the lack of a closed-form solution.

2.4 Using discrete stable distributions for differential expression analysis

Differential expression analysis for miRNAs aims to identify miRNAs whose expression levels change notably under different experimental conditions, such as various cancer subtypes, healthy versus diseased states, or pre- and post-treatment stages. The objective is to understand the roles of these molecular regulators in biological processes, disease mechanisms, and as potential therapeutic targets. This analysis involves comparing miRNA expression profiles to discover specific miRNAs that might be crucial for driving biological changes. These changes can involve overexpression (increased expression levels) or underexpression (reduced expression levels), providing critical information on disease progression and treatment responses. The primary challenge in differential expression analysis is to distinguish genuine biological expression changes from the inherent variability of the underlying nature of miRNA expression. Extremely large values often observed in miRNA expression have traditionally been considered experimental noise, leading to their elimination through normalization and preprocessing procedures. However, this heavy-tailed nature might represent the underlying distribution rather than mere experimental noise, as discussed in the Introduction section. Therefore, the development of alternative methods for differential expression analysis directly handling raw data becomes crucial. At the outset, we highlight the availability of popular robust non-parametric tests, recommended for data containing extreme values (e.g., Staudte and Sheather 2011; Wilcox 2022). For example, the rank-based Brunner-Munzel test (e.g., Kume et al. 2017) for the stochastic equality of two populations. This test accommodates ties and does not assume equal variances between the groups. However, if the null hypothesis is rejected, it implies differences between distributions, but it provides limited insight into the nature of these differences. Additionally, the Moody’s median test can be applied to test for differences in medians. Although the median is robust against heavy-tailed distributions, a limitation arises as our interest extends beyond merely determining median based central tendency. That is, the median serves as a singular metric, and no significant differences in medians does not necessarily imply that heavy-tailed distributions are identical or vice versa. However, we recommend using both tests as additional tools in the context of heavy-tailed differential expression analysis. In R, the Brunner-Munzel test can be applied by using the package brunnermunzel (Ara 2020) and the median test via the package agricolae (Felipe de Mendiburu 2023). In this paper, we propose that examining differences in the parameters α and λ of the discrete stable distribution could be more sensitive and suitable for discovering group differences, considering that miRNAs are derived from heavy-tailed distributions. For each miRNA, let λ ̂ n > 0 and α ̂ n ∈ ( 0,1 ] be the estimates of the parameters of the discrete stable distribution, fitted to the raw expressions of the miRNA sequence belonging to a group of interest. We denote the (true) parameters of some group A by λ _A , α _A and corresponding estimates by ( λ ̂ n ) A , ( α ̂ n ) A . To measure the difference between the parameter estimates, for each miRNA, calculate (1) the ratios of the parameter estimates of λ of the reference A to some comparison group B: ( λ ̂ n ) A / ( λ ̂ n ) B and (2) the differences of the corresponding tail index estimates: ( α ̂ n ) A − ( α ̂ n ) B . (we do not use absolute value, as the sign indicates the direction of the tail heaviness. However, for illustrative purposes, the absolute difference may be used in some situations). This approach allows for a visual representation of the differences in estimates in high dimensions: assuming that the x-axis represents ( α ̂ n ) A − ( α ̂ n ) B and the y-axis represents ( λ ̂ n ) A / ( λ ̂ n ) B . All miRNAs with minimal differences are located around zero. MiRNAs with similar tail index estimates but differences in parameter lambda estimates align around a vertical line drawn from zero, parallel to the y-axis. In general, differences in both α and λ are of interest. Differences in the tail index α estimates indicate a possible difference in population tails, which affects the proneness to extreme values. However, in a differential expression analysis context, the parameter λ might be more interesting. Since λ impacts both the center and concentration of the primary mass of the distribution in discrete stable distributions (see the discussion of Figure 1 in Section 1.1), it might provide a more comprehensive understanding of differences between distributions across groups. However, as we see from formulas (7), differences in the tail index estimates α complicate the direct interpretation of the related difference (or lack thereof) in the estimates of λ. On the other hand, minor variations in the tail index estimates α indicate comparable tail characteristics, potentially facilitating a more straightforward comparison of the λ estimates. In practical terms, a difference of up to 0.1 could be considered small, indicating similar tail index. Importantly, as differences are identified from parameter estimates, the incorporation of a formal test proves to be beneficial. For similar tail indexes, differences in λ estimates might signal differences in typical patient expressions: higher λ estimates indicate a tendency toward elevated miRNA expressions by shifting and/or spreading the distribution’s mass over a wider range.

Regarding the significance of differences in the parameter λ, a permutation test (e.g., Wilcox 2022) based on a statistic such as:

can be used to assess the significance of differences in λ between groups A and B by comparing the statistic computed from the original sample with those obtained from permuted samples. Permutation tests involve multiple random reassignments of observations, simulating the null hypothesis distribution to evaluate the statistical significance of observed differences. However, testing all possible permutations can be computationally intensive. Hence, in practice, random sampling from permutations is typically adopted. The corresponding p-value based on i = 1, …, n _perm values D(i) is calculated as:

where D(i) represents the difference based on the ith sample permutation, D _o is the observed value of the original data, and n _perm is the number of permutations.

3.1 Exploratory data analysis for heavy tails

The summary of some classic and robust statistics for the TCGA-BRCA data is given in Table 2 and for the NOWAC-LUCA data in Table 3.

Both Tables 2 and 3 show misalignment between the classic and robust central tendency and dispersion measures, and both tables suggest that the underlying probability distribution might be skewed to the right and heavy tailed. We point out that even if the theoretical mean does not exist due to the heavy-tailed nature of the distribution, the sample average (sample mean), as presented in Tables 2 and 3, can still be calculated from the observed data. However, in this case, the sample mean will not serve as an estimate of the non-existent theoretical mean. With regard to the popular assumption that miRNAs follow a Poisson distribution, we also draw attention to the inconsistency with this distribution. Theoretically, for the Poisson distribution, the variance-to-mean ratio, the product of mean and skewness squared, and the product of mean and excess kurtosis are all equal to 1. A summary of these measures for the TCGA-BRCA data is given in Table 2 and, correspondingly, for the NOWAC-LUCA data in Table 3. From both Tables 2 and 3, it follows that the miRNA-seq count data are overdispersed, overskewed, and have heavier tails, relative to the Poisson distribution. However, note that sample skewness and excess kurtosis can be severely biased in finite samples (e.g., Joanes and Gill 1998), so these results should be treated with care.

The importance of exploratory data analysis to understand heavy-tailed distributions and the properties of underlying discrete probability distributions is highly emphasized; However, the identification of heavy tails often lacks formal criteria and relies instead on visual inspection (e.g., Embrechts et al. 2013, Chapter 6). Common methods in exploratory data analysis for heavy tails, such as probability and QQ plots, the mean excess function plots, and Gumbel’s method of exceedance, tend to be better suited for a single-variable or few-variables analysis (e.g., Embrechts et al. 2013, Section 6.2), and might become impractical in high-dimensional settings. A simple tool for exploring heavy tails and the finiteness of moments is the ratio of maximum to sum (e.g., Embrechts et al. 2013, Section 6.2.6). Suppose X ₁, …, X _n are iid random variables and for any positive r define the quantity R _n (r) = M _n (r)/S _n (r), where S n ( r ) = ∣ X 1 ∣ r + ⋯ + ∣ X n ∣ r , M n ( r ) = max ∣ X 1 ∣ r , … , ∣ X n ∣ r , r > 0, n ≥ 1. The following equivalence holds: R _n (r) → ^a.s. 0 if and only if E∣X∣ ^r < ∞. Therefore, R _n (r) should be small for large n provided that E∣X∣ ^r < ∞, otherwise it indicates heavy tails. The advantage of using R _n (r) is its ability to visually combine many variables into a plot, which is especially useful in high-dimensional data analysis. For r = 1, 2, 3, 4, the ratios of maximum to sum R _n (r) for the TCGA-BRCA data are presented in Figure A.1-F.1 and for the NOWAC-LUCA data in Figure A.3-F.1. From both Figures A.1-F.1 and A.3-F.1, the results suggest that most of the analyzed miRNAs may not have moments of order r = 1, 2, 3, 4. In summary, exploratory data analysis for heavy tails of miRNA expressions indicates that the underlying distributions are likely to be heavy tailed.

4.1 TCGA-BRCA: estimates of discrete stable distribution parameters

The discrete stable distributions were fitted to p = 400 miRNAs, both for data from all samples (all, n = 1076) and for data from the molecular subtypes, normal-like (Normal, n = 40), luminal A (LumA, n = 568), luminal B (LumB, n = 202), basal-like (Basal, n = 185), HER2-amplified (positive) (Her2, n = 81). A summary of the corresponding parameter estimates α ̂ n and λ ̂ n is given in Table 4. Table 4 displays estimates derived from various samples. ‘All samples’ encompass all subtypes together, while samples from specific subtypes are selected as subsets. The parameters are then estimated separately for each of these subsets.

From Table 4, the tail index α of each of the p = 400 miRNAs was estimated to be less than 1, which means that the underlying distribution was estimated to be heavy tailed. In particular, for data from all samples, half of the miRNAs are estimated to have a tail index less than 0.65. Focusing on the subgroups, also here all tail index estimates was less than 1. However, the estimates are slightly larger (i.e., slightly lighter tails) than the ones for all data. Estimates of the parameter λ range from almost zero to hundreds of thousands, both for samples from all data and from molecular subtypes. When comparing the subtypes, the λ estimates are similar to those for all samples, except for the mean and maximum. Indeed, the maximum of the λ estimates is affected by assigning data entries to subgroups, which accordingly affects the mean. However, in all groups, the median of λ ̂ n is around 10, while at least a quarter of λ estimates are more than 50. More analysis of the differences between the groups, at the miRNA level, is given in Sections 3.3 and 3.4.

4.2 TCGA-BRCA: comparison of the goodness-of-fit

The goodness-of-fit of the Poisson, negative binomial, and discrete stable distributions are compared using AIC. For comparison, for every miRNA, the AIC ratio of the estimated discrete stable to the estimated Poisson distribution was evaluated, and the corresponding ratio of the estimated discrete stable to the estimated negative binomial was evaluated. To summarize these results, box plots of the ratios are presented in Figure 3.

Figure 3:

TCGA-BRCA: for p = 400 miRNAs, the ratios of AICs of discrete stable models to Poisson and to negative binomial for data from all samples (n = 1076) and from samples of molecular subtypes, normal-like (n = 40), LumA (n = 568), LumB (n = 202), Basal (n = 185), Her2 (n = 81).

The ratios are given both for the models fitted to the data from all samples together (labelled ‘all’) and for data from the molecular subtypes (labelled correspondingly). If the ratio is less than 1, then the discrete stable distribution gives a better fit. The smaller the ratio, the better the fit of the discrete stable distribution (in comparison to the Poisson and negative binomial distribution). From the left side of Figure 3, it can be seen that for all miRNAs, the discrete stable distribution gives a much better fit than the Poisson, both for the data from all samples together and for the data from the molecular subgroups. In all samples, for around 200 miRNAs the AIC of the discrete stable distribution is approximately 90 percent reduced, and for around 300 miRNAs the values are at least 50 percent reduced. For no miRNA, the Poisson distribution was a better fit than the discrete stable distribution. From the right side of Figure 3, it follows that for all samples together, the discrete stable distribution gives a better fit than the negative binomial for around 344 (86 %) of the miRNAs. For the remaining 56 (14 %) miRNAs, the ratios are close to 1, meaning that the goodness-of-fit is similar. For data from the LumA, LumB, Basal, Her2 and normal-like subtypes, the discrete stable distribution gives a better fit than negative binomial for 58 %, 69 %, 63 %, 73 % and 58 % of the miRNAs, respectively. For a more thorough investigation of goodness-of-fit for a selected set of miRNAs, see Section 3.4.

4.3 TCGA-BRCA: differential expression via discrete stable distributions

In this section, differential expression analysis is performed through estimates of the parameters of discrete stable distributions within the molecular subtypes of breast cancer. In particular, the estimates for the subtypes LumA (n = 568), LumB (n = 202), HeR2 (n = 81) and Basal (n = 186) are compared with the ones for the normal-like subgroup (n = 40). The general summary of the summed expression for all reads aligned per 400 miRNA (measured in primary breast cancer tissue) is given in Table 5.

Certain tendencies can be observed from Table 5. For example, the minimum expression in each group is 0, and the quartiles are also broadly similar. However, for all quartiles, the values of the other subtypes are slightly lower (i.e., underexpressed) compared to the normal-like subgroup, which is in agreement with previous findings (e.g., Fontana et al. 2021, and references therein.) Extremely large maximum values are present in all groups, but no obvious differences appear.

The general summary of the parameter estimates of the fitted discrete stable distributions of samples of molecular subtypes is given in Table 4. Next, for each miRNA the differential expression analysis was performed between normal-like and the other molecular subtypes, as explained in Section 1.4.

For each miRNA, we calculated (1) the ratios of estimates for the parameter λ of the normal-like subgroup and the other subgroups and (2) the differences in estimates for the tail index λ of the normal-like subgroup and the other subgroups. For each miRNA, both calculations are plotted in Figure 4.

Figure 4:

TCGA-BRCA: for p = 400 miRNAs, the ratios of estimates for parameter λ and the differences in estimates for the tail index α, between the normal-like subgroup and the other breast cancer subgroups, labeled correspondingly.

Note that in Figure 4, the ratios of the estimates of λ larger than 1 represent underexpressions, while the ratios less than 1 represent overexpressions, compared to the normal-like subgroup.

In Figure 4, compared to the normal-like subgroup, the five miRNAs that are most underexpressed (irrespective of differences in the tail index α estimates) in each subtype are as follows:

1. LumA: hsa-mir-21, hsa-mir-10a, hsa-mir-199b, hsa-mir-199a-2, hsa-mir-182;

2. LumB: hsa-mir-143, hsa-mir-21, hsa-mir-199b, hsa-mir-93, hsa-mir-199a-2;

3. Basal: hsa-mir-199a-2, hsa-mir-199b, hsa-mir-199a-1, hsa-mir-126, hsa-mir-100;

4. Her2: hsa-mir-199b, hsa-mir-183, hsa-mir-199a-2, hsa-mir-151a, hsa-mir-182.

Although the aforementioned miRNAs have the highest ratios of λ estimates, they also have large differences in the tail index estimates, making comparison of the estimates of λ less straightforward. As such, we can conclude that for these miRNAs, the discrete stable distribution parameters were estimated to be different. However, as explained in Section 1.4, in Figure 4 we are even more interested in the cases where the difference in estimates for the tail index α is relatively small, while the corresponding ratio of estimates for the parameter λ is different from one. In that case, the interpretation of λ differences is more straightforward. For that purpose, a zoom into Figure 4, where only the part with difference in tail index α estimates between the groups of interest is less than 0.1, is presented in Figure 5. In Figure 5, the ratios within the LumA subtype appear to fluctuate less around the value of 1 than within other subtypes. This agrees with previous studies in which miRNA expressions of the LumA subtype have been found to be similar to the normal-like subgroup (e.g., Fontana et al. 2021).

Figure 5:

TCGA-BRCA: a zoom to Figure 4, where the absolute differences in estimates for tail index α is less than 0.1.

The ratios of λ estimates that stand out are all larger than 1, representing underexpressions compared to the normal-like subgroup. The top 10 ratios (larger than 1) from Figure 5, with the corresponding estimated differences in tail index, are given in Table 6.

In particular, in Table 6 the miRNAs with ratios larger than 3.5 (i.e., most underexpressed), are as follows: (1) for LumA: none; (2) LumB: hsa-mir-1262, hsa-mir-202, hsa-mir-138-1; (3) Basal: hsa-mir-26a-1, hsa-mir-26a-2, hsa-mir-29c, hsa-mir-202, hsa-mir-101-1; (4) Her2: hsa-mir-135a-2, hsa-mir-202, hsa-mir-93.

For all the miRNAs that stand out with especially large differences in the discrete stable distribution parameter estimates, a further analysis needs to be done for exploring them as potential biomarkers for the (molecular subtypes of) breast cancer. As an example, for the miR-200 family, a more detailed analysis with robust testing is provided in Section 3.4.

4.4 TCGA-BRCA miR-200: an illustrative example for the miR-200 family

The miR-200 family is one of the most frequent groups of miRNAs whose expression is altered in (breast) cancer (e.g., Cavallari et al. 2021; Fontana et al. 2021; Wen et al. 2021; Ye et al. 2014). The family consists of five miRNAs: miR-141, miR-200a, miR-200b, miR-200c and miR-429. Despite comprehensive studies, the potential role of these miRNAs as biomarkers in cancer has not yet been completely understood. In this section, an illustrative example of differential expression analysis via discrete stable distributions is given for the miR-200 family for the PAM50 breast cancer molecular subtypes normal-like, LumA, LumB, Her, and Basal.

5.1 NOWAC-LUCA: estimates of discrete stable distribution parameters

Discrete stable distributions were fitted to 198 variables of miRNA expressions for data from all samples (n = 240) and to data from ‘case’ (n = 124) and ‘control’ (n = 116) subgroups. The summary of parameter estimates of the discrete stable distribution is given in Table 9.

For all miRNAs, the underlying distributions were estimated to be heavy tailed (i.e., α ̂ n < 1 ). This is true for all samples together, and when stratified by case-control status. In particular, for 50 % of miRNAS from all samples tail index is estimated less than 0.59. The estimates of λ for the data from all samples vary from 6 to 66 while around half are between 8 and 25. The minimum and 1st quartile of λ estimates for the data from cases and controls are similar while median, mean, 3rd quartile and maximum are quite different in comparison to the ones based on all samples.

5.2 NOWAC-LUCA: comparison of the goodness-of-fit

The goodness-of-fit of the Poisson, the negative binomial and the discrete stable distributions are compared using AIC. To compare goodness-of-fit, the ratios of corresponding AIC values were evaluated as before. If the ratio is less than 1, then the discrete stable distribution gives better fit. The ratios for p = 198 miRNAs are presented in Figure 6.

Figure 6:

NOWAC-LUCA: for p = 198 miRNAs, the ratios of AICs of discrete stable models to Poisson and negative binomial models from all samples (n = 240) and from samples grouped to case (n = 116) and control (n = 124).

From the left side of Figure 6 it can be seen that the discrete stable distribution was a better fit than the Poisson distribution for all miRNAs. In fact, more than half of the miRNAs have an AIC value for the Poisson distribution more than 10 times bigger than those for the discrete stable distributions. From the right side of Figure 6, for the data from ‘all’, ‘control’ and ‘case’ samples, the discrete stable distribution give a better fit than the negative binomial distribution for 58 %, 73 %, 54 % of the miRNAs, respectively. For the remaining miRNAs on the right side of Figure 6, the goodness-of-fit is similar.

5.3 NOWAC-LUCA: differential expression via discrete stable distribution

In this section, the differential expression analysis via the estimates of the parameters of discrete stable distributions for the p = 198 miRNA read counts is performed. First, the summary of data is given in Table 10. From this table, the quantiles of reads of the data from case and control samples are quite similar, except for the maximums.

The summary of the discrete stable distribution parameter estimates is given in Table 9. Although the quartiles of the estimates of α are quite similar, they hint that the tail index for the data from the ‘case’ samples are slightly smaller than for the data from the ‘control’ samples, that is, the miRNA-seq counts in the ‘case’ group have heavier tails. The quartiles of estimates of λ are lower (indicating underexpression) for the data from the ‘case’ samples, especially the 3rd quartiles and maximums. Next, for each miRNA, the differential expression analysis between cases and controls was carried out, as explained in Section 1.4. For each miRNA, the ratio of the λ estimates and corresponding difference in α estimates are presented in Figure 7.

Figure 7:

NOWAC-LUCA: for p = 198 miRNAs, the ratios of estimates for parameter λ and the differences in estimates for tail index α, between the ‘control’ and ‘case’ groups, labeled accordingly.

On the left side of Figure 7, the ratios that stand out represent underexpressions of certain miRNAs in the ‘case’ subgroup as compared to the controls. In particular, the top 5 underexpressed miRNAs (i.e., the largest ratios of λ estimates) in decreasing order are hsa-miR-10b-5p, hsa-miR-186-5p, hsa-miR-128-3p, hsa-miR-222-3p, hsa-miR-101-3p, hsa-miR-146b-5p, hsa-miR-98-5p, hsa-miR-342-3p. However, while these miRNAs have clear differences in the λ estimates, they also have the biggest differences in the tail index α estimates, which means the parameter λ differences should be treated with care. For more insights, a zoom into the left side of Figure 7, where only the part with differences in tail index estimates is less than 0.1, is presented on the right side of Figure 7. In Table 11, the miRNAs with ratios of λ estimates most different from 1, from the right side of Figure 7, and the corresponding differences in the tail index estimates, are listed.

For all these miRNAs we observe an underexpression of miRNAs among cases compared to controls, except for hsa-miR-142-3p, for which we observe an overexpression. In addition, results from an even more conservative approach regarding the tail differences, is presented in Table 12. In particular, from the right side of Figure 7, the top 20 ratios of λ estimates with constrained differences in tail index ≤0.01, is listed in Table 12.

Here, all ratios of the λ estimates represent underexpression of miRNAs in the ‘case’ subgroup. Although the ratios of the λ’s in both Tables 11 and 12 are not too remarkable, the listed miRNAs could be considered as potential biomarkers. Similarly, from the left side of Figure 7, miRNAS with the largest relative differences in λ estimates as well as differences in α estimates should be studied further as potential biomarkers for lung cancer.

In (Nøst et al. 2023) nine candidate miRNAs were proposed as potential early markers of lung cancer; miR-320d, miR-320c, miR-320b, miR-92b-3p,miR-130b-3p, miR-200c-3p, miR-375-3p, miR-335-5p,and miR-323a-3p. This study found miR-320c and miR-320 elevated in pre-diagnostic specimen, i.e., they can be suggested as early markers of lung cancer. Based on Table 12, we would in addition highlight hsa-miR-335-5p (and hsa-miR-320a-3p from the miR-320 family).