Decision variables to be discovered in modelling high-dimensional omics data for cancer studies *

Abstract

High-dimensional omics data are often contaminated by sources of unwanted variations caused by platforms, batches, or other external factors. These interferences and noise can obscure critical signals related to cancer. Contaminated data are modeled as a combination of variables derived from the phenotype of interest (POI) and confounding factors. To identify these variables, a novel method called Decision Variable Analysis (DVA) is proposed. The novelty of DVA is to iteratively extract independent decisive variables for modeling the data. Specifically, a priori knowledge introduced as the definite variable linked with POI is removed from data through a residual operation. The number of variables is estimated from the residual matrix based on the zero gradient of singular values, rather than relying on random matrix theory or principal components analysis, which can produce unreliable results when the number of features exceeds the number of samples. Applications of DVA to both synthetic and real data demonstrate superior performance in identifying variables compared to conventional approaches. Improvements offered by DVA are illustrated across high-dimensional omics datasets, particularly those with smaller sample sizes relative to the number of features on different platforms. The results indicate that DVA is an effective method for dissecting sources of variation in high-dimensional data with disturbances.

Keywords

High-dimensional data decision variable analysis confounding factors deconvolution

1. Introduction

One of the most important issues for high-dimensional omics data analysis is the identification of meaningful markers.^1,2 For instance, differentially expressed genes or differentially methylated CpG sites were identified between normal and cancer samples from the Cancer Genome Atlas (TCGA).^3,4 However, these large-scale profiling data are subjected to heterogeneity caused by confounding factors (CFs), resulting from both technical and biological sources, e.g. ambient conditions during sample preparation and processing, sample storage conditions, or any other technical discrepancy between or within ‘batches’, which describe some samples being processed differently than others.^5–7

Confounding factors are often presented against a background of biological variables, which are not the interesting factors in previous studies. A typical example is the gene expression analysis, which involves identifying differential expressed genes under various conditions, such as the cancer state. Not only changes in expression are related to the cancer state, but also the age, sex or body mass index of the individuals. Some genes are observed to have differential expression^8–10 pertaining to cancer,^11,12 some corresponding to these unwanted factors.^14,13 In the context of DNA methylation studies, many epigenome-wide association studies (EWAS) are performed using blood, as it is easily accessible. However, blood is a heterogeneous tissue that consists of several different cell types. It is critical to account for cell composition as it could influence the differential methylation calling.^16,15 Other genetic, demographic, and environmental factors have also been presented to have an impact on the omics data profiling.^18,19,17 If these factors are not excluded in the large-scale profiling analysis when identifying differential expression or methylation, this may introduce additional spurious signals resulting from the effect of these factors or induce ‘regular’ noise in the data.²⁰

Normalization techniques are usually used for adjusting ‘regular’ noise or systematic variation resulting from laboratory and technical sources. Some studies indicate that unknown variations remain unremoved after normalizing.²¹ Thereby, a few methods have been proposed to solve the problem of unwanted variation in large profiling studies. Conventional techniques are linear modeling and the empirical Bayes (EB) framework, which deem the confounding factors as additional variables. ComBat is a method to assume the factors inducing the unwanted variation are known, especially for small data samples.^22,23 However, many factors can hardly be obtained in practice, and it is difficult for researchers to select all factors into the model even if the factors (variables) are known because there are too many of them. Moreover, factors may interact with each other and complicate the influence on high-dimensional omics data. Another method assumes that the sources of the unwanted variations are unknown. For instance, remove unwanted variation (RUV), 2-step (RUV-2) or RUV-like method uses negative control features that are not associated with phenotype of interest (e.g. cancer, age, or sex) but are affected by unknown variations to estimate the unwanted factors.¹⁴ Negative controls play a crucial role in these methods. Negative controls are features (e.g. genes) that must be known $a priori$ as unassociated factors with the phenotype of interest (POI). However, in reality, the negative controls are generally partially known or utterly ignorant.

The unwanted or unknown factors can be derived from the data using tricks of component analysis, e.g. singular value decomposition (SVD), principal components analysis (PCA),²⁴ and independent variables analysis (ICA).²⁵ Surrogate variable analysis (SVA) has been proposed to model the potential factors for feature selection in omics data based on SVD.²⁶ Whereafter, independent surrogate variable analysis (ISVA) has been proposed using random matrix theory (RMT) based on PCA to deconvolve CFs.²⁷ Whereas, the result is not reliable if the omics data table with rows (features) greater than columns (samples) is processed by PCA-based methods.²⁸ Furthermore, there are no convincing explanations for the meanings of principal components or surrogate variables associated with the sources of variations, whatever unmeasured or measured factors resulted from both biological and technical sources. To solve the small sample size problem when dealing with high-dimensional (features) omics data and model the CFs better, a novel framework called decision variable analysis (DVA) is proposed by eliminating the influence of the main POI corresponding to the decision variable.

The motivation of DVA is that a priori knowledge is closely related to the primary POI, which can be transformed into a decision variable in a linear regression equation. After removing the effect of the known decision variable on high-dimensional omics data, the residual matrix is obtained, and the residual part is decomposed into statistically independent variables related to phenotypes or potential variables. The proposed algorithm improves the SVA/ISVA in two aspects: (i) The nonlocal gradient is used to estimate the number of components instead of RMT based on PCA in ISVA because PCA can hardly obtain reliable results when processing high-dimensional samples of small size. (ii) Considering the biological significance of the data matrix, the independent decisive variables (DVs) of CFs is progressively extracted against POI interference and model misspecification instead of uncorrelated surrogate variables (SVs) as SVA does.

The rest of this article is organized as follows. Section 2 introduces the concepts and method of DVA, along with simulated and factual data. In Section 3, for the effectiveness of proposed method, DVA is implemented on the simulated data and real datasets in the context of gene expression and DNA methylation, respectively. The results show that the inferred decision variables proposed by this method can capture confounding factors more efficiently when compared to other methods. Section 4 summarizes that DVA can find decisive variables for cancer study and hold for different high-dimensional omics data.

2. Method

2.1. Review of SVA/ISVA

A biological data matrix X is composed of m rows and n columns, where m and n denote the number of features (genes, CpGs, etc.) and samples, respectively. The number of features is greater than the number of samples (m > n). To estimate confounding factors ( $u_{i}$ ), the model of SVA uses this form:

x_{i j} = u_{i} + f_{i} (y_{i}) + ε_{i j}

(1)

where,

x_{i j}

is the element in matrix X,

f_{i} (y_{i})

