Shixiang Wang

>上士闻道
勤而行之

UCSCXenaTools: Retrieve Gene Expression and Clinical Information from UCSC Xena for Survival Analysis

Shixiang Wang · 2019-08-25

Categories: bioinformatics  
Tags: r   software   package   bioinformatics   data-access   survival-analysis   UCSCXenaTools  

I thank the edition made by Stefanie Butland.

The UCSC Xena platform provides an unprecedented resource for public omics data from big projects like The Cancer Genome Atlas (TCGA), however, it is hard for users to incorporate multiple datasets or data types, integrate the selected data with popular analysis tools or homebrewed code, and reproduce analysis procedures. To address this issue, we developed an R package UCSCXenaTools for enabling data retrieval, analysis integration and reproducible research for omics data from the UCSC Xena platform1.

In this technote we will outline how to use the UCSCXenaTools package to pull gene expression and clinical data from UCSC Xena for survival analysis.

For general usage of UCSCXenaTools, please refer to the package vignette. Any bug or feature request can be reported in GitHub issues.

Installation

UCSCXenaTools is available from CRAN:

install.packages("UCSCXenaTools")

Alternatively, the latest development version can be downloaded from GitHub:

remotes::install_github("ropensci/UCSCXenaTools", build_vignettes = TRUE, dependencies = TRUE)

How it works

Before actually pulling data, understanding how UCSCXenaTools works (see Figure 1) will help users locate the most important function to use.

Generally,

Figure 1. The UCSCXenaTools pipeline

We will provide an example illustrating how to use UCSCXenaTools to study the effect of expression of the KRAS gene on prognosis of Lung Adenocarcinoma (LUAD) patients. KRAS is a known driver gene in LUAD. We retrieve expression data for the KRAS gene and survival status data for LUAD patients from the TCGA and use these as input to a survival analysis, frequently used in cancer research.

Download data

First we get information on all datasets in the TCGA LUAD cohort and store as luad_cohort object.

suppressMessages(library(UCSCXenaTools))
suppressMessages(library(dplyr))
luad_cohort = XenaData %>% 
  filter(XenaHostNames == "tcgaHub") %>% # select TCGA Hub
  XenaScan("TCGA Lung Adenocarcinoma")   # select LUAD cohort

luad_cohort
#> # A tibble: 23 x 17
#>    XenaHosts XenaHostNames XenaCohorts XenaDatasets SampleCount DataSubtype
#>    <chr>     <chr>         <chr>       <chr>              <int> <chr>      
#>  1 https://… tcgaHub       TCGA Lung … RABIT/separ…         467 pathway ac…
#>  2 https://… tcgaHub       TCGA Lung … RABIT/separ…         120 pathway ac…
#>  3 https://… tcgaHub       TCGA Lung … TCGA.LUAD.s…         151 DNA methyl…
#>  4 https://… tcgaHub       TCGA Lung … TCGA.LUAD.s…         492 DNA methyl…
#>  5 https://… tcgaHub       TCGA Lung … TCGA.LUAD.s…         516 copy numbe…
#>  6 https://… tcgaHub       TCGA Lung … TCGA.LUAD.s…         237 protein ex…
#>  7 https://… tcgaHub       TCGA Lung … TCGA.LUAD.s…         576 gene expre…
#>  8 https://… tcgaHub       TCGA Lung … TCGA.LUAD.s…          60 miRNA matu…
#>  9 https://… tcgaHub       TCGA Lung … TCGA.LUAD.s…         576 gene expre…
#> 10 https://… tcgaHub       TCGA Lung … TCGA.LUAD.s…         576 gene expre…
#> # … with 13 more rows, and 11 more variables: Label <chr>, Type <chr>,
#> #   AnatomicalOrigin <chr>, SampleType <chr>, Tags <chr>, ProbeMap <chr>,
#> #   LongTitle <chr>, Citation <chr>, Version <chr>, Unit <chr>, Platform <chr>

Download clinical dataset

