
==== Front
BMC Bioinformatics
BMC Bioinformatics
BMC Bioinformatics
1471-2105
BioMed Central London

5897
10.1186/s12859-024-05897-1
Research
A novel approach to the analysis of Overall Survival (OS) as response with Progression-Free Interval (PFI) as condition based on the RNA-seq expression data in The Cancer Genome Atlas (TCGA)
Lin Bo 1
Wang Kaipeng 2
Yuan Yuan 34
Wang Yueguo 4
Liu Qingyuan 5
Wang Yulan 4
Sun Jian 4
Wang Wenwen 4
Wang Huanli 6
Zhou Shusheng 4
Jin Kui 4
Zhang Mengping 1
Lai Yinglei Laiyinglei@ustc.edu.cn

17
1 grid.59053.3a 0000000121679639 School of Mathematical Sciences, University of Science and Technology of China, Hefei, 230026 Anhui China
2 https://ror.org/00xp9wg62 grid.410579.e 0000 0000 9116 9901 School of Mathematics and Statistics, Nanjing University of Science and Technology, Nanjing, 210094 Jiangsu China
3 https://ror.org/01f8qvj05 grid.252957.e 0000 0001 1484 5512 Graduate School of Bengbu Medical College, Bengbu, Anhui China
4 https://ror.org/04c4dkn09 grid.59053.3a 0000 0001 2167 9639 Department of Emergency Medicine, The First Affiliated Hospital of USTC, Division of Life Sciences and Medicine, University of Science and Technology of China, Hefei, 230001 Anhui China
5 https://ror.org/0108wjw08 grid.440647.5 0000 0004 1757 4764 School of Mathematics and Physics, Anhui Jianzhu University, Hefei, 230009 Anhui China
6 https://ror.org/04c4dkn09 grid.59053.3a 0000 0001 2167 9639 Department of Information Center, The First Affiliated Hospital of USTC, Division of Life Sciences and Medicine, University of Science and Technology of China, Hefei, 230001 Anhui China
7 https://ror.org/00y4zzh67 grid.253615.6 0000 0004 1936 9510 Department of Statistics, The George Washington University, Washington, DC USA
13 9 2024
13 9 2024
2024
25 30013 4 2024
9 8 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
Background

Overall Survival (OS) and Progression-Free Interval (PFI) as survival times have been collected in The Cancer Genome Atlas (TCGA). It is of biomedical interest to consider their dependence in pathway detection and survival prediction. We intend to develop novel methods for integrating PFI as condition based on parametric survival models for identifying pathways associated with OS and predicting OS.

Results

Based on the framework of conditional probability, we developed a family of frailty-based parametric-models for this purpose, with exponential or Weibull distribution as baseline. We also considered two classes of existing methods with PFI as a covariate. We evaluated the performance of three approaches by analyzing RNA-seq expression data from TCGA for lung squamous cell carcinoma and lung adenocarcinoma (LUNG), brain lower grade glioma and glioblastoma multiforme (GBMLGG), as well as skin cutaneous melanoma (SKCM). Our focus was on fourteen general cancer-related pathways. The 10-fold cross-validation was employed for the evaluation of predictive accuracy. For LUNG, p53 signaling and cell cycle pathways were detected by all approaches. Furthermore, three approaches with the consideration of PFI demonstrated significantly better predictive performance compared to the approaches without the consideration of PFI. For GBMLGG, ten pathways (e.g., Wnt signaling, JAK-STAT signaling, ECM-receptor interaction, etc.) were detected by all approaches. Furthermore, three approaches with the consideration of PFI demonstrated better predictive performance compared to the approaches without the consideration of PFI. For SKCM, p53 signaling pathway was detected only by our Weibull-baseline-based model. And three approaches with the consideration of PFI demonstrated significantly better predictive performance compared to the approaches without the consideration of PFI.

Conclusions

Based on our study, it is necessary to incorporate PFI into the survival analysis of OS. Furthermore, PFI is a survival-type time, and improved results can be achieved by our conditional-probability-based approach.

Keywords

TCGA
RNA-seq
Overall Survival
Progression-Free Interval
Pathway
Conditional
http://dx.doi.org/10.13039/501100018530 Major Science and Technology Projects in Anhui Province 2020b07050001 Jin Kui National Natural Science Foundation of China12126604 Zhang Mengping R&D project of Pazhou Lab (Huangpu)Grant2023K0609 Zhang Mengping issue-copyright-statement© BioMed Central Ltd., part of Springer Nature 2024
==== Body
pmcBackground

Overall Survival (OS) and Progression-Free Interval (PFI) as survival times have been widely employed to demonstrate the clinical impact of treatments, covariates, or biomarkers [1, 2]. In this study, Overall Survival (OS), defined as the duration from diagnosis to death from any cause, and Progression-Free Interval (PFI), defined as the time from diagnosis to the occurrence of a new tumor event were selected as survival times. In biomedical research, the Cox proportional hazards model [3] is predominantly utilized to analyze OS or PFI separately. The principal objective of these studies is to understand the marginal effects of certain covariates on a chosen survival time and to facilitate prediction. The detection of pathways is pivotal in revealing the intricate biological processes and interactions that significantly impact patient outcomes [4, 5]. Furthermore, it is evident that there exists some dependence between OS and PFI, and their dependence may play a critical role in unraveling the mechanisms underlying disease progression [6–8]. Therefore, it is necessary to consider their dependence in pathway detection and survival prediction. In this study, we intend to develop novel methods for integrating PFI as condition based on parametric survival models for identifying pathways associated with OS and predicting OS.

These methods can be developed within the framework of a bivariate survival analysis, focusing on how to obtain an appropriate conditional distribution for the second survival time T2, when observing a survival time T1 (whether censored or not). Although current methodologies were available for analyzing the related bivariate survival problem, they often exhibit limitations in their applicability and flexibility [9, 10]. Day et al. [9] explored a similar problem, but they assumed that times were not sequential. In contrast, Henderson et al. [10] regarded times as sequential, yet their interest was confined to instances where the observed survival time was not censored. However, Day et al. [9] and Henderson et al. [10] did not consider the situation that observed survival time was censored. In this study, we proposed an approach to accommodate this situation.

Based on the framework of conditional probability distribution, we developed a family of novel frailty-based methods [11–13] for this purpose. Our frailty-based methods can accommodate the situation when the observed survival time as condition is censored and also capture certain nonlinear relationships. In our frailty models, we consider a gamma distribution for the frailty. We developed two shared-frailty-model based methods, with exponential or Weibull as baseline hazard function. The estimators are obtained by the maximum likelihood method. We derived the related likelihood ratio test to identify covariates associated with the future survival event when the observed survival event is considered as condition. Furthermore, we achieved the prediction of future survival event.

In the literature [10], there are two classes of existing methods related to this analysis: both with T1 as a covariate; T2 was modeled by traditional parametric regression models or Cox regression models. In the Cox regression model, T1 can be used as a time-varying/time-dependent effect [14]. Both methods mirrored the approach suggested by Gail et al. [15] in their analysis of multiple tumor times.

We evaluated the performance of three approaches by analyzing RNA-seq expression data from TCGA for lung squamous cell carcinoma and lung adenocarcinoma (LUNG), brain lower grade glioma and glioblastoma multiforme (GBMLGG), as well as skin cutaneous melanoma (SKCM). Our focus was on fourteen general cancer-related pathways (Table 2). We utilized the existing single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm [16], and each pathway is summarized as a covariate. The 10-fold cross-validation approach was employed for the evaluation of predictive accuracy, which was based on the time-dependent Receiver Operating Characteristic (ROC) methodology [17–19]. In Fig. 1, we present the workflow diagram of our framework, which highlights the key steps and processes. To further illustrate our frailty-based methods, we evaluated the impact of censoring on likelihood ratio test and parameter estimation (In the Supplementary Materials). Then, we also demonstrated the advantages of frailty-based methods over traditional regression models (In the Supplementary Materials).Fig. 1 The work flow diagram of our framework

Methods

We first provide an overview of regression models used in survival analysis, beginning with the Cox proportional hazards model and traditional regression models, and then extending to shared gamma frailty models.

Background: Cox proportional hazards model

The Cox proportional hazards model [3] is a widely used method in survival analysis that relates the time of an event to covariates without making specific assumptions about the underlying hazard function. The model can be described by the following hazard function:λ(t|x)=λ0(t)exp(βx),

where λ(t|x) is the hazard function, λ0(t) is the baseline hazard, β represents the coefficients, and x represents the covariates.

Background: time-varying Cox model

The Cox model can be extended to handle time-varying/time-dependent effects [20], where covariates change over time or coefficients change over time. The hazard function for a Cox model with time-varying covariates or time-varying coefficients can be written as:λ(t|x(t))=λ0(t)exp(βx(t)),λ(t|x)=λ0(t)exp(β(t)x),

where x(t) denotes the covariates that change over time, β(t) denotes the coefficients that change over time. This extension allows for more flexibility in modeling scenarios where the risk factors change as time progresses.

Background: baseline hazard functions and regression approach

