start_time <- Sys.time()

References

Cox Proportional-Hazards Model reference used was (Kassambara 2023):

  • (Kassambara 2023);

Comparison between KM and Cox:

  • (Marin and Hamadani 2020a);

Complete Survival Implementation in R:

  • (Marin and Hamadani 2020b).

Load Libraries

library(survival)
library(survminer)
library(dplyr)
library(tibble)
library(kableExtra)
library(reshape)
library(expss)

Cox Regression vs. Kaplan Meier analysis

Survival analysis as evaluated in the Kaplan-Meier (KM) method, allows for a simple interpretation of the results and evaluation of a parameter such as a gene set score effect in the survival rates or survival function. However, KM does not allow the calculation or estimation of risk or hazards associated with a phenotype. For this reason we are required to use another method to calculate Hazards Ratio (HR) of a phenotype.

Kaplan Meier Cox Regression
Pros Simple interpretation Hazards fluctuate with time
Can estimate survival function S(t) Can estimate H.R.
Handles categorical and continuous variables
Cons No functional form
Can not handle H.R. Can not estimate S(t)
Only includes few categorical variables

In KM analyis we evaluated association of phenotype with survival rates. In survival analyis, the method we use to estimate HR of a phenotype is the Cox Regression model. The Cox Regression model can convey univariate or multivariate analysis results.

Biological pathways that interfere with survival outcomes

  • Immunotherapies are mediated by monoclonal antibodies or other entities that interact with the host immune system and are responsible for “enhancing” the immune response;

  • One enhancement mechanism occurs through the Antibody Dependent Cellular Cytotoxicity (ADCC) response: host effector cells (immune cells) recognize, bind and destroy tumor cells “decorated” with antibodies.

Figure 4.1. Antibody Dependent Cellular Cytotoxicity (ADCC) mediates immune attacks of the host innate immune system against tumor cells. Figure is based on Lo Nigro et al., 2019.

Biological pathways that interfere with survival outcomes (continued)

  • There must be an antigen recognized by the antibody;

  • The antibody then is recognized in its Fc fragment, by the Fc receptor of effector cells such as macrophages, NK cells, B lymphocytes, follicular dendritic cells, neutrophils, eosinophils, basophils, human platelets and mast cells.

Figure 4.2. Fc receptor in effector cells, recognize Fc fragment in antibody bound to target cell (e.g.: pathogen or tumor cell).

ADRN to MES transition (AMT) increases resistance to immunotherapy

  • ADRN to MES transition (AMT) drives neuroblastoma aggressiveness;

  • AMT is thought to drive aggressiveness by means of chemotherapy and immunotherapy resistance.

  • The figure bellow, depicts the AMT involvement in immunotherapy resistance.

Figure 4.3. ADRN to MES transition described in Mabe et al., 2022.

ADRN to MES transition (AMT) increases resistance to immunotherapy

  • Given that hypoxia can be one of the “natural” stimuli that causes AMT, how are HIF1A and hypoxia involved in chemotherapy resistance mechanisms?

  • Chemotherapy resistance mechanisms might involve quiescence and exit or arrest of the cell cycle.

  • HIF1A is involved in both immunotherapy and chemotherapy resistance.

1) Survival Analysis: Cox Regression

1.1) Multivariate Cox Regression analysis

In the code bellow, we will evaluate how HIF1A and other hypoxia gene markers, correlate with the survival information stored in the R2 GSE62564 RNA-Seq data-frame.

3.1) Using the survivalAnalysis R package

  • The first reference we try uses the survivalAnalysis package in R

  • This is the reference we followed:

https://cran.r-project.org/web/packages/survivalAnalysis/vignettes/multivariate.html

3.2) Using the Ggforest R package (a relative to the ggplot package?)

https://rpkgs.datanovia.com/survminer/reference/ggforest.html

The link above, shows how to use the GGforest R package (Alboukadel Kassambara 2023).

We begin by showing that we can repeat the code used by others (Alboukadel Kassambara 2023):

3.3) Kocak data

1.2) Make variables of interest, categorical

r2_gse62564_GSVA_Metadata <- readRDS("../results/r2_gse62564_GSVA_Metadata.rds")

Gene_GSVA_Metadata_DF        <- r2_gse62564_GSVA_Metadata

library(maditr)
## OS Status must be made numeric; check the neuroblastoma-r2-plots-GSE62564.html notebook
Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
  dplyr::mutate(os_bin_numeric = ifelse(Gene_GSVA_Metadata_DF$os_bin == "event", "1", "0"))
                
## Change the os_bin_numeric variable to numeric
Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
  dplyr::mutate(os_bin_numeric = as.numeric(os_bin_numeric))

## Change the OS time variable to numeric
Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
  dplyr::rename(os_days = os_day)%>%
  dplyr::mutate(os_days = as.numeric(os_days))
  
## EFS Status must be numeric
Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
  dplyr::mutate(efs_bin_numeric =
                  ifelse(Gene_GSVA_Metadata_DF$efs_bin == "event", "1", "0"))

## Change the OS time variable to numeric
Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
  dplyr::rename(efs_days = efs_day)%>%
  dplyr::mutate(efs_days = as.numeric(efs_days))

Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
  ## Mutate current risk category to high_risk
  dplyr::rename(high_risk = high_risk)

cross_cases(Gene_GSVA_Metadata_DF, high_risk)
 #Total 
 high_risk 
   no  322
   yes  176
   #Total cases  498 
Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
  dplyr::rename(mycn_status = mycn_status)

cross_cases(Gene_GSVA_Metadata_DF, high_risk)
 #Total 
 high_risk 
   no  322
   yes  176
   #Total cases  498 
## Construct the high-risk bin with 1/0
## Strangely, I was not able to do this in the previous block
Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
  dplyr::mutate(high_risk = ifelse(Gene_GSVA_Metadata_DF$high_risk == "yes", "1", "0"))

cross_cases(Gene_GSVA_Metadata_DF, high_risk)
 #Total 
 high_risk 
   0  322
   1  176
   #Total cases  498 
