DeepPIG: deep neural network architecture with pairwise connected layers and stochastic gates using knockoff frameworks for feature selection

TL;DR

DeepPIG is a new deep learning model that improves feature selection in machine learning, especially when feature signals are weak.

Contribution

DeepPIG introduces a novel deep neural network architecture with pairwise connected layers and stochastic gates for robust feature selection.

Findings

DeepPIG outperformed baseline and recent models in synthetic data without violating the FDR level.

It showed superior classification performance in real-world cancer prognosis and microbiome datasets.

The model is effective even when feature signals are weak.

Abstract

Selecting relevant feature subsets is essential for machine learning applications. Among the feature selection techniques, the knockoff filter procedure proposes a unique framework that minimizes false discovery rates (FDR). However, employing a deep neural network architecture for a knockoff filter framework requires higher detection power. Using the knockoff filter framework, we present a Deep neural network with PaIrwise connected layers integrated with stochastic Gates (DeepPIG) for the feature selection model. DeepPIG exhibited better detection power in synthetic data than the baseline and recent models such as Deep feature selection using Paired-Input Nonlinear Knockoffs (DeepPINK), Stochastic Gates (STG), and SHapley Additive exPlanations (SHAP) while not violating the preselected FDR level, especially when the signal of the features were weak. The selected features determined by DeepPIG demonstrated superior classification performance compared with the baseline model in real-world data analyses, including the prediction of certain cancer prognosis and classification tasks using microbiome and single-cell datasets. In conclusion, DeepPIG is a robust feature selection approach even when the signals of features are weak. Source code is available at https://github.com/DMCB-GIST/DeepPIG.

Figures

The architecture of DeepPIG and training strategy. (a) DeepPIG utilized the feature selection layer incorporating original and knockoff variables to select important features related to responses. (b) Detailed view of feature selection unit and their training scheme.

Simulation study results. (a–d) DeepPIG, DeepPINK, STG, and SHAP were applied to the synthetic dataset for feature selection. Empirical power and FDR of DeepPIG and DeepPINK were obtained using the knockoff+ threshold. Powers are illustrated as solid lines and FDRs as dashed lines. The preselected FDR target is 0.2, as shown in black dashed lines.

Full text

Introduction

Since the era of big data, revolutionary improvements have been made in various fields. A deep neural network (DNN) is a plausible approach for treating complex data. While DNNs offer remarkable predictive abilities in various tasks, their “black box” nature was most concerning to many experts who needed to understand data features used to make such decisions^1,2^. Furthermore, most datasets usually contain features irrelevant to the responses of interest, leading to suboptimal training or overfitting^3,4^. In this context, identifying the crucial features contributing to a specific response and reducing the feature dimensions are essential^5^.

Various feature selection and importance scoring methods have been proposed for statistics and machine learning^6^. Feature selection methods should ideally control the rate of selecting irrelevant features while maintaining high power to identify relevant features. Traditional approaches, including the Benjamini and Hochberg procedures^7–9^, use p-values that reflect feature importance. While these methods are effective for simple models, they face challenges with complex models such as DNNs. In such cases, the generation of meaningful p-values that reflect feature importance becomes ambiguous^10^. Moreover, high-dimensional data incur high costs for computing the p-values.

The model-X knockoff framework for feature selection was proposed to bypass the usage of p-values without violating the false discovery rate (FDR) above a preselected level^11^. The knockoff framework starts by generating knockoff variables that mimic the arbitrary dependence structure among the original features without looking at the responses. Knockoff variables have been used as controls in feature selection by comparing the importance of the original features and their knockoff counterparts.

Recently, feature selection approaches using modified layer architectures from vanilla neural networks have been proposed, such as stochastic gates (STG)^12^. Although they achieved a high detection power in many applications, they often did not consider FDR control explicitly. In addition, most procedures for identifying significant features depend on empirical thresholds, which make them less deterministic. Meanwhile, a DNN architecture suitable for knockoff frameworks, such as DeepPINK, was proposed^13^. DeepPINK introduces a pairwise connected filter layer for the original and knockoff variables. However, it often fails to select a single feature when the signal of important features is dim.