Two baseline hazard functions are considered in the analysis: the exponential distribution and the Weibull distribution (Table 1). In practice, traditional regression model is often used to analysis ordered survival events. The probability density function of the traditional regression model is [10]f(t2|x,t1)=e-t2δλeβx+γlog(t1)eβx+γlog(t1)λδt2δ-1.

where δ, λ, β and γ represent unknown parameter, x represents covariate, t1 and t2 represent two survival times.Table 1 Two baseline hazard functions for shared gamma frailty models and regression approach

Distribution	λ0(t)	
Exponential	λ,λ>0	
Weibull	λρtρ-1,λ>0,ρ>0	
Regression approach	Density function f(t2)=e-t2δλeβx+γlog(t1)eβx+γlog(t1)λδt2δ-1	

Conditional models

Background: shared gamma frailty models

Shared frailty models are widely used models in survival analysis [21, 22]. In these models, we assume that the two time components are T1 and T2. In shared frailty models, we assume that T1 and T2 are conditionally independent given a shared unobservable random effect, namely the frailty Z. We consider that the frailty Z follows a gamma distribution (Z∼Γ(θ,θ),θ>0). Given the gamma frailty Z with the observed covariate x, the conditional hazard function and the conditional cumulative hazard function of T1 and T2 are as followλi(t|Z,x)=Zλi0(t)eβixi=1,2;Λi(t|Z,x)=ZΛi0(t)eβixi=1,2.

where λ10(t) and λ20(t) are baseline hazard functions, Λ10(t) and Λ20(t) are baseline cumulative hazard functions.

The conditional survival functions of T1 and T2 are as followS1(t|Z,x)=e-ZΛ10(t)eβ1x;S2(t|Z,x)=e-ZΛ20(t)eβ2x.

Because T1 and T2 are conditionally independent given the gamma frailty Z, the bivariate survival function of T1and T2 given the gamma frailty Z isS(t1,t2|Z,x)=e-Z(Λ10(t1)eβ1x+Λ20(t2)eβ2x).

Then, the bivariate density function of T1 and T2 is (The derivation details are described in the Supplementary Materials.)f(t1,t2|x)=λ10(t1)λ20(t2)e(β1+β2)x(θ+1)θθ+1(θ+Λ10(t1)eβ1x+Λ20(t2)eβ2x)θ+2.

Conditional distribution

The conditional density function of T2 given T1 is (The derivation details are described in the Supplementary Materials.)P(T2=t2|T1=t1,X=x)=λ20(t2)eβ2x(θ+1)(θ+Λ10(t1)eβ1x)θ+1(θ+Λ10(t1)eβ1x+Λ20(t2)eβ2x)θ+2.

In practice, it is often to observe right censored T1 and T2. The corresponding conditional probabilities are as follow (The derivation details are described in the Supplementary Materials)P(T2=t2|T1≥t1,X=x)=λ20(t2)θeβ2x(θ+Λ10(t1)eβ1x)θ(θ+Λ10(t1)eβ1x+Λ20(t2)eβ2x)θ+1;P(T2≥t2|T1=t1,X=x)=(θ+Λ10(t1)eβ1x)θ+1(θ+Λ10(t1)eβ1x+Λ20(t2)eβ2x)θ+1;P(T2≥t2|T1≥t1,X=x)=(θ+Λ10(t1)eβ1x)θ(θ+Λ10(t1)eβ1x+Λ20(t2)eβ2x)θ.

It can be seen from the above formulas that our models actually do not assume that T1≤T2. However, in practice, it is often to observe T1≤T2. Based on the above relationships, we find that the conditional probability can maintain an increasingly monotone relationship between T1 and T2 (regardless of whether T1 is censored or not), which is commonly observed in practice. For example, in the TCGA datasets, if the tumor recurrence could be delayed (PFI), then the survival probability (OS) could be higher.

Model interpretation

The hazard functions of conditional model are as followλ(t2|T1=t1,X=x)=∂[-log(∫t2+∞P(T2=t|T1=t1)dt)]∂t2=(θ+1)λ20(t2)eβ2xθ+Λ10(t1)eβ1x+Λ20(t2)eβ2x=(θ+1)λ20(t2)θe-β2x+Λ10(t1)e(β1-β2)x+Λ20(t2);λ(t2|T1≥t1,X=x)=∂[-log(∫t2+∞P(T2=t|T1≥t1)dt)]∂t2=θλ20(t2)eβ2xθ+Λ10(t1)eβ1x+Λ20(t2)eβ2x=θλ20(t2)θe-β2x+Λ10(t1)e(β1-β2)x+Λ20(t2).

It can be seen from the above formulas that the predictive score of evaluating individual hazard is θe-β2x+Λ10(t1)e(β1-β2)x. We can obtain the relationship of the individual hazard and covariate x as follow β2≥0,β1≤β2.

The conditional hazard function increases with respect to x;

β2≥0,β1>β2.

The conditional hazard function with respect to x increases at [0,1β1logβ2θΛ10(t1)(β1-β2)] and decreases at (1β1logβ2θΛ10(t1)(β1-β2),+∞];

β2<0,β1≥β2.

The conditional hazard function decreases with respect to x;

β2<0,β1<β2.

The conditional hazard function with respect to x decreases at [0,1β1logβ2θΛ10(t1)(β1-β2)] and increases at (1β1logβ2θΛ10(t1)(β1-β2),+∞].

Clearly, censored t1 does not affect the relationship above.

Theoretical results

We have proven the five lemmas required for our models in the Supplementary Materials. In Lemma3 and Lemma4, we generalize two baseline hazard functions to any positive function and generalize the frailty distribution to the power-variance-function (PVF) family [11].

Likelihood-based estimation

We assume that the sample data were (t1i,t2i,δ1i,δ2i,xii=1⋯n), where δ1i=0 represents censored t1i, δ1i=1 represents uncensored t1i, δ2i=0 represents censored t2i and δ2i=1 represents uncensored t2i.

In the exponential-baseline-based model, the baseline hazard functions areλ10(t1)=λ1,λ20(t2)=λ2,Λ10(t1)=λ1t1,Λ20(t2)=λ2t2,λ1>0,λ2>0.

The likelihood function L(θ,λ1,λ2,β1,β2) is∏i=1n(λ2eβ2xi(θ+1)(θ+λ1t1ieβ1xi)θ+1(θ+λ1t1ieβ1xi+λ2t2ieβ2xi)θ+2)δ1iδ2i(λ2eβ2xiθ(θ+λ1t1ieβ1xi)θ(θ+λ1t1ieβ1xi+λ2ti2eβ2xi)θ+1)(1-δ1i)δ2i((θ+λ1t1ieβ1xi)θ+1(θ+λ1t1ieβ1xi+λ2t2ieβ2xi)θ+1)δ1i(1-δ2i)((θ+λ1t1ieβ1xi)θ(θ+λ1t1ieβ1xi+λ2t2ieβ2xi)θ)(1-δ1i)(1-δ2i).

In the Weibull-baseline-based model, the baseline hazard functions areλ10(t1)=λ1ρ1t1ρ1-1,λ20(t2)=λ2ρ2t2ρ2-1,Λ10(t1)=λ1t1ρ1,Λ20(t2)=λ2t2ρ2,λ1>0,λ2>0,ρ1>0,ρ2>0.

The likelihood function L(θ,λ1,λ2,β1,β2,ρ1,ρ2) is∏i=1n(λ2ρ2t2iρ2-1eβ2xi(θ+1)(θ+λ1t1iρ1eβ1xi)θ+1(θ+λ1t1iρ1eβ1xi+λ2t2iρ2eβ2xi)θ+2)δ1iδ2i(λ2ρ2t2iρ2-1eβ2xiθ(θ+λ1t1iρ1eβ1xi)θ(θ+λ1t1iρ1eβ1xi+λ2ti2ρ2eβ2xi)θ+1)(1-δ1i)δ2i((θ+λ1t1iρ1eβ1xi)θ+1(θ+λ1t1iρ1eβ1xi+λ2t2iρ2eβ2xi)θ+1)δ1i(1-δ2i)((θ+λ1t1iρ1eβ1xi)θ(θ+λ1t1iρ1eβ1xi+λ2t2iρ2eβ2xi)θ)(1-δ1i)(1-δ2i).

The maximum likelihood estimators asymptotically follow normal distributions with mean being the true value of the parameters and covariance matrix being the inverse of the n-fold Fisher information matrix by the Lemma1 (In the Supplementary Materials). We calculated the maximum likelihood estimates by the optim function in R or the NLMIXED Procedure in SAS.

Variance and confidence interval

In Section Likelihood-based Estimation, we have stated that the estimated variances of the maximum likelihood estimates of β1 and β2 can be calculated by the inverse of the n-fold Fisher information matrix and we can calculate the estimated variances of the maximum likelihood estimates of β1 and β2 by the optim function in R or the NLMIXED Procedure in SAS. The maximum likelihood estimates of β1 and β2 asymptotically follow normal distributions. Therefore, the 95% confidence interval of β1 and β2 are as follow: assume the estimated standard deviations of β1 and β2 are respectively s1^,s2^, the estimated n-fold Fisher information matrix is H^ and the estimation of β1 and β2 are respectively β1^,β2^, the 95% confidence interval of β1 and β2 are respectively[β1^-1.96∗s1^,β1^+1.96∗s1^]=[β1^-1.96∗H^-1(1,1),β1^+1.96∗H^-1(1,1)];[β2^-1.96∗s2^,β2^+1.96∗s2^]=[β2^-1.96∗H^-1(2,2),β2^+1.96∗H^-1(2,2)].

