The microarray data and clinical information for Control, Healthy obese, NAFLD, NASH, and HCC samples were obtained from the NCBI Gene Expression Omnibus (GEO) (Barrett et al., 2013) and Array Express (Parkinson et al., 2007) database. This dataset comprised a total of 132 samples with GEO accession numbers: GSE48452 with platform ID - GPL11532 (Affymetrix Human Gene 1.1 ST Array); GSE25097 with platform ID - GPL10687 (Rosetta/Merck Human RSTA Affymetrix 1.0 microarray). Details of datasets and respective sample information were provided as Supplementary Material S1. Following the merging of data, the R package “imputeTS” (version 3.3) (Moritz and Bartz-Beielstein, 2017) was utilized to impute missing values. This imputation process was only carried out on less than 5% of randomly missing values. Afterward, the “Limma” (version 3.57.3) package (Ritchie et al., 2015) was utilized to remove batch effects specific to each study, and data normalization was performed utilizing the robust multichip averaging (RMA) method (Bioconductor, 2023). Figure 1 presents a schematic diagram representing the methodology.

The entire dataset, encompassing samples from various disease stages ranging from NAFLD to HCC, is utilized to identify a suitable machine-learning technique with the primary objective of effectively categorizing the data into distinct groups based on disease stages. Each of the following machine learning (ML) models is applied individually to the entire pre-processed dataset using 10-fold cross-validation: DISCR (Discriminant Analysis), NB (Naive Bayes), RF (Random Forest), DT (Decision Tree), KNN (K-Nearest Neighbors), SVM (Support Vector Machine), and ANN (Artificial Neural Network). This procedure guarantees a thorough assessment and reduces the risk of overfitting by training and verifying the models on distinct subsets of the data. Through this screening, the ML technique with the highest accuracy values and other pertinent performance indicators is chosen for further investigation.

Twelve different feature selection methods are applied to the entire dataset. Each of these methods individually selects subsets of genes interactively, comprising varying quantities, such as the top 20, 30, 40, 50, and beyond. Subsequently, each subset generated by the different methods is utilized in the chosen machine learning technique, employing 10-fold cross-validation to calculate the classification accuracy. This guarantees that the performance of the chosen characteristics is comprehensively assessed. Furthermore, each feature selection method yields a subset of genes based on the highest level of accuracy attained. This technique enables the identification of the most pertinent genes that make a major contribution to classification, ensuring a strong and dependable selection of features. Among the feature selection techniques used, the ReliefF algorithm evaluates each feature’s significance based on its ability to distinguish between similar cases. It iteratively selects a sample, comparing it with the nearest samples from both the same and different classes. Features that effectively differentiate between classes receive higher weights, while those that distinguish within the same class receive lower weights (Urbanowicz et al., 2018). The process is repeated across multiple instances to reliably estimate feature importance, making ReliefF effective for noisy, multi-class datasets. The basic equation for updating the feature weights in ReliefF is: W A = W A − 1 m ∑ i = 1 k diff A , R , Hit i + 1 m ∑ C ≠ class R P C 1 − P class R ∑ i = 1 k diff A , R , Miss i C

Where,

W [A] This represents the updated weight or importance score of feature AAA.

W A on the right-hand side: This term indicates the current weight of feature A.

1 m ∑ i = 1 k diff A , R , Hit i This part calculates the average difference in feature A values between the Randomly decided on instance R and its K nearest neighbors that belong to the same class as R (denoted as Hit i ).

Cramér’s V is a filter-based technique that measures the correlation between two nominal variables. It involves creating a contingency table to compute the Chi-Squared statistic, which quantifies the strength of the association between the variables (Kearney, 2017). To account for the bias, the Chi-Squared estimate is reformulated as a Phi-Squared value. The degrees of freedom are subtracted from both the row and column counts to obtain the denominator of the Cramér’s V formula. Taking the square root of the quotient resulting from dividing the corrected Phi-Squared value by the denominator yields Cramér’s V in the range from 0 for no association up to 1 for the perfect association. This measure helps determine categorical data as it offers an idea of the strength with which variables are correlated without the assumption of any linear relation.