In this study, we designed a novel architecture of DNN and a unique training scheme called pairwise connected layers and stochastic gates (DeepPIG). We combined the core layer architecture of DeepPINK and STG to achieve higher detection power while preventing FDR violations. Furthermore, we developed distinctive training algorithms to compute feature importance using the original and corresponding knockoff variables. From the experimental results, we observed enhanced feature selection performance compared with baseline models using synthetic datasets. DeepPIG exhibited good sensitivity by finding significant features with higher power, especially when the signals of the true features were weak, while ensuring that the FDR was not violated. We also conducted the real data analysis, including cancer survival prediction tasks and the previously used microbiome and single-cell datasets in the baseline model study. The features selected by DeepPIG exhibited better classification performance and robustness than that of DeepPINK. Finally, we report the identified cancer prognostic genes, frequently identified as significant genes for classifying long-term survivors of kidney, liver, and pancreatic cancers. DeepPIG selected several prognostic genes at higher frequencies than the baseline model, highlighting its robustness. Taken together, DeepPIG provides robust feature selection by enhancing selection power and maintaining strict FDR control, making it applicable to various biological datasets.

Methods

Knockoff framework

The knockoff filter procedure was introduced as a variable selection method that controls the FDR^11^. The knockoff filter method is well known for providing accurate FDR control while bypassing p-values. The knockoff procedure has two main parts: constructing knockoff variables that imitate the original variables and defining the knockoff statistics that can be taken as feature importance scores. When generating knockoff variables, the responses of interest must not be associated. Formally, knockoff variables for a set of original random variables \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf{x}}=(x_1,...x_d )^T$$\end{document} are defined as a new set of random variables \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{{\widetilde{\hbox {x}}}}=({\widetilde{x}}_1,...{\widetilde{x}}_d )^T$$\end{document} that satisfy the following properties:

1. For any subset \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S \subset \left\{ 1,...p \right\} ,(x^T,{\widetilde{x}}^T )_{\text {swap}(S)} {\mathop {=}\limits ^{d}} (x ^T,{\widetilde{x}}^T )$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(x^T,{\widetilde{x}}^T )_{\text {swap}(S)}$$\end{document} is obtained by swapping the components \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_j$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{x}}_j$$\end{document} in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(x^T,{\widetilde{x}}^T )$$\end{document} for each \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j\in S$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathop {=}\limits ^{d}}$$\end{document} denotes equal in distribution;
2. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{x}}\perp y\ |\ x$$\end{document}

Among the methods for constructing knockoff variables, one promising approach is to use DNN-based models such as DeepLINK^14,15^. DeepLINK takes advantage of an autoencoder with flexible nonlinear factor modeling power. Because the feature vectors generated from the autoencoder are nonlinear, one can generate knockoff variables without assuming a joint distribution of x, such as Gaussian.

Next, knockoff variables were used as controls for the original variables; therefore, original variables with a significantly stronger relationship with the response than their corresponding knockoffs were considered important features. For each feature index \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1,...,d$$\end{document} , we defined \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_j$$\end{document} as the knockoff statistic to measure the importance of the j-th original feature. A large positive value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_j$$\end{document} provides evidence that the jth original feature is important, whereas small magnitudes around zero of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_j$$\end{document} are expected to be null features. Formally, knockoff statistics \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_j$$\end{document} is a function of the augmented data matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[{\textbf {X}},\varvec{{\widetilde{\hbox {X}}}}]$$\end{document} and the response vector y with a function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_j$$\end{document} which satisfies the “sign-flip” property:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} k_j([{\textbf {X}},\varvec{{\widetilde{\hbox {X}}}}]_{\text {swap}(S)},{\textbf {y}}) = \left\{ \begin{matrix} k_j([{\textbf {X}},\varvec{{\widetilde{\hbox {X}}}}],{\textbf {y}}), j\notin S\\ -k_j([{\textbf {X}},\varvec{{\widetilde{\hbox {X}}}}],{\textbf {y}}), j \in S \end{matrix}\right. \end{aligned}$$\end{document}

where S denotes any subset of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\{ 1,...,d\right\}$$\end{document} . A threshold that does not violate the target FDR is required to select significant features using the constructed knockoff statistics. The set of important features is selected as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widehat{S}}=\left\{ j:K_j\ge t \right\}$$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=T$$\end{document} or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=T_+$$\end{document} , where T is the knockoff threshold, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_+$$\end{document} is the knockoff+ threshold. For target FDR level q, the knockoff thresholds are defined as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} T= & {} \min \left\{ t>0:\frac{\left| \left\{ j:K_j\le -t \right\} \right| }{\max \left\{ \left| j:K_j\ge t \right| ,1\right\} } \le q \right\} \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} T_+= & {} \min \left\{ t>0:\frac{1+\left| \left\{ j:K_j\le -t \right\} \right| }{\max \left\{ \left| j:K_j\ge t \right| ,1\right\} } \le q \right\} \end{aligned}$$\end{document}

