Rank‐based Bayesian variable selection for genome‐wide transcriptomic analyses

Variable selection is crucial in high‐dimensional omics‐based analyses, since it is biologically reasonable to assume only a subset of non‐noisy features contributes to the data structures. However, the task is particularly hard in an unsupervised setting, and a priori ad hoc variable selection is still a very frequent approach, despite the evident drawbacks and lack of reproducibility. We propose a Bayesian variable selection approach for rank‐based unsupervised transcriptomic analysis. Making use of data rankings instead of the actual continuous measurements increases the robustness of conclusions when compared to classical statistical methods, and embedding variable selection into the inferential tasks allows complete reproducibility. Specifically, we develop a novel extension of the Bayesian Mallows model for variable selection that allows for a full probabilistic analysis, leading to coherent quantification of uncertainties. Simulation studies demonstrate the versatility and robustness of the proposed method in a variety of scenarios, as well as its superiority with respect to several competitors when varying the data dimension or data generating process. We use the novel approach to analyze genome‐wide RNAseq gene expression data from ovarian cancer patients: several genes that affect cancer development are correctly detected in a completely unsupervised fashion, showing the usefulness of the method in the context of signature discovery for cancer genomics. Moreover, the possibility to also perform uncertainty quantification plays a key role in the subsequent biological investigation.


Introduction
Over the last decades, technological advances have generated an explosion of data in a variety of fields. High-throughput technologies have made -omics data more accessible to researchers, allowing for a systematic exploration of the genetic and epigenetic basis of cancer. However, -omics data pose several significant analytic challenges, including high-dimensionality, high complexity and nonnormality and thus, standard statistical approaches are not applicable. To overcome these compli-cations in high-dimensional omics-based analyses, it is considered biologically reasonable to assume that only a small piece of information is relevant for prediction or subtyping [3,14,15], so that a lowdimensional solution can also lead to better results interpretability. However, a frequent approach to variable selection for -omics applications consists in performing a priori ad hoc selection of genes according to in-sample statistics or outer literature information, a strategy that heavily affects the solidity and reproducibility of results [16]. Moreover, it has been established in the last decades, with the explosion of high-dimensional data methods, that variable selection should be performed jointly with the main inferential tasks, to properly estimate the uncertainties and the impact of selection on inference results [8]. Among unsupervised methods for variable selection in the context of -omics data, Bayesian model-based approaches are highly suitable, due to their capability of propagating all uncertainties. While unsupervised methods may lead to inaccurate estimates because of the low signal-to-noise ratio [17], especially for high-throughput -omics data, the Bayesian approach allows for incorporating prior biological knowledge when estimating relevant genes from the data. Current methods for Bayesian variable selection for high-dimensional unsupervised problems include probabilistic sparse principal component analysis [2,7], and Bayesian graphical structure learning [22,24], all with applications to -omics data. In the case of ranking data, there currently exists no Bayesian unsupervised methods capable of tackling the dimensions of -omics data.
Converting continuous data to rankings and then performing rank-based analyses is a practice that has gained much popularity in statistical genomics [11]. Different -omics data layers can have incomparable absolute scales, and the use of ranks leads to robust inference [1]. Moreover, analysing rankings instead of continuous variables allows for easy integration of multiple heterogeneous data, and enhances the reproducibility and interpretability of results. The analysis of rank and preference data has recently seen a growing interest in the machine learning community [4], with the Plackett-Luce [30,33] and Mallows model [31] being among the most commonly used models. The Mallows model is an exponential family-type model on the permutation space, where the probability density depends on a distance between rankings. Compared to Plackett-Luce, the Mallows model shows greater flexibility in the choice of the distance between permutations, and it is also more versatile in adapting to different kinds of data. As a probabilistic model for rank data, the Mallows model enjoys great interpretability, model compactness, inference and computational efficiency. Our current work is based on the Bayesian Mallows Model proposed in [38], whose implementing solution, BayesMallows, is described in [35]: BayesMallows already provides an inferential computationally feasible approach for the most important choices of distance among permutations, and it shows good accuracy when compared to competitors on datasets of moderate size [28]. However, for this method to be applicable to the common data dimensions in -omics applications, variable selection is crucial. Moreover, from a purely modeling perspective, imposing the Mallows data generating process on the full (genome-wide) ranked list of genes seems highly unrealistic: RNAseq measurements of moderately expressed genes tend to be very noisy, thus making a low-dimensional solution more suited to this kind of data. Imposing a ranking model only on a selection of genes allows focusing on consistently highly expressed genes, since those will be selected as highly ranked items, thus allowing to detect genome-wide high-quality signals in a more robust way. More generally, this leads to an unsupervised variable selection procedure that can learn from the available rankings to distinguish relevant items from background ones in a high-dimensional setting, as well as provide an estimate of the position in the ranking of the selected relevant items.
The aim of this paper is to develop a novel extension of the Bayesian Mallows model, the lowerdimensional Bayesian Mallows Model (lowBMM), providing a low-dimensional solution for rankbased unsupervised variable selection in high-dimensional settings. To the best of our knowledge, the only proposal similar in scope to ours is [42], proposing a partition-Mallows model (PAMA) for rank aggregation that can distinguish relevant items, as well as ranking the relevant ones. This method can handle partial lists, incorporate any available covariate information and can estimate the quality of rankers, as it is based on BARD [11], which assigns aggregated ranks based on the posterior probability that a certain item is relevant. However, PAMA only considers the Kendall distance for the Mallows model, thus allowing the partition function to have a closed form, while the Bayesian Mallows model can use any distance. Additionally, PAMA is unable to perform individuallevel inference, a typical advantage of the Bayesian Mallows model, and it is tested on a relatively complete selection of simulated scenarios, however lacking large-data examples. Concerning rankbased variable selection in -omics analyses, the recent work in [9] converts microbiome data to rankings and develops a Bayesian variable selection model using ranks to identify the best covariates for modeling the outcome. However, this latter method is fully supervised, and moreover, common for both of these methods is the lack of applicability to high-dimensional settings: all data examples are of much smaller dimensions compared to the typical dimension of -omics data (> 10 4 ). Other Bayesian implementations to analyse rank data have been proposed to reduce computational burden [21,40,41], to merge rankings from multiple sources [11,20], and to combine rank data with other types of data [5]. There exist other methods which are not directly related to lowBMM, but can handle rank data, such as the Markov chain-based methods in [26,27] as well as the Mallows based methods in [18,25]. Simulation studies show that the proposed approach is superior to existing methods in high dimensions, and has no competitor in ultra-high dimensions.
The remainder of the paper is laid out as follows: in Section 2 we describe the novel lowBMM for variable selection, and its computational implementation. Several simulation studies are described in Section 3, aimed at demonstrating the versatility and robustness of the novel method under different scenarios. Included here is also a sensitivity study on the tuning parameters involved in the model, and a comparison with existing methods. In Section 4, we present an application of this modeling approach to RNAseq gene expression data from ovarian cancer patients. We conclude with a brief summary and discussion in Section 5. the form where α > 0 is a positive scale parameter, ρ ∈ P n is the latent consensus ranking, Z n (α, ρ) is the model partition function, 1 S (·) is the indicator function of the set S, and finally d n (·, ·) : P n × P n → [0, +∞) is a right-invariant 2 distance function between two rankings. We have here chosen to let the distance function explicitly depend on n, so that d n (·, ·) denotes the distance computed on the complete set of items. Several possibilities for choosing this distance function exist, such as the footrule distance, the Spearman distance, and the Kendall distance [12]. In this paper we choose to use the footrule distance, defined as d n (R, ρ) = n i=1 |R i − ρ i |, the equivalent of an 1 measure between rankings. The choice of this distance is motivated by its greater computational efficiency compared to its competitors, such as the Kendall distance, which is more computationally intensive [35]. Moreover, the popularity of the Kendall lies in the fact that a closed form of the model partition function Z n (α, ρ) exists in this case; however, we have previously proposed computational solutions to be able to use other distances [35], while also showing that the Kendall is often less accurate [28]. The Spearman distance is also problematic in our case, as it tends to largely penalize extreme data, being the equivalent of an 2 measure between rankings. On the other hand, the equivalent of the 1 between continuous data seems more suited to -omics applications.
Since the footrule distance is right-invariant, the partition function Z n (α, ρ) is independent of ρ and only dependent on α. Moreover, we assume in the remainder of the paper that α is fixed and known, and we provide strategies to tune this parameter off-line. This simplification is done to avoid having to compute the normalizing constant Z n (α) at each iteration of the Markov Chain Monte Carlo (MCMC) algorithm (see subsection 2.3), thus greatly lightening the computational burden of the algorithm, at the cost of the tuning of a single scalar parameter, which proves to be straightforwardly carried out in practice (see subsection 2.4 for further details). Nonetheless, the performance of the approximated approaches to estimate the normalizing constant Z n (α) introduced in [38] greatly depends on the size of n, with the performance degrading for decreasing n, and with the approximation becoming increasingly computationally demanding when n → ∞. Thus, fixing α simplifies the inferential procedure quite substantially.
Given the model in (1) and the choices specified above, the likelihood associated to the observed rankings R 1 , . . . , R N under the Mallows model takes the form: In order to use the Bayesian version of this model as introduced in [38], we need to specify a prior for ρ. In this paper, we set a uniform prior π(ρ) = 1 n! 1 Pn (ρ) in the space P n of n-dimensional permutations. The posterior distribution for ρ then becomes To perform inference, [38] proposed a MCMC algorithm based on a Metropolis-Hastings (MH) scheme (see [35] for details on the implementation).

