Feature screening for survival trait with application to TCGA high-dimensional genomic data

Data structure and methods partial review

We consider a study with $\mathrm{n}$ independent subjects. For a subject $\mathrm{i}$, suppose that there are $\mathrm{p}$ genes expression ${\left(x_{i1},\dots ,x_{ip}\right)}^{{}^{\prime }}$ related to clinical survival outcomes $T_i$. Note that the number of the genes is far greater than the sample size $\mathrm{n}$. Usually, the survival outcome is subject to censoring, so we define $C_i$ as censoring time, and use ${\delta }_{\mathrm{i}}$ as the indicator of whether the survival time of subject $\mathrm{i}$ is censored, then define $V_i=\mathrm{m}\mathrm{i}\mathrm{n}\left(T_i,C_i\right)$ as observed survival time.

We list the common and effective survival feature screening methods, and provide readers with an overview summarized in Table 1. Edelmann et al. (2020) provided an R package “MVS” that can be downloaded from https://github.com/thomashielscher/MVS, which contains the first six screening methods discussed in Table 1, and R codes for the last two screening methods discussed in Table 1 are available at https://figshare.com/articles/software/CODE/16677070, laying the foundation for meaningful comparison.

TCGA cancer data source

TCGA RNA-Seq expression data and phenotypic data including survival time and censoring status data can be downloaded from the R package ‘TCGAbiolinks’ (Colaprico et al., 2016), or ‘UCSCXenaTools’ (Wang & Liu, 2019). The TCGA ESCA, PAAD, LUAD, and BRCA genomic data with survival traits analyzed during this study are all available at Figshare: https://figshare.com/articles/dataset/DATA/16677619. The TCGA HNSCC genomic data can be downloaded from the R package “GEInter” (Wu, Qin & Ma, 2021).

Evaluation performance in the simulation study

In performance measurement, we report the percentiles of the minimum model size (MMS) statistics among 200 replications through violin plot to view the distribution of MMS data and its summary statistics, where MMS is the minimum size of a selected model, including underlying effective predictors. MMS measures the complexity of the selected model and reflects the accuracy of the screening process; a smaller MMS value indicates the higher accuracy of feature screening. We note that a violin chart can be constructed through the “vioplot” R package (Adler & Kelly, 2021). We also performed additional simulation studies to investigate the survival time prediction errors by giving the average number of c-index (Harrell, Lee & Mark, 1996) among 200 replications as a function of the number of selected features for each method. The c-index metric compares the subjects’ predicted survival time rankings with their real survival time rankings. A larger c-index indicates better prediction accuracy.

For our comparative study, we assess the similarity of different screening methods in terms of the list of detected features. To this end, for each simulation run, the Jaccard index (J(A, B) = |A∩B|/|A∪B|) of the true feature list is calculated by selection of by two methods under a specific model size of 500. The average Jaccard index of 200 simulation repetitions is used as a measure of similarity between the two methods; in this way, the similarity matrix of all the screening methods under consideration can be constructed and visualized, for example, using the “corrplot” R package (Wei & Simko, 2017).

In addition, we calculate the overlap coefficient (O(A, B) = |A∩B|/min (|A|, |B|)) set similarity analysis of features selected by each method with a ground truth set of predictors. The average number of overlap coefficient index among 200 replications as a function of the number of selected features for each method is used as a measure of similarity between the feature screening method and the ground truth set of predictors. A larger overlap coefficient index indicates highly similarity with a ground truth set of predictors.

Simulation scenarios

For the first simulated settings, we follow the simulation settings of Song et al. (2014), and generate a cohort of 500 subjects. Each subject’s survival time follows the linear transformation model

$H\left(T_i\right)=-x_i^{{}^{\prime }}{\beta }_0+\epsilon ,$where $H\left(t\right)=\mathrm{l}\mathrm{o}\mathrm{g}\left(0.5\left(e^{2t}-1\right)\right)$, the covariates $\mathrm{x}$ jointly follow a 2,000-dimensional multivariate standard normal distribution with the first-order autoregressive (AR(1)) structure that is $\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}\left(x_{.j},x_{.k}\right)={0.5}^{\left|j-k\right|}$. The distribution of the error term $\epsilon$ follows a standard extreme value distribution, which corresponds to a proportional hazards (PH) model.