## Here, os_day and efs_day refer to how the variables (columns) are named in the Kocak (n=498) dataset.
Gene_GSVA_Metadata_DF["os_days"]                        <- as.numeric(Gene_GSVA_Metadata_DF$os_days)  ## os_days variable is defined above, in line 143
Gene_GSVA_Metadata_DF["efs_days"]                       <- as.numeric(Gene_GSVA_Metadata_DF$efs_days) ## Not sure if os_days and efs_days need to be here
Gene_GSVA_Metadata_DF["os_bin_numeric"]                 <- as.numeric(Gene_GSVA_Metadata_DF$os_bin_numeric)
Gene_GSVA_Metadata_DF$efs_bin_numeric                   <- as.numeric(Gene_GSVA_Metadata_DF$efs_bin_numeric)
Gene_GSVA_Metadata_DF["HIF1A"]                          <- as.numeric(Gene_GSVA_Metadata_DF$HIF1A)
Gene_GSVA_Metadata_DF["T cell exhaustion"]              <- as.numeric(Gene_GSVA_Metadata_DF$`T cell exhaustion`)
Gene_GSVA_Metadata_DF["ELVIDGE_HIF1A_TARGETS_UP"]       <- as.numeric(Gene_GSVA_Metadata_DF$ELVIDGE_HIF1A_TARGETS_UP)
Gene_GSVA_Metadata_DF["HALLMARK_INFLAMMATORY_RESPONSE"] <- as.numeric(Gene_GSVA_Metadata_DF$HALLMARK_INFLAMMATORY_RESPONSE)
Gene_GSVA_Metadata_DF["Quiescence Up Cheung"]           <- as.numeric(Gene_GSVA_Metadata_DF$`Quiescence Up Cheung`)
Gene_GSVA_Metadata_DF["Quiescence Down Cheung"]         <- as.numeric(Gene_GSVA_Metadata_DF$`Quiescence Down Cheung`)
Gene_GSVA_Metadata_DF["ST8SIA1"]                        <- as.numeric(Gene_GSVA_Metadata_DF$ST8SIA1)

library(dplyr)
Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
  mutate(hif1a_qtle       = ntile(HIF1A, 3))%>%
  mutate(hif1a_qtle       = as.factor(hif1a_qtle))%>%
  mutate(hallm_hypo_qtle  = ntile(HALLMARK_HYPOXIA, 3))%>%
  mutate(hallm_hypo_qtle  = as.factor(hallm_hypo_qtle))%>%
  mutate(t_cell_qtle      = ntile(`T cell exhaustion`, 3))%>%
  mutate(t_cell_qtle      = as.factor(t_cell_qtle))%>%
  mutate(elvige_Up_qtle   = ntile(ELVIDGE_HIF1A_TARGETS_UP, 3))%>%
  mutate(elvige_Up_qtle   = as.factor(elvige_Up_qtle))%>%
  mutate(hall_inflm_qtle  = ntile(HALLMARK_INFLAMMATORY_RESPONSE, 3))%>%
  mutate(hall_inflm_qtle  = as.factor(hall_inflm_qtle))%>%
  mutate(qui_up_cheung    = ntile(`Quiescence Up Cheung`, 3))%>%
  mutate(qui_up_cheung    = as.factor(qui_up_cheung))%>%
  mutate(qui_down_cheung  = ntile(`Quiescence Down Cheung`, 3))%>%
  mutate(qui_down_cheung  = as.factor(qui_down_cheung))%>%
  mutate(st8sia1_qtle     = ntile(ST8SIA1, 3))%>%
  mutate(st8sia1_qtle     = as.factor(st8sia1_qtle))

1.3) Survival analysis of the hypoxia gene sets

1.3.1) OS

In this part, it is best to code the gene sets of interest first considering the OS analysis to then evaluate the EFS analysis.

The advantage here is that we can use the code and only modify the gene set of interest.

Also, we can follow the procedure from above, where the ntile() function was used to construct 3 quantiles, or construct 2 quantiles as required.

1.3.1.1) HIF1A

We begin visualizing the quatiles that we constructed, in this case, with the values of the expression of HIF1A.

library(tidyverse)

Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
    mutate(quantile = ntile(HIF1A, 2))

Gene_GSVA_Metadata_DF$quantile <- as.factor(Gene_GSVA_Metadata_DF$quantile )

qplot(HALLMARK_HYPOXIA, HIF1A, 
      data = Gene_GSVA_Metadata_DF,
      colour=quantile,
      ylab = "HIF1A",
      xlab = "Hallmark Hypoxia")

library(survival)
fit_os_HIF1A <- survfit(Surv(os_days, os_bin_numeric) ~ Gene_GSVA_Metadata_DF$quantile,
                   data = Gene_GSVA_Metadata_DF)

library(survminer)
ggsurvplot(fit_os_HIF1A, 
           data = Gene_GSVA_Metadata_DF, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

1.3.1.2) HALLMARK_HYPOXIA

library(tidyverse)

Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
    mutate(quantile = ntile(HALLMARK_HYPOXIA, 2))

Gene_GSVA_Metadata_DF$quantile <- as.factor(Gene_GSVA_Metadata_DF$quantile )

qplot(HALLMARK_HYPOXIA, HALLMARK_HYPOXIA, 
      data = Gene_GSVA_Metadata_DF,
      colour=quantile,
      ylab = "HALLMARK_HYPOXIA",
      xlab = "Hallmark Hypoxia")

library(survival)
fit_os_geneSetInterest <- survfit(Surv(os_days, os_bin_numeric) ~ Gene_GSVA_Metadata_DF$quantile,
                   data = Gene_GSVA_Metadata_DF)

library(survminer)
ggsurvplot(fit_os_geneSetInterest, 
           data = Gene_GSVA_Metadata_DF, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

1.3.1.3) HALLMARK_INFLAMMATORY_RESPONSE

library(tidyverse)

Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
    mutate(quantile = ntile(HALLMARK_INFLAMMATORY_RESPONSE, 2))

Gene_GSVA_Metadata_DF$quantile <- as.factor(Gene_GSVA_Metadata_DF$quantile )

qplot(HALLMARK_HYPOXIA, HALLMARK_INFLAMMATORY_RESPONSE, 
      data = Gene_GSVA_Metadata_DF,
      colour=quantile,
      ylab = "HALLMARK_INFLAMMATORY_RESPONSE",
      xlab = "Hallmark Hypoxia")

library(survival)
fit_os_geneSetInterest <- survfit(Surv(os_days, os_bin_numeric) ~ Gene_GSVA_Metadata_DF$quantile,
                   data = Gene_GSVA_Metadata_DF)

library(survminer)
ggsurvplot(fit_os_geneSetInterest, 
           data = Gene_GSVA_Metadata_DF, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

1.3.1.4) T cell exhaustion

In this part it is necessary to solve the spaces problem in the names of the variables.

