‘create_se_object’

Expression and survival data are stored in a SummarizedExperiment object. The function ‘create_se_object’ reads a file with expression data and a file with survival data, and returns a SummarizedExperiment object.

# Read names of the built-in gene expression data file and survival data file
fdat = system.file("extdata", "example_genes.txt", package = "kmcut")
sdat = system.file("extdata", "survival_data.txt", package = "kmcut")

# Create a SummarizedExperiment object 'se'
se = create_se_object(efile = fdat, sfile = sdat)

‘km_opt_pcut’

This function uses each distinct value of every feature from the dataset as a stratification cutoff to select a cutoff that results in the maximum separation of the Kaplan-Meier survival curves, and then estimates the statistical significance of this optimal cutoff by means of the permutation test. A detailed description of the steps implemented in this function is provided below and is also available in our original publication (Wei et al, 2018).

First, an ordered list C of all distinct values of a given feature observed in the group of test subjects (samples) is created, all values in the list being sorted from smallest to largest. Then, each value from the list is used as a stratification cutoff: samples with the feature below or equal to the cutoff are labeled as ‘low’ and above the cutoff as ‘high’. To avoid edge effects, if for a particular cutoff the size of low or high sub-group is smaller than ‘min_fraction’ of the total number of samples, this cutoff is discarded from the list. The log-rank test is applied to compare survival distributions between the low and high groups for each stratification cutoff, and the value of the test statistic is recorded. After all stratification cutoffs from the list are tested, the vector of the observed values of the test statistic is created, O=(o₁, o₂ ,…, o_n), where o_i is the observed value of the test statistic for the cutoff c_i in the list C=(c₁, c₂ ,…, c_n). The cutoff that results in the largest value of the test statistic, O_MAX, is selected as the optimal cutoff, C_OPT. If two or more cutoffs result in the same observed value of the test statistic, the cutoff closest to the median is selected.

Additionally, the shape of the plot of the observed values of the test statistic is compared to the expected “ideal” plot. The empirical assumption behind this comparison is that the “ideal” optimization should produce a plot of the observed values of the test statistic with exactly one peak, meaning that the values monotonically increase before the peak and monotonically decrease after the peak. Initially, three possible expected plots are defined: E₁, E₂, and E₃. Each plot consists of cutoff points from C, has a single peak, with the height of this peak being equal to O_MAX. In the plot E₁ the peak is located at the 25th percentile of all distinct values of the feature, in the plot E₂ at the 50th percentile, and in the plot E₃ at the 75th percentile. The final expected plot, E, is selected from (E₁, E₂, E₃) to minimize the distance from the location of its peak to C_OPT. The similarity of is quantified by the Spearman rank correlation between the vectors of the observed and expected values, r(O, E).

The statistical significance of the stratification cutoff C_OPT obtained from the optimization procedure for a given feature is estimated by means of a random permutation test run for ‘n_iter’ iterations. On each iteration i, the sample labels in the survival data and in the feature data are randomly shuffled, the same optimization procedure is applied to the randomized data, and the largest value of the observed test statistic O_RAND(i) and the Spearman rank correlation between the observed and expected plot r_RAND(i) are recorded. After all ‘n_iter’ iterations are completed, the p-value is calculated as follows:

The equation used to calculate the p-value.

where d[O_RAND(i) O_MAX AND r_RAND(i) r(O, E)] = 1 if on random iteration i the largest value of the observed test statistic is equal to or greater than O_MAX and the Spearman rank correlation between the observed and the expected plot is equal to or greater than r(O, E). Otherwise,
d[O_RAND(i) O_MAX AND r_RAND(i) r(O, E)] = 0.

The random permutation test described above is time-consuming. It is recommended that the users first run the ‘km_opt_scut’ function to identify a relatively small sub-set of candidate genes, and then apply ‘km_opt_pcut’ to this sub-set.