Now we download clinical datasets of TCGA LUAD cohort and load them into R.

cli_query = luad_cohort %>% 
  filter(DataSubtype == "phenotype") %>%  # select clinical dataset
  XenaGenerate() %>%  # generate a XenaHub object
  XenaQuery() %>% 
  XenaDownload()
#> This will check url status, please be patient.
#> All downloaded files will under directory /var/folders/bj/nw1w4g1j37ddpgb6zmh3sfh80000gn/T//RtmpLpcBQu.
#> The 'trans_slash' option is FALSE, keep same directory structure as Xena.
#> Creating directories for datasets...
#> Downloading TCGA.LUAD.sampleMap/LUAD_clinicalMatrix
#> Downloading survival/LUAD_survival.txt.gz

cli = XenaPrepare(cli_query)

NOTE: from 2019-12-06, UCSC Xena stores phenotype data in two datasets, one is survival data, the other is other phenotype data. So cli is a list here.

Keep survival data here, we only need it.

cli = cli$LUAD_survival.txt.gz

# See a few rows
head(cli)
#> # A tibble: 6 x 11
#>   sample `_PATIENT`    OS OS.time   DSS DSS.time   DFI DFI.time   PFI PFI.time
#>   <chr>  <chr>      <dbl>   <dbl> <dbl>    <dbl> <dbl>    <dbl> <dbl>    <dbl>
#> 1 TCGA-… TCGA-05-4…     0       0     0        0    NA       NA     0        0
#> 2 TCGA-… TCGA-05-4…     0    1523     0     1523    NA       NA     0     1523
#> 3 TCGA-… TCGA-05-4…     1     121    NA      121    NA       NA     0      121
#> 4 TCGA-… TCGA-05-4…     0     607     0      607     1      334     1      334
#> 5 TCGA-… TCGA-05-4…     0     426     0      426    NA       NA     1      183
#> 6 TCGA-… TCGA-05-4…     0    1369     0     1369    NA       NA     0     1369
#> # … with 1 more variable: Redaction <lgl>

Download KRAS gene expression

To download gene expression data, first we need to select the right dataset.

ge = luad_cohort %>% 
  filter(DataSubtype == "gene expression RNAseq", Label == "IlluminaHiSeq")
ge
#> # A tibble: 1 x 17
#>   XenaHosts XenaHostNames XenaCohorts XenaDatasets SampleCount DataSubtype Label
#>   <chr>     <chr>         <chr>       <chr>              <int> <chr>       <chr>
#> 1 https://… tcgaHub       TCGA Lung … TCGA.LUAD.s…         576 gene expre… Illu…
#> # … with 10 more variables: Type <chr>, AnatomicalOrigin <chr>,
#> #   SampleType <chr>, Tags <chr>, ProbeMap <chr>, LongTitle <chr>,
#> #   Citation <chr>, Version <chr>, Unit <chr>, Platform <chr>

Now we fetch KRAS gene expression values.

# You can pass gene symbols to 'identifiers' option
# to obtain their values in a dataset.
# A matrix will be returned by 'fetch_dense_values' function
# with rows corresponding to genes,
# so here we extract the first row.
KRAS = fetch_dense_values(host = ge$XenaHosts, 
                          dataset = ge$XenaDatasets, 
                          identifiers = "KRAS",
                          use_probeMap = TRUE) %>% 
  .[1, ]
#> -> Checking identifiers...
#> -> use_probeMap is TRUE, skipping checking identifiers...
#> -> Done.
#> -> Checking samples...
#> -> Done.
#> -> Checking if the dataset has probeMap...
#> -> Done. ProbeMap is found.

head(KRAS)
#> TCGA-69-7978-01 TCGA-62-8399-01 TCGA-78-7539-01 TCGA-50-5931-11 TCGA-73-4658-01 
#>           10.25           10.29           10.82           10.29           10.36 
#> TCGA-44-6775-01 
#>           10.03

Merge expression data and survival status