The way the Gene_GSVA_Metadata_DF dataframe was constructed does not allow the “T cell exhaustion” dataframe to be processed.

We need to use the strategy used in notebook “06-neuroblastoma-target-plots.Rmd” when we used the function str_replace_all from package stringr.

In this strategy, we replace spaces with “_“. The difference here is that we’re not changing the dataframe name.

library(stringr)
names(Gene_GSVA_Metadata_DF)<-str_replace_all(names(Gene_GSVA_Metadata_DF), c(" " = "_"))

Then we can use the T_cell_exhaustion as the variable name:

library(tidyverse)

Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
    mutate(quantile = ntile(T_cell_exhaustion, 2))

Gene_GSVA_Metadata_DF$quantile <- as.factor(Gene_GSVA_Metadata_DF$quantile )

qplot(HALLMARK_HYPOXIA, T_cell_exhaustion, 
      data = Gene_GSVA_Metadata_DF,
      colour=quantile,
      ylab = "T cell exhaustion",
      xlab = "Hallmark Hypoxia")

library(survival)
fit_os_geneSetInterest <- survfit(Surv(os_days, os_bin_numeric) ~ Gene_GSVA_Metadata_DF$quantile,
                   data = Gene_GSVA_Metadata_DF)

library(survminer)
ggsurvplot(fit_os_geneSetInterest, 
           data = Gene_GSVA_Metadata_DF, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

1.3.1.5) ELVIDGE_HIF1A_TARGETS_UP

library(tidyverse)

Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
    mutate(quantile = ntile(ELVIDGE_HIF1A_TARGETS_UP, 2))

Gene_GSVA_Metadata_DF$quantile <- as.factor(Gene_GSVA_Metadata_DF$quantile )

qplot(HALLMARK_HYPOXIA, ELVIDGE_HIF1A_TARGETS_UP, 
      data = Gene_GSVA_Metadata_DF,
      colour=quantile,
      ylab = "ELVIDGE_HIF1A_TARGETS_UP",
      xlab = "Hallmark Hypoxia")

library(survival)
fit_os_geneSetInterest <- survfit(Surv(os_days, os_bin_numeric) ~ Gene_GSVA_Metadata_DF$quantile,
                   data = Gene_GSVA_Metadata_DF)

library(survminer)
ggsurvplot(fit_os_geneSetInterest, 
           data = Gene_GSVA_Metadata_DF, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

1.3.1.6) Quiescence Up Cheung

library(tidyverse)

Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
    mutate(quantile = ntile(Quiescence_Up_Cheung, 2))

Gene_GSVA_Metadata_DF$quantile <- as.factor(Gene_GSVA_Metadata_DF$quantile )

qplot(HALLMARK_HYPOXIA, Quiescence_Up_Cheung, 
      data = Gene_GSVA_Metadata_DF,
      colour=quantile,
      ylab = "Quiescence Up Cheung",
      xlab = "Hallmark Hypoxia")

library(survival)
fit_os_geneSetInterest <- survfit(Surv(os_days, os_bin_numeric) ~ Gene_GSVA_Metadata_DF$quantile,
                   data = Gene_GSVA_Metadata_DF)

library(survminer)
ggsurvplot(fit_os_geneSetInterest, 
           data = Gene_GSVA_Metadata_DF, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

1.3.1.7) Quiescence Down Cheung

library(tidyverse)

Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
    mutate(quantile = ntile(Quiescence_Down_Cheung, 2))

Gene_GSVA_Metadata_DF$quantile <- as.factor(Gene_GSVA_Metadata_DF$quantile )

qplot(HALLMARK_HYPOXIA, Quiescence_Down_Cheung, 
      data = Gene_GSVA_Metadata_DF,
      colour=quantile,
      ylab = "Quiescence Down Cheung",
      xlab = "Hallmark Hypoxia")

library(survival)
fit_os_geneSetInterest <- survfit(Surv(os_days, os_bin_numeric) ~ Gene_GSVA_Metadata_DF$quantile,
                   data = Gene_GSVA_Metadata_DF)

library(survminer)
ggsurvplot(fit_os_geneSetInterest, 
           data = Gene_GSVA_Metadata_DF, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

1.3.1.8) ST8SIA1/GD2 surface expression

library(tidyverse)

Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
    mutate(quantile = ntile(ST8SIA1, 2))

Gene_GSVA_Metadata_DF$quantile <- as.factor(Gene_GSVA_Metadata_DF$quantile )

qplot(HALLMARK_HYPOXIA, ST8SIA1, 
      data = Gene_GSVA_Metadata_DF,
      colour=quantile,
      ylab = "ST8SIA1",
      xlab = "Hallmark Hypoxia")

library(survival)
fit_os_geneSetInterest <- survfit(Surv(os_days, os_bin_numeric) ~ Gene_GSVA_Metadata_DF$quantile,
                   data = Gene_GSVA_Metadata_DF)

library(survminer)
ggsurvplot(fit_os_geneSetInterest, 
           data = Gene_GSVA_Metadata_DF, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

1.3.2) EFS

Now we start the EFS, whose analysis code is similar to OS.

1.3.2.1) HIF1A

library(tidyverse)

Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
    mutate(quantile = ntile(HIF1A, 2))

Gene_GSVA_Metadata_DF$quantile <- as.factor(Gene_GSVA_Metadata_DF$quantile )

qplot(HALLMARK_HYPOXIA, HIF1A, 
      data = Gene_GSVA_Metadata_DF,
      colour=quantile,
      ylab = "HIF1A",
      xlab = "Hallmark Hypoxia")

library(survival)
fit_efs_HIF1A <- survfit(Surv(efs_days, efs_bin_numeric) ~ Gene_GSVA_Metadata_DF$quantile,
                   data = Gene_GSVA_Metadata_DF)

library(survminer)
ggsurvplot(fit_efs_HIF1A, 
           data = Gene_GSVA_Metadata_DF, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

1.3.2.2) HALLMARK_HYPOXIA

library(tidyverse)

Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
    mutate(quantile = ntile(HALLMARK_HYPOXIA, 2))

Gene_GSVA_Metadata_DF$quantile <- as.factor(Gene_GSVA_Metadata_DF$quantile )

qplot(HALLMARK_HYPOXIA, HALLMARK_HYPOXIA, 
      data = Gene_GSVA_Metadata_DF,
      colour=quantile,
      ylab = "HALLMARK_HYPOXIA",
      xlab = "Hallmark Hypoxia")

