Connection to the cBioPortal Server

Before starting a TCGAretriever pipeline, it’s recommended to check whether the remote data server can be reached. Thos is achieved via the check_TCGAretriever_connection() function. This should return TRUE.

library(TCGAretriever)
library(reshape2)
library(ggplot2)

# Check connection
con_alive <- check_TCGAretriever_connection()
con_alive

## [1] TRUE

Cancer Study IDs

Here we retrieve the list of available studies and select all tcga_pub entries. Each of the values included in the studyId column is a study identifier and can be passed as csid argument to other functions such as get_genetic_profiles() or get_case_lists().

# If connection is alive, obtain a list of cancer studies from cBio
if (con_alive) {
  
  all_studies <- get_cancer_studies()

  # Find published TCGA datasets
  keep <- grepl("tcga_pub$", all_studies[,'studyId'])
  tcga_studies <- all_studies[keep, ]

  # Show results
  utils::head(tcga_studies[, c(11, 1, 4)])
  
} else {
  
  print('The server cannot be reached. All code below won\'t be executed.')

}

##     readPermission         studyId
## 33            TRUE   blca_tcga_pub
## 62            TRUE   brca_tcga_pub
## 135           TRUE   hnsc_tcga_pub
## 153           TRUE   laml_tcga_pub
## 155           TRUE   kich_tcga_pub
## 163           TRUE lgggbm_tcga_pub
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          description
## 33                                                                                                                                                                                                                                                                                                                                                                                                            Whole-exome sequencing of 131 high-grade muscle-invasive urothelial bladder carcinomas
## 62  Whole-exome sequencing (510 samples with matched normals), genomic DNA copy number arrays, DNA methylation, messenger RNA arrays, microRNA sequencing and reverse-phase protein arrays analysis on 825 primary breast cancer samples. <a href="http://cancergenome.nih.gov/">The Cancer Genome Atlas (TCGA)</a> Breast Invasive Carcinoma project. 825 cases.<br><i>Nature 2012.</i> <a href="https://tcga-data.nci.nih.gov/docs/publications/brca_2012/">Raw data via the TCGA Data Portal</a>.
## 135                                                                                                                                                                                                                                                                                             Whole-exome sequencing and/or whole-genome sequencing of 279 head and neck squamous cell carcinoma tumor/normal pairs. The Cancer Genome Atlas (TCGA) head and neck squamous cell carcinoma Project.
## 153                                                                                                                                                                                                            Whole-genome or whole-exome sequencing analysis of 200 adult de novo AML tumor/normal pairs. TCGA Acute Myeloid Leukemia Project; raw data at the <A HREF="https://tcga-data.nci.nih.gov/tcga/tcgaCancerDetails.jsp?diseaseType=LAML&diseaseName=Acute%20Myeloid%20Leukemia">NCI</A>.
## 155                                                                                                                                                                                                                                                                                                                                                                                                       Whole-exome sequencing of 66 chromophobe renal cell carcinoma (ChRCCs) tumor/normal pairs.
## 163                                                                                                                                                                                                                                                                                                                                                                                                                                  Whole exome sequencing of 1,122 LGG and GBM tumor/normal pairs.

Genetic Profiles (Assays) and Case Lists

It is possible to retrieve genetic profiles and case lists associated to each study of interest. Genetic profiles indicate what kind of data are available for a certain study (e.g., RNA-seq, micro-array, DNA mutation…). Typically, a data type of interest may only be available for a fraction of patients in the cohort. Therefore, it is also important to obtain the identifier of the corresponding case list of interest.

# Define the cancer study id: brca_tcga_pub
my_csid <- "brca_tcga_pub"

# Obtain genetic profiles
brca_pro <- get_genetic_profiles(csid = my_csid)
utils::head(brca_pro[, c(7, 1)])

##   showProfileInAnalysisTab                 molecularProfileId
## 1                     TRUE               brca_tcga_pub_gistic
## 2                    FALSE           brca_tcga_pub_linear_CNA
## 3                    FALSE     brca_tcga_pub_methylation_hm27
## 4                    FALSE                brca_tcga_pub_mirna
## 5                    FALSE brca_tcga_pub_mirna_median_Zscores
## 6                    FALSE                 brca_tcga_pub_mrna

# Obtain cases 
brca_cas <- get_case_lists(csid = my_csid)
utils::head(brca_cas[, c(4, 1)])