The ‘km_opt_pcut’ function can be run in parallel on multiple processors by setting the nproc argument to the desired number of processors. The graph below shows how the run time of the ‘km_opt_pcut’ function decreases when the number of processors is increased. The results were obtained by using a data set of 200 genes and 215 samples, with 10000 iterations per gene. Benchmarking utilized the NIH High Performance Computing cluster Biowulf.

The run time vs. the number of processors.

An example of how to run ‘km_opt_pcut’ on 1 processor with the data files included in the package:

# Read names of the built-in gene expression data file and survival data file
fdat = system.file("extdata", "example_genes.txt", package = "kmcut")
sdat = system.file("extdata", "survival_data.txt", package = "kmcut")

# Create SummarizedExperiment object
se = create_se_object(efile = fdat, sfile = sdat)

# Run the permutation test for 10 iterations for each gene, use 1 processor
km_opt_pcut(obj = se, bfname = "test", n_iter = 10, wlabels = TRUE,
                wpdf = FALSE, verbose = FALSE, nproc = 1)

## Running on 1 CPU(s)

Two graphs are created for each gene, MYCN and MYH2.

One graph shows the observed optimization plot (blue circles) and the expected optimization plot (green triangles). The optimal stratification cutoff is highlighted by the red circle. The Spearman rank correlation between the observed and the expected plots is printed after gene name. For instance, in the case of MYCN the correlation is high, R = 0.909, indicating good optimization. In the case of MYH2, the correlation is low and negative, R = -0.043, indicating poor optimization.
The other graph shows the Kaplan-Meier survival curves for groups with low and high gene expression groups stratified using the optimal cutoff. The value of this optimal cutoff and the p-value are printed after gene name. For instance, in the case of MYCN gene, the optimal cutoff is 5.84464 and the p-value is 0 (note that in the figure the p-value is shown as exactly 0 because only 10 iterations of the permutation test were used in this example).

The MYCN oncogene expression is known to be a strong predictor of survival outcome for neuroblastoma patients - low expression levels correspond to better survival, high expression levels correspond to poor survival (Norris et al, 1997). The significant results of the optimization for MYCN confirm this information.

The MYH2 gene encodes the myosin heavy chain 2, which is a protein found in the muscle tissue and the level of its expression has nothing to do with neuroblastoma survival (Smerdu et al, 1994). The statistically not significant results of the optimization for MYH2 confirm this information.

Additionally, this run will create the following two output files in the current working directory (the names of output files will be created automatically by adding the run information to the base name ‘bfname’).

a) Tab-delimited text file with the results
“test_KMoptp_minf_0.10_iter_10.txt”

Table 3. An illustration of the file format for tab-delimited text file with the results for two genes. 1st column – gene id, 2nd – optimal stratification cutoff, 3rd – test statistic calculated for the optimal cutoff, 4th – number of samples in low-expression sub-group, 5th - number of samples in high-expression sub-group, 6th – permutation test p-value, 7th – FDR-adjusted p-value.

b) CSV file with low/high sample labels
“test_KMoptp_minf_0.10_iter_10_labels.csv”

Table 4. An illustration of the file format for the CVS file with low/high labels for two genes and five samples. 1st column – sample ids, all subsequent columns contain low/high labels for each gene, where 1 and 2 correspond to low- and high-expression sub-groups, respectively (‘low’ means below the cutoff and ‘high’ means above the cutoff).

‘km_opt_scut’

This function uses each distinct value of a given feature observed in the dataset as a stratification cutoff to select a cutoff that results in the maximum separation of the Kaplan-Meier survival curves, but does not use the permutation test to estimate the statistical significance of this optimal cutoff. ‘km_opt_scut’ produces graphs and output files virtually identical to the output of ‘km_opt_pcut’ (except for the permutation test p-value) and is meant to be used as a fast exploratory alternative to ‘km_opt_pcut’. An example of how to use ‘km_opt_scut’ with the data files included in the package:

# Read names of the built-in gene expression data file and survival data file
fdat = system.file("extdata", "example_genes.txt", package = "kmcut")
sdat = system.file("extdata", "survival_data.txt", package = "kmcut")

# Create SummarizedExperiment object
se <- create_se_object(efile = fdat, sfile = sdat)
 
# Search for optimal cutoffs
km_opt_scut(obj = se, bfname = "test", wpdf = FALSE, verbose = FALSE)

Two graphs are created for each gene, MYCN and MYH2.

Additionally, this run will create the following two output files in the current working directory (the names of output files will be created automatically by adding the run information to the base name ‘bfname’).

a) Tab-delimited text file with the results “test_KMopt_minf_0.10.txt

b) CSV file with low/high sample labels “test_KMopt_minf_0.10_labels.csv”

The format of these files is identical to the output files described for function ‘km_opt_pcut’. The only difference is that the tab-delimited file with the results in column ‘P’ and the PDF file with Kaplan-Meier curves contain the p-value obtained from the log-rank test performed with the optimal cutoff. This p-value is provided for information purposes only and should not be treated as an actual p-value because it does not reflect the calculation of multiple test statistics involved in the optimization procedure applied to select the optimal cutoff.

‘km_qcut’

For each feature, this function uses the cutoff supplied as a quantile (from 0 to 100) to stratify samples into 2 groups (below and above this quantile), plots Kaplan-Meier survival curves for these two groups and calculates the log-rank test p-value. Since for each feature only one cutoff corresponding to the specified quantile is used, it is equivalent to performing one log-rank test per feature (that is, no optimization is performed). An example of how to use ‘km_qcut’ with the data files included in the package:

# Read names of the built-in gene expression data file and survival data file
fdat = system.file("extdata", "example_genes.txt", package = "kmcut")
sdat = system.file("extdata", "survival_data.txt", package="kmcut")

# Create SummarizedExperiment object
se <- create_se_object(efile = fdat, sfile = sdat)
 
# Use the 50th quantile (the median) to stratify the samples
km_qcut(obj = se, bfname = "test", quant = 50, wpdf = FALSE)

A graph with Kaplan-Meier curves will be created for each gene, MYCN and MYH2.

Additionally, this run will create the following two output files in the current working directory (the names of output files will be created automatically by adding the run information to the base name ‘bfname’).

a) Tab-delimited text file with the results “test_KM_quant_50.txt”

b) CSV file with low/high sample labels “test_KM_quant_50_labels.csv”

The format of these two files is identical to the output files described for function ‘km_opt_pcut’.

‘km_ucut’

For each feature, this function applies the user-supplied fixed value as a cutoff to stratify samples into 2 groups (below/above this cutoff), plots Kaplan-Meier survival curves for these two groups, and calculates the log-rank test p-value. Since for each feature the same fixed cutoff is used, it is equivalent to performing one log-rank test per feature (no optimization is performed). An example of how to use ‘km_ucut’ with the data files included in the package:

# Read names of the built-in gene expression data file and survival data file
fdat = system.file("extdata", "example_genes.txt", package = "kmcut")
sdat = system.file("extdata", "survival_data.txt", package = "kmcut")

# Create SummarizedExperiment object
se <- create_se_object(efile = fdat, sfile = sdat)
 
# Use the cutoff = 5 to stratify the samples and remove features that have 
# less than 90% unique values (this removes the MYH2 gene from the analysis)
km_ucut(obj = se, bfname = "test", cutoff = 5, min_uval = 90, wpdf = FALSE)

A graph with Kaplan-Meier curves is created for gene (there is only MYCN here).

Additionally, this run will create the following two output files in the current working directory (the names of output files will be created automatically by adding the run information to the base name ‘bfname’).

a) Tab-delimited text file with the results “test_KM_ucut_5.txt”

b) CSV file with low/high sample labels “test_KM_ucut_5_labels.csv”

