Comparison of Normalization Methods for Selecting a Non-Small Cell Lung Cancer Marker Panel from Circulating miRNA RT-qPCR Data

1. Introduction

MiRNAs are small non-coding RNAs that play a key role in regulating gene expression and, consequently, biological processes. The stability of miRNAs in biological fluids and the correlation of extracellular miRNA expression levels in these fluids with pathological conditions make them promising biomarkers for non-invasive or minimally invasive diagnostics of various diseases, including oncological, cardiovascular, neurological, etc [1,2,3]. Among the drawbacks of using miRNAs as biomarkers are the variability of their expression depending on sample type and physiological state, as well as the use of indirect methods to assess their concentrations (reverse transcription polymerase chain reaction, RT-qPCR), which implies the need for normalization of individual miRNA expression data. Inadequate normalization can lead to incorrect assessment of expression levels and false conclusions [4]. Thus, the reliability of data regarding the content of individual miRNAs depends on the correctness of the normalization applied, and different normalization strategies can significantly affect the analysis outcomes of miRNA expression and the identification of meaningful biological differences [5,6,7]. Several normalization methods are widely used to assess individual miRNA levels, namely: normalization to stably expressed RNA molecules (e.g., miR-16, U6, miR-191), global normalization to the mean level of all miRNAs studied in the experiment, pairwise normalization of individual miRNAs, as well as TMM (trimmed mean of M-values) and VSN (variance stabilization normalization) methods [4,5,8,9]. Despite the abundance of methods, there is no universal "best" normalization technique for assessment of miRNA levels; the choice depends on data type, measurement platform, and research objective [4,5,10]. Each normalization method has advantages and disadvantages. For example, normalization to “housekeeping” miRNAs like miR-16 is often used as a simple internal control, but the stability of miR-16 expression varies across tissues, pathological and physiological states. Variants of global normalization (to the overall mean or similar) and functional normalization (to functionally important miRNAs) consider more complex biological relationships and can provide information on the contribution of studied miRNAs to the investigated phenomenon but require larger datasets and more complex computations. Quantile normalization is a method aimed at making signal intensity distributions identical across all samples. For this, values in each sample are sorted, averaged by rank, and returned to their original order, ensuring the same empirical distribution across all data. The method is widely used in microarray data analysis and similar data arrays obtained by other methods and effectively removes technical variation without considering a function describing data distribution [5]. The TMM method is based on calculating the trimmed mean of log-fold changes between pairs of samples, excluding extreme values and accounting for technical biases such as total expression level [11]. TMM normalization effectively balances differences in sequencing library depth and is often used to analyze RNA expression profiles after massively parallel sequencing. Another normalization method, VSN, transforms data so that variance minimally depends on the mean, thereby stabilizing expression variability and improving data comparability [12]. VSN uses models based on assumptions about noise and signal intensity distribution, increasing the reliability of subsequent analysis [5]. However, despite their advantages, quantile normalization, TMM, and VSN are generally not applied to normalize RT-qPCR data, including miRNA analysis [5]. The cycle threshold (Ct) value measured by RT-qPCR reflects the expression level of individual miRNAs rather than the distribution of signal intensities across the entire panel simultaneously, as in microarrays or RNA-Seq. Since quantile normalization aims to equalize value distributions across the entire sample set, its application to Ct values, which reflect individual rather than global characteristics, is inappropriate. Similarly, the TMM method, designed to align library size and remove extreme values in RNA-Seq data, loses relevance for RT-qPCR, where the number of miRNAs studied is small. VSN, intended to stabilize variance in large continuous datasets, is ineffective for RT-qPCR Ct values as its statistical assumptions do not fit the nature of discrete RT-qPCR data. As a result, RT-qPCR data normalization commonly uses stable reference genes (e.g., miR-16 or snRNA U6), the mean Ct of several stable miRNAs, or the ΔCt method to correct technical variation without distorting biological data [4].

