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 expression) 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("clusterProfiler")
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")])
##   project_id
## 1 TARGET-NBL
## 2 GENIE-GRCC
## 3 GENIE-DFCI
## 4  GENIE-NKI
## 5 GENIE-VICC
## 6  GENIE-UHN
##                                                                  name
## 1                                                       Neuroblastoma
## 2         AACR Project GENIE - Contributed by Institut Gustave Roussy
## 3    AACR Project GENIE - Contributed by Dana-Farber Cancer Institute
## 4    AACR Project GENIE - Contributed by Netherlands Cancer Institute
## 5 AACR Project GENIE - Contributed by Vanderbilt-Ingram Cancer Center
## 6 AACR Project GENIE - Contributed by Princess Margaret Cancer Centre

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")
## $file_count
## [1] 16965
## 
## $data_categories
##   file_count case_count               data_category
## 1       2282        376       Copy Number Variation
## 2       2502        377            Sequencing Reads
## 3       5256        375 Simple Nucleotide Variation
## 4       1290        377             DNA Methylation
## 5        423        377                    Clinical
## 6       1698        376     Transcriptome Profiling
## 7       1634        377                 Biospecimen
## 8        184        184          Proteome Profiling
## 9       1696        371        Structural Variation
## 
## $case_count
## [1] 377
## 
## $file_size
## [1] 2.671516e+13

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
  data.type = "Gene Expression Quantification",
  experimental.strategy = "RNA-Seq",
  #workflow.type = "HTSeq - Counts"
  workflow.type = "STAR - 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] "id"                        "data_format"              
##  [3] "cases"                     "access"                   
##  [5] "file_name"                 "submitter_id"             
##  [7] "data_category"             "type"                     
##  [9] "file_size"                 "created_datetime"         
## [11] "md5sum"                    "updated_datetime"         
## [13] "file_id"                   "data_type"                
## [15] "state"                     "experimental_strategy"    
## [17] "version"                   "data_release"             
## [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] "Primary Tumor"       "Solid Tissue Normal" "Primary Tumor"

sample_type is a categorical variable, which can be better visualized with the table function, so that we can check how many different categories exist in the data.

table(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 Tumors). For simplicity, we will ignore the small class of recurrent tumors. Therefore, we will redo the query as