Lower-dimensional Bayesian Mallows model (lowBMM) for variable selection
We here introduce a generalization of the BMM presented in Section 2.1 to the scopes of variable selection, which in the context of models for rankings assumes the form of a selection of the relevant items, i.e., the items worth being ranked. We are thus thinking of generalizing BMM to make it suitable to situations in which n is large or ultra-large, and therefore it becomes unrealistic to assume that all n items can be completely ranked. We define A * = {A i1 , ..., A i n * } as an n * −dimensional subset of the original set of items A, with n * << n and A * ⊂ A. When defining the lower-dimensional Bayesian Mallows model, the underlying assumption is that only a portion of the data follows the Mallows distribution, while the rest is unranked: this is a very realistic assumption, especially in the large n setting. We then formally assume that items in A * are such that R j | A * , j = 1, . . . , N, follows a Mallows model over the permutation space of dimension n * , where with the notation R j | S we indicate the restriction of an n-dimensional ranking to only the items belonging to the set S, which is also a ranking in dimension |S|. The remaining items in A \ A * are simply irrelevant to the scopes of the analysis, and therefore show no specific pattern in the data. This is equivalent to assuming that R j | A\A * , j = 1, . . . , N, is uniformly distributed over the space of permutations of dimension n − n * . A * is then the novel model parameter characterising the lower-dimensional Bayesian Mallows model, the variable selection setting for the Bayesian Mallows model. Not only, lowBMM is assuming the data follow a Mallows model only on the lower-dimensional set of items, so that the consensus ranking parameter is now defined as ρ ∈ P n * . Therefore, the likelihood of lowBMM is defined on a lower-dimensional space of dimension n * , and it takes the form: . This is the same distance as d n (·, ·) but restricted to the n * −dimensional set of items included in A * . The vector R j | A * is the restriction of R j to the set A * : it thus defines a partial ordering from which an n * −dimensional permutation can be derived. Moreover, U S is the uniform distribution over the domain S, meaning that the term U P n−n * (R j | A\A * ) refers to the assumption that R j | A\A * includes only noisy unranked data.
Remark. The way the modeling framework is defined does not result in a variable selection procedure that selects the top-ranked items only. Indeed, the items in A * do not need to be those that are most often top-n * ranked in the data, even if this is the most intuitive solution to a lower-dimensional ranking problem. The relevant items could also show a "consistency" pattern, i.e., be often ranked in the same respective order across assessors: this latter possibility constitutes a valid data generating process for lowBMM, as proved in the simulation study in Section 3.1, as the only requirement for lowBMM to be a suitable model for the data is that the ranks of the items follow a Mallows model in the lower-dimensional space P n * . In other words, the model is able to pick out the n * items on which the assessors agree the most, as the model assigns the largest probability to the items whose ranks show the smallest distance to the consensus, which is constrained to be an n * -dimensional permutation. In conclusion, the results from lowBMM will depend on the specific dataset: if there exists a top-rank solution in the data, such a solution will be estimated by lowBMM as the one with the largest posterior. On the other hand, if there exists a pattern in the data showing certain items consistently ranked in a specific order (not necessarily top-ranked), lowBMM will give those items the largest probability in the marginal posterior distribution of A * , and rank them accordingly in the marginal posterior distribution of ρ.
Since the inferential approach is Bayesian, we have to decide priors for all parameters. As we did previously, we set a uniform prior on ρ, however restricted to the n * −dimensional space of permutations of elements of A * : π(ρ|A * ) = 1 n * ! 1 P n * (ρ). We also set a uniform prior for A * over C, π(A * ) = 1 |C| 1 C (A * ), where we define C as the collection of all n n * possible sets of items of dimension n * chosen from a set of dimension n. The posterior distribution of the variable selection for the Bayes Mallows model can then be written as which, by removing the terms not depending on any of the model parameters and by assuming the data are proper rankings, can be simplified to Note that, from the posterior distribution in (3), marginal posterior distributions of both model parameters, ρ and A * , can be easily derived. Often one is interested in computing posterior summaries of such distributions. By inspecting the marginal posterior distribution of A * estimated by lowBMM, one can for instance a posteriori select the items in A that maximise this marginal posterior under all possible model reductions of dimension n * , thus practically attaining variable selection via the marginal posterior mode (i.e., the maximum a posteriori -MAP) of A * . However, computing the MAP might not be the best way of summarizing the marginal posterior distributions of both ρ and A * , in addition to being a computationally intensive procedure for larger n and n * . Instead, we introduce novel ways for providing posterior summaries of both the consensus ranking of the relevant items, calledρ A * , and of the set of relevant items, calledÂ * ; these approaches are described in Section 2.3.1.
Furthermore, the inspection of the marginal posterior distributions of ρ and A * has much wider implications than simply obtaining posterior summaries: for example, by inspecting the marginal posterior distribution of A * one can decide whether variable selection is appropriate to the data at hand (few items are clearly assigned rankings while the rest shows high uncertainty), or rather items are a posteriori ranked in blocks (several items share the same rank with large probability, and only a respective ordering between groups of items can be properly estimated). Scope of the simulation studies in Section 3 is to assess the correct behaviour of the model in these different scenarios, together with its gain in efficiency in the large n context.
Remark. Note that we have assumed n * to be fixed and known. Estimating n * together with the other model parameters would imply a model selection step, which would make the inferential task substantially more challenging, or require to completely change our inferential approach, for instance towards the use of stick-breaking priors [19]. Both possibilities are completely out of the scope of the present paper and are left to future speculations. On the other hand, the choice of fixing n * is typically not problematic in practice, as one can often tune it in connection to the real application (e.g. in genomics, biologists are often interested in a specific "number" of relevant genes, even more if we consider that each gene-regulating biological pathway includes on average a hundred genes [34]).
Another practical possibility is to simply choose n * as large as allowed by the computing resources at hand. This latter approach is also supported by the simulation studies (see Section 3.2), where we realized that choosing n * larger than necessary would result in a solution including and automatically detecting the lower dimensional relevant subset.