Chi-Squared Statistic χ 2 : χ 2 = ∑ O i − E i 2 E i where O i is the observed frequency and E i is the expected frequency.

Phi-Squared ( ϕ 2 ): ϕ 2 = χ 2 n where n is the total number of observations.

Corrected Phi-Squared ϕ corr 2 : ϕ corr 2 = max 0 , ϕ 2 − k − 1 r − 1 n − 1 where r and k are the number of rows and columns in the contingency table.

Corrected Row and Column Counts r corr = r − r − 1 2 n − 1 k corr = k − k − 1 2 n − 1

Denominator: denominator = min k corr − 1 , r corr − 1

Cramér’s V: V = ϕ corr 2 denominator

Pearson’s correlation coefficient (r) quantifies the strength and direction of a linear relationship between two variables (Nasir et al., 2020), ranging from −1 to +1. A value of +1 indicates a perfect positive correlation, where both variables increase together; −1 signifies a perfect negative correlation, where one variable increases as the other decreases; and 0 denotes no linear correlation.

Pearson’s correlation coefficient ( r ) is defined as: r = n ∑ x y − ∑ x ∑ y n ∑ x 2 − ∑ x 2 n ∑ y 2 − ∑ y 2

Where:

Kendall’s Tau is a non-parametric measure that quantifies the ordinal correlation between two variables (Valencia et al., 2019), focusing on the direction and magnitude of their association. While Pearson’s correlation coefficient measures only linear relationships, Kendall’s Tau is a measure that is very useful in cases where observations do not meet the assumptions of normality or linearity. It counts concordant and discordant pairs of observations for its calculation. Values range from −1 (perfect inverse correlation) up to +1 for perfect direct correlation. A value of 0 indicates no association.

Kendall’s Tau τ can be defined as: τ = Number   of   concordant   pairs − Number   of   discordant   pairs n 2

Where:

Another method, Spearman’s rank correlation coefficient (ρ) is a statistical measure that quantifies the degree and direction of the relationship between two variables that have been ranked (Schober and Schwarte, 2018). It measures how well the relationship between two variables can be represented as a monotonic function. This means that the Spearman rank correlation varies from Pearson’s correlation, as it considers the variables’ ranks while trying to figure out their relationship and not the variable’s actual values. Hence, it is applicable in assessing ordinal data that deviate from the assumptions of linearity and normalcy.

Spearman’s rank correlation coefficient ρ is defined as: ρ = 1 − 6 ∑ d i 2 n n 2 − 1

Where:

On the other hand, Conditional Infomax Feature Extraction (CIFE) selects features by maximizing a score J(S) (Lin and Tang, 2006). Features are incrementally added to the set S based on their significance and redundancy relative to already chosen features. In this way, the algorithm continues till it gets some pre-defined features or threshold scores. CIFE intends to offer an informative feature set that is concise in length, by managing relevance and redundancy towards better facilitation of subsequent ML models.

The CIFE can be expressed as J S = ∑ X i ∈ S I X i ; Y − β ∑ X i , X j ∈ S I X i ; X j ∣ Y

Where:

Among mutual information-based feature selection methods, the Joint Mutual Information (JMI) criterion is used to identify the most relevant features for a specific task. JMI selects features based on their high mutual information with the target variable while considering the combined information among selected features (Bennasar et al., 2015). This way of approaching features ensures that they are individually relevant and also collectively informative to reduce redundancy while enhancing the predictive capability of the model. JMI helps build an efficient and effective feature set, resulting in balancing relevance with redundancy.

The JMI criterion can be defined as: J X i ; Y , S = I X i ; Y + ∑ X j ∈ S I X i ; X j ∣ Y

Where:

The goal is to maximize J X i ; Y , S to select features that contribute the most information about the target variable while considering the redundancy with already selected features.

