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

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 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).

Related project

XenaShiny, a Shiny project based on UCSCXenaTools, is under development by my friends and me.

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.