Metropolis-Hastings algorithm for inference in the lowBMM
In order to obtain samples from the posterior in equation (2), we set up a MH-MCMC scheme. The algorithm iterates between two steps: in one step the consensus ranking ρ given the current set A * is updated, and in the other the set A * is updated given the current consensus ranking.
In the first step of the algorithm, we propose a new consensus ranking ρ ∈ P n * using the "leapand-shift" proposal distribution described in [38]. The acceptance probability for updating ρ in the MH algorithm is In equation (4) above, P l denotes the probability mass function associated to moving from one n *dimensional ranking vector to another according to the leap-and-shift proposal. The parameter l denotes the number of items perturbed in the consensus ranking ρ to get a new proposal ρ , and is used to tune the acceptance probability. To compute the distances restricted to A * , the R j 's needs to be updated so that they correspond to the items in the current set A * . The second step of the algorithm updates the set A * . We propose a new set A * prop by perturbing L ∈ {1, ..., n * } elements in the current A * , selected with uniform probability. The L items are swapped with L items from the set A \ A * , again uniformly. Therefore, we can formally write the proposal distribution as q( where we define C in L as the collection of all n−n * L possible sets of items of dimension L chosen from a set of dimension n − n * , to be brought "in" to the reduced set A * in order to obtain A * prop . Likewise, we define C out L as the collection of all n * L possible sets of items of dimension L chosen from a set of dimension n * , to be brought "out" of the reduced set A * in order to obtain A * prop . Note that such proposal distribution is by definition symmetrical. The move from A * to A * prop is then accepted with probability: .
Since the prior for A * is uniform, the prior terms in the ratio simplify, as do the symmetrical proposal distributions. We are then left with a ratio of likelihoods, that can be written as the distance is computed on the restricted set of the current A * prop , and R j is also restricted to the same set. Therefore, one needs to decide how to evaluate ρ on A * prop . There are multiple ways of doing this, the most simple being to randomly perturb the current ρ from the first step of the algorithm to accommodate the new L items proposed for A * prop . That is, given the current ρ, we discard the ranks corresponding to the L items that were left out of A * prop , and we randomly re-assign those same ranks to the new items.
The MH-MCMC algorithm described above is summarized in Algorithm 1.