library(survival)
fit_efs_geneSetInterest <- survfit(Surv(efs_days, efs_bin_numeric) ~ Gene_GSVA_Metadata_DF$quantile,
                   data = Gene_GSVA_Metadata_DF)

library(survminer)
ggsurvplot(fit_efs_geneSetInterest, 
           data = Gene_GSVA_Metadata_DF, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

1.3.2.3) HALLMARK_INFLAMMATORY_RESPONSE

library(tidyverse)

Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
    mutate(quantile = ntile(HALLMARK_INFLAMMATORY_RESPONSE, 2))

Gene_GSVA_Metadata_DF$quantile <- as.factor(Gene_GSVA_Metadata_DF$quantile )

qplot(HALLMARK_HYPOXIA, HALLMARK_INFLAMMATORY_RESPONSE, 
      data = Gene_GSVA_Metadata_DF,
      colour=quantile,
      ylab = "HALLMARK_INFLAMMATORY_RESPONSE",
      xlab = "Hallmark Hypoxia")

library(survival)
fit_efs_geneSetInterest <- survfit(Surv(efs_days, efs_bin_numeric) ~ Gene_GSVA_Metadata_DF$quantile,
                   data = Gene_GSVA_Metadata_DF)

library(survminer)
ggsurvplot(fit_efs_geneSetInterest, 
           data = Gene_GSVA_Metadata_DF, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

1.3.2.4) T cell exhaustion

library(tidyverse)

Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
    mutate(quantile = ntile(T_cell_exhaustion, 2))

Gene_GSVA_Metadata_DF$quantile <- as.factor(Gene_GSVA_Metadata_DF$quantile )

qplot(HALLMARK_HYPOXIA, T_cell_exhaustion, 
      data = Gene_GSVA_Metadata_DF,
      colour=quantile,
      ylab = "T cell exhaustion",
      xlab = "Hallmark Hypoxia")

library(survival)
fit_efs_geneSetInterest <- survfit(Surv(efs_days, efs_bin_numeric) ~ Gene_GSVA_Metadata_DF$quantile,
                   data = Gene_GSVA_Metadata_DF)

library(survminer)
ggsurvplot(fit_efs_geneSetInterest, 
           data = Gene_GSVA_Metadata_DF, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

1.3.2.5) ELVIDGE_HIF1A_TARGETS_UP

library(tidyverse)

Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
    mutate(quantile = ntile(ELVIDGE_HIF1A_TARGETS_UP, 2))

Gene_GSVA_Metadata_DF$quantile <- as.factor(Gene_GSVA_Metadata_DF$quantile )

qplot(HALLMARK_HYPOXIA, ELVIDGE_HIF1A_TARGETS_UP, 
      data = Gene_GSVA_Metadata_DF,
      colour=quantile,
      ylab = "ELVIDGE_HIF1A_TARGETS_UP",
      xlab = "Hallmark Hypoxia")

library(survival)
fit_os_geneSetInterest <- survfit(Surv(efs_days, efs_bin_numeric) ~ Gene_GSVA_Metadata_DF$quantile,
                   data = Gene_GSVA_Metadata_DF)

library(survminer)
ggsurvplot(fit_efs_geneSetInterest, 
           data = Gene_GSVA_Metadata_DF, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

1.3.2.6) Quiescence Up Cheung

library(tidyverse)

Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
    mutate(quantile = ntile(Quiescence_Up_Cheung, 2))

Gene_GSVA_Metadata_DF$quantile <- as.factor(Gene_GSVA_Metadata_DF$quantile )

qplot(HALLMARK_HYPOXIA, Quiescence_Up_Cheung, 
      data = Gene_GSVA_Metadata_DF,
      colour=quantile,
      ylab = "Quiescence Up Cheung",
      xlab = "Hallmark Hypoxia")

library(survival)
fit_os_geneSetInterest <- survfit(Surv(efs_days, efs_bin_numeric) ~ Gene_GSVA_Metadata_DF$quantile,
                   data = Gene_GSVA_Metadata_DF)

library(survminer)
ggsurvplot(fit_efs_geneSetInterest, 
           data = Gene_GSVA_Metadata_DF, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

1.3.2.7) Quiescence Down Cheung

library(tidyverse)

Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
    mutate(quantile = ntile(Quiescence_Down_Cheung, 2))

Gene_GSVA_Metadata_DF$quantile <- as.factor(Gene_GSVA_Metadata_DF$quantile )

qplot(HALLMARK_HYPOXIA, Quiescence_Down_Cheung, 
      data = Gene_GSVA_Metadata_DF,
      colour=quantile,
      ylab = "Quiescence Down Cheung",
      xlab = "Hallmark Hypoxia")

library(survival)
fit_os_geneSetInterest <- survfit(Surv(efs_days, efs_bin_numeric) ~ Gene_GSVA_Metadata_DF$quantile,
                   data = Gene_GSVA_Metadata_DF)

library(survminer)
ggsurvplot(fit_efs_geneSetInterest, 
           data = Gene_GSVA_Metadata_DF, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

1.3.2.8) ST8SIA1/GD2 surface expression

library(tidyverse)

Gene_GSVA_Metadata_DF <- Gene_GSVA_Metadata_DF %>%
    mutate(quantile = ntile(ST8SIA1, 2))

Gene_GSVA_Metadata_DF$quantile <- as.factor(Gene_GSVA_Metadata_DF$quantile )

qplot(HALLMARK_HYPOXIA, ST8SIA1, 
      data = Gene_GSVA_Metadata_DF,
      colour=quantile,
      ylab = "ST8SIA1",
      xlab = "Hallmark Hypoxia")

library(survival)
fit_os_geneSetInterest <- survfit(Surv(efs_days, efs_bin_numeric) ~ Gene_GSVA_Metadata_DF$quantile,
                   data = Gene_GSVA_Metadata_DF)

library(survminer)
ggsurvplot(fit_efs_geneSetInterest, 
           data = Gene_GSVA_Metadata_DF, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

In this part we may encounter this error for the Kocak dataset:

Error in contrasts<-(*tmp*, value = contr.funs[1 + isOF[nn]]) : contrasts can be applied only to factors with 2 or more levels

I found this explanation:

This error occurs when you attempt to fit a regression model using a predictor variable that is either a factor or character and only has one unique value.

This was the reference used:

https://www.statology.org/contrasts-applied-to-factors-with-2-or-more-levels/

This happened because high_risk variable in the Kocak dataset, is set to be yes or no, not 1 or 2, as may have happened with the other datasets.

Check the levels in the high risk variable, to make sure it came with the expected number of levels.

cross_cases(Gene_GSVA_Metadata_DF, high_risk)
 #Total 
 high_risk 
   0  322
   1  176
   #Total cases  498

2) Forest Plots