The format of these two files is identical to the output files described for function ‘km_opt_pcut’.

‘km_val_cut’

This function creates Kaplan-Meier survival curves for each feature in a validation data file by using a file with previously determined stratification cutoffs (one cutoff per feature) and calculates the log-rank test p-value. Since for each feature only one previously determined cutoff is used, it is equivalent to performing one log-rank test per feature (no optimization is performed). The file with previously determined stratification cutoffs is a tab-delimited file with the table that contains one or more feature and a stratification cutoff for each feature (the table is generated by functions ‘km_opt_scut’, ‘km_opt_pcut’, ‘km_qcut’ or ‘km_ucut’). It must have first two columns named as ‘tracking_id’ and ‘CUTOFF’. The ‘tracking_id’ column contains gene (feature) names, the ‘CUTOFF’ column contains stratification cutoff for each gene (feature). An example of how to use ‘km_val_cut’ to validate stratification cutoffs determined by ‘km_qcut’:

# Read names of training (fdat1) and validation (fdat2) gene expression data
# files and survival data file (sdat).
fdat1 <- system.file("extdata", "expression_data_1.txt", package = "kmcut")
fdat2 <- system.file("extdata", "expression_data_2.txt", package = "kmcut")
sdat <- system.file("extdata", "survival_data.txt", package = "kmcut")

# Create SummarizedExperiment object with training data
se1 <- create_se_object(efile = fdat1, sfile = sdat)

# Step 1: Run 'km_qcut' on the training data in 'se1'
km_qcut(obj = se1, bfname = "training_data", quant = 50, min_uval = 40)

Step 1 will create three output files in the current working directory:

1. PDF file with plots “training_data_KM_quant_50.pdf”
   
2. Tab-delimited text file with the cutoffs “training_data_KM_quant_50.txt”
   
3. CSV file with low/high sample labels “training_data_KM_quant_50.csv”

The format of these three files is identical to the output files described for function ‘km_opt_pcut’. The only difference is that the PDF file does not contain observed vs. expected optimization plots.

# Create SummarizedExperiment object with test data
se2 <- create_se_object(efile = fdat2, sfile = sdat)
 
# Step 2: Validate the thresholds from "training_data_KM_quant_50.txt" on
# the test data in 'se2'.
km_val_cut(infile = "training_data_KM_quant_50.txt", obj = se2, 
             bfname = "test", wpdf = TRUE, min_uval = 40)

Step 2 will create three output files in the current working directory:

1. PDF file with plots “test_KM_val.pdf”
2. Tab-delimited text file with the validation results “test_KM_val.txt”
   
3. CSV file with low/high sample labels “test_KM_val_labels.csv”

The format of these three files is identical to the output files described for function ‘km_opt_pcut’. The only difference is that the PDF file does not contain observed vs. expected optimization plots.

‘ucox_batch’

This function fits a univariate Cox proportional hazard regression model and performs the likelihood ratio test for each feature in the dataset. An example of how to use ‘ucox_batch’ with the data files included in the package:

# Read names of the built-in gene expression data file (fdat) and
# survival data file (sdat)
fdat = system.file("extdata", "example_genes.txt", package = "kmcut")
sdat = system.file("extdata", "survival_data.txt", package = "kmcut")

# Create SummarizedExperiment object
se <- create_se_object(efile = fdat, sfile = sdat)

# Perform the regression on the data in 'se'
ucox_batch(obj = se, bfname = "test")

## Processing 1 of 2

## Processing 2 of 2

This will create in the current working directory a tab-delimited text file with results “test_ucoxbatch.txt”

Table 5. An illustration of the file format for the tab-delimited text file with Cox regression summary. 1st column – gene id, 2nd – concordance coefficient, 3rd – hazard ratio, 4th – likelihood ratio test p-value, 5th – FDR-adjusted p-value.

‘ucox_pred’