. Postprocessing of the MCMC results
In order to derive posterior summaries of the consensus ranking of the relevant items and of the set of relevant items, namedρ A * andÂ * respectively, we inspect the marginal posterior distribution of ρ and A * (after burn-in). Suppose M posterior samples are obtained: For the computation of the posterior summaries, we first need the following definitions. Definition 1. Given a vector of real numbers (x 1 , ..., x n ) ∈ R n , its corresponding rank vector is obtained as follows: Definition 2. For n * ≤ k ≤ n, the Highest Probability Set (HPS) of A * is defined as follows: for all We now quantify the two posterior summaries of ρ and A * as follows: Note that in all the figures in the paper reporting heatplots of the marginal posterior distribution of ρ, the items on the x-axis will be ordered according to their respective ranking inρ A * .
In some cases, it might be of interest to obtain a "top probability selection within the selection" A * top ∈Â * for further downstream analysis (for instance in cases when not only n but also n * is quite large). The strategy for computingÂ * top consists in inspecting the probability distributions P (A i ∈ top-K, i = 1, ..., n | A i ∈Â * ) for varying K ≤ n * . From these probability distributions, we would like to choose a K that discriminates the items that have a high probability of being selected compared to the rest of the items. Once the K is specified, the probability distribution will by construction be bimodal, and it will be possible to find a cut-off point c that discriminates its two peaks well (this has explicitly been explained in the context of the case study presented in Section 4, see Figure 11). The cut-off c is then used as a lower probability bound for being included in

Off-line estimation of α
To simplify inference by avoiding the computation of the normalizing constant, we decided to fix α. Nonetheless, we would still need an estimation of its most reasonable value when we are given a new ranking dataset. We propose one possible method for estimating such a value, inspired by [29].
As α describes how closely the rankings are distributed around the common consensus, or in other words, how similar the individual rankings are to each other, we can assume that a higher α indicates a higher "similarity" between assessors. Therefore, when assuming that the individual rankings of the assessors are drawn from a Mallows with parameter α, a higher α results in a lower mean distance between the assessors. To estimate α, we assume that the mean distance among N assessors for a real ranking dataset should resemble the mean distance of a simulated ranking dataset of the same dimension, generated by drawing independent samples from the Mallows distribution. Precisely, suppose we want to estimate α for a real ranking dataset with n items ranked by N assessors: we can then generate several ranking datasets of dimension (N, n) over a grid of α 0 values, R α0 1 , . . . , R α0 N (the chosen consensus ranking ρ 0 does not matter). Then, we calculate the mean pairwise distance between assessors for each simulated dataset (i.e. for every α 0 ), defined as: Finally, we compute the same mean distance between assessors as in (7) for the dataset we are interested in, and choose the α 0 value where the two distance measures intersect. The approach to estimate α described so far works for a complete ranking dataset where we assume that all n items are ranked according to independent draws from the Mallows distribution. However, when using the lowBMM, we will be assuming that only n * items are ranked according to a Mallows model. Hence, when using the method described above we will need to re-scale the estimated α from dimension (N, n) to dimension (N, n * ). Letα n be the optimal α 0 value estimated from the datasets of dimension (N, n), as described above. We can then re-scaleα n to obtain the optimalα n * value in dimension n * by matching the terms in the exponent in the Mallows model: where max dn (respectively max d n * ) is the maximum attainable footrule distance between two rankings in dimension n (respectively n * ). The method is assessed in a simulation study described in Section 3.5.

Simulation experiments
This section describes how lowBMM performs under different data generating processes, how its performance is affected by the various parameters involved, and how it compares to other methods. The variable selection procedure provides a posterior distribution for the consensus ranking in a reduced dimension, and the type of distribution obtained in turn depends on the data generating process. Therefore, inspecting the solution provided by lowBMM under different scenarios for generating the data allows us to understand the data structure and patterns in real-life situations. First, we consider a top-rank data generating process, where only a subset of the items is ranked, and assigned the top ranks in reduced dimension, while the rest is noisy: this scenario results in a top-rank solution, i.e., lowBMM correctly selects the top-items following a rank model. Second, we generate data where items are consistently awarded ranks in a specific order: this results in a rank-consistency solution, i.e., items that can be ordered are selected and top-ranked. A study on how lowBMM works under different data generating processes is presented in Section 3.1, and a sensitivity analysis on the tuning parameters involved in lowBMM can be found in Section 3.2. We also test the method's robustness to the noise level in the data in Section 3.3, and we provide an extensive study on how lowBMM compares to existing methods in Section 3.4. Finally, we present a performance assessment of the off-line estimation of α in Section 3.5.

