GeneCompete: an integrative tool of a novel union algorithm with various ranking techniques for multiple gene expression data

Win-loss method

The win-loss method finds the ratio of the number of wins to the number of matches attended. Let n_ij be the number of games played between player i and player j, and w_ij be the number of times player i wins player j. The ranking of player i can be computed as: (2) $$r_w\left(i\right)=\sum _i{\displaystyle \frac{w_{ij}}{n_{ij}}}$$

For the intersection strategy, all genes participate in the same of games, with n_ij =k(| S_int |−1), ∀i,j. Then, the ranking can be calculated based on the total number of wins: $r_w\left(i\right)=\sum _i{w}_{ij}$. The maximum value is k(| S_int |−1) in the case of winning all players in all matches, and the minimum is 0 when losing every game. In the case of the union strategy, the term w_ij/n_ij only occurs when player i competes with player j. Therefore, the more datasets to which gene i is connected, the more opponents the gene has. Thus, achieving a larger winning percentage with more occurrence in datasets leads to higher gene scores. The advantage of this algorithm lies in its simple concept and low computational time.

Massey’s least squares method

Massey algorithm was originally proposed by Massey (1997) to rank football teams. Let $\tau =\sum _i{N}_i$ and X be the $\tau \times N$ matrix present the outcome of games, $$X_{ti}=\{\begin{array}{c}1ifteamiwin{t}^{th}game\\ -1ifteamilost{t}^{th}game\\ 0otherwise\end{array}$$

Each row in X_ti contains just two non-zero elements: 1 for the winner and −1 for the loser. In the Massey algorithm, explicit opponent indices are not used; instead, it relies on game outcomes (wins and losses) to compute team rankings. This method considers how teams perform against different opponents, ultimately assigning ratings or rankings. While the precise mathematical procedures can differ in various Massey algorithm implementations, the opponent's identity is typically inferred directly from the game results.

The Massey matrix is defined as M = X^TX. It can be expressed in terms of number of games as $M_{ij}=\{\begin{array}{c}{N}_{i}ifi=j\\ -n_{ij}ifi\ne j\end{array}$, where $N_i$ is the number of games played by player i. This definition is explained by Massey as ${\left(X^{T}X\right)}_{ij}={\mathit{x}}_i\cdot {\mathit{x}}_j$, where ${\mathit{x}}_i=\left\{1,2\dots ,N\right\}$ is the column vector of X. For the diagonal elements, let’s consider ${\mathit{x}}_i\cdot {\mathit{x}}_i$, in cases where the game is played, $x_i$ can be 1 or −1, resulting in the summation of the number of games played by player i. When i ≠ j, the term ${\mathit{x}}_i\cdot {\mathit{x}}_j$ can be non-zero only when there is a match between player i and j, which has exact values, i.e., 1 and −1 are multiplied together and summed for every game, resulting in $-n_{ij}$.

Let y be the vector of point differentials, where the t^th component of y is the point difference in the t^th game, and p = X^Ty. Since the Massey rating r_ms can be calculated from Xr_ms = y, then X^TXr_ms = X^TXy. Subsequently, the simple Massey linear equation is given as: (3) $$Mr{}_m{}_s=p$$

Notice that the last row of M is replaced by a vector of ones, and the last row of p is replaced by zeros row because M is s singular matrix, and Eq. (3) cannot be solved. Consequently, the addition of Massey scores of all players equals zero, $\sum _i{r}_{ms}\left(i\right)=0$.

Notably, the win-loss method simply counts the number of wins and losses for each team without considering the margin of victory or defeat. It's a straightforward way to assess performance based on win-loss records. On the other hand, the Massey algorithm takes into account the margin of victory or defeat. It does not just treat all wins and losses equally. Teams are ranked based on a more sophisticated assessment of their performance, which can provide a more accurate representation of team strengths. In the intersection strategy, both Ni and nij (number of games played and number of games won) are equal for all genes. The rating is primarily based on the win-loss record, which is similar to the win-loss method. It focuses on the number of games won and lost without considering the margin of victory or defeat. In contrast, the union strategy assesses performance based on the percentage of games won. This strategy gives more weight to how convincingly teams win games, considering the margin of victory. It can provide a more fine-grained evaluation of team strength, rewarding teams that not only win but also do so decisively. Key distinction lies in how the algorithms handle the margin of victory or defeat. The win-loss method and intersection strategy primarily focus on the number of wins and losses, whereas the Massey algorithm and union strategy consider the margin of victory, providing a more nuanced and accurate assessment of team performance.

Colley’s least squares method

Colley (2002) discovered a ranking model that applies Laplace’s rule of succession in a linear model. Let $W_i=\sum _j{w}_{ij}$ and $L_i=\sum _j{l}_{ij}$ be the number of wins and losses for team i. First, the Colley matrix has a high connection with Massey matrix, C = M + 2I, where I is the N × N identity matrix, or it can be defined as $C_{ij}=\{\begin{array}{c}{N}_i+2ifi=j\\ -n_{ij}ifi\ne j\end{array}$.

Let $b_i=1+{\displaystyle \frac{W_i-L_i}{2}}$ be the difference between the number of wins and losses, which is derived from the modified winning percentage ${\displaystyle \frac{W_i+1}{N_i+2}}$ to start the rating at 0.5. The derivation of b_i from the modified winning percentage is presented in Data S1. Then, the rating r_c is computed by solving the equation (4) $$Cr{}_c=b$$

Although, the Colley and Massey algorithms stem from different motivations (Devlin & Treloar, 2018), the two matrices are similar. Hence, they have led to nearly the same results since our application considers the score of matches as a win-loss record in Massey.

Keener’s method

Keener (1993) proposed a eigen-based ranking model by applying the Perron Frobenius eigenvector. The concept behind constructing the Keener matrix is to differentiate between a dense number of players with a win probability around 0.5. This model uses a non-linear skewing function, $h\left(x\right)=0.5+0.5\left(sgn\left(x-0.5\right)\sqrt{\left|2x-1\right|}\right)$. Next, the probability of winning, based on Laplace’s rule of succession, is represented as $a_{ij}={\displaystyle \frac{W_{ij}+1}{W_{ij}+W_{ji}+2}}$. Then, the Keener matrix is defined as K_ij = h(a_ij). Definitively, the Keener rating r_k is obtained by solving: (5) $$Kr{}_k=\lambda r_k$$where λ is the largest eigenvalue, and r_k is the corresponding eigenvector.

As mentioned, the intersection strategy considers the same number of games for all genes. The rating result is based on the winning percentage, similar to Colley, but it has the advantage of distinguishing near-zero probability, which can lead to more accurate results than Massey and Colley. For the union strategy, more candidate genes are considered; however, Keener requires more time to solve for eigenvectors and eigenvalues.

Elo’s ranking method

Elo (1978)’s system was first established for ranking chess players. Player rankings are based on their previous performances, changes in ratings occur through iterations. Elo’s rating for team i can be computed as: (6) $$r_E^{new}\left(i\right)=r_E^{old}\left(i\right)+f\left({\kappa }_{ij}-{\mu }_{ij}\right)$$where

${\kappa }_{ij}=\{\begin{array}{c}1ifplayeriwinplayerj\\ 0ifplayerilossplayerj\\ 0.5ifplayeriandplayerjaretie\end{array}$ represents the actual outcome of the game.

The constant f is set to 10, and ${\mu }_{ij}={\displaystyle \frac{1}{1+{10}^{\left[r_E^{old}\left(j\right)-r_E^{old}\left(i\right)\right]/400}}}$ is the expected probability of player i winning against player j, constructed using a logistic function. Elo’s method requires the initial ratings of each player as input; this work uses equal values for all players, with $r_E^{old}=1500$. When considering the first pair of players, the player who wins the game gains a higher rating, while the loser’s rating decreases. The process is iterative until the last pair of players is considered.

Markov method

The Markov chain can be used for ranking by considering a voting process, where the stronger alternative is voted for by the weaker alternative (Von Hilgers & Langville, 2006). Nowadays, the Markov chain has been developed in various forms; for example, the (1,∞) variant is proposed to reduce the sensitivity of the model (Vaziri, Yih & Morin, 2018). Generally, the original Markov method is constructed using the (0,1) voting matrix $V_{ij}=\{\begin{array}{c}1ifplayeriwinplayerj\\ 0otherwise\end{array}$.

Next, the transition probability matrix P is obtained by normalizing the voting matrix or dividing each element by its row summation. Then, the Markov rating vector r_mk is obtained by solving: (7) $$r{}_{m}k=Pr{}_m{}_k.$$

The Markov ranking method takes into account both the opponents and their level of strength. However, this method has displayed sensitivity to small changes in data and also requires a long computational time, making it more suitable for solving problems with small number of players.

PageRank method

PageRank was first proposed by Google’s founders, Larry Page and Sergey Brin, to rank web pages (Brin & Page, 1998). Unlike the previous methods, this model is constructed using a network. First, a directed graph G is constructed by considering each node as a player, with directed edge pointing from the losing player to the winning player. A higher number of in-degrees for a node indicates a stronger opponent. Let B_u be the set of neighboring nodes pointing to u, and |v| be the outgoing degree from node v. Then, the PageRank score r_p(u) of player u is defined as: (8) $$r_p\left(u\right)=\sum _{v\in B_u}{\displaystyle \frac{r_p\left(v\right)}{\left|v\right|}}$$

Alternatively, the power method is applied to quickly solve for the PageRank rating. Let A be the N×N adjacency matrix of the graph G. A is normalized by row summation to obtain Z, and G is a Google matrix computed as G = αZ + (1− α) E, where E is an entirely 1/N matrix, and α is a damping factor, usually set to 0.85. Thus, the PageRank score r_p(u) can be computed as: (9) $$r{}_p=r{}_{0}G{}^c$$where r₀ is the initial vector, which is set to 1/N if no initial vector is provided, and c is the number of iterations needed to reach convergence. To apply PageRank in the union strategy, genes in S_up are treated as nodes in the network for the up-regulated case, and S_down for down-regulated case. Genes that participate many games and frequently win against strong opponents tend to have high PageRank scores.

Bi-directional PageRank method

The improved model of PageRank is developed for sport ranking (Zhou et al., 2022). Bi-directional PageRank (BiPageRank) considers both win and loss whereas PageRank computes only the winning score. This work shows the outperforming results of BiPageRank when compared with PageRank using both synthetic data and application of four sports: soccer, basketball, ice hockey, and baseball. The BiPageRank can be computed as: (10) $$r{}_s=r{}_p-r{}_q$$where r_p is the PageRank score, and r_q is the backward propagation of PageRank score, which can be calculated as $r_q\left(u\right)=\sum _{v\in Q_u}{\displaystyle \frac{r_q\left(v\right)}{|v{|}^{in}}}$, where Qu is the in-neighbor of node u, and |v|ⁱⁿ is the in-degree of node v. PageRank r_p assigns a higher score to players who frequently win against strong teams. In contrast, r_q assigns a higher value to players who lose to low-rated teams. Thus, the BiPageRank score improves upon PageRank by considering both the wins and losses of the players.

Gene expression data of hypertrophic cardiomyopathy

Different experimental types, sample origins, and platforms used to collect data on HCM provide varying information. To ensure the reliability of results and confirm the importance of genes to the disease, it is crucial to include a larger number of samples in the model. In this study, we gathered nine datasets of HCM gene expression data collected from the Gene Expression Omnibus (GEO) database. The data comprise four microarray datasets: GSE36961, GSE32453, GSE68316 (Yang et al., 2015), and GSE1145 and five RNA-Seq datasets; GSE89714, GSE130036 (Liu et al., 2019), GSE160997 (Maron et al., 2021), GSE180313 (Ranjbarvaziri et al., 2021), and GSE141910. In total, there are 464 samples, consisting of 213 cases and 251 controls. The characteristics of the HCM gene expression data are provided in Table S1.

Microarray quality control project

The United States Food and Drug Administration (FDA) provides data of Microarray Quality Control (MAQC) and Sequencing Quality Control (SEQC). MAQC was first developed to evaluate agreement across microarray data and is provided in GSE5350 (Li et al., 2014; MAQC Consortium, 2006; Su et al., 2014; Wen et al., 2010). With the emergence of next-generation sequencing technologies, SEQC was introduced to access RNA-Seq performance, and it is available in in GSE56457 (MAQC Consortium, 2014), GSE47774 (Su et al., 2014), and GSE48016 (Munro et al., 2014; Wang et al., 2014). From the four types of samples provided by the United States Food and Drug Administration (FDA), we have selected two types of RNA samples: A (Universal Human Reference RNA) and B (Human Brain Reference RNA). We gathered nine datasets from GEO database to obtain gene expression data from 1442 samples, with 721 samples for each type, as provided in Table S2.

Prioritization techniques with up- and down-regulation genes

LOOCV involves leaving one dataset as the testing set while combining the others using ranking methods. Higher performance indicates a stronger relationship with the DEGs of the test set.

To compare the ranking performance with the original method, we consider two cases: up-regulation and down-regulation. The original or classical method for identifying DEGs based on applying specific criteria: a logFC > 1 and an adj.p.val. < 0.05 for up-regulation genes, and a logFC < −1 and an adj.p.val. < 0.05 for down-regulation genes. Subsequently, each gene’s count score is calculated by summing the total number of datasets meeting the criteria of logFC > 1 and logFC < −1, for each up or down cases. The average logFC (Avg_logFC) is then directly computed as the mean logFC across datasets. The ROC curve becomes visible after a single iteration of leave-one-out cross-validation. The ROC curves for union strategy in up-regulation and down-regulation can be found in Datas S2 and S3, while the results for intersections are presented in Datas S4 and S5. Fig. 3 shows that the original method exhibits the low predictive power for DEGs in both up-regulation and down-regulation cases. The count score method improves upon the original approach, suggesting that genes present in more datasets with a high absolute logFC are more likely to predict DEGs. Consequently, ranking scores that consider the number of datasets can be valuable for the prediction.

Notice that many cases demonstrate the union strategy yielding higher performance than the intersection approach across many ranking methods. In the case of the win-loss method, the union strategy, which considers the number of datasets a gene has joined, demonstrates better performance than the intersection strategy, with the same number of datasets. Massey’s approach applied the winning score as input, showing similar performance for both methods. Colley and Elo, utilizing a probability of win, can reduce the effectiveness of the union ranking strategy due to genes that appear in a greater number of datasets having a higher probability of score reduction, unlike the intersection strategy that maintains a consistent score regardless of dataset count. A modification of the Keener matrix did not improve rankings for this task and achieved a similar performance to the win-loss method. The Markov method shows high sensitivity in union ranking, especially in cases of up-regulation. This implies that minor data changes can lead to substantial ranking differences; for instance, a lower-ranked player defeating the highest-ranked player might result in a significant increase in the former’s ranking. Both PageRank and BiPageRank exhibit similar behaviors, though BiPageRank displays slightly lower performance. By the concept of PageRank, genes that wins other important genes tend to have a higher score. Moreover, we observe that genes that are presented in a higher number of datasets play a higher number of matches and have a higher chance to receive a PageRank score. Among PageRank, Win-loss, Keener, and BiPageRank, as illustrated in Fig. 3, most instances demonstrate that PageRank and BiPageRank exhibit superior performance in terms of both the AUC of the ROC curves and the AUPR of the precision-recall curves. The common thread shared by PageRank, Win-loss, Keener's Algorithm, and BiPageRank is their focus on ranking genes within their respective relationship networks. While PageRank and BiPageRank underscore the significance of connections and links to other nodes, Keener's Algorithm takes into account both local and global influences from others. Conversely, Win-loss simplifies ranking by relying on binary outcomes in competitive scenarios.

It is worth noting that in both the intersection and union strategies, genes can be categorized based on their positive and negative logFC values. Comparisons between the separated and non-separated versions can be found in Fig. S1. Interestingly, the outcomes appear more favorable for the separated approach. However, it is important to highlight that the separated approach yields only 23 candidate genes for up-regulation and 105 for down-regulation.

Winner genes with the best top ranking

The PageRank method method stands out as the best approach for intersection strategy. Tables S3 and S4 present the top 10 ranking genes for both up-regulated and down-regulated cases. The logFC values of genes in many datasets are lower than 1 in up-regulated case and greater than −1 in down-regulated case. This suggests the gene expression differences in HCM are not reaching a two-fold changes compared to normal patients.

Tables S5 and S6 display the top 10 ranking genes obtained by each method. Interestingly, the top genes identified by the win-loss, Keener, PageRank, and BiPageRank methods are quite similar. In Table 2, our method identifies the top-ranking genes using PageRank, which is supported by existing literature evidence. HCM is closely associated with dilated cardiomyopathy (DCM), a type of heart muscle disease that can lead to heart failure and life-threatening arrhythmia (Tobita et al., 2018).

Among the top 10 up-regulated genes, SLITRK4 has been identified as a promising biomarker for HCM gene tests (Zheng et al., 2021). It is commonly differentially expressed in the datasets related to the study of HCM patients, such as GSE130036 and GSE36961 (Cui et al., 2022). SFRP4 is known to be involved in cardiac development and various cardiovascular diseases (Zeng et al., 2019). It has been identified as a hub gene in the HCM key module (Ma et al., 2021b) and as an up-regulated DEGs in HCM (Ren et al., 2016). In addition, SFRP4 is associated with ischemic cardiomyopathy, a type of heart muscle disease (Alimadadi et al., 2020), and has been verified as a hub gene associated with heart failure (HF) (Zhou et al., 2020).

For CA3, an increase in its expression has been confirmed by immunohistochemistry as a myocardial protein (Coats et al., 2018). Moreover, CA3 expression levels were significantly higher in the plasma of heart failure patients than in control patients (Su et al., 2021). FRZB has been identified as a hub gene in the HCM key module (Ma et al., 2021b) and hub biomarkers for dilated cardiomyopathy (DCM) (Fang et al., 2022). In addition, FRZB has been recognized as a potential immune-related key genes involved in ischemic cardiomyopathy through random forest analysis and nomogram (Zheng et al., 2023).

MXRA5 has been identified as a key gene with prognostic value in left-sided HF (Zhou et al., 2020). It is extracellular-associated proteins included in the top 500 genes in the HF consensus signature (Ramirez Flores et al., 2021). SMOC2 has been defined as a real hub gene of HCM due to its high intramodular connectivity values (Jiang et al., 2021). The protein encoded by the differentially expressed methylated gene SMOC2 was found to be upregulated in Chagas disease cardiomyopathy (Shi et al., 2022). THBS4 is implicated in severe HCM and heart failure pathogenesis (Tsoutsman et al., 2013). It is also predicted to play a role in the development of DCM (Zhao et al., 2018). THBS4 expression has been associated with hypertrophic cardiac disease (Peisker et al., 2022). FNDC1 was among the 10 most up-regulated transcripts in patients undergoing repair of tetralogy of Fallot heart tissue, compared with right ventricle donor tissue (Brayson et al., 2020). Both FNDC1 and MXRA5 have been identified as novel extracellular matrix (ECM) biomarkers in calcified valves, making them potential targets in the development and progression of aortic stenosis (Bouchareb et al., 2021).

FMOD has been identified as upregulated DEGs in heart failure (Kolur et al., 2021). It is a type of fibromodulin that is upregulated in clinical and experimental heart failure (Andenæs et al., 2018). DIO2 is a direct transcriptional target of the FoxO1 protein, which is involved in relative hypertrophic growth of neonatal cardiomyocytes in vitro and in vivo (Ferdous et al., 2020). It has been reported that DIO2 is up-regulate in the hearts of DCM mice (Wang et al., 2010).

For the top 10 of down-regulated genes, FCN3 is a key dysfunctional gene. It was identified by studying the network of differentially expressed genes between HCM and healthy controls (Cui et al., 2022). Additionally, FCN3 is associated with the development of HF (Jiang, Zhang & Zhao, 2022). CORIN was identified as a downregulated mRNAs in the myocardial tissues of patients with HCM (Cao & Yuan, 2022). It was reported to be a cardiac protease that activates natriuretic peptides, the expression of which has been examined and studied in the activity of mouse and human failing hearts (Chen et al., 2010). Regarding HOPX, the relationship between HOPX gene variations and HCM was investigated. The results suggest that HOPX may cause pathogenesis or manifestation of HCM (Güleç et al., 2014). A study showed that HOPX expression is reduced and completely absent in severe heart failure (Trivedi et al., 2011). In addition, the HOPX gene plays an adjusted role in HCM pathogenesis through SRF-dependent genes (Alkanli & Ay, 2019). It was reported that the expression of MYH6 is dominant in human cardiac atria and plays roles in cardiac muscle contraction, including the composition of the cardiac muscle thick filament (Razmara & Garshasbi, 2018). Moreover, MYH6 mutations were evaluated in HCM phenotypes (Hsieh et al., 2022). The study showed that mutations in the MYH6 gene result in the abnormal development of cardiac muscle cells, which can lead to HCM.

SERPINA3 is significantly perturbed in heart failure proteins shared between two studies (Chen et al., 2022). It was reported to be downregulated in HCM compared to healthy controls (Chen et al., 2018). It is a common down-regulated DEGs in GSE130036 and GSE36961 (Cui et al., 2022). TUBA3E was identified as an HCM hub gene in the negative module (Jiang et al., 2021). It is also a common down-regulated DEGs in GSE130036 and GSE36961 (Cui et al., 2022). TUBA3E is included in a list of down-regulated genes expressed in patients with both HCM and DCM (Chaffin et al., 2022).

A study suggested that the potential function of CD163 macrophages is in supporting the homeostasis of cardiac tissue (Zhang et al., 2021). CD163 plays a key role in the pathogenesis of HCM (Zhao et al., 2016). It is a common down-regulated DEGs in GSE130036 and GSE36961 (Cui et al., 2022). A study reported that SMTNL2 is a down-regulated HF gene (Kolur et al., 2021). Regarding CCL2, it was reported that the CCL2-CCR2 signaling pathways are associated with the development and progression of cardiovascular disease (Zhang et al., 2022). RARRES1 expression was observed to be absent in the HCM samples in many of the fibroblast populations (Larson et al., 2020). It is one of the top three DCM down-regulated genes (Ma et al., 2021a).

To confirm the biological relevance of each ranking, we performed a gene set enrichment analysis. The results of up-regulation and down-regulation are shown in Tables S7 and S8, respectively. For up-regulation, genes in each method are involved in similar pathways, such as the extracellular region (GO:0005576), extracellular region part (GO:0044421), extracellular matrix (GO:0031012), extracellular space (GO:0005615), proteinaceous extracellular matrix (GO:0005578) and neurogenesis (GO:0022008). For down-regulation, genes were enriched in immune response (GO:0006955), extracellular region (GO:0005576), and defense response (GO:0006952).