FDR control was achieved if knockoff statistics of the null features were symmetrically distributed.

Proposed model

Figure 1. The architecture of DeepPIG and training strategy. (a) DeepPIG utilized the feature selection layer incorporating original and knockoff variables to select important features related to responses. (b) Detailed view of feature selection unit and their training scheme.

Here, we designed a novel DNN architecture DeepPIG for constructing knockoff statistics to improve the detection power and maintain the FDR control property of knockoff filter methods. We integrated the architectures of STG and DeepPINK^12,13^. The original and knockoff variables were combined and fed into the feature selection unit, followed by the hidden layers (Fig 1a). The final output layer produced the response variable y, which is the target output of the network. A detailed view of the feature selection unit and the training strategy are illustrated in Fig. 1b. In the feature selection unit, the stochastic gates are attached to each variable and utilized in a linear combination. STG was employed to utilize the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell 0$$\end{document} -norm of features or determine the number of selected features during DNN training. Since exact \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell 0$$\end{document} regularization can be computationally expensive and intractable for high dimensions, the stochastic gate component was designed to relax the Bernoulli distribution for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell 0$$\end{document} -norm with a continuous probabilistic distribution. Stochastic gates were attached to each input feature, where the trainable parameter, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _j$$\end{document} , and the random noise, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon _j$$\end{document} , regulate the probability of the jth gate being active; these are called relaxed Bernoulli variables and are given as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} g_j=\text {max}(0,\text {min}(0,\mu _j+\epsilon _j)), \epsilon _j \sim N(0,\sigma ^2) \end{aligned}$$\end{document}

where N denotes the normal distribution with fixed variance \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document} . The relaxed Bernoulli variables were clipped, mean-shifted, and random Gaussian. Given a loss L, the stochastic gate model is trained by minimizing the empirical risk:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \min _{\theta ,\mu }\mathbb {{\widehat{E}}}_{X,Y}{\mathbb {E}}_{G}[L(f_{\theta }({\textbf{X}}\odot {\textbf{G}},Y)+\lambda \left\| {\textbf{G}} \right\| _{0} )], \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{\theta }$$\end{document} is a model parameterized by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf{G}}$$\end{document} is a random vector with D independent variables \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_j$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j\in [D]$$\end{document} , Y is a response vector and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\odot$$\end{document} denotes element-wise multiplication. STG considers features to be important if their gate probability values are high, such as 1. Although STG achieved noteworthy performance in finding important features in various experiments, it did not explicitly consider FDR control, which led to the selection of too many features and failure to control FDR in several experimental settings.

Next, we paired the original knockoff variable that successfully passed through the stochastic gate with its corresponding knockoff variable by the plug-in filter layer. The output of this layer is a linear combination of the weighted input variables:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {output}_{filter}=g_kz_kx_k + \tilde{g_k}\tilde{z_k}\tilde{x_k} \end{aligned}$$\end{document}

Through this design, the filter weights connected to the original and knockoff features compete with each other during training. The filter weights corresponding to the original features were much larger if the original features were significant for the response vectors, thereby providing evidence for the selection of important features.

To train DeepPIG, we applied two strategies: (1) a pre-training effect by masking knockoff variables into null vectors in the early stage and (2) a training-stopping criterion using the paired t-test results of the knockoff statistics K. See Algorithm 1 for the pseudocode of the training scheme. In the first phase, we only fed the original features to the model, replacing the knockoff variables with vectors in which all the elements were set to zero. In addition, gating probabilities \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{\mu }$$\end{document} were frozen with the open state (Algorithm 1 lines 2–5). When the validation loss stabilized, the model weights were restored to a point corresponding to half of the epoch of the stopping point, and the second phase began (Algorithm 1 line 7). For example, if the validation loss stabilized at epoch 10, all model weights were restored to epoch 5. After resetting, the numeric values of the filter weights z, connected to the original variables, were copied to the filter weights \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{z}}$$\end{document} which were connected to knockoff variables (Algorithm 1 line 8). Simultaneously, the gating probabilities of both variables ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{\mu }$$\end{document} ) were set to 0–1 scaled absolute values of the filter weights of the original variables z. This procedure allows the model to roughly identify probable features.

