1 Intro

This exercise will show how to obtain clinical and genomic data from the Cancer Genome Atlas (TGCA) and to perform classical analysis important for clinical data.

These include:

  1. Download the data (clinical and expresion) from TGCA
  2. Processing of the data (normalization) and saving it locally using simple table formats.
  3. Unsupervised analysis includes differential expression, PCA and clustering.
  4. Build a machine learning model (classifier) to predict cancer.
  5. Perform survival analysis of molecular markers detected in previous analysis.

First, we start by loading all libraries necessary for this exercise. Please check their documentation if you want to know more.

# Load packages
library("TCGAbiolinks")
library("limma")
library("edgeR")
library("glmnet")
library("factoextra")
library("FactoMineR")
library("caret")
library("SummarizedExperiment")
library("gplots")
library("survival")
library("survminer")
library("RColorBrewer")
library("gProfileR")
library("genefilter")

2 TCGA data

In this tutorial, we will focus on Liver Hepatocellular Carcinoma, which is identified in TCGA as LIHC. For LIHC, TCGA provides data for 377 patients including: clinical, expression, DNA methylation and genotyping data. In this tutorial, we will work with clinical and expression data (RNA-seq). Go to https://portal.gdc.cancer.gov/ and search for TCGA-LIHC if you want to understand the data deposited in TCGA. You can also try to find your way through the previous data to look for other data sets of your interest.

We will make use of the TCGAbiolinks library, which allows us to query, prepare and download data from the TCGA portal. TCGAbiolinks provides important functionality as matching data of same the donors across distinct data types (clinical vs expression) and provides data structures to make its analysis in R easy.

To download TCGA data with TCGAbiolinks, you need to follow 3 steps. First, you will query the TCGA database through R with the function GDCquery. This will allow you to investigate the data available at the TCGA database. Next, we use GDCdownload to download raw version of desired files into your computer. Finally GDCprepare will read these files and make R data structures so that we can further analyse them.

Before we get there however we need to know what we are searching for. We can check all the available projects at TCGA with the command bellow. Since there are many lets look at the first 6 projects using the command head().

GDCprojects = getGDCprojects()

head(GDCprojects[c("project_id", "name")])

As a general rule in R (and especially if you are working in RStudio) whenever some method returns some value or table you are not familiar with, you should check its structure and dimensions. You can always use functions such as head() to only show the first entries and dim() to check the dimension of the data.

We already know that Liver Hepatocellular Carcinoma has as id TCGA-LIHC. We can use the following function to get details on all data deposited for TCGA-LIHC.

TCGAbiolinks:::getProjectSummary("TCGA-LIHC")
## $case_count
## [1] 377
## 
## $file_count
## [1] 11578
## 
## $file_size
## [1] 1.780846e+13
## 
## $data_categories
##   case_count file_count               data_category
## 1        376       2122     Transcriptome Profiling
## 2        375       3032 Simple Nucleotide Variation
## 3        377       1634                 Biospecimen
## 4        377        423                    Clinical
## 5        377        430             DNA Methylation
## 6        376       2283       Copy Number Variation
## 7        377       1654            Sequencing Reads

Of note, not all patients were measured for all data types. Also, some data types have more files than samples. This is the case when more experiments were performed per patient, i.e. transcriptome profiling was performed both in mRNA and miRNA, or that data have been analysed by distinct computational strategies.

Let us start by querying all RNA-seq data from LIHC project. When using GDCquery we always need to specify the id of the project, i.e. “TCGA_LIHC”, and the data category we are interested in, i.e. “Transcriptome Profiling”. Here, we will focus on a particular type of data summarization for mRNA-seq data (workflow.type), which is based on raw counts estimated with HTSeq.

Note that performing this query will take a few of minutes

query_TCGA = GDCquery(
  project = "TCGA-LIHC",
  data.category = "Transcriptome Profiling", # parameter enforced by GDCquery
  experimental.strategy = "RNA-Seq",
  workflow.type = "HTSeq - Counts")

To visualize the query results in a more readable way, we can use the command getResults.