Another mutual information-based method, The Maximum Relevance (MIM) criterion is a technique for selecting features that have the most mutual information with the target variable (Che et al., 2017). The main objective of MIM is to identify those attributes with the highest relevance to predict the outcome, such that it improves the performance of the machine learning models. This technique lessens the dimensionality of the data without losing the most informative attributes. MIM especially comes in handy in applications where the dataset comprises several features, and the identification of the most relevant ones can greatly enhance model efficiency and accuracy. The MIM criterion can be mathematically defined as: M I M X i , Y = I X i ; Y

Where:

Mutual information I X ;   Y measures the amount of information obtained about one variable through another variable and is defined as: I X i ; Y = ∑ x i ∈ X i ∑ y ∈ Y p x i , y log p x i , y p x i p y

where:

Normalized Mutual Information Feature Selection (NMIFS) is a feature selection technique that aims to find and rank features based on their normalized mutual information in relation to the target variable (Estévez et al., 2009). Mutual information quantifies interdependence between variables, showing how much one can inform about another. By standardizing this metric, NMIFS accounts for variable scales and distributions, helping to identify the most informative features for predicting the target variable and enhancing the effectiveness and precision of machine learning models. The NMIFS score for a feature X i with respect to the target variable Y can be defined as: NMIFS X i , Y = I X i ; Y H X i + H Y

Where:

Lasso, or Least Absolute Shrinkage and Selection Operator, is a regularization technique in linear regression that enhances feature selection and controls model complexity (Muthukrishnan and Rohini, 2017). By adding a penalty based on the absolute values of coefficients, Lasso reduces the coefficients of less important features to zero, promoting simpler models and avoiding overfitting. It is widely used in data science and machine learning to streamline models and focus on the most relevant features. The objective function minimized by Lasso is given by: min β 1 2 N ∑ i = 1 N y i − ∑ j = 1 p x i j β j 2 + α ∑ j = 1 p β j

Where:

Ridge regression is a regularization method employed in linear regression models to mitigate overfitting and enhance generalization (Paul and Drineas, 2016). Ridge regression employs a L 2 penalty on the regression coefficients, in contrast to Lasso regression which utilizes a L 1 penalty. The penalty term is determined by a regularization parameter α , which determines the trade-off between accurately fitting the data and punishing big coefficients. min β 1 2 N ∑ i = 1 N y i − ∑ j = 1 p x i j β j 2 + α ∑ j = 1 p β j 2

N is the number of samples.

p is the number of features.

y i is the target variable for sample i .

x i j is the value of feature j e for sample i .

β j is the coefficient of feature j .

α is the regularization parameter that controls the strength of the penalty term ∑ j = 1 p β j 2

Gradient Boosting is an advanced ensemble learning method that iteratively combines decision trees. Unlike Random Forests, which build trees independently, Gradient Boosting constructs trees sequentially, with each tree correcting errors made by its predecessors (Otchere et al., 2022). It uses feature importance to identify significant features, revealing how much each feature contributes to accurate predictions. High feature importances allow for more contribution toward the overall impact of model performance while allowing some form of implicit feature selection, especially with complex datasets, making the model perform better and be more interpretable. The method includes several steps, such as -

Initialization of the model with a Constant:

Initialization of the model with a constant value, typically the mean of the target values for a regression problem. This can be represented as: F 0 x = arg min c ∑ i = 1 n L y i , c

Where L y , c is the loss function, y i are the true target values, and c is a constant.

Iterative tree building:

For m = 1 to M (where M is the number of trees to be built):

a. Computing the Pseudo-Residuals: r i m = − ∂ L y i , F m − 1 x i ∂ F m − 1 x i

Combining the Trees:

The final model F x after M iterations is: F x = F 0 x + ∑ m = 1 M η h m x

Next, Feature selection in Gradient Boosting is typically achieved by examining feature importances derived from the model. Each feature’s importance is calculated based on its contribution to reducing the model’s prediction error. Feature importance can be computed as follows:

Mathematically, the importance of the feature j can be represented as: Importance j = ∑ m = 1 M ∑ t = 1 T m Δ L j , m t

Where, T m is the number of nodes in the tree, and is the reduction in the loss function due to the split on feature j at node t in tree m .