is the relationship between POI and the feature

i

ε_{i j}

is the random noise.

SVA perform a SVD on the residual matrix R which is defined by $R = X - Y, Y$ is the regression estimation for known factors:

R = U Λ V^{T}

(2)

where,

U

is an m

\times

m orthogonal matrix,

Λ

is a diagonal n

\times

n matrix with non-negative real numbers on the diagonal,

V

is an n

\times

n orthogonal matrix, and

V^{T}

is the transpose of

V

The singular vectors of SVD represent the variation associated with the factors, including some biological factors and experimental factors. These factors are mixed and constitute the potential CFs. To untangle the mixed effects in surrogate variables derived from singular vectors, ISVA suggests using independent component analysis (ICA) on the residual matrix R. The model of ISVA based on ICA can be expressed as:

x = A s(t) + ϵ = \sum_{i = 1}^{m} a_{i} s_{i} + ε_{i j}

(3)

where,

x = (x_{1}, x_{2}, \dots, x_{n})^{T}

denotes the n-dimensional biological observations,

A

is the mixed matrix,

a_{i}

is the i column of

A

, and

s_{i}

is the i row of

s(t)

s(t) = (s_{1}, s_{2}, \dots, s_{m})^{T}

is the statistically independent source signal, and

ϵ

ε_{i j}

is the random noise.

Introducing ICA is to estimate the mixed matrix $A$ and the source $s(t)$ , whatever the technical or biological source. This is solved by adaptively computing the $w$ (weight vector) and establishing a cost function that minimizes the mutual information or maximizes the non-Gaussianity of the calculated $S = W^{T} X$ . The original source $S$ then can be reconstructed by multiplying the observations $X$ with the inverse of matrix $A$ . The fastICA algorithm usually uses a nonquadratic nonlinear function to measure non-Gaussianity.²⁹ An R-package is already available and has many successful applications. ISVA is an estimation of K components (dimensions) imported in ICA based on RMT.³⁰ RMT estimates the number of components of a data covariance matrix by comparing the eigenvalues of the data covariance matrix derived from PCA. The number of eigenvalues of the data covariance matrix larger than the analytical maximum is the estimation of K significant components (dimensions). Subsequently, ICA is implemented in this K-dimensional subspace of the data matrix.

2.2. Decision variable analysis

To determine the factors as variables related to the cancer, the high-dimensional profiling data can be deemed as a fusion pattern by different factors: the first one is definite POI, the second ones are known or unknown factors as confounding factors (CFs), and the third is random or stochastic fluctuation. First, the data should be resolved into (almost) independent components. For instance, if samples are obtained from two related but distinct cohorts and profiled on microarrays, the array and cohort can be considered as two potential and independent factors. Whereas. these factors are modeled as ‘surrogate variables’ only derived from SVD in SVA. These ‘surrogate variables’ are linearly uncorrelated and not independent, which is inconsistent with the reality. Thus, these factors are need to be deconvolved into statistically independent components using ICA. In view of ICA, identifying different impact factors can be deemed as ‘blind source signal separation’. To trace different variables, the proposed method should be designed on the basis of independent non-Gaussian components corresponding to the POI and CFs. Although ISVA also uses ICA to deconvolve the confounding factors, it applies ICA over the K-dimensional subspace of the data matrix reduced by RMT based on PCA. PCA is a crucial step in the dimensionality estimation of ISVA. In this step, after normalizing the data matrix, PCA is conducted using the eigenvalue decomposition of the data covariance matrix. It should be noted that the meaning of the biological data columns represents samples in normalization. When the sample number is smaller than the feature number, PCA cannot obtain reliable results. Thus, ISVA cannot achieve effective results in real biological data, whatever DNA methylation data or gene expression data; this is proved by the experiments in the following sections. For extensive adaptability to the special biological data, the decisive variables (DVs) of CFs are progressively purified against interference. The flow chart of the proposed method for purifying DVs, called decision variable analysis (DVA), is shown in Figure 1.

(1)
The data are correlated to known variables (if there are $a priori$ knowledge) with the smallest P-value which means remarkable statistical significance as the POI.
(2)
The effect of known variable (y) is removed from the data by $R = X - L (y)$ , where $L (y)$ means the linear part of the data matrix $X$ affected by y and $R$ is residual matrix.
(3)
The number of uncovered variables is estimated when the non-local gradient^31,32 of the singular value approximates zero (0.01).
(4)
Independent component analysis (ICA) is performed on the residual matrix $R$ in (2): $R = S A$ , where $A$ is the mixed matrix of residual matrix and $S$ is the source matrix of residual matrix.
(5)
The mixed matrix of residual matrix is correlated to the data to find the most related contents in the data.
(6)
The data matrix is clipped by features that are significantly correlated ( $P < 0.05$ ) with the mixed matrix of residual matrix: $C \subseteq X$ , where $X$ is data matrix and $C$ is a clipped subset of $X$ .
(7)
Again, ICA is conducted on $C$ in (6): $C = S_{c} A_{c}$ , where $A_{c}$ is the mixed matrix of $C$ and $S_{c}$ is the source matrix of $C$ .
(8)
The vectors of $A_{c}$ are selected as decision variables which have the largest Pearson correlations with the mixed matrix of residual matrix, which means the effects of CFs are purified again.

Figure 1.
Flowchart of the proposed decision variable analysis.
2.3. Simulated dataset

We generate the simulated data with 50 samples representing simulated units based on Lee et al.²⁶ and Teschendorff et al.²⁷ for comparing to other methods with the same standard. The simulated data contain three variables ( $V_{1}, V_{2}, V_{3}$ ) corresponding to a POI and two confounding factors (CF1 and CF2), respectively. Let $V_{1}$ a binary variable. $V_{1} = 0$ denote 25 randomly selected samples that are categorized into one class, and $V_{1} = 1$ represents the other 25 samples that are classified into another category (Figure 2(B)). Similarly, two CFs are also assumed to be binary variables. Half of the samples are assigned randomly with $V_{2} = 1$ or $V_{3} = 1$ and the other half with $V_{2} = 0$ or $V_{3} = 0$ (Figure 2(C) and (D)). Then, given one sample s, a three-tuple of variables $V_{s} = (V_{1}, V_{2}, V_{3})$ is built. Since each variable in triad $V_{s}$ just takes value 0 or 1, there are 8 different states of $V_{s}$ . The influences of variables are formulated as follows.

