References

The survival analysis reference used was writen by Zabor

  • (Zabor 2023).

Transfer libraries

The survival analysis reference used was Jaeger.

  • (Jaeger 2015).

To add p-values to plots, I used the Kassambara 2017 and 2020 references

  • (Kassambara 2017)

Load libraries

library(maditr)
## Warning: package 'maditr' was built under R version 4.1.2
## 
## To modify variables or add new variables:
##              let(mtcars, new_var = 42, new_var2 = new_var*hp) %>% head()

1) Continued processing of the R2 dataset

Load the GSVA object back to the R environment. This object now contains all the metadata information.

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

In neuroblastoma, the main risk classification system classifies patients using features such as age, stage, histology and MYCN amplification status.

MYCN is a gene that encodes a transcription factor that is very important in tumor prognosis, the association between tumor features and patient survival.

Let’s take a look at the MYCN expression compared by the neuroblastoma risk group.

First, let’s see which type of variable is the the MYCN status and the MYCN expression

class(r2_gse62564_GSVA_Metadata$mycn_status)
## [1] "character"
class(r2_gse62564_GSVA_Metadata$MYCN)
## [1] "character"

They are both of type “character.” In order to plot, we have to make the MYCN expression to be of type numeric

r2_gse62564_GSVA_Metadata$MYCN <- as.numeric(r2_gse62564_GSVA_Metadata$MYCN)
class(r2_gse62564_GSVA_Metadata$MYCN)
## [1] "numeric"

Check the age at diagnosis

r2_gse62564_GSVA_Metadata$age_at_diagnosis <- as.numeric(r2_gse62564_GSVA_Metadata$age_at_diagnosis)
age_diagnosis_days_gse62564_498 <- mean(r2_gse62564_GSVA_Metadata$age_at_diagnosis)

2) First Plots of the Neuroblastoma R2 Dataset

2.1) MYCN Amplification status and MYCN expression

Let’s try and see the boxplot of the MYCN expression vs. MYCN amplification status.

Using the boxplot function that we tried before, we will have an error.

boxplot(data=r2_gse62564_GSVA_Metadata, r2_gse62564_GSVA_Metadata$mycn_status, r2_gse62564_GSVA_Metadata$MYCN)

Maybe we could try the boxplot function making the mycn_status a numeric variable.

boxplot(data=r2_gse62564_GSVA_Metadata, as.numeric(r2_gse62564_GSVA_Metadata$mycn_status), r2_gse62564_GSVA_Metadata$MYCN)

Let’s try and see if ggplot, which we have also used before, takes the mycn_status variable as categorical as it is right now.

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.1.2
ggplot(r2_gse62564_GSVA_Metadata, aes(x=mycn_status, y=MYCN, fill=mycn_status)) + 
    geom_boxplot()

In this figure, we see that group with clinical MYCN amplification, also shows high MYCN expression

What is the difference, in measurement, between MYCN amplification and MYCN expression?

How are the MYCN amplification and MYCN gene expression correlated?

For now, I am not bothering neither with formating of the X and Y labels nor with the N/A values of the MYCN amplification status.

2.2) High Risk status and MYCN expression

We can plot the high risk status vs. MYCN amplification and leave the fill argument as mycn_status

library(ggplot2)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=MYCN, fill=mycn_status)) + 
    geom_boxplot()

Or we can make it the high_risk status instead

library(ggplot2)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=MYCN, fill=high_risk)) + 
    geom_boxplot()

Here we see that the high_risk status correlates positively with MYCN expression.

3) Survival Analysis

Make the OS (overall survival) status numeric, instead of character. Do the same for the EFS.

## OS
r2_gse62564_GSVA_Metadata <- r2_gse62564_GSVA_Metadata %>%
  dplyr::mutate(os_bin_numeric =
                  ifelse(r2_gse62564_GSVA_Metadata$os_bin == "event", "1", "0"))

r2_gse62564_GSVA_Metadata["os_bin_numeric"] <- as.numeric(r2_gse62564_GSVA_Metadata$os_bin_numeric)