2.1.1) OS

model_kocak <- coxph( Surv(os_days, os_bin_numeric) ~ high_risk + 
                        hif1a_qtle +
                        hallm_hypo_qtle +
                        hall_inflm_qtle +
                        t_cell_qtle +
                        elvige_Up_qtle +
                        qui_up_cheung +
                        qui_down_cheung +
                        st8sia1_qtle,
                      data = Gene_GSVA_Metadata_DF )
## Model Plot
ggforest(model_kocak)
## Warning in .get_data(model, data = data): The `data` argument is not provided.
## Data will be extracted from model fit.

model_kocak <- coxph( Surv(os_days, os_bin_numeric) ~ hif1a_qtle,
                data = Gene_GSVA_Metadata_DF )
## Model Plot
ggforest(model_kocak)
## Warning in .get_data(model, data = data): The `data` argument is not provided.
## Data will be extracted from model fit.

2.1.2) EFS

model_kocak <- coxph( Surv(efs_days, efs_bin_numeric) ~ high_risk + 
                        mycn_status +
                        inss_stage + # Something is not allowing plot of inss_stage
                        hif1a_qtle +
                        hallm_hypo_qtle +
                        hall_inflm_qtle +
                        t_cell_qtle +
                        elvige_Up_qtle+
                        qui_up_cheung +
                        qui_down_cheung +
                        st8sia1_qtle,
                      data = Gene_GSVA_Metadata_DF )
## Model Plot
ggforest(model_kocak)
## Warning in .get_data(model, data = data): The `data` argument is not provided.
## Data will be extracted from model fit.

model_kocak <- coxph( Surv(efs_days, efs_bin_numeric) ~ hif1a_qtle,
                data = Gene_GSVA_Metadata_DF )
## Model Plot
ggforest(model_kocak)
## Warning in .get_data(model, data = data): The `data` argument is not provided.
## Data will be extracted from model fit.

3) Cumulative incidence Plots

  • Reference

https://www.emilyzabor.com/tutorials/survival_analysis_in_r_tutorial.html

Kocak

3.1) INSS stage

OS

library(tidycmprsk)
## Warning: package 'tidycmprsk' was built under R version 4.1.2
Gene_GSVA_Metadata_DF$os_bin_numeric <- as.factor(Gene_GSVA_Metadata_DF$os_bin_numeric)

Gene_GSVA_Metadata_DF$os_day <- as.numeric(Gene_GSVA_Metadata_DF$os_day)

library(ggsurvfit)
## Warning: package 'ggsurvfit' was built under R version 4.1.2
cuminc(Surv(os_day, os_bin_numeric) ~ mycn_status, data = Gene_GSVA_Metadata_DF) %>% 
  ggcuminc(outcome = c()) + # This line did not work
  # ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

Gene_GSVA_Metadata_DF$high_risk <- as.factor(Gene_GSVA_Metadata_DF$high_risk)

cross_cases(Gene_GSVA_Metadata_DF, high_risk)
 #Total 
 high_risk 
   0  322
   1  176
   #Total cases  498 
cuminc(Surv(os_day, os_bin_numeric) ~ high_risk, data = Gene_GSVA_Metadata_DF) %>% 
  ggcuminc(outcome = c()) + # This line did not work
  # ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

3.2) hif1a_qtle

3.2.1) OS

library(tidycmprsk)
# r2_gse62564_GSVA_Metadata <- readRDS("../results/r2_gse62564_GSVA_Metadata.rds")

Gene_GSVA_Metadata_DF$os_bin_numeric <- as.factor(Gene_GSVA_Metadata_DF$os_bin_numeric)

Gene_GSVA_Metadata_DF$os_day <- as.numeric(Gene_GSVA_Metadata_DF$os_day)

## HALLMARK_HYPOXIA has not been defined as numeric for the Gene_GSVA_Metadata_DF dataframe
Gene_GSVA_Metadata_DF$HALLMARK_HYPOXIA <- as.numeric(Gene_GSVA_Metadata_DF$HALLMARK_HYPOXIA)

qplot(data = Gene_GSVA_Metadata_DF, colour = Gene_GSVA_Metadata_DF$hif1a_qtle,
        HALLMARK_HYPOXIA,
        HIF1A,
        xlab = "HALLMARK_HYPOXIA",
        ylab = "hif1a_qtle",
      main = "Hallmark Hypoxia vs. HIF1A Quantiles")+
       labs(colour = "hif1a_qtle")+
  theme_bw()

library(ggsurvfit)
cuminc(Surv(os_day, os_bin_numeric) ~ hif1a_qtle, data = Gene_GSVA_Metadata_DF) %>% 
  # ggcuminc(outcome = c(0, 1)) + # This line did not work 
  ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

3.2.2) EFS

library(tidycmprsk)
# r2_gse62564_GSVA_Metadata <- readRDS("../results/r2_gse62564_GSVA_Metadata.rds")

Gene_GSVA_Metadata_DF$efs_bin_numeric <- as.factor(Gene_GSVA_Metadata_DF$efs_bin_numeric)

Gene_GSVA_Metadata_DF$efs_days <- as.numeric(Gene_GSVA_Metadata_DF$efs_days)

library(ggsurvfit)
cuminc(Surv(efs_days, efs_bin_numeric) ~ hif1a_qtle, data = Gene_GSVA_Metadata_DF) %>% 
  # ggcuminc(outcome = c(0, 1)) + # This line did not work 
  ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

3.3) hallm_hypo_qtle

3.3.1) OS

library(tidycmprsk)
# r2_gse62564_GSVA_Metadata <- readRDS("../results/r2_gse62564_GSVA_Metadata.rds")

Gene_GSVA_Metadata_DF$os_bin_numeric <- as.factor(Gene_GSVA_Metadata_DF$os_bin_numeric)

Gene_GSVA_Metadata_DF$os_day <- as.numeric(Gene_GSVA_Metadata_DF$os_day)

## HALLMARK_HYPOXIA has not been defined as numeric for the Gene_GSVA_Metadata_DF dataframe
Gene_GSVA_Metadata_DF$HALLMARK_HYPOXIA <- as.numeric(Gene_GSVA_Metadata_DF$HALLMARK_HYPOXIA)