lihc_res = getResults(query_TCGA) # make results as table
# head(lihc_res) # data of the first 6 patients.
colnames(lihc_res) # columns present in the table
##  [1] "file_id"                   "md5sum"                   
##  [3] "file_name"                 "data_type"                
##  [5] "data_category"             "data_format"              
##  [7] "created_datetime"          "id"                       
##  [9] "updated_datetime"          "data_release"             
## [11] "submitter_id"              "version"                  
## [13] "experimental_strategy"     "access"                   
## [15] "cases"                     "state"                    
## [17] "type"                      "file_size"                
## [19] "project"                   "analysis_id"              
## [21] "analysis_state"            "analysis_submitter_id"    
## [23] "analysis_workflow_link"    "analysis_workflow_type"   
## [25] "analysis_workflow_version" "sample_type"              
## [27] "is_ffpe"                   "cases.submitter_id"       
## [29] "sample.submitter_id"

One interesting question is the tissue type measured at an experiment (normal, solid tissue, cell line). This information is present at column “sample_type”.

head(lihc_res$sample_type) # first 6 types of tissue.
## [1] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [4] "Solid Tissue Normal" "Primary Tumor"       "Primary Tumor"

We can visualize it better with the summary function.

summary(factor(lihc_res$sample_type)) # summary of distinct tissues types present in this study
##       Primary Tumor     Recurrent Tumor Solid Tissue Normal 
##                 371                   3                  50

As you see, there are 50 controls (Solid Tissue Normal) and 371 cancer samples (Primary Solid Tumors). For simplicity, we will ignore the small class of recurrent solid tumors. Therefore, we will redo the query as

query_TCGA = GDCquery(
  project = "TCGA-LIHC",
  data.category = "Transcriptome Profiling", # parameter enforced by GDCquery
  experimental.strategy = "RNA-Seq",
  workflow.type = "HTSeq - Counts",
  sample.type = c("Primary Tumor", "Solid Tissue Normal"))

Next, we need to download the files from the query. Before, be sure that you set your current working directory to the place you want to save your data. TCGA will save the data in a directory structure startign with a directory “GDCdata”.

Let us now download the files specified in the query.

GDCdownload(query = query_TCGA)

Given that you need to download many files, the previous operation might take some time. Often the download fails for one or another file. You can re-run the previous command until no error message is given. The method will recognize that the data has already been downloaded and won’t download the data again.

Finally, lets load the actual RNASeq data into R. Remember that the output directory set must be the same to where you downloaded the data.

tcga_data = GDCprepare(query_TCGA)

We can then check the size of the object with the command.

dim(tcga_data)
## [1] 56602   421

There are 3 functions that allow us to access to most important data present in this object, these are: colData(), rowData(), assays(). Use the command ?SummarizedExperiment to find more details. colData() allows us to access the clinical data associated with our samples. The functions colnames() and rownames() can be used to extract the column and rows names from a given table respectively.

# In R (and other programming languages) you can often
# chain functions to save time and space
colnames(colData(tcga_data))
##  [1] "barcode"                     "patient"                    
##  [3] "sample"                      "shortLetterCode"            
##  [5] "definition"                  "sample_submitter_id"        
##  [7] "sample_type_id"              "days_to_collection"         
##  [9] "sample_type"                 "state"                      
## [11] "pathology_report_uuid"       "is_ffpe"                    
## [13] "tissue_type"                 "initial_weight"             
## [15] "submitter_id"                "sample_id"                  
## [17] "oct_embedded"                "tumor_grade"                
## [19] "ajcc_pathologic_stage"       "primary_diagnosis"          
## [21] "progression_or_recurrence"   "age_at_diagnosis"           
## [23] "year_of_diagnosis"           "ajcc_pathologic_m"          
## [25] "treatments"                  "ajcc_staging_system_edition"
## [27] "tissue_or_organ_of_origin"   "diagnosis_id"               
## [29] "prior_treatment"             "site_of_resection_or_biopsy"
## [31] "days_to_diagnosis"           "morphology"                 
## [33] "tumor_stage"                 "classification_of_tumor"    
## [35] "days_to_last_follow_up"      "ajcc_pathologic_n"          
## [37] "icd_10_code"                 "last_known_disease_status"  
## [39] "ajcc_pathologic_t"           "prior_malignancy"           
## [41] "synchronous_malignancy"      "height"                     
## [43] "alcohol_history"             "exposure_id"                
## [45] "weight"                      "bmi"                        
## [47] "year_of_death"               "gender"                     
## [49] "demographic_id"              "days_to_birth"              
## [51] "race"                        "age_at_index"               
## [53] "ethnicity"                   "year_of_birth"              
## [55] "vital_status"                "days_to_death"              
## [57] "bcr_patient_barcode"         "name"                       
## [59] "disease_type"                "project_id"                 
## [61] "released"                    "primary_site"               
## [63] "releasable"                  "sample.aux"