query_TCGA = GDCquery(
  project = "TCGA-LIHC",
  data.category = "Transcriptome Profiling", # parameter enforced by GDCquery
  data.type = "Gene Expression Quantification",
  experimental.strategy = "RNA-Seq",
  workflow.type = "STAR - 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 starting 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] 60660   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. But be careful if you have large objects. The colnames() and rownames() functions return all column and row names, even if you have tens of thousands. The combination header(colnames()) can make more sense in these cases.

# 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"              "oct_embedded"               
##  [9] "sample_id"                   "submitter_id"               
## [11] "state"                       "is_ffpe"                    
## [13] "sample_type"                 "tissue_type"                
## [15] "tumor_descriptor"            "composition"                
## [17] "days_to_collection"          "initial_weight"             
## [19] "pathology_report_uuid"       "synchronous_malignancy"     
## [21] "ajcc_pathologic_stage"       "days_to_diagnosis"          
## [23] "treatments"                  "last_known_disease_status"  
## [25] "tissue_or_organ_of_origin"   "days_to_last_follow_up"     
## [27] "age_at_diagnosis"            "primary_diagnosis"          
## [29] "prior_malignancy"            "year_of_diagnosis"          
## [31] "prior_treatment"             "ajcc_staging_system_edition"
## [33] "ajcc_pathologic_t"           "morphology"                 
## [35] "ajcc_pathologic_n"           "ajcc_pathologic_m"          
## [37] "classification_of_tumor"     "diagnosis_id"               
## [39] "icd_10_code"                 "site_of_resection_or_biopsy"
## [41] "tumor_grade"                 "progression_or_recurrence"  
## [43] "alcohol_history"             "exposure_id"                
## [45] "race"                        "gender"                     
## [47] "ethnicity"                   "vital_status"               
## [49] "age_at_index"                "days_to_birth"              
## [51] "year_of_birth"               "demographic_id"             
## [53] "year_of_death"               "days_to_death"              
## [55] "bcr_patient_barcode"         "primary_site"               
## [57] "project_id"                  "disease_type"               
## [59] "name"                        "releasable"                 
## [61] "released"                    "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$ajcc_pathologic_stage)
## 
##    Stage I   Stage II  Stage III Stage IIIA Stage IIIB Stage IIIC   Stage IV 
##        189         97          6         73          8         10          3 
##  Stage IVA  Stage IVB 
##          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] 60660   421
head(assay(tcga_data)[,1:10]) # expression of first 6 genes and first 10 samples
##                    TCGA-DD-AACJ-01A-11R-A41C-07 TCGA-DD-A4NE-01A-11R-A27V-07
## ENSG00000000003.15                         3311                         3132
## ENSG00000000005.6                             2                            0
## ENSG00000000419.13                         1216                         1596
## ENSG00000000457.14                          758                         1241
## ENSG00000000460.17                          156                          530
## ENSG00000000938.13                           91                           66
##                    TCGA-ED-A8O6-01A-11R-A36F-07 TCGA-MI-A75I-01A-11R-A32O-07
## ENSG00000000003.15                         7510                         5443
## ENSG00000000005.6                             0                            0
## ENSG00000000419.13                         1904                         1217
## ENSG00000000457.14                          719                          432
## ENSG00000000460.17                          233                          138
## ENSG00000000938.13                          106                          133
##                    TCGA-FV-A3I0-11A-11R-A22L-07 TCGA-G3-AAV4-01A-11R-A38B-07
## ENSG00000000003.15                         3957                          839
## ENSG00000000005.6                             0                           18
## ENSG00000000419.13                          873                         1475
## ENSG00000000457.14                          316                          568
## ENSG00000000460.17                           57                          201
## ENSG00000000938.13                          278                           64
##                    TCGA-BC-4072-01B-11R-A155-07 TCGA-2Y-A9H6-01A-11R-A39D-07
## ENSG00000000003.15                         4823                         8519
## ENSG00000000005.6                             0                            5
## ENSG00000000419.13                         2426                          937
## ENSG00000000457.14                          773                          536
## ENSG00000000460.17                          220                          122
## ENSG00000000938.13                          958                          343
##                    TCGA-DD-A1EC-11A-11R-A131-07 TCGA-CC-5260-01A-01R-A131-07
## ENSG00000000003.15                         3414                         4241
## ENSG00000000005.6                             5                            1
## ENSG00000000419.13                          819                          786
## ENSG00000000457.14                          416                          535
## ENSG00000000460.17                          466                          200
## ENSG00000000938.13                          274                          507
head(rowData(tcga_data))     # ensembl id and gene id of the first 6 genes.
## DataFrame with 6 rows and 10 columns
##                      source     type     score     phase            gene_id
##                    <factor> <factor> <numeric> <integer>        <character>
## ENSG00000000003.15   HAVANA     gene        NA        NA ENSG00000000003.15
## ENSG00000000005.6    HAVANA     gene        NA        NA  ENSG00000000005.6
## ENSG00000000419.13   HAVANA     gene        NA        NA ENSG00000000419.13
## ENSG00000000457.14   HAVANA     gene        NA        NA ENSG00000000457.14
## ENSG00000000460.17   HAVANA     gene        NA        NA ENSG00000000460.17
## ENSG00000000938.13   HAVANA     gene        NA        NA ENSG00000000938.13
##                         gene_type   gene_name       level     hgnc_id
##                       <character> <character> <character> <character>
## ENSG00000000003.15 protein_coding      TSPAN6           2  HGNC:11858
## ENSG00000000005.6  protein_coding        TNMD           2  HGNC:17757
## ENSG00000000419.13 protein_coding        DPM1           2   HGNC:3005
## ENSG00000000457.14 protein_coding       SCYL3           2  HGNC:19285
## ENSG00000000460.17 protein_coding    C1orf112           2  HGNC:25565
## ENSG00000000938.13 protein_coding         FGR           2   HGNC:3697
##                             havana_gene
##                             <character>
## ENSG00000000003.15 OTTHUMG00000022002.2
## ENSG00000000005.6  OTTHUMG00000022001.2
## ENSG00000000419.13 OTTHUMG00000032742.2
## ENSG00000000457.14 OTTHUMG00000035941.6
## ENSG00000000460.17 OTTHUMG00000035821.9
## ENSG00000000938.13 OTTHUMG00000003516.3

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

Basics

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

Details

In our pipeline function, we use the package limma. We will select a particular clinical feature of the data to use as class for grouping the samples as either tumor vs normal tissue. This data is available under the column definition for tcga_data, but needs the use of the function colData to be accessed. In addition, limma requires this data to be a factor, so we convert it as such:

clinical_data = colData(tcga_data)
group = factor(clinical_data$definition)

As seen before, we have 2 distinct groups of tissues defined in this column, Solid Tissue Normal (our control samples) and Primary solid Tumor (our samples of interest). In this factor, we also want to define Solid Tissue Normal as being our base or reference level.

group = relevel(group, ref="Solid Tissue Normal")

Next, we need to create a design matrix, which will indicate the conditions to be compared by the DE analysis. The ~ symbol represents that we are constructing a formula.

design = model.matrix(~group)
head(design)
##   (Intercept) groupPrimary solid Tumor
## 1           1                        1
## 2           1                        1
## 3           1                        1
## 4           1                        1
## 5           1                        0
## 6           1                        1

Before performing DE analysis, we remove genes, which have low amount of counts. We transform our tcga_data object as DGEList, which provides functions for filtering. By default genes with counts with less than 10 reads are removed.

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

# filtering
keep = filterByExpr(dge,design) # defining which genes to keep
dge = dge[keep,,keep.lib.sizes=FALSE] # filtering the dge object
rm(keep) #  use rm() to remove objects from memory if you don't need them anymore

Before we fit a model to our data, we normalize the data to minimize batch effects and technical variation with the TMM (trimmed mean of M-values) normalization method. Moreover, to apply limma on RNA-seq, we need to convert the data to have a similar variance as arrays. This is done with the VOOM method.

dge = calcNormFactors(dge, method = "TMM")
v = voom(dge, design, plot=TRUE)