(1)
$V_{s} = (0, 0, 0)$ means that there are null features with no influence of variables
$p (x | V_{s}) \sim N (0, 1)$
(4)
(2)
$V_{s} = (0, V_{2}, 0)$ means that there are null features influenced by the CF1
$p (x | V_{s}) \sim δ_{V_{s}, x 1 z} N (e_{2}, 1) + δ_{V_{s}, x 0 z} N (0, 1)$
(5)
(3)
$V_{s} = (0, 0, V_{3})$ means that there are null features influenced by the CF2
$p (x | V_{s}) \sim δ_{V_{s}, x y 1} N (e_{3}, 1) + δ_{V_{s}, x y 0} N (0, 1)$
(6)
(4)
$V_{s} = (0, V_{2}, V_{3})$ means that there are null features influenced by two CFs
$p (x | V_{s}) \sim δ_{V_{s}, x 11} N (e_{2} + e_{3}, 1) + δ_{V_{s}, x 01} N (e_{3}, 1) + δ_{V_{s}, x 10} N (e_{2}, 1) + δ_{V_{s}, x 00} N (0, 1)$
(7)
(5)
$V_{s} = (V_{1}, 0, 0)$ means that there are true positives with no influence of CFs
$p (x | V_{s}) \sim δ_{V_{s}, 0 y z} N (0, 1) + δ_{V_{s}, 1 y z} (π_{- 1} N (- e_{1}, 1) + π_{1} N (e_{1}, 1))$
(8)
(6)
$V_{s} = (V_{1}, V_{2}, 0)$ means that there are true positives influenced by CF1
$\begin{aligned} p (x | V_{s}) & \sim δ_{V_{s}, 00 z} N (0, 1) + δ_{V_{s}, 01 z} N (e_{2}, 1) + δ_{V_{s}, 10 z} (π_{- 1} N (- e_{1}, 1) + π_{1} N (e_{1}, 1)) \\ + δ_{V_{s}, 11 z} (π_{- 1} N (- e_{1} + e_{2}, 1) + π_{1} N (e_{1} + e_{2}, 1)) \end{aligned}$
(9)
(7)
$V_{s} = (V_{1}, 0, V_{3})$ means that there are true positives influenced by CF2
$\begin{aligned} p (x | V_{s}) & \sim δ_{V_{s}, 0 y 0} N (0, 1) + δ_{V_{s}, 0 y 1} N (e_{3}, 1) + δ_{V_{s}, 1 y 0} (π_{- 1} N (- e_{1}, 1) + π_{1} N (e_{1}, 1)) \\ + δ_{V_{s}, 1 y 1} (π_{- 1} N (- e_{1} + e_{3}, 1) + π_{1} N (e_{1} + e_{3}, 1)) \end{aligned}$
(10)
(8)
$V_{s} = (V_{1}, V_{2}, V_{3})$ means that there are true positives influenced by all CFs
$\begin{aligned} p (x | V_{s}) & \sim δ_{V_{s}, 000} N (0, 1) + δ_{V_{s}, 001} N (e_{3}, 1) + δ_{V_{s}, 010} N (e_{2}, 1) + δ_{V_{s}, 011} N (e_{2} + e_{3}, 1) \\ + δ_{V_{s}, 101} (π_{- 1} N (- e_{1} + e_{3}, 1) + π_{1} N (e_{1} + e_{3}, 1)) + δ_{V_{s}, 110} (π_{- 1} N (- e_{1} + e_{2}, 1) \\ + π_{1} N (e_{1} + e_{2}, 1)) + δ_{V_{s}, 111} (π_{- 1} N (- e_{1} + e_{2} + e_{3}, 1) + π_{1} N (e_{1} + e_{2} + e_{3}, 1)) \end{aligned}$
(11)
where $δ_{x^{^{'}} y^{^{'}} z^{^{'}}, x y z}$ means the triple Kronecker delta: $δ_{x^{^{'}} y^{^{'}} z^{^{'}}, x y z} = 1$ if and only if $x^{^{'}} = x, y^{^{'}} = y$ , and $z^{^{'}} = z$ , otherwise $δ_{x^{^{'}} y^{^{'}} z^{^{'}}, x y z} = 0$ and $π_{- 1} = π_{1} = 0.5$ .

Figure 2.
An example of simulated data. (A) A heatmap of a simulated data containing 2,000 genes measured on 50 samples. (B) Genes 1801–2000 in this data are affected by the POI. The value of POI is taken either 0 or 1 for each sample, which indicates the sample belongs to one or the other class. (C) Genes 1201–1800 in this data are affected by both two confounding factors (CFs), i.e. CF1 and CF2. These two CFs are also binary variables that take value either 0 or 1. (D) Genes 601–1200 in this data are affected by one confounding factor (CF1 or CF2). Genes 1-600 are not affected by any POI or CFs.

Base on Teschendorff et al.²⁷ and McAlpine et al.,³³ ten percent of features derived from equations (4) to (11) of triad $V_{s}$ are assumed to be true positives (TPs) which can be differentiated between two classes of POI. Two CFs are modeled as follows: each CF affects 50 features independently, and 50 features are affected by both CFs. The noise is modeled by a Gaussian of the mean of zero and the unit variance of N(0, 1). Moreover, $e_{1}, e_{2}$ , and $e_{3}$ can be used to represent the relative effect scope of POI, CF1, and CF2 on the simulated data, respectively. Let $e_{2} = e_{3}$ for simplicity. By and large, the whole generated simulated data can be visualized in Figure 2.
2.4. Real biological datasets: DNA methylation

Four DNA methylation datasets produced by Illumina’s Infinium Human Methylation 450k Beadchips are used in this study. The Illumina Human Methylation 450 BeadChip measures the methylation status of CpGs at 485,577 sites across the human genome. The data is expressed as $x_{i j}$ where i denotes the CpG and j represents the sample. The methylation values of the CpGs are normalized between 0 (unmethylated locus) and 1(methylated) which obey an approximate $β$ -valued distribution based on Xie et al.³⁴ The $β$ is defined as the ratio of methylated to the combined intensity values because it is easier to interpret the related biological process:³⁵

β = \frac{m a x (M, 0)}{m a x (U, 0) + m a x (M, 0) + 100}

(12)

where, U and M are the unmethylated and methylated intensity values of the probe, and 100 is a correction term to regularize probes of low total signal intensity (i.e. probes of

U + M \approx 0

after subtracting background signal). The

x_{i j}

then can be expressed as

x_{i j} = β_{i j}

, where

β_{i j}

is the normalized methylation value at the

i

CpG locus of

j

sample.

