Network-Constrained Eigen-Single-Cell Profile Estimation for Uncovering Crucial Immunogene Regulatory Systems in Human Bone Marrow

Abstract

We focus on characterizing cell lines from young and aged-healthy and -AML (acute myeloid leukemia) cell lines, and our goal is to identify the key markers associated with the progression of AML. To characterize the age-related phenotypes in AML cell lines, we consider eigenCell analysis that effectively encapsulates the primary expression level patterns across the cell lines. However, earlier investigations utilizing eigenGenes and eigenCells analysis were based on linear combination of all features, leading to the disturbance from noise features. Moreover, the analysis based on a fully dense loading matrix makes it challenging to interpret the results of eigenCells analysis. In order to address these challenges, we develop a novel computational approach termed network-constrained eigenCells profile estimation, which employs a sparse learning strategy. The proposed method estimates eigenCell based on not only the lasso but also network constrained penalization. The use of the network-constrained penalization enables us to simultaneously select neighborhood genes. Furthermore, the hub genes and their regulator/target genes are easily selected as crucial markers for eigenCells estimation. That is, our method can incorporate insights from network biology into the process of sparse loading estimation. Through our methodology, we estimate sparse eigenCells profiles, where only critical markers exhibit expression levels. This allows us to identify the key markers associated with a specific phenotype. Monte Carlo simulations demonstrate the efficacy of our method in reconstructing the sparse structure of eigenCells profiles. We employed our approach to unveil the regulatory system of immunogenes in both young/aged-healthy and -AML cell lines. The markers we have identified for the age-related phenotype in both healthy and AML cell lines have garnered strong support from previous studies. Specifically, our findings, in conjunction with the existing literature, indicate that the activities within this subnetwork of CD79A could be pivotal in elucidating the mechanism driving AML progression, particularly noting the significant role played by the diminished activities in the CD79A subnetwork. We expect that the proposed method will be a useful tool for characterizing disease-related subsets of cell lines, encompassing phenotypes and clones.

1. INTRODUCTION

Acute myeloid leukemia (AML) is a fast-progressing cancer that can lead to rapid fatality if not promptly addressed, particularly affecting people of older age (Jayavelu et al., 2022; Zjablovskaja and Florian, 2019). The dysregulation of immune responses, including innate immune and inflammatory pathways, is strongly associated with the development of abnormalities in hematopoietic stem and progenitor cells, which play a crucial role in the pathogenesis of myelodysplastic syndrome (MDS) and AML (Putnam et al., 2023). Moreover, the immune surveillance mechanisms, which encompass the adaptive and innate immune systems, are inherently crafted to eradicate the development of AML (Taghiloo and Asgarian-Omran, 2021).

Single-cell technologies have considerably advanced, and single-cell RNA-seq widely used in various fields of research. To effectively analyze the single-cell dataset, various computational strategies have been developed and used. Many popular packages, such as, Seurat, Monocle, Scanpy, and orchestrating single-cell analysis, have been proposed for analysis of single-cell RNA-seq (Stuart et al., 2019; Wolf et al., 2018; Trapnell et al., 2014; Amezquita et al., 2020). Dohmen et al. (2022) developed a machine learning pipeline for distinguishing tumor cells from normal cells at the single-cell levels. Furthermore, supervised machine learning strategies, such as, logistic regression, support vector machines, random forests, neural networks, and light gradient boosting machine, have been widely used to classify cell types using single-cell RNA-seq datasets (Le et al., 2022).

In this study, we focus on the immunogene regulatory system in both healthy and AML cell lines because immune dysfunction is important for AML pathogenesis, and understanding of immunobiology of AML is critical for immunotherapy (Khaldoyanidi et al., 2021). Especially, we consider highlighting the unique characteristics specific to young and aged cell lines. To characterize disease-related phenotype, eigenGene and eigenCells analysis has been utilized for various purposes, including the characterization of disease-related phenotypes, gene expression clustering, marker identification, and pathways reconstruction (Mansfeldt et al., 2015; Singh et al., 2020; Alter and Golub, 2005). We employ eigenCells analysis to characterize cell lines of both young/aged-healthy and those with AML. Although the eigenGenes and eigenCells analysis has been widely used, the previous research has predominantly utilized a fully dense loading matrix in singular value decomposition (SVD). However, fully dense loadings may be affected by noise features, as the loading matrix is estimated by a linear combination of all features, which unavoidably includes elements of noise. In addition, the use of fully dense loadings estimated from all features can make it difficult to interpret the results of eigenCell analysis.

To address the issues, we adopt sparse loading estimation in accordance with sparse principal component analysis (sparse PCA) (Zou et al., 2006). Furthermore, we draw attention to a shortcoming in current machine learning methodologies, which are predominantly designed from computational and algorithmic perspectives. This can result in a lack of biologically meaningful interpretations. To ensure biologically reliable results in eigenCells analysis, we introduce a novel computational approach called network-constrained eigenCells profile estimation (nwEigenCELLs). This method incorporates knowledge from network biology, recognizing that genes connected within networks may share similar biological functions, and that hub genes, which are highly connected in the estimated gene network, play a crucial role in gene regulation and biological processes during the sparse eigenCell profile estimation procedure. With the proposed nwEigenCELLs method, we are able to simultaneously select linked genes. Additionally, hub genes and their regulator/target genes tend to have relatively large coefficients and/or are easily selected as crucial markers for eigenCells profile estimation. Consequently, we can interpret the results of eigenCells analysis by focusing on the sub-network (i.e., specific pathway) comprising hub genes, along with their target and regulator genes. It implies that our strategy produces biologically meaningful outcomes.

The effectiveness of our proposed method is illustrated through Monte Carlo simulations. The nwEigenCELLs method demonstrates remarkable performance in recovering the sparse structure of the loading matrix and producing stable results. We apply our strategy to the publicly accessible single-cell gene expression data (Triana et al., 2021) to uncover the immune system traits of young/aged-healthy and -AML cell lines. We use “3 healthy young and 3 healthy old bone marrow donors” and “15 leukemic bone marrow donors” datasets. The immunogenes annotated in the immune system process pathway of Gene Ontology term analysis are extracted. We then uncover eigenCell profiles of the immunogene for four phenotypes and their molecular interplays. The AML immune markers we’ve identified are supported by robust evidence as potential biomarkers for AML. Specifically, the loss of activities of the identified CD79 subnetwork, which includes CD79A, CD79B, CD20, IGHM, IGHD, and IGKC, are likely markers of AML progression. Furthermore, S100A family’s hubness is regarded as a characteristic feature of young cell lines. The results of our scRNA-seq analysis, combined with existing literature, suggest that the reduction in activities within the CD79A subnetwork may play a pivotal role in the progression of AML. Moreover, activators of this subnetwork may offer essential insights for AML prevention.

The remainder of this paper is organized as follows: The methods section introduces the proposed strategy, nwEiegenCELLs. We also introduce computation strategies for gene network estimation and regularization parameter selection in the methods section. Then, we show the results of simulation studies in the Monte Carlo simulations section. Finally, we describe the results of young/aged-healthy and -AML-specific immunogene regulatory system identification. The Discussion section offers the conclusions.

2. NETWORK-CONSTRAINED SINGULAR VALUE DECOMPOSITION

The expression matrix of p genes across n cell lines is represented as $X = (x_{1}, x_{2}, ..., x_{p}) \in ℝ^{n \times p}$ , where $x_{j} = {(x_{1 j}, x_{2 j}, \dots, x_{n j})}^{T}$ for $j = 1, ..., p$ . The SVD of X is given as $X \approx U D V^{T} = \sum_{k = 1}^{K} d_{k} u_{k} v_{k}^{T},$ (1)where $U = [u_{1}, u_{2}, ..., u_{K}] \in ℝ^{n \times K}$ and $V = [v_{1}, v_{2}, ..., v_{K}] \in ℝ^{p \times K}$ are orthogonal matrices, and $D = diag (d_{1}, d_{2}, ..., d_{K})$ with positive singular values $d_{1} \geq d_{2} \geq \dots \geq d_{K}$ on its diagonal. The SVD enables us to perform a linear transformation of the expression levels in X , moving from the space of “n cells” × “p genes” to the reduced space of “K-eigenGenes” × “K-eigenCells”, with $K = min {n, p}$ (Alter et al., 2000). The transformation matrices U and V delineate the expression levels of K-eigenGenes in n cell lines and the expression levels of p genes in K-eigenArrays, respectively (Alter et al., 2000).

We characterize a specific phenotype by the expression levels of eigenCells. This means that the characteristic of a group of cells is described by the eigenCells module, $E C = v_{1} .$ (2)

The EC captures vital information regarding the expression profiles of n cell lines within a specific phenotype. That is, the eigenCells serve as a summary of the principal pattern of expression levels of cell lines, and thus can be considered as pivotal modules for characterizing the phenotype.

Although the SVD is a useful tool to determine the expression levels of eigenGenes and eigenCells that characterize the overall expression of a specific phenotype, previous studies have conducted eigen-Gene and -Cell analysis using the fully dense loading matrix U and V . However, the fully dense loadings exhibit the following limitations:

Disturbance by noise features: The utilization of fully dense loadings in eigenCells analysis may result in erroneous outcomes, as it involves estimating expression levels based on a linear combination of all features, which inevitably includes noisy features.

Difficulty in interpretation: Interpreting the extracted expression levels of all p genes in K-eigenCells is a challenging task.

In order to address this issue, we consider the following sparse PCA to estimate sparse expression levels of eigenCells (Zou et al., 2006). $(\hat{A}, \hat{B}) = \underset{A, B}{arg min} {\sum_{i = 1}^{n} ‖ x_{i} - A B^{T} x_{i} ‖^{2} + λ_{1} \sum_{k = 1}^{K} ‖ β_{k} ‖_{1} + λ_{2} \sum_{k = 1}^{K} ‖ β_{k} ‖_{2}^{2}}, subject to A^{T} A = I_{K},$ (3)where $λ_{1}, λ_{2} > 0, β_{k}$ is a p-dimensional column vector and $B = [β_{1}, β_{2}, ..., β_{K}]$ , and $A = [a_{1}, a_{2}, ..., a_{K}]$ and $a_{k}$ is a p-dimensional vector called k^th principal component loading; that is, the projection of the data ( $X a_{k}$ ) is called the k^th principal component. Then, ${\hat{β}}_{k} \propto v_{k}$ describing the expression levels in k^th eigenCells. Utilizing the sparse loading estimation technique, we are able to derive expression levels of crucial genes in K-eigenCells. In essence, only critical genes display expression profiles within eigenCells. This allows us to efficiently discern the key markers for each phenotype.