## EFS
r2_gse62564_GSVA_Metadata <- r2_gse62564_GSVA_Metadata %>%
  dplyr::mutate(efs_bin_numeric =
                  ifelse(r2_gse62564_GSVA_Metadata$efs_bin == "event", "1", "0"))

r2_gse62564_GSVA_Metadata["efs_bin_numeric"] <- as.numeric(r2_gse62564_GSVA_Metadata$efs_bin_numeric)

3.1) Survival Plots

library(survival)
## Warning: package 'survival' was built under R version 4.1.2
fit_os_high_risk <- survfit(Surv(as.numeric(os_day), os_bin_numeric) ~ r2_gse62564_GSVA_Metadata$high_risk,
                   data = r2_gse62564_GSVA_Metadata)
library(survminer)
## Loading required package: ggpubr
## Warning: package 'ggpubr' was built under R version 4.1.2
## 
## Attaching package: 'survminer'
## The following object is masked from 'package:survival':
## 
##     myeloma
ggsurvplot(fit_os_high_risk, 
           data = r2_gse62564_GSVA_Metadata, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(2000, 1),
           )

ggsurvplot(fit_os_high_risk, 
           data = r2_gse62564_GSVA_Metadata, 
           # risk.table = TRUE,
           risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(2000, 1),
           )

  • Have students generate the EFS plot

4) The Hypoxia Effect “Controversy”

Let’s now establish a relationship between neuroblastoma high risk status and hypoxia, one of the novelties of our work.

4.1) HIF1A

r2_gse62564_GSVA_Metadata$HIF1A <- as.numeric(r2_gse62564_GSVA_Metadata$HIF1A)

compare_means(data = r2_gse62564_GSVA_Metadata, HIF1A ~ high_risk, method = "t.test")
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=HIF1A, fill=high_risk)) + 
    geom_boxplot()+
  stat_compare_means(method = "t.test", label.x = 1.4)+
  theme_linedraw()

4.2) HALLMARK_HYPOXIA

We also must make HALLMARK_HYPOXIA, numeric

r2_gse62564_GSVA_Metadata$HALLMARK_HYPOXIA <- as.numeric(r2_gse62564_GSVA_Metadata$HALLMARK_HYPOXIA)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=HALLMARK_HYPOXIA, fill=high_risk)) + 
    geom_boxplot()

Include statistics analysis

In the next chunk, we construct a data-frame with statistics calculations for our plot.

The reference used here was the Kassambara 2023 article:

  • (Kassambara 2023)
library(rstatix)
## Warning: package 'rstatix' was built under R version 4.1.2
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
# Statistical test
stat.test <- r2_gse62564_GSVA_Metadata %>%
  t_test(HALLMARK_HYPOXIA ~ high_risk) %>%
  add_significance()
stat.test

In the next chunk, we create the object that will show the statistical results. After, we use the stat_pvalue_manual function.

We also change the theme of the ggplot2 figure plot.

library(ggpubr)

r2_gse62564_GSVA_Metadata$high_risk <- as.factor(r2_gse62564_GSVA_Metadata$high_risk)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=HALLMARK_HYPOXIA)) + 
    geom_boxplot()+
  stat_pvalue_manual(stat.test, label = "p", y.position = 0.6) + # y.position is mandatory
  scale_y_continuous(expand = expansion(mult = c(0.05, 0.1)))

r2_gse62564_GSVA_Metadata$high_risk <- as.factor(r2_gse62564_GSVA_Metadata$high_risk)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=HALLMARK_HYPOXIA, fill="high_risk")) + 
    geom_boxplot()+
  stat_pvalue_manual(stat.test, label = "p", y.position = 0.6) + # y.position is mandatory
  scale_y_continuous(expand = expansion(mult = c(0.05, 0.1)))

compare_means(data = r2_gse62564_GSVA_Metadata, HALLMARK_HYPOXIA ~ high_risk, method = "t.test")
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=HALLMARK_HYPOXIA, fill=high_risk)) + 
    geom_boxplot()+
  stat_compare_means(method = "t.test")+
  theme_linedraw()