The true regression coefficient vector is sparse:

For the second simulated settings, we follow the simulation settings of Edelmann et al. (2020). The simulated settings are the same as the first simulated settings except for the relationships of true parameter vector, i.e.,

$H\left(T_i\right)=-z_i^{{}^{\prime }}{\beta }_0+\epsilon ,$where ${\mathrm{z}}_{\mathrm{i}}^{{}^{\prime }}=\left(g_1\left(x_{i1}\right),\dots ,g_p\left(x_{ip}\right)\right)$, we assume $g_1\left(x\right)={\beta }_1\left|x\right|,g_2\left(x\right)={\beta }_2\left|x\right|,g_4\left(x\right)={\beta }_4{x}^2,$, $,g_5\left(x\right)={\beta }_{5}1(x>0)$, and all other j, we assume $g_j\left(x\right)={\beta }_{j}x$.

For the third simulated settings, we follow the simulation settings of Wang & Chen (2021). The simulated settings are the same as the first simulated settings except for the correlation structure of the variables and the true regression coefficient. We generate the covariates jointly following a 2,000-dimensional multivariate standard normal distribution with different network structures, including “hub”, “band”, “cluster”, and “scale-free”. These network structures can be generated by the “huge” package with huge.generator function (Zhao et al., 2012) and the true regression coefficient vector is defined as

For each simulation scenario, we perform 200 replications to investigate the numerical performances of different methods.

Survival prediction measure in real data application

To evaluate the performance of survival prediction, we report three measures of prediction accuracy: the c-index, deviance, and the number of selected features (NOSF) metrics. Larger c-index/smaller deviance and number of selected features indicates better prediction accuracy. Since the c-index metric can be used to compare a subject’s predicted survival time ranking with the true survival time ranking. And, the deviance metric is defined as

$D=-2\left(logL\left(\hat{\mathbf{\beta }}\right)-logL\left(0\right)\right),$where $logL\left(\hat{\mathbf{\beta }}\right)$ is the log partial likelihood function of Cox’s model from the testing set of the data, $\hat{\mathbf{\beta }}$ is an estimator of the penalized Cox’s regression with the MCP penalty parameter in a prediction model obtained from the training set of the data, and $logL\left(0\right)$ is the log partial likelihood function of Cox’s null model from the testing set of the data, where all predictors are assumed not related to the true survival time.

The deviance metric can therefore be considered as a comparison between the survival prediction model and null model (no predictors considered). The deviance metric is also a suitable survival prediction criterion.

Finally, we choose the number of selected features metric as a precision criterion; the reason is that the feature screening approach is used to find parsimonious precision models. Meaning, if we are modeling two sets of features with the same predictive accuracy, we want to choose the model that uses the features, as a smaller number of selected features are easier to interpret/evaluate in follow-up studies about biological function.

Simulation studies

In simulation studies, a series of simulation studies are conducted to investigate the performance of the existing feature screening methods for survival trait in identifying true associated predictors and survival prediction errors.

The simulation results are summarized in Figs. 1–6 and Figs. S1–S12. Note that Figs. 1–3 indicate the performances for the simulation study 1 with AR(1) structure; Figs. 4–6 indicate the performances for the simulation study 2 with nonlinear structure; Figs. S1–S3 indicate the performances for the simulation study 3 with band structure; Figs. S4–S6 indicate the performances for the simulation study 3 with hub structure; Figs. S7–S9 indicate the performances for the simulation study 3 with cluster structure; and Figs. S10–S12 indicate the performances for the simulation study 3 with scale-free structure. We note that the “IPCW-tau (MB)” method is also an IPCW Kendall tau method of network adjustment, but it uses the original predictors without using nonparanormal procedure to transform them.