##                                 name           sampleListId
## 1                        All samples      brca_tcga_pub_all
## 2                        PAM50 Basal    brca_tcga_pub_basal
## 3                  PAM50 Claudin low  brca_tcga_pub_claudin
## 4              Samples with CNA data      brca_tcga_pub_cna
## 5 Samples with mutation and CNA data   brca_tcga_pub_cnaseq
## 6                   Complete samples brca_tcga_pub_complete

Retrieve Genomic Data

Once we identified a genetic profile of interest and a case list of interest, we can obtain clinical data via the TCGAretriever::get_clinical_data(). Additionally, molecular data can be obtained via the TCGAretriever::get_molecular_data() function (non-mutation data) or the get_mutation_data() function (mutation data).

# Define a set of genes of interest
q_csid <- 'brca_tcga_pub'
q_genes <- c("TP53", "HMGA1", "E2F1", "EZH2")
q_cases <- "brca_tcga_pub_complete"
rna_prf <- "brca_tcga_pub_mrna"
mut_prf <- "brca_tcga_pub_mutations"

# Download Clinical Data
brca_cli <- get_clinical_data(csid = q_csid, case_list_id = q_cases)

# Download RNA
brca_RNA <- get_molecular_data(case_list_id = q_cases, 
                               gprofile_id = rna_prf, 
                               glist = q_genes)

• NOTE 1: the resulting data.frame includes ENTREZ_GENE_IDs and OFFICIAL_SMBOLs as first and second column. The third column indicates the gene/assay type.

• NOTE 2: up to n=500 genes can be queried via the get_molecular_data() function.

# Show results
brca_RNA[, 1:5]

##   entrezGeneId hugoGeneSymbol           type TCGA-A1-A0SD-01 TCGA-A1-A0SE-01
## 1         1869           E2F1 protein-coding         -1.4535         -1.3525
## 2         2146           EZH2 protein-coding         -1.7928         -1.5652
## 3         3159          HMGA1 protein-coding         -2.5195         -1.9663
## 4         7157           TP53 protein-coding         -1.0107         -0.3385

# Set SYMBOLs as rownames
# Note that you may prefer to use the tibble package for this
rownames(brca_RNA) <- brca_RNA$hugoGeneSymbol
brca_RNA <- brca_RNA[, -c(1, 2, 3)]

# Round numeric vals to 3 decimals
for (i in 1:ncol(brca_RNA)) {
  brca_RNA[, i] <- round(brca_RNA[, i], digits = 3)
}

# Download mutations
brca_MUT <- get_mutation_data(case_list_id = q_cases, 
                              gprofile_id = mut_prf, 
                              glist = q_genes)

# Identify Samples carrying a TP53 missense mutation
tp53_mis_keep <- brca_MUT$hugoGeneSymbol == 'TP53' &
  brca_MUT$mutationType == 'Missense_Mutation' &
  !is.na(brca_MUT$sampleId)
tp53_mut_samples <- unique(brca_MUT$sampleId[tp53_mis_keep])

# Show results
keep_cols <- c('sampleId', 'hugoGeneSymbol', 'mutationType',  'proteinChange')
utils:::head(brca_MUT[, keep_cols])

##          sampleId hugoGeneSymbol      mutationType proteinChange
## 1 TCGA-A1-A0SK-01           TP53 Missense_Mutation         M133K
## 2 TCGA-A1-A0SO-01           TP53 Missense_Mutation         Y220C
## 3 TCGA-A2-A04P-01           TP53   Frame_Shift_Del    G108Vfs*15
## 4 TCGA-A2-A04R-01           E2F1 Missense_Mutation         Y168C
## 5 TCGA-A2-A04W-01           TP53 Missense_Mutation         V272M
## 6 TCGA-A2-A0CL-01           TP53 Missense_Mutation         R273L

NOTE: Mutation data are returned in a different format (molten data format) compared to data corresponding to other genetic profiles. Results are formatted as molten data.frames where each row is a DNA variant. Each sample / case may have 0 or more variants for each of the query genes.

Retrieve Data for all available genes

It is possible to download data corresponding to all available genes using the fetch_all_tcgadata() function. This function takes the same arguments as the get_profile_data() function but the glist argument. Gene lists are obtained automatically via the get_gene_identifiers() function. Each job is automatically split into multiple batches. Results are automatically aggregated before being returned. If mutation data are being queried, results are returned as a molten data.frame. The mutations argument is used to guide the query type as well as the output format.

