Problems with survival analysis data

[ghost]

Hi everyone.

Hope you are all doing well.

I recently posted about my difficulty in locating the p-values generated by tcga_surv_plot() in the discussion thread here. Following a suggestion, I attempted to perform a coxph analysis with the tcga_surv_get() but encountered a new problem, I can't seem to replicate the p-values from the plots.

You can follow my steps and compare the results below:
Step 1: Retrieving Data and Generating Plot
I started by retrieving the clinical data and generating the survival plot:

#Variables used within both functions
gene_ID <- "ZNF714"
gene_var <- "mRNA"
tumor <- "UCEC"

# Getting clinical data 
data_surv <- tcga_surv_get(
  item = gene_ID,
  TCGA_cohort = tumor, 
  profile = gene_var,
  TCGA_cli_data = dplyr::full_join(load_data("tcga_clinical"), load_data("tcga_surv"),
                                   by = "sample"))

#Generating plot 
plot_surv <- tcga_surv_plot(
  data_surv,
  time = "OS.time",
  status = "OS",
  cutoff_mode = "Auto"
  cutpoint = c(50, 50),
  cnv_type = c("Normal"),
  profile = gene_var,
  palette = "aaas")

#Printing plot 
print(plot_surv)

Step 2: Cox Proportional Hazards Model
Next, I attempted to fit the Cox proportional hazards model and check assumptions:

#Prepare the data for Cox analysis
data_surv$OS.time <- as.numeric(data_surv$OS.time)
data_surv$OS <- as.numeric(data_surv$OS)
data_surv$value <- as.numeric(data_surv$value)

#Fit the Cox proportional hazards model
cox_model <- coxph(Surv(OS.time, OS) ~ value, data = data_surv)

# Check assumptions
cox_zph <- cox.zph(cox_model)
cox_zph # Print the results

Coxzph output:

Plot output:

Any suggestions or insights into what might be causing the discrepancies in p-values would be highly appreciated.

Furthermore, the plots don't show the log-rank statistics (even though they have "Log-rank" written in them). I'm unsure if this omission is intentional or if there's a way to include log-rank statistics in the plots, similar to what's displayed in the Xena Browser Kaplan-Meier plot:
E.g.:

[ShixiangWang]

@Banirelia You should use the following code to set up the 50-50 cutoff.

plot_surv <- tcga_surv_plot(
  data_surv,
  time = "OS.time",
  status = "OS",
  cutoff_mode = "Custom",
  cutpoint = c(50, 50),
  profile = gene_var,
  palette = "aaas")

For the result difference, it should be explained by the original data difference. The function used ?UCSCXenaShiny::query_pancan_value to query TPM data from Toil datahub for TCGA analysis. To note, the tcga_surv_plot included more data data_surv |> dplyr::filter(!is.na(value), !is.na(OS.time), !is.na(OS)) |> nrow() (188) than your plot using UCSC Xena browser (175).

In summary, the following reasons should be considered to explain the difference:

• Data hub (Toil)
• Normalized method (TPM)
• Other hidden preprocessing differences between the two tools

As the tcga_surv_get is open-source, you can debug the code debug(tcga_surv_get) to check if there are any problems in the preprocessing.

For the stats, I cannot find an option in the http://rpkgs.datanovia.com/survminer/reference/ggsurvplot.html to enable this, so I cannot support it currently. However, you can obtain the stats by yourself. An example is given.

library(survival)

> summary(coxph(Surv(time, status) ~ ph.ecog, data=lung))
Call:
coxph(formula = Surv(time, status) ~ ph.ecog, data = lung)

  n= 227, number of events= 164 
   (因为不存在，1个观察量被删除了)

          coef exp(coef) se(coef)     z  Pr(>|z|)    
ph.ecog 0.4759    1.6095   0.1134 4.198 0.0000269 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

        exp(coef) exp(-coef) lower .95 upper .95
ph.ecog      1.61     0.6213     1.289      2.01

Concordance= 0.604  (se = 0.024 )
Likelihood ratio test= 17.57  on 1 df,   p=0.00003
Wald test            = 17.62  on 1 df,   p=0.00003
Score (logrank) test = 17.89  on 1 df,   p=0.00002

[ghost]

Thank you for the feedback provided earlier. I now understand the limitations regarding the log-rank statistics and the differences between the Xena browser data.

However, I'm still facing challenges understanding why there's a significant difference between the p-values generated by the Cox proportional hazards model (coxph) and those produced by the tcga_surv_plot function, even when using the cutoff_mode = "Custom" argument and the same data from tcga_surv_get.

Below is the code snippet I've been using, along with the outputs demonstrating the discrepancies:

#Variables used within both functions
gene_ID <- "GLIS3"
gene_var <- "mRNA"
tumor <- "LGG"

# Getting clinical data 
data_surv <- tcga_surv_get(
  item = gene_ID,
  TCGA_cohort = tumor, 
  profile = gene_var,
  TCGA_cli_data = dplyr::full_join(load_data("tcga_clinical"), load_data("tcga_surv"),
                                   by = "sample"))

#Generating plot 
plot_surv <- tcga_surv_plot(
  data_surv,
  time = "OS.time",
  status = "OS",
  cutoff_mode = "Custom",
  cutpoint = c(50, 50),
  profile = gene_var,
  palette = "aaas")

#Printing plot 
print(plot_surv)

#Prepare the data for Cox analysis
data_surv$OS.time <- as.numeric(data_surv$OS.time)
data_surv$OS <- as.numeric(data_surv$OS)
data_surv$value <- as.numeric(data_surv$value)

#Fit the Cox proportional hazards model
cox_model <- coxph(Surv(OS.time, OS) ~ value, data = data_surv)

# Check assumptions
cox_zph <- cox.zph(cox_model)
cox_zph # Print the results

Outputs, how the coxph p-value and the tcga_surv_plot differ significantly, even using the same data from tcga_surv_get:

I'm aiming for at least close enough results between the Coxph model and the plots generated for an accurate batch analysis. Any insights or assistance on resolving this discrepancy would be greatly appreciated.

Thank you for your help!

[lishensuo]

The tcga_surv_plot() function is utilized to perform KM (Kaplan-Meier) survival analysis based on log-rank test and the model is different with Cox regression model. You can manually calculate the p.value of KM analysis after getting data_surv data.frame.

data_surv <- data_surv %>%
  dplyr::arrange(.data$value) %>%
  dplyr::mutate(per_rank = 100 / nrow(.) * (1:nrow(.)))
## dividing into two groups as the 'group' column
cutpoint = c(50, 50)
data_surv <- data_surv %>%
  dplyr::mutate(group = dplyr::case_when(
    .data$per_rank > !!cutpoint[2] ~ "High",
    .data$per_rank <= !!cutpoint[1] ~ "Low",
    TRUE ~ NA_character_
  ))

surv_diff <- survdiff(Surv(OS.time, OS) ~ group, data = data_surv)
p.val = 1 - pchisq(surv_diff$chisq, length(surv_diff$n) - 1)
p.val
# [1] 0.001764663

[ghost]

This is perfect, thank you very much!!