ELVIDGE_HIF1A_TARGETS_UP

r2_gse62564_GSVA_Metadata$ELVIDGE_HIF1A_TARGETS_UP <- as.numeric(r2_gse62564_GSVA_Metadata$ELVIDGE_HIF1A_TARGETS_UP)

compare_means(data = r2_gse62564_GSVA_Metadata, ELVIDGE_HIF1A_TARGETS_UP ~ high_risk, method = "t.test")
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=ELVIDGE_HIF1A_TARGETS_UP, fill=high_risk)) + 
    geom_boxplot()+
  stat_compare_means(method = "t.test")+
  theme_linedraw()

GROSS_HYPOXIA_VIA_HIF1A_UP

r2_gse62564_GSVA_Metadata$GROSS_HYPOXIA_VIA_HIF1A_UP <- as.numeric(r2_gse62564_GSVA_Metadata$GROSS_HYPOXIA_VIA_HIF1A_UP)

compare_means(data = r2_gse62564_GSVA_Metadata, GROSS_HYPOXIA_VIA_HIF1A_UP ~ high_risk, method = "t.test")
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=GROSS_HYPOXIA_VIA_HIF1A_UP, fill=high_risk)) + 
    geom_boxplot()+
  stat_compare_means(method = "t.test")+
  theme_linedraw()

Hypoxia trial

r2_gse62564_GSVA_Metadata$GROSS_HYPOXIA_VIA_HIF1A_DN <- as.numeric(r2_gse62564_GSVA_Metadata$GROSS_HYPOXIA_VIA_HIF1A_DN)

compare_means(data = r2_gse62564_GSVA_Metadata, GROSS_HYPOXIA_VIA_HIF1A_DN ~ high_risk, method = "t.test")
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=GROSS_HYPOXIA_VIA_HIF1A_DN, fill=high_risk)) + 
    geom_boxplot()+
  stat_compare_means(method = "t.test")+
  theme_linedraw()

4.3) Hypoxia Down-Regulation

No statistical analysis

r2_gse62564_GSVA_Metadata$ADRN_Norm_vs_Hypo_Down_635.txt <- as.numeric(r2_gse62564_GSVA_Metadata$ADRN_Norm_vs_Hypo_Down_635.txt)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=ADRN_Norm_vs_Hypo_Down_635.txt, fill=high_risk)) + 
    geom_boxplot()

Statistical analysis

compare_means(data = r2_gse62564_GSVA_Metadata, ADRN_Norm_vs_Hypo_Down_635.txt ~ high_risk, method = "t.test")
r2_gse62564_GSVA_Metadata$ADRN_Norm_vs_Hypo_Down_635.txt <- as.numeric(r2_gse62564_GSVA_Metadata$ADRN_Norm_vs_Hypo_Down_635.txt)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=ADRN_Norm_vs_Hypo_Down_635.txt, fill=high_risk)) + 
    geom_boxplot()+
  stat_compare_means(method = "t.test")+
  theme_linedraw()

4.4) T cell exhaustion

Reference paper for the use of the T cell exhaustion gene set in this section was Chen 2020:

  • (Baotao Chen 2020)
r2_gse62564_GSVA_Metadata$`T cell exhaustion` <- as.numeric(r2_gse62564_GSVA_Metadata$`T cell exhaustion`)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=`T cell exhaustion`, fill=high_risk)) + 
    geom_boxplot()

Statistical analysis

The compare_means function can not deal with spaces in the T cell exhaustion variable name.

We will have to use a trick to deal with the spaces in “T cell exhaustion.”

We will use the package stringr to insert dots (.) replacing the spaces in the r2_gse62564_GSVA_Metadata dataframe.

To start, we have to create the r2_gse62564_GSVA_Metadata_Spaces dataframe, which will have structure names that contained spaces, changed to one of these new characteres: “.” , “,” or ““.