Varying the data generating process: top-rank and rank consistency experiments
To test the method's applicability to different data structures, we consider two data generating processes: top-rank and rank-consistency. In the case of top-rank we generate data in the following way: n * items were given ranks sampled from the Mallows model, .., n * ) while the rest of the items were randomly assigned ranks, R j | A\A * ∼ U(P n * +1,...,n ).
Here the top items are ranked according to the Mallows model, and thus result in a top-rank solution.
In the rank consistency experiment on the other hand, the aim is to see whether lowBMM can select consistently ranked items even if not necessarily top-ranked. The setting was the following: items in A * were given a respective ranking sampled from a Mallows model in dimension n * ; then, for obtaining each assessor-specific ranking, the same items were assigned a random ranking in dimension n in their respective positions; finally, the rest of the items had randomly assigned ranks. Therefore, items in the true A * were always ranked in a certain order for each assessor, however their specific ranks were not consistent. For the top-rank data generating process, we performed two data simulation experiments: first a toy example with n = 20 items, and then a larger example with n = 1000 items.
In the small example, we considered a group of N = 50 assessors, and n * = 8 relevant items. For what concerns the tuning parameters, we set L = 1 and l = round(n * /5) as suggested by the tuning parameter study in Section 3.2. We used M = 5 · 10 3 MCMC iterations, and the algorithm took 10 seconds to run. The left panel in Figure 1 displays the marginal posterior distribution of ρ where the items on the x-axis have been ordered according toρ A * (see Section 2.3.1 for the computation of posterior summaries). The results show that the method is able to correctly select the relevant items to be included in the selection set, as the items whose ranking was simulated from the Mallows model are all top-ranked. However, there is uncertainty associated to their specific rankings: some items seem to "compete" for ranks, e.g. item 6 (true rank 2) and item 18 (true rank 3) both show a large posterior marginal probability of being ranked 2; at the same time, their final ordering is correct (the ordering on the x-axis is consistent with the rainbow-grid representing the true ordering ρ A * ). This is not an unexpected behaviour, as the variability in the items' true ranks is dependent on the scale parameter, α. As we do not estimate α in the current version of lowBMM, we inspected how the solution changed when varying this parameter. The α parameter regulates how much the assessors agree with each other on the ranking of the items: a large α means a better agreement between the assessors, and a lower α means less agreement. So, when α is large, the items distributed according to the Mallows are easier to pick out, making it less likely that other items are accepted in the proposal set. This is evident in the distribution of the bars on the top of the plots in Figure S1 in the supplementary material, where the larger α value shows less variability in the selected items. Moreover, when the data is "easier", the algorithm converges much more quickly to the correct set. The bottom items (whose ranks are essentially noise) are clearly distributed on the bottom for α = 10, while there is a bit more uncertainty for α = 3, as expected.  18 1 10 7 20 16 17 8 11 14 13 5 2 9 3 4 12 18 1 10 7 20 16 2 3 4 5 8 9 11 12 13 14 15   In order to test the method in a more realistic scenario, we simulated the data in a dimension closer to that often encountered in -omics applications. This larger experiment aims to investigate the method performance and its computing capabilities on a more demanding dataset when using the same data generating process as before (top-rank). We set the total number of items to n = 1000, the number of relevant items to n * = 50 and the number of assessors to N = 50. We also set the tuning parameters to L = 1 and l = round(n * /5), as previously done. We ran the MCMC for 5 · 10 4 iterations in this larger experiment, and the algorithm used 155 seconds to run. Note that we used a larger M in the larger examples to ensure good space exploration and convergence of the chains, and the computing times were all scaling well with the increasing dimension. The posterior of ρ can be seen in the left plot of Figure 2: more than 75% of the items in the true selection are assigned top ranks with large marginal posterior probability, while the uncertainty associated to the items' rankings increases moving towards the bottom, for the items not in the top-n * . The trace plot on the right in Figure 2 demonstrates convergence, with most top-ranked items clearly converging to one final rank, and only a few items competing for the same rank throughout the chain. This "competing" behaviour should to some degree be expected, however we noticed it to be enhanced by increasing the tuning parameter L (see Figure S9 in the supplementary material), suggesting to keep L small (this is confirmed in the sensitivity study on the tuning parameters in Section 3.2). For the rank consistency data generating process, we carried out a simulation study similar to the small top-rank example, also here focusing on a small experiment with n = 20, n * = 8 and N = 50. The right panel in Figure 1 presents the marginal posterior distribution of ρ where the items have been ordered according toρ A * on the x-axis. Here the performance of lowBMM is evident: the correct items (indicated by the rainbow grid on the top) are given the highest marginal posterior probability in A * (as shown by the bar plot on top), and the ordering of the items is also consistent with the true ρ A * . This shows that lowBMM is able to handle varying data structures, and quite interestingly, is performing slightly better on a dataset with consistently ranked items, compared to the top-rank data example. The method also seems to be somewhat more robust to the choice of the tuning parameters for this type of data generating process (see Section 3.2, Figure 3 for details). The rank consistency simulation experiment can be considered a more realistic example of how ranking data could be structured in real-life situations, particularly in genomics: a subset of genes might be consistently ranked in a certain pattern across samples (e.g. consistently more or less expressed), while their individual rankings do not really matter as much. At the same time, such a pattern would be much more interesting to detect than the actual true ranking estimation. Furthermore, experiments making use of this data generating process were overall more challenging to handle for all the other competing methods considered for comparison (see Section 3.4).
In the rank consistency simulation experiment, we also experimented by reducing the number of items being given the consistent ranks, n * , as well as the level of precision in the data sampled from the Mallows, α. Decreasing both n * and α makes the estimation via lowBMM more difficult, as less items show a detectable pattern, and the noise level in the data is higher. Figure 3 shows the results obtained with α = 0.5: as expected, we see a much larger uncertainty in the ranks with respect to what observed in Figure 1 (right) for α = 10. It is interesting to note the performance of the method: despite the much larger uncertainty, still most weight in the marginal posterior of ρ is given to the correct items (looking at the bar plots at the top), and the ordering of the items is consistent with the true ordering (rainbow grid on top). This shows that the accuracy of lowBMM is unaffected by decreasing α or n * , even in the case of a rank consistency data generating process. Note that in this experiment we have used the same n * both in the data generation and in the analysis, as the main point here was assessing the robustness of lowBMM to a less detectable pattern. In the next Section 3.2 we will instead verify the method robustness to the misspecification of n * .