The aim of this study is a comparative analysis of normalization methods for circulating miRNA expression data obtained by quantitative reverse transcription PCR to identify diagnostic miRNA markers, select the best ones, and form diagnostic panels based on them. Accordingly, we focused on comparing several normalization techniques: pairwise normalization, Tres- and Quadro-normalization (to the mean of pairs and triplets of miRNAs), normalization to the mean level of studied miRNAs, and specific functional normalization considering the division of miRNAs into pools of tumor suppressor and oncogenic miRNAs. The choice of these approaches is based on previously reported data regarding the efficiency of different miRNA expression normalization methods for solving the stated task and gaining insight into the biological role of miRNA markers. In this study, normalization methods were evaluated using model quality metrics, computational resources, and stability of feature selection on RT-qPCR miRNA level data from blood extracellular vesicles (EVs) of patients with non-small cell lung cancer (NSCLC) and healthy donors (HD).

It was shown that pairwise normalization, Tres-, and Quadro-normalizations are optimal for RT-qPCR data normalization aimed at discovering potential lung cancer markers and marker panels, which may improve NSCLC diagnostic accuracy using circulating miRNAs analyzed by RT-qPCR.

2. Materials and Methods

2.1. Samples of Biological Fluids and Their Analysis

This study utilized plasma EVs samples from 27 HD and 19 NSCLC patients, obtained from the National Medical Research Center named after Academician Meshalkin and the Federal State Budgetary Institution "Research Institute of Pulmonology" of the Federal Medical-Biological Agency of Russia.

All patients and donors provided informed voluntary consents to participate in the research. The study was approved by the Ethics Committees of the Federal State Budgetary Institution Institute of Chemical Biology and Fundamental Medicine of the Siberian Branch of the Russian Academy of Sciences (protocol No. 10 dated December 22, 2008) and the Federal State Budgetary Institution "Research Institute of Pulmonology" (protocol LEC No. 04-24 dated June 24, 2024). The donors' ages and disease stages according to the TNM classification are provided in Table 1.

Blood processing and isolation of EVs with subsequent miRNA extraction were performed as described in [13]. Real-time reverse transcription PCR (RT-qPCR) to assess expression levels of 17 microRNAs (primer and probe sequences are listed in Table 2) was conducted as described in [13] in two setups (miR-19b, -374a, -324, -22, -222, -133b, -144, -425 and miR-205, -660, -30e, -125b, -92a, -378a, -375, -27b, -31).

*Conditionally, since these normalization types include all investigated miRNAs.

As an internal spike-in control, cel-miR-39 was added to each sample before extraction, and its amount was subsequently quantified by real-time RT-qPCR (primer and probe sequences are given in the table) and used as a criterion for sample inclusion/exclusion in the analysis.

2.2. Methods of miRNA Expression Normalization and Statistical Analysis

The following normalization methods were used for miRNA expression analysis:

• Pairwise normalization (PN): For each unique miRNA pair, the difference between their Ct values was calculated. This accounts for relative expression of two miRNAs, reducing the influence of common technical variability.
• Tres normalization (TN): ΔCt values were calculated for each triplet of microRNAs by subtracting the mean Ct of the other two miRNAs in the triplet from the Ct of each miRNA. This method extends the concept of pairwise normalization by considering relationships among three markers.
• Quadro normalization (QN): ΔCt values were calculated for each quadruplet of miRNAs by subtracting the mean Ct of the other three miRNAs in the quadruplet from the Ct of each miRNA. This further extends pairwise normalization by considering relationships within groups of four markers.
• Normalization to the arithmetic mean (overall mean, OM): Ct values for each miRNA were normalized by subtracting the mean Ct of all miRNAs in the panel. This approach standardizes variability of overall expression in the sample.
• Normalization to the exclusive mean (EM): For each miRNA, Ct was normalized by subtracting the mean Ct of all other miRNAs in the panel, excluding the one being normalized.
• Normalization to mean considering expression level (NME): miRNAs were classified as tumor suppressors or oncogenes based on average Ct differences between NSCLC patients and HD. Tumor suppressor miRNAs (upregulated in HD and downregulated in NSCLC) were normalized against the mean of oncogenic miRNAs (upregulated in NSCLC and downregulated in HD), and vice versa. This approach accounts for biological roles of miRNAs and allows correcting mutual influence of functional groups.
• Normalization to mean considering function (NMF): miRNAs were classified into tumor suppressor and oncogenic groups based on two criteria: a) expression differences between NSCLC and HD; b) only miRNAs with statistically significant intergroup differences were included. Tumor suppressor miRNAs were normalized to the mean of oncogenic, and vice versa. This method, like the previous one, incorporates biological function but applies stricter group selection criteria.