r2_gse62564_GSVA_Metadata_Spaces <- r2_gse62564_GSVA_Metadata
library(stringr)
## Warning: package 'stringr' was built under R version 4.1.2
names(r2_gse62564_GSVA_Metadata_Spaces)<-str_replace_all(names(r2_gse62564_GSVA_Metadata_Spaces), c(" " = "." , "," = "" ))

Then, instead of using the r2_gse62564_GSVA_Metadata dataframe, we use r2_gse62564_GSVA_Metadata_Spaces. Note that now we will use the variable name as T.cell.exhaustion, without the spaces.

compare_means(data = r2_gse62564_GSVA_Metadata_Spaces, T.cell.exhaustion ~ high_risk, method = "t.test")
r2_gse62564_GSVA_Metadata_Spaces$T.cell.exhaustion <- as.numeric(r2_gse62564_GSVA_Metadata_Spaces$T.cell.exhaustion)
ggplot(r2_gse62564_GSVA_Metadata_Spaces, aes(x=high_risk, y=T.cell.exhaustion, fill=high_risk)) + 
    geom_boxplot()+
  stat_compare_means(method = "t.test")+
  theme_linedraw()

4.5) Quiescence

Hypoxia induces quiescence, the cell cycle arrest at the G0 phase (Barbara Muz 2015).

We found a quiescence signature on the work that discussed signatures of quiescent stem cells (Cheung and Rando 2013).

HIF1A activation in cancer cells augments T cell exhaustion and reduces T cell killing (Pilar Baldominos 2022).

In another example, we interrogate about the Quiescence gene signatures.

4.5.1) Quiescence Up-Regulated

r2_gse62564_GSVA_Metadata_Spaces$Quiescence.Up.Cheung <- as.numeric(r2_gse62564_GSVA_Metadata_Spaces$Quiescence.Up.Cheung)

compare_means(data = r2_gse62564_GSVA_Metadata_Spaces, Quiescence.Up.Cheung ~ high_risk, method = "t.test")
r2_gse62564_GSVA_Metadata_Spaces$Quiescence.Up.Cheung <- as.numeric(r2_gse62564_GSVA_Metadata_Spaces$Quiescence.Up.Cheung)
ggplot(r2_gse62564_GSVA_Metadata_Spaces, aes(x=high_risk, y=Quiescence.Up.Cheung, fill=high_risk)) + 
    geom_boxplot()+
  stat_compare_means(method = "t.test")+
  theme_linedraw()

4.5.2) Quiescence Down-Regulated

r2_gse62564_GSVA_Metadata_Spaces$Quiescence.Down.Cheung <- as.numeric(r2_gse62564_GSVA_Metadata_Spaces$Quiescence.Down.Cheung)

compare_means(data = r2_gse62564_GSVA_Metadata_Spaces, Quiescence.Down.Cheung ~ high_risk, method = "t.test")
r2_gse62564_GSVA_Metadata_Spaces$Quiescence.Down.Cheung <- as.numeric(r2_gse62564_GSVA_Metadata_Spaces$Quiescence.Down.Cheung)
ggplot(r2_gse62564_GSVA_Metadata_Spaces, aes(x=high_risk, y=Quiescence.Down.Cheung, fill=high_risk)) + 
    geom_boxplot()+
  stat_compare_means(method = "t.test")+
  theme_linedraw()

4.6) HALLMARK_INFLAMMATORY_RESPONSE

No statistical analysis

r2_gse62564_GSVA_Metadata$HALLMARK_INFLAMMATORY_RESPONSE <- as.numeric(r2_gse62564_GSVA_Metadata$HALLMARK_INFLAMMATORY_RESPONSE)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=HALLMARK_INFLAMMATORY_RESPONSE, fill=high_risk)) + 
    geom_boxplot()

Statistical analysis