To attain more biologically reliable results in sparse eigenCells profile estimation, we develop a novel computational strategy that integrates the following knowledge of network biology into the process of sparse loading estimation. •

Genes that are interconnected in the networks may share similar biological functions.

We encourage similarity in gene coefficients for linked genes when estimating the expression levels of eigenCells. •

Hub genes, which are connected to a multitude of other genes, play a pivotal role in gene regulation and various biological processes (Yu et al., 2017).

In the estimation of the eigenCells profile, we regard hub genes as crucial markers and encourage that the hubgenes are easily selected.

The estimated gene network can be given as a weighted graph $G = (V, E, W)$ . Let $V = {1, ..., p}$ be the set of vertices, each corresponding to one of the p genes. Furthermore, $E \in V \times V$ represents the set of edges, where (i, j) signifies a link between vertices i and j. Lastly, $W = (w_{i j}), (i, j) \in E$ represents the edge weights, and hereafter referred to as an adjacency matrix. From the adjacency matrix W , the normalized Laplacian matrix L for the graph can be computed as follows (Li and Li, 2008; Sun et al., 2014), $L = l_{i j} = {\begin{array}{l} 1 - \frac{w_{i j}}{d_{i}} & if i = j and d_{i} \neq 0, \\ - \frac{w_{i j}}{\sqrt{d_{i} d_{j}}} & if (i, j) \in E, \\ 0 & otherwise \end{array}$ (4)where $d_{i} = \sum_{i \sim j}^{w_{i j}}$ is the degree (i.e., number of edges of the i^th gene or directly connected genes with the i^th gene) of each gene; that is, the large value of d_i indicates that the i^th gene is a hubgene in the gene regulatory network.

We develop a novel machine learning approach, nwEigenCELLs, which integrates the network biology knowledge using network-constrained regularization (Sun et al., 2014), $(\hat{A}, \hat{B}) = \underset{A, B}{arg min} {\sum_{i = 1}^{n} ‖ x_{i} - A B^{T} x_{i} ‖^{2} + λ_{1} \sum_{k = 1}^{K} ‖ β_{r} ‖_{1} + λ_{2} \sum_{k = 1}^{K} β_{k}^{T} L^{k} β_{k}}, = \underset{A, B}{arg min} {\sum_{i = 1}^{n} ‖ x_{i} - A B^{T} x_{i} ‖_{F}^{2} + λ_{1} \sum_{k = 1}^{K} ‖ β_{k} ‖_{1} + λ_{2} \sum_{k = 1}^{K} \sum_{q = 1}^{p} \sum_{j = 1}^{p} {(\frac{sgn (β_{k q}) β_{k q}}{\sqrt{d_{q}}} - \frac{sgn (β_{k j}) β_{k j}}{\sqrt{d_{j}}})}^{2} w_{q j}}, subject to A^{T} A = I_{K},$ (5)where $L^{k} = {S^{k}}^{T} L S^{k}$ with $S^{k} = diag (sgn ({\hat{β}}_{k 1}), ..., sgn ({\hat{β}}_{k p}))$ . Then, ${\hat{β}}_{k} \propto v_{k}$ . Through the last term in the objective function in equation (5), we are able to locally smooth the network and promote the collective selection of linked genes. It was pointed out that the network constrained penalty; that is, $\sum_{k = 1}^{K} β_{k}^{T} L β_{k}$ cannot perform well when neighboring genes have opposite signs of regression coefficients, because the penalty lead to similar values of coefficient of neighboring genes. On the other hand, the sign adjusted network constrained penalty used in our method; that is, $\sum_{k = 1}^{K} β_{k}^{T} L^{k} β_{k} = \sum_{q = 1}^{p} \sum_{j = 1}^{p} {(\frac{sgn (β_{k q}) β_{k q}}{\sqrt{d_{q}}} - \frac{sgn (β_{k j}) β_{k j}}{\sqrt{d_{j}}})}^{2} w_{q j}$ can effectively locally smooth the network and encourage the simultaneous selection of neighboring genes, even though neighboring genes have opposite coefficient signs. Furthermore, a relatively small penalty is imposed to the hub genes by scaling the coefficients using the square root of their degrees. Consequently, the hub genes and their regulator/target genes tend to have relatively large coefficients and are more easily identified as essential markers for eigenCells profile estimation. In summary, our method can integrate knowledge from network biology and thus we can estimate the expression levels of eigenCells using crucial subnetworks composed of hub genes and their regulator and target genes. Consequently, we can interpret the results of eigenCell analysis based on the genes within a specific pathway, yielding biologically reliable interpretations.

2.1. Implementation

The algorithm of nwEigenCELLs is given as follows:

Step 1 Compute the adjacency matrix and the normalized Laplacian matrix.

Step 2 Start A at $V = [v_{1}, v_{2}, ..., v_{K}]$ , which is from ordinary SVD of X .

Step 3 Given a fixed $A = [a_{1}, a_{2}, ..., a_{K}]$ , estimate $β_{k}, k = 1, 2, ..., K$ , ${\hat{β}}_{k} = \underset{β_{k}}{arg min} {‖ z_{k} - X β_{k} ‖^{2}} + λ_{1} ‖ β_{k} ‖_{1} + λ_{2} \sum_{q = 1}^{p} \sum_{j = 1}^{p} {(\frac{sgn (β_{k q}) β_{k q}}{\sqrt{d_{q}}} - \frac{sgn (β_{k j}) β_{k j}}{\sqrt{d_{j}}})}^{2} w_{q j}$ where $z_{k} = X a_{k}$ . Update $\hat{B} = [{\hat{β}}_{1}, {\hat{β}}_{2}, ..., {\hat{β}}_{K}]$ .

Step 4 For a fixed $\hat{B}$ , compute the SVD of $X^{T} X \hat{B} = U Γ V^{T}$ and update $\hat{A} = U V^{T}$ . The solution is based on Theorem 4 of Zou et al. (2006).

Step 5 Repeat Steps 3–4, until convergence with the convergence condition $m a x_{j} | {\hat{β}}_{k j}^{(r)} - {\hat{β}}_{k j}^{(r + 1)} | < τ$ , where ${\hat{β}}_{k j}^{(r)}$ is the r^th iteration of the estimated parameter β_kj.

Step 6 The expression levels of eigenCells are given by scale ${\hat{β}}_{k}$ to have unit l₂-norm, that is, $v_{k} = \frac{{\hat{β}}_{k}}{‖‖ {\hat{β}}_{k} ‖‖}$ for $k = 1, ..., K$ .

The input and output of the algorithm are expression levels of all genes and expression levels of eigenCells, respectively.

2.2. Regularization Parameter Selection

Selecting the appropriate regularization parameters is a critical issue of sparse eigenCell analysis, as the outcomes of the proposed nwEigenCELLs heavily depend on the parameters λ₁ and λ₂. We employ the Bayesian Information Criterion presented in Guo et al. (2010) to choose the appropriate regularization parameters. $BIC = ‖ X - XBA ‖_{F}^{2} / σ^{2} + log (n) d f$ (6)where $σ^{2} = \frac{1}{n} | | X - X \hat{A} {\hat{A}}^{T} {| |}_{F}^{2}$ and df is the degree of freedom of the nwEigenCELLs. We employ the number of non-zero elements in B as an estimate for the degrees of freedom in line with Zou et al. (2006) and Guo et al. (2010).

3. MONTE CARLO Simulations

We employ Monte Carlo simulations to investigate the performance of the proposed nwEigenCELLs. The simulation scenarios are designed to be consistent with the benchmarks established in previous studies (Zou et al., 2006; Guo et al., 2010). We evaluate that the methods can identify crucial features in loading estimation, that is, true positive and true negative in sparse loading estimation.

3.1. Simulation 1

Scenario 1: The hidden vectors $V_{1}, V_{2}, V_{3}$ for representing principal components are generated using the following procedure: $\begin{array}{l} V_{1} \sim N (0, 250) \\ V_{2} \sim N (0, 350) \\ V_{3} = - 0.3 V_{1} + 0.65 V_{2} + ε_{V_{3}}, ε_{V_{3}} \sim N (0, 1) \end{array}$

The n × p data matrix $X = (x_{1}, x_{2}, ..., x_{p})$ is subsequently generated using the following process: $\begin{array}{l} x_{i j} = V_{i 1} + e_{i j}, & for 1 \leq j < 5, \\ x_{i j} = V_{i 2} + e_{i j}, & for 5 \leq j < 9, \\ x_{i j} = V_{i 3} + e_{i j}, & for j = 9, 10. \end{array}$ where $\begin{array}{l} e_{i j} \sim N (0.7 \times e_{i 1}, σ_{e}^{2}) & for 2 < j \leq 4, \\ e_{i j} \sim N (0.7 \times e_{i 5}, σ_{e}^{2}) & for 6 < j \leq 8, \\ e_{i j} \sim (0, 1) & for 8 < j \leq 10, \end{array}$ where $e_{i j} \sim N (0, 1)$ for j = 1, 5. The features $x_{1}, x_{2}, x_{3}, x_{4}$ form the derived first component (i.e., V₂), while the features $x_{5}, x_{6}, x_{7}, x_{8}$ and $x_{9}, x_{10}$ form the derived second and third components (i.e., V₁ and V₃), respectively. Ideally, the features $x_{1}, x_{2}, x_{3}, x_{4}$ have non-zero loadings for the first eigenvalue, while the features $x_{5}, x_{6}, x_{7}, x_{8}$ have non-zero loadings for the second eigenvalue. We suppose that the $1^{s t}$ and $5^{t h}$ genes are transcription factors regulating $2^{n d} \sim 4^{t h}$ genes and $6^{t h} \sim 8^{t h}$ genes, respectively.

Scenario 2: The hidden vectors $V_{1}, V_{2}, V_{3}$ for representing principal components are generated using the following procedure: $\begin{array}{l} V_{1} \sim N (0, 350) \\ V_{2} \sim N (0, 250) \\ V_{3} = - 0.3 V_{1} + 0.65 V_{2} + ε_{V_{3}}, ε_{V_{3}} \sim N (0, 1) \end{array}$

The n × p data matrix $X = (x_{1}, x_{2}, ..., x_{p})$ is subsequently generated using the following process: $\begin{array}{l} x_{i j} = V_{i 1} + e_{i j}, & for 1 \leq j < 6, \\ x_{i j} = V_{i 2} + e_{i j}, & for 6 \leq j < 9, \\ x_{i j} = V_{i 3} + e_{i j}, & for j = 9, 10. \end{array}$ where $\begin{array}{l} e_{i j} = N (0.7 \times e_{i 1}, σ_{e}^{2}) & for 1 < j \leq 5, \\ e_{i j} = N (0.7 \times e_{i 6}, σ_{e}^{2}) & for 6 < j \leq 8, \\ e_{i j} = N (0, 1) & for 8 < j \leq 10, \end{array}$ where $e_{i j} \sim N (0, 1)$ for j = 1, 6. The features $x_{1}, x_{2}, x_{3}, x_{4}, x_{5}$ form the derived first component (i.e., V₂), while the features $x_{6}, x_{7}, x_{8}$ and x₉, x₁₀ form the derived second and third components (i.e., V₁ and V₃), respectively. Ideally, the features $x_{1}, x_{2}, x_{3}, x_{4}, x_{5}$ have non-zero loadings for the first eigenvalue, while the features $x_{6}, x_{7}, x_{8}$ have non-zero loadings for the second eigenvalue. We suppose that the $1^{s t}$ and $6^{t h}$ genes are transcription factors regulating $2^{n d} \sim 5^{t h}$ genes and $7^{t h} \sim 8^{t h}$ genes, respectively.

We set n = 20, and the simulations are replicated 500 times. In simulation study, we first compute covariance of X to describe the gene regulatory network. Our method (nwEigenCELLs) is evaluated through a comparative analysis with results from both ordinary sparse PCA (SPCA) and fused loading-based sparse PCA (fSPCA) (Guo et al., 2010). Tables 1 and 2 show the simulation outcomes for both $σ_{e}^{2} = 0.51, 1$ and $σ_{e}^{2} = 3, 5$ , respectively, specifically focusing on the medians of the first two loading vectors. The median of loading values indicate accuracy of loading estimation in the viewpoint of true positive and true negative. The values in parentheses represent the median absolute deviation of the estimated loadings across 500 replications, where the bold numbers indicate more stable results among the SPCA and nwEigenCELLs. As shown in Table 1, the nwEigenCELLs method successfully reconstructs the sparse structure of eigenCell expression levels. On the other hand, the fused loading-based sparse PCA exhibits limitations in recovering these sparse structures. While the ordinary sparse PCA does exhibit effectiveness in identifying the sparse structure of the loading matrix, it tends to yield less stable results compared to our proposed method, in overall. This is evidenced by the relatively high median absolute deviation of the estimated loadings observed in the ordinary sparse PCA (see the values in parentheses). Furthermore, the ordinary sparse PCA cannot perform well for reconstructing the sparse structure for large values of noise, that is, $σ_{e}^{2} = 3, 5$ as shown in Table 2. In scenario 1, the ordinary sparse PCA cannot consistently reconstruct sparse structure of loading of PC2, while the nwEigenCELLs shows effective performance (see the values in parentheses). Especially, the performance of the ordinary sparse PCA becomes worse in $σ_{e}^{2} = 5$ , that is, the method cannot identify the crucial variables for the loading vector estimation of PC1. To conclude, the nwEigenCELLs method stands out as an efficient tool for uncovering the sparse expression levels of eigenCells.

Table 1.
Simulation Results

$σ_{e}^{2}$ Loadings Scenario 1 Scenario 2

PC1 PC2 PC1 PC2

PCA SPCA FSPCA nwEigen PCA SPCA FSPCA nwEigen PCA SPCA FSPCA nwEigen PCA SPCA FSPCA nwEigen

0.51 1 0.001 (0.16) 0 (0) −0.005 (0.01) 0 (0) −0.246 (0.25) −0.371 (0.25) −0.500 (0) −0.397 (0.22) −0.235 (0.21) −0.444 (0.07) −0.447 (0) −0.426 (0.07) 0.007 (0.11) 0 (0) 0 (0) 0 (0)

2 −0.004 (0.15) 0 (0) −0.005 (0.01) 0 (0) −0.240 (0.26) −0.362 (0.22) −0.500 (0) −0.390 (0.19) −0.238 (0.21) −0.418 (0.07) −0.447 (0) −0.427 (0.07) 0.010 (0.11) 0 (0) 0 (0) 0 (0)

3 0.000 (0.15) 0 (0) −0.005 (0.01) 0 (0) −0.244 (0.25) −0.382 (0.22) −0.500 (0) −0.399 (0.19) −0.232 (0.21) −0.436 (0.07) −0.447 (0) −0.436 (0.07) 0.007 (0.11) 0 (0) 0 (0) 0 (0)

4 −0.001 (0.15) 0 (0) −0.005 (0.01) 0 (0) −0.241 (0.25) −0.396 (0.22) −0.500 (0) −0.406 (0.20) −0.236 (0.21) −0.426 (0.06) −0.447 (0) −0.424 (0.07) 0.010 (0.11) 0 (0) 0 (0) 0 (0)

5 0.096 (0.37) 0.289 (0.29) 0.361 (0.12) 0.291 (0.29) −0.012 (0.15) 0 (0) 0 (0) 0 (0) −0.237 (0.21) −0.422 (0.07) −0.447 (0) −0.433 (0.07) 0.006 (0.11) 0 (0) 0 (0) 0 (0)

6 0.100 (0.36) 0.288 (0.29) 0.361 (0.12) 0.324 (0.32) −0.011 (0.15) 0 (0) 0 (0) 0 (0) −0.006 (0.12) 0 (0) 0 (0) 0 (0) −0.053 (0.47) −0.233 (0.34) −0.006 (0.54) −0.318 (0.32)

7 0.093 (0.37) 0.310 (0.31) 0.361 (0.12) 0.333 (0.29) −0.010 (0.15) 0 (0) 0 (0) 0 (0) −0.007 (0.12) 0 (0) 0 (0) 0 (0) −0.047 (0.47) −0.242 (0.34) −0.006 (0.54) −0.306 (0.31)

8 0.098 (0.36) 0.300 (0.30) 0.361 (0.12) 0.312 (0.31) −0.006 (0.15) 0 (0) 0 (0) 0 (0) −0.004 (0.12) 0 (0) 0 (0) 0 (0) −0.051 (0.47) −0.329 (0.33) −0.006 (0.54) −0.337 (0.34)

9 0.048 (0.25) 0 (0) 0.173 (0.17) 0 (0) 0.012 (0.12) 0 (0) 0 (0) 0 (0) 0.000 (0.16) 0 (0) 0 (0) 0 (0) 0.009 (0.28) 0 (0) −0.003 (0.16) 0 (0)

10 0.053 (0.25) 0 (0) 0.174 (0.17) 0 (0) 0.011 (0.12) 0 (0) 0 (0) 0 (0) 0.001 (0.16) 0 (0) 0 (0) 0 (0) 0.001 (0.28) 0 (0) −0.003 (0.15) 0 (0)

1 1 0.001 (0.16) 0 (0) −0.005 (0.01) 0 (0) −0.246 (0.25) −0.371 (0.25) −0.500 (0) −0.397 (0.22) −0.235 (0.21) −0.444 (0.07) −0.447 (0) −0.426 (0.07) 0.007 (0.11) 0 (0) 0 (0) 0 (0)

2 −0.004 (0.15) 0 (0) −0.005 (0.01) 0 (0) −0.240 (0.26) −0.362 (0.22) −0.500 (0) −0.390 (0.19) −0.238 (0.21) −0.418 (0.07) −0.447 (0) −0.427 (0.07) 0.010 (0.11) 0 (0) 0 (0) 0 (0)

3 0.000 (0.15) 0 (0) −0.005 (0.01) 0 (0) −0.244 (0.25) −0.382 (0.22) −0.500 (0) −0.399 (0.19) −0.232 (0.21) −0.436 (0.07) −0.447 (0) −0.436 (0.07) 0.007 (0.11) 0 (0) 0 (0) 0 (0)

4 −0.001 (0.15) 0 (0) −0.005 (0.01) 0 (0) −0.241 (0.25) −0.396 (0.22) −0.500 (0) −0.406 (0.20) −0.236 (0.21) −0.426 (0.06) −0.447 (0) −0.424 (0.07) 0.010 (0.11) 0 (0) 0 (0) 0 (0)

5 0.096 (0.37) 0.289 (0.29) 0.361 (0.12) 0.291 (0.29) −0.012 (0.15) 0 (0) 0 (0) 0 (0) −0.237 (0.21) −0.422 (0.07) −0.447 (0) −0.433 (0.07) 0.006 (0.11) 0 (0) 0 (0) 0 (0)

6 0.100 (0.36) 0.288 (0.29) 0.361 (0.12) 0.324 (0.32) −0.011 (0.15) 0 (0) 0 (0) 0 (0) −0.006 (0.12) 0 (0) 0 (0) 0 (0) −0.053 (0.47) −0.233 (0.34) −0.006 (0.54) −0.318 (0.32)

7 0.093 (0.37) 0.310 (0.31) 0.361 (0.12) 0.333 (0.29) −0.010 (0.15) 0 (0) 0 (0) 0 (0) −0.007 (0.12) 0 (0) 0 (0) 0 (0) −0.047 (0.47) −0.242 (0.34) −0.006 (0.54) −0.306 (0.31)

8 0.098 (0.36) 0.300 (0.30) 0.361 (0.12) 0.312 (0.31) −0.006 (0.15) 0 (0) 0 (0) 0 (0) −0.004 (0.12) 0 (0) 0 (0) 0 (0) −0.051 (0.47) −0.329 (0.33) −0.006 (0.54) −0.337 (0.34)

9 0.048 (0.25) 0 (0) 0.173 (0.17) 0 (0) 0.012 (0.12) 0 (0) 0 (0) 0 (0) 0.000 (0.16) 0 (0) 0 (0) 0 (0) 0.009 (0.28) 0 (0) −0.003 (0.16) 0 (0)

10 0.053 (0.25) 0 (0) 0.174 (0.17) 0 (0) 0.011 (0.12) 0 (0) 0 (0) 0 (0) 0.001 (0.16) 0 (0) 0 (0) 0 (0) 0.001 (0.28) 0 (0) −0.003 (0.15) 0 (0)

$σ_{e}^{2}$	Loadings	Scenario 1	Scenario 2
0.51	1	0.001 (0.16)	0 (0)	−0.005 (0.01)	0 (0)	−0.246 (0.25)	−0.371 (0.25)	−0.500 (0)	−0.397 (0.22)	−0.235 (0.21)	−0.444 (0.07)	−0.447 (0)	−0.426 (0.07)	0.007 (0.11)	0 (0)	0 (0)	0 (0)
2	−0.004 (0.15)	0 (0)	−0.005 (0.01)	0 (0)	−0.240 (0.26)	−0.362 (0.22)	−0.500 (0)	−0.390 (0.19)	−0.238 (0.21)	−0.418 (0.07)	−0.447 (0)	−0.427 (0.07)	0.010 (0.11)	0 (0)	0 (0)	0 (0)
3	0.000 (0.15)	0 (0)	−0.005 (0.01)	0 (0)	−0.244 (0.25)	−0.382 (0.22)	−0.500 (0)	−0.399 (0.19)	−0.232 (0.21)	−0.436 (0.07)	−0.447 (0)	−0.436 (0.07)	0.007 (0.11)	0 (0)	0 (0)	0 (0)
4	−0.001 (0.15)	0 (0)	−0.005 (0.01)	0 (0)	−0.241 (0.25)	−0.396 (0.22)	−0.500 (0)	−0.406 (0.20)	−0.236 (0.21)	−0.426 (0.06)	−0.447 (0)	−0.424 (0.07)	0.010 (0.11)	0 (0)	0 (0)	0 (0)
5	0.096 (0.37)	0.289 (0.29)	0.361 (0.12)	0.291 (0.29)	−0.012 (0.15)	0 (0)	0 (0)	0 (0)	−0.237 (0.21)	−0.422 (0.07)	−0.447 (0)	−0.433 (0.07)	0.006 (0.11)	0 (0)	0 (0)	0 (0)
6	0.100 (0.36)	0.288 (0.29)	0.361 (0.12)	0.324 (0.32)	−0.011 (0.15)	0 (0)	0 (0)	0 (0)	−0.006 (0.12)	0 (0)	0 (0)	0 (0)	−0.053 (0.47)	−0.233 (0.34)	−0.006 (0.54)	−0.318 (0.32)
7	0.093 (0.37)	0.310 (0.31)	0.361 (0.12)	0.333 (0.29)	−0.010 (0.15)	0 (0)	0 (0)	0 (0)	−0.007 (0.12)	0 (0)	0 (0)	0 (0)	−0.047 (0.47)	−0.242 (0.34)	−0.006 (0.54)	−0.306 (0.31)
8	0.098 (0.36)	0.300 (0.30)	0.361 (0.12)	0.312 (0.31)	−0.006 (0.15)	0 (0)	0 (0)	0 (0)	−0.004 (0.12)	0 (0)	0 (0)	0 (0)	−0.051 (0.47)	−0.329 (0.33)	−0.006 (0.54)	−0.337 (0.34)
9	0.048 (0.25)	0 (0)	0.173 (0.17)	0 (0)	0.012 (0.12)	0 (0)	0 (0)	0 (0)	0.000 (0.16)	0 (0)	0 (0)	0 (0)	0.009 (0.28)	0 (0)	−0.003 (0.16)	0 (0)
10	0.053 (0.25)	0 (0)	0.174 (0.17)	0 (0)	0.011 (0.12)	0 (0)	0 (0)	0 (0)	0.001 (0.16)	0 (0)	0 (0)	0 (0)	0.001 (0.28)	0 (0)	−0.003 (0.15)	0 (0)
1	1	0.001 (0.16)	0 (0)	−0.005 (0.01)	0 (0)	−0.246 (0.25)	−0.371 (0.25)	−0.500 (0)	−0.397 (0.22)	−0.235 (0.21)	−0.444 (0.07)	−0.447 (0)	−0.426 (0.07)	0.007 (0.11)	0 (0)	0 (0)	0 (0)
2	−0.004 (0.15)	0 (0)	−0.005 (0.01)	0 (0)	−0.240 (0.26)	−0.362 (0.22)	−0.500 (0)	−0.390 (0.19)	−0.238 (0.21)	−0.418 (0.07)	−0.447 (0)	−0.427 (0.07)	0.010 (0.11)	0 (0)	0 (0)	0 (0)
3	0.000 (0.15)	0 (0)	−0.005 (0.01)	0 (0)	−0.244 (0.25)	−0.382 (0.22)	−0.500 (0)	−0.399 (0.19)	−0.232 (0.21)	−0.436 (0.07)	−0.447 (0)	−0.436 (0.07)	0.007 (0.11)	0 (0)	0 (0)	0 (0)
4	−0.001 (0.15)	0 (0)	−0.005 (0.01)	0 (0)	−0.241 (0.25)	−0.396 (0.22)	−0.500 (0)	−0.406 (0.20)	−0.236 (0.21)	−0.426 (0.06)	−0.447 (0)	−0.424 (0.07)	0.010 (0.11)	0 (0)	0 (0)	0 (0)
5	0.096 (0.37)	0.289 (0.29)	0.361 (0.12)	0.291 (0.29)	−0.012 (0.15)	0 (0)	0 (0)	0 (0)	−0.237 (0.21)	−0.422 (0.07)	−0.447 (0)	−0.433 (0.07)	0.006 (0.11)	0 (0)	0 (0)	0 (0)
6	0.100 (0.36)	0.288 (0.29)	0.361 (0.12)	0.324 (0.32)	−0.011 (0.15)	0 (0)	0 (0)	0 (0)	−0.006 (0.12)	0 (0)	0 (0)	0 (0)	−0.053 (0.47)	−0.233 (0.34)	−0.006 (0.54)	−0.318 (0.32)
7	0.093 (0.37)	0.310 (0.31)	0.361 (0.12)	0.333 (0.29)	−0.010 (0.15)	0 (0)	0 (0)	0 (0)	−0.007 (0.12)	0 (0)	0 (0)	0 (0)	−0.047 (0.47)	−0.242 (0.34)	−0.006 (0.54)	−0.306 (0.31)
8	0.098 (0.36)	0.300 (0.30)	0.361 (0.12)	0.312 (0.31)	−0.006 (0.15)	0 (0)	0 (0)	0 (0)	−0.004 (0.12)	0 (0)	0 (0)	0 (0)	−0.051 (0.47)	−0.329 (0.33)	−0.006 (0.54)	−0.337 (0.34)
9	0.048 (0.25)	0 (0)	0.173 (0.17)	0 (0)	0.012 (0.12)	0 (0)	0 (0)	0 (0)	0.000 (0.16)	0 (0)	0 (0)	0 (0)	0.009 (0.28)	0 (0)	−0.003 (0.16)	0 (0)
10	0.053 (0.25)	0 (0)	0.174 (0.17)	0 (0)	0.011 (0.12)	0 (0)	0 (0)	0 (0)	0.001 (0.16)	0 (0)	0 (0)	0 (0)	0.001 (0.28)	0 (0)	−0.003 (0.15)	0 (0)

Median values of the first two loading vectors, where the values in parentheses is the median absolute deviation over 500 replications ( $σ_{e}^{2} = 0.51$ and $σ_{e}^{2} = 0.1$ ).

FSPCA, fused loading-based sparse PCA; PCA, principal component analysis; SPCA, ordinary sparse PCA.

Table 2.

Simulation Results

$σ_{e}^{2}$	Loadings	Scenario 1								Scenario 2
		PC1				PC2				PC1				PC2
		PCA	SPCA	FSPCA	nwEigen	PCA	SPCA	FSPCA	nwEigen	PCA	SPCA	FSPCA	nwEigen	PCA	SPCA	FSPCA	nwEigen
3	1	0.015 (0.15)	0 (0)	−0.005 (0.01)	0 (0)	−0.052 (0.43)	−0.352 (0.26)	−0.500 (0)	−0.364 (0.22)	−0.111 (0.33)	−0.435 (0.06)	−0.447 (0)	−0.433 (0.06)	−0.014 (0.12)	0 (0)	0 (0)	0 (0)
	2	0.016 (0.15)	0 (0)	−0.005 (0.01)	0 (0)	−0.053 (0.43)	−0.373 (0.25)	−0.500 (0)	−0.392 (0.22)	−0.112 (0.33)	−0.425 (0.06)	−0.447 (0)	−0.423 (0.06)	−0.014 (0.12)	0 (0)	0 (0)	0 (0)
	3	0.013 (0.15)	0 (0)	−0.005 (0.01)	0 (0)	−0.054 (0.43)	−0.368 (0.24)	−0.500 (0)	−0.387 (0.21)	−0.118 (0.33)	−0.441 (0.07)	−0.447 (0)	−0.436 (0.07)	−0.015 (0.12)	0 (0)	0 (0)	0 (0)
	4	0.013 (0.15)	0 (0)	−0.005 (0.01)	0 (0)	−0.054 (0.43)	−0.380 (0.25)	−0.500 (0)	−0.397 (0.22)	−0.108 (0.34)	−0.423 (0.07)	−0.447 (0)	−0.424 (0.06)	−0.015 (0.12)	0 (0)	0 (0)	0 (0)
	5	0.061 (0.40)	0.319 (0.30)	0.382 (0.10)	0.303 (0.30)	0.009 (0)	0 (0.16)	0 (0)	0 (0)	−0.108 (0.34)	−0.430 (0.06)	−0.447 (0)	−0.436 (0.05)	−0.016 (0.12)	0 (0)	0 (0)	0 (0)
	5	0.067 (0.39)	0.335 (0.31)	0.382 (0.10)	0.332 (0.30)	0.008 0.00	0 (0.16)	0 (0)	0 (0)	0.003 (0.13)	0 (0)	0 (0.01)	0 (0)	0.201 (0.33)	−0.319 (0.32)	−0.426 (0.15)	−0.380 (0.36)
	7	0.059 (0.40)	0.302 (0.30)	0.382 (0.10)	0.299 (0.30)	0.012 0.00	0 (0.16)	0 (0)	0 (0)	0.004 (0.12)	0 (0)	0 (0.01)	0 (0)	0.197 (0.33)	−0.292 (0.29)	−0.426 (0.15)	−0.347 (0.35)
	8	0.065 (0.39)	0.323 (0.32)	0.382 (0.10)	0.338 (0.31)	0.004 0.00	0 (0.16)	0 (0)	0 (0)	0.002 (0.13)	0 (0)	0 (0.01)	0 (0)	0.191 (0.34)	−0.367 (0.37)	−0.426 (0.15)	−0.388 (0.37)
	9	0.016 (0.28)	0 (0)	0.183 (0.17)	0 (0)	0.008 0.00	0 (0.11)	0 (0)	0 (0)	−0.004 (0.17)	0 (0)	0.010 (0)	0 (0)	0.103 (0.23)	0 (0)	−0.062 (0.16)	0 (0)
	10	0.020 (0.28)	0 (0)	0.183 (0.17)	0 (0)	0.009 0.00	0 (0.12)	0 (0)	0 (0)	−0.006 (0.17)	0 (0)	0.010 (0)	0 (0)	0.093 (0.24)	0 (0)	−0.062 (0.16)	0 (0)
5	1	−0.005 (0.16)	0 (0)	−0.006 (0.01)	0 (0)	−0.005 (0.46)	−0.382 (0.24)	−0.500 (0)	−0.396 (0.20)	0.298 (0.15)	−0.437 (0.06)	−0.447 (0)	−0.434 (0.06)	0.005 (0.12)	0 (0)	0 (0)	0 (0)
	2	−0.005 (0.16)	0 (0)	−0.006 (0.01)	0 (0)	−0.005 (0.46)	−0.364 (0.26)	−0.500 (0)	−0.386 (0.21)	0.296 (0.15)	−0.435 (0.07)	−0.447 (0)	−0.436 (0.06)	0.007 (0.12)	0 (0)	0 (0)	0 (0)
	3	0.000 (0.16)	0 (0)	−0.006 (0.01)	0 (0)	0.001 (0.46)	−0.388 (0.24)	−0.500 (0)	−0.391 (0.21)	0.300 (0.15)	−0.437 (0.07)	−0.447 (0)	−0.427 (0.06)	0.004 (0.12)	0 (0)	0 (0)	0 (0)
	4	−0.002 (0.16)	0 (0)	−0.006 (0.01)	0 (0)	0.005 (0.46)	−0.364 (0.25)	−0.500 (0)	−0.387 (0.22)	0.298 (0.15)	−0.428 (0.07)	−0.447 (0)	−0.430 (0.06)	0.007 (0.12)	0 (0)	0 (0)	0 (0)
	5	0.061 (0.40)	0 (0.48)	0.294 (0.19)	0.014 (0.47)	0.019 (0.17)	0 (0)	0 (0)	0 (0)	0.306 (0.14)	−0.425 (0.07)	−0.447 (0)	−0.427 (0.07)	0.004 (0.12)	0 (0)	0 (0)	0 (0)
	6	0.060 (0.40)	0 (0.48)	0.294 (0.19)	0.033 (0.47)	0.015 (0.17)	0 (0)	0 (0)	0 (0)	−0.014 (0.13)	0 (0)	0 (0.01)	0 (0)	0.325 (0.20)	−0.346 (0.35)	−0.478 (0.09)	−0.382 (0.35)
	7	0.061 (0.40)	0 (0.48)	0.294 (0.19)	0.079 (0.45)	0.019 (0.16)	0 (0)	0 (0)	0 (0)	−0.013 (0.13)	0 (0)	0 (0.01)	0 (0)	0.332 (0.20)	−0.336 (0.34)	−0.478 (0.09)	−0.379 (0.37)
	8	0.061 (0.40)	0 (0.48)	0.294 (0.19)	0.033 (0.46)	0.019 (0.16)	0 (0)	0 (0)	0 (0)	−0.014 (0.13)	0 (0)	0 (0.01)	0 (0)	0.324 (0.20)	−0.334 (0.33)	−0.478 (0.09)	−0.378 (0.33)
	9	0.016 (0.27)	0 (0)	0.167 (0.16)	0 (0)	0.004 (0.12)	0 (0)	0 (0)	0 (0)	−0.046 (0.16)	0 (0)	0.009 (0)	0 (0)	0.148 (0.19)	0 (0)	−0.088 (0.13)	0 (0)
	10	0.030 (0.26)	0 (0)	0.167 (0.16)	0 (0)	0.009 (0.12)	0 (0)	0 (0)	0 (0)	−0.047 (0.16)	0 (0)	0.009 (0)	0 (0)	0.139 (0.19)	0 (0)	−0.088 (0.13)	0 (0)

Median values of the first two loading vectors, where the values in parentheses is the median absolute deviation over 500 replications ( $σ_{e}^{2} = 3$ and $σ_{e}^{2} = 5$ ).

FSPCA, fused loading-based sparse PCA; PCA, principal component analysis; SPCA, ordinary sparse PCA.

3.2. Simulation 2

We also consider the simulation situation with the number of features p = 50 (i.e., p > n).

Scenario 1: The n × p data matrix $X = (x_{1}, x_{2}, ..., x_{p})$ is subsequently generated using the following process: $\begin{array}{l} x_{i j} = V_{i 1} + e_{i j}, & for 1 \leq j \leq 20, \\ x_{i j} = V_{i 2} + e_{i j}, & for 20 < j \leq 40, \\ x_{i j} = V_{i 3} + e_{i j}, & for 40 < j \leq 50. \end{array}$ where $e_{i j} \sim N (0, 1)$ for $j = 1, 11, 21, 31$ , and $\begin{array}{l} e_{i j} = N (0.7 \times e_{i 1}, σ_{e}^{2}) & for 1 \leq j \leq 10, \\ e_{i j} = N (0.7 \times e_{i 11}, σ_{e}^{2}) & for 12 < j \leq 20, \\ e_{i j} = N (0.7 \times e_{i 21}, σ_{e}^{2}) & for 22 < j \leq 30, \\ e_{i j} = N (0.7 \times e_{i 31}, σ_{e}^{2}) & for 32 < j \leq 40. \end{array}$

The features $x_{1} \sim x_{20}$ form the derived first component (i.e., V₂), while the features $x_{21} \sim x_{40}$ and $x_{41} \sim x_{50}$ form the derived second and third components (i.e., V₁ and V₃), respectively. Ideally, the features $x_{1} \sim x_{20}$ have non-zero loadings for the first eigenvalue, while the features $x_{21} \sim x_{40}$ have non-zero loadings for the second eigenvalue. We suppose that the $1^{s t}, 11^{t h}, 21^{t h}, 31^{t h}$ genes are transcription factors regulating $2^{n d} \sim 10^{t h}$ genes, $12^{t h} \sim 20^{t h}$ genes, $22^{t h} \sim 30^{t h}$ genes and $32^{t h} \sim 40^{t h}$ genes, respectively.

Scenario 2: The n × p data matrix $X = (x_{1}, x_{2}, ..., x_{p})$ is subsequently generated using the following process: $\begin{array}{l} x_{i j} = V_{i 1} + e_{i j}, & for 1 \leq j \leq 30, \\ x_{i j} = V_{i 2} + e_{i j}, & for 30 < j \leq 40, \\ x_{i j} = V_{i 3} + e_{i j}, & for 40 < j \leq 50. \end{array}$ where $e_{i j} \sim N (0, 1)$ for $j = 1, 11, 21, 31$ , and $\begin{array}{l} e_{i j} = N (0.7 \times e_{i 1}, σ_{e}^{2}) & for 1 \leq j \leq 10, \\ e_{i j} = N (0.7 \times e_{i 11}, σ_{e}^{2}) & for 12 < j \leq 20, \\ e_{i j} = N (0.7 \times e_{i 21}, σ_{e}^{2}) & for 22 < j \leq 30, \\ e_{i j} = N (0.7 \times e_{i 31}, σ_{e}^{2}) & for 32 < j \leq 40. \end{array}$

The features $x_{1} \sim x_{30}$ form the derived first component (i.e., V₂), while the features $x_{31} \sim x_{40}$ and $x_{41} \sim x_{50}$ form the derived second and third components (i.e., V₁ and V₃), respectively. Ideally, the features $x_{1} \sim x_{30}$ have non-zero loadings for the first eigenvalue, while the features $x_{31} \sim x_{40}$ have non-zero loadings for the second eigenvalue. We suppose that the $1^{s t}, 11^{t h}, 21^{t h}$ and $31^{t h}$ genes are transcription factors regulating $2^{n d} \sim 10^{t h}$ genes, $12^{t h} \sim 20^{t h}$ genes, $22^{t h} \sim 30^{t h}$ genes and $32^{t h} \sim 40^{t h}$ genes, respectively.

The median values of the first two loading vectors for $σ_{e}^{2} = 0.51$ and $σ_{e}^{2} = 1$ based on 500 replications, are displayed in Figures 1 and 2, respectively. The sparse PCA and nwEigenCELLs methods also demonstrates effectiveness in recovering the sparse structure of the loading matrix for high dimensional situations, that is, $x_{1} \sim x_{20}$ ( $x_{1} \sim x_{30}$ ) have non-zero loadings for the first eigenvalue, while the features $x_{21} \sim x_{40}$ ( $x_{31} \sim x_{40}$ ) have non-zero loadings for the second eigenvalue in scenarios 1 and 2, respectively. On the other hands, the fSPCA cannot effectively perform for sparse loading identification.

FIG. 1.

The median the first two loadings for $σ_{e}^{2} = 0.51$ (p = 50).

FIG. 2.

The median the first two loadings for $σ_{e}^{2} = 1$ (p = 50).

3.3. Simulation 3

The simulation studies for p = 150 (i.e., p > n) are also conducted.

Scenario 1: The n × p data matrix $X = (x_{1}, x_{2}, ..., x_{p})$ is subsequently generated using the following process: $\begin{array}{l} x_{i j} = V_{i 1} + e_{i j}, & for 1 \leq j \leq 60, \\ x_{i j} = V_{i 2} + e_{i j}, & for 60 < j \leq 120, \\ x_{i j} = V_{i 3} + e_{i j}, & for 120 < j \leq 150. \end{array}$ where $e_{i j} \sim N (0, 1)$ for $j = 1, 31, 61, 91$ , and $\begin{array}{l} e_{i j} = N (0.7 \times e_{i 1}, σ_{e}^{2}) & for 1 \leq j \leq 30, \\ e_{i j} = N (0.7 \times e_{i 31}, σ_{e}^{2}) & for 32 < j \leq 60, \\ e_{i j} = N (0.7 \times e_{i 61}, σ_{e}^{2}) & for 62 < j \leq 90, \\ e_{i j} = N (0.7 \times e_{i 91}, σ_{e}^{2}) & for 92 < j \leq 120. \end{array}$

The features $x_{1} \sim x_{60}$ form the derived first component (i.e., V₂), while the features $x_{31} \sim x_{120}$ and $x_{121} \sim x_{150}$ form the derived second and third components (i.e., V₁ and V₃), respectively. Ideally, the features $x_{1} \sim x_{60}$ have non-zero loadings for the first eigenvalue, while the features $x_{61} \sim x_{120}$ have non-zero loadings for the second eigenvalue. We suppose that the $1^{s t}, 31^{t h}, 61^{t h}, 91^{t h}$ genes are transcription factors regulating $2^{n d} \sim 30^{t h}$ genes, $32^{t h} \sim 60^{t h}$ genes, $62^{t h} \sim 90^{t h}$ genes and $92^{t h} \sim 120^{t h}$ genes, respectively.

Scenario 2: The n × p data matrix $X = (x_{1}, x_{2}, ..., x_{p})$ is subsequently generated using the following process: $\begin{array}{l} x_{i j} = V_{i 1} + e_{i j}, & for 1 \leq j \leq 90, \\ x_{i j} = V_{i 2} + e_{i j}, & for 90 < j \leq 120, \\ x_{i j} = V_{i 3} + e_{i j}, & for 120 < j \leq 150. \end{array}$ where $e_{i j} \sim N (0, 1)$ for $j = 1, 31, 61, 91$ , and $\begin{array}{l} e_{i j} = N (0.7 \times e_{i 1}, σ_{e}^{2}) & for 1 \leq j \leq 30, \\ e_{i j} = N (0.7 \times e_{i 31}, σ_{e}^{2}) & for 12 < j \leq 90, \\ e_{i j} = N (0.7 \times e_{i 61}, σ_{e}^{2}) & for 22 < j \leq 120, \\ e_{i j} = N (0.7 \times e_{i 91}, σ_{e}^{2}) & for 32 < j \leq 150. \end{array}$

The features $x_{1} \sim x_{90}$ form the derived first component (i.e., V₂), while the features $x_{91} \sim x_{120}$ and $x_{121} \sim x_{150}$ form the derived second and third components (i.e., V₁ and V₃), respectively. Ideally, the features $x_{1} \sim x_{90}$ have non-zero loadings for the first eigenvalue, while the features $x_{91} \sim x_{120}$ have non-zero loadings for the second eigenvalue. We suppose that the $1^{s t}, 31^{t h}, 61^{t h}$ and $91^{t h}$ genes are transcription factors regulating $2^{n d} \sim 30^{t h}$ genes, $32^{t h} \sim 60^{t h}$ genes, $62^{t h} \sim 90^{t h}$ genes and $92^{t h} \sim 120^{t h}$ genes, respectively.

Figures 3 and 4 show the median values of the first two loading vectors for $σ_{e}^{2} = 0.51$ and $σ_{e}^{2} = 1$ , respectively. Similar to situation p = 50, our method also shows effective results. Tt can be considered through simulation results of various scenarios that the nwEigenCELLs method has the potential to be a powerful tool for recovering the sparse structure of eigenCells profiles, specifically in the crucial marker identification. Furthermore, our method incorporate knowledge of network biology and thus can provide biologically reliable results.

FIG. 3.

The median the first two loadings for $σ_{e}^{2} = 0.51$ (p = 150).

FIG. 4.

The median the first two loadings for $σ_{e}^{2} = 1$ (p = 150).

Table 3 shows the computational time required for network-constrained eigenCell profile estimation by the proposed nwEigenCELLS. The run time indicates the total time required to estimate network-constrained eigenCell profile via regularization parameters selection.

Table 3.

Running Timings for Network-Constrained eigenCell Profile Estimation via R 4.2.1

p	Scenario 1			Scenario 2
p	10	50	150	10	50	150
Times	29.67	270.30	2713.83	15.25	107.52	2106.80

Unit: minute.

4. UNCOVERING IMMUNOGENE REGULATORY SYSTEM OF BONE MARROW BASED ON SINGLE-CELL RNA-SEQUENCING PROFILING ANALYSIS

We focus on exploring the immunogene regulatory system in AML, with a special focus on age-specific gene regulation, particularly in young and aged cell lines. The proposed method was employed to uncover immunogene regulatory networks in the human bone marrow. We use a publicly available single-cell gene expression data that has been sourced from Triana (Triana et al., 2021). The main focus of our study is the investigation of immunogene interactions within young- and aged-healthy cell lines. Our analysis is based on a dataset that encompasses expression levels for 461 genes across 49,507 cell lines, sourced from “3 healthy young and 3 healthy old bone marrow donors”. We also consider age-specific immunogene regulatory systems in AML, where the networks are estimated by “15 leukemic bone marrow donors” dataset consisting of expression levels for 458 genes in 31586 cell lines. The distribution of cell line development stages is presented in Table 4.

Table 4.
Number of Cell Lines of the Datasets “3 Healthy Young and 3 Healthy Old Bone Marrow Donors” and “15 Leukemic Bone Marrow Donors”

Data BM type $♯$ cells Development stage $♯$ cells

Healthy Young 21,026 25-year-old human stage 9,751

26-year-old human stage 11,317

34-year-old human stage 6,963

Old 28,031 59-year-old human stage 4,071

65-year-old human stage 10,639

69-year-old human stage 6,316

AML 31,586 32-year-old human stage 545

44-year-old human stage 4,149

45-year-old human stage 603

46-year-old human stage 687

51-year-old human stage 306

54-year-old human stage 527

55-year-old human stage 11,941

56-year-old human stage 1,396

61-year-old human stage 7,547

62-year-old human stage 1,341

66-year-old human stage 492

70-year-old human stage 880

71-year-old human stage 565

78-year-old human stage 607

Data	BM type	$♯$ cells	Development stage	$♯$ cells
Healthy	Young	21,026	25-year-old human stage	9,751
26-year-old human stage	11,317
34-year-old human stage	6,963
Old	28,031	59-year-old human stage	4,071
65-year-old human stage	10,639
69-year-old human stage	6,316
AML		31,586	32-year-old human stage	545
44-year-old human stage	4,149
45-year-old human stage	603
46-year-old human stage	687
51-year-old human stage	306
54-year-old human stage	527
55-year-old human stage	11,941
56-year-old human stage	1,396
61-year-old human stage	7,547
62-year-old human stage	1,341
66-year-old human stage	492
70-year-old human stage	880
71-year-old human stage	565
78-year-old human stage	607

AML, acute myeloid leukemia; BM, bone marrow.

In the dataset of “15 leukemic bone marrow donors”, we define cell lines with a development stage of $\leq 56$ as young-AML cell lines (comprising 20,154 cell lines), and cell lines with a development stage of $\geq 61$ as aged-AML cell lines (amounting to 11,432 cell lines). We consider the immune system process pathway of Gene Ontology term analysis, and then extract immnogenes annotated genes within this pathway. As a result, 165 and 164 immnogenes are identified for young/aged-healthy and young/aged-AML cell lines, respectively. In short, expression levels of 164 immnogenes in 20,154 and 11,432 cell lines are used for young- and aged-healthy-specific network-constrained eigenCell profile estimation, respectively, while the network-constrained eigenCell profile of young and aged AML cell lines are estimated by the expression levels of 165 immnogenes in 21026 and 28031 cell lines, respectively. We consider four phenotypes, young-healthy, aged-health, young-AML, and aged-AML, and the dataset of the four phenotypes are used to each phenotypes-specific gene network, i.e., young-healthy, aged-health, young-AML, and aged-AML-specific inmmunogene networks. We then uncover the immunogene regulatory systems that can highlight the unique characteristics specific to young/aged-healthy and young/aged-AML cell lines.

In order to estimate gene network, various computational strategies have been developed and used in the various research on gene network analysis, such as, gene networks from database, co-expression network based on correlation, non-linear regression model with bayesian network, gaussian graphical model with L₁-type regularization, and the graphical lasso (Friedman et al., 2008). Our strategy can apply to gene networks estimated by various network estimation strategies. In this study, we consider linear regression framework to describe gene regulatory system. The estimation of the gene regulatory network is carried out within the context of a linear regression framework, where, $x_{j} = \sum_{l \neq j} γ_{j l} x_{l} + ϵ_{j}, j = 1, ..., p,$ (7)where γ_jl is the regression coefficient describing the effect of l^th regulator gene on the $j^{t h}$ target gene and $ϵ_{j} = {(ϵ_{j 1}, ..., ϵ_{j n})}^{T}$ is a random error vector, assumed to be independently and identically distributed with a mean 0 and variance $σ^{2}$ . We employ L₁-type regularization methods, notably lasso (Tibshirani, 1996), to estimate the gene network, $γ_{j} = \underset{γ_{j}}{arg min} {\sum_{i = 1}^{n} {(x_{i j} - \sum_{l \neq j}^{p} x_{i l} γ_{j l})}^{2} + λ \sum_{l \neq j}^{p} | γ_{j l} |},$ (8)where $λ > 0$ is the regularization parameters of $γ$ and is selected by 3-fold cross validation. We then compute the adjacency matrix as follow, $W = w_{i j} = \frac{| γ_{i j} | + | γ_{j i} |}{2} .$ (9)

Because our method shows effective performance in accurately recovering the sparse structure of eigenCell profile and stable estimation results in simulation studies, we implement the proposed method on each group of cell lines, which includes young/aged-healthy as well as young/aged-AML cell lines. Following this, we identify the eigenCell profiles and crucial immnogenes for each of the four groups of cell lines. Table 5 shows the computational time required for network-constrained eigenCell profile estimation for uncovering crucial immunogene regulatory systems.

Table 5.

Running Timings Network-Constrained eigenCell Profile Estimation for Uncovering Crucial Immunogene Regulatory Systems via R 4.2.1

	AML		Healthy
	Young	Aged	Young	Aged
$♯$ Genes	164		165
$♯$ Cells	20,154	11,432	21,026	28,031
Times	439.57	253.14	671.22	451.467

Unit: minute.

AML, acute myeloid leukemia.

We extract the regulatory systems of the crucial immnogenes. The immnogenes regulatory networks in both young/aged-healthy and young/aged-AML cell lines are given in Figure 5, where the edges with an absolute value of edge weight greater than 0.1, that is, $| γ_{j l} | > 0.1$ are only considered. As depicted in Figure 5, there are marked differences in the immnogene regulatory systems between the young and aged cell lines. The subnetwork encompassing CD79A, comprising CD20, CD79B, IGKC, IGHM, and IGHD, may be seen as immunogene interactions that are specific to healthy conditions. The hubness of NKG7 is regarded as a feature specific to the aged group. The interactions associated with the S100A family, namely S100A8 and S100A9, are regarded as a network specific to the young group, and their hubness becomes disappear in aged cell lines.

FIG. 5.

Immunogene regulatory network of bone marrow for healthy and AML, where color indicate sign of edges (i.e., red:− and blue:+), arrow indicates the regulator and target genes (i.e., regulator → target). The yollow marked genes are crucial markers that shows phenotype (Yong/Aged and Health/AML) -specific molecular interplays. AML, acute myeloid leukemia.

The identified markers and their corresponding references related to AML are presented in Table 6. The markers that have been identified exhibit compelling evidence as markers for AML, and their mechanisms in AML have been elucidated in previous studies.

Table 6.

The Identified Markers from Immunogene Regulatory Networks for Young/Old-Healthy and Young/Old-AML Cell Lines

Cell lines	Genes	References
Healthy	CD79A	Kozlov et al., 2005; Luger et al., 2013; Paiva et al., 2018; Wozniak et al., 2008; Zhang et al., 2018
	CD79B
	CD20(MS4A1)	Bolouri et al., 2021
	IGHM	Xia et al., 2022; Qiu et al., 2021; Yu et al., 2020
	IGHD
	IGKC
Old	NKG7	Tang et al., 2023; Fujii et al., 2022; Zheng et al., 2021
Young	S100A8	Mondet et al., 2021; Karjalainen et al., 2019; Yang et al., 2020; Laouedj et al., 2017
Young	S100A9

4.1. Cd79 family: CD79A, CD79B

The observed level of CD79A on myeloid cells is contingent on the specific antibody clone utilized, indicating that approximately 30%–45% of AML cases show positivity for CD79a (Luger et al., 2013). Paiva et al. (2018) revealed that anti-MPO, CD79A, and TdT serve as excellent cell markers for the diagnosis and classification of acute leukemias. The expression profile of CD13, CD33, and MPO, coupled with the absence of CD19, CD79A, CD3, and TdT, is indicative of a typical phenotype for AML. Specifically, this phenotype, characterized by the expression of CD13, CD33, CD34, CD117, and MPO, and the absence of CD19, CD7, CD3, CD79A, CD14, and CD64, signifies AML originating from early progenitors and excludes monocytoid or lymphoid proliferation (Wozniak and Szlezak, 2008). The study conducted by Zhang et al. (2018) revealed that CD79B is the most prominently upregulated gene in both MDS and AML, and they uncovered novel genes, such as RAG1, CD79B, VPREB1, and C18orf1, which had not been previously associated with MDS and AML.

4.2. Cd20

Bolouri et al. (2021) have suggested that bromodomain inhibitors and immune therapies targeting CD19, CD20, CD22, and CD79A could potentially be effective in treating infant AML.

4.3. Immunoglobulin Gene Family: IGHM, IGHD, IGKC

Through next-generation sequencing, Xia et al. (2022) demonstrated that dysregulation of immunoglobulin gene family members (IGHA, IGHD, IGHE, IGHG, IGHM, IGK, and IGL) occurs in AML. The study conducted by Qiu et al. (2021) showed an elevated mRNA level of IGHDs in AML patients and emphasized the significance of IGHDs as a marker for assessing the overall survival of AML patients. Yu et al. (2017) identified immune genes such as CCL5, CXCL9, CXCL13, CD80, CD21, CD8A, and IGKC as potent prognostic indicators for AML.

4.4. NKG7

Tang et al. (2023) reported a substantial increase in the expression of cytotoxicity-related genes, GZMB, NKG7, PRF1, CST7, CTSW, GNLY, GZMH, KIR3DL2, KIR2DL3, and KIR2DS4, in decitabine-or azacitidine-treated CD44v6 CAR-T cells. CD44v6 chimeric antigen receptor T (CD44v6 CAR-T) cells demonstrate significant anti-tumor activity and have been shown to be safe for use in AML. Six tumor cytotoxicity-related molecules (GZMA, GZMB, GZMH, IFNG, FASLG, and NKG7) originating from T cells demonstrated markedly higher levels in the bone marrow in comparison to the peripheral blood in cases of AML (Fujii et al., 2022).

4.5. S100A Family: S100A8, S100A9

In AML, a co-occurrence of venetoclax resistance was noted alongside the overexpression of S100A8 and S100A9, as well as TET2 mutations (P <0.05) (Karjalainen et al., 2019). Yang et al. (2020) demonstrated that the expression levels of S100A4, S100A6, S100A8, S100A9, S100A10, S100A12, and S100B were elevated in AML compared to normal blood samples. AML is characterized by high expression levels of S100A8 and S100A9, with S100A8 expression specifically associated with an unfavorable prognosis (Laouedj et al., 2017). The presence of S100A9 promotes differentiation of AML cells, while S100A8 counteracts the differentiation induced by S100A9, thereby preserving the immature phenotype of AML (Laouedj et al., 2017). Laouedj et al. (2017) revealed that treatment with recombinant S100A9 proteins led to increased maturation of AML cells, resulting in growth arrest and prolonged survival in an AML mouse model. It has been suggested that S100A8 and S100A9 play a regulatory role in myeloid differentiation in leukemia and may hold therapeutic potential, particularly in cases of myelomonocytic and monocytic AML (Laouedj et al., 2017).

Figure 6 illustrates the expression levels of the identified markers for healthy, aged, and young cell lines. CD79A, CD79B, CD20, IGHM, IGHD, and IGKC, which are markers specific to healthy cells, display notably heightened activities in both young- and aged-healthy cell lines. Conversely, in young- and aged-AML cell lines, their activities are notably reduced. The markers specific to young cells, S100A9 and S100A9.1, demonstrate higher expression levels in young cell lines compared to aged cell lines. In brief, the markers identified for each group of cell lines display elevated activities, encompassing both molecular interplays and expression levels within the cell lines. Through our results and in line with the literature (Wozniak and Szlezak, 2008), we suggest that the activities within the subnetwork containing CD79A, CD20, CD79B, IGKC, IGHM, and IGHD may yield crucial information for understanding the mechanism behind AML progression. The decrease in activities within the CD79A subnetwork may hold pivotal roles in the progression of AML, and activators of the subnetwork could provide crucial clues for preventing AML.

FIG. 6.

Expression levels of the identified dominant markers for Healthy, Aged and Young cell lines in the four phenotypes, cell lines of Yong-health (color: cyan), Aged-Healthy (color: dark cyan), Young-AML (color: red) and Aged-AML (color: dark red). AML, acute myeloid leukemia.

5. CONCLUSIONS

This study presents a novel computational approach for eigenCell analysis. In order to effectively estimate eigenCell profiles and determine the dominant markers for a specific phenotype, we have introduced a sparse learning approach. This approach estimates the expression levels of eigenCells using only the critical markers; that is, the eigenCells estimation is based only on the genes having non-zero loadings. Additionally, we have integrated insights from network biology into the sparse eigenCell profile estimation process. This ensures that our method yields biologically reliable results and interpretations.

Monte Carlo simulations have illustrated that the proposed method exhibits exceptional performance in accurately recovering the sparse structure of eigenCell profiles (i.e., estimating the sparse loading matrix within singular value decomposition). Our method is employed to unveil essential immune markers and their regulatory systems in immunogens in young/aged-healthy and -AML cell lines. Although not the immunogenes-focused but all genes-based analysis lead to a compressive understanding of the mechanism of AML, identification of crucial immunogenes and their molecular interplay is one of goal of our study, because immune dysfunction plays a crucial role in the pathogenesis of AML. We estimate immnunogenes regulatory networks for young- and aged-healthy cell lines and young- and aged-AML cell lines. We apply the developed strategy to identify essential immune markers and their interplays in the four regions of cell lines. The subnetwork of CD79A, encompassing CD20, CD79B, IGKC, IGHM, and IGHD, is identified as immune markers exclusive to healthy cells, while the S100A family is uncovered as markers specific to young cells. The findings from our study suggest that the activities within the subnetwork may be crucial in uncovering the mechanism underlying AML progression, particularly noting the significant role played by the diminished activities in the CD79A subnetwork. Therefore, it can also be proposed that activators of the subnetwork may hold crucial clues for preventing AML.

As shown in Tables 2 and 4, the running time of the nwEigenCELLs heavily rely on the number of genes. The computational complexity is one of the limitations of our method. Furthermore, incorrectly estimated network may lead to a false result in the eigenCell profile estimation, even though incorporating knowledge of network biology (i.e., network constraint penalty) can provide biologically reliable results. That is, the gene regulatory network estimation is also crucial issue to implement our method, and the use of biological knowledge-based gene networks (e.g., protein–protein interaction) can be considered for the network-constrained eigenCell profile estimation.

The gene networks of four phenotypes, young/aged-healthy and young/aged-AML, are estimated by cell lines of each phenotype; that is, there are differences in group sample sizes. It leads to bias the output toward seeing differences and comparing molecular interplays between phenotypes. Thus, interpreting the networks needs to be handled with great care, and we consider comparing networks of phenotypes with different group sample sizes as one of the future works of our study.

Although our strategy uncovered immunogene interplays in young/aged-healthy and young/aged-AML cell lines, and the crucial marker identification was validated via literature, there are probably potentially unidentified crucial genes. Future work can involve analysis of the well-known markers and uncover molecular interplays of the markers in our results.

Footnotes

ACKNOWLEDGMENT

This research used the computational resources of the Super Computer System, Human Genome Center, Institute of Medical Science, University of Tokyo.

AUTHORS’ CONTRIBUTIONS

H.P. developed the method, performed the analysis, and drafted the article. S.M. supervised the work. All authors have read and approved the final version of the article.

AUTHOR DISCLOSURE STATEMENT

The authors declare they have no conflicting financial interests.

FUNDING INFORMATION

This work was supported by the Sungshin Women’s University Research Grant of 2024 (Grant Recipient: Heewon Park). The work of Heewon Park was supported by NRF (RS-2023–00276559). This research was also supported by AMED under Grant Number 23tk0124003h0001, 24tk0124003h0002 and JSPS KAKENHI Grant Number JP24H00009, JP22H03692.

References

Alter

, Brown

, Botstein

. Singular value decomposition for genome-wide expression data processing and modeling. Proc Natl Acad Sci U S A, 2000; 97(18):10101–10106; doi: 10.1073/pnas.97.18.10101

Alter

, Golub

. Reconstructing the pathways of a cellular system from genome-scale signals by using matrix and tensor computations. Proc Natl Acad Sci U S A, 2005; 102(49):17559–17564; doi: 10.1073/pnas.0509033102

Amezquita

, Lun

ATL

, Becht

, et al. Orchestrating single cell analysis with Bioconductor. Nat Methods, 2020; 17(2):137–145; doi: 10.1038/s41592-019-0654-x

Bolouri

, Ries

, Pardo

, et al. A B-cell developmental gene regulatory network is activated in infant AML. PLoS One, 2021; 16(11):e0259197; doi: 10.1371/journal.pone.0259197

Dohmen

, Baranovskii

, Ronen

, et al. Identifying tumor cells at the single-cell level using machine learning. Genome Biol, 2022; 23(1):123; doi: 10.1186/s13059-022-02683-1

Friedman

, Hastie

, Tibshirani

. Sparse inverse covariance estimation with the graphical lasso. Biostat, 2008; 9(3):432–441; doi: 10.1093/biostatistics/kxm045

Fujii

, Kawamata

, Shimizu

, et al. Reinvigoration of innate and adaptive immunity via therapeutic cellular vaccine for patients with AML. Molec Ther—Oncol, 2022; 27(15):315–332; doi: 10.1016/j.omto.2022.09.001

Guo

, James

, Levina

, et al. Principal Component Analysis With Sparse Fused Loadings. J Comput Graph Stat, 2010; 19(4):930–946; doi: 10.1198/jcgs.2010.08127

Jayavelu

, Wolf

, Buettner

, et al. The proteogenomic subtypes of acute myeloid leukemia. Cancer Cell, 2022; 40(3):301–317.e12; doi: 10.1016/j.ccell.2022.02.006

10.

Karjalainen

, Liu

, Kumar

, et al. Elevated expression of S100A8 and S100A9 correlates with resistance to the BCL-2 inhibitor venetoclax in AML. Leukemia, 2019; 33(10):2548–2553; doi: 10.1038/s41375-019-0504-y

11.

Khaldoyanidi

, Nagorsen

, Stein

, et al. Immune Biology of Acute Myeloid Leukemia: Implications for Immunotherapy. J Clin Oncol, 2021; 39(5):419–432; doi: 10.1200/JCO.20.00475

12.

Kozlov

, Beason

, Yu

, et al. CD79a expression in acute myeloid leukemia t(8;21) and the importance of cytogenetics in the diagnosis of leukemias with immunophenotypic ambiguity. Cancer Genet Cytogenet, 2005; 163(1):62–67; doi: 10.1016/j.cancergencyto.2005.06.002

13.

Laouedj

, Tardif

, Gil

, et al. S100A9 induces differentiation of acute myeloid leukemia cells through TLR4. Blood, 2017; 129(14):1980–1990; doi: 10.1182/blood-2016-09-738005

14.

, Peng

, Uy

, et al. Machine learning for cell type classification from single nucleus RNA sequencing data. PLoS One, 2022; 17(9):e0275070; doi: 10.1371/journal.pone.0275070

15.

, Li

. Network-constrained regularization and variable selection for analysis of genomic data. Bioinformatics, 2008; 24(9):1175–1182; doi: 10.1093/bioinformatics/btn081

16.

Luger

, Yang

, Raviv

, et al. Expression of the B-cell receptor component CD79a on immature myeloid cells contributes to their tumor promoting effects. PLoS One, 2013; 8(10):e76115; doi: 10.1371/journal.pone.0076115

17.

Mansfeldt

, Logsdon

, Debs

, et al. SPINE: SParse eIgengene NEtwork linking gene expression clusters in Dehalococcoides mccartyi to perturbations in experimental conditions. PLoS One, 2015; 10(2):e0118404; doi: 10.1371/journal.pone.0118404

18.

Mondet

, Chevalier

, Mossuz

. Pathogenic Roles of S100A8 and S100A9 Proteins in Acute Myeloid and Lymphoid Leukemia: Clinical and Therapeutic Impacts. Molecules, 2021; 26(5):1323; doi: 10.3390/molecules26051323

19.

Paiva

, Paiva

, Cavalcanti

, et al. Importance of Detection of Intracellular Myeloperoxidase, CD13, CD79a, CD22, CD3 and Terminal Deoxynucleotidyl Transferase By Flow Cytometry Diagnosis of Acute Leukemias. Blood, 2018; 132(Supplement 1):5194–5194; doi: 10.1182/blood-2018-99-118880

20.

Putnam

, Kondeti

, Kesler

MBA

, et al. Modulating the immune system as a therapeutic target for myelodysplastic syndromes and acute myeloid leukemia. Biochem Cell Biol, 2023; 101(6):481–495; doi: 10.1139/bcb-2022-0374

21.

Qiu

, Zhang

, et al. Development of a Prognostic Nomogram for Acute Myeloid Leukemia on IGHD Gene Family. Int J Gen Med, 2021; 14:4303–4316; doi: 10.2147/IJGM.S317528

22.

Singh

, Rai

, Dohare

, et al. Network-based identification of signature genes KLF6 and SPOCK1 associated with oral submucous fibrosis. Mol Clin Oncol, 2020; 12(4):299–310; doi: 10.3892/mco.2020.1991

23.

Stuart

, Butler

, Hoffman

, et al. Comprehensive Integration of Single-Cell Data. Cell, 2019; 177(7):1888–1902.e21; doi: 10.1016/j.cell.2019.05.031

24.

Sun

, Lin

, Feng

, et al. Network-regularized high-dimensional cox regression for analysis of genomic data. Stat Sin, 2014; 24(3):1433–1459; doi: 10.5705/ss.2012.317

25.

Taghiloo

, Asgarian-Omran

. Immune evasion mechanisms in acute myeloid leukemia: A focus on immune checkpoint pathways. Crit Rev Oncol Hematol, 2021; 157:103164; doi: 10.1016/j.critrevonc.2020.103164

26.

Tang

, Kong

, Wang

, et al. Demethylating therapy increases cytotoxicity of CD44v6 CAR-T cells against acute myeloid leukemia. Front Immunol, 2023; 14:1145441; doi: 10.3389/fimmu.2023.1145441

27.

Tibshirani

. Regression shrinkage and selection via the lasso. J Roy Stat Soc Ser B, 1996; 58(1):267–288; doi: 10.1111/j.2517-6161.1996.tb02080.x

28.

Trapnell

, Cacchiarelli

, Grimsby

, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol, 2014; 32(4):381–386; doi: 10.1038/nbt.2859

29.

Triana

, Vonficht

, Jopp-Saile

, et al. Single-cell proteo-genomic reference maps of the hematopoietic system enable the purification and massive profiling of precisely defined cell states. Nat Immunol, 2021; 22(12):1577–1589; doi: 10.1038/s41590-021-01059-0

30.

Wolf

, Angerer

, Theis

. SCANPY: Large-scale single-cell gene expression data analysis. Genome Biol, 2018; 19(1):15; doi: 10.1186/s13059-017-1382-0

31.

Wozniak

, Szlezak

. Standard immunophenotyping of leukemia cells in acute myeloid leukemia (AML). Central Europ J Immun, 2008; 33(1).

32.

Xia

, Wu

, Sun

, et al. Next-Generation Sequencing Revealed a Distinct Immunoglobulin Repertoire with Specific Mutation Hotspots in Acute Myeloid Leukemia. Biology (Basel), 2022; 11(2):161; doi: 10.3390/biology11020161

33.

Yang

, Jin

, Huang

, et al. Expression and clinical significance profile analysis of S100 family members in human acute myeloid leukemia. Eur Rev Med Pharmacol Sci, 2020; 24(13):7324–7334; doi: 10.26355/eurrev_202007_21896

34.

, Lim

, Wang

, et al. Enhanced construction of gene regulatory networks using hub gene information. BMC Bioinformatics, 2017; 18(1):186; doi: 10.1186/s12859-017-1576-1

35.

, Xin

, Dong

, et al. Weighted gene coexpression network analysis identifies critical genes in different subtypes of acute myeloid leukaemia. Biotech Biotech Equip, 2020; 34(1):925–936; doi: 10.1080/13102818.2020.1811767

36.

Zhang

, Zhao

, Wei

, et al. Integrated bioinformatic analysis of microarray data reveals shared gene signature between MDS and AML. Oncol Lett, 2018; 16(4):5147–5159; doi: 10.3892/ol.2018.9237

37.

Zheng

, Cai

, Ren

, et al. Identification of immune related cells and crucial genes in the peripheral blood of ankylosing spondylitis by integrated bioinformatics analysis. Peer J, 2021; 9:e12125; doi: 10.7717/peerj.12125

38.

Zjablovskaja

, Florian

. Acute Myeloid Leukemia: Aging and Epigenetics. Cancers (Basel), 2019; 12(1):103; doi: 10.3390/cancers12010103

39.

Zou

, Hastie

, Tibshirani

. Sparse principal component analysis. J Comp Grap Stat, 2006; 15(2):265–286; doi: 10.1198/106186006X113430