Next, the training was continued with both the original and knockoff variables to drop insignificant features until the model encountered the stopping criterion. The stopping criterion considers the knockoff statistics, which we have revised as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} K_j=Z_j-{\widetilde{Z}}_j,j=1,...,d, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_j=\mu _j (z_j w_j )^2,\ {\widetilde{Z}}_j=\widetilde{\mu }_j ({\widetilde{z}}_j w_j )^2$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w\in R^{d}$$\end{document} are products of fully connected layer weights in the reshaped dimension. Because the importances of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_j$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{Z}}_j$$\end{document} are paired and assuming that the distribution of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_j$$\end{document} is greater than that of its counterpart, we applied a paired t-test between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_j$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{Z}}_j$$\end{document} , where the alternative hypothesis is that the mean of the original importance \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_j$$\end{document} is greater than that of the knockoff importance \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{Z}}_j$$\end{document} . Outliers within two standard deviations of the mean were excluded from this test (Algorithm 1 lines 14–17). This stopping strategy was intended for FDR control because it depends on the assumption that the original and knockoff importance scores of the null features are symmetrically distributed. As the training progressed, the knockoff importance \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{Z}}_j$$\end{document} increased, and the difference between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_j$$\end{document} decreased. The training was stopped when 1) the p-value of the paired t-test was no longer significant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p > 0.05)$$\end{document} and 2) the validation loss stabilized (Algorithm 1 lines 18–20). Finally, feature selection was conducted with a knockoff filter procedure using knockoff statistics.

Algorithm 1Pseudocode for the DeepPIG training scheme

Results

Simulation studies

Synthetic data

Mirroring the simulation studies by DeepLINK, we designed the simulation experiment settings as follows: linear factor model and logistic factor model (Eqs. 8, 9).

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} x_i= & {} \Lambda f_i+\epsilon _i \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} x_{ij}= & {} \frac{c_j}{1+\exp ([1,{{\textbf {f}}}_i^T] \varvec{\lambda }_j )}+\epsilon _{ij}, j=1,...,d \end{aligned}$$\end{document}

Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_i=(f_i^1,f_i^2,f_i^3 )^T$$\end{document} is the latent factor vector, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\lambda }_j$$\end{document} are the factor loading parameters of the desirable dimensions, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_j$$\end{document} ’s are constants, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon$$\end{document} denotes random noise. All the parameters were drawn independently from the standard normal distribution N(0, 1). The response vector \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y=(y_1,...,y_n)^T$$\end{document} is assumed to depend on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_i$$\end{document} via the following linear and nonlinear link functions (Eqs. 10, 11).

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} y_\text {linear}= & {} {\textbf {x}}^T\varvec{\beta } \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} y_\text {nonlinear}= & {} \sin ({\textbf {x}}^T \varvec{\beta })\exp ({\textbf {x}}^T\varvec{\beta }) \end{aligned}$$\end{document}

where the coefficient vector \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }=(\beta _1,...,\beta _d )^T$$\end{document} . The locations of the true features were randomly selected, and the corresponding \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _j$$\end{document} was set to the amplitude of nonzero A, either positive or negative, with equal probability. The remaining features were considered to be null, and their corresponding \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _j$$\end{document} was set to zero.

Here, we set the sample size n to 1000, feature dimension d to 500, and true feature size s to 10. The values of the amplitude A varied from three to 25. For all settings, we conducted experiments 100 times with the target FDR q as 0.2.

Simulation results

The model’s feature selection performances on the synthetic datasets were determined using metrics such as power and FDR. Power is defined as the expectation of the true discovery proportion (TDP).

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {Power}:= {\mathbb {E}}[\text {TDP}] \text { with TDP}:= \frac{|S \bigcap S_0|}{|S_0|} \end{aligned}$$\end{document}

where S denotes the subset of selected features and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_0$$\end{document} denotes the subset of true features. Power can be interpreted as recall as well. In contrast, FDR is formally defined as the expectation of the false discovery proportion (FDP).

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {FDR}:= {\mathbb {E}}[\text {FDP}] \text { with FDP}:= \frac{|S \bigcap S_0^c|}{\text {max}\{|S|,1\}} \end{aligned}$$\end{document}

FDR can be interpreted as 1-precision.