Sensitivity study on tuning parameters
A sensitivity study was carried out to verify the effect of the tuning parameters L and l on lowBMM, as well as the effect of misspecifying the number of relevant items for the selection, n * . We first focused on studying the effect of L and l and therefore we set: n = 100, N = 10, n * = 10 and α = 5. The data was generated according to the top-rank data generating process described in Section 3.1: n * items were given ranks sampled from the Mallows model, while the rest of the items were ranked randomly lower. We experimented by varying the two tuning parameters introduced in Section 2.3: L, i.e., the number of items perturbed in A * for proposing a new set A * prop , and l, i.e., the number of items perturbed in the current ρ for proposing a new consensus ρ . Values for L and l were explored over a grid (L, l) ∈ [1, 2, 3, 4, 5] × [1,2,3,4,5]. Results were collected after running the method on 20 simulated datasets, when using M = 5000 MCMC iterations, which always showed to be more than enough for convergence.
To evaluate the results we computed the footrule distance between the true and the estimated consensus ranking, normalized by the number of correctly selected items n corr = |A * Â * |: d norm (ρ A * ,ρ A * ) = d(ρ A * ,ρ A * )/n corr (see Section 2.3.1 for details on the computation of the posterior summariesρ A * andÂ * ). For brevity, we will sometimes indicate d norm (ρ A * ,ρ A * ) simply with d norm .
The left panel in Figure 4 shows d norm (ρ A * ,ρ A * ) over the 20 runs and for varying L and the left panel in Figure S2 in the supplementary material displays the proportion of correctly selected items in A * for each L. The results showed a clear indication that a smaller L results in a lower error overall. Similarly, the middle panel in Figure 4 (and also the middle panel in Figure S2 in the supplementary material) indicates that the leap size l has a smaller effect on the error, even if l = round(n * /5) yields a slightly better result, both in terms of the footrule distance (d norm ) and in terms of the selection of items. Similar conclusions can be made by inspecting the heatplots of the marginal posterior distributions of ρ for varying values of the tuning parameters (see Figure S3 in the supplementary material for α = 3, and Figure S4 in the supplementary material for α = 10), where it is evident that varying l does not affect the error (but does contribute slightly to the degree of uncertainty in the ranks), while varying L does seem to affect the error to a slight higher degree compared to l. Furthermore, as noticeable in the barplots in Figures S3 and S4, increasing L contributes to worsening the degree of exploration of the items space when estimating the posterior distribution for A * , thus supporting the decision to keep L small. It is also worth exploring variations in the acceptance probability of the two parameters estimated in the MCMC, ρ and A * , when varying the tuning parameters L and l. From inspection of Figure  5 (α = 3, n * = 8, n = 20 and N = 50 in this scenario), it seems that the acceptance probability for ρ increases with increasing L, while it decreases with increasing l; the latter behaviour is expected, while the former is less trivial to understand. For what concerns A * , the acceptance probability decreases greatly with increasing L, and it seems to not be much affected by the leap size l, which are both quite natural behaviours. Therefore, we generally recommend to keep L = 1 to explore the space appropriately when estimating A * , and to tune l by inspecting the acceptance probability of ρ. The tuning of l in the BMM has been extensively studied, resulting in a recommendation to use l = round(n/5) as a default value to ensure a good balance between space exploration and accuracy (see [38] for further discussions). However, the respective values of n * and n, which constitute the novelty of lowBMM with respect to the original BMM, might make it challenging to choose a combination of tuning parameters that results in a chain that explores the space well. Therefore, even if l = round(n * /5) will generally be a solid choice, we do recommend to try some values of l before setting it.
Since the number of relevant items n * is in general unknown, we also tested the performance of lowBMM when assuming a different n * than what was used when generating the data. We did a systematic study with n * guess ∈ [6,8,10,12,20], and with the truth set to n * = 10. The results from the study are displayed in the right panels in Figure 4 and in Figure S2 in the supplementary material, both suggesting that setting n * guess > n * increases the accuracy. This is not surprising, as it allows the algorithm to search for a good selection in a larger permutations space. It is also interesting to evaluate how the consensus ranking ρ changes when n * guess differs from the true value: the marginal posterior distribution of ρ obtained when n * guess < n * is displayed in Figure 6 on the left, and when n * guess > n * on the right. By inspecting the left panel we see that assuming n * guess < n * yields a clustered structure, thus implying that the best lower-dimensional solution of ranked data in a higher-dimension is consistent with a partial ordering, also called "bucket solution" [13]. On the other hand, assuming n * guess > n * (Figure 6, right panel) yields a top-rank solution, as the algorithm seems Posterior probability Figure 6: Results from the sensitivity study on the tuning parameters as described in Section 3.2, here when setting n * guess differently from the true n * = 8. Both panels display heatplots of the marginal posterior distribution of ρ with the items ordered according toρ A * on the x-axis; left: n * guess = 4, right: n * guess = 12. The rainbow grid indicates the true ranking ρ A * , and the bar plot indicates the proportion of times the items were selected in A * over all MCMC iterations. n = 20, N = 50, α = 3, L = 1 and l = round(n * guess /5).

Robustness to noise experiment
In the simulated data experiments conducted so far, we only tested a varying degree of agreement among the rankings provided by the assessors by varying α. It could then be argued that both data generating processes were somewhat easy to handle for lowBMM, since the noise structure was only induced by controlling a single model parameter. To test the robustness of lowBMM to increasingly large noise levels, regardless the value of α, we performed another simulation experiment: data were generated as previously, but an increasingly large random perturbation of the rankings was then enforced in the data. Precisely, we first simulated the data according to the top-rank data generating process described in Section 3.1 (with n * = 8, n = 20 and N = 25), and then we iteratively swapped the rankings of items in the top n * (sampled from the Mallows) with those of items ranked below. In order to perform this iterative rank swapping, we followed the following procedure: at level 1, we take the "bottom item" in the true set A * (the item whose rank is n * in the true ρ A * ), and we swap the rank for that item with the rank of an item outside of the true set, for 90% of the assessors. Then, at any level i > 1, we repeat the same procedure for the items in the true set A * whose ranks in ρ A * are n * − i. So for instance, if we repeat the swapping procedure four times, at level 4 the ranks of four items in the true set have been swapped with the ranks of four items from the outside. The tuning parameters were distributed over the grid (L, l) ∈ [1, 2, 3] × [1, 2, 3].
The results can be seen in Figure 7, displaying boxplots of d norm (ρ A * ,ρ A * ). The results generally show a very good performance of lowBMM even with an increasing level of noise, i.e., several levels in the iterative perturbation/swapping. The figure also suggests that a lower L is more robust to noise, while the leap size l has very little effect: therefore, studying lowBMM's robustness to noise also confirmed the conclusions about the tuning of L and l reached in Section 3.2.