Each normalization approach was implemented in R. Subsequently, statistical filtering of variables was performed based on distribution analysis using the Shapiro-Wilk test (if normal distribution, t-test was used to identify intergroup differences; if not, the non-parametric Kruskal-Wallis test was applied). Only features with significant intergroup differences (p < 0.05) were included in the final dataset.

Next steps included feature space expansion and cleaning, class balancing via bootstrapping, and 100-fold bootstrap on training data to select optimal regularization parameter λ (each time the sample was randomly split into training (30%) and testing (70%) subsets). The median λ from the 100 repeats was used to build the final LASSO model on the training data. Informative features were selected based on LASSO coefficients and used to build a diagnostic model based on Random Forests (RF) with 100-fold bootstrapping (each time the sample was randomly split into training (30%) and testing (70%) subsets, averaged results for test sets are presented in the Results section).

This combined analysis approach showed the highest effectiveness in a previous comparative study of machine learning methods [14]. The script also included data visualization: average ROC curve from 100 repeats, mean SHAP (Shapley Additive exPlanations, importance ranking) values for all donors (overall feature importance for the model).

The models were evaluated by metrics averaged over 100 bootstrap iterations:

• Accuracy (mean ± SD): proportion of correct predictions among all observations; overall classification quality. Calculated as the sum of true positives and true negatives divided by total number of cases; range 0-1.
• F1 score (mean ± SD): harmonic mean between Precision and Sensitivity. Balances the ability to detect positive diagnoses while minimizing false positives; range 0-1 where 1 indicates perfect balance.
• AUC (mean ± SD): area under the ROC curve. Measures model’s ability to distinguish classes at different probability thresholds; ranges from 0.5 (random chance) to 1 (perfect).
• Precision (mean ± SD): proportion of correctly classified positive cases among all cases predicted as positive by the model; important when false positives are risky; range 0-1.
• Sensitivity (mean ± SD): proportion of true positives among all actual positives (true patients); range 0-1.
• Specificity (mean ± SD): proportion of correctly identified negatives among all negatives; higher specificity means fewer false positives; range 0-1.
• OverfitGap (mean ± SD): metric showing model overfitting degree, e.g., difference between performance on training and testing sets; lower is better; expressed as fraction or percent.
• MCC (mean ± SD) (Matthews Correlation Coefficient): balanced metric considering all confusion matrix elements; ranges from -1 (complete error) to +1 (perfect); 0 means random chance.
• Jaccard Index: feature selection stability, reflecting similarity between feature sets in different repetitions or model variations; used to assess reliability and reproducibility of feature selection in diagnostic and prognostic models; defined as the intersection size divided by union size of two feature sets; ranges 0 (no overlap) to 1 (perfect overlap).

These metrics jointly assess the accuracy and reliability of diagnostic models and help select the most optimal depending on clinical conditions and data analysis tasks. Models with the highest values for all listed metrics except OverfitGap, which should be minimized, possess the greatest diagnostic significance. Additional evaluation metrics included training and inference time, i.e., the time spent on data processing and obtaining results.

3. Results

The analysis included data where the Ct values of investigated miRNAs were below 40; all samples met this criterion. The Ct values for the spike-in control cel-miR-39 across all samples averaged 25 with a variability of ±1. Its consistent Ct across samples indicates uniform RNA extraction and reverse transcription efficiency, validating the reliability of microRNA quantifications in the study.

A comparative analysis of various miRNA normalization methods was performed, assessing their effectiveness using key model quality metrics and computational parameters. No statistically significant differences were found between the two variants of mean normalization methods (OM and EM); the number of features in the final model was identical (5 out of 17). The only difference was a higher average training time for the normalization method using the exclusive mean of miRNAs (EM), likely due to computational characteristics, but it did not significantly affect model quality. Therefore, both variants are further considered a single method.