DeepPIG showed better detection power than DeepLINK while controlling the FDR to less than 0.2 in various settings, as illustrated in Fig. 2. When the link function was linear, the signal of the true features increased as the amplitude increased. By contrast, the signal and amplitude are no longer monotonic to the nonlinear link function, as shown in Fig. 2b, d. Notably, DeepPIG showed higher detection power even when the signal amplitudes were weak. Statistically, we conducted the paired t-test to compare the obtained powers by DeepPIG and DeepPINK when the amplitudes were three to eight combined. The p-values were 1.61 e−08, 1.12 e−14, 5.56 e−08 and 0.695 for the settings described in Fig. 2a–d. For the last case, the p-value was 1.20e-08 when amplitudes were 20 to 25 combined. The full scales of the results are depicted in Fig. S1. Further, DeepPIG exhibited higher F1 scores than baseline models when signals were weak, as illustrated in Fig. S2. Based on these observations, it is expected that DeepPIG will be effective in identifying inconspicuous features. STG failed to control the FDR when the response vector was generated using a nonlinear link function, and the amplitudes were large.

We compared the performance of our model with other methods such as SHapley Additive exPlanations (SHAP)^16^, a representative method in the explainable AI area. We utilized the SHAP method on basic DNN models. Furthermore, we applied linear regression and recursive feature elimination (RFE) on Lasso regression for conventional feature selection approaches. We selected the top 10 features based on their shapley values or coefficients, respectively. We observed that these methods were effective in scenarios where the significant features were conspicuous but failed to control FDR when the signals were weak (Fig. 2 and Fig. S3). When the link function was linear, it was difficult to identify important features when the amplitudes were low because the amplitudes and signal were monotonic. SHAP, linear regression, and RFE showed relatively high FDRs when the amplitudes of the features were low. Conversely, when the link function was nonlinear, these methods also exhibited high FDRs for relatively high amplitudes, as the relationship between the amplitudes and signal was not monotonic.

We conducted additional experiments for hyperparameter analysis for the simulation study results. We experimented with the restore epoch that was set to 30% and 90% in addition to 50% to demonstrate the effect of weight transfer timing. The later the weight transfer occurred, DeepPIG tended to select more features. DeepPIG failed to control FDR when the restore epoch was set too late, meaning the knockoff variables were not trained enough. We observed that the optimal performance was achieved when the restore epoch was set to 50% (Fig. S4a).

Next, various ranges of regularization coefficients were assessed to determine their effects on the model parameter. A higher regularization coefficient typically resulted in a controlled FDR with decreased power, leading to selection of fewer features (Fig. S4b). To verify the robustness of our findings, we conducted the experiments again using a synthetic dataset with a different system that had different random seeds and achieved equivalent results (Fig. S5).Figure 2. Simulation study results. (a–d) DeepPIG, DeepPINK, STG, and SHAP were applied to the synthetic dataset for feature selection. Empirical power and FDR of DeepPIG and DeepPINK were obtained using the knockoff+ threshold. Powers are illustrated as solid lines and FDRs as dashed lines. The preselected FDR target is 0.2, as shown in black dashed lines.

Real data analysis

Transcriptomic markers of cancer prognosis

Predicting cancer prognosis is challenging in the field of cancer therapy. In this study, we used DeepPIG to identify cancer prognostic genes. Transcriptomic profiles and prognostic information of patients with kidney, liver, and pancreatic cancer were collected from The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) databases^17–20^. For each cancer type, we categorized the patients as long-term survivors (LTS) or short-term survivors (non-LTS) based on their survival duration and death event occurrences. Patients were labeled as LTS if their survival times were larger than a specific threshold regardless of their death event occurrence, whereas non-LTS were labeled if their survival times were less than the threshold and the event had occurred. The properties of each dataset are listed in Table 1.

The distance correlation was applied to all three datasets as a screening step before they were fed into the feature selection models^21^. The prediction performance is reportedly poor when these datasets are employed without screening steps, suggesting that the feature space must be reduced before training the models^15^. Using ICGC datasets for screening steps, prognostic genes were ranked by the distance correlation between the gene expression level and LTS status. TCGA datasets were filtered using screened genes and employed in feature selection models.

We applied DeepPIG and DeepPINK 100 times to select the features and compared the number of selected features. For DeepPIG, we set the minimum epoch for weight transfer to increase the sensitivity. For each repetition, we randomly split each dataset into 80% training and 20% testing. To ensure fairness, each repetition used identical knockoff variables and a training-test split for DeepPIG and DeepPINK. Further, repetitions that could not select any features were denoted as “empty repetitions” and excluded from the 100 repetitions when deriving prediction performance. After selection, their prediction abilities were measured using independent vanilla DNNs, and the area under the curve (AUC) and classification errors were determined as the performance metrics. For further analysis, we used all the screened features and the same number of randomly selected features as those selected by DeepPIG for the same repetition. We also computed the precision, recall, and F1 scores and have reported them in Table S7.