compare_means(data = r2_gse62564_GSVA_Metadata, HALLMARK_INFLAMMATORY_RESPONSE ~ high_risk, method = "t.test")
r2_gse62564_GSVA_Metadata$HALLMARK_INFLAMMATORY_RESPONSE <- as.numeric(r2_gse62564_GSVA_Metadata$HALLMARK_INFLAMMATORY_RESPONSE)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=HALLMARK_INFLAMMATORY_RESPONSE, fill=high_risk)) + 
    geom_boxplot()+
  stat_compare_means(method = "t.test")+
  theme_linedraw()

4.7) ST8SIA1/GD2

No Statistical Analysis

Reference paper for this section is Mabe 2022:

  • (Nathaniel W Mabe 2022)
r2_gse62564_GSVA_Metadata$ST8SIA1 <- as.numeric(r2_gse62564_GSVA_Metadata$ST8SIA1)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=ST8SIA1, fill=high_risk)) + 
    geom_boxplot()

Statistical analysis

compare_means(data = r2_gse62564_GSVA_Metadata, ST8SIA1 ~ high_risk, method = "t.test")
r2_gse62564_GSVA_Metadata$ST8SIA1 <- as.numeric(r2_gse62564_GSVA_Metadata$ST8SIA1)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=ST8SIA1, fill=high_risk)) + 
    geom_boxplot()+
  stat_compare_means(method = "t.test")+
  theme_linedraw()

5) Survival analysis of the hypoxia gene sets

5.1) HIF1A

OS

library(tidyverse)

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

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

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

library(survival)
fit_os_HIF1A <- survfit(Surv(as.numeric(os_day), os_bin_numeric) ~ r2_gse62564_GSVA_Metadata_survival$quantile,
                   data = r2_gse62564_GSVA_Metadata_survival)
library(survminer)
ggsurvplot(fit_os_HIF1A, 
           data = r2_gse62564_GSVA_Metadata_survival, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

EFS

library(tidyverse)

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

library(survival)
fit_efs_HIF1A <- survfit(Surv(as.numeric(efs_day), efs_bin_numeric) ~ r2_gse62564_GSVA_Metadata_survival$quantile,
                   data = r2_gse62564_GSVA_Metadata_survival)
library(survminer)
ggsurvplot(fit_efs_HIF1A, 
           data = r2_gse62564_GSVA_Metadata_survival, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

5.2) ELVIDGE_HIF1A_TARGETS_UP

library(tidyverse)

r2_gse62564_GSVA_Metadata_survival <- r2_gse62564_GSVA_Metadata %>%
    mutate(quantile = ntile(ELVIDGE_HIF1A_TARGETS_UP, 4))

library(survival)
fit_os_ELVIDGE_HIF1A_TARGETS_UP <- survfit(Surv(as.numeric(os_day), os_bin_numeric) ~ r2_gse62564_GSVA_Metadata_survival$quantile,
                   data = r2_gse62564_GSVA_Metadata_survival)
library(survminer)
ggsurvplot(fit_os_ELVIDGE_HIF1A_TARGETS_UP, 
           data = r2_gse62564_GSVA_Metadata_survival, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

Another solution to this is using the ifelse statement to generate the gene set score as a categorical vector.

Once the categorical vector is generated, it can be used to be plotted in the Kaplan Meyer Survival curve.

5.3) T.cell.exhaustion (OS)

library(tidyverse)

r2_gse62564_GSVA_Metadata_Spaces_survival <- r2_gse62564_GSVA_Metadata_Spaces %>%
    mutate(quantile = ntile(T.cell.exhaustion, 4))

library(survival)
fit_os_T.cell.exhaustion <- survfit(Surv(as.numeric(os_day), os_bin_numeric) ~ r2_gse62564_GSVA_Metadata_Spaces_survival$quantile,
                   data = r2_gse62564_GSVA_Metadata_Spaces_survival)
library(survminer)
ggsurvplot(fit_os_T.cell.exhaustion, 
           data = r2_gse62564_GSVA_Metadata_Spaces_survival, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

5.4) T.cell.exhaustion (EFS)

library(tidyverse)

r2_gse62564_GSVA_Metadata_Spaces_survival <- r2_gse62564_GSVA_Metadata_Spaces %>%
    mutate(quantile = ntile(T.cell.exhaustion, 2))

library(survival)
fit_efs_T.cell.exhaustion <- survfit(Surv(as.numeric(efs_day), efs_bin_numeric) ~ r2_gse62564_GSVA_Metadata_Spaces_survival$quantile,
                   data = r2_gse62564_GSVA_Metadata_Spaces_survival)
library(survminer)
ggsurvplot(fit_efs_T.cell.exhaustion, 
           data = r2_gse62564_GSVA_Metadata_Spaces_survival, 
           risk.table = TRUE,
           # risk.table = "abs_pct",
           pval = TRUE, 
           pval.coord = c(4000, 1),
           )

6) HIF1A relationship to age at diagnosis