This link provides a basic explanation about tcga_data. Note that both clinical and expression data are present in this object.

Lets look at some potentially interesting features. The table() function (in this context) produces a small summary with the sum of each of the factors present in a given column.

table(tcga_data@colData$vital_status)
## 
##        Alive         Dead Not Reported 
##          255          164            2
table(tcga_data@colData$tumor_stage)
## 
## not reported      stage i     stage ii    stage iii   stage iiia   stage iiib 
##           32          189           97            6           73            8 
##   stage iiic     stage iv    stage iva    stage ivb 
##           10            3            1            2
table(tcga_data@colData$definition)
## 
## Primary solid Tumor Solid Tissue Normal 
##                 371                  50
table(tcga_data@colData$tissue_or_organ_of_origin)
## 
## Liver 
##   421
table(tcga_data@colData$gender)
## 
## female   male 
##    143    278
table(tcga_data@colData$race)
## 
## american indian or alaska native                            asian 
##                                2                              164 
##        black or african american                     not reported 
##                               24                               13 
##                            white 
##                              218

Is there a particular column (feature) that allows you to distinguish tumor tissue from normal tissue?

What about the RNA-seq data? We can use the assay function to obtain the RNA-seq count matrices and rowData to see gene mapping information. Can you tell how many genes and how many samples are included there?

dim(assay(tcga_data))     # gene expression matrices.
## [1] 56602   421
head(assay(tcga_data)[,1:10]) # expression of first 6 genes and first 10 samples
##                 TCGA-DD-AAE7-01A-11R-A41C-07 TCGA-DD-AADJ-01A-11R-A41C-07
## ENSG00000000003                         2422                         3658
## ENSG00000000005                            0                            1
## ENSG00000000419                          380                          571
## ENSG00000000457                          302                          365
## ENSG00000000460                           84                           68
## ENSG00000000938                          147                          124
##                 TCGA-UB-A7MD-01A-12R-A352-07 TCGA-BD-A2L6-11A-21R-A213-07
## ENSG00000000003                         2151                         4219
## ENSG00000000005                            4                            4
## ENSG00000000419                         1268                         1070
## ENSG00000000457                          553                          565
## ENSG00000000460                          861                          133
## ENSG00000000938                          356                         1493
##                 TCGA-CC-5259-01A-31R-A213-07 TCGA-HP-A5MZ-01A-21R-A27V-07
## ENSG00000000003                         3324                        10319
## ENSG00000000005                            1                            2
## ENSG00000000419                          341                         1137
## ENSG00000000457                          170                          687
## ENSG00000000460                           43                          226
## ENSG00000000938                           32                          382
##                 TCGA-DD-AAVS-01A-11R-A41C-07 TCGA-CC-A7IH-01A-11R-A33J-07
## ENSG00000000003                         8408                         9419
## ENSG00000000005                            0                            0
## ENSG00000000419                         1502                         1404
## ENSG00000000457                         1196                         1650
## ENSG00000000460                          436                         1585
## ENSG00000000938                          277                           81
##                 TCGA-FV-A3I0-01A-11R-A22L-07 TCGA-MI-A75H-01A-11R-A32O-07
## ENSG00000000003                        16772                         2814
## ENSG00000000005                            0                            2
## ENSG00000000419                         1422                          733
## ENSG00000000457                          662                          489
## ENSG00000000460                          164                          326
## ENSG00000000938                          190                          234
head(rowData(tcga_data))     # ensembl id and gene id of the first 6 genes.
## DataFrame with 6 rows and 3 columns
##                 ensembl_gene_id external_gene_name original_ensembl_gene_id
##                     <character>        <character>              <character>
## ENSG00000000003 ENSG00000000003             TSPAN6       ENSG00000000003.13
## ENSG00000000005 ENSG00000000005               TNMD        ENSG00000000005.5
## ENSG00000000419 ENSG00000000419               DPM1       ENSG00000000419.11
## ENSG00000000457 ENSG00000000457              SCYL3       ENSG00000000457.12
## ENSG00000000460 ENSG00000000460           C1orf112       ENSG00000000460.15
## ENSG00000000938 ENSG00000000938                FGR       ENSG00000000938.11