Datasets 1 $\sim$ 4 are DNA methylation datasets from the Cancer Genome Atlas (TCGA) (https://tcga-data.nci.nih.gov/tcga).

Dataset 1: the DNAm dataset is composed of 395 samples from patients with stomach adenocarcinoma and normal controls (a total of 137 women and 258 men).

Dataset 2: this DNAm dataset consists of 472 samples from patients with skin cutaneous melanoma and normal people (a total of 179 women and 293 men).

Dataset 3: this dataset consists of 310 samples from a total of 141 women and 169 men (38 healthy people and 272 patients with colon adenocarcinoma).

Dataset 4: this dataset consists of 102 samples from 7 healthy people and 95 patients with rectum adenocarcinoma, a total of 47 women and 55 men.

In these Illumina Infinium Human Methylation 450k arrays, three main CFs are taken into account: sex, age and race variations. It is noteworthy that the variations caused by these factors have not been removed by quantile normalization procedures. Specifically, sex, age, and race are taken as confounding factors in Datasets 1

\sim

3 and sex and age as confounding factors in Dataset 4 because there are too many missing values in the race factor of this dataset.

Another readily obtaining DNA samples for methylation studies is derived from the whole blood. However, blood contains many functionally distinct cell populations in various proportions. Such confounding due to cellular heterogeneity might affect the interpretation of methylation. To validate the identifiability of the proposed algorithm to the cellular heterogeneity of blood samples, the DNA methylation data (eighteen samples from adult donors) from GEO (GSE77797, Dataset 5) is used. Cancer status is used as the phenotype of interest. The proportions of B cells, CD4T, CD8T, natural killer cells, purified granulocytes, and monocytes in mixed white blood cells and whole blood are used as confounding factors. The proposed algorithm is expected to identify these factors for the deconvolution of adult human whole blood in the DNA methylation study.

2.5. Real biological datasets: gene expression

Datasets 6 $\sim$ 9 are gene expression datasets from TCGA (https://tcga-data.nci.nih.gov/tcga).

Dataset 6: this dataset consists of 16861 genes and 529 samples from 471 patients with lung adenocarcinoma and 58 healthy people (a total of 289 women and 240 men).

Dataset 7: this dataset is composed of 16,633 genes and 587 samples from 515 patients with kidney renal clear cell carcinoma and 72 normal, a total of 202 women and 385 men, respectively.

Dataset 8: this gene expression dataset contains 16,990 genes and 448 samples from a total of 158 women and 290 men (35 healthy people and 413 patients with stomach adenocarcinoma).

Dataset 9: this dataset includes 16,161 genes and 399 samples from 50 healthy people and 349 patients with liver hepatocellular carcinoma, a total of 136 women and 263 men.

Likewise, gender, age, and race are taken as the potential confounding factors in these datasets.

2.6. Algorithm comparison using the same evaluation indexes

Using the same statistic tests and evaluation indexes, the comparison is performed to competing DVA, SVA, and standard linear regression (SLR) with no adjustment for CFs on the same simulated and real biological datasets, respectively. The statistical tests and evaluation indexes are as follows.

Let $H_{0}$ denote null features and $H_{A}$ indicate the number of features following the alternative hypothesis. Let $H_{AA}$ represent the number of features following the alternative hypothesis and affected by the CFs ( $H_{AA} < H_{A}$ ). DVA, SVA and SLR are compared without adjustment for CFs on their performance to detect the number of true positives at a given false discovery rate (FDR) threshold of 0.05, the positive predictive value $PPV = TP / H_{A}$ (where TP is the number of true positives passing the same FDR threshold of 0.05), the false negative rate $FNR = FN / H_{A}$ (where FN is the number of false negatives at the same FDR threshold), and the sensitivity to detect true positives affected by the CFs, that is, SEA = $T P_{A} / H_{AA}$ (where $T P_{A}$ is the number of features affected by the CFs passing the FDR threshold of 0.05). The identifiability of derived surrogate or independent decision variables is evaluated with CFs using the average of the best Pearson correlations (CORR) between them. Moreover, the Kolmogorov–Smirnov test (KST) is implemented to test the P-values of nulls following the uniform distribution between 0 and 1.

3. Results

3.1. Simulated data

We generate the simulated data, including 2000 features and 50 samples with a primary binary variable of POI and two CFs (also binary variables) in Section 2. It is assumed that 10% of features are TPs distinguishing for POI and have a 25% overlap for each CF. Using the same FDR threshold of 0.05 and the relative effect size, five methods widely used so far (ComBat-seq,²³ DVA, RUV-4,¹⁴ SLR, and SVA^36,37) have repeatedly run 100 times, respectively. Thereupon, performance indexes of these algorithms are illustrated in Figure 3.

Figure 3.

Performance indexes of different algorithms over 100 runs of simulated data. The algorithms for comparison are Com-Bat-seq, DVA, RUV4, SLR, and SVA. Using the same given false discovery rate (FDR) threshold of 0.05, the metrics of TP are compared passing the same FDR threshold, PPV, FNR, SEA, CORR, and KST (Please see Section 2 for more details). (A) Relative effect size = 1. (B) Effect size of CFs is set to be 3 $\times$ the effect size of the phenotype of interest.

Figure 4.

(A) The variation of 14 significant singular vectors used to estimate variables, which measured relative to the total variation in the data. The number of variables can be estimated when the parameter of non-local gradient is set to 0.01. (B) Heatmap of P-values of relationships between 14 significant singular vectors and the factors (cancer, age, race, and sex). P-values are estimated using linear model. Color codes: $P < 1 \times 10^{- 10}$ (darkred), $P < 1 \times 10^{- 5}$ (red), $P < 0.01$ (orange), $P < 0.05$ (pink), $P > 0.05$ (white).

In Figure 3(A), it is observed that DVA outperforms other methods in all indexes except KST. In terms of KST, DVA show almost the same performance with ComBat-seq, RUV4, SVA except SLR. This indicate that the null gene P-values of DVA and other methods except SLR exhibit agreement with the uniform distribution. From Figure 3, it is noted that DVA remarkably reduces the FNR and increases the PPV and SEA with lower variations compared to other algorithms. As expected, the SLR has the lowest performance because it cannot take potential confounding factors into account. SVA, RUV4, and ComBat-seq have almost the same largest variations in terms of TP, PPV, FNR, and SEA, suggesting that the performance of these methods is unstable. Moreover, the CORR of DVA is the highest among these algorithms, which indicates that the inferred decision variables (DVs) have the closest relationships with the CFs than others. Using DVs as surrogates for CFs is a better choice of which variables to contain as covariates.

In generating the simulated data, the relative effect sizes of POI and CFs are also taken into account, as previous studies suggested.²⁷ All methods run 100 times with the effect size changing from 1 to 4. Figure 3(A) and (B) illustrate the results of the effect size of CFs being assumed to be 1 and 3 times that of the POI, respectively. From Figure 3(B), the number of TP detected by DVA is the closest to the true number of TP with the smallest deviation. For most indexes, DVA has the best performance among the algorithms (Figure 3(B)). Meanwhile, DVA almost outperforms other methods in all metrics, and the performance at a large effect size of CFs is better than the counterpart at a small effect size (Figure 3(A) and (B)). The variation of performance metrics (expect KST and CORR) of DVA becomes smaller when increasing the effect sizes of CFs. Note that the variation of metrics (expect KST and CORR) of SVA has significantly decreased. This suggests that it is easier to identify the CFs for the algorithm when the effects of CFs become larger. When the effects of CFs are small, DVA can have more steady performance than others (Figure 3(A)). In another word, DVA is more stable or reliable than other methods in modeling the CFs when the confounders have a small influence on the data. From the experiments with different effect sizes, we found that as decreased effect size, the interferences stemming from CFs become smaller, and the advantages of DVA in identifying CFs are more obvious, since it is necessary to consider the stability of the algorithm in deconvolving confounding factors.

3.2. DNA methylation datasets

We implemented the proposed DVA on four DNA methylation datasets from TCGA (Datasets 1-4 as detailed in Section 2.4). In the procedure of DVA (Section 2.2), the data is firstly correlated with different variables (e.g. cancer status, age, race, and sex). Cancer status has been selected as the POI with the smallest P-value using linear regression. The number of variables is estimated by the number of significant singular vectors (SSV) decomposed by SVD from the data matrix. In Figure 4(A), it is observed that the variation of SSV decreases as the number of SSV increases. The large variation of SSV indicates that unnegligible information is contained in the dimension. Thus, the number of variables can be estimated when the variation of SSV does not change. From Figure 4(A), the slope of the variation curve of singular vectors is (almost) flat when the parameter of the non-local gradient is set to 0.01.

The heatmap of P-values of correlation between singular vectors and the factors is drawn in Figure 4(B). It is important to note that many singular vectors are significantly related to factors other than cancer, especially sex. Based on the above experiments, it is found that SVA almost has the same performance with ComBat and RUV-like methods. Furthermore, SVA-like methods are more popular than ComBat and ruv-like methods based on literature investigation.^36,37 Thereupon, SVA is applied to the same four methylation datasets using cancer status as the POI. The derived surrogate variable (SV) and decision variable (DV), which are the best correlated with the CFs based on the relationships and P-value statistics, are recorded in Table 1. It is found that all selected DVs are significantly correlated with the sex factor. The derived DVs are more related to sex factors than others, whereas SVs are not. In most cases, low correlations are found between inferred SVs and the factors ( $R < 0.3$ ); the correlation values between DVs and the factors are larger than the corresponding values between SVs and the factors. These suggest that DVA compares favorably to SVA, especially in modeling sex factor.

Table 1.
The performance of DVA and SVA on DNA methylation datasets 1-4.

Datasets Factors SVA-P SVA-R DVA-P DVA-R

Stomach adenocarcinoma Sex 0.361 0.05(1) $1 \times 10^{- 167}$ 0.94(9)*

Age $1 \times 10^{- 7}$ 0.27(2) $1 \times 10^{- 5}$ 0.23(7)

Race 0.018 0.12(1) $2 \times 10^{- 10}$ 0.33(5)

Skin cutaneous melanoma Sex 0.036 0.10(2) $4 \times 10^{- 128}$ 0.85(10)

Age 0.005 0.13(2) $3 \times 10^{- 8}$ 0.26(9)

Race 0.757 0.01(2) 0.035 0.10(8)

Colon adenocarcinoma Sex $1 \times 10^{- 83}$ 0.86(4) $1 \times 10^{- 81}$ 0.86(13)

Age 0.017 0.14(3) 0.006 0.16(7)

Race 0.254 0.07(5) 0.140 0.09(8)

Rectum adenocarcinoma Sex – – $3 \times 10^{- 40}$ 0.91(2)

Age – – 0.044 0.20(10)

Blood CD4T – – 0.095 0.41(5)

CD8T – – $8 \times 10^{- 6}$ 0.85 (7)

B cells – – $5 \times 10^{- 5}$ 0.81 (8)

Natural killer cells – – $9 \times 10^{- 7}$ 0.89 (8)

Monocytes – – 0.001 0.69 (8)

Granulocytes – – $2 \times 10^{- 11}$ 0.97 (8)

Datasets	Factors	SVA-P	SVA-R	DVA-P	DVA-R
Stomach adenocarcinoma	Sex	0.361	0.05(1)	$1 \times 10^{- 167}$	0.94(9)*
	Age	$1 \times 10^{- 7}$	0.27(2)	$1 \times 10^{- 5}$	0.23(7)
	Race	0.018	0.12(1)	$2 \times 10^{- 10}$	0.33(5)
Skin cutaneous melanoma	Sex	0.036	0.10(2)	$4 \times 10^{- 128}$	0.85(10)
	Age	0.005	0.13(2)	$3 \times 10^{- 8}$	0.26(9)
	Race	0.757	0.01(2)	0.035	0.10(8)
Colon adenocarcinoma	Sex	$1 \times 10^{- 83}$	0.86(4)	$1 \times 10^{- 81}$	0.86(13)
	Age	0.017	0.14(3)	0.006	0.16(7)
	Race	0.254	0.07(5)	0.140	0.09(8)
Rectum adenocarcinoma	Sex	–	–	$3 \times 10^{- 40}$	0.91(2)
	Age	–	–	0.044	0.20(10)
Blood	CD4T	–	–	0.095	0.41(5)
	CD8T	–	–	$8 \times 10^{- 6}$	0.85 (7)
	B cells	–	–	$5 \times 10^{- 5}$	0.81 (8)
	Natural killer cells	–	–	$9 \times 10^{- 7}$	0.89 (8)
	Monocytes	–	–	0.001	0.69 (8)
	Granulocytes	–	–	$2 \times 10^{- 11}$	0.97 (8)

$^{*}$ The best SV or DV among inferred SVs/DVs related to the confounding factors is shown in the bracket.

Figure 5.

Comparison of DVA with SVA in identifying sex effect in Dataset 1. The weights (y-axis) in the two decision variables and surrogate variables the most correlated with the sex factor are plotted against genders (x-axis). The F-statistics of a linear ANOVA with sex as the independent variable are also provided for comparing the identifiability of sex effect.

Considering a single DV or SV may not capture all possible effects, the second best DV/SV correlated with the factor is included in the analysis. The first and second best DVs/SVs correlated with sex factors are shown in Figures 5 and 6, respectively. There are bad correlations ( $R < 0.2$ ) between the second best DV/SV and the factor, although DVs are more relevant to the sex factor than SVs. This indicates that the second best DVs/SVs are related to some uncertainties or measurement errors. The corresponding F-statistics of DVA are larger than those of SVA, implying that the inferred DVs grasp more variability in gender differences (Figures 5 and 6). Nevertheless, DVA captures the CF fairly well ( $R > 0.85$ ), even if the CF is subject to measurement error or uncertainties, but not so for SVA. The derived DVs accurately model the effect that unknown or known CFs have on the data. Using DVs as covariates can achieve better results than using one or two factors or SVs as covariates.

Figure 6.

Comparison of DVA with SVA in identifying sex effect in Dataset 2. The weights (y-axis) in the two decision variables and surrogate variables the most correlated with the sex factor are plotted against genders (x-axis). The F-statistics of a linear ANOVA with sex as the independent variable are also provided for comparing the identifiability of sex effect.

In addition, it indicates that DVA can capture the factor of sex in the dataset of rectum adenocarcinoma, which contains 102 samples, whereas SVA cannot (Table 1). SVA cannot achieve results on blood dataset either. In blood dataset, pertaining to whole-blood DNA methylation, there are only eighteen samples. In this case of small samples, our proposed method can achieve promising results and model the CFs of cellular heterogeneity in blood (Table 1 and Figure 7).

Figure 7.

The performance of inferred DVs in modeling the proportions of different cell types. The values of DV (y-axis) which is the most correlated with cell populations are plotted against the proportions of cell types (x-axis).

In most cases, the inferred DVs from whole-blood methylation data have high correlations with cell populations in varying proportions with small P-values. Since cell compositions in blood as confounding factors convolve the interpretation of DNA methylation, these highly related DVs can be readily used as covariates in modeling. Besides, for Dataset 3, appreciable differences cannot be observed between DVA and SVA. This is not surprising when CFs are subject to no measurement error or uncertainties (at least in the way they affect the data). Nonetheless, the performance metrics achieved by DVA are better than those by SVA.

3.3. Gene expression datasets

We further apply DVA and SVA to four gene expression datasets (Datasets 6-9, Section 2.5). These gene expression datasets all have a low signal-to-noise ratio (SNR), and potential CFs are sex, age, and race. The derived SV/DV is best correlated with the CFs and P-values, which are listed in Table 2. It suggests that the selected DVs are all significantly related to the sex factor. The derived SVs/DVs are more correlated with the gender factor than the other two factors. A low correlation between SVs (from SVA) and sex factor in Dataset 7 (stomach adenocarcinoma). The correlation values between DVs and sex factor are all larger than those between SVs and sex, which also indicates that DVA models the factor better than SVA-like method.

Table 2.
The performance of DVA and SVA on gene expression datasets 6-9.

Datasets Factors SVA-P SVA-R DVA-P DVA-R

Lung adenocarcinoma Sex $3 \times 10^{- 50}$ 0.62(5) 0 0.98(16)*

Age $5 \times 10^{- 3}$ 0.13(38) $3 \times 10^{- 4}$ 0.17 (2)

Race $8 \times 10^{- 5}$ 0.18 (20) $4 \times 10^{- 2}$ 0.09 (12)

Kidney renal clear cell carcinoma Sex $3 \times 10^{- 44}$ 0.54 (9) $1 \times 10^{- 224}$ 0.91 (11)

Age $2 \times 10^{- 5}$ 0.18 (9) $8 \times 10^{- 4}$ 0.14 (11)

Race $1 \times 10^{- 7}$ 0.22 (17) $1 \times 10^{- 3}$ 0.13 (11)

Stomach adenocarcinoma Sex $1 \times 10^{- 12}$ 0.35 (10) $5 \times 10^{- 77}$ 0.77 (4)

Age $7 \times 10^{- 4}$ 0.17 (19) $2 \times 10^{- 3}$ 0.16 (10)

Race $7 \times 10^{- 6}$ 0.23 (27) $1 \times 10^{- 2}$ 0.13 (13)

Liver hepatocellular carcinoma Sex $3 \times 10^{- 28}$ 0.52 (10) $3 \times 10^{- 105}$ 0.84 (9)

Age $4 \times 10^{- 8}$ 0.27 (15) $1 \times 10^{- 7}$ 0.26 (4)

Race $4 \times 10^{- 9}$ 0.29 (15) $2 \times 10^{- 5}$ 0.21 (9)

Datasets	Factors	SVA-P	SVA-R	DVA-P	DVA-R
Lung adenocarcinoma	Sex	$3 \times 10^{- 50}$	0.62(5)	0	0.98(16)*
	Age	$5 \times 10^{- 3}$	0.13(38)	$3 \times 10^{- 4}$	0.17 (2)
	Race	$8 \times 10^{- 5}$	0.18 (20)	$4 \times 10^{- 2}$	0.09 (12)
Kidney renal clear cell carcinoma	Sex	$3 \times 10^{- 44}$	0.54 (9)	$1 \times 10^{- 224}$	0.91 (11)
	Age	$2 \times 10^{- 5}$	0.18 (9)	$8 \times 10^{- 4}$	0.14 (11)
	Race	$1 \times 10^{- 7}$	0.22 (17)	$1 \times 10^{- 3}$	0.13 (11)
Stomach adenocarcinoma	Sex	$1 \times 10^{- 12}$	0.35 (10)	$5 \times 10^{- 77}$	0.77 (4)
	Age	$7 \times 10^{- 4}$	0.17 (19)	$2 \times 10^{- 3}$	0.16 (10)
	Race	$7 \times 10^{- 6}$	0.23 (27)	$1 \times 10^{- 2}$	0.13 (13)
Liver hepatocellular carcinoma	Sex	$3 \times 10^{- 28}$	0.52 (10)	$3 \times 10^{- 105}$	0.84 (9)
	Age	$4 \times 10^{- 8}$	0.27 (15)	$1 \times 10^{- 7}$	0.26 (4)
	Race	$4 \times 10^{- 9}$	0.29 (15)	$2 \times 10^{- 5}$	0.21 (9)

$^{*}$ The best SV or DV among inferred SVs/DVs related to the confounding factors is shown in the bracket.

The similar improvements in DVA across four gene expression datasets can ben found in Table 2. The inferred DVs model the CFs better because DVs are mutually independent statistically, whereas SVs are irrelevant variables and not independent variables. The independent variables are preferred when modeling the potential effects of the data. Furthermore, independent DVs model the confounders directly from the data rather than measuring different phenotypes, which may be subject to measurement error or uncertainties. The performance of DVs can remain stable when the data has measurement errors or uncertain interference. In other words, when phenotypes are hard to measure or subject to measurement error or uncertainties, using DVs as covariates is an advisable choice.

4. Discussion and conclusion

A novel framework for DVA has been proposed, which takes advantage of SVA and ISVA, to model the confounding factors and use inferred DVs as covariates. One simulated and nine real datasets, including two different data types (DNA methylation and gene expression), are used for comparing the performance of the algorithms. In the procedure of DVA, the data is firstly correlated with different variables to find the variable that has the smallest P-value as POI using linear regression. The effect of POI is then removed by subtracting the regression part from the data matrix. The dimension of DVs is estimated using a non-local gradient of singular value, indicating the variation of information content, instead of using random matrix theory or PCA as SVA/ISVA do. The threshold, which indicates the non-local gradient of the singular value approaching zero, is set to 0.01. After dimension estimating, the fastICA is used on the residual matrix and the correlated data matrix with the mixed matrix of the residual matrix. The original data matrix is reduced by the indexes of statistically significant correlated features. The reduced data matrix thus contains information with respect to the CFs other than the POI. The fastICA is secondarily applied to this reduced matrix and the vectors of the mixed matrix of the reduced matrix with the best correlation with the mixed matrix of the residual matrix are selected as the decision variables. The results across simulated and real datasets suggest that the proposed method improves the ability to identify the CFs over other methods, and DVA is a more stable algorithm for measurement error or uncertainties. The inferred DVs have the highest correlation with the CFs in all datasets, which indicates that the specific biological signature in each dataset is clearly distinguished by using DVA.

The improvement of modeling confounders by DVA is derived from the optimized combination of SVD and ICA on the residual matrix obtained by removing the effect of POI from the original data. Although the difference between SVA/ISVA and DVA is mainly on the step of dimension estimation, it has been tested that the effectiveness of random matrix theory based on PCA, which ISVA used, relies on the data that has more samples than features. Furthermore, using the same arrangement-based dimension estimation algorithm as SVA, DVA still better identifies the confounding factors because of the blind source separation technique introduced in DVA.

Comparative experiments on simulated data and distinct real datasets suggest that DV is more relevant to the confounder (especially sex) than other surrogates. There is no high correlation ( $R > 0.9$ ) between inferred SVs by SVA and the CFs. In most cases, the correlation values between SVs and CFs are lower than 0.55 ( $P < 0.05$ ). This indicates that there are other mixed components in SV, that is, SV is not pure because SVA is based on irrelevant variables and not on independent variables. In these datasets, age, race, and sex are statistically independent variables, which cause model misspecification of SVA. Indeed, SVA is a tool for discovering a subspace of features or dimension reduction rather than identifying CFs. Since ISVA and DVA both include ICA as a step in data processing, ISVA may identify the CFs by which the subspace spanned may be similar to that spanned by DVs. However, they are different at the dimensionality of the spanned subspaces. ISVA uses random matrix theory based on PCA to estimate dimensionality. PCA is usually conducted by using eigenvalue decomposition of data covariance (or correlation) matrix after normalization. It is worth noting that the rows of the biological data matrix represent features and the columns represent samples when normalizing each feature. PCA cannot achieve reliable results when there are more features than samples. On the contrary, the proposed estimation of dimensionality using a non-local gradient of singular value approximating zero suggests that the cutoff dimensions contain little additional information rather than random matrix theory. Furthermore, since the inferred SVs/ISVs are not pure and even related to the POI because the heterogeneity of data is not uncommon in reality, including SVs/ISVs as covariates in the model disturbs the real biological signal. Although the simulated and real data are highly heterogeneous, the proposed method can achieve better correlations between DVs and CFs. Since DVA models the CFs more accurately and better, it provides a more robust framework for improving the identifiability of CFs for heterogeneous high-dimensional data.

The identified SVs/ISVs by using SVA/ISVA are problematic when they are both correlated with POI and CFs. In this case, distinguishing the effect of the primary signal from that caused by CFs is critical. In our method, the primary biological signal (e.g. related to cancer) of samples can be deemed as a priori knowledge. A priori knowledge correlated with the POI is first removed by solving the residual matrix of the original data matrix. The downstream processing of DVA (e.g. finding significant singular vectors and calculating independent components) is conducted on this residual matrix. This ensures DV is highly related to CF (e.g. sex, $R = 0.91 \sim 0.98$ , $P ≪ 0.05$ ) in more than half of the data sets.

In the case of DNA samples from whole blood, cellar components in the blood determine the interpretation of DNA methylation. Compared to other techniques, the proposed method shows adaptability in identifying different cell types in real whole-blood methylation data with a small sample size. Since the variation of blood cell composition as confounders disturbs the methylation calling, these confounding factors correlated with DVs created by our method can be thought of as covariates in the downstream regression analysis.

In conclusion, the proposed DVA method works well in identifying confounding factors in gene expression and DNA methylation datasets (including whole-blood methylation data). This novel method provides a flexible framework for statistical inference and makes high-latency sample processing more reasonable and accurate. It is expected to be a reliable data mining tool for discovering decisive variables resulting in cancers with known or hidden interference and disturbance.

Footnotes

Acknowledgements

The authors express their gratitude to the anonymous reviewers for the constructive suggestions.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported in part by the National Natural Science Foundation of China under Grants 62301353 and 12272064, in part by Key R&D Program of Hunan (2022SK2104) and Major Project of Basic Science (Natural Science) Research in Colleges and Universities of Jiangsu Province (23KJA580001), in part by Leading plan for scientific and technological innovation of high-tech industries of Hunan (2022GK4010), in part by National Key R&D Program of China (2021YFF0900602).

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

ORCID iD

Jie Xie

References

Zhou

Zhu

. A grouping feature selection method based on feature interaction. Intell Data Anal 2023; 27: 361–377.

Chen

Ryu

Vinyard

, et al. SIMBA: single-cell embedding along with features. Nat Methods 2024; 21: 1003–1013.

Zhao

Ruan

Koestler

, et al. Epigenome-wide scan identifies differentially methylated regions for lung cancer using pre-diagnostic peripheral blood. Epigenetics 2022; 17: 460–472.

Lin

Guo

, et al. Multi-omic characterization of genome-wide abnormal DNA methylation reveals diagnostic and prognostic markers for esophageal squamous-cell carcinoma. Signal Transduction Targeted Ther 2022; 7: 13.

Varghese

Barefoot

Jain

, et al. Integrative analysis of DNA methylation and microrna expression reveals mechanisms of racial heterogeneity in hepatocellular carcinoma. Front Genet 2021; 12: 15.

Luo

Vogeli

. Reclassification of kidney clear cell carcinoma based on immune cell gene-related DNA CpG Pairs. Biomedicines 2021; 9: 24.

Chen

Fujimoto

, et al. Immunogenomic intertumor heterogeneity across primary and metastatic sites in a patient with lung adenocarcinoma. J Exp Clin Cancer Res 2022; 41: 13.

Grant

Wang

Kumari

, et al. Characterising sex differences of autosomal DNA methylation in whole blood using the Illumina EPIC array. Clin Epigenetics 2022; 14: 16.

Nitsche

Vedire

Kannisto

, et al. Visceral obesity in non-small cell lung cancer. Cancers 2022; 14: 3450.

10.

Khomicheva

Vityaev

Ananko

, et al. ExpertDiscovery system application for the hierarchical analysis of eukaryotic transcription regulatory regions based on DNA codes of transcription. Intell Data Anal 2008; 12: 481–494.

11.

Bartish

Abraham

Goncalves

, et al. The role of eiF4F-driven mRNA translation in regulating the tumour microenvironment. Nat Rev Cancer 2023; 23: 408–425.

12.

Chen

Higgins

, et al. DNA methylation differences in noncoding regions in ER-negative breast tumors between black and white women. Front Oncol 2023; 13: 11.

13.

Simonea

Sergioa

Diegob

, et al. An unsupervised clustering approach for leukaemia classification based on DNA micro-arrays data. Intell Data Anal 2007; 11: 175–188.

14.

Salim

Molania

Wang

, et al. RUV-III-NB: normalization of single cell RNA-seq data. Nucleic Acids Res 2022; 50: e96.

15.

González

Guerra

Robles

, et al. CliDaPa: a new approach to combining clinical data with DNA microarrays. Intell Data Anal 2010; 14: 207–223.

16.

Nguyen

Yao

, et al. Associations between DNA methylation and BMI vary by metabolic health status: a potential link to disparate cardiovascular outcomes. Clin Epigenetics 2021; 13: 12.

17.

Jukka

Martti

. Comparing identification methods for DNA investigations of crimes and accidents. Intell Data Anal 2008; 12: 409–423.

18.

Katsaounou

Nicolaou

Vogazianos

, et al. Colon cancer: from epidemiology to prevention. Metabolites 2022; 12: 26.

19.

Boniolo

Dorigatti

Ohnmacht

, et al. Artificial intelligence in early drug discovery enabling precision medicine. Expert Opin Drug Discov 2021; 16: 991–1007.

20.

Chen

, et al. scGMAI: a Gaussian mixture model for clustering single-cell RNA-seq data based on deep autoencoder. Brief Bioinformatics 2021; 22: 1–10.

21.

Dressler

Brägelmann

Reischl

, et al. Normics: proteomic normalization by variance and data-inherent correlation structure. Mol Cell Proteomics MCP 2022; 21: 100269.

22.

Radua

Vieta

Shinohara

, et al. Increased power by harmonizing structural MRI site differences with the ComBat batch method in ENIGMA. Neuroimage 2020; 218: 14.

23.

Zhang

Parmigiani

Johnson

. ComBat-seq: batch effect adjustment for RNA-seq count data. NAR Genomics Bioinf 2020; 2: lqaa078.

24.

Raula

Pedroa

Monicab

. Independent component analysis algorithms for microarray data analysis. Intell Data Anal 2010; 14: 193–206.

25.

Xue

Yazar

Neavin

, et al. Pitfalls and opportunities for applying latent variables in single-cell eQTL analyses. Genome Biol 2023; 24: 11.

26.

Lee

Sun

Wright

, et al. An improved and explicit surrogate variable analysis procedure by coefficient adjustment. Biometrika 2017; 2: 303–316.

27.

Teschendorff

Zhuang

Widschwendter

. Independent surrogate variable analysis to deconvolve confounding factors in large-scale microarray profiling studies. Bioinformatics 2011; 27: 1496–1505.

28.

Ricardoa

Gildardoa

Marco

, et al. Improving pattern classification of DNA microarray data by using PCA and logistic regression. Intell Data Anal 2016; 20: S53–S67.

29.

Hyvarinen

Oja

. Independent component analysis: algorithms and applications. Neural Netw 2000; 13: 411–430.

30.

Plerou

Gopikrishnan

Rosenow

, et al. Random matrix approach to cross correlations in financial data. Phys Rev E 2002; 65: 18.

31.

Zhu

Zhang

, et al. Free vibration of self-powered nanoribbons subjected to thermal-mechanical-electrical fields based on a nonlocal strain gradient theory. Appl Math Model 2022; 110: 583–602.

32.

Jin

Sui

Zhu

, et al. Axial free vibration of rotating FG piezoelectric nano-rods accounting for nonlocal and strain gradient effects. J Vib Eng Technol 2023; 11: 537–549.

33.

McAlpine

Chiu

Nout

, et al. Evaluation of treatment effects in patients with endometrial cancer and POLE mutations: an individual patient data meta-analysis. Cancer 2021; 127: 2409–2422.

34.

Xie

Leung

Chen

, et al. Differential methylation values in differential methylation analysis. Bioinformatics 2019; 35: 1094–1097.

35.

Reho

Saez-Atienzar

Ruffo

, et al. Differential methylation analysis in neuropathologically confirmed dementia with Lewy bodies. Commun Biol 2024; 7.

36.

Lawlor

Marquez

Lee

, et al. V-SVA: an R shiny application for detecting and annotating hidden sources of variation in single-cell RNA-seq data. Bioinformatics 2020; 36: 3582–3584.

37.

Xiao

Jiang

, et al. Genes associated with inflammation for prognosis prediction for clear cell renal cell carcinoma: a multi-database analysis. Transl Cancer Res 2023; 12: 2629–2645.