Analysis of survival data for Ewing sarcoma

For this analysis, we will start by getting the data from GEO, processing it into a more suitable object for a survival analysis.

Setting up the environment

These are packages you will need for this notebook. For exact versions used, please refer to the session info at the bottom of this notebook.

library('tidyverse')
library('GEOquery')
library('ggplot2')
library('RColorBrewer')
library('survival')
library('survminer')

I usually set a couple of directories of where to store things, or where I have located things we may need. Below, generalDatasets is just where we have stored more global data, such as a general use database. The baseWorkspace is where we will store my coding files that will usually end up getting pushed to GitHub. The last is baseRepository. This is set up the same as baseWorkspace but is not pushed to GitHub, so will store larger files (such as the data we download from GEO).

generalDatasets = 'C:/Users/chris/OneDrive/Documents/bccrc/projectsRepository/generalDatasets'
baseWorkspace = 'C:/Users/chris/OneDrive/Documents/bccrc/protocolsDryLab/relatedToPublishedData'
baseRepository = 'C:/Users/chris/OneDrive/Documents/bccrc/protocolsRepository'

Data processing

This data set is from Savola et al. from GSE17679 that was published in PMID 22084725. There appears to be two sets of replicates here, with a bunch of other cell line and other tissue samples analyzed as well that may need filtering downstream. If you wanted to analyze different data here, replace with your own GEO identifier. The destdir portion of the command here is telling getGEO where we want to save the data. Here we are using the ‘paste’ function to combine my baseRepository from above with the name of the folder where we want to save it. The sep = ’’ just tells it to leave no spaces between the individual elements (it is a file path, we don’t want spaces).

##########################################################################################
esArrayData = getGEO('GSE17679',
             destdir = paste(baseRepository, '/getGeoDataForSurvivalAnalysis/', sep = ''),
             GSEMatrix = TRUE)

Found 1 file(s)
GSE17679_series_matrix.txt.gz
Using locally cached version: C:/Users/chris/OneDrive/Documents/bccrc/protocolsRepository/getGeoDataForSurvivalAnalysis//GSE17679_series_matrix.txt.gz
Parsed with column specification:
cols(
  .default = col_double(),
  ID_REF = col_character()
)
See spec(...) for full column specifications.
Using locally cached version of GPL570 found here:
C:/Users/chris/OneDrive/Documents/bccrc/protocolsRepository/getGeoDataForSurvivalAnalysis//GPL570.soft 
62 parsing failures.
  row     col           expected    actual         file
54614 SPOT_ID 1/0/T/F/TRUE/FALSE --Control literal data
54615 SPOT_ID 1/0/T/F/TRUE/FALSE --Control literal data
54616 SPOT_ID 1/0/T/F/TRUE/FALSE --Control literal data
54617 SPOT_ID 1/0/T/F/TRUE/FALSE --Control literal data
54618 SPOT_ID 1/0/T/F/TRUE/FALSE --Control literal data
..... ....... .................. ......... ............
See problems(...) for more details.

This command returns a list like object, so we will extract it to be an expressionSet.

##########################################################################################
esData = esArrayData[[1]]

We now have an ‘expressionSet’ that contains different types of information. We can look at some of the data to see what we have using different commands. Below, pData gives us the ‘phenotype data’ - this is usually the sample details provided by the submitter and is often where the clinical information is found.

##########################################################################################
head(pData(esData))

The exprs gives us the expression data for each sample analyzed, and all probes.

##########################################################################################
head(exprs(esData)[,1:5])

          GSM439886 GSM439887 GSM439888 GSM439889 GSM439890
1007_s_at  9.151022  8.160573  9.016025  8.609008  9.199717
1053_at    6.085073  5.719322  6.775956  4.984402  6.582644
117_at     4.108801  4.166347  4.740777  4.149291  4.202501
121_at     7.121038  7.437739  6.235114  6.537358  6.519501
1255_g_at  2.735556  2.889045  3.078523  2.681159  4.497911
1294_at    6.426440  6.113979  5.224356  6.830285  5.732191

The fData will give us ‘feature data’, which is usually related to the array itself, so things like probe identifiers and annotation.

##########################################################################################
head(fData(esData))

First we just want to annotate the expression data to include gene symbols (right now it is just probe identifiers), as these are what we will use for the survival analysis downstream. We can use the feature data to get this annotation. Below, we extract the expression data into the ‘esExpression’ object.

##################################################################################
esExpression = as_tibble(exprs(esData)) %>%
  mutate(ID = row.names(exprs(esData)))