Finally, we can use some core functionality of R to save the TCGA_data as a .RDS file. This is faster than repeating the previous operations and useful if you need to work in your data for several days.

# Save the data as a file, if you need it later, you can just load this file
# instead of having to run the whole pipeline again
saveRDS(object = tcga_data,
        file = "tcga_data.RDS",
        compress = FALSE)

The data can be loaded with the following command

tcga_data = readRDS(file = "tcga_data.RDS")

3 RNASeq Normalization

3.1 Defining a pipeline

A typical task on RNA-Seq data is differential expression (DE) analysis, based on some clinical phenotypes. This, in turn, requires normalization of the data, as in its raw format it may have batch effects and other artifacts.

A common approach to such complex tasks is to define a computational pipeline, performing several steps in sequence, allowing the user to select different parameters.

We will now define and run one such pipeline, through the use of an R function.

The function is called limma_pipeline(tcga_data, condition_variable, reference_group), where tcga_data is the data object we have gotten from TCGA and condition_variable is the interesting variable/condition by which you want to group your patient samples. You can also decide which one of the values of your conditional variable is going to be the reference group, with the reference_group parameter.

This function returns a list with three different objects:

  • A complex object, resulting from running voom, this contains the TMM+voom normalized data;
  • A complex object, resulting from running eBayes, this contains the the fitted model plus a number of statistics related to each of the probes;
  • A simple table, resulting from running topTable, which contains the top 100 differentially expressed genes sorted by p-value. This is how the code of this function:
limma_pipeline = function(
  tcga_data,
  condition_variable,
  reference_group=NULL){

  design_factor = colData(tcga_data)[, condition_variable, drop=T]

  group = factor(design_factor)
  if(!is.null(reference_group)){group = relevel(group, ref=reference_group)}

  design = model.matrix(~ group)

  dge = DGEList(counts=assay(tcga_data),
                 samples=colData(tcga_data),
                 genes=as.data.frame(rowData(tcga_data)))

  # filtering
  keep = filterByExpr(dge,design)
  dge = dge[keep,,keep.lib.sizes=FALSE]
  rm(keep)

  # Normalization (TMM followed by voom)
  dge = calcNormFactors(dge)
  v = voom(dge, design, plot=TRUE)

  # Fit model to data given design
  fit = lmFit(v, design)
  fit = eBayes(fit)

  # Show top genes
  topGenes = topTable(fit, coef=ncol(design), number=100, sort.by="p")

  return(
    list(
      voomObj=v, # normalized data
      fit=fit, # fitted model and statistics
      topGenes=topGenes # the 100 most differentially expressed genes
    )
  )
}

Please read the Workflow in Details tab for a step by step explanation of what this pipeline does. For now, copy the limma_pipeline function in your R console, so we can start using it on the TCGA data.

With the following command, we can obtain the DE analysis comparing Primary solid Tumor samples against Solid Tissue Normal. This will be used in the next section, on the classification task.

limma_res = limma_pipeline(
  tcga_data=tcga_data,
  condition_variable="definition",
  reference_group="Solid Tissue Normal"
)

Let’s save this object to file, like we did with tcga_data:

# Save the data as a file, if you need it later, you can just load this file
# instead of having to run the whole pipeline again
saveRDS(object = limma_res,
        file = "limma_res.RDS",
        compress = FALSE)

As an example, we also show here how we can use limma_pipeline to perform DE analysis by grouping patients by gender instead of by tissue type.

gender_limma_res = limma_pipeline(
  tcga_data=tcga_data,
  condition_variable="gender",
  reference_group="female"
)

3.1.1 Visualization

We have also prepared a function that produces PCA plots given the voom object created by the limma_pipeline function. You can inspect this function and try to figure out how it works.

plot_PCA = function(voomObj, condition_variable){
  group = factor(voomObj$targets[, condition_variable])
  pca = prcomp(t(voomObj$E))
  # Take PC1 and PC2 for the plot
  plot(pca$x[,1:2],col=group, pch=19)
  # include a legend for points
  legend("bottomleft", inset=.01, levels(group), pch=19, col=1:length(levels(group)))
  return(pca)
}

By calling the function plot_PCA, we get a plot of the first two principal components:

res_pca = plot_PCA(limma_res$voomObj, "definition")