Comparison with other methods
To further assess the accuracy and solidity of the lowBMM method, it is also important to compare its performance with competitor approaches. In order to carry out a fair comparison, we set up simulation experiments with varying computational demands, taking inspiration from the simulation study in [42]. The first experiment is a toy example with a small number of items, n = 20, so that all compared approaches can manage to provide a solution in a small setting. Two additional experiments with increasing dimensions, instead, are used to test how the various methods are able to handle a larger number of items, specifically: n = 100 and n = 1000. The methods included in the comparison were: four Markov chain-based methods MC 1 , MC 2 , MC 3 [26], and Cross Entropy Monte Carlo (CEMC) [27]; three Mallows-based methods: the original Mallows Model (MM) [31], the Extended Mallows Model (EMM) [25], and the Partition Mallows model (PAMA) [42]. Where possible, a comparison with the original BMM [38] was also included. Finally, we included BORDA [10] as a reference, as it is a simple method and quite commonly used as a "baseline" comparison. The TopKLists R package was used for the MC methods, CEMC and BORDA; PerMallows was used for MM; ExtMallows was used for EMM; the PAMA package was used for PAMA; and finally the BayesMallows package was used for BMM. All methods were implemented for all simulation experiments, however some methods did not manage to complete/converge by the set time limit.
We computed three performance measures in order to compare and evaluate the methods. Let ρ A * be the true consensus ranking of the relevant items, andρ A * be its estimate (the point estimate for the frequentist methods, and the default posterior summary for the Bayesian ones). Similarly, let A * be the true set of relevant items, andÂ * be its estimate, with |Â * | = |A * | = n * . We then consider the following measures of performance: 1. the footrule distance between the true and the estimated consensus ranking, normalized by the number of correctly selected items n corr , denoted as d norm (this is the same measure as described in Section 3.2).
3. the recovery distance d R = d τ (ρ A * ,ρ A * ) + (n * − n corr ) × (n + n * + 1)/2 as defined in [42], where d τ is the Kendall distance between two rankings, and n * −n corr is the number of wrongly selected items.  The results from the first experiment with n = 20 items are shown in the upper portion of Table 1, and in the left panel of Figure 8. BORDA, EMM and lowBMM all perform similarly, outperforming the other methods in terms of accuracy, both in the selection of the items, and also in their respective rankings. Similar results were obtained for the experiment with n = 100 items, where however we were only able to run the methods displayed in the right panel of Figure 8 in the set time limit of 72hrs (see also the lower portion of Table 1). Finally, concerning the large experiment with n = 1000 items, we simulated data according to the two data generating processes (top-rank and rank consistency) described in Section 3.1. Results are reported in Table 2: only the Mallows-based methods, and not the PAMA, were able to converge in these scenarios. The results for the top-rank simulation experiment are similar to those obtained for the smaller set of items, however with EMM performing slightly better in terms of d R andp. On the other hand, EMM was by far the slowest method in this computationally demanding scenario, as the algorithm is based on an iterative optimization procedure and not on MCMC-sampling. In the rank consistency experiment, lowBMM is performing slightly better than EMM for d R andp (see Figure 9, and the right portion of Table 2). These results also show that BMM and MM are not even close to converging to the correct solution in the rank consistency scenario. It is also worth noting that lowBMM and BMM are Bayesian methods that provide uncertainty quantification, while EMM and MM are methods that result in point estimates.

Comparison with the Bayesian Mallows Model (BMM)
Among the performance comparisons of lowBMM with alternative methods, a special role is played by BMM: indeed BMM, as introduced in [38] and subsequently implemented in the BayesMallows R package [35], is the original full dimensional version of lowBMM, so we also would like to verify that lowBMM is consistent with the full model in a limiting situation. To this aim, we ran the full BMM as well as the variable selection lowBMM on the same simulated dataset with n = 20, N = 50, n * = 8. For a comparison of the marginal posterior distribution of ρ obtained on the reduced set of items with the lowBMM to the one obtained with BMM on the complete set of items, see Figure  S10 in the supplementary material, right and left panel, respectively. Both models show consistent results for the top items, with the full model being slightly more certain in its top ranks. The noisy items seem to be randomly distributed on the bottom in the full model, while in lowBMM they are assigned the bottom ranks. We then compared BMM and lowBMM on a larger simulated scenario (n = 1000, N = 50, n * = 50), when using both data generating processes: Figure 10 demonstrates that lowBMM outperforms BMM in terms of accuracy and stability in both scenarios. In the top-rank simulation experiment, lowBMM converges to the correct solution with fewer iterations, while BMM reaches the solution at a later stage. This shows the consistency of lowBMM, as we obtain the same solution as for the full model, however in the correct reduced dimension and with lower computational demands. For the rank-consistency simulation experiment, lowBMM clearly outperforms BMM: while lowBMM converges to the solution after 3000 iterations, BMM exhibits an oscillating behaviour, and it never manages to reach convergence. This last simulation example further shows the flexibility in the variable selection procedure under different data generating procedures, even in a high-dimensional setting with basically no possible (Bayesian) competitor.

Off-line estimation of α
Finally, we also aimed at testing the approach described in Section 2.4 for the off-line estimation of α, in a simulation setting with the same top-rank data generating process as previously used. We therefore generated datasets with varying dimensions (N, n), where n * = round(n/3) items were sampled from the Mallows model with a fixed α true , and the rest of the items were sampled from a uniform distribution. For each dataset in a given dimension, we ran the estimation procedure of Section 2.4 by generating datasets of equal dimension over a grid of α 0 values, and computed the mean distance between assessors in each dataset. Finally, we obtained the final estimateα n * by determining the intersection between the distance measure for the simulated datasets and the distance measure for the dataset generated from α true , in each dimension. Performance of the estimation procedure can be inspected in Table 3, where the method stability is demonstrated even for increasing n. The estimated valuesα n * are consistently slightly smaller than α true , which however does not affect the method performance, but simply makes the posterior distributions for ρ and A * slightly more peaked to absorb the unnecessary uncertainty. We remark that a smaller α than necessary is in theory better for estimation, as it allows for more model flexibility and for better space exploration in the MCMC.  Table 3: Results from simulation experiments described in Section 3.5, for testing the off-line estimation of α for varying (N, n) as described in Section 2.4. α true = 3 for all scenarios.