Next, we join the two data.frame by sample and keep necessary columns. Here we focus on ‘Primary Tumor’ for simplicity.

merged_data = tibble(sample = names(KRAS),
                     KRAS_expression = as.numeric(KRAS)) %>% 
  left_join(cli, by = "sample") %>% 
  filter(substr(sample, 14, 15) == "01") %>%  # Keep only 'Primary Tumor'
  select(sample, KRAS_expression, OS.time, OS) %>% 
  rename(time = OS.time, 
         status = OS)
  

head(merged_data)
#> # A tibble: 6 x 4
#>   sample          KRAS_expression  time status
#>   <chr>                     <dbl> <dbl>  <dbl>
#> 1 TCGA-69-7978-01           10.2    134      0
#> 2 TCGA-62-8399-01           10.3   2696      0
#> 3 TCGA-78-7539-01           10.8    791      0
#> 4 TCGA-73-4658-01           10.4   1600      1
#> 5 TCGA-44-6775-01           10.0    705      0
#> 6 TCGA-44-2655-01            9.75  1324      0

Survival analysis

To study the effect of KRAS gene expression on prognosis of LUAD patients, we show two approaches:

  1. use Cox model to determine the effect when KRAS gene expression increases
  2. use Kaplan-Meier curve and log-rank test to observe the difference in different ofKRAS gene expression status, i.e. high or low

We will use package survival and survminer to create models and plot survival curves, respectively.

library(survival)
library(survminer)
#> Loading required package: ggplot2
#> Loading required package: ggpubr

Cox model

fit = coxph(Surv(time, status) ~ KRAS_expression, data = merged_data)
fit
#> Call:
#> coxph(formula = Surv(time, status) ~ KRAS_expression, data = merged_data)
#> 
#>                   coef exp(coef) se(coef)     z       p
#> KRAS_expression 0.2922    1.3394   0.1020 2.864 0.00418
#> 
#> Likelihood ratio test=7.64  on 1 df, p=0.0057
#> n= 506, number of events= 183 
#>    (9 observations deleted due to missingness)

We can find that patients with higher KRAS gene expression have higher risk (34% increase per KRAS gene expression unit increase), and the effect of KRAS gene expression is statistically significant (p<0.05).

If you know little about survival analysis, two blogs are recommended to read:

  1. Survival Analysis Basics
  2. Cox Proportional-Hazards Model

Risk between expression groups

We can also divide patients into two groups using KRAS median as a cutoff.

merged_data = merged_data %>% 
  mutate(group = case_when(
    KRAS_expression > quantile(KRAS_expression, 0.5) ~ 'KRAS_High',
    KRAS_expression < quantile(KRAS_expression, 0.5) ~ 'KRAS_Low',
    TRUE ~ NA_character_
  ))

fit = survfit(Surv(time, status) ~ group, data = merged_data)

Then we can plot the survival curves for each group.

ggsurvplot(fit, pval = TRUE)
Figure 2. Kaplan-Meier curve. Survival probability vs Time (days)

Figure 1: Figure 2. Kaplan-Meier curve. Survival probability vs Time (days)

The Kaplan-Meier plot shows what percent of patients are alive at a time point. We can clearly see that patients in ‘KRAS_Low’ group have better survival than patients in ‘KRAS_High’ group because the survival probability of ‘KRAS_High’ group is always lower than ‘KRAS_Low’ group over time (the unit is ‘day’ here). The difference between the two groups is statistically significant (p<0.05 by log-rank test).

Acknowledgements

We thank Christine Stawitz and Carl Ganz for their constructive comments. This package is reviewed by rOpenSci at https://github.com/ropensci/software-review/issues/315.


  1. Wang et al., (2019). The UCSCXenaTools R package: a toolkit for accessing genomics data from UCSC Xena platform, from cancer multi-omics to single-cell RNA-seq. Journal of Open Source Software, 4(40), 1627, https://doi.org/10.21105/joss.01627↩︎