The genes obtained from different feature selection approaches are used as subsets. These subsets were then used to create an ensemble by categorizing the genes by assigning values 12 to 1. A gene is classified as 12 if it appears in all twelve subsets, 11 if it is present in eleven of the twelve subsets, and so on. Subsequently, genes that fell into at least 6 subsets were selected for Cox regression-based survival analysis. This decision guarantees that a minimum of 50% of the feature selection methods will support the inclusion of these genes, hence strengthening the reliability subsequent to analysis.

The ensemble approach yielded a consolidated list of genes, then annotated using DAVID tool version 6.8 (https://david.ncifcrf.gov/) (Sherman et al., 2022). The criteria for conducting Gene Ontology (GO) and KEGG pathway enrichment analysis were set as a p-value below 0.05 and a false discovery rate (FDR) below 0.05. The GO enrichment analysis was utilized to ascertain the biological activities of these genes. GO ontologies are divided into three categories: molecular function (MF), cellular component (CC), and biological process (BP). The KEGG pathway enrichment analysis identified metabolic pathways that showed a significant enrichment of genes, as compared to the total genome background. The SRPlot online toolkit (http://www.bioinformatics.com.cn/srplot) (Tang et al., 2023) was used to display the findings of the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment studies.

Survival analysis was performed using the TCGA-LIHC (Liver Hepatocellular Carcinoma) dataset to investigate the association between the expression level of the genes and overall survival. First, clinical data and gene expression data were retrieved from the TCGA-LIHC cohort using the “TCGAbiolinks” package (Colaprico et al., 2016). The gene expression data is preprocessed, and a Variance Stabilizing Transformation (VST) is applied using DESeq2 (Love et al., 2014). The expression levels of the genes were extracted, and a median value was calculated to stratify the samples into “HIGH” and “LOW” expression groups. Next, the clinical data is merged with the gene expression data, and a Cox proportional hazards regression model is fitted to estimate the hazard ratio (HR) and its 95% confidence interval (CI) using the survival package (Therneau, 2021). Finally, a Kaplan-Meier survival curve is plotted using the “survminer” package (Pawar et al., 2022).

The biomarkers identified through survival analysis were used as input in two separate databases i.e., DGIdb (The Drug Gene Interaction Database) (Cannon et al., 2024) and Drug repurposing encyclopedia (Li et al., 2023). Additionally, clue.io COMMAND app web tool (Xie et al., 2022) was utilized to screen potential drug candidates. The DGIdb database uses a combination of expert curation and text-mining approaches to mine drug-gene interactions mined from DrugBank, PharmGKB, ChEMBL, Drug Target Commons, and others. On the other hand, The Drug Repurposing Encyclopedia utilizes the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/), and consensus drug profiles from DREIMT (http://www.dreimt.org/), which are derived from the Connectivity Map (CMap) LINCS gene expression dataset (https://clue.io/). The command app also utilizes connectivity map (CMap) analysis to screen possible drugs.

After data preparation, seven well-known machine learning techniques such as DISCR, NB, RF, DT, KNN, SVM, and ANN, were tested on the entire dataset. Figure 2 illustrates the accuracy scores of these machine learning techniques, each conducted independently on the whole preprocessed dataset using 10-fold cross-validation. Based on the accuracy values, DISCR outperformed other machine learning techniques, achieving the highest accuracy score of 0.90.

After selecting DISCR as the optimal classification method, twelve feature selection models were applied to identify the most significant gene subsets. These models included ReliefF, Cramér’s V, Kendall’s Tau, Pearson’s correlation, Spearman’s correlation coefficient, CIFE, JMI, MIM, NMIFS, LASSO, RIDGE, and Gradient Boosting. Iteratively, each of the twelve methods separately selects different subsets of genes like top 30, 40, 50, and so on. DISCR with 10-fold cross-validation is used for each subset of the different methods separately to check the classification accuracy. The accuracy details for different subsets of the genes, such as 30, 40, 50 to 150 genes, are reported in Table 1. Out of all the filter-based strategies, Relief F achieved the greatest accuracy of 0.70297 using a set of 30 characteristics. Similarly, the accuracy reached by Cramér’s V, Kendall Tau, Pearson’s correlation, and Spearman’s correlation coefficient was 0.623762, 0.683168, 0.722772, and 0.673267, respectively. These accuracies were obtained using 40, 20, 20, and 10 features. Among the approaches that use mutual information, CIFE, JMI, MIM, and NMIFS achieved accuracies of 0.643564, 0.752475, 0.475248, and 0.712871 using 40, 30, 70, and 20 features, respectively. Within the set of embedded methods, the LASSO algorithm achieved the highest accuracy score of 1 while using 40 features. Ridge and elastic-net algorithms achieved an accuracy of 0.752475 and 0.712871, respectively, using a total of 30 features. Figure 3 illustrates the performance matrices of different feature selection models (A-D) (Performance metrics for different feature selection models (E-L) were provided as Supplementary Material S2).

The subsets of genes obtained from different feature selection approaches are subsequently merged to create an ensemble set of genes. These genes are then categorized from 12 to 1 using a consensus approach with those scoring at least 6 being considered as most important features. Through this ensemble approach a set of ten genes were identified. These genes include C8B (Complement C8 Beta Chain), APOF (Apolipoprotein F), FBXL3 (F-Box and Leucine Rich Repeat Protein 3), ABAT (4-Aminobutyrate Aminotransferase), ZFP1 (ZFP1 Zinc Finger Protein), MBTPS1(Membrane Bound Transcription Factor Peptidase, Site 1), CENPV (Centromere Protein V), METTL23 (Methyltransferase 23, Arginine), RPL9 (Ribosomal Protein L9), and ABCB11 (ATP Binding Cassette Subfamily B Member 11).

Kaplan-Meier survival analysis was used to evaluate the predicted survival probability over time in liver cancer patients, comparing normal and HCC samples based on the expression levels of key genes identified through an ensemble approach. Out of the ten identified genes, eight genes were found to have the worst overall survival rate for HCC. The identified genes were ABAT (HR: 1.69), C8B (HR: 1.69), FBXL3 (HR: 1.43), ZFP1 (HR: 1.35), ABCB11(HR: 1.33), MBTPS1(HR: 1.21), CENPV (HR: 1.2) and APOF (HR: 1.18). In the case of ABAT gene, the average hazard ratio stands at 1.69. This suggests that individuals with high gene expression face a 69% elevated risk of death compared to those with low gene expression, regardless of the circumstance as depicted in Figure 4. Similarly, for C8B and FBXL3, high expression corresponds to a 69% and 43% increased risk of death, respectively, in comparison to scenarios with low expression of these genes. HR value greater than 1 indicates a higher hazard rate in the reference group compared to the comparison group, suggesting that genes with higher expression in the comparison group are associated with an increased risk of death. Therefore, these genes could serve as predictive markers for poor prognosis. Survival analysis results of the rest of the genes are provided in Supplementary Material S3.

The gene ontology analysis using the ten identified genes from feature selection and based on the selected identifier - “OFFICIAL GENE SYMBOL,” and p-value < 0.05 and FDR <0.05 cut-offs yielded significant GO terms, depicted in Figure 5, for the marker genes. The enriched GO terms related to biological processes found to be negative regulation of gamma-aminobutyric acid secretion (GO:0014053), histone H3-R17 methylation (GO:0034971), gamma-aminobutyric acid biosynthetic process (GO:0009449), positive regulation of catecholamine metabolic process (GO:0045915), regulation of cholesterol biosynthetic process (GO:0045540), regulation of sterol biosynthetic process (GO:0106118), etc. The results suggest the involvement of metabolism and epigenetic processes, indicating dynamic changes in the disease progression. The top GO terms associated with cellular components terms identified as membrane attack complex (GO:0005579), low-density lipoprotein particle (GO:0034362), pore complex (GO:0046930), plasma lipoprotein particle (GO:0034358), etc; suggesting the involvement of membrane structure changes, lipid transport, etc, that are key characteristics of changes due to lipid accumulation and oxidative damage in hepatic cells in cellular dysfunction. The GO terms associated with molecular functions are histone-arginine N-methyltransferase activity (GO:0008469), carbon-sulfur lyase activity (GO:0016846), cholesterol binding (GO:0015485), pyridoxal phosphate binding (GO:0030170), etc., indicating lipid interaction and metabolic processes as important factors in this disease progression.

The genes were also subjected to KEGG pathway analysis to determine their association with their corresponding biological pathways. A total of 13 pathways-related KEGG terms were obtained from the database. The criteria of p-value less than 0.05 and FDR less than 0.05 were used for analysis. The top enriched pathways and the respective gene counts are depicted in Figure 6. The KEGG analysis results showed that the genes are enriched in Butanoate metabolism (hsa00650), Propanoate metabolism (hsa00640), Alanine aspartate and glutamate metabolism (hsa00250), Complement and coagulation cascades (hsa04610), Protein processing in the endoplasmic reticulum (hsa04141), etc.

The integrated approach of drug repurposing, which utilized connectivity map analysis and drug-gene interaction text-mining, etc, facilitated the creation of a library of potential repurposed drugs. This resulted in the identification of 81 potential candidate drugs. All these candidates have been utilized for molecular docking to analyze further interactions (details of the drugs are given in Supplementary Material S4). Among these 81 drugs, DGIdb (The Drug Gene Interaction Database) based screening identified a total of 19 drugs. The drugs with the highest interaction scores (IS) were Vigabatrin (IS- 11.79601078), Divalproex sodium (IS- 8.847008084), and Pyruvic Acid (IS-5.89800539). Screening using the Drug Repurposing Encyclopedia (DRE) database identified a total of 55 drugs, with Scopolamine (Enrichment score (ES) - 0.973), Amiloride (ES - 0.971), Damnacanthal (ES - 0.967), and Esculin (ES - 0.962) showing the highest enrichment scores. In addition, the search conducted by the Command app revealed a total of 6 drugs, including Glibenclamide, Phenelzine, etc.

Docking results showed good interaction between ZFP1, C8B, MBPTS1, CENPV, ABCB11, and Diosmin with docking scores of −11.6821 kcal/mol, −11.134 kcal/mol, 10.4712 kcal/mol, −10.2391 kcal/mol, and −9.85606 kcal/mol respectively (Figures 7A–E). There are three hydrogen bonds between ZFP1 and Diosmin at THR112, ASN116 and GLU133 positions (Figure 7A) and eight hydrogen bonds between C8B and Diosmin at positions–ARG82, GLN65, CYS79, THR423, ASP424, LEU259, TYR166 and TYR141 (Figure 7B). The hydrogen bond interactions between MBPTS1 and Diosmin found at VAL54, TRP556, MET353, ARG386, and ASN515 (Figure 7C) and the interactions between CENPV and Diosmin included Hydrogen bonds at HIS210, ARG225, SER254 positions and PI-PI interactions at TYR561 and TRP256 positions (Figure 7D). The details about the interactions between top protein-ligand complexes are provided in Supplementary Material S5. ABAT and Esculin with a docking score of −4.30418 kcal/mol (Figure 7H). APOF and Lapatinib with a docking score of −6.90333 kcal/mol (Figure 7G), indicating a strong interaction between the protein and ligand. A satisfactory interaction has also been observed between FBX23 and Phenelzine with a docking score of −6.30952 kcal/mol (Figure 7F). The MM-GBSA analysis supported the docking results, showing strong binding affinities. CENPV-Diosmin had a binding free energy of −100.71 kcal/mol, with Van der Waals and electrostatic contributions of −49.50 and −46.93 kcal/mol, respectively. Similarly, C8B-Diosmin exhibited a binding free energy of −74.58 kcal/mol, with Van der Waals and electrostatic contributions of −37.40 and −52.24 kcal/mol (The detailed result of MM-GBSA analysis of the top ligand-protein complexes were provided as Supplementary Material S6). These results highlight a significant binding affinity driven by a balanced interplay of electrostatic, van der Waals, and lipophilic interactions.