Now, we extract the feature data, keep the first gene symbol in the ‘Gene Symbol’ column, select only the columns we want to keep with the ‘select’ function, and then merge it with the esExpression object we just created.

##################################################################################
esExpressionAnnotated = as_tibble(fData(esData)) %>%
  mutate(symbol = sapply(strsplit(`Gene Symbol`, ' '), '[', 1)) %>%
  dplyr::select(ID, GB_ACC, symbol) %>%
  left_join(esExpression)

Joining, by = "ID"

So now we have an expression object with gene symbols and all of the associated data. Now I want to look at the phenotype data. There appears to be PNET cases in here, cell lines, and some other tissues. I will keep only things classified as Ewing. Most of the command below is just text parsing to extract the information we want. This is the purpose of the ‘sub’ function. I know it can look a bit confusing, but if you want to learn more about it, you should look into regular expressions or ask here and I can do my best to explain it.

##################################################################################
esPhenotype = as_tibble(pData(esData)) %>%
  dplyr::select(geo_accession, characteristics_ch1.1:characteristics_ch1.8) %>%
  filter(grepl('Ewing', characteristics_ch1.1)) %>%
  mutate(type = sub('state\\: (.*)$', '\\1', characteristics_ch1.3)) %>%
  mutate(age = sub('age\\: (.*)$', '\\1', characteristics_ch1.4)) %>%
  mutate(sex = sub('sex\\: (.*)$', '\\1', characteristics_ch1.5)) %>%
  mutate(efs = as.double(sub('efs \\(months\\)\\: (.*)$', '\\1', characteristics_ch1.6))) %>%
  mutate(ovs = as.double(sub('ovs \\(months\\)\\: (.*)$', '\\1', characteristics_ch1.7))) %>%
  mutate(status = sub('status\\: (.*)$', '\\1', characteristics_ch1.8)) %>%
  dplyr::select(geo_accession, type:status)

Now we assign a numerical value to alive and dead, where dead = 2.

##################################################################################
esPhenotype$status = ifelse(esPhenotype$status == 'Dead', 2, 1)

Now we can save the two data sets, just in case we want to use them later on for something else and don’t want to redo this processing.

##################################################################################
saveRDS(esExpressionAnnotated, paste(baseRepository, '/getGeoDataForSurvivalAnalysis/dataset_savolaExpression.rds', sep = ''))
saveRDS(esPhenotype, paste(baseRepository, '/getGeoDataForSurvivalAnalysis/dataset_savolaPhenotype.rds', sep = ''))

Now we are ready to move on and do the survival analysis. For simplicity and to avoid messing up my original files, I am first going to read in the data we just saved to and start fresh from here with new names.

##########################################################################################
esExpression = readRDS(paste(baseRepository, '/getGeoDataForSurvivalAnalysis/dataset_savolaExpression.rds', sep = ''))
esPhenotype = readRDS(paste(baseRepository, '/getGeoDataForSurvivalAnalysis/dataset_savolaPhenotype.rds', sep = ''))

Now everything is set to do a survival analysis. First, we can just do this simply for a single gene. I find the easiest way to do this is to write a function as we want it to be reproducible if we decide to do it for many genes. For this, I like to create a new .R file, write the function in it, and call it into this session. So, if you are in RStudio, click File>New File>R Script. In my case, I have provided this file on GitHub and it is called getGeoDataForSurvivalAnalysisFunctions.R. You can open this file with any text editor, or RStudio. So, create the file, and read it in to this R session.

##########################################################################################
source(paste(baseWorkspace, '/getGeoDataForSurvivalAnalysisFunctions.R', sep = ''))

In this file, you will find a function called ‘survivalAnalysis’. This function will take your gene of interest, your phenotype and expression data, as well as inputs as to whether you would like to generate a plot or not. The function expects the data to be in the specific format as we have created here, so be careful of that in your own use. Expression data:

##########################################################################################
head(esExpression)

NA

And the phenotype data:

##########################################################################################
head(esPhenotype)

Now we can call our function on a gene we are interested in. By default, this function will write plot and text files to our baseRepository directory. The ‘writeDirectory’ option only allows you to change folders within the baseRepository location. You do not have to specify ‘writeDirectory’ if the baseRepository is suitable for you.

##########################################################################################
geneOfInterest = 'YBX1'
singleSurvival = survivalAnalysis(geneOfInterest = geneOfInterest,
                                expressionData = esExpression, 
                                phenotypeData = esPhenotype, 
                                survivalPlot = TRUE, 
                                survivalData = TRUE,
                                writeDirectory = '/getGeoDataForSurvivalAnalysis')