Hypothesis testing

When we consider T1, the relationship between T2 and the covariate x is not only related to β2, it may also be related to β1. According to Lemma2 (In the Supplementary Materials), β1=β2=0 is equivalent to that T2 is independent of the covariate x given T1 in three shared frailty models.

Therefore, we can consider the following hypothesis testH0:β1=β2=0.vsH1:β1≠0orβ2≠0.

For the above hypothesis test, we used the likelihood ratio test.

The likelihood ratio statistics are as followΛ=sup(θ,λ1,λ2,β1,β2)∈R5L(θ,λ1,λ2,β1,β2)sup(θ,λ1,λ2,β1,β2)∈(R,R,R,0,0)L(θ,λ1,λ2,β1,β2);Λ=sup(θ,λ1,λ2,ρ1,ρ2,β1,β2)∈R7L(θ,λ1,λ2,ρ1,ρ2,β1,β2)sup(θ,λ1,λ2,ρ1,ρ2,β1,β2)∈(R,R,R,R,R,0,0)L(θ,λ1,λ2,ρ1,ρ2,β1,β2);

According to Lemma5 (In the Supplementary Materials), the above likelihood ratio test statistics have asymptotical distributions, i.e. the double log-likelihood ratio test statistics asymptotically follow chi-square distribution with two degrees of freedom (χ2(2)). Therefore, we test the covariate x by the likelihood ratio test. The above likelihood ratio test statistics Λ can be calculated using the optim function in R or the NLMIXED Procedure in SAS.

Single-sample gene set enrichment analysis (ssGSEA)

We utilized the R package GSVA to perform ssGSEA [16] at the individual sample level, using pathways as gene sets. The pathway information was obtained from the Molecular Signatures Database (MSigDB) [23, 24].

Kaplan–Meier (K–M) curve

We utilized the R package survival to create a Kaplan–Meier (K–M) [25] curves. The K–M curves estimates and visualizes the survival probability over time.

Time-dependent receiver operating characteristic (ROC)

We use the time-dependent ROC methodology (a modified version of the conventional ROC methodology) to evaluate the predictive accuracy of our models [17–19]. The areas under the time-dependent ROC curves (AUC(t)) are often used to compare the predictive accuracy of several prediction models. Since Heagerty [17, 18] proposed the time-dependent ROC methodology and some definitions of cases and controls, there are many methods to estimate the time-dependent ROC curve. We choose a nonparametric estimator of the time-dependent ROC curve using the inverse probability of censoring weighting (IPCW) approach [19]. We plot AUC(t) of IPCW estimators by the R package timeROC.

Results

TCGA data