We observed that DeepPIG outperformed DeepPINK in terms of the performance metrics and the number of selected features in all three datasets, as summarized in Table 1. Notably, DeepPIG has fewer empty repetitions and greater robustness in feature selection. Despite searching for appropriate hyperparameters, such as the learning rate and regularization coefficients, DeepPINK showed empty repetitions with high chances. The top 10 most frequently selected genes and their selection ratios, i,e., the number of selection times when it was not an empty repetition, are reported in Table 2. DeepPIG showed higher selection ratios for top-ranked genes than DeepPINK, suggesting the robustness of DeepPIG during repetitions. The top 100 genes and their selection ratios are listed in Tables S1–S3.

Finally, we investigated the biological roles and associations of the top-ranked genes with cancers. COL11A1 is essential for bone development and collagen fiber assembly and acts as a prognostic marker in many solid cancers, including renal carcinoma^22–24^. Hepatocyte growth factor (HGF) is a pleiotropic factor that is crucial for tubular repair, regeneration after acute renal injury, renal development, and the maintenance of normal adult kidney structure and function^25,26^. Increased PLOD2 expression is often found in advanced tumors and is correlated with a poor prognosis in patients with hepatocellular carcinoma^27^. It was reported that the overexpression of Stanniocalcin 2 (STC2) was correlated with tumor growth, invasion, metastasis, and prognosis associated with many types of cancers, including liver cancer^28,29^. C15orf48 is highly expressed in pancreatic cancer and is significantly associated with the prognosis of pancreatic adenocarcinoma^30^. Furthermore, the role of RRAD in the occurrence of ferroptosis in pancreatic cancer has been previously reported^31^. Table 1. Long-term survivor classification results for cancer datasets.TissueDataset \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\text{a}}$$\end{document}

of screened featuresMean # of selected featuresEmpty repetitionsMean ± SD test AUC (mean ± SD test classification error) $\documentclass[12pt]{minimal}

			\usepackage{amsmath}
			\usepackage{wasysym} 
			\usepackage{amsfonts} 
			\usepackage{amssymb} 
			\usepackage{amsbsy}
			\usepackage{mathrsfs}
			\usepackage{upgreek}
			\setlength{\oddsidemargin}{-69pt}
			\begin{document}$$^{\text{b}}$$\end{document}$ 

MainScreeningDeepPIGDeepPINKDeepPIGDeepPINKDeepPIGDeepPINKAll screened featuresRandom features \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\text{c}}$$\end{document} KidneyTCGA-KIRC, 222 LTS (48), 82 nonLTS (24)RECA-EU, 43 LTS (60) 18 nonLTS (18)4009.001.920 / 10056 / 1000.651 ± 0.087 (0.259 ± 0.044)0.587 ± 0.087 (0.266 ± 0.030)0.596 ± 0.076 (0.308 ± 0.046)0.600 ± 0.089 (0.272 ± 0.032)LiverTCGA-LIHC, 91 LTS (36), 104 nonLTS (36)LIRI-JP, 19 LTS (48), 17 nonLTS (12)2008.272.869 / 10043 / 1000.631 ± 0.094 (0.404 ± 0.080)0.596 ± 0.109 (0.419 ± 0.085)0.661 ± 0.084 (0.374 ± 0.076)0.570 ± 0.104 (0.455 ± 0.079)PancreasTCGA-PAAD, 66 LTS (18), 66 nonLTS (18)PAAD-CA, 19 LTS (48), 43 nonLTS (12)1005.541.2616 / 10073 / 1000.573 ± 0.106 (0.455 ± 0.090)0.537 ± 0.103 (0.479 ± 0.092)0.628 ± 0.091 (0.427 ± 0.092)0.563 ± 0.111 (0.461 ± 0.091) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\rm{a}}$$\end{document} Dataset name, the number of long-term survivors (LTS) and short-term survivors (nonLTS), and survival time threshold in the month in parentheses. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\rm b}$$\end{document} Empty repetitions were excluded. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\rm c}$$\end{document} The same number of features as those of DeepPIG was randomly selected. Table 2. Top 10 frequently selected prognostic genes with selection ratio. (Ratio of times selected to non-empty repetitions).KidneyLiverPancreasRankDeepPIGDeepPINKDeepPIGDeepPINKDeepPIGDeepPINK1COL11A10.88IFI440.36PLOD20.87PLOD20.33C15orf480.8C15orf480.482HGF0.88HGF0.27STC20.67ADAM90.25RRAD0.7RRAD0.263IFI440.78COL11A10.23TMX10.6GCLM0.23PSMB80.52PSMB80.264LYPD6B0.7NTM0.23ADAM90.49STC20.23GPBAR10.45USP220.195KCNE50.67KCNE50.18PARD30.41PIK3IP10.19MAP1LC3B0.43RCOR10.196BCAT10.63NKAIN40.14IFI60.34MRPL30.16PPP1R100.39FAM19A50.197PGC0.5LYPD6B0.14MERTK0.31ZWINT0.16TGFBR30.27PPP1R100.198C16orf890.38PGC0.14GCLM0.3PARD30.14UGT2B150.19HIST1H2AC0.159HSPB70.25BCAT10.11IGFBP30.29GTPBP40.12SH3BP40.19SOX90.1510APOD0.24IRS40.11C50.26MEX3D0.12UGT2B170.18GPBAR10.11