Application to RNA-seq data from ovarian cancer patients
In this section we describe the application of the lowBMM method on RNAseq data from TCGA ovarian cancer patients. The gene level expression data was available at the TCGA Data Portal (https://tcga-data.nci.nih.gov/tcga/) 3 , and the selected ovarian cancer patients were the same as previously analysed [15].
Serous ovarian tumors share similar genetic origin and features with basal-like breast cancers [6]: both are difficult-to-treat cancer types, and they are highly commonly characterized at the molecular level, especially for what concerns mutations (with TP53 being the most frequent example). Massive whole-genome data analyses are therefore key to understanding individual tumor biology for these cancer types. This is the main motivation for focusing on RNAseq data from serous ovarian tumors, as selecting basal-like breast cancer patients would require an additional subtyping preprocessing step that we would like to avoid. Additionally, since our method assumes complete data, preprocessing steps included keeping only the items (genes) with less than 50% missing values: all of these missing values were imputed using k-nearest neighbor averaging [37]. This resulted in a final dataset of ultra-high dimension with n = 15348 genes and N = 265 patients.
For running our lowBMM method, we set n * = 500, l = 20 and L = 1 (tuning parameters were set based on the conclusions from the simulation studies in Section 3.2, and from some trial-anderror). We also performed an off-line estimation of α on the experimental dataset, on a search grid α ∈ [10 −15 , 10 −5 , 1, 10, 50, 100] that was set as such to ensure adequate coverage over several orders of magnitude, as described in Section 2.4, resulting in α = 10. Due to the high-dimensionality of the data, we ran two chains, each with M = 5 × 10 6 MCMC iterations on a big memory HPC server. We discarded the first 5×10 4 iterations as burn-in in both chains before merging them for post-processing (see Section 2.3.1 for details on the post-processing of the results). As A * is a probabilistic selection, we will get slightly different posterior summaries from different runs of lowBMM. Moreover, A * is a discrete parameter defined in a space whose dimension grows exponentially with n, and in this particular case the space will be exceedingly large. This will in turn result in small variations in each run of the algorithm, as we are limited by computational capacities in how much we can realistically explore this space. Furthermore, this behaviour can also be affected by aspects related to the current problem setting: particularly as the n * value used in the analysis might be too low compared to n. We remark that the sensitivity study on the tuning parameters in Section 3.2 suggested that running lowBMM with a smaller n * compared to the "truth" might make it more difficult for the method to converge to the true ρ A * . Moreover, during the off-line estimation of α, we estimated the n * that best fits the experimental data to n * ≈ 2400. However, we kept the final n * = 500 value as it is more suited to allow results interpretability, and as it also allows the MCMC to perform a good exploration of the state space given the current computational capacity. The computations were performed on resources provided by UNINETT Sigma2, the National Infrastructure for High-Performance Computing and Data Storage in Norway.
To gain insight on the results provided by lowBMM on the experimental dataset, we computed a "top probability selection" as described in Section 2.3.1, to be used in a Gene Set Enrichment Analysis (GSEA) [32,36]. This analysis was done to verify whether the method could identify common genes related to cancer. To this aim, we computed various top-K sets, with a threshold for the lower items: genes that had probability less than 0.1 of being in the different top-K sets were dropped from the plot. This resulted in the violin plots shown in Figure 11, left panel: this plot indicates that K = 75 is a good choice for allowing many genes to have a quite large posterior probability of being top-ranked. We therefore further inspected the probability distribution of the top-75 probabilities across the genes, by drawing a histogram of the same probabilities ( Figure 11, right panel): this second plot allowed us to select a cut-off point of c = 0.77 for the top-75 probability set, thus resulting in a final selection of |Â * top | = 63 genes whereÂ * Figure 12 shows that the items in the refined selectionÂ * top are selected with very large posterior probability, even if they sometimes "compete" for some of the available ranks: this effect had already been observed in the context of the simulation studies in Section 3.2, and it might be due to n * being smaller than n * true . There is overall more variability in the rankings of the items as compared to the simulation studies, which is to some degree expected in experimental data with a higher degree of noise. Higher variability in the ranks was also typical for the simulation studies where n * guess << n * true (see Figure S5 in the supplementary material). The GSEA indicated that the lowBMM identified known genes involved in cell differentiation and cancer development. The top-ranked genes selected inÂ * top were enriched in pathways such as cell differentiation, morphogenesis and developmental processes (Table 4), suggesting a role in the deregulation of cellular identity essential for cancer development.   Ranking

Discussion
In this paper, we have developed a novel rank-based unsupervised variable selection method for high-dimensional data. The method extends the applicability of the Bayesian Mallows ranking model to unprecedented data dimensions and provides a much better modeling approach. Indeed, an important advantage of the variable selection procedure is the ability to distinguish relevant genes from the background ones, and to provide an estimate of the relative importance of the selected relevant genes. Furthermore, the variable selection procedure is able to work in highdimensional settings where the number of items is n >> 10 3 , with no competitor model for ranking data capable of scaling to such setting. The simulation studies described in Section 3 showed that the method performs well on datasets of varying sizes and under varying data generating procedures. Additionally, lowBMM is superior in terms of accuracy and computational time compared to existing methods in high dimensions, and has no competitor in ultra-high dimensions. We applied the proposed model to RNAseq data of ovarian cancer samples from TCGA (Section 4). Although we consider a cancer genomics application, the methodological contribution given in this paper is very general, and can thus be applied to any high-dimensional setting with heterogeneous rank data.
Throughout this work, we have assumed that the number of relevant items to be selected, n * , is known. This is reasonable in many practical problems where a specific n * can be tuned according to the research questions/demands or has to be decided according to computational capacity. Simulation studies showed that assuming n * guess << n * true might pose problems in high-dimensional settings, and therefore the joint estimation of this parameter as part of the existing hierarchical framework constitutes a very interesting direction for future research. Nonetheless, as mentioned in the simulation studies while commenting upon the marginal posterior distribution of ρ obtained when n * guess < n * true , our method can be seen as an alternative to the constrained median bucket order technique introduced in [13]. It would then be interesting to fully compare the respective estimates obtained via lowBMM and the median bucket order, in the light of the relationship between the choice of permutation space (strong vs weak rankings), distance (any distance vs the Kemeny only), and estimation procedure (Bayesian vs frequentist). This comparison is clearly out of the scopes of the present paper, but makes up a very interesting future research direction.
The current version of lowBMM works under several assumptions: no handling of missing data, no clustering, fixed scale parameter α, all aspects which we plan to implement in future versions of the method. The somewhat challenging procedure of choosing the tuning parameters l and L might be solved with an adaptive MCMC and would allow the algorithm to explore the space more thoroughly. We leave this task for future consideration.
To our knowledge, no rank-based unsupervised variable selection procedure that is capable of scaling to the common data dimensions in -omics applications currently exists.