Normalization by functional group (NMF) - separating oncogenic and tumor suppressor miRNAs, using a pre-selection of miRNAs with significant intergroup differences - resulted in selection of only miR-324 with modest diagnostic performance: AUC 0.73, sensitivity 83%, specificity 40%. Due to low diagnostic efficiency, this method was excluded from further analysis. The comparative results of other normalization methods with evaluation of model quality metrics are presented in Table 4.

The final number of features in diagnostic panels ranged from 8 to 3 (Table 2, Table 3 and Table 4, Figure 1). Information on the set and representation of miRNAs in the final diagnostic models is shown in Table 2 and Table 3.

The pairs selected in the PN diagnostic model included 10 miRNAs (8 features) (Table 2). Twelve miRNAs were part of the 8 features in the TN diagnostic model. Thirteen microRNAs were included among the 7 features in the QN diagnostic model.

The final OM and EM normalization models included 5 features: miR-30, miR-205, miR-660, miR-222, and miR-125. Three of these microRNAs also appeared in the functional group-based normalization (NME) diagnostic panel: miR-30, miR-660, miR-125.

None of these normalization methods' final models included miR-31 or miR-375, and miR-200 and miR-324 appeared only once (Table 2).

The AUC metric reached its maximum in TN and QN methods, with PN showing a marginal difference (Figure 2). AUC values dropped notably for normalization to means (OM and EM) and even more for normalization considering miRNA functions (NME). Similar patterns were found for recall, sensitivity, specificity, and MCC (Table 4).

The overfit gap metric was minimal in PN, TN, and QN, indicating absence of overfitting in diagnostic panels based on these normalizations, but it increased substantially across OM, EM, and NME.

All examined methods demonstrate good control of overfitting with low overfit gap values, with the lowest values observed in PN, TN, and QN, despite the increased number of markers analyzed.

Training and inference times were lowest for pairwise normalization (0.02 s), increasing sequentially for TN, QN, functional group, and mean normalizations. Memory consumption was minimal and did not differ significantly among the methods.

Feature selection stability, measured by the Jaccard index, was highest in normalization accounting for functional groups, associated with a reduced number of selected features. This was followed by TN, PN and QN, whereas stability was lower for mean normalizations. The low stability observed for OM and EM normalizations likely indicates averaging of feature importance under these methods, resulting in different features being selected across bootstraps but with overall similar diagnostic significance of the resulting models.

Thus, normalization methods based on one or a few miRNAs (PN, TN, and QN) provide the highest accuracy and balance in models. Though similar in performance, TN achieves the best diagnostic results.

4. Discussion

The question of optimal normalization of miRNA data obtained by RT-qPCR remains unresolved. Current studies largely focus on identifying single miRNA normalizers [15,16] or several miRNAs for this purpose [17]. Rarely, normalization to initial sample amount, total nucleic acid concentration, or volume is used [18]. This study presents a comparative analysis of several miRNA normalization methods applied to RT-qPCR data aiming to optimize the search for NSCLC diagnostic markers.

Our results demonstrate high effectiveness of normalization methods using 1-3 miRNA normalizers for building NSCLC diagnostic models. High AUC, accuracy, sensitivity, and specificity values achieved with PN indicate a superior capacity to reduce technical variability and enhance biological signal. These findings align with previous reports; Boeri et al. [8] proposed pairwise normalization to increase RT-qPCR sensitivity and specificity in oncological microRNA analysis. Similarly, Schwarzenbach et al. [5] and Meyer et al. [4] reported that using multiple miRNAs as normalizers better compensates technical variation and stabilizes models compared to normalization with universal housekeeping miRNAs.

Normalization to the mean and to the mean considering classification of miRNAs into oncogenic and tumor suppressor groups showed comparatively poorer classification performance and lower feature selection stability than PN, TN, and QN methods. This likely reflects over-averaging of biologically meaningful variation or insufficient adaptation to clinical data specifics - a conclusion supported by literature noting that excessive averaging during normalization can obscure biologically important miRNA differences and reduce diagnostic power [7].