To further investigate the significance of the selected prognostic genes, we conducted a univariate Cox proportional hazards analysis on the corresponding TCGA cohorts, including all patients whose prognostic information was available. We observed that the top-ranked genes frequently exhibited significant p-values, showing their association with survival in cancer patients, as summarized in Table 3Table 3. Survival analysis of top 10 frequently selected prognostic genes.KidneyLiverPancreasGeneHazard ratioCoxPH p-valueGeneHazard ratioCoxPH p-valueGeneHazard ratioCoxPH p-valueCOL11A11.0180.0001235PLOD21.0275.43E–06C15orf481.0190.000525HGF1.0170.00024STC21.0250.0002739RRAD1.0070.000857IFI441.0266.39E–12TMX11.0670.0056695PSMB81.0282.94E-05LYPD6B1.0170.3040778ADAM91.0883.83E–05GPBAR10.9190.032268KCNE51.0160.0399017PARD31.0546.63E–05MAP1LC3B1.0030.865728BCAT11.0673.66E–08IFI610.6973015PPP1R100.9320.000291PGC1.0190.0828003MERTK1.0110.214183TGFBR30.9940.847321C16orf890.9970.6961147GCLM1.0171.43E–05UGT2B151.0040.323858HSPB70.9940.7014953IGFBP31.0010.5672436SH3BP41.0160.289861APOD1.0160.0052317C50.9960.0007997UGT2B170.960.351064p-values under 0.05 are presented in bold.

Microbiome and single-cell datasets

We compared the performances of DeepPIG and DeepPINK using the same datasets and procedures described in the DeepLINK study^15^. The datasets used were the human microbiome dataset from a colorectal cancer (CRC) study^32,33^, the human single-cell dataset from a glioblastoma study^34^, and the murine single-cell dataset from a lipopolysaccharide (LPS)-stimulated transcriptomic effect study^35^. Human microbiome datasets were utilized to identify important microbial species related to colorectal cancer by classifying 184 individuals (91 patients with colorectal cancer and 93 healthy controls). The human single-cell dataset contained 632 cells (580 tumor cells and 52 surrounding peripheries) from patients with glioblastoma, and these single-cell gene expression data were employed to investigate the differential gene expression between both cell types. The murine single-cell dataset was collected to investigate the effect of LPS-stimulated nuclear factor- \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B (NF- \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B) on gene expression. Classification of 580 cells (202 unstimulated cells and 368 LPS-stimulated cells) based on their condition revealed significantly differentially expressed genes under both conditions.

Similar to the previous section, we applied screening steps and applied DeepPIG and DeepPINK 100 times to select the features. As shown in Table 4, DeepPIG selected more features than DeepPINK and exhibited better test errors. The frequently selected features are listed in Tables S4–S6.