qplot(data = Gene_GSVA_Metadata_DF, colour = Gene_GSVA_Metadata_DF$hallm_hypo_qtle,
        HALLMARK_HYPOXIA,
        HALLMARK_HYPOXIA,
        xlab = "HALLMARK_HYPOXIA",
        ylab = "HALLMARK_HYPOXIA",
      main = "Hallmark Hypoxia vs. HIF1A Quantiles")+
       labs(colour = "hallm_hypo_qtle")+
  theme_bw()

library(ggsurvfit)
cuminc(Surv(os_day, os_bin_numeric) ~ hallm_hypo_qtle, data = Gene_GSVA_Metadata_DF) %>% 
  # ggcuminc(outcome = c(0, 1)) + # This line did not work 
  ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

3.3.2) EFS

library(tidycmprsk)
# r2_gse62564_GSVA_Metadata <- readRDS("../results/r2_gse62564_GSVA_Metadata.rds")

Gene_GSVA_Metadata_DF$efs_bin_numeric <- as.factor(Gene_GSVA_Metadata_DF$efs_bin_numeric)

Gene_GSVA_Metadata_DF$efs_days <- as.numeric(Gene_GSVA_Metadata_DF$efs_days)

library(ggsurvfit)
cuminc(Surv(efs_days, efs_bin_numeric) ~ hallm_hypo_qtle, data = Gene_GSVA_Metadata_DF) %>% 
  # ggcuminc(outcome = c(0, 1)) + # This line did not work 
  ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

3.4) HALLMARK_INFLAMMATORY_RESPONSE

3.4.1) OS

library(tidycmprsk)
# r2_gse62564_GSVA_Metadata <- readRDS("../results/r2_gse62564_GSVA_Metadata.rds")

Gene_GSVA_Metadata_DF$os_bin_numeric <- as.factor(Gene_GSVA_Metadata_DF$os_bin_numeric)

Gene_GSVA_Metadata_DF$os_days <- as.numeric(Gene_GSVA_Metadata_DF$os_days)

## HALLMARK_HYPOXIA has not been defined as numeric for the Gene_GSVA_Metadata_DF dataframe
Gene_GSVA_Metadata_DF$HALLMARK_HYPOXIA <- as.numeric(Gene_GSVA_Metadata_DF$HALLMARK_HYPOXIA)

qplot(data = Gene_GSVA_Metadata_DF, colour = Gene_GSVA_Metadata_DF$hall_inflm_qtle,
        HALLMARK_HYPOXIA,
        HALLMARK_INFLAMMATORY_RESPONSE,
        xlab = "HALLMARK_HYPOXIA",
        ylab = "hall_inflm_qtle",
      main = "Hallmark Hypoxia vs. HIF1A Quantiles")+
       labs(colour = "hall_inflm_qtle")+
  theme_bw()

library(ggsurvfit)
cuminc(Surv(os_days, os_bin_numeric) ~ hall_inflm_qtle, data = Gene_GSVA_Metadata_DF) %>% 
  # ggcuminc(outcome = c(0, 1)) + # This line did not work 
  ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

3.4.2) EFS

library(tidycmprsk)
# r2_gse62564_GSVA_Metadata <- readRDS("../results/r2_gse62564_GSVA_Metadata.rds")

Gene_GSVA_Metadata_DF$efs_bin_numeric <- as.factor(Gene_GSVA_Metadata_DF$efs_bin_numeric)

Gene_GSVA_Metadata_DF$efs_days <- as.numeric(Gene_GSVA_Metadata_DF$efs_days)

library(ggsurvfit)
cuminc(Surv(efs_days, efs_bin_numeric) ~ hall_inflm_qtle, data = Gene_GSVA_Metadata_DF) %>% 
  # ggcuminc(outcome = c(0, 1)) + # This line did not work 
  ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

3.5) T-cell Exhaustion

3.5.1) OS

library(tidycmprsk)
# r2_gse62564_GSVA_Metadata <- readRDS("../results/r2_gse62564_GSVA_Metadata.rds")

Gene_GSVA_Metadata_DF$os_bin_numeric <- as.factor(Gene_GSVA_Metadata_DF$os_bin_numeric)

Gene_GSVA_Metadata_DF$os_days <- as.numeric(Gene_GSVA_Metadata_DF$os_days)

## HALLMARK_HYPOXIA has not been defined as numeric for the Gene_GSVA_Metadata_DF dataframe
Gene_GSVA_Metadata_DF$T_cell_exhaustion <- as.numeric(Gene_GSVA_Metadata_DF$T_cell_exhaustion)

qplot(data = Gene_GSVA_Metadata_DF, colour = Gene_GSVA_Metadata_DF$t_cell_qtle,
        HALLMARK_HYPOXIA,
        T_cell_exhaustion,
        xlab = "HALLMARK_HYPOXIA",
        ylab = "t_cell_qtle",
      main = "Hallmark Hypoxia vs. HIF1A Quantiles")+
       labs(colour = "t_cell_qtle")+
  theme_bw()

library(ggsurvfit)
cuminc(Surv(os_days, os_bin_numeric) ~ t_cell_qtle, data = Gene_GSVA_Metadata_DF) %>% 
  # ggcuminc(outcome = c(0, 1)) + # This line did not work 
  ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

3.5.2) EFS

library(tidycmprsk)
# r2_gse62564_GSVA_Metadata <- readRDS("../results/r2_gse62564_GSVA_Metadata.rds")

Gene_GSVA_Metadata_DF$efs_bin_numeric <- as.factor(Gene_GSVA_Metadata_DF$efs_bin_numeric)

Gene_GSVA_Metadata_DF$efs_days <- as.numeric(Gene_GSVA_Metadata_DF$efs_days)

library(ggsurvfit)
cuminc(Surv(efs_days, efs_bin_numeric) ~ t_cell_qtle, data = Gene_GSVA_Metadata_DF) %>% 
  # ggcuminc(outcome = c(0, 1)) + # This line did not work 
  ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

3.6) ELVIDGE_HIF1A_TARGETS_UP

3.6.1) OS

library(tidycmprsk)
# r2_gse62564_GSVA_Metadata <- readRDS("../results/r2_gse62564_GSVA_Metadata.rds")

Gene_GSVA_Metadata_DF$os_bin_numeric <- as.factor(Gene_GSVA_Metadata_DF$os_bin_numeric)