Normalization by miRNA functional groups is conceptually attractive as it incorporates biological roles and interactions of markers. However, our results showed relatively low efficiency and stability for this approach compared to PN, TN, and QN. MiRNAs regulate numerous genes and pathways, with one miRNA targeting multiple genes and vice versa [19]. In tumors, miRNAs can compensate for each other’s expression changes, creating complex nonlinear interactions that challenge simplified normalization models based on group averaging. Reduced accuracy, AUC, sensitivity, and specificity suggest that functional group averaging dilutes important individual expression differences, introducing biological ambiguity. Methods in the PN-QN group better capture individual relative changes. This aligns with evidence that miRNA biological heterogeneity and complex regulation limit simplified mean-based normalizations [4,5]. Nonetheless, such approaches may enhance interpretability and provide functionally sound normalization in specific research contexts, especially in pathogenetic and therapeutic studies (e.g., suppressing oncogenic miRNAs and restoring tumor suppressors).

The computational efficiency depends on normalization method; PN requires the shortest training and inference times, consistent with the fundamental algorithmic load differences reported by Mestdagh et al. [9], representing a practical advantage for clinical and research applications necessitating resource minimization.

The number of miRNAs included in diagnostic models ranged from 3 to 8 (out of 17 studied), balancing biological signal completeness and model parsimony to enhance stability and interpretability. The PN and TN models included 8 features each, while the QN model had 7. The total unique miRNAs increased from PN to QN, reflecting higher-dimensional feature space expansion. The PN model comprised 10 miRNAs; miR-19, -425, -125, -205, -92, and -660 predominated, forming diagnostic pairs (8 features total). TN increased unique miRNAs to 12 (adding miR-133, -222, -27; dropping miR-30), with miR-660 (4 features) and miR-144 (3 features) dominating, suggesting diagnostic relevance and likely involvement in disease-associated biological processes. QN altered the panel by adding miR-324 and -200, totaling 14 miRNAs forming 7 diagnostic features. Overall, PN, TN, and QN identify overlapping marker miRNAs; adding new miRNAs does not always yield additional diagnostic features. Considering feasibility, simplicity, performance, and cost, pairwise normalization is optimal, though TN and QN can validate marker selection or improve diagnostic metrics and provide insights into miRNA regulatory interactions implicated in pathology.

Most studied miRNAs featured in PN, TN, and QN diagnostic signatures, appearing in numerator, denominator, or with varying SHAP importance signs, complicating cross-method data analysis but allowing comparisons. For example, miR-19 appears 2, 2, and 3 times in PN, TN, and QN respectively, always in the numerator with positive SHAP importance, similarly for miR-22, suggesting lower expression in patients relative to HD. This partly aligns with literature describing both oncogenic (miR-19 [20,21]) and tumor suppressor (miR-22 [22,23,24]) roles. MiR-425 is largely in denominators with mostly positive SHAP values (except 1 of 7 unique markers), indicating higher expression in patients, consistent with its role in promoting NSCLC growth via the PTEN/PI3K/AKT pathway [25,26].

All miRNAs from mean normalization models (miR-30, miR-205, miR-660, miR-222, miR-125) also appear in PN, TN, and QN panels; miR-30 exclusively in PN, miR-222 in TN and QN. MiR-125 and miR-660 feature in all five diagnostic panels (including NME); miR-205 in four, highlighting their robustness as NSCLC biomarker candidates confirmed in literature [27,28,29,30,31,32,33,34,35,36].

MiRNAs appearing only once or twice in models were few (two in each case), suggesting significant overlap between normalization variants and indicating lower biomarker potential for those (miR-27, -133, -200, -324). MiR-31 and miR-375 did not appear in any model, denoting low diagnostic value and limited utility in NSCLC marker panels.

In summary, PN, TN, and QN normalization methods effectively analyze RT-qPCR data to identify potential cancer markers. Other methods, such as normalization to overall or exclusive means (OM, EM) and mean normalization considering oncogenic and tumor suppressor marker grouping, are less effective for marker discovery and diagnostic system building but provide supplementary insights into disease mechanisms and help identify biologically significant, stably represented candidate miRNA markers in biological fluids. These methods are recommended as complementary tools, at least in early marker panel selection stages.