[1] "Performing analysis for YBX1"

A couple of things of note about how this function operates:

• For a given input gene, it will first parse the expression data to retain only probes for that gene. Because there are usually many probes for a gene, it ranks them and keeps only the highest expressing probe.
• After this, the expression values across all patients for this probe are ranked, and the top and bottom 25% of patients are assigned as ‘high’ and ‘low’ expression cases.
• Then, it performs survival analysis based on overall survival using only these high and low expression cases.
• The plot command tells it whether to output a survival plot, and the data command whether to output the text of the results.
• The variable assignment here to ‘testSurvival’ returns the coxZscore for that gene, allowing you to get results for many genes to compare.

If we want to do a list of a bunch of genes, we can do that as below. Keep in mind that if you feed it 1,000 genes, it will make 1,000 plots and text files if you tell it to, so make sure you write to a directory where you can easily make sense of this if you need to (perhaps a new folder holding just these output files).

##########################################################################################
geneOfInterest = as.list(c('YBX1','EIF4A1','PODXL'))
multiSurvival = lapply(geneOfInterest, survivalAnalysis, 
                       expressionData = esExpression, 
                       phenotypeData = esPhenotype, 
                       survivalPlot = TRUE, 
                       survivalData = TRUE, 
                       writeDirectory = '/getGeoDataForSurvivalAnalysis')

[1] "Performing analysis for YBX1"
[1] "Performing analysis for EIF4A1"
[1] "Performing analysis for PODXL"

This is the end of this document. As a final command, we will run sessionInfo to ensure others who look at our work know exactly how our environment was set up.

sessionInfo()

R version 3.6.3 (2020-02-29)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18363)

Matrix products: default

locale:
[1] LC_COLLATE=English_Canada.1252  LC_CTYPE=English_Canada.1252    LC_MONETARY=English_Canada.1252 LC_NUMERIC=C                   
[5] LC_TIME=English_Canada.1252    

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] survminer_0.4.7     ggpubr_0.3.0        survival_3.1-12     RColorBrewer_1.1-2  GEOquery_2.54.1     Biobase_2.46.0      BiocGenerics_0.32.0
 [8] forcats_0.5.0       stringr_1.4.0       dplyr_0.8.5         purrr_0.3.4         readr_1.3.1         tidyr_1.0.3         tibble_3.0.1       
[15] ggplot2_3.3.1       tidyverse_1.3.0    

loaded via a namespace (and not attached):
 [1] httr_1.4.1        jsonlite_1.6.1    splines_3.6.3     carData_3.0-4     modelr_0.1.8      assertthat_0.2.1  blob_1.2.1        cellranger_1.1.0 
 [9] yaml_2.2.1        pillar_1.4.4      backports_1.1.7   lattice_0.20-38   glue_1.4.1        limma_3.42.2      digest_0.6.25     ggsignif_0.6.0   
[17] rvest_0.3.5       colorspace_1.4-1  Matrix_1.2-18     pkgconfig_2.0.3   broom_0.5.6       haven_2.2.0       xtable_1.8-4      scales_1.1.1     
[25] km.ci_0.5-2       openxlsx_4.1.5    rio_0.5.16        KMsurv_0.1-5      farver_2.0.3      generics_0.0.2    car_3.0-8         ellipsis_0.3.1   
[33] withr_2.2.0       cli_2.0.2         magrittr_1.5      crayon_1.3.4      readxl_1.3.1      fs_1.4.1          fansi_0.4.1       nlme_3.1-144     
[41] rstatix_0.5.0     xml2_1.3.2        foreign_0.8-75    tools_3.6.3       data.table_1.12.8 hms_0.5.3         lifecycle_0.2.0   munsell_0.5.0    
[49] reprex_0.3.0      zip_2.0.4         compiler_3.6.3    rlang_0.4.6       grid_3.6.3        rstudioapi_0.11   labeling_0.3      gtable_0.3.0     
[57] abind_1.4-5       DBI_1.1.0         curl_4.3          R6_2.4.1          gridExtra_2.3     zoo_1.8-8         lubridate_1.7.8   knitr_1.28       
[65] survMisc_0.5.5    stringi_1.4.6     Rcpp_1.0.4.6      vctrs_0.3.0       dbplyr_1.4.4      tidyselect_1.1.0  xfun_0.14