From simulation studies 1, 2, and 3 with band structure, we see that the IPCW-tau (NPN-MB) method outperforms all alternative methods in terms of the overlap coefficient (top three panels of Figs. 1, 4, S1), and MMS measure (Figs. 2, 5, S2). Overall, the variable selection accuracy of the IPCW-tau (NPN-MB) method is substantially better than all other methods. From bottom three panels of Figures 1, 4, S1, the FAST method has a higher c-index when the number of selected variables is larger, and the IPCW-tau (NPN-MB) method outperforms most alternative methods when the number of selected variables is medium or small.

In simulation study 3 with hub structure, we see that the IPCW-tau (NPN-MB) method performs the best when the number of selected variables is larger or medium in terms of the overlap coefficient (top three panels of Fig. S4), then the CINDEX method performs the best when the number of selected variables is small. The c-index measure patterns (bottom three panels of Fig. S4) are similar to these in the previous simulation studies. On the MMS measure metric (Fig. S5), the IPCW-tau (NPN-MB) method performs the best.

In simulation study 3 with cluster and scale-free structures, we see that the CINDEX method performs the best when the censoring rate is high in terms of the overlap coefficient (top three panels of Figs. S7, S10), then the IPCW-tau (NPN-MB) method performs the best when the censoring rate is medium or small. The c-index measure patterns (bottom three panels of Figs. S7, S10) are similar to these in the previous simulation studies. In the MMS measure metric (Figs. S8, S11), the performance of IPCW-tau (NPN-MB) method is better than other methods.

According to the average of Jaccard similarity index (Figs. 3, 6, S3, S6, S9, S12), we find that screening methods belonging to the same category have higher similarity than screening methods not belonging to the same category.

A further simulation study is conducted under the scenario where the survival time follows a proportional odds (PO) model (i.e., the error term follows a standard logistic distribution) and all other settings are the same as those in the previous simulation study. We still obtain similar numeric results and the proposed IPCW-tau (NPN-MB) method has a better performance than the alternative methods at most gene structures. These correspondence figures are all available at Figshare: https://figshare.com/articles/figure/Survival_Feature_Screening_based_on_proportional_odds_model/17089013.

Real data application with TCGA ESCA data

After excluding patients with missing survival time data, our analysis is focused on the subset of the TCGA ESCA data with 368 patients and 20,501 gene expression variables. The censoring rate in the data is about 58%. As the number of disease-associated biomarkers is not expected to be large, the top 2,000 genes with the smallest p-values based on marginal Cox’s model are selected for downstream analysis. We take five random splits of the whole data into 294:74 training/test sets of the data to evaluate the performance of all methods for survival prediction in the TCGA ESCA data.

According to the procedure of Wang & Chen (2021), we apply eight screening methods, “PL”, “SIS”, “FAST”, “RCDCS”, “CRCDCS”, “CINDEX” “IPCW-tau”, “IPCW-tau (NPN-MB)”, to the TCGA ESCA data. After grid search from the top 10 to the top 300 ranked genes, the best overall prediction performance of all methods is attained by using the top 150 genes, so the top-ranked 150 predictors are selected as the candidate covariates for each method and the Cox’s regression model with the candidate covariates and the MCP penalty (Zhang, 2010) is applied to the training data to establish the final prediction model. Besides, the MCP-penalized Cox model with the top 2,000 genes selected by the univariate Cox’s test is applied to the training data to build the prediction model. We also take the published biomarker genes (DNAJB1, BNIP1, VAMP7, TBK1) related to ESCA (Du et al., 2021) as a survival prediction model to make comparisons.

The prediction accuracy performances for different methods are evaluated and the numerical results are provided in Table 2 that reports the median of the survival prediction results among five treatments. Overall, we can see that the proposed IPCW-tau (NPN-MB) method outperforms the alternative methods for survival prediction in the TCGA ESCA test data.

In addition, we apply the proposed IPCW-tau (NPN-MB) method for whole data to identify several important biomarker genes and estimate the correspondence parameters by penalized Cox’s regression model with the MCP penalty. Please see Table 3 for the list of selected associated predictors with their correspondence weights. We identify ten genes and find the two genes (GFPT1, HNRNPC) genes that are related to ESCA in the literature (Zhang et al., 2020b; and Xu, Pan & Pan, 2020).