6.1) HIF1A, age at diagnosis and HR

r2_gse62564_GSVA_Metadata_Spaces_survival$age_at_diagnosis <- as.numeric(r2_gse62564_GSVA_Metadata_Spaces_survival$age_at_diagnosis)

max(r2_gse62564_GSVA_Metadata_Spaces_survival$age_at_diagnosis)
## [1] 8983
min(r2_gse62564_GSVA_Metadata_Spaces_survival$age_at_diagnosis)
## [1] 0
r2_gse62564_GSVA_Metadata_Spaces_survival_age <- r2_gse62564_GSVA_Metadata_Spaces_survival %>%
    mutate(age_qtls = ntile(age_at_diagnosis, 4))

r2_gse62564_GSVA_Metadata_Spaces_survival_age$age_qtls <- as.factor(r2_gse62564_GSVA_Metadata_Spaces_survival_age$age_qtls)
ggplot(r2_gse62564_GSVA_Metadata_Spaces_survival_age, aes(x=age_qtls, y=HIF1A, fill=high_risk)) + 
    geom_boxplot()+
  stat_compare_means(method = "t.test")+
  theme_linedraw()

6.2) HIF1A, age at diagnosis and MYCN status

r2_gse62564_GSVA_Metadata_Spaces_survival_age <- r2_gse62564_GSVA_Metadata_Spaces_survival %>%
    mutate(age_qtls = ntile(age_at_diagnosis, 4))

r2_gse62564_GSVA_Metadata_Spaces_survival_age$age_qtls <- as.factor(r2_gse62564_GSVA_Metadata_Spaces_survival_age$age_qtls)
ggplot(r2_gse62564_GSVA_Metadata_Spaces_survival_age, aes(x=age_qtls, y=HIF1A, fill=mycn_status)) + 
    geom_boxplot()+
  stat_compare_means(method = "anova")+
  theme_linedraw()

Remove the NAs from the mycn_status variable

r2_gse62564_GSVA_Metadata_Spaces_survival_age <- r2_gse62564_GSVA_Metadata_Spaces_survival_age[which(
  r2_gse62564_GSVA_Metadata_Spaces_survival_age$mycn_status != "na"  ), ]

r2_gse62564_GSVA_Metadata_Spaces_survival_age <- r2_gse62564_GSVA_Metadata_Spaces_survival_age %>%
    mutate(age_qtls = ntile(age_at_diagnosis, 4))

r2_gse62564_GSVA_Metadata_Spaces_survival_age$age_qtls <- as.factor(r2_gse62564_GSVA_Metadata_Spaces_survival_age$age_qtls)
ggplot(r2_gse62564_GSVA_Metadata_Spaces_survival_age, aes(x=age_qtls, y=HIF1A, fill=mycn_status)) + 
    geom_boxplot()+
  stat_compare_means(method = "anova")+
  theme_linedraw()

7) References