Gene_GSVA_Metadata_DF$os_days <- as.numeric(Gene_GSVA_Metadata_DF$os_days)

## HALLMARK_HYPOXIA has not been defined as numeric for the Gene_GSVA_Metadata_DF dataframe
Gene_GSVA_Metadata_DF$ELVIDGE_HIF1A_TARGETS_UP <- as.numeric(Gene_GSVA_Metadata_DF$ELVIDGE_HIF1A_TARGETS_UP)

qplot(data = Gene_GSVA_Metadata_DF, colour = Gene_GSVA_Metadata_DF$elvige_Up_qtle,
        HALLMARK_HYPOXIA,
        ELVIDGE_HIF1A_TARGETS_UP,
        xlab = "HALLMARK_HYPOXIA",
        ylab = "elvige_Up_qtle",
      main = "Hallmark Hypoxia vs. ELVIDGE_HIF1A_TARGETS_UP")+
       labs(colour = "elvige_Up_qtle")+
  theme_bw()

library(ggsurvfit)
cuminc(Surv(os_days, os_bin_numeric) ~ elvige_Up_qtle, data = Gene_GSVA_Metadata_DF) %>% 
  # ggcuminc(outcome = c(0, 1)) + # This line did not work 
  ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

3.6.2) EFS

library(tidycmprsk)
# r2_gse62564_GSVA_Metadata <- readRDS("../results/r2_gse62564_GSVA_Metadata.rds")

Gene_GSVA_Metadata_DF$efs_bin_numeric <- as.factor(Gene_GSVA_Metadata_DF$efs_bin_numeric)

Gene_GSVA_Metadata_DF$efs_days <- as.numeric(Gene_GSVA_Metadata_DF$efs_days)

library(ggsurvfit)
cuminc(Surv(efs_days, efs_bin_numeric) ~ elvige_Up_qtle, data = Gene_GSVA_Metadata_DF) %>% 
  # ggcuminc(outcome = c(0, 1)) + # This line did not work 
  ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

3.7) Quiescence Up

3.7.1) OS

library(tidycmprsk)
# r2_gse62564_GSVA_Metadata <- readRDS("../results/r2_gse62564_GSVA_Metadata.rds")

Gene_GSVA_Metadata_DF$os_bin_numeric <- as.factor(Gene_GSVA_Metadata_DF$os_bin_numeric)

Gene_GSVA_Metadata_DF$os_days <- as.numeric(Gene_GSVA_Metadata_DF$os_days)

## HALLMARK_HYPOXIA has not been defined as numeric for the Gene_GSVA_Metadata_DF dataframe
Gene_GSVA_Metadata_DF$Quiescence_Up_Cheung <- as.numeric(Gene_GSVA_Metadata_DF$Quiescence_Up_Cheung)

qplot(data = Gene_GSVA_Metadata_DF, colour = Gene_GSVA_Metadata_DF$qui_up_cheung,
        HALLMARK_HYPOXIA,
        Quiescence_Up_Cheung,
        xlab = "HALLMARK_HYPOXIA",
        ylab = "qui_up_cheung",
      main = "Hallmark Hypoxia vs. Quiescence Up Cheung")+
       labs(colour = "qui_Up_cheung")+
  theme_bw()

library(ggsurvfit)
cuminc(Surv(os_days, os_bin_numeric) ~ qui_up_cheung, data = Gene_GSVA_Metadata_DF) %>% 
  # ggcuminc(outcome = c(0, 1)) + # This line did not work 
  ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

3.7.2) EFS

library(tidycmprsk)
# r2_gse62564_GSVA_Metadata <- readRDS("../results/r2_gse62564_GSVA_Metadata.rds")

Gene_GSVA_Metadata_DF$efs_bin_numeric <- as.factor(Gene_GSVA_Metadata_DF$efs_bin_numeric)

Gene_GSVA_Metadata_DF$efs_days <- as.numeric(Gene_GSVA_Metadata_DF$efs_days)

library(ggsurvfit)
cuminc(Surv(efs_days, efs_bin_numeric) ~ qui_up_cheung, data = Gene_GSVA_Metadata_DF) %>% 
  # ggcuminc(outcome = c(0, 1)) + # This line did not work 
  ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

3.7) Quiescence Down

3.7.1) OS

library(tidycmprsk)
# r2_gse62564_GSVA_Metadata <- readRDS("../results/r2_gse62564_GSVA_Metadata.rds")

Gene_GSVA_Metadata_DF$os_bin_numeric <- as.factor(Gene_GSVA_Metadata_DF$os_bin_numeric)

Gene_GSVA_Metadata_DF$os_days <- as.numeric(Gene_GSVA_Metadata_DF$os_days)

## HALLMARK_HYPOXIA has not been defined as numeric for the Gene_GSVA_Metadata_DF dataframe
Gene_GSVA_Metadata_DF$Quiescence_Down_Cheung <- as.numeric(Gene_GSVA_Metadata_DF$Quiescence_Down_Cheung)

qplot(data = Gene_GSVA_Metadata_DF, colour = Gene_GSVA_Metadata_DF$qui_down_cheung,
        HALLMARK_HYPOXIA,
        Quiescence_Down_Cheung,
        xlab = "HALLMARK_HYPOXIA",
        ylab = "qui_down_cheung",
      main = "Hallmark Hypoxia vs. HIF1A Quantiles")+
       labs(colour = "qui_down_cheung")+
  theme_bw()

library(ggsurvfit)
cuminc(Surv(os_days, os_bin_numeric) ~ qui_down_cheung, data = Gene_GSVA_Metadata_DF) %>% 
  # ggcuminc(outcome = c(0, 1)) + # This line did not work 
  ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

3.7.2) EFS

library(tidycmprsk)
# r2_gse62564_GSVA_Metadata <- readRDS("../results/r2_gse62564_GSVA_Metadata.rds")

Gene_GSVA_Metadata_DF$efs_bin_numeric <- as.factor(Gene_GSVA_Metadata_DF$efs_bin_numeric)

Gene_GSVA_Metadata_DF$efs_days <- as.numeric(Gene_GSVA_Metadata_DF$efs_days)

library(ggsurvfit)
cuminc(Surv(efs_days, efs_bin_numeric) ~ qui_down_cheung, data = Gene_GSVA_Metadata_DF) %>% 
  # ggcuminc(outcome = c(0, 1)) + # This line did not work 
  ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

3.8) ST8SIA1/GD2 Surface Expression