This function fits a univariate Cox proportional hazard regression model for each feature in the training dataset and then uses the models to calculate risk scores for the same features in the test dataset. An example of how to use ‘ucox_pred’ with the data files included in the package:

# Read names of the built-in training (fdat1) and test (fdat2) 
# gene expression data files and survival data file (sdat)
fdat1 = system.file("extdata", "expression_data_1.txt", package = "kmcut")
fdat2 = system.file("extdata", "expression_data_2.txt", package = "kmcut")
sdat = system.file("extdata", "survival_data.txt", package = "kmcut")

# Create SummarizedExperiment object with training data
se1 <- create_se_object(efile = fdat1, sfile = sdat)
# Create SummarizedExperiment object with test data
se2 <- create_se_object(efile = fdat2, sfile = sdat)

# Fit Cox model on the training data in 'se1' and use it to calculate the risk
# scores for the test data in 'se2'.
ucox_pred(obj1 = se1, obj2 = se2, bfname = "demo", min_uval = 40)

## Processing 1 of 2

## Processing 2 of 2

This will create three output files in the current working directory:

a) Tab-delimited text file with Cox regression summary for the training data
“demo_cox_train_sum.txt”

Table 6. An illustration of the file format for the tab-delimited text file with Cox regression summary. 1st column – gene id, 2nd – concordance coefficient, 3rd – hazard ratio, 4th – likelihood ratio test p-value, 5th – FDR-adjusted p-value.

b) Tab-delimited text file with the risk scores for the training data
“demo_train_score.txt”
In this file, rows are genes (features) and columns are samples.

c) Tab-delimited text file with the risk scores for the test data
“demo_test_score.txt”
In this file, rows are genes (features) and columns are samples.

‘extract_rows’

This function extracts a sub-set of rows (such as a group of genes) from a data table. All columns will be preserved. Names of the rows to be extracted must be in a text file, one name per line (the exact names, case-sensitive, no extra symbols are allowed). An example of how to use the function with the data files included in the package:

# Read the name of the built-in gene expression data file with 2 genes (2 rows)
fdat = system.file("extdata", "example_genes.txt", package = "kmcut")
# Read the name of the built-in list file that contains one gene id (MYCN)
idlist = system.file("extdata", "rowids.txt", package = "kmcut")

# Run the function
extract_rows(fnamein = fdat, fids = idlist,
            fnameout = "example_genes_subset.txt")

This will create a tab-delimited text file “example_genes_subset.txt” with one row “MYCN” in the current working directory.

‘extract_columns’

This function extracts a sub-set of columns (such as a group of samples) from a data table. All rows will be preserved. Names of the columns to be extracted must be in a text file, one name per line (the exact names, case-sensitive, no extra symbols are allowed). An example of how to use the function with the data files included in the package:

# Read the name of the built-in gene expression data file with 2 genes (2 rows)
fdat = system.file("extdata", "example_genes.txt", package = "kmcut")
# Read the name of the built-in list file that contains a sub-set of 
# column (sample) ids
idlist = system.file("extdata", "columnids.txt", package = "kmcut")

# Run the function
extract_columns(fnamein = fdat, fids = idlist,
                    fnameout = "example_samples_subset.txt")

This will create a tab-delimited text file “example_samples_subset.txt” in the current working directory. This file will contain columns (samples) from the list.

‘transpose_table’

This function transposes a data table (that is, converts rows to columns and columns to rows). Row names will become column names, and column names will become row names. An example of how to use the function with the data file included in the package:

# Read the name of the built-in gene expression data file.
# In this file, genes are rows and samples are columns.
fdat = system.file("extdata", "example_genes.txt", package = "kmcut")

# Run the function
transpose_table(fnamein = fdat, fnameout = "example_genes_transposed.txt")

This will create a tab-delimited text file “example_genes_295_transposed.txt” with the transposed table in the current working directory. In this file, genes (features) are columns and samples are rows.