Baotao Chen, Mengyuan Li, Lin Li. 2020. Hif1a Expression Correlates with Increased Tumor Immune and Stromal Signatures and Aggressive Phenotypes in Human Cancers. https://pubmed.ncbi.nlm.nih.gov/32488852/.
Barbara Muz, Feda Azab, Pilar de la Puente. 2015. The Role of Hypoxia in Cancer Progression, Angiogenesis, Metastasis, and Resistance to Therapy. https://pubmed.ncbi.nlm.nih.gov/27774485/.
Cheung, Tom H., and Thomas A. Rando. 2013. Molecular Regulation of Stem Cell Quiescence. https://pubmed.ncbi.nlm.nih.gov/23698583/.
Jaeger, T. Florian. 2015. Transfer Installed Libraries to a Different Version of r. https://hlplab.wordpress.com/2012/06/01/transferring-installed-packages-between-different-installations-of-r/.
Kassambara, Alboukadel. 2017. Add p-Values and Significance Levels to Ggplots. http://www.sthda.com/english/articles/24-ggpubr-publication-ready-plots/76-add-p-values-and-significance-levels-to-ggplots/.
———. 2023. How to Add p-Values onto Basic GGPLOTS. https://www.datanovia.com/en/blog/how-to-add-p-values-onto-basic-ggplots/.
Nathaniel W Mabe, Guillermo N Dalton, Min Huang. 2022. Transition to a Mesenchymal State in Neuroblastoma Confers Resistance to Anti-Gd2 Antibody via Reduced Expression of St8sia1. https://pubmed.ncbi.nlm.nih.gov/35817829/.
Pilar Baldominos, Olga Barreiro, Alex Barbera-Mourelle. 2022. Quiescent Cancer Cells Resist t Cell Attack by Forming an Immunosuppressive Niche. https://pubmed.ncbi.nlm.nih.gov/35447074/.
Zabor, Emily C. 2023. Survival Analysis in r. https://www.emilyzabor.com/tutorials/survival_analysis_in_r_tutorial.html.

8) 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] lubridate_1.9.2 forcats_1.0.0   dplyr_1.1.2     purrr_1.0.1    
##  [5] readr_2.1.4     tidyr_1.3.0     tibble_3.2.1    tidyverse_2.0.0
##  [9] stringr_1.5.0   rstatix_0.7.2   survminer_0.4.9 ggpubr_0.6.0   
## [13] survival_3.5-5  ggplot2_3.4.2   maditr_0.8.3   
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.10       lattice_0.21-8    zoo_1.8-12        digest_0.6.32    
##  [5] utf8_1.2.3        R6_2.5.1          backports_1.4.1   evaluate_0.21    
##  [9] highr_0.10        pillar_1.9.0      rlang_1.1.1       rstudioapi_0.14  
## [13] data.table_1.14.8 car_3.1-2         jquerylib_0.1.4   Matrix_1.5-1     
## [17] rmarkdown_2.22    labeling_0.4.2    splines_4.1.1     munsell_0.5.0    
## [21] gridtext_0.1.5    broom_1.0.5       compiler_4.1.1    xfun_0.39        
## [25] pkgconfig_2.0.3   htmltools_0.5.5   ggtext_0.1.2      tidyselect_1.2.0 
## [29] gridExtra_2.3     km.ci_0.5-6       fansi_1.0.4       tzdb_0.4.0       
## [33] withr_2.5.0       commonmark_1.9.0  grid_4.1.1        jsonlite_1.8.7   
## [37] xtable_1.8-4      gtable_0.3.3      lifecycle_1.0.3   magrittr_2.0.3   
## [41] KMsurv_0.1-5      scales_1.2.1      stringi_1.7.12    cli_3.6.1        
## [45] cachem_1.0.8      carData_3.0-5     farver_2.1.1      ggsignif_0.6.4   
## [49] xml2_1.3.4        bslib_0.5.0       survMisc_0.5.6    generics_0.1.3   
## [53] vctrs_0.6.3       tools_4.1.1       glue_1.6.2        markdown_1.7     
## [57] hms_1.1.3         abind_1.4-5       fastmap_1.1.1     yaml_2.3.7       
## [61] timechange_0.2.0  colorspace_2.1-0  knitr_1.43        sass_0.4.6