3.8.1) OS

library(tidycmprsk)
# r2_gse62564_GSVA_Metadata <- readRDS("../results/r2_gse62564_GSVA_Metadata.rds")

Gene_GSVA_Metadata_DF$os_bin_numeric <- as.factor(Gene_GSVA_Metadata_DF$os_bin_numeric)

Gene_GSVA_Metadata_DF$os_days <- as.numeric(Gene_GSVA_Metadata_DF$os_days)

## HALLMARK_HYPOXIA has not been defined as numeric for the Gene_GSVA_Metadata_DF dataframe
Gene_GSVA_Metadata_DF$ST8SIA1 <- as.numeric(Gene_GSVA_Metadata_DF$ST8SIA1)

qplot(data = Gene_GSVA_Metadata_DF, colour = Gene_GSVA_Metadata_DF$st8sia1_qtle,
        HALLMARK_HYPOXIA,
        ST8SIA1,
        xlab = "HALLMARK_HYPOXIA",
        ylab = "st8sia1_qtle",
      main = "Hallmark Hypoxia vs. HIF1A Quantiles")+
       labs(colour = "st8sia1_qtle")+
  theme_bw()

library(ggsurvfit)
cuminc(Surv(os_days, os_bin_numeric) ~ st8sia1_qtle, data = Gene_GSVA_Metadata_DF) %>% 
  # ggcuminc(outcome = c(0, 1)) + # This line did not work 
  ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

3.8.2) EFS

library(tidycmprsk)
# r2_gse62564_GSVA_Metadata <- readRDS("../results/r2_gse62564_GSVA_Metadata.rds")

Gene_GSVA_Metadata_DF$efs_bin_numeric <- as.factor(Gene_GSVA_Metadata_DF$efs_bin_numeric)

Gene_GSVA_Metadata_DF$efs_days <- as.numeric(Gene_GSVA_Metadata_DF$efs_days)

library(ggsurvfit)
cuminc(Surv(efs_days, efs_bin_numeric) ~ st8sia1_qtle, data = Gene_GSVA_Metadata_DF) %>% 
  # ggcuminc(outcome = c(0, 1)) + # This line did not work 
  ggcuminc() +
  ylim(c(0, 0.7)) + 
  labs(
    x = "Days"
  )
## Plotting outcome "1".

4) References

Alboukadel Kassambara, Przemyslaw Biecek, Marcin Kosinski. 2023. Forest Plot for Cox Proportional Hazards Model. https://rpkgs.datanovia.com/survminer/reference/ggforest.html.
Kassambara, Alboukadel. 2023. Cox Proportional-Hazards Model. http://www.sthda.com/english/wiki/cox-proportional-hazards-model.
Marin, Mike, and Ladan Hamadani. 2020a. Kaplan Meier Vs. Exponential Vs. Cox Proportional Hazards (Pros & Cons). https://www.youtube.com/watch?v=K7bmmbD7KIg.
———. 2020b. Survival Analysis | Concepts and Implementation in r. https://www.youtube.com/playlist?list=PLqzoL9-eJTNDdnKvep_YHIwk2AMqHhuJ0.

5) Session Info

sessionInfo()
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggsurvfit_0.3.0  tidycmprsk_0.2.0 lubridate_1.9.2  forcats_1.0.0   
##  [5] stringr_1.5.0    purrr_1.0.1      readr_2.1.4      tidyr_1.3.0     
##  [9] tidyverse_2.0.0  expss_0.11.4     maditr_0.8.3     reshape_0.8.9   
## [13] kableExtra_1.3.4 tibble_3.2.1     dplyr_1.1.2      survminer_0.4.9 
## [17] ggpubr_0.6.0     ggplot2_3.4.2    survival_3.5-5  
## 
## loaded via a namespace (and not attached):
##  [1] ggtext_0.1.2         cmprsk_2.2-11        matrixStats_1.0.0   
##  [4] webshot_0.5.5        httr_1.4.6           tools_4.1.1         
##  [7] backports_1.4.1      bslib_0.5.0          utf8_1.2.3          
## [10] R6_2.5.1             colorspace_2.1-0     withr_2.5.0         
## [13] tidyselect_1.2.0     gridExtra_2.3        compiler_4.1.1      
## [16] cli_3.6.1            rvest_1.0.3          gt_0.9.0            
## [19] htmlTable_2.4.1      xml2_1.3.4           labeling_0.4.2      
## [22] sass_0.4.6           scales_1.2.1         checkmate_2.2.0     
## [25] survMisc_0.5.6       commonmark_1.9.0     systemfonts_1.0.4   
## [28] digest_0.6.32        rmarkdown_2.22       svglite_2.1.1       
## [31] pkgconfig_2.0.3      htmltools_0.5.5      fastmap_1.1.1       
## [34] highr_0.10           htmlwidgets_1.6.2    rlang_1.1.1         
## [37] rstudioapi_0.14      jquerylib_0.1.4      generics_0.1.3      
## [40] farver_2.1.1         zoo_1.8-12           jsonlite_1.8.7      
## [43] car_3.1-2            magrittr_2.0.3       huxtable_5.5.2      
## [46] Matrix_1.5-1         Rcpp_1.0.10          munsell_0.5.0       
## [49] fansi_1.0.4          abind_1.4-5          lifecycle_1.0.3     
## [52] stringi_1.7.12       yaml_2.3.7           gtsummary_1.7.1     
## [55] carData_3.0-5        plyr_1.8.8           grid_4.1.1          
## [58] crayon_1.5.2         lattice_0.21-8       cowplot_1.1.1       
## [61] splines_4.1.1        gridtext_0.1.5       hms_1.1.3           
## [64] knitr_1.43           pillar_1.9.0         markdown_1.7        
## [67] ggsignif_0.6.4       glue_1.6.2           evaluate_0.21       
## [70] broom.helpers_1.13.0 data.table_1.14.8    vctrs_0.6.3         
## [73] tzdb_0.4.0           gtable_0.3.3         km.ci_0.5-6         
## [76] assertthat_0.2.1     cachem_1.0.8         xfun_0.39           
## [79] xtable_1.8-4         broom_1.0.5          rstatix_0.7.2       
## [82] viridisLite_0.4.2    hardhat_1.3.0        KMsurv_0.1-5        
## [85] timechange_0.2.0     ellipsis_0.3.2
end_time <- Sys.time()
end_time - start_time
## Time difference of 2.934798 mins