The RNA sequencing expression data from The Cancer Genome Atlas (TCGA) database were considered. Liu et al. [26] recently curated and standardized genomic and survival information for the TCGA dataset into the TCGA Pan-Cancer Clinical Data Resource (TCGA-CDR). Therefore, we downloaded genomic and survival data for the TCGA datasets from UCSC Xena Browser (https://xenabrowser.net/datapages/). In this study, Overall Survival (OS), defined as the duration from diagnosis to death from any cause, and Progression-Free Interval (PFI), defined as the time from diagnosis to the occurrence of a new tumor event were selected as endpoints because these were deemed to be relatively accurate endpoints by Liu et al. But since PFI contains information about OS, we did not consider the case where PFI equals OS. Then the analysis was conducted with the PFI as T1 and OS as T2. Time is measured in days. Among 186 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, we primarily focus on the analytical findings for a subset specifically associated with cancer. There exists a total of sixteen pathways in KEGG but fourteen of them are catalogued within the MSigDB. In this study, ssGSEA was performed on these fourteen pathways (Table 2) to calculate cancer-related predictive scores.

In this study, we focused on cancer datasets with relatively large sample sizes. After excluding missing values and removing cases where PFI equals OS, there are three datasets that can be used for patient survival analysis. Therefore, our subsequent analyses are based on lung squamous cell carcinoma and lung adenocarcinoma (LUNG; n=306), brain lower grade glioma and glioblastoma multiforme (GBMLGG; n=276), as well as skin cutaneous melanoma (SKCM; n=268).

LUNG data

Data description Fig. 2 The scatter plot of PFI and OS. The solid points represent the data with both PFI and OS not censored. The empty points represent the data with PFI or OS censored. The x-axis represents PFI while y-axis represents OS. a LUNG Dataset. b GBMLGG Dataset

Table 2 Fourteen cancer-related pathways and analysis results of six models for OS in LUNG dataset

KEGG Pathway	Cox1	Coxtime	Regression1	Regression2	Exponential	Weibull	
Wnt signaling	0.177(−)	0.359(−)	0.119(−)	0.545(−)	0.241(−)	0.262(−)	
TGF-beta signaling	0.614(−)	0.462(−)	0.489(−)	0.639(−)	0.558(−)	0.610(−)	
JAK-STAT signaling	0.405(−)	0.462(−)	0.432(−)	0.279(−)	0.843(−)	0.976(−)	
MAPK signaling	0.121(−)	0.107(−)	0.056(−)	0.063(−)	0.366(−)	0.467(−)	
ECM-receptor interaction	0.901(−)	0.743(−)	0.688(−)	0.890(−)	0.858(−)	0.920(−)	
Cytokine-cytokine receptor interaction	0.416(−)	0.342(−)	0.426(−)	0.171(−)	0.434(+)	0.983(−)	
Focal adhesion	0.919(−)	0.907(−)	0.682(−)	1.000(−)	0.456(−)	0.992(−)	
Adherens junction	0.604(−)	0.583(−)	0.513(−)	0.435(−)	0.934(−)	0.968(−)	
mTOR signaling	0.845(−)	0.717(−)	0.956(−)	0.733(−)	0.452(−)	0.920(+)	
PPAR signaling	0.298(−)	0.303(−)	0.275(−)	0.363(−)	0.688(−)	0.669(−)	
VEGF signaling	0.288(−)	0.674(−)	0.217(−)	0.500(−)	0.577(−)	0.928(−)	
Apoptosis	1.000(+)	0.394(−)	0.876(+)	0.244(−)	0.952(+)	0.999(−)	
p53 signaling	0.014(+)	0.001(+)	0.023(+)	0.001(+)	0.003(+)	0.000(+)	
Cell cycle	0.009(+)	0.001(+)	0.008(+)	0.000(+)	0.000(+)	0.000(+)	
In the table, “Cox1” represents the Cox regression model (consider corresponding pathway as only covariate) hypothesis test significance, “Coxtime” represents the Cox regression model (consider corresponding pathway as a covariate, but PFI as a time-varying effect, using time-varying coefficient) hypothesis test significance, “Regression1” represents the traditional regression model (consider corresponding pathway as only covariate) hypothesis test significance,“Regression2” represents the traditional regression model (consider corresponding pathway and PFI as covariates) hypothesis test significance, “exponential” represents the exponential-baseline-based model hypothesis test significance, “Weibull” represents the Weibull-baseline-based model hypothesis test significance. The level of 0.05 is considered for claiming significance. The symbols in parentheses represent the positive or negative dependence of the pathway on the OS. All model coefficients are provided in Table S1–S3 (In the Supplementary Materials)

From Fig. 2, we can observe that PFI and OS are mostly concentrated in the interval [0, 4000]. It is evident from Fig. 2 that PFI and OS are highly positively correlated. Among these patients, the censoring proportion for PFI is 0.0% while the censoring proportion for OS is 35.6%. An intuitive overview of gene set activity across the lung cancer samples is provided by the heatmap of the GSVA output in the LUNG dataset, shown in Figure S1 (In the Supplementary Materials).

Pathway detection

Six different models were performed to identify pathways significantly associated with OS: “Cox1” represents the Cox regression model (consider corresponding pathway as only covariate), “Coxtime” represents the Cox regression model (consider corresponding pathway as a covariate, but PFI as a time-varying/time-dependent effect, using time-varying coefficient), “Regression1” represents the traditional regression model (consider corresponding pathway as only covariate), “Regression2” represents the traditional regression model (consider corresponding pathway and PFI as covariates), the exponential-baseline-based model, the Weibull-baseline-based model. Table 2 presents the analysis results of six models:Notably, these six models all detected p53 signaling pathway and cell cycle pathway significantly positively correlated with OS.

It is well-known that p53 signaling pathway and cell cycle pathway are closely linked to the development of lung squamous cell carcinoma (LUSC) and lung adenocarcinoma (LUAD) [27–30].Fig. 3 AUC(t2) curve within 5 years or survival curve for t2 in LUNG dataset. a AUC(t2) curve. b Survival curve for t2. In plot, x-axis represents time t2 when y-axis represents AUC(t2) or survival probability. “Weibull” represents the Weibull-baseline-based model, “Coxtime” represents the Cox regression model (consider pathway scores that are not highly correlated as covariates, but PFI as a time-varying effect, using time-varying coefficient), “Regression2” represents the traditional regression model (consider pathway scores that are not highly correlated and PFI as covariates). The solid lines represent AUC(t2) curves or survival curve for t2. The dashed lines represent the pointwise confidence intervals of AUC(t2) curves or survival curve for t2

Predictive performance

In our study, we employed a 10-fold cross-validation method to compare the predictive performance of the models. Initially, the entire dataset was divided into ten subsets of equal size. This method involved ten independent cycles of training and validation. In each cycle, a unique subset was selected as the validation set, while the remaining nine subsets were combined to form the training set. This approach ensured that each subset served as a validation set once and as part of the training set at other times. Following each training and validation cycle, the predictive score of models on the current validation set was assessed. Upon completion of all ten cycles, the predictive scores from the ten validation sets (encompassing the entire dataset) were aggregated to derive the predictive score of models for the full dataset.

We compared the predictive performance of six models. For each fold (10-fold cross-validation), based on the training data, the association p-values were first computed; the pairwise pearson correlations were then calculated. A network adjacency matrix was constructed accordingly based on the absolute correlation value (>0.6); for each connected component in the network, only the pathway with the smallest associated p-value was selected. (Please see Table S13–S15 in the supplementary materials for the related details, including the Cox1 model, the Coxtime model, the Regression1 model, the Regression2 model, the exponential-baseline-based model, the Weibull-baseline-based model.) The predictive performance was evaluated using the time-dependent ROC methodology. The predictive scores of the Cox1 model, the Coxtime model, the Regression1 model and the Regression2 model are all based on linear predictors. The predictive scores of the exponential-baseline-based model and the Weibull-baseline-based model are based on θe-β2Tx+Λ10(t1)e(β1-β2)Tx (see details in Section Model Interpretation).

In Fig. 3 and Figure S5 (In the Supplementary Materials), we can observe that our Weibull-baseline-based model, the Coxtime model and the Regression2 model demonstrate significantly better predictive performance in comparison to the Cox1 model and the Regression1 model in [0, 3000]. And we also can observe that the exponential-baseline-based model, the Cox1 model and the Regression1 model demonstrate poor predictive capabilities, with the AUC values being approximately 0.5.

GBMLGG data

Data description

From Fig. 2, we can similarly observe that PFI and OS are mostly concentrated in the interval [0, 4000]. It is evident from Fig. 2 that PFI and OS are highly positively correlated. Among these patients, the censoring proportion for PFI is 0.0% while the censoring proportion for OS is 34.1%. An intuitive overview of gene set activity across the glioma cancer samples is provided by the heatmap of the GSVA output in the GBMLGG dataset, shown in Figure S2 (In the Supplementary Materials).Table 3 Fourteen cancer-related pathways and analysis results of six models for OS in GBMLGG dataset

KEGG Pathway	Cox1	Coxtime	Regression1	Regression2	Exponential	Weibull	
Wnt signaling	0.000(−)	0.000(−)	0.000(−)	0.000(−)	0.000(−)	0.000(−)	
TGF-beta signaling	0.019(−)	0.515(−)	0.006(−)	0.068(−)	0.069(+)	0.115(−)	
JAK-STAT signaling	0.000(+)	0.000(+)	0.000(+)	0.000(+)	0.001(+)	0.002(+)	
MAPK signaling	0.001(−)	0.081(−)	0.000(−)	0.102(−)	0.114(−)	0.043(−)	
ECM-receptor interaction	0.000(+)	0.000(+)	0.000(+)	0.000(+)	0.000(+)	0.000(+)	
Cytokine-cytokine receptor interaction	0.000(+)	0.000(+)	0.000(+)	0.000(+)	0.000(+)	0.009(+)	
Focal adhesion	0.000(+)	0.000(+)	0.000(+)	0.000(+)	0.000(+)	0.000(+)	
Adherens junction	0.320(+)	0.619(+)	0.427(+)	0.690(+)	0.840(+)	0.911(−)	
mTOR signaling	0.000(−)	0.000(−)	0.000(−)	0.000(−)	0.000(−)	0.000(−)	
PPAR signaling	0.061(+)	0.084(+)	0.020(+)	0.118(+)	0.378(−)	0.592(−)	
VEGF signaling	0.000(+)	0.001(+)	0.000(+)	0.001(+)	0.081(−)	0.007(+)	
Apoptosis	0.000(+)	0.000(+)	0.000(+)	0.000(+)	0.001(+)	0.009(+)	
p53 signaling	0.000(+)	0.000(+)	0.000(+)	0.000(+)	0.000(+)	0.000(+)	
Cell cycle	0.000(+)	0.000(+)	0.000(+)	0.000(+)	0.000(+)	0.000(+)	
In the table, “Cox1” represents the Cox regression model (consider corresponding pathway as only covariate) hypothesis test significance, “Coxtime” represents the Cox regression model (consider corresponding pathway as a covariate, but PFI as a time-varying effect, using time-varying coefficient) hypothesis test significance, “Regression1” represents the traditional regression model (consider corresponding pathway as only covariate) hypothesis test significance,“Regression2” represents the traditional regression model (consider corresponding pathway and PFI as covariates) hypothesis test significance, “exponential” represents the exponential-baseline-based model hypothesis test significance, “Weibull” represents the Weibull-baseline-based model hypothesis test significance. The level of 0.05 is considered for claiming significance. The symbols in parentheses represent the positive or negative dependence of the pathway on the OS. All model coefficients are provided in Table S4–S6 (In the Supplementary Materials)

Pathway detection

Similar six different models were performed to identify pathways significantly associated with OS. Table 3 presents the analysis results of six models:The TGF-beta signaling pathway exhibited a significantly negative correlation with OS in the Cox1 and the Regression1 models. In contrast, the Coxtime model and the Regression2 model detected that TGF-beta signaling pathway is not significant correlated with OS. Similarly, the Weibull model detected that TGF-beta signaling pathway is not significant correlated with OS.

The Cox1 model and the Regression1 model both detected MAPK signaling pathway significantly negatively correlated with OS. Similarly, the Weibull model detected that MAPK signaling pathway is significantly negatively correlated with OS. But the Coxtime model and the Regression2 models detected that MAPK signaling pathway is not significant correlated with OS.

The Regression1 model detected PPAR signaling pathway significantly positively correlated with OS. The Weibull model detected that PPAR signaling pathway is not significant correlated with OS.

Notably, the Cox1 model, the Coxtime model, the Regression1 model, the Regression2 model and the Weibull model all detected that eight pathways (JAK-STAT signaling, ECM-receptor interaction, cytokine-cytokine receptor interaction, focal adhesion, VEGF signaling, apoptosis, p53 signaling, cell cycle) are significantly positively correlated with OS, two pathways (Wnt signaling, mTOR signaling) are significantly negatively correlated with OS.

For brain lower grade glioma (LGG) and glioblastoma multiforme (GBM), many pathways have been extensively studied. For instance, it has been studied that Wnt signaling pathway, JAK-STAT signaling pathway, ECM-receptor interaction pathway, cytokine-cytokine receptor interaction pathway and focal adhesion pathway are related to the growth and survival of tumors [31–35].Fig. 4 AUC(t2) curve within 5 years or survival curve for t2 in GBMLGG dataset. a AUC(t2) curve. b Survival curve for t2. In plot, x-axis represents time t2 when y-axis represents AUC(t2) or survival probability. “Weibull” represents the Weibull-baseline-based model, “Coxtime” represents the Cox regression model (consider pathway scores that are not highly correlated as covariates, but PFI as a time-varying effect, using time-varying coefficient), “Regression2” represents the traditional regression model (consider pathway scores that are not highly correlated and PFI as covariates). The solid lines represent AUC(t2) curves or survival curve for t2. The dashed lines represent the pointwise confidence intervals of AUC(t2) curves or survival curve for t2

Predictive performance

Similarly, we employed a 10-fold cross-validation method to compare the predictive performance of the models. We compared the predictive performance of six models. For each fold (10-fold cross-validation), based on the training data, the association p-values were first computed; the pairwise pearson correlations were then calculated. A network adjacency matrix was constructed accordingly based on the absolute correlation value (>0.6); for each connected component in the network, only the pathway with the smallest associated p-value was selected. (Please see Table S16–S18 in the supplementary materials for the related details, including the Cox1 model, the Coxtime model, the Regression1 model, the Regression2 model, the exponential-baseline-based model, the Weibull-baseline-based model.) The predictive performance was evaluated using the time-dependent ROC methodology. The predictive scores of the Cox1 model, the Coxtime model, the Regression1 model and the Regression2 model are all based on linear predictors. The predictive scores of the exponential-baseline-based model and the Weibull-baseline-based model are based on θe-β2Tx+Λ10(t1)e(β1-β2)Tx (see details in Section Model Interpretation).

In Fig. 4 and Figure S6 (In the Supplementary Materials), we can observe that our Weibull-baseline-based model, the Coxtime model and the Regression2 model demonstrate better predictive performance in comparison to the Cox1 model and the Regression1 model. And we can observe that the exponential-baseline-based model demonstrate poor predictive capabilities, with the AUC values being approximately 0.5.

SKCM data

Data description Fig. 5 The scatter plot of PFI and OS in the SKCM dataset. The solid points represent the data with both PFI and OS not censored. The empty points represent the data with PFI or OS censored. The x-axis represents PFI while y-axis represents OS

From Fig. 5, we can similarly observe that PFI and OS are mostly concentrated in the interval [0, 7000]. It is evident from Fig. 5 that PFI and OS are highly positively correlated. Among these patients, the censoring proportion for PFI is 0.0% while the censoring proportion for OS is 39.6%. An intuitive overview of gene set activity across the skin cancer samples is provided by the heatmap of the GSVA output in the SKCM dataset, shown in Figure S3 (In the Supplementary Materials).Table 4 Fourteen cancer-related pathways and analysis results of six models for OS in SKCM dataset

KEGG Pathway	Cox1	Coxtime	Regression1	Regression2	Exponential	Weibull	
Wnt signaling	0.265(+)	0.009(+)	0.254(+)	0.013(+)	0.001(+)	0.000(+)	
TGF-beta signaling	0.811(+)	0.400(+)	0.767(+)	0.403(+)	0.585(+)	0.613(+)	
JAK-STAT signaling	0.003(−)	0.001(−)	0.002(−)	0.000(−)	0.013(−)	0.000(−)	
MAPK signaling	0.122(−)	0.227(−)	0.140(−)	0.183(−)	0.823(+)	0.820(−)	
ECM-receptor interaction	0.631(+)	0.972(−)	0.620(+)	0.876(−)	0.235(+)	0.627(+)	
Cytokine-cytokine receptor interaction	0.002(−)	0.000(−)	0.002(−)	0.000(−)	0.000(−)	0.000(−)	
Focal adhesion	0.544(+)	0.554(+)	0.527(+)	0.674(+)	0.316(+)	0.428(+)	
Adherens junction	0.025(+)	0.000(+)	0.019(+)	0.001(+)	0.060(+)	0.000(+)	
mTOR signaling	0.118(+)	0.005(+)	0.126(+)	0.008(+)	0.048(+)	0.006(+)	
PPAR signaling	0.621(+)	0.876(+)	0.650(+)	0.892(+)	0.552(+)	0.981(+)	
VEGF signaling	0.692(+)	0.546(+)	0.679(+)	0.750(+)	0.138(+)	0.513(+)	
Apoptosis	0.010(−)	0.003(−)	0.006(−)	0.000(−)	0.045(−)	0.005(−)	
p53 signaling	0.400(+)	0.851(−)	0.389(+)	0.802(−)	0.958(+)	0.000(+)	
Cell cycle	0.049(+)	0.016(+)	0.045(+)	0.014(+)	0.362(+)	0.032(+)	
In the table, “Cox1” represents the Cox regression model (consider corresponding pathway as only covariate) hypothesis test significance, “Coxtime” represents the Cox regression model (consider corresponding pathway as a covariate, but PFI as a time-varying effect, using time-varying coefficient) hypothesis test significance, “Regression1” represents the traditional regression model (consider corresponding pathway as only covariate) hypothesis test significance,“Regression2” represents the traditional regression model (consider corresponding pathway and PFI as covariates) hypothesis test significance, “exponential” represents the exponential-baseline-based model hypothesis test significance, “Weibull” represents the Weibull-baseline-based model hypothesis test significance. The level of 0.05 is considered for claiming significance. The symbols in parentheses represent the positive or negative dependence of the pathway on the OS. All model coefficients are provided in Tables S7–S9 (In the Supplementary Materials)

Pathway detection

Similar six different models were performed to identify pathways significantly associated with OS. Table 4 presents the analysis results of six models:The Wnt signaling pathway and mTOR signaling pathway exhibited a significantly positive correlation with OS in the Coxtime model and the Regression2 model. Similarly, the Weibull model detected that Wnt signaling pathway and mTOR signaling pathway are significantly positive correlated with OS. In contrast, the Cox1 and the Regression1 models detected that Wnt signaling pathway and mTOR signaling pathway are not significantly correlated with OS.

The Weibull model detected p53 signaling pathway significantly positively correlated with OS, but this correlation was not detected by other models.

Notably, the Cox1 model, the Coxtime model, the Regression1 model, the Regression2 model and the Weibull model all detected that two pathways (adherens junction, cell cycle) are significantly positively correlated with OS, three pathways (JAK-STAT signaling, cytokine-cytokine receptor interaction, apoptosis) are significantly negatively correlated with OS.

For skin cutaneous melanoma (SKCM), many pathways have been extensively studied. For instance, it has been shown that JAK-STAT signaling pathway is related to the disease development, as well as cell cycle pathway, p53 signaling pathway and apoptosis pathway [36]; it also well-known that mTOR signaling pathway plays an important role in tumor progression [37].Fig. 6 AUC(t2) curve within 5 years or survival curve for t2 in SKCM dataset. a AUC(t2) curve. b Survival curve for t2. In plot, x-axis represents time t2 when y-axis represents AUC(t2) or survival probability. “Weibull” represents the Weibull-baseline-based model, “Coxtime” represents the Cox regression model (consider pathway scores that are not highly correlated as covariates, but PFI as a time-varying effect, using time-varying coefficient), “Regression2” represents the traditional regression model (consider pathway scores that are not highly correlated and PFI as covariates). The solid lines represent AUC(t2) curves or survival curve for t2. The dashed lines represent the pointwise confidence intervals of AUC(t2) curves or survival curve for t2

Predictive performance

Similarly, we employed a 10-fold cross-validation method to compare the predictive performance of the models. We compared the predictive performance of six models. For each fold (10-fold cross-validation), based on the training data, the association p-values were first computed; the pairwise pearson correlations were then calculated. A network adjacency matrix was constructed accordingly based on the absolute correlation value (>0.6); for each connected component in the network, only the pathway with the smallest associated p-value was selected. (Please see Table S19–S21 in the supplementary materials for the related details, including the Cox1 model, the Coxtime model, the Regression1 model, the Regression2 model, the exponential-baseline-based model, the Weibull-baseline-based model.) The predictive performance was evaluated using the time-dependent ROC methodology. The predictive scores of the Cox1 model, the Coxtime model, the Regression1 model and the Regression2 model are all based on linear predictors. The predictive scores of the exponential-baseline-based model and the Weibull-baseline-based model are based on θe-β2Tx+Λ10(t1)e(β1-β2)Tx (see details in Section Model Interpretation).

In Fig. 6 and Figure S7 (In the Supplementary Materials), we can observe that our Weibull-baseline-based model, the Coxtime model and the Regression2 model demonstrate significantly better predictive performance in comparison to the Cox1 model and the Regression1 model in [0, 4500]. And we also can observe that the exponential-baseline-based model, the Cox1 model and the Regression1 model demonstrate poor predictive capabilities, with the AUC values being approximately 0.5.

Discussion and conclusions

We presented some comparisons of several models using TCGA datasets for LUNG, GBMLGG and SKCM. For LUNG, p53 signaling and cell cycle pathways were detected significantly positively correlated with OS by all approaches. Our Weibull-baseline-based model, the Coxtime model and the Regression2 model demonstrate significantly better predictive performance in comparison to the approaches without the consideration of PFI in [0, 3000]. For GBMLGG, eight pathways (JAK-STAT signaling, ECM-receptor interaction, cytokine-cytokine receptor interaction, focal adhesion, VEGF signaling, apoptosis, p53 signaling, cell cycle) were detected significantly positively correlated with OS by all approaches, two pathways (Wnt signaling, mTOR signaling) were detected significantly negatively correlated with OS by all approaches. Our Weibull-baseline-based model, the Coxtime model and the Regression2 model demonstrate better predictive performance in comparison to the approaches without the consideration of PFI. For SKCM, p53 signaling pathway was detected significantly positively correlated with OS by our Weibull-baseline-based model but not by the other approaches. Our Weibull-baseline-based model, the Coxtime model and the Regression2 model demonstrate significantly better predictive performance in comparison to the approaches without the consideration of PFI in [0, 4500]. In the Supplementary Materials, we present additional analysis based on a cohort of breast invasive carcinoma (BRCA) samples (n=111). For BRCA, MAPK signaling pathway was detected significantly negatively correlated with OS by our Weibull-baseline-based model but not by the other approaches. Our Weibull-baseline-based model demonstrate significantly better predictive performance in comparison to the approaches without the consideration of PFI in [0, 3000]. Based on our study, it is necessary to incorporate PFI into a survival analysis so that additional pathways associated with OS may be discovered, and the predictive performance for OS can be improved. PFI is a survival-type time. Based on the conditional probability framework, further improved results can be obtained by our proposed approaches. We also performed Cox regression with time-varying/time-dependent covariates for prediction performance evaluation [38], and our conclusion still remains, please see Figs. S46–S49 (In the Supplementary Materials) for more details. In practice, it may be generally difficult to obtain PFI in advance. In some instances, the PFI may be observed in advance [39, 40]. If PFI could be observed in advance, it would enhance the prediction of OS. This would encourage us to collect information about PFI as much as possible. Due to the correlations between PFI and OS, the prediction performance on OS can be significantly enhanced. In the analysis of bivariate correlated survival data, it is common to construct frailty models. This approach helps account for the correlation and unobserved heterogeneity between bivariate survival times.

Due to the complexity of frailty models, parameter estimation through maximum likelihood optimization may be sensitive to the choice of initial values. Therefore, we recommend exploring multiple initial values during model application to ensure appropriate parameter estimation. The baseline hazard function for Weibull-baseline based data exhibits a strictly increasing trend. The exponential-baseline-based model and the piecewise-exponential-baseline-based model (In the Supplementary Materials) are less likely to capture this trend. So in the Weibull-baseline based data, we observed that the likelihood ratio test distribution does not asymptotically follow χ2(2) for the exponential-baseline-based model and piecewise-exponential-baseline-based model (In the Supplementary Materials). And the corresponding parameter estimation and hypothesis testing are abnormal. Furthermore, when identifying associated covariates, excessive zero values for covariate may result in the failure of the likelihood ratio test. In this case, appropriate transformations to “x” (covariate) can solve this issue. In our study, we use bulk RNA-seq data. When analyzing bulk RNA-seq data, mixed cell types can lead to confounding results. Deconvolution can separate the gene expression levels of each cell type, thereby avoiding the confusion caused by cell mixture. We obtained a ranked list of genes through the deconvolution method (NMF), and then used Gene Set Enrichment Analysis (gseapy module in Python) to assess the significance of fourteen cancer-related pathways. The specific results are shown in Table S37–S40.

For the piecewise-exponential-baseline-based model (In the Supplementary Materials), there are several methods to divide the interval, such as quantiles or empirical hazard plots [41]. In our simulation studies (In the Supplementary Materials), the intervals of the piecewise-exponential-baseline-based model match the optimal intervals (consistent with the data generation process). In this case, its performance is also not clearly higher to the Weibull-baseline-based model. We primarily focused on the frailty following a gamma distribution. Other distributions like log-normal or inverse Gaussian can also be considered [42], but the computing will become complicated. Our models primarily utilize shared frailty models. However, it is worth noting that they can also be expanded to the multivariable frailty models, as demonstrated in prior works [43–46]. Additionally, another generalization of shared frailty models can be achieved through the incorporation of copula models, as explored in previous studies [47–49]. When confronted with high censoring rates, our models can be optimized by using the multiple imputation method, as indicated in previous studies [50, 51]. Moreover, our models exhibit the versatility to accommodate not only random right censoring but also random left censoring or random interval censoring, as explored in the literature [52–55]. To enhance our model adaptability, the baseline hazard function can be approximated with spline functions [8]. The spline functions provide a more flexible approach to capture the shape of the hazard function.

Additional file

Supplementary file 1

Abbreviations

OS Overall Survival

PFI Progression-Free Interval

TCGA The Cancer Genome Atlas

LUNG Lung squamous cell carcinoma and lung adenocarcinoma

GBMLGG Brain lower grade glioma and glioblastoma multiforme

SKCM Skin cutaneous melanoma

ssGSEA Single-sample Gene Set Enrichment Analysis

ROC Receiver Operating Characteristic

PVF Power-variance-function

MSigDB Molecular Signatures Database

K-M Kaplan-Meier

IPCW Inverse probability of censoring weighting

TCGA-CDR TCGA Pan-Cancer Clinical Data Resource

KEGG Kyoto Encyclopedia of Genes and Genomes

GSEA Gene Set Enrichment Analysis

LUSC Lung squamous cell carcinoma

LUAD Lung adenocarcinoma

LGG Brain lower grade glioma

GBM Glioblastoma multiforme

Supplementary Information

The online version contains supplementary material available at 10.1186/s12859-024-05897-1.

Acknowledgements

During the preparation of this work the authors used ChatGPT in order to improve language and readability. After using this tool/service, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.

Author Contributions

BL, KJ, MZ, and YL designed the research; BL and YL developed the methods; BL, KW, YW (Yueguo Wang), QL, YW (Yulan Wang), JS, WW, YY, HW, SZ, KJ and YL contributed to the acquisition, analysis, and interpretation of the data; BL and YL drafted the manuscript; KJ, MZ and YL revised the manuscript; All authors read and approved the final manuscript.

Funding

YL was partially supported by a start-up fund from the University of Science and Technology of China. KJ was partially supported by Anhui Major Science and Technology Project under grant number 2020b07050001. MZ was partially supported by the National Natural Science Foundation of China under grant number 12126604 and R &D project of Pazhou Lab (Huangpu) under Grant2023K0609.

Availability of data and materials

For each type of tumor, we downloaded TCGA genomic data and survival data based on the following dataset ID from UCSC Xena Browser (https://xenabrowser.net/datapages/).

∙ Lung Cancer: TCGA.LUNG.sampleMap/HiSeqV2_PANCAN and survival/LUNG_survival.txt; ∙ Lower grade glioma and glioblastoma: TCGA.GBMLGG.sampleMap/HiSeqV2_PANCAN and survival/GBMLGG_survival.txt; ∙ Melanoma: TCGA.SKCM.sampleMap/HiSeqV2_PANCAN and survival/SKCM_survival.txt; ∙ Breast Cancer: TCGA.BRCA.sampleMap/HiSeqV2_PANCAN and survival/BRCA_survival.txt. The code used in this study is available at https://github.com/1LinBo/Survival/tree/main.

Declarations

Ethics Approval and Consent to Participate

This study adhered to the established guidelines of TCGA and UCSC, consequently, the requirement for ethical approval and patient informed consent was waived (https://xenabrowser. net/datapages/).

Consent for Publication

Not applicable.

Conflict of interest

The authors declare that they have no competing financial interests.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Hudis CA Barlow WE Costantino JP Proposal for standardized definitions for efficacy end points in adjuvant breast cancer trials: the STEEP system J Clin Oncol 2007 25 15 2127 2132 10.1200/JCO.2006.10.3523 17513820
Hudis CA, Barlow WE, Costantino JP. Proposal for standardized definitions for efficacy end points in adjuvant breast cancer trials: the STEEP system. J Clin Oncol. 2007;25(15):2127–32.17513820 10.1200/JCO.2006.10.3523
2. Punt CJA Buyse M Köhne C Endpoints in adjuvant treatment trials: a systematic review of the literature in colon cancer and proposed definitions for future trials J Natl Cancer Inst 2007 99 13 998 1003 10.1093/jnci/djm024 17596575
Punt CJA, Buyse M, Köhne C. Endpoints in adjuvant treatment trials: a systematic review of the literature in colon cancer and proposed definitions for future trials. J Natl Cancer Inst. 2007;99(13):998–1003.17596575 10.1093/jnci/djm024
3. David Cox R David R Regression models and life tables J Roy Stat Soc 1972 34 2 187 220 10.1111/j.2517-6161.1972.tb00899.x
David Cox R, David R. Regression models and life tables. J Roy Stat Soc. 1972;34(2):187–220.10.1111/j.2517-6161.1972.tb00899.x
4. Xie HY Wang WJ Sun FY Proteomics analysis to reveal biological pathways and predictive proteins in the survival of high-grade serous ovarian cancer Sci Rep 2017 7 1 9896 10.1038/s41598-017-10559-9 28852147
Xie HY, Wang WJ, Sun FY, et al. Proteomics analysis to reveal biological pathways and predictive proteins in the survival of high-grade serous ovarian cancer. Sci Rep. 2017;7(1):9896.28852147 10.1038/s41598-017-10559-9
5. Lyudovyk O Shen YF Tatonetti NP Pathway analysis of genomic pathology tests for prognostic cancer subtyping J Biomed Inform 2019 98 103286 10.1016/j.jbi.2019.103286 31499184
Lyudovyk O, Shen YF, Tatonetti NP, et al. Pathway analysis of genomic pathology tests for prognostic cancer subtyping. J Biomed Inform. 2019;98: 103286.31499184 10.1016/j.jbi.2019.103286
6. Michiels S Le MA Buyse M Surrogate endpoints for overall survival in locally advanced head and neck cancer: meta-analyses of individual patient data Lancet Oncol 2009 10 4 341 350 10.1016/S1470-2045(09)70023-3 19246242
Michiels S, Le MA, Buyse M, et al. Surrogate endpoints for overall survival in locally advanced head and neck cancer: meta-analyses of individual patient data. Lancet Oncol. 2009;10(4):341–50.19246242 10.1016/S1470-2045(09)70023-3
7. Rondeau V Pignon JP Michiels S A joint model for the dependence between clustered times to tumour progression and deaths: a meta-analysis of chemotherapy in head and neck cancer Stat Methods Med Res 2015 24 6 711 729 10.1177/0962280211425578 22025414
Rondeau V, Pignon JP, Michiels S. A joint model for the dependence between clustered times to tumour progression and deaths: a meta-analysis of chemotherapy in head and neck cancer. Stat Methods Med Res. 2015;24(6):711–29.22025414 10.1177/0962280211425578
8. Emura T Nakatochi M Murotani K A joint frailty-copula model between tumour progression and death for meta-analysis Stat Methods Med Res 2017 26 6 2649 2666 10.1177/0962280215604510 26384516
Emura T, Nakatochi M, Murotani K, et al. A joint frailty-copula model between tumour progression and death for meta-analysis. Stat Methods Med Res. 2017;26(6):2649–66.26384516 10.1177/0962280215604510
9. Day R Bryant J Lefkopoulou M Adaptation of bivariate frailty models for prediction, with application to biological markers as prognostic indicators Biometrika 1997 84 1 45 56 10.1093/biomet/84.1.45
Day R, Bryant J, Lefkopoulou M. Adaptation of bivariate frailty models for prediction, with application to biological markers as prognostic indicators. Biometrika. 1997;84(1):45–56.10.1093/biomet/84.1.45
10. Henderson R Prince H Choice of conditional models in bivariate survival Stat Med 2000 19 4 563 574 10.1002/(SICI)1097-0258(20000229)19:4<563::AID-SIM356>3.0.CO;2-K 10694736
Henderson R, Prince H. Choice of conditional models in bivariate survival. Stat Med. 2000;19(4):563–74.10694736 10.1002/(SICI)1097-0258(20000229)19:4<563::AID-SIM356>3.0.CO;2-K
11. Duchateau L Janssen P The frailty model 2008 Berlin Springer
Duchateau L, Janssen P. The frailty model. Berlin: Springer; 2008.
12. Manda SOM A nonparametric frailty model for clustered survival data Commun Stat-Theory Methods 2011 40 5 863 875 10.1080/03610920903480882
Manda SOM. A nonparametric frailty model for clustered survival data. Commun Stat-Theory Methods. 2011;40(5):863–75.10.1080/03610920903480882
13. Gasperoni F Ieva F Paganoni AM Non-parametric frailty Cox models for hierarchical time-to-event data Biostatistics 2020 21 3 531 544 10.1093/biostatistics/kxy071 30590499
Gasperoni F, Ieva F, Paganoni AM, et al. Non-parametric frailty Cox models for hierarchical time-to-event data. Biostatistics. 2020;21(3):531–44.30590499 10.1093/biostatistics/kxy071
14. Rogoz B, de l’Aulnoit A H, Duhamel A, et al. Thirty-year trends of survival and time-varying effects of prognostic factors in patients with metastatic breast cancer—a single institution experience. Clin Breast Cancer. 2018;18(3):246–53.
15. Gail MH Santner TJ Brown CC An analysis of comparative carcinogenesis experiments based on multiple times to tumor Biometrics 1980 1 255 266 10.2307/2529977
Gail MH, Santner TJ, Brown CC. An analysis of comparative carcinogenesis experiments based on multiple times to tumor. Biometrics. 1980;1:255–66.10.2307/2529977
16. Hänzelmann S Castelo R Guinney J GSVA: gene set variation analysis for microarray and RNA-seq data BMC Bioinf 2013 1 1
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf. 2013;1:1.
17. Heagerty PJ Lumley T Pepe MS Time-dependent ROC curves for censored survival data and a diagnostic marker Biometrics 2000 56 2 337 344 10.1111/j.0006-341X.2000.00337.x 10877287
Heagerty PJ, Lumley T, Pepe MS. Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics. 2000;56(2):337–44.10877287 10.1111/j.0006-341X.2000.00337.x
18. Heagerty PJ Zheng Y Survival model predictive accuracy and ROC curves Biometrics 2005 61 1 92 105 10.1111/j.0006-341X.2005.030814.x 15737082
Heagerty PJ, Zheng Y. Survival model predictive accuracy and ROC curves. Biometrics. 2005;61(1):92–105.15737082 10.1111/j.0006-341X.2005.030814.x
19. Blanche P Dartigues JF Jacqmin-Gadda H Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks Stat Med 2013 32 30 5381 5397 10.1002/sim.5958 24027076
Blanche P, Dartigues JF, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat Med. 2013;32(30):5381–97.24027076 10.1002/sim.5958
20. Fisher LD Lin DY Time-dependent covariates in the Cox proportional-hazards regression model Annu Rev Public Health 1999 20 1 145 157 10.1146/annurev.publhealth.20.1.145 10352854
Fisher LD, Lin DY. Time-dependent covariates in the Cox proportional-hazards regression model. Annu Rev Public Health. 1999;20(1):145–57.10352854 10.1146/annurev.publhealth.20.1.145
21. Mwikali Muli A, Houwing-Duistermaat J, Gusnanto A. Use of shared gamma frailty model in analysis of survival data in twins. Use of Shared Gamma Frailty Model in Analysis of Survival Data in Twins, 2021; pp. 45–58.
22. Esayas LM Akessa GM Kifle DD Application of parametric shared frailty models to analyze time-to-death of gastric cancer patients J Gastrointest Cancer 2023 54 1 104 116 10.1007/s12029-021-00775-y 35064523
Esayas LM, Akessa GM, Kifle DD. Application of parametric shared frailty models to analyze time-to-death of gastric cancer patients. J Gastrointest Cancer. 2023;54(1):104–16.35064523 10.1007/s12029-021-00775-y
23. Mootha VK Lindgren CM Eriksson K PGC-1α-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes Nat Genet 2003 34 3 267 273 10.1038/ng1180 12808457
Mootha VK, Lindgren CM, Eriksson K. PGC-1-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet. 2003;34(3):267–73.12808457 10.1038/ng1180
24. Subramanian A Tamayo P Mootha V Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles Proc Natl Acad Sci 2005 102 43 15545 15550 10.1073/pnas.0506580102 16199517
Subramanian A, Tamayo P, Mootha V. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci. 2005;102(43):15545–50.16199517 10.1073/pnas.0506580102
25. Kaplan Edward L Meier Paul Nonparametric estimation from incomplete observations J Am Stat Assoc 1958 53 282 457 481 10.1080/01621459.1958.10501452
Kaplan Edward L, Meier Paul. Nonparametric estimation from incomplete observations. J Am Stat Assoc. 1958;53(282):457–81.10.1080/01621459.1958.10501452
26. Liu JF Lichtenberg T Hoadley KA An integrated TCGA pan-cancer clinical data resource to drive high-quality survival outcome analytics Cell 2018 173 2 400 416 10.1016/j.cell.2018.02.052 29625055
Liu JF, Lichtenberg T, Hoadley KA. An integrated TCGA pan-cancer clinical data resource to drive high-quality survival outcome analytics. Cell. 2018;173(2):400–16.29625055 10.1016/j.cell.2018.02.052
27. Kaiser AM Gatto A Hanson KJ p53 governs an AT1 differentiation programme in lung cancer suppression Nature 2023 619 7971 851 859 10.1038/s41586-023-06253-8 37468633
Kaiser AM, Gatto A, Hanson KJ. p53 governs an AT1 differentiation programme in lung cancer suppression. Nature. 2023;619(7971):851–9.37468633 10.1038/s41586-023-06253-8
28. Bretz AC Gittler MP Charles JP ΔNp63 activates the Fanconi anemia DNA repair pathway and limits the efficacy of cisplatin treatment in squamous cell carcinoma Nucl Acids Res 2016 44 7 3204 3218 10.1093/nar/gkw036 26819410
Bretz AC, Gittler MP, Charles JP. Np63 activates the Fanconi anemia DNA repair pathway and limits the efficacy of cisplatin treatment in squamous cell carcinoma. Nucl Acids Res. 2016;44(7):3204–18.26819410 10.1093/nar/gkw036
29. Liu W Du Q Mei T Comprehensive analysis the prognostic and immune characteristics of mitochondrial transport-related gene SFXN1 in lung adenocarcinoma BMC Cancer 2024 24 1 94 10.1186/s12885-023-11646-z 38233752
Liu W, Du Q, Mei T. Comprehensive analysis the prognostic and immune characteristics of mitochondrial transport-related gene SFXN1 in lung adenocarcinoma. BMC Cancer. 2024;24(1):94.38233752 10.1186/s12885-023-11646-z
30. Ye G Luo H Zhang T Knockdown of RNF183 suppressed proliferation of lung adenocarcinoma cells via inactivating the STAT3 signaling pathway Cell Cycle 2022 21 9 948 960 10.1080/15384101.2022.2035617 35104174
Ye G, Luo H, Zhang T. Knockdown of RNF183 suppressed proliferation of lung adenocarcinoma cells via inactivating the STAT3 signaling pathway. Cell Cycle. 2022;21(9):948–60.35104174 10.1080/15384101.2022.2035617
31. Precilla SD Biswas I Kuduvalli SS Crosstalk between PI3K/AKT/mTOR and WNT/β-Catenin signaling in GBM-Could combination therapy checkmate the collusion? Cell Signal 2022 95 110350 10.1016/j.cellsig.2022.110350 35525406
Precilla SD, Biswas I, Kuduvalli SS. Crosstalk between PI3K/AKT/mTOR and WNT/-Catenin signaling in GBM-Could combination therapy checkmate the collusion? Cell Signal. 2022;95: 110350.35525406 10.1016/j.cellsig.2022.110350
32. Heynckes S Daka K Franco P Crosslink between Temozolomide and PD-L1 immune-checkpoint inhibition in glioblastoma multiforme BMC Cancer 2019 19 1 7 10.1186/s12885-019-5308-y 30606139
Heynckes S, Daka K, Franco P. Crosslink between Temozolomide and PD-L1 immune-checkpoint inhibition in glioblastoma multiforme. BMC Cancer. 2019;19:1–7.30606139 10.1186/s12885-019-5308-y
33. Jiang Y He J Guo Y Identification of genes related to low-grade glioma progression and prognosis based on integrated transcriptome analysis J Cell Biochem 2020 121 5–6 3099 3111 10.1002/jcb.29577 31886582
Jiang Y, He J, Guo Y. Identification of genes related to low-grade glioma progression and prognosis based on integrated transcriptome analysis. J Cell Biochem. 2020;121(5–6):3099–111.31886582 10.1002/jcb.29577
34. Zan X Li L Construction of lncRNA-mediated ceRNA network to reveal clinically relevant lncRNA biomarkers in glioblastomas Oncol Lett 2019 17 5 4369 4374 30944630
Zan X, Li L. Construction of lncRNA-mediated ceRNA network to reveal clinically relevant lncRNA biomarkers in glioblastomas. Oncol Lett. 2019;17(5):4369–74.30944630
35. Alowaidi F Hashimi SM Alqurashi N Cripto-1 overexpression in U87 glioblastoma cells activates MAPK, focal adhesion and ErbB pathways Oncol Lett 2019 18 3 3399 3406 31452820
Alowaidi F, Hashimi SM, Alqurashi N. Cripto-1 overexpression in U87 glioblastoma cells activates MAPK, focal adhesion and ErbB pathways. Oncol Lett. 2019;18(3):3399–406.31452820
36. Zhang W Zhao H Chen J Mining database for the expression and gene regulation network of JAK2 in skin cutaneous melanoma Life Sci 2020 253 117600 10.1016/j.lfs.2020.117600 32234492
Zhang W, Zhao H, Chen J. Mining database for the expression and gene regulation network of JAK2 in skin cutaneous melanoma. Life Sci. 2020;253: 117600.32234492 10.1016/j.lfs.2020.117600
37. Jiang Y Hu X Wang Z RPTOR mutation: a novel predictor of efficacious immunotherapy in melanoma Invest New Drugs 2024 42 1 60 69 10.1007/s10637-023-01413-z 38071684
Jiang Y, Hu X, Wang Z. RPTOR mutation: a novel predictor of efficacious immunotherapy in melanoma. Invest New Drugs. 2024;42(1):60–9.38071684 10.1007/s10637-023-01413-z
38. Therneau T Crowson C Atkinson E Using time dependent covariates and time dependent coefficients in the Cox model Surv Vignettes 2017 2 3 1 25
Therneau T, Crowson C, Atkinson E. Using time dependent covariates and time dependent coefficients in the Cox model. Surv Vignettes. 2017;2(3):1–25.
39. Hamada T Nakai Y Isayama H Progression-free survival as a surrogate for overall survival in first-line chemotherapy for advanced pancreatic cancer Eur J Cancer 2016 65 11 20 10.1016/j.ejca.2016.05.016 27451020
Hamada T, Nakai Y, Isayama H, et al. Progression-free survival as a surrogate for overall survival in first-line chemotherapy for advanced pancreatic cancer. Eur J Cancer. 2016;65:11–20.27451020 10.1016/j.ejca.2016.05.016
40. Buyse M Burzykowski T Carroll K Progression-free survival is a surrogate for survival in advanced colorectal cancer J Clin Oncol 2007 25 33 5218 5224 10.1200/JCO.2007.11.8836 18024867
Buyse M, Burzykowski T, Carroll K, et al. Progression-free survival is a surrogate for survival in advanced colorectal cancer. J Clin Oncol. 2007;25(33):5218–24.18024867 10.1200/JCO.2007.11.8836
41. Ma Y. Flexible isotonic regression in survival data analysis. PhD thesis, The George Washington University, 2010.
42. Hanagal DD Modeling survival data using frailty models 2011 Berlin Springer
Hanagal DD. Modeling survival data using frailty models. Berlin: Springer; 2011.
43. Yashin AI Vaupel JW Iachine IA Correlated individual frailty: an advantageous approach to survival analysis of bivariate data Math Popul Stud 1995 5 2 145 159 10.1080/08898489509525394 12290053
Yashin AI, Vaupel JW, Iachine IA. Correlated individual frailty: an advantageous approach to survival analysis of bivariate data. Math Popul Stud. 1995;5(2):145–59.12290053 10.1080/08898489509525394
44. Wienke A. Frailty models in survival analysis. Chapman and Hall/CRC, 2010.
45. Martins A Aerts M Hens N Correlated gamma frailty models for bivariate survival time data Stat Methods Med Res 2019 28 10–11 3437 3450 10.1177/0962280218803127 30319043
Martins A, Aerts M, Hens N, et al. Correlated gamma frailty models for bivariate survival time data. Stat Methods Med Res. 2019;28(10–11):3437–50.30319043 10.1177/0962280218803127
46. Ng SK Tawiah R Mclachlan GJ Joint frailty modeling of time-to-event data to elicit the evolution pathway of events: a generalized linear mixed model approach Biostatistics 2021 1 1
Ng SK, Tawiah R, Mclachlan GJ, et al. Joint frailty modeling of time-to-event data to elicit the evolution pathway of events: a generalized linear mixed model approach. Biostatistics. 2021;1:1.
47. Romeo JS Meyer R Gallardo DI Bayesian bivariate survival analysis using the power variance function copula Lifetime Data Anal 2018 24 2 355 383 10.1007/s10985-017-9396-1 28536818
Romeo JS, Meyer R, Gallardo DI. Bayesian bivariate survival analysis using the power variance function copula. Lifetime Data Anal. 2018;24(2):355–83.28536818 10.1007/s10985-017-9396-1
48. Emura T, Matsui S, Rondeau V. Survival analysis with correlated endpoints: Joint Frailty-Copula models. Springer; 2019.
49. Sofeu CL Emura T Rondeau V A joint frailty-copula model for meta-analytic validation of failure time surrogate endpoints in clinical trials Biom J 2021 63 2 423 446 10.1002/bimj.201900306 33006170
Sofeu CL, Emura T, Rondeau V. A joint frailty-copula model for meta-analytic validation of failure time surrogate endpoints in clinical trials. Biom J. 2021;63(2):423–46.33006170 10.1002/bimj.201900306
50. Pan W A multiple imputation approach to Cox regression with interval-censored data Biometrics 2000 56 1 199 203 10.1111/j.0006-341X.2000.00199.x 10783796
Pan W. A multiple imputation approach to Cox regression with interval-censored data. Biometrics. 2000;56(1):199–203.10783796 10.1111/j.0006-341X.2000.00199.x
51. Lam KF Xu Y Cheung TL A multiple imputation approach for clustered interval-censored survival data Stat Med 2010 29 6 680 693 10.1002/sim.3835 20069624
Lam KF, Xu Y, Cheung TL. A multiple imputation approach for clustered interval-censored survival data. Stat Med. 2010;29(6):680–93.20069624 10.1002/sim.3835
52. Goggins WB Finkelstein DM A proportional hazards model for multivariate interval-censored failure time data Biometrics 2000 56 3 940 943 10.1111/j.0006-341X.2000.00940.x 10985240
Goggins WB, Finkelstein DM. A proportional hazards model for multivariate interval-censored failure time data. Biometrics. 2000;56(3):940–3.10985240 10.1111/j.0006-341X.2000.00940.x
53. Kim MY Xue X The analysis of multivariate interval-censored survival data Stat Med 2002 21 23 3715 3726 10.1002/sim.1265 12436466
Kim MY, Xue X. The analysis of multivariate interval-censored survival data. Stat Med. 2002;21(23):3715–26.12436466 10.1002/sim.1265
54. Bellamy SL Li Y Ryan LM Analysis of clustered and interval censored data from a community-based study in asthma Stat Med 2004 23 23 3607 3621 10.1002/sim.1918 15534894
Bellamy SL, Li Y, Ryan LM, et al. Analysis of clustered and interval censored data from a community-based study in asthma. Stat Med. 2004;23(23):3607–21.15534894 10.1002/sim.1918
55. Wong MCM Lam KF Lo ECM Bayesian analysis of clustered interval-censored data J Dent Res 2005 84 9 817 821 10.1177/154405910508400907 16109990
Wong MCM, Lam KF, Lo ECM. Bayesian analysis of clustered interval-censored data. J Dent Res. 2005;84(9):817–21.16109990 10.1177/154405910508400907