# Download all brca_pub mutation data (complete samples)
all_brca_MUT <- fetch_all_tcgadata(case_list_id = q_cases, 
                                   gprofile_id = mut_prf, 
                                   mutations = TRUE)
  
# Download all brca_pub RNA expression data (complete samples)
all_brca_RNA <- fetch_all_tcgadata(case_list_id = q_cases, 
                                   gprofile_id = rna_prf, 
                                   mutations = FALSE)

Relationship between E2F1 and EZH2 in BRCA

Here we use the data retrieved before to analyze the relationship between E2F1 and EZH2 expression in TCGA breast cancer samples.

# Visualize the correlation between EZH2 and E2F1
df <- data.frame(sample_id = colnames(brca_RNA), 
                 EZH2 = as.numeric(brca_RNA['EZH2', ]), 
                 E2F1 = as.numeric(brca_RNA['E2F1', ]), 
                 stringsAsFactors = FALSE)

ggplot(df, aes(x = EZH2, y = E2F1)) +
  geom_point(color = 'gray60', size = 0.75) +
  theme_bw() +
  geom_smooth(method = 'lm', color = 'red2', 
              size=0.3, fill = 'gray85') +
  ggtitle('E2F1-EZH2 correlation in BRCA') + 
  theme(plot.title = element_text(hjust = 0.5))

Scatterplot - showing E2F1 RNA expression with respect to EZH2 RNA expression across TCGA breast cancer samples.

E2F1 Expression vs. Binned EZH2 Expression in BRCA

This is an alternative visualization that illustrates the same relationship between E2F1 and EZH2 expression. Here, we make use of the make_groups() function to discretize the numeric vector corresponding to EZH2 expression counts. This function takes a numeric vector as its num_vector argument (e.g., a gene expression vector), performs binning, and returns the results as a data.frame. Here we show the distribution of E2F1 expression counts with respect to binned EZH2 expression.

# Bin samples according to EZH2 expression
EZH2_bins <- make_groups(num_vector = df$EZH2, groups = 5) 
utils::head(EZH2_bins, 12)

##     value rank
## 1  -1.793    2
## 2  -1.565    2
## 3  -2.495    1
## 4  -2.136    1
## 5   0.529    5
## 6  -1.439    3
## 7  -0.426    5
## 8  -1.165    4
## 9  -0.894    4
## 10 -0.998    4
## 11 -1.132    4
## 12 -1.483    3

# attach bin to df
df$EZH2_bin <- EZH2_bins$rank

# build Boxplot
ggplot(df, aes(x = as.factor(EZH2_bin), y = E2F1)) +
  geom_boxplot(outlier.shape = NA, fill = '#fed976') +
  geom_jitter(width = 0.1, size = 1) +
  theme_bw() +
  xlab('EZH2 Bins') +
  ggtitle('E2F1 Expression vs. Binned EZH2 Expression') +
  theme(plot.title = element_text(face = 'bold', hjust = 0.5))

Boxplot - showing E2F1 RNA expression with respect to binned EZH2 RNA expression across TCGA breast cancer samples.

Relationship between HMGA1 and TP53 with respect to TP53 mutation status in BRCA

We use the data retrieved before to analyze the relationship between HMGA1 and TP53 expression with respect to the TP53 WT status in TCGA breast cancer samples.

# Coerce to data.frame with numeric features 
mol_df <- data.frame(sample_id = colnames(brca_RNA), 
                     HMGA1 = as.numeric(brca_RNA['HMGA1', ]),
                     TP53 = as.numeric(brca_RNA['TP53', ]),
                     stringsAsFactors = FALSE)

mol_df$TP53.status = factor(ifelse(mol_df$sample_id %in% tp53_mut_samples, 
                                   '01.wild_type', '02.mutated'))
  
# Visualize the correlation between EZH2 and E2F1
ggplot(mol_df, aes(x = TP53, y = HMGA1)) +
  geom_point(color = 'gray60', size = 0.75) +
  facet_grid(cols = vars(TP53.status)) +
  theme_bw() +
  geom_smooth(mapping = aes(color = TP53.status), 
              method = 'lm', size=0.3, fill = 'gray85') +
  ggtitle('HMGA1-TP53 correlation in BRCA') + 
  theme(plot.title = element_text(hjust = 0.5))

Scatterplot - showing HGMA1 RNA expression with respect to TP53 RNA expression across TCGA breast cancer samples. Samples with wild type (left) and mutant (right) TP53 status are shown in the two panels.