The CPSM R-package is an advanced computational pipeline designed to predict the survival probability of cancer patients with precision and efficiency. This package automates critical steps in survival analysis, including data preprocessing, normalization, and splitting datasets into training and testing subsets. It identifies significant features using univariate survival analysis and LASSO COX-regression, and also generates LASSO-based prognostic index (PI) scores to optimize predictive performance. CPSM builds robust survival prediction models using selected clinical, molecular, and integrated feature sets. To support data interpretation and clinical applications, the package includes powerful visualization tools, such as survival curves, bar plots for predicted mean and median survival times, and nomograms. Designed for multi-omics data, CPSM simplifies complex workflows, empowering researchers to uncover novel biomarkers and advance precision oncology.
Figure: The workflow of the CPSM package represents different steps performed by various functions of the CPSM package.
#Step1: First Install remote package
install.packages("remotes")
#load remotes package
library("remotes")
#Step2: install CPSM package
remotes::install_github("hks5august/CPSM", local = TRUE, , dependencies=TRUE)
#or use the following command
#remotes::install_github("hks5august/CPSM", ref = "v1.0.0", dependencies=TRUE)
# Check if package get installed, load package
library("CPSM") #Load CPSM packages
library(CPSM)
#Load other required packages
library(preprocessCore)
library(ggfortify)
library(survival)
library(survminer)
library(ggplot2)
library(MASS)
library(MTLR)
library(dplyr)
library(SurvMetrics)
library(pec)
library(glmnet)
library(reshape2)
library(rms)
library(Matrix)
library(Hmisc)
library(survivalROC)
library(ROCR)
#set seed
set.seed(7)
Example Input data: "Example_TCGA_LGG_FPKM_data" is a tab separated file. It contains Samples (184 LGG Cancer Samples) in the rows and Features in the columns. Gene Expression is available in terms of FPKM values in the data. Features information: In the data there are 11 clinical + demographic, 4 types survival with time and event information and 19,978 protein coding genes. Clinical and demographic features: Clinical demographic features that are present in this example data include Age, subtype, gender, race, ajcc_pathologic_tumor_stage, histological_type, histological_grade, treatment_outcome_first_course, radiation_treatment_adjuvant, sample_type, type. Types of Survival: 4 types of Survival include OS (overall survival), PFS (progression-free survival), DSS (disease-specific survival), DFS (Disease-free survival). In the data, column names OS, PFS, DSS and DFS represent event information, while OS.time, PFS.time, DSS.time and DFS.time indicate survival time in days.
This function converts OS time (in days) into months and then removes samples where OS/OS.time information is missing. Here, we need to provide input data in tsv or txt format. Further, we needs to define col_num (column number at which clinical/demographic and survival information ends,e.g. 20, surv_time (name of column which contain survival time (in days) information, e.g. OS.time ) and output file name, e.g. “New_data.txt”
library(SummarizedExperiment)
#load data
data(Example_TCGA_LGG_FPKM_data, package="CPSM")
#check data
Example_TCGA_LGG_FPKM_data
#preprpocess data
New_data <- data_process_f(assays(Example_TCGA_LGG_FPKM_data)$expression,
col_num=20, surv_time = "OS.time")
str(New_data[1:10])
After data processing, we got a new output file “New_data”, which contains 176 samples. Thus, data_process_f function removes 8 samples where OS/OS time information is missing. Besides, here is a new 21st column in the data with column name “OS_month” where OS time is available in months.
Before proceeding further, we need to split our data into training and test subset for the purpose of feature selection and model development. Here, we need output from the previous step as an input ( which was “New_data.txt”). Next we need to define the fraction (e.g. 0.9) by which we want to split data into training and test. Thus, fraction=0.9 will split data into 90% training and 10% as test set. Besides, we also need to provide training and set output names (e.g. train_FPKM.txt,test_FPKM.txt )
data(New_data, package = "CPSM")
# Call the function
result <- tr_test_f(data = assays(New_data)$expression, fraction = 0.9)
# Access the train and test data
train_FPKM <- result$train_data
str(train_FPKM[1:10])
test_FPKM <- result$test_data
str(test_FPKM[1:10])
After the train-test split, we got a two new outputs: “train_FPKM”, “test_FPKM”, where, train_FPKM contains 158 samples and test_FPKM contains 18 samples. Thus, tr_test_f function splits data into 90:10 ratio.
Next to select features and develop ML models, data must be normalized. Since, expression is available in terms of FPKM values. Thus, train_test_normalization_f function will first convert FPKM value into log scale [log2(FPKM+1) followed by quantile normalization using the “preprocessCore” package. Here, training data will be used as a target matrix for quantile normalization. Here, we need to provide training and test datasets (that we obtained from the previous step of Train/Test Split). Further, we need to provide column number where clinical information ends (e.g. 21) in the input datasets. Besides, we also need to provide output files names (train_clin_data (which contains only Clinical information of training data), test_clin_data (which contains only Clinical information of training data), train_Normalized_data_clin_data (which contains Clinical information and normalized values of genes of training samples), test_Normalized_data_clin_data (which contains Clinical information and normalized values of genes of test samples).
# Normalize the training and test data sets
data(train_FPKM, package = "CPSM")
data(test_FPKM, package = "CPSM")
Result_N_data <- train_test_normalization_f(train_data = train_FPKM,
test_data = test_FPKM,
col_num = 21)
# Access the Normalized train and test data
Train_Clin <- Result_N_data$Train_Clin
Test_Clin <- Result_N_data$Test_Clin
Train_Norm_data <- Result_N_data$Train_Norm_data
Test_Norm_data <- Result_N_data$Test_Norm_data
str(Train_Clin[1:10])
str(Train_Norm_data[1:10])
After, running the function, we obtained 4 outputs: Train_Clin - Contains only Clinical features, Test_Clin - contains only Clinical features of Test samples; Train_Norm_data - Clinical features with normalized values of genes for training samples; Test_Norm_data - Clinical features with normalized values of genes for test samples.
Next to create a survival model, we will create a Prognostic Index (PI) Score. PI score is calculated based on the expression of the features selected by the LASSO regression model and their beta coefficients. For instance, 5 features (G1, G2, G3, G4, and G5 and their coefficient values are B1, B2, B3, B4, and B5, respectively) selected by the LASSO method. Then PI score will be computed as following:
PI score = G1B1 + G2B2 + G3 * B3 + G4B4+ G5B5
Here, we need to provide Normalized training (Train_Norm_data) and test data (Test_Norm_data)as input data that we have obtained from the previous function “train_test_normalization_f”. Further, we need to provide col_num n column number at which clinical features ends (e.g. 21), nfolds (number of folds e.g. 5) for the LASSO regression method to select features. We implemented LASSO using the “glmnet” package. Further, we need to provide surv_time (name of column containing survival time in months, e.g. OS_month) and surv_event (name of column containing survival event information, e.g. OS) information in the data. Besides, we also need to provide names and training and test output file names to store data containing LASSO genes and PI values.
# Step 4 - Lasso PI Score
data(Train_Norm_data, package = "CPSM")
data(Test_Norm_data, package = "CPSM")
Result_PI <- Lasso_PI_scores_f(train_data = Train_Norm_data,
test_data = Test_Norm_data,
nfolds=5,
col_num=21,
surv_time = "OS_month",
surv_event = "OS")
# Extract list of selected features with beta coeff values
Train_Lasso_key_variables <- Result_PI$Train_Lasso_key_variables
# Train data with PI
Train_PI_data <- Result_PI$Train_PI_data
# Test data with PI score
Test_PI_data <- Result_PI$Test_PI_data
#view train and test data structures
str(Train_PI_data[1:10])
str(Test_PI_data[1:10])
# Lambda plot
plot(Result_PI$cvfit)
#Coefficient regression plot
plot(Result_PI$cvfit$glmnet.fit, "lambda", label=TRUE)
Thus, Lasso_PI_scores_f gave us following outputs:
- Train_Lasso_key_variables: List of features selected by LASSO and their beta coefficient values
- Train_Cox_Lasso_Regression_lamda_plot: Lasso Regression Lambda plot.
- Train_PI_data: It contains expression of genes selected by LASSO and PI score in the last column for training samples.
- Test_PI_data: It contains expression of genes selected by LASSO and PI score in the last column for test samples.
Figure: Lasso Regression Lambda plot
Besides PI score, with the “Univariate_sig_features_f” function of CPSM package, we can select significant (p-value <0.05) features based on univariate survival analysis. These features are selected based on their capability to stratify high-risk and low-risk survival groups using the cut off value of their median expression.
Here, we need to provide Normalized training (Train_Norm_data.txt) and test data (Test_Norm_data.txt)as input data that we have obtained from the previous function “train_test_normalization_f”. Further, we need to provide a “col_num” (e.g 21)column number at which clinical features ends. Further, we need to provide surv_time (name of column containing survival time in months, e.g. OS_month) and surv_event (name of column containing survival event information, e.g. OS) information in the data. Besides, we also need to provide names and training and test output file names to store data containing expression of selected genes.
#Step 4b - Univariate Survival Significant Feature Selection.
# load train and test data
data(Train_Norm_data, package = "CPSM")
data(Test_Norm_data, package = "CPSM")
Result_Uni <- Univariate_sig_features_f(train_data = Train_Norm_data,
test_data = Test_Norm_data,
col_num=21,
surv_time = "OS_month" ,
surv_event = "OS")
#significant genes with HR, P-val
Univariate_Survival_Significant_genes_List <- Result_Uni$Univariate_Survival_Significant_genes_List
#training data with significant genes
Train_Uni_sig_data <- Result_Uni$Train_Uni_sig_data
#test data with significant genes
Test_Uni_sig_data <- Result_Uni$Test_Uni_sig_data
str(Univariate_Survival_Significant_genes_List[1:10])
#significant clinical features with HR, P-val
Univariate_Survival_Significant_clin_List <- Result_Uni$Univariate_Survival_Significant_clin_List
Train_Uni_sig_clin_data <- Result_Uni$Train_Uni_sig_clin_data
Test_Uni_sig_clin_data <- Result_Uni$Test_Uni_sig_clin_data
str(Univariate_Survival_Significant_clin_List)
Thus, Univariate_sig_features_f gave us following outputs: Univariate_Survival_Significant_genes_List: a table of univariate significant genes along with their corresponding coefficient values, HR value, P-values, C-Index values. Train_Uni_sig_data: It contains expression of significant genes selected by univariate survival analysis for training samples. Test_Uni_sig_data: It contains expression of significant genes selected by univariate survival analysis for test samples.
After selecting significant or key features using LASSO or Univariate survival analysis, next we want to develop an ML prediction model to predict survival probability of patients. MTLR_pred_model_f function of CPSM give us multiple options to develop models including Only Clinical features (Model_type=1), PI score (Model_type=2), PI Score + Clinical features (Model_type=3), Significant Univariate features (Model_type=4), Significant Univariate features Clinical features (Model_type=5) using MTLR package. Further, here, we were interested in developing a model based on PI score. Thus, we need to provide following inputs: (1) Training data with only clinical features, (2) Test data with only clinical features, (3) Model type (e.g. 2, since we want to develop model based on PI score), (4) Training data with PI score , (5) Test data with PI score, (6) Clin_Feature_List (e.g. Key_PI_list.txt), a list of features which will be used to build model . Furthermore, we also need to provide surv_time (name of column containing survival time in months, e.g. OS_month) and surv_event (name of column containing survival event information, e.g. OS) information in the clinical data
data(Train_Clin, package = "CPSM")
data(Test_Clin, package = "CPSM")
data(Key_Clin_feature_list, package = "CPSM")
Result_Model_Type1 <- MTLR_pred_model_f(train_clin_data = Train_Clin,
test_clin_data = Test_Clin,
Model_type = 1,
train_features_data = Train_Clin,
test_features_data = Test_Clin,
Clin_Feature_List = Key_Clin_feature_list,
surv_time = "OS_month",
surv_event = "OS")
survCurves_data <- Result_Model_Type1$survCurves_data
mean_median_survival_time_data <- Result_Model_Type1$mean_median_survival_time_data
survival_result_based_on_MTLR <- Result_Model_Type1$survival_result_based_on_MTLR
Error_mat_for_Model <- Result_Model_Type1$Error_mat_for_Model
data(Train_Clin, package = "CPSM")
data(Test_Clin, package = "CPSM")
data(Train_PI_data, package = "CPSM")
data(Test_PI_data, package = "CPSM")
data(Key_PI_list, package = "CPSM")
Result_Model_Type2 <- MTLR_pred_model_f(train_clin_data = Train_Clin,
test_clin_data = Test_Clin,
Model_type = 2,
train_features_data = Train_PI_data ,
test_features_data = Test_PI_data ,
Clin_Feature_List = Key_PI_list,
surv_time = "OS_month",
surv_event = "OS")
survCurves_data <- Result_Model_Type2$survCurves_data
mean_median_survival_time_data <- Result_Model_Type2$mean_median_survival_time_data
survival_result_based_on_MTLR <- Result_Model_Type2$survival_result_based_on_MTLR
Error_mat_for_Model <- Result_Model_Type2$Error_mat_for_Model
data(Train_Clin, package = "CPSM")
data(Test_Clin, package = "CPSM")
data(Train_PI_data, package = "CPSM")
data(Test_PI_data, package = "CPSM")
data(Key_Clin_features_with_PI_list, package = "CPSM")
Result_Model_Type3 <- MTLR_pred_model_f(train_clin_data = Train_Clin,
test_clin_data = Test_Clin,
Model_type = 3,
train_features_data = Train_PI_data,
test_features_data = Test_PI_data,
Clin_Feature_List = Key_Clin_features_with_PI_list,
surv_time = "OS_month",
surv_event = "OS")
survCurves_data <- Result_Model_Type3$survCurves_data
mean_median_survival_time_data <- Result_Model_Type3$mean_median_survival_time_data
survival_result_based_on_MTLR <- Result_Model_Type3$survival_result_based_on_MTLR
Error_mat_for_Model <- Result_Model_Type3$Error_mat_for_Model
data(Train_Clin, package = "CPSM")
data(Test_Clin, package = "CPSM")
data(Train_Uni_sig_data, package = "CPSM")
data(Test_Uni_sig_data, package = "CPSM")
data(Key_univariate_features_with_Clin_list, package = "CPSM")
Result_Model_Type5 <- MTLR_pred_model_f(train_clin_data = Train_Clin,
test_clin_data = Test_Clin,
Model_type = 4,
train_features_data = Train_Uni_sig_data,
test_features_data = Test_Uni_sig_data,
Clin_Feature_List =Key_univariate_features_with_Clin_list,
surv_time = "OS_month",
surv_event = "OS")
survCurves_data <- Result_Model_Type5$survCurves_data
mean_median_survival_time_data <- Result_Model_Type5$mean_median_survival_time_data
survival_result_based_on_MTLR <- Result_Model_Type5$survival_result_based_on_MTLR
Error_mat_for_Model <- Result_Model_Type5$Error_mat_for_Model
After, implementing MTLR_pred_model_f function , we got following outputs:
- Model_with_PI.RData : Model on training data
- survCurves_data : Table containing predicted survival probability of each patient at different time points. This data can be further used to plot the survival curve of patients.
- mean_median_survival_time_data : Table containing predicted mean and median survival time of each patient in the test data. This data can be further used for bar plots.
- Error_mat_for_Model : Table containing performance parameters obtained on test data based on prediction model. It contains IBS score (Integrated Brier Score) =0.192, C-Index =0.81.
Next to visualize survival of patients, we will plot survival curve plots using the surv_curve_plots_f function based on the data “survCurves_data ” that we obtained from the previous step after running the MTLR_pred_model_f function. Further, the surv_curve_plots_f function also allows highlighting a specific patient on the curve. Thus the function needs only two inputs: 1) Surv_curve_data, (2) Sample ID of a specific patient (e.g. TCGA-DB-A4XF-01) that needs to be highlighted.
#Create Survival curves/plots for individual patients
data(survCurves_data, package = "CPSM")
plots <- surv_curve_plots_f(Surv_curve_data = survCurves_data,
selected_sample = "TCGA-TQ-A7RQ-01")
# Print the plots
# Survival curves for all patients in test data
print(plots$all_patients_plot)
# Survival curve for selected patient in test data
print(plots$highlighted_patient_plot)
Here, we obtained two output plots:
- Survival curves for all patients in the test data with different colors
- Survival curves for all patients (in black) and highlighted patient (red) in the test data
Figure: Survival curves for all patients in the test data.
Figure: Survival curves for all patients (in black/grey) and highlighted patient (red).
Next, to visualize the predicted survival time of patients, we will plot the barplot for mean/median using “mean_median_surv_barplot_f” function based on the data that we obtained from step 5 after running the MTLR_pred_model_f function. Further, the mean_median_surv_barplot_f function also allows highlighting a specific patient on the curve. Thus the function needs only two inputs: 1) surv_mean_med_data, (2) Sample ID of a specific patient (e.g. TCGA-DB-A4XF-01) that needs to be highlighted.
data(mean_median_survival_time_data, package = "CPSM")
plots_2 <- mean_median_surv_barplot_f(surv_mean_med_data =
mean_median_survival_time_data,
selected_sample = "TCGA-TQ-A7RQ-01")
# Print the plots
# Barplot for mean and median survival time of all patients in Test data
print(plots_2$mean_med_all_pat)
# Barplot for mean and median survival time of all patients (grey) in Test data and selected patient (colored)
print(plots_2$highlighted_selected_pat)
Here, we obtained two output plots:
- Barplot for all patients in the test data, where the red color bar represents mean survival and cyan/green color bar represents median survival time.
- Barplot for all patients with a highlighted patient (dashed black outline) in the test data. It shows this patient has a predicted mean and median survival is 81.58 and 75.50 months.
Figure: Barplot for mean and median survival time of all patients in the data, where the red color bar represents mean survival and cyan/green color bar represents median survival time.
Figure: Barplot for mean and median survival time of all patients (grey color) with a highlighted patient (colored) in the test data.
Next, the Nomogram_generate_f function of CPSM will allow users to generate a nomogram plot for their data (training data containing all samples) based on user-defined clinical and other features in their data. For instance, we will generate a nomogram based on 6 features (Age, gender, race, histological_type, sample_type, PI). Here, we will provide data containing all the features (Samples in rows and features in columns) (e.g. Train_Data_Nomogram_input) and a list of features (feature_list_for_Nomogram) based on which we want to generate a nomogram. Further, we also need to provide surv_time (name of column containing survival time in months, e.g. OS_month) and surv_event (name of column containing survival event information, e.g. OS) information in the data.
data(Train_Data_Nomogram_input, package = "CPSM")
data(feature_list_for_Nomogram, package = "CPSM")
Result_Nomogram <- Nomogram_generate_f(data = Train_Data_Nomogram_input,
Feature_List = feature_list_for_Nomogram,
surv_time = "OS_month",
surv_event = "OS")
C_index_mat <- Result_Nomogram$C_index_mat
Here, we will get a Nomogram based on features that we provide. This nomogram can predict Risk (Event risk, eg, Death), 1-year, 3-year, 5-year and 10 years survival of patients.
As last part of this document, we call the function "sessionInfo()", which reports the version numbers of R and all the packages used in this session. It is good practice to always keep such a record as it will help to trace down what has happened in case that an R script ceases to work because the functions have been changed in a newer version of a package.
sessionInfo()
CPSM: R-package of an Automated Machine Learning Pipeline for Predicting the Survival Probability of Single Cancer Patient. (2024) Harpreet Kaur, Pijush Das, Kevin Camphausen, Uma T Shankavaram. bioRxiv 2024.11.14.623597; doi: https://doi.org/10.1101/2024.11.14.623597
- Kuhn, Max (2008). “Building Predictive Models in R Using the caret Package.” Journal of Statistical Software, 28(5), 1–26. doi:10.18637/jss.v028.i05, https://www.jstatsoft.org/index.php/jss/article/view/v028i05.
- Bolstad B (2024). preprocessCore: A collection of pre-processing functions. R package version 1.66.0, https://github.com/bmbolstad/preprocessCore.
- Horikoshi M, Tang Y (2018). ggfortify: Data Visualization Tools for Statistical Analysis Results. https://CRAN.R-project.org/package=ggfortify.
- Therneau T (2024). A Package for Survival Analysis in R. R package version 3.7-0, https://CRAN.R-project.org/package=survival.
- Terry M. Therneau, Patricia M. Grambsch (2000). Modeling Survival Data: Extending the Cox Model. Springer, New York. ISBN 0-387-98784-3.
- Kassambara, A., Kosinski, M., Biecek, P., & Scheipl, F. (2021). survminer: Drawing survival curves using 'ggplot2' (Version 0.4.9) [R package]. CRAN. https://doi.org/10.32614/CRAN.package.survminer
- Haider, H. (2019). MTLR: Survival Prediction with Multi-Task Logistic Regression (Version 0.2.1) [R package]. CRAN. https://doi.org/10.32614/CRAN.package.MTLR
- Wickham H, François R, Henry L, Müller K, Vaughan D (2023). dplyr: A Grammar of Data Manipulation. R package version 1.1.4, https://github.com/tidyverse/dplyr, https://dplyr.tidyverse.org.
- Wickham H (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4, https://ggplot2.tidyverse.org.
- Zhou, H., Cheng, X., Wang, S., Zou, Y., & Wang, H. (2022). SurvMetrics: Predictive Evaluation Metrics in Survival Analysis (Version 0.5.0) [R package]. CRAN. https://doi.org/10.32614/CRAN.package.SurvMetrics
- Simon N, Friedman J, Tibshirani R, Hastie T (2011). “Regularization Paths for Cox's Proportional Hazards Model via Coordinate Descent.” Journal of Statistical Software, 39(5), 1–13. doi:10.18637/jss.v039.i05.
- Gerds TA (2023). pec: Prediction Error Curves for Risk Prediction Models in Survival Analysis. R package version 2023.04.12, https://CRAN.R-project.org/package=pec.
- Heagerty, P. J., & Saha-Chaudhuri, P. (2022). survivalROC: Time-Dependent ROC Curve Estimation from Censored Survival Data (Version 1.0.3.1) [R package]. CRAN. https://doi.org/10.32614/CRAN.package.survivalROC
- Harrell, F. E. Jr. (2024). rms: Regression Modeling Strategies (Version 6.8-1) [R package]. CRAN. https://doi.org/10.32614/CRAN.package.rms
- Sing T, Sander O, Beerenwinkel N, Lengauer T (2005). “ROCR: visualizing classifier performance in R.” Bioinformatics, 21(20), 7881. http://rocr.bioinf.mpi-sb.mpg.de.
- Bates, D., Maechler, M., Jagan, M., Davis, T. A., Karypis, G., Riedy, J., Oehlschlägel, J., & R Core Team. (2024). Matrix: Sparse and Dense Matrix Classes and Methods (Version 1.7-0) [R package]. CRAN. https://doi.org/10.32614/CRAN.package.Matrix
- Harrell, F. E. Jr., & Dupont, C. (2024). Hmisc: Harrell Miscellaneous (Version 5.1-3) [R package]. CRAN. https://doi.org/10.32614/CRAN.package.Hmisc Wickham H (2007). “Reshaping Data with the reshape Package.” Journal of Statistical Software, 21(12), 1–20. http://www.jstatsoft.org/v21/i12/.