The frequently selected features are similar to those in the analysis in the DeepLINK study. However, some genes were selected more frequently by DeepPIG, suggesting its sensitive detection capability. For instance, for the human microbiome dataset analysis, Parvimonas micra was selected 94 times out of 100 repetitions by DeepPIG, whereas it was selected 52 times by DeepPINK. Several recent studies have reported the biological effects of P. micra on CRC. Parvimonas micra promotes the development of CRC and can be considered as a predictor of poor outcomes in patients with CRC^36,37^. It also influences proliferation, wound healing, and inflammation in CRC cell lines^38^. For human single-cell dataset analysis, B2M and C1R were selected 47 and 42 times, respectively, using DeepPIG, compared to 6 and 3 times using DeepPINK. Some studies have reported that B2M has a significant relationship with the tumor-immune microenvironment and plays a critical role in tumor progression, patient prognosis, and immunotherapy of gliomas^39–41^. Furthermore, the expression level of C1R was associated with immune cell infiltration and prognosis of glioblastoma^42^. Table 4. Performance comparison with microbiome and single-cell datasets.Dataset# of screened featuresMean # of selected featuresEmpty repetitionsMean ± SD test classification error \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\text{a}}$$\end{document} DeepPIGDeepPINKDeepPIGDeepPINKDeepPIGDeepPINKAll screened featuresRandom features \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\text{b}}$$\end{document} Human microbiome3013.6910.660/1000/1000.275 ± 0.0750.287 ± 0.0750.298 ± 0.0720.342 ± 0.083Human single-cell20028.9317.360/10011/1000.048 ± 0.0260.053 ± 0.0320.035 ± 0.0220.068 ± 0.029Murine single-cell20050.9634.480/1000/1000.011 ± 0.0130.015 ± 0.0170.009 ± 0.0120.038 ± 0.031 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\text{a}}$$\end{document} Empty repetitions were excluded. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\text{b}}$$\end{document} The same number of features as those of DeepPIG was randomly selected.

Discussion

Various explainable AI (XAI) models such as SHAP^16^ are used for measuring feature importance and understanding which feature contribute significantly to the output of the neural network models. Although the knockoff models and XAI techniques both employ the weights of parameters to compute the feature importance, their main aspects are somewhat different. XAI techniques are mainly applied to interpret the results of trained models that are considered “black box” models that users can not catch how the model comes to the specific results. On the contrary, the knockoff framework is designed to select significant features among a large number of features while keeping the FDR not violated.

One advantage of the knockoff framework over conventional feature selection approaches or SHAP is its ability to determine the threshold using knockoff variables. In other methods, it is necessary to specify the number of features to be selected. We examined the performance of existing methods with the top 10 features, as we were aware that 10 significant features existed within the synthetic datasets. The knockoff framework is useful for real dataset analysis, especially when it is unclear how many features should be selected.

The motivation for designing DeepPIG is that the knockoff filter method often fails to select a single feature in real data analysis. We focused on constructing a sensitive feature-selection model for low-amplitude signal cases. Cancer prognosis prediction is an example of this, as genes related to survival are uncommon. DeepPIG exhibited better detection power when the amplitude of the features was weak, as demonstrated by synthetic and real data.

Some components of the proposed training strategy were determined empirically. For example, in selecting the specific epoch for weight transfer, the late transfer failed in FDR control. If DeepPIG was “overcooked” with the original variable only, the knockoff variables had an insufficient influence on the model and could not overcome the score differences between the original and knockoff variables. In this context, the “half point epoch” criterion was decided empirically based on the simulation study. Although an epoch earlier than half is acceptable, we recommend staying within half of the epoch where the validation loss is stabilized.

Additionally, the 0.05 p-value of the paired t-test for the stopping criterion was adjusted. This is because the knockoff statistics within non-outlier regions do not necessarily indicate that they are actual null features. Since a p-value of 0.05 is considered the general criterion in statistical fields, we stuck with it as a criterion for our study. We utilized paired t-tests as a “gadget” to verify how the knockoff importance scores catch up with that of the original. An alternative method for determining the time for weight transfer and testing the symmetry of the original and knockoff scores should be explored in future studies.

Conclusion

In this study, we present DeepPIG, a DNN architecture, and a training scheme for feature selection. We integrated the key structures of DeepPINK and STG to improve detection power while keeping FDR under a preselected level. Using synthetic data, we achieved a higher power, especially when the amplitudes of the features were weak. We applied DeepPIG to renal carcinoma, hepatocellular carcinoma, and pancreatic carcinoma datasets to classify patients with cancer as LTS or non-LTS. DeepPIG robustly selects several prognostic genes at high frequencies. Furthermore, we compared DeepPIG with the baseline model DeepPINK using the human microbiome, human single-cell, and murine single-cell datasets, which were employed in the baseline model study. It was observed that DeepPIG selected a greater number of features and had superior prediction capacities. In conclusion, DeepPIG is a robust feature selection approach even when the signals of features are weak.

Supplementary Information

Supplementary Information.