Real data application with TCGA PAAD data

After excluding patients with missing survival time data, our analysis is focused on the subset of the TCGA PAAD data with 178 patients and 20,501 gene expression variables. The censoring rate in the data is about 48%. As the number of disease-associated biomarkers is not expected to be large, the top 2,000 genes with the smallest p-values based on marginal Cox’s model are selected for downstream analysis. We take five random splits of the whole data into 142:36 training/test sets of the data to evaluate the performance of all methods for survival prediction in the TCGA PAAD data.

According to the procedure of Wang & Chen (2021), we apply eight screening methods, “PL”, “SIS”, “FAST”, “RCDCS”, “CRCDCS”, “CINDEX” “IPCW-tau”, “IPCW-tau (NPN-MB)”, to the TCGA PAAD data. After grid search from the top 10 to the top 300 ranked genes, the best overall prediction performance of all methods is attained by using the top 20 genes, so the top-ranked 20 predictors are selected as the candidate covariates for each method, and the Cox’s regression model with the candidate covariates and the MCP penalty (Zhang, 2010) is applied to the training data to establish the final prediction model. Besides, the MCP-penalized Cox model with the top 2,000 genes selected by the univariate Cox’s test is applied to the training data to build the prediction model. We also take the published biomarker genes (CDKN2A, TP53, TTN, KRAS) related to PAAD (Baek & Lee, 2020) as a survival prediction model to make comparisons.

The prediction accuracy performances for different methods are evaluated and the numerical results are provided in Table 4 that reports the median of the survival prediction results among five folds. Overall, we can see that the proposed IPCW-tau (NPN-MB) method outperforms the alternative methods for survival prediction in the TCGA PAAD test data.

In addition, we apply the proposed IPCW-tau (NPN-MB) method for whole data to identify several important biomarker genes and estimate the correspondence parameters by penalized Cox’s regression model with the MCP penalty. Please see Table 5 for the list of selected associated predictors with their correspondence weights. We identify two genes (MET, ZMAT1) and find the MET gene that is related to PAAD in the literature (Li et al., 2021; Huang et al., 2021; Wu et al., 2019a; Vanderwerff et al., 2019; and Li et al., 2019).

Finally, the analysis results for TCGA HNSCC, TCGA LUAD, and TCGA BRCA data are provided in supplementary materials. Note that we take the published biomarker genes (GIMAP6, SELL, TIFAB, KCNA3, CCR4) related to HNSCC (Ran et al., 2021); (ALK, BRAF, EGFR, ROS1) related to LUAD (Chen et al., 2021); (TMEM190, TUBA3D, LYVE1, LILBR5, CD209) related to BRCA (Liu et al., 2019) as a survival prediction model to make comparisons. We identify nine genes and find the four genes (PITPNM3, MXD4, ABCB1, BATF) that are related to HNSCC in the literature (Aravind et al., 2021; Wu et al., 2019b; da Silva et al., 2021; Duz & Karatas, 2021; Wang et al., 2020; and Wen et al., 2015). We identify fifteen genes and find the seven genes (EPB41L5, INPP5J, KRT16, MS4A1, MYLIP, PEBP1, SFTPB) that are related to LUAD in the literature (Li et al., 2020a; Zhang et al., 2020a; Yuanhua et al., 2019; Song et al., 2020; Liu et al., 2021b; Li et al., 2020b; Zhang et al., 2021; Cao et al., 2021; and Zhang et al., 2019). We identify ten genes and find the four genes (EDA2R, PCMT1, QPRT, SKP1) that are related to BRCA in the literature (Liu, Kain & Wang, 2012; Kyritsis et al., 2021; Liu et al., 2021a; and Tian et al., 2020). The proposed IPCW-tau (NPN-MB) method consistently performs well in these cancer datasets (refer to Tables S2, S4, and S6).