Downloading data from TCGA data portal

The Bioconductor package TCGAbiolinks (Colaprico et al. 2016) has three main functions GDCquery, GDCdownload and GDCprepare that should sequentially be used to respectively search, download and load the data as an R object.

GDCquery uses GDC API to search the data for a given project and data category and filters the results by samples, sample type, file type and others features if requested by the user. This function returns an object with a summary table with the results found (samples, files and other useful information) and the arguments used in the query. The most important GDCquery arguments are project which receives a GDC project (TCGA-USC, TCGA-LGG, TARGET-AML, etc), data.category which receives a data category (Transcriptome Profiling, Copy Number Variation, DNA methylation, Gene expression, etc), data.type which receives a data type (Gene expression quantification, Isoform Expression Quantification, miRNA Expression Quantification, Copy Number Segment, Masked Copy Number Segment, etc), workflow.type, which receives a GDC workflow type (HTSeq - Counts, HTSeq - FPKM-UQ, HTSeq - FPKM), legacy, which selects to use the legacy database or the harmonized database, file.type, which receives a file type for the searches in the legacy database (hg18.seg, hg19.seg, nocnv_,hg18.seg, nocnv_hg19.seg, rsem.genes.results, rsem.genes.normalized_results, etc) and platform, which receives a platform for the searches in the legacy database (HumanMethylation27, Genome_Wide_SNP_6, IlluminaHiSeq_RNASeqV2, etc). A complete list of possible entries for arguments can be found in the TCGAbiolinks vignette. Listing 1 shows an example of this function.

After the search step, the user will be able to download the data using the GDCdownload function which can use either the GDC API to download the samples, or the gdc client tools. The downloaded data will be saved in a directory with the project name and a sub-folder with the data.category, for example “TCGA-GBM/DNA_methylation”.

Finally, GDCprepare transforms the downloaded data into a summarizedExperiment object (Huber et al. 2015) or a data frame. If SummarizedExperiment is set to TRUE, TCGAbiolinks will add to the object sub-type information, which was defined by The Cancer Genome Atlas (TCGA) Research Network reports (the full list of papers can be seen in TCGAquery_subtype section in TCGAbiolinks vignette), and clinical information. Listing 1 shows how to use these functions to download DNA methylation and gene expression data from the GDC legacy database and 2 shows how to download copy number variation from harmonized data portal. Other examples, that access the harmonized data can be found in the TCGAbiolinks vignette.

library(TCGAbiolinks)
# Obs: The data in the legacy database has been aligned to hg19
query.met.gbm <- GDCquery(project = "TCGA-GBM", 
                          legacy = TRUE,
                          data.category = "DNA methylation",
                          platform = "Illumina Human Methylation 450", 
                          barcode = c("TCGA-76-4926-01B-01D-1481-05", "TCGA-28-5211-01C-11D-1844-05"))
GDCdownload(query.met.gbm)

met.gbm.450 <- GDCprepare(query = query.met.gbm,
                          save = TRUE, 
                          save.filename = "gbmDNAmet450k.rda",
                          summarizedExperiment = TRUE)

query.met.lgg <- GDCquery(project = "TCGA-LGG", 
                          legacy = TRUE,
                          data.category = "DNA methylation",
                          platform = "Illumina Human Methylation 450",
                          barcode = c("TCGA-HT-7879-01A-11D-2399-05", "TCGA-HT-8113-01A-11D-2399-05"))
GDCdownload(query.met.lgg)
met.lgg.450 <- GDCprepare(query = query.met.lgg,
                          save = TRUE, 
                          save.filename = "lggDNAmet450k.rda",
                          summarizedExperiment = TRUE)
met.gbm.lgg <- SummarizedExperiment::cbind(met.lgg.450, met.gbm.450)


query.exp.lgg <- GDCquery(project = "TCGA-LGG", 
                          legacy = TRUE,
                          data.category = "Gene expression",
                          data.type = "Gene expression quantification",
                          platform = "Illumina HiSeq", 
                          file.type = "results",
                          sample.type = "Primary solid Tumor")
GDCdownload(query.exp.lgg)
exp.lgg <- GDCprepare(query = query.exp.lgg, save = TRUE, save.filename = "lggExp.rda")

query.exp.gbm <- GDCquery(project = "TCGA-GBM", 
                          legacy = TRUE,
                          data.category = "Gene expression",
                          data.type = "Gene expression quantification",
                          platform = "Illumina HiSeq", 
                          file.type = "results",
                          sample.type = "Primary solid Tumor")
GDCdownload(query.exp.gbm)
exp.gbm <- GDCprepare(query = query.exp.gbm, save = TRUE, save.filename = "gbmExp.rda")
exp.gbm.lgg <- SummarizedExperiment::cbind(exp.lgg, exp.gbm)

#-----------------------------------------------------------------------------
#                   Data.category: Copy number variation aligned to hg38
#-----------------------------------------------------------------------------
query <- GDCquery(project = "TCGA-ACC",
                  data.category = "Copy Number Variation",
                  data.type = "Copy Number Segment",
                  barcode = c( "TCGA-OR-A5KU-01A-11D-A29H-01", "TCGA-OR-A5JK-01A-11D-A29H-01"))
GDCdownload(query)
data <- GDCprepare(query)

query <- GDCquery("TCGA-ACC",
                  "Copy Number Variation",
                  data.type = "Masked Copy Number Segment",
                  sample.type = c("Primary solid Tumor")) # see the barcodes with getResults(query)$cases
GDCdownload(query)
data <- GDCprepare(query)

If a SummarizedExperiment object was chosen, the data can be accessed with three different accessors: assay for the data information, rowRanges to gets the range of values in each row and colData to get the sample information (patient, batch, sample type, etc) (Huber et al. 2015; H. J. Morgan M Obenchain V and H., n.d.). An example is shown in listing below.

library(SummarizedExperiment)

# Load object from TCGAWorkflowData package
# THis object will be created in the further sections,
data(GBMIllumina_HiSeq) 

# get expression matrix
data <- assay(gbm.exp)
datatable(data[1:10,], 
          options = list(scrollX = TRUE, keys = TRUE, pageLength = 5), 
          rownames = TRUE)

# get genes information
genes.info <- rowRanges(gbm.exp)
genes.info

## GRanges object with 21022 ranges and 4 metadata columns:
##                          seqnames              ranges strand |      gene_id
##                             <Rle>           <IRanges>  <Rle> |  <character>
##                   A1BG|1    chr19   58856544-58864865      - |         A1BG
##                    A2M|2    chr12     9220260-9268825      - |          A2M
##                   NAT1|9     chr8   18027986-18081198      + |         NAT1
##                  NAT2|10     chr8   18248755-18258728      + |         NAT2
##              SERPINA3|12    chr14   95058395-95090983      + | RP11-986E7.7
##                      ...      ...                 ...    ... .          ...
##     NCRNA00182|100302692     chrX   73183790-73513409      - |          FTX
##   TMED7-TICAM2|100302736     chr5 114914339-114961858      - | TMED7-TICAM2
##   TMED7-TICAM2|100302736     chr5 114949205-114968689      - |        TMED7
##   TMED7-TICAM2|100302736     chr5 114914339-114961876      - |       TICAM2
##   LOC100303728|100303728     chrX 118599997-118603061      - |  SLC25A5-AS1
##                          entrezgene ensembl_gene_id
##                           <numeric>     <character>
##                   A1BG|1          1 ENSG00000121410
##                    A2M|2          2 ENSG00000175899
##                   NAT1|9          9 ENSG00000171428
##                  NAT2|10         10 ENSG00000156006
##              SERPINA3|12         12 ENSG00000273259
##                      ...        ...             ...
##     NCRNA00182|100302692  100302692 ENSG00000230590
##   TMED7-TICAM2|100302736  100302736 ENSG00000251201
##   TMED7-TICAM2|100302736  100302736 ENSG00000134970
##   TMED7-TICAM2|100302736  100302736 ENSG00000243414
##   LOC100303728|100303728  100303728 ENSG00000224281
##                                                                    transcript_id.transcript_id_TCGA-28-1753-01A-01R-1850-01
##                                                                                                                 <character>
##                   A1BG|1                                                                              uc002qsd.3,uc002qsf.1
##                    A2M|2                                                                   uc001qvj.1,uc001qvk.1,uc009zgk.1
##                   NAT1|9 uc003wyq.2,uc003wyr.2,uc003wys.2,uc003wyt.2,uc003wyu.2,uc003wyv.2,uc010ltc.2,uc010ltd.2,uc011kyl.1
##                  NAT2|10                                                                                         uc003wyw.1
##              SERPINA3|12                       uc001ydo.3,uc001ydp.2,uc001ydq.2,uc001ydr.2,uc001yds.2,uc010avf.1,uc010avg.2
##                      ...                                                                                                ...
##     NCRNA00182|100302692                                                                              uc004ebr.1,uc010nlq.1
##   TMED7-TICAM2|100302736                                                                              uc003krd.2,uc003kre.2
##   TMED7-TICAM2|100302736                                                                              uc003krd.2,uc003kre.2
##   TMED7-TICAM2|100302736                                                                              uc003krd.2,uc003kre.2
##   LOC100303728|100303728                                                                              uc004ere.1,uc004erg.1
##   -------
##   seqinfo: 24 sequences from an unspecified genome; no seqlengths

# get sample information
sample.info <- colData(gbm.exp)
datatable(as.data.frame(sample.info), 
          options = list(scrollX = TRUE, keys = TRUE, pageLength = 5), 
          rownames = FALSE)

The clinical data can be obtained using TCGAbiolinks through two methods. The first one will download only the indexed GDC clinical data which includes diagnoses (vital status, days to death, age at diagnosis, days to last follow up, days to recurrence), treatments (days to treatment, treatment id, therapeutic agents, treatment intent type), demographic (gender, race, ethnicity) and exposures (cigarettes per day, weight, height, alcohol history) information. This indexed clinical data can be obtained using the function GDCquery_clinical which can be used as described in listing below. This function has two arguments project (“TCGA-GBM”,“TARGET-AML”,etc) and type (“Clinical” or “Biospecimen”). The second method will download the XML files with all clinical data for the patient and retrieve the desired information from it. This will give access to all clinical data available which includes patient (tumor tissue site, histological type, gender, vital status, days to birth, days to last follow up, etc), drug (days to drug therapy start, days to drug therapy end, therapy types, drug name), radiation (days to radiation therapy start, days to radiation therapy end, radiation type, radiation dosage ), new tumor event (days to new tumor event after initial treatment, new neoplasm event type, additional pharmaceutical therapy), follow up (primary therapy outcome success, follow up treatment success, vital status, days to last follow up, date of form completion), stage event (pathologic stage), admin (batch number, project code, disease code, Biospecimen Core Resource).

# get indexed clinical patient data for GBM samples
gbm_clin <- GDCquery_clinic(project = "TCGA-GBM", type = "Clinical")

# get indexed clinical patient data for LGG samples
lgg_clin <- GDCquery_clinic(project = "TCGA-LGG", type = "Clinical")

# Bind the results, as the columns might not be the same,
# we will will plyr rbind.fill, to have all columns from both files
clinical <- plyr::rbind.fill(gbm_clin,lgg_clin)

datatable(clinical[1:10,], options = list(scrollX = TRUE, keys = TRUE), rownames = FALSE)

# Fetch clinical data directly from the clinical XML files.
# if barcode is not set, it will consider all samples.
# We only set it to make the example faster
query <- GDCquery(project = "TCGA-GBM",
                  file.type = "xml",
                  data.category = "Clinical",
                  barcode = c("TCGA-08-0516","TCGA-02-0317")) 
GDCdownload(query)
clinical <- GDCprepare_clinic(query, clinical.info = "patient")

datatable(clinical, options = list(scrollX = TRUE, keys = TRUE), rownames = FALSE)

clinical.drug <- GDCprepare_clinic(query, clinical.info = "drug")

datatable(clinical.drug, options = list(scrollX = TRUE, keys = TRUE), rownames = FALSE)

clinical.radiation <- GDCprepare_clinic(query, clinical.info = "radiation")

datatable(clinical.radiation, options = list(scrollX = TRUE,  keys = TRUE), rownames = FALSE)

clinical.admin <- GDCprepare_clinic(query, clinical.info = "admin")

datatable(clinical.admin, options = list(scrollX = TRUE, keys = TRUE), rownames = FALSE)

Mutation information is stored in two types of Mutation Annotation Format (MAF): Protected and Somatic (or Public) MAF files, which are derived from the GDC annotated VCF files. Annotated VCF files often have variants reported on multiple transcripts whereas the protected MAF (*protected.maf) only reports the most critically affected one and the Somatic MAFs (*somatic.maf) are further processed to remove low quality and potential germline variants. To download Somatic MAFs data using TCGAbiolinks, GDCquery_maf function is provided (see listing below).

LGGmut <- GDCquery_Maf(tumor = "LGG", pipelines = "mutect2")

data(mafMutect2LGGGBM)
datatable(LGGmut[1:10,], options = list(scrollX = TRUE, keys = TRUE), rownames = FALSE)

Finally, the Cancer Genome Atlas (TCGA) Research Network has reported integrated genome-wide studies of various diseases (ACC (Zheng et al. 2016), BRCA (C. G. A. Network and others 2012b), COAD (C. G. A. Network and others 2012a), GBM (Ceccarelli et al. 2016), HNSC (C. G. A. Network and others 2015a), KICH (Davis et al. 2014), KIRC (Network and others 2013), KIRP (Network and others 2016), LGG (Ceccarelli et al. 2016), LUAD (Network and others 2014c), LUSC (Network and others 2012), PRAD (Network and others 2015), READ (C. G. A. Network and others 2012a), SKCM (C. G. A. Network and others 2015b), STAD (Network and others 2014a), THCA (Network and others 2014d) and UCEC (Network and others 2014b)) which classified them in different subtypes. This classification can be retrieved using the TCGAquery_subtype function or by accessing the samples information in the SummarizedExperiment object that created by the GDCprepare function, which automatically incorporates it into the object.

gbm.subtypes <- TCGAquery_subtype(tumor = "gbm")

datatable(gbm.subtypes[1:10,], options = list(scrollX = TRUE, keys = TRUE), rownames = FALSE)

Downloading data from Broad TCGA GDAC

The Bioconductor package RTCGAToolbox (Samur 2014) provides access to Firehose Level 3 and 4 data through the function getFirehoseData. The following arguments allow users to select the version and tumor type of interest:

• dataset - Tumor to download. A complete list of possibilities can be view with getFirehoseDatasets function.

• runDate - Stddata run dates. Dates can be viewed with getFirehoseRunningDates function.

• gistic2_Date - Analyze run dates. Dates can viewed with getFirehoseAnalyzeDates function.

These arguments can be used to select the data type to download: RNAseq_Gene, Clinic, miRNASeq_Gene, ccRNAseq2_Gene_Norm, CNA_SNP, CNV_SNP, CNA_Seq, CNA_CGH, Methylation, Mutation, mRNA_Array , miRNA_Array, and RPPA.

By default, RTCGAToolbox allows users to download up to 500 MB worth of data. To increase the size of the download, users are encouraged to use fileSizeLimit argument. An example is found in listing below. The getData function allows users to access the downloaded data (see lines 22-24 of listing below as an S4Vector object.

library(RTCGAToolbox)

# Get the last run dates
lastRunDate <- getFirehoseRunningDates()[1]

# get DNA methylation data, RNAseq2 and clinical data for GBM
gbm.data <- getFirehoseData(dataset = "GBM",
                            runDate = lastRunDate, 
                            gistic2Date = getFirehoseAnalyzeDates(1),
                            Methylation = FALSE,  
                            clinical = TRUE,
                            RNASeq2GeneNorm  = FALSE, 
                            Mutation = TRUE,
                            fileSizeLimit = 10000)

gbm.mut <- getData(gbm.data,"Mutation")
gbm.clin <- getData(gbm.data,"clinical")

Finally, RTCGAToolbox can access level 4 data, which can be handy when the user requires GISTIC results. GISTIC is used to identify genes targeted by somatic copy-number alterations (SCNAs) (Mermel et al. 2011).

# Download GISTIC results
lastanalyzedate <- getFirehoseAnalyzeDates(1)
gistic <- getFirehoseData("GBM",gistic2Date = lastanalyzedate)

# get GISTIC results
gistic.allbygene <- getData(gistic, type = "GISTIC", platform = "AllByGene")
gistic.thresholedbygene <- getData(gistic, type = "GISTIC", platform = "ThresholdedByGene")

data(GBMGistic)
datatable(gistic.allbygene,
          filter = 'top',
          options = list(scrollX = TRUE, keys = TRUE, pageLength = 5), 
          rownames = FALSE)

datatable(gistic.thresholedbygene,
          filter = 'top',
          options = list(scrollX = TRUE, keys = TRUE, pageLength = 5), 
          rownames = FALSE)

Identification of recurrent CNV in cancer

Cancer-related CNV has to be present in many of the analyzed genomes. The most significant recurrent CNV were identified using GAIA (Morganella, Pagnotta, and Ceccarelli, n.d.), an iterative procedure where a statistical hypothesis framework is extended to take into account within-sample homogeneity. GAIA is based on a conservative permutation test allowing the estimation of the probability distribution of the contemporary mutations expected for non-driver markers.

Segmented data retrieved from TCGA were used to generate a matrix including all needed information about the observed aberrant regions. Furthermore, GAIA requires genomic probes metadata (specific for each CNV technology), that can be downloaded from the broadinstitute website.

# -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=--=-==--=--==---=-=-=-=-=-=-=-=-=-=-=-=-=-=-=--=-==--=--==--
# Retrieve probes meta file from broadinstitute website for hg19
# For hg38 analysis please take a look on:
# https://gdc.cancer.gov/about-data/data-harmonization-and-generation/gdc-reference-files
# File: SNP6 GRCh38 Liftover Probeset File for Copy Number Variation Analysis
# -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=--=-==--=--==---=-=-=-=-=-=-=-=-=-=-=-=-=-=-=--=-==--=--==--
gdac.root <- "ftp://ftp.broadinstitute.org/pub/GISTIC2.0/hg19_support/"
file <- paste0(gdac.root, "genome.info.6.0_hg19.na31_minus_frequent_nan_probes_sorted_2.1.txt")
if(!file.exists(basename(file))) downloader::download(file, basename(file))
markersMatrix <-  readr::read_tsv(basename(file), col_names = FALSE, col_types = "ccn", progress = FALSE)
save(markersMatrix, file = "markersMatrix.rda", compress = "xz")

cancer <- "GBM"
message(paste0("Starting ", cancer))

# get objects created above
data(GBMnocnvhg19)
data(markersMatrix)

# Add label (0 for loss, 1 for gain)
cnvMatrix <- cbind(cnvMatrix,Label=NA)
cnvMatrix[cnvMatrix[,"Segment_Mean"] < -0.3,"Label"] <- 0
cnvMatrix[cnvMatrix[,"Segment_Mean"] > 0.3,"Label"] <- 1
cnvMatrix <- cnvMatrix[!is.na(cnvMatrix$Label),]

# Remove "Segment_Mean" and change col.names
cnvMatrix <- cnvMatrix[,-6]
colnames(cnvMatrix) <- c("Sample.Name", "Chromosome", "Start", "End", "Num.of.Markers", "Aberration")

# Substitute Chromosomes "X" and "Y" with "23" and "24"
cnvMatrix[cnvMatrix$Chromosome == "X","Chromosome"] <- 23
cnvMatrix[cnvMatrix$Chromosome == "Y","Chromosome"] <- 24
cnvMatrix$Chromosome <- as.integer(cnvMatrix$Chromosome)

# Recurrent CNV identification with GAIA
colnames(markersMatrix) <- c("Probe.Name", "Chromosome", "Start")
unique(markersMatrix$Chromosome)

##  [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14"
## [15] "15" "16" "17" "18" "19" "20" "21" "22" "X"  "Y"

markersMatrix[markersMatrix$Chromosome == "X","Chromosome"] <- 23
markersMatrix[markersMatrix$Chromosome == "Y","Chromosome"] <- 24
markersMatrix$Chromosome <- as.integer(markersMatrix$Chromosome)
markerID <- paste(markersMatrix$Chromosome,markersMatrix$Start, sep = ":")
# Removed duplicates
markersMatrix <- markersMatrix[!duplicated(markerID),]
# Filter markersMatrix for common CNV
markerID <- paste(markersMatrix$Chromosome,markersMatrix$Start, sep = ":")

file <- "ftp://ftp.broadinstitute.org/pub/GISTIC2.0/hg19_support/CNV.hg19.bypos.111213.txt"
if(!file.exists(basename(file))) downloader::download(file, basename(file))
commonCNV <- readr::read_tsv(basename(file), progress = FALSE)
commonID <- paste(commonCNV$Chromosome,commonCNV$Start, sep = ":")
markersMatrix_fil <- markersMatrix[!markerID %in% commonID,]

library(gaia)
set.seed(200)
markers_obj <- load_markers(as.data.frame(markersMatrix_fil))
nbsamples <- length(unique(cnvMatrix$Sample))
cnv_obj <- load_cnv(cnvMatrix, markers_obj, nbsamples)
suppressWarnings({
  results <- runGAIA(cnv_obj,
                     markers_obj,
                     output_file_name = paste0("GAIA_",cancer,"_flt.txt"),
                     aberrations = -1,  # -1 to all aberrations
                     chromosomes = 9, # -1 to all chromosomes
                     approximation = TRUE, # Set to TRUE to speed up the time requirements
                     num_iterations = 5000, # Reduced to speed up the time requirements
                     threshold = 0.25)
})
# Set q-value threshold
# Use a smalled value for your analysis. We set this as high values
# due to the small number of samples which did not reproduced
# results with smaller q-values
threshold <- 0.3

# Plot the results
RecCNV <- t(apply(results,1,as.numeric))
colnames(RecCNV) <- colnames(results)
RecCNV <- cbind(RecCNV, score = 0)
minval <- format(min(RecCNV[RecCNV[,"q-value"] != 0,"q-value"]), scientific = FALSE)
minval <- substring(minval,1, nchar(minval) - 1)
RecCNV[RecCNV[,"q-value"] == 0,"q-value"] <- as.numeric(minval)
RecCNV[,"score"] <- sapply(RecCNV[,"q-value"],function(x) -log10(as.numeric(x)))
RecCNV[RecCNV[,"q-value"] == as.numeric(minval),]

##   Chromosome Aberration Kind Region Start [bp] Region End [bp]
## 1          9               0          18259823        18387361
## 2          9               0          19290143        19485455
## 3          9               0          19535498        26257238
## 4          9               0          26433917        26725609
## 5          9               0          26834232        27009100
##   Region Size [bp]    q-value    score
## 1           127539 0.00875417 2.057785
## 2           195313 0.00875417 2.057785
## 3          6721741 0.00875417 2.057785
## 4           291693 0.00875417 2.057785
## 5           174869 0.00875417 2.057785

gaiaCNVplot(RecCNV,threshold)

save(results, RecCNV, threshold, file = paste0(cancer,"_CNV_results.rda"))

Recurrent amplifications and deletions were identified for GBM, and represented in chromosomal overview plots by a statistical score (\(-log_{10}\) corrected p-value for amplifications and \(log_{10}\) corrected p-value for deletions). Genomic regions identified as significantly altered in copy number (corrected p-value < \(10^{-4}\)) were then annotated to report amplified and deleted genes potentially related to cancer.

Gene annotation of recurrent CNV

The aberrant recurrent genomic regions in cancer, as identified by GAIA, have to be annotated to verify which genes are significantly amplified or deleted. Using biomaRt we retrieved the genomic ranges of all human genes and we compared them with significant aberrant regions to select full length genes.

library(GenomicRanges)
# Get gene information from GENCODE using biomart
genes <- TCGAbiolinks:::get.GRCh.bioMart(genome = "hg19") 
genes <- genes[genes$external_gene_id != "" & genes$chromosome_name %in% c(1:22,"X","Y"),]
genes[genes$chromosome_name == "X", "chromosome_name"] <- 23
genes[genes$chromosome_name == "Y", "chromosome_name"] <- 24
genes$chromosome_name <- sapply(genes$chromosome_name,as.integer)
genes <- genes[order(genes$start_position),]
genes <- genes[order(genes$chromosome_name),]
genes <- genes[,c("external_gene_id", "chromosome_name", "start_position","end_position")]
colnames(genes) <- c("GeneSymbol","Chr","Start","End")
genes_GR <- makeGRangesFromDataFrame(genes,keep.extra.columns = TRUE)
save(genes_GR,genes,file = "genes_GR.rda", compress = "xz")

##############################
## Recurrent CNV annotation ## 
##############################
# Get gene information from GENCODE using biomart
data(genes_GR) # downloaded in the previous step (available in TCGAWorkflowData)

load(paste0(cancer,"_CNV_results.rda"))
sCNV <- RecCNV[RecCNV[,"q-value"] <= threshold,c(1:4,6)]
sCNV <- sCNV[order(sCNV[,3]),]
sCNV <- sCNV[order(sCNV[,1]),]
colnames(sCNV) <- c("Chr","Aberration","Start","End","q-value")
sCNV_GR <- makeGRangesFromDataFrame(sCNV,keep.extra.columns = TRUE)

hits <- findOverlaps(genes_GR, sCNV_GR, type = "within")
sCNV_ann <- cbind(sCNV[subjectHits(hits),],genes[queryHits(hits),])
AberrantRegion <- paste0(sCNV_ann[,1],":",sCNV_ann[,3],"-",sCNV_ann[,4])
GeneRegion <- paste0(sCNV_ann[,7],":",sCNV_ann[,8],"-",sCNV_ann[,9])
AmpDel_genes <- cbind(sCNV_ann[,c(6,2,5)],AberrantRegion,GeneRegion)
AmpDel_genes[AmpDel_genes[,2] == 0,2] <- "Del"
AmpDel_genes[AmpDel_genes[,2] == 1,2] <- "Amp"
rownames(AmpDel_genes) <- NULL

save(RecCNV, AmpDel_genes, file = paste0(cancer,"_CNV_results.rda"))

Visualizing multiple genomic alteration events

In order to visualize multiple genomic alteration events, we recommend using maftools plot which is provided by Bioconductor package maftools (Mayakonda and Koeffler 2016). The listing below shows how to download mutation data using GDCquery_maf (line 4) and prepare it to use with maftools.

library(maftools)
# Download Mutation Annotation Format (MAF) files
LGGmut <- GDCquery_Maf(tumor = "LGG", pipelines = "mutect2")
GBMmut <- GDCquery_Maf(tumor = "GBM", pipelines = "mutect2")

# Merge them 
mut <- plyr::rbind.fill(LGGmut, GBMmut)
save(maf,file ="mafMutect2LGGGBM.rda", compress = "xz")

The function read.maf is used to prepare the MAF data to be used with maftools. We also added clinical information that will be used in survival plots.

library(maftools)
# recovering data from TCGAWorkflowData package.
data(mafMutect2LGGGBM)

# To prepare for maftools we will also include clinical data
# For a mutant vs WT survival analysis 
# get indexed clinical patient data for GBM samples
gbm_clin <- GDCquery_clinic(project = "TCGA-GBM", type = "Clinical")
# get indexed clinical patient data for LGG samples
lgg_clin <- GDCquery_clinic(project = "TCGA-LGG", type = "Clinical")
# Bind the results, as the columns might not be the same,
# we will will plyr rbind.fill, to have all columns from both files
clinical <- plyr::rbind.fill(gbm_clin,lgg_clin)
colnames(clinical)[1] <- "Tumor_Sample_Barcode"
clinical$Overall_Survival_Status <- 1
clinical$Overall_Survival_Status[which(clinical$vital_status != "dead")] <- 0
clinical$time <- clinical$days_to_death
clinical$time[is.na(clinical$days_to_death)] <- clinical$days_to_last_follow_up[is.na(clinical$days_to_death)]

# Create object to use in maftools
maf <- read.maf(maf = mut, clinicalData = clinical, isTCGA = T)

We can plot a MAF summary.

plotmafSummary(maf = maf,
               rmOutlier = TRUE,
               addStat = 'median',
               dashboard = TRUE)

We can draw oncoplot with the top 20 most mutated genes and add metadata information such as molecular subtypes information.

gbm.subtypes <- TCGAquery_subtype(tumor = "gbm")
lgg.subtypes <- TCGAquery_subtype(tumor = "lgg")
subtype.info <- plyr::rbind.fill( lgg.subtypes, gbm.subtypes)
colnames(subtype.info)[1] <- "Tumor_Sample_Barcode" # required
oncoplot(maf = maf,
         top = 20,
         legendFontSize = 8,
         clinicalFeatures = c("Supervised.DNA.Methylation.Cluster","Original.Subtype"), # which columns should be plotted
         annotationDat = subtype.info) # data frame with metadata

We can also perform survival analysis by grouping samples from MAF based on mutation status of given gene(s).

plot <- mafSurvival(maf = maf,
                   genes = "IDH1",
                   fn = NULL,
                   time = 'time',
                   Status = 'Overall_Survival_Status',
                   isTCGA = TRUE)

## IDH1 
##  413 
##     Group medianTime   N
## 1: Mutant        748 412
## 2:     WT        372 485

Overview of genomic alterations by circos plot

Genomic alterations in cancer, including CNV and mutations, can be represented in an effective overview plot named circos. We used circlize CRAN package to represent significant CNV (resulting from GAIA analysis) and recurrent mutations (selecting curated genetic variations retrieved from TCGA that are identified in at least two tumor samples) in LGG. Circos plot can illustrate molecular alterations genome-wide or only in one or more selected chromosomes.

LGGmut <- GDCquery_Maf(tumor = "LGG", pipelines = "mutect2")

###############################################
## Genomic aberration overview - Circos plot ## 
###############################################
# Retrieve curated mutations for selected cancer (e.g. "LGG") 
data(mafMutect2LGGGBM)
# Select only potentially damaging mutations
LGGmut <- LGGmut[LGGmut$Variant_Classification %in% c("Missense_Mutation",
                                                      "Nonsense_Mutation",
                                                      "Nonstop_Mutation",
                                                      "Frame_Shift_Del",
                                                      "Frame_Shift_Ins"),]
# Select recurrent mutations (identified in at least two samples)
mut.id <- paste0(LGGmut$Chromosome,":",LGGmut$Start_Position,"-",LGGmut$End_Position,"|",LGGmut$Reference_Allele)
mut <- cbind(mut.id, LGGmut)
# Prepare selected mutations data for circos plot
s.mut <- mut[mut$mut.id %in% unique(mut.id[duplicated(mut.id)]),]
s.mut <- s.mut[,c("Chromosome","Start_Position","End_Position","Variant_Classification","Hugo_Symbol")]
s.mut <- unique(s.mut)
typeNames <- unique(s.mut[,4])
type <- c(4:1)
names(type) <- typeNames[1:4]
Type <- type[s.mut[,4]]
s.mut <- cbind(s.mut,Type)
s.mut <- s.mut[,c(1:3,6,4,5)]

# Load recurrent CNV data for selected cancer (e.g. "LGG")
load("GBM_CNV_results.rda")
# Prepare selected sample CNV data for circos plot
s.cnv <- as.data.frame(RecCNV[RecCNV[,"q-value"] <= threshold,c(1:4,6)])
s.cnv <- s.cnv[,c(1,3,4,2)]
s.cnv[s.cnv$Chromosome == 23,"Chromosome"] <- "X"
s.cnv[s.cnv$Chromosome == 24,"Chromosome"] <- "Y"
Chromosome <- paste0("chr",s.cnv[,1])
s.cnv <- cbind(Chromosome, s.cnv[,-1])
s.cnv[,1] <- as.character(s.cnv[,1])
s.cnv[,4] <- as.character(s.cnv[,4])
s.cnv <- cbind(s.cnv,CNV=1)
colnames(s.cnv) <- c("Chromosome","Start_position","End_position","Aberration_Kind","CNV")

library(circlize)
# Draw genomic circos plot
par(mar=c(1,1,1,1), cex=1)
circos.initializeWithIdeogram()
# Add CNV results
colors <- c("forestgreen","firebrick")
names(colors)  <- c(0,1)
circos.genomicTrackPlotRegion(s.cnv,  ylim = c(0,1.2),
                              panel.fun = function(region, value, ...) {
                                circos.genomicRect(region, value, ytop.column = 2, ybottom = 0,
                                                   col = colors[value[[1]]], 
                                                   border="white")
                                cell.xlim = get.cell.meta.data("cell.xlim")
                                circos.lines(cell.xlim, c(0, 0), lty = 2, col = "#00000040")
                              })
# Add mutation results
colors <- c("blue","green","red","gold")
names(colors)  <- typeNames[1:4]
circos.genomicTrackPlotRegion(s.mut, ylim = c(1.2,4.2),
                              panel.fun = function(region, value, ...) {
                                circos.genomicPoints(region, value, cex = 0.8, pch = 16, col = colors[value[[2]]], ...)
                              })

circos.clear()

legend(-0.2, 0.2, bty="n", y.intersp=1, c("Amp","Del"), pch=15, 
       col=c("firebrick","forestgreen"), title="CNVs", text.font=1, cex=0.4, title.adj=0)
legend(-0.2, 0, bty="n", y.intersp=1, names(colors), pch=16, 
       col=colors, title="Mutations", text.font=1, cex=0.4, title.adj=0)

par(mar=c(1,1,1,1),cex=1.5)
circos.par("start.degree" = 90, canvas.xlim = c(0, 1), canvas.ylim = c(0, 1), 
           gap.degree = 270, cell.padding = c(0, 0, 0, 0), track.margin = c(0.005, 0.005))
circos.initializeWithIdeogram(chromosome.index = "chr17")
circos.par(cell.padding = c(0, 0, 0, 0))
# Add CNV results
colors <- c("forestgreen","firebrick")
names(colors)  <- c(0,1)
circos.genomicTrackPlotRegion(s.cnv,  ylim = c(0,1.2),
                              panel.fun = function(region, value, ...) {
                                circos.genomicRect(region, value, ytop.column = 2, ybottom = 0,
                                                   col = colors[value[[1]]], 
                                                   border="white")
                                cell.xlim = get.cell.meta.data("cell.xlim")
                                circos.lines(cell.xlim, c(0, 0), lty = 2, col = "#00000040")
                              })

# Add mutation results representing single genes
genes.mut <- paste0(s.mut$Hugo_Symbol,"-",s.mut$Type)
s.mutt <- cbind(s.mut,genes.mut)
n.mut <- table(genes.mut)
idx <- !duplicated(s.mutt$genes.mut)
s.mutt <- s.mutt[idx,]
s.mutt <- cbind(s.mutt,num=n.mut[s.mutt$genes.mut])
genes.num <- paste0(s.mutt$Hugo_Symbol," (",s.mutt$num.Freq,")")
s.mutt <- cbind(s.mutt[,-c(6:8)],genes.num)
s.mutt[,6] <- as.character(s.mutt[,6])
s.mutt[,4] <- s.mutt[,4]/2
s.mutt$num.Freq <- NULL
colors <- c("blue","green","red","gold")
names(colors)  <- typeNames[1:4]
circos.genomicTrackPlotRegion(s.mutt, ylim = c(0.3,2.2), track.height = 0.05,
                              panel.fun = function(region, value, ...) {
                                circos.genomicPoints(region, value, cex = 0.4, pch = 16, col = colors[value[[2]]], ...)
                              })

circos.genomicTrackPlotRegion(s.mutt, ylim = c(0, 1), track.height = 0.1, bg.border = NA)
i_track = get.cell.meta.data("track.index")

circos.genomicTrackPlotRegion(s.mutt, ylim = c(0,1),
                              panel.fun = function(region, value, ...) {
                                circos.genomicText(region, value, 
                                                   y = 1, 
                                                   labels.column = 3,
                                                   col = colors[value[[2]]],
                                                   facing = "clockwise", adj = c(1, 0.5),
                                                   posTransform = posTransform.text, cex = 0.4, niceFacing = TRUE)
                              }, track.height = 0.1, bg.border = NA)

circos.genomicPosTransformLines(s.mutt,
                                posTransform = function(region, value)
                                  posTransform.text(region, 
                                                    y = 1, 
                                                    labels = value[[3]],
                                                    cex = 0.4, track.index = i_track+1),
                                direction = "inside", track.index = i_track)

circos.clear()

legend(0.25, 0.2, bty="n", y.intersp=1, c("Amp","Del"), pch=15, 
       col=c("firebrick","forestgreen"), title="CNVs", text.font=1, cex=0.4, title.adj=0)
legend(0, 0.2, bty="n", y.intersp=1, names(colors), pch=16, 
       col=colors, title="Mutations", text.font=1, cex=0.4, title.adj=0)

Pre-Processing Data

The LGG and GBM data used for following transcriptomic analysis were downloaded using TCGAbiolinks. We downloaded only primary solid tumor (TP) samples, which resulted in 516 LGG samples and 156 GBM samples, then prepared it in two separate rse objects (RangedSummarizedExperiment) saving them as an R object with a filename including both the cancer name and the name of the platform used for gene expression data (see listing below).

query <- GDCquery(project = "TCGA-GBM",
                  data.category = "Gene expression",
                  data.type = "Gene expression quantification",
                  platform = "Illumina HiSeq", 
                  file.type  = "results", 
                  sample.type = c("Primary solid Tumor"),
                  legacy = TRUE)
# We will use only 20 samples to make the example faster
query$results[[1]] <-  query$results[[1]][1:20,]                  
GDCdownload(query)
gbm.exp <- GDCprepare(query, 
                      save = TRUE, 
                      summarizedExperiment = TRUE, 
                      save.filename = "GBMIllumina_HiSeq.rda")

query <- GDCquery(project = "TCGA-LGG",
                  data.category = "Gene expression",
                  data.type = "Gene expression quantification",
                  platform = "Illumina HiSeq", 
                  file.type  = "results", 
                  sample.type = c("Primary solid Tumor"),
                  legacy = TRUE)
# We will use only 20 samples to make the example faster
query$results[[1]] <-  query$results[[1]][1:20,]
GDCdownload(query)
lgg.exp <- GDCprepare(query, 
                      save = TRUE, 
                      summarizedExperiment = TRUE, 
                      save.filename = "LGGIllumina_HiSeq.rda")

To pre-process the data, first, we searched for possible outliers using the TCGAanalyze_Preprocessing function, which performs an Array Array Intensity correlation AAIC (lines 14-17 and 26-29 of listing below). In this way, we defined a square symmetric matrix of pearson correlation among all samples in each cancer type (LGG or GBM). This matrix found 0 samples with low correlation (cor.cut = 0.6) that can be identified as possible outliers.

Second, using the TCGAanalyze_Normalization function, which encompasses the functions of the EDASeq package, we normalized mRNA transcripts.

This function implements Within-lane normalization procedures to adjust for GC-content effect (or other gene-level effects) on read counts: loess robust local regression, global-scaling, and full-quantile normalization (Risso et al. 2011) and between-lane normalization procedures to adjust for distributional differences between lanes (e.g., sequencing depth): global-scaling and full-quantile normalization (Bullard et al. 2010).

data("LGGIllumina_HiSeq")
data("GBMIllumina_HiSeq")

dataPrep_LGG <- TCGAanalyze_Preprocessing(object = lgg.exp,
                                          cor.cut = 0.6,    
                                          datatype = "raw_count",
                                          filename = "LGG_IlluminaHiSeq_RNASeqV2.png")

dataPrep_GBM <- TCGAanalyze_Preprocessing(object = gbm.exp,
                                          cor.cut = 0.6, 
                                          datatype = "raw_count",
                                          filename = "GBM_IlluminaHiSeq_RNASeqV2.png")

dataNorm <- TCGAanalyze_Normalization(tabDF = cbind(dataPrep_LGG, dataPrep_GBM),
                                      geneInfo = TCGAbiolinks::geneInfo,
                                      method = "gcContent") #18323   672

dataFilt <- TCGAanalyze_Filtering(tabDF = dataNorm,
                                  method = "quantile",
                                  qnt.cut =  0.25)  # 13742   672

save(dataFilt, file = paste0("LGG_GBM_Norm_IlluminaHiSeq.rda"))

dataFiltLGG <- subset(dataFilt, select = substr(colnames(dataFilt),1,12) %in% lgg_clin$bcr_patient_barcode)
dataFiltGBM <- subset(dataFilt, select = substr(colnames(dataFilt),1,12) %in% gbm_clin$bcr_patient_barcode)

dataDEGs <- TCGAanalyze_DEA(mat1 = dataFiltLGG,
                            mat2 = dataFiltGBM,
                            Cond1type = "LGG",
                            Cond2type = "GBM",
                            fdr.cut = 0.01 ,
                            logFC.cut = 1,
                            method = "glmLRT")

# Number of differentially expressed genes (DEG)
nrow(dataDEGs)

[1] 2535

EA: enrichment analysis

In order to understand the underlying biological process of DEGs we performed an enrichment analysis using TCGAanalyze_EA_complete function.

#-------------------  4.2 EA: enrichment analysis             --------------------
ansEA <- TCGAanalyze_EAcomplete(TFname="DEA genes LGG Vs GBM", RegulonList = rownames(dataDEGs))
TCGAvisualize_EAbarplot(tf = rownames(ansEA$ResBP),
                        filename = NULL,
                        GOBPTab = ansEA$ResBP, 
                        nRGTab = rownames(dataDEGs),
                        nBar = 20)

TCGAvisualize_EAbarplot(tf = rownames(ansEA$ResBP),
                        filename = NULL,
                        GOCCTab = ansEA$ResCC,
                        nRGTab = rownames(dataDEGs),
                        nBar = 20)

TCGAvisualize_EAbarplot(tf = rownames(ansEA$ResBP),
                        filename = NULL,
                        GOMFTab = ansEA$ResMF, 
                        nRGTab = rownames(dataDEGs),
                        nBar = 20)

TCGAvisualize_EAbarplot(tf = rownames(ansEA$ResBP),
                        filename = NULL,
                        PathTab = ansEA$ResPat,
                        nRGTab = rownames(dataDEGs),
                        nBar = 20)

The plot shows canonical pathways significantly overrepresented (enriched) by the DEGs (differentially expressed genes) with the number of genes for the main categories of three ontologies (GO:biological process, GO:cellular component, and GO:molecular function, respectively). The most statistically significant canonical pathways identified in DEGs list are listed according to their p-value corrected FDR (-Log) (colored bars) and the ratio of list genes found in each pathway over the total number of genes in that pathway (ratio, red line).]

TCGAanalyze_EAbarplot outputs a bar chart as shown in figure with the number of genes for the main categories of three ontologies (i.e., GO:biological process, GO:cellular component, and GO:molecular function).

The Figure shows canonical pathways significantly overrepresented (enriched) by the DEGs. The most statistically significant canonical pathways identified in the DEGs are ranked according to their p-value corrected FDR (-Log10) (colored bars) and the ratio of list genes found in each pathway over the total number of genes in that pathway (ratio, red line).

PEA: Pathways enrichment analysis

To verify if the genes found have a specific role in a pathway, the Bioconductor package pathview (Luo and Brouwer 2013) can be used. Listing below shows an example how to use it. It can receive, for example, a named vector of genes with their expression level, the pathway.id which can be found in KEGG database, the species (’hsa’ for Homo sapiens) and the limits for the gene expression.

library(SummarizedExperiment)
GenelistComplete <- rownames(assay(lgg.exp,1))

# DEGs TopTable
dataDEGsFiltLevel <- TCGAanalyze_LevelTab(dataDEGs,"LGG","GBM",
                                          dataFilt[,colnames(dataFiltLGG)],
                                          dataFilt[,colnames(dataFiltGBM)])

dataDEGsFiltLevel$GeneID <- 0

library(clusterProfiler)
# Converting Gene symbol to geneID
eg = as.data.frame(bitr(dataDEGsFiltLevel$mRNA,
                        fromType="SYMBOL",
                        toType="ENTREZID",
                        OrgDb="org.Hs.eg.db"))
eg <- eg[!duplicated(eg$SYMBOL),]

dataDEGsFiltLevel <- dataDEGsFiltLevel[dataDEGsFiltLevel$mRNA %in% eg$SYMBOL,]

dataDEGsFiltLevel <- dataDEGsFiltLevel[order(dataDEGsFiltLevel$mRNA,decreasing=FALSE),]
eg <- eg[order(eg$SYMBOL,decreasing=FALSE),]

# table(eg$SYMBOL == dataDEGsFiltLevel$mRNA) should be TRUE
all(eg$SYMBOL == dataDEGsFiltLevel$mRNA)

[1] TRUE

dataDEGsFiltLevel$GeneID <- eg$ENTREZID

dataDEGsFiltLevel_sub <- subset(dataDEGsFiltLevel, select = c("GeneID", "logFC"))
genelistDEGs <- as.numeric(dataDEGsFiltLevel_sub$logFC)
names(genelistDEGs) <- dataDEGsFiltLevel_sub$GeneID
library(pathview)
# pathway.id: hsa05214 is the glioma pathway
# limit: sets the limit for gene expression legend and color
hsa05214 <- pathview::pathview(gene.data  = genelistDEGs,
                               pathway.id = "hsa05214",
                               species    = "hsa",
                               limit = list(gene=as.integer(max(abs(genelistDEGs)))))

The red genes are up-regulated and the green genes are down-regulated in the LGG samples compared to GBM.

Pathways enrichment analysis: glioma pathway. Red defines genes that are up-regulated and green defines genes that are down-regulated.

Inference of gene regulatory networks

Starting with the set of differentially expressed genes, we infer gene regulatory networks using the following state-of-the-art inference algorithms: ARACNE (Margolin et al. 2006), CLR (Faith et al. 2007), MRNET (Meyer et al. 2007) and C3NET (Altay and Emmert-Streib 2010). These methods are based on mutual inference and use different heuristics to infer the edges of the network. These methods have been made available via Bioconductor/CRAN packages (MINET (Meyer, Lafitte, and Bontempi 2008) and c3net, (Altay and Emmert-Streib 2010) respectively). Many gene regulatory interactions have been experimentally validated and published. These ’known’ interactions can be accessed using different tools and databases such as BioGrid (Stark et al. 2006) or GeneMANIA (Montojo et al. 2010). However, this knowledge is far from complete and in most cases only contains a small subset of the real interactome. The quality assessment of the inferred networks can be carried out by comparing the inferred interactions to those that have been validated. This comparison results in a confusion matrix as presented in Table below.

Different quality measures can then be computed such as the false positive rate \[fpr=\frac{FP}{FP+TN},\] the true positive rate (also called recall) \[tpr=\frac{TP}{TP+FN}\] and the precision \[p=\frac{TP}{TP+FP}.\] The performance of an algorithm can then be summarized using ROC (false positive rate versus true positive rate) or PR (precision versus recall) curves.

A weakness of this type of comparison is that an edge that is not present in the set of known interactions can either mean that an experimental validation has been tried and did not show any regulatory mechanism or (more likely) has not yet been attempted.
In the following, we ran the nce on i) the 2,901 differentially expressed genes identified in Section “Transcriptomic analysis”.

Retrieving known interactions

We obtained a set of known interactions from the BioGrid database.

There are 3,941 unique interactions between the 2,901 differentially expressed genes.

Using differentially expressed genes from TCGAbiolinks workflow

We start this analysis by inferring one gene set for the LGG data.

### read biogrid info (available in TCGAWorkflowData as "biogrid")
### Check last version in https://thebiogrid.org/download.php 
file <- "http://thebiogrid.org/downloads/archives/Release%20Archive/BIOGRID-3.4.146/BIOGRID-ALL-3.4.146.tab2.zip"
if(!file.exists(gsub("zip","txt",basename(file)))){
  downloader::download(file,basename(file))
  unzip(basename(file),junkpaths =TRUE)
}

tmp.biogrid <- read.csv(gsub("zip","txt",basename(file)), header=TRUE, sep="\t", stringsAsFactors=FALSE)

### plot details (colors & symbols)
mycols <- c('#e41a1c','#377eb8','#4daf4a','#984ea3','#ff7f00','#ffff33','#a65628')

### load network inference libraries
library(minet)
library(c3net)

### deferentially identified genes using TCGAbiolinks
# we will use only a subset (first 50 genes) of it to make the example faster
names.genes.de <- rownames(dataDEGs)[1:30]

data(biogrid)
net.biogrid.de <- getAdjacencyBiogrid(tmp.biogrid, names.genes.de)

mydata <- dataFiltLGG[names.genes.de, ]

### infer networks
t.mydata <- t(mydata)
net.aracne <- minet(t.mydata, method = "aracne")
net.mrnet <- minet(t.mydata)
net.clr <- minet(t.mydata, method = "clr")
net.c3net <- c3net(mydata)

### validate compared to biogrid network
tmp.val <- list(validate(net.aracne, net.biogrid.de), 
                validate(net.mrnet, net.biogrid.de),
                validate(net.clr, net.biogrid.de), 
                validate(net.c3net, net.biogrid.de))

### plot roc and compute auc for the different networks
dev1 <- show.roc(tmp.val[[1]],cex=0.3,col=mycols[1],type="l")
res.auc <- auc.roc(tmp.val[[1]])
for(count in 2:length(tmp.val)){
  show.roc(tmp.val[[count]],device=dev1,cex=0.3,col=mycols[count],type="l")
  res.auc <- c(res.auc, auc.roc(tmp.val[[count]]))
}

legend("bottomright", legend=paste(c("aracne","mrnet","clr","c3net"), signif(res.auc,4), sep=": "),
       col=mycols[1:length(tmp.val)],lty=1, bty="n" )
# Please, uncomment this line to produce the pdf files.
# dev.copy2pdf(width=8,height=8,device = dev1, file = paste0("roc_biogrid_",cancertype,".pdf"))

ROC with corresponding AUC for inferred GBM networks compared to BioGrid interactions

In Figure above, the obtained ROC curve and the corresponding area under the curve (AUC) are presented. It can be observed that CLR and MRNET perform best when comparing the inferred network with known interactions from the BioGrid database.

Visualizing the mean DNA methylation of each patient

It should be highlighted that some pre-processing of the DNA methylation data was done. The DNA methylation data from the 450k platform has three types of probes cg (CpG loci) , ch (non-CpG loci) and rs (SNP assay). The last type of probe can be used for sample identification and tracking and should be excluded for differential methylation analysis according to the ilumina manual. Therefore, the rs probes were removed (see listing below lines 68). Also, probes in chromosomes X, Y were removed to eliminate potential artifacts originating from the presence of a different proportion of males and females (Marabita et al. 2013). The last pre-processing steps were to remove probes with at least one NA value (see listing below lines 65).

After this pre-processing step and using the function TCGAvisualize_meanMethylation function, we can look at the mean DNA methylation of each patient in each group. It receives as argument a SummarizedExperiment object with the DNA methylation data, and the arguments groupCol and subgroupCol which should be two columns from the sample information matrix of the SummarizedExperiment object (accessed by the colData function) (see listing below lines 70-74).

#----------------------------
# Obtaining DNA methylation
#----------------------------
# Samples
lgg.samples <- matchedMetExp("TCGA-LGG", n = 10)
gbm.samples <- matchedMetExp("TCGA-GBM", n = 10)
samples <- c(lgg.samples,gbm.samples)

#-----------------------------------
# 1 - Methylation
# ----------------------------------
# For methylation it is quicker in this case to download the tar.gz file
# and get the samples we want instead of downloading files by files
query <- GDCquery(project = c("TCGA-LGG","TCGA-GBM"),
                  data.category = "DNA methylation",
                  platform = "Illumina Human Methylation 450",
                  legacy = TRUE, 
                  barcode = samples)
GDCdownload(query)
met <- GDCprepare(query, save = FALSE)

# We will use only chr9 to make the example faster
met <- subset(met,subset = as.character(seqnames(met)) %in% c("chr9"))
# This data is avaliable in the package (object elmerExample)

data(elmerExample)
#----------------------------
# Mean methylation
#----------------------------
# Plot a barplot for the groups in the disease column in the
# summarizedExperiment object

# remove probes with NA (similar to na.omit)
met <- subset(met,subset = (rowSums(is.na(assay(met))) == 0))

TCGAvisualize_meanMethylation(met,
                              groupCol = "project_id",
                              group.legend  = "Groups",
                              filename = NULL,
                              print.pvalue = TRUE)

##     groups      Mean    Median       Max       Min
## 1 TCGA-GBM 0.4755915 0.4755397 0.5039685 0.4284987
## 2 TCGA-LGG 0.5271248 0.5283638 0.5661392 0.4441902
##             TCGA-GBM    TCGA-LGG
## TCGA-GBM          NA 0.001033341
## TCGA-LGG 0.001033341          NA

The figure above illustrates a mean DNA methylation plot for each sample in the GBM group (140 samples) and a mean DNA methylation for each sample in the LGG group. Genome-wide view of the data highlights a difference between the groups of tumors.

Searching for differentially methylated CpG sites

The next step is to define differentially methylated CpG sites between the two groups. This can be done using the TCGAanalyze_DMR function (see listing below). The DNA methylation data (level 3) is presented in the form of beta-values that uses a scale ranging from 0.0 (probes completely unmethylated ) up to 1.0 (probes completely methylated).

To find these differentially methylated CpG sites, first, the function calculates the difference between the mean DNA methylation (mean of the beta-values) of each group for each probe. Second, it tests for differential expression between two groups using the Wilcoxon test adjusting by the Benjamini-Hochberg method. Arguments of TCGAanalyze_DMR was set to require a minimum absolute beta-values difference of 0.15 and an adjusted p-value of less than \(0.05\).

After these tests, a volcano plot (x-axis: difference of mean DNA methylation, y-axis: statistical significance) is created to help users identify the differentially methylated CpG sites and return the object with the results in the rowRanges.

#------- Searching for differentially methylated CpG sites     ----------
met <- TCGAanalyze_DMR(met,
                       groupCol = "project_id", # a column in the colData matrix
                       group1 = "TCGA-GBM", # a type of the disease type column
                       group2 = "TCGA-LGG", # a type of the disease column
                       p.cut = 0.05,
                       diffmean.cut = 0.15,
                       save = FALSE,
                       legend = "State",
                       plot.filename = "LGG_GBM_metvolcano.png",
                       cores = 1 # if set to 1 there will be a progress bar
)

The figure below shows the volcano plot produced by listing below. This plot aids the user in selecting relevant thresholds, as we search for candidate biological DMRs.

Volcano plot: searching for differentially methylated CpG sites (x-axis:difference of mean DNA methylation, y-axis: statistical significance)

To visualize the level of DNA methylation of these probes across all samples, we use heatmaps that can be generated by the Bioconductor package complexHeatmap (Z., n.d.). To create a heatmap using the complexHeatmap package, the user should provide at least one matrix with the DNA methylation levels. Also, annotation layers can be added and placed at the bottom, top, left side and right side of the heatmap to provide additional metadata description. The listing below shows the code to produce the heatmap of a DNA methylation datum.

#--------------------------
# DNA Methylation heatmap
#-------------------------
library(ComplexHeatmap)
clinical <- plyr::rbind.fill(gbm_clin,lgg_clin)

# get the probes that are Hypermethylated or Hypomethylated
# met is the same object of the section 'DNA methylation analysis'
status.col <- "status.TCGA.GBM.TCGA.LGG"
sig.met <- met[values(met)[,status.col] %in% c("Hypermethylated","Hypomethylated"),]


# top annotation, which sampples are LGG and GBM
# We will add clinical data as annotation of the samples
# we will sort the clinical data to have the same order of the DNA methylation matrix
clinical.order <- clinical[match(substr(colnames(sig.met),1,12),clinical$bcr_patient_barcode),]
ta = HeatmapAnnotation(df = clinical.order[,c("disease","gender","vital_status","race")],
                       col = list(disease = c("LGG" = "grey", "GBM" = "black"),
                                  gender = c("male" = "blue","female" = "pink")))

# row annotation: add the status for LGG in relation to GBM
# For exmaple: status.gbm.lgg Hypomethyated means that the
# mean DNA methylation of probes for lgg are hypomethylated
# compared to GBM ones.
ra = rowAnnotation(df = values(sig.met)[status.col],
                   col = list("status.TCGA.GBM.TCGA.LGG" = 
                                c("Hypomethylated" = "orange",
                                  "Hypermethylated" = "darkgreen")),
                   width = unit(1, "cm"))

heatmap  <- Heatmap(assay(sig.met),
                    name = "DNA methylation",
                    col = matlab::jet.colors(200),
                    show_row_names = FALSE,
                    cluster_rows = TRUE,
                    cluster_columns = FALSE,
                    show_column_names = FALSE,
                    bottom_annotation = ta,
                    column_title = "DNA Methylation") 
# Save to pdf
png("heatmap.png",width = 600, height = 400)
draw(heatmap, annotation_legend_side =  "bottom")
dev.off()

Motif analysis

Motif discovery is the attempt to extract small sequence signals hidden within largely non-functional intergenic sequences. These small sequence nucleotide signals (6-15 bp) might have a biological significance as they can be used to control the expression of genes. These sequences are called Regulatory motifs. The Bioconductor package rGADEM (Droit et al. 2015; Li 2009) provides an efficient de novo motif discovery algorithm for large-scale genomic sequence data.

The user may be interested in looking for unique signatures in the regions defined by ‘differentially methylated’ to identify candidate transcription factors that could bind to these elements affected by the accumulation or absence of DNA methylation. For this analysis we use a sequence of 100 bases before and after the probe location (See lines 6-8 in the listing below). An object will be returned which contains all relevant information about your motif analysis (i.e., sequence consensus, PWM, chromosome, p-value, etc).

Using Bioconductor package motifStack (Ou et al. 2013) it is possible to generate a graphic representation of multiple motifs with different similarity scores.

library(rGADEM)
library(BSgenome.Hsapiens.UCSC.hg19)
library(motifStack)

probes <- rowRanges(met)[values(met)[,status.col] %in% c("Hypermethylated" ,"Hypomethylated") ,]
# Get hypo/hyper methylated probes and make a 200bp window 
# surrounding each probe.
sequence <- RangedData(space = as.character(seqnames(probes)),
                       IRanges(start = start(ranges(probes)) - 100,
                               end = start(ranges(probes)) + 100), strand = "*")
#look for motifs
gadem <- GADEM(sequence, verbose = FALSE, genome = Hsapiens)

top 3 4, 5-mers: 20 40 60 top 3 4, 5-mers: 20 40 60 top 3 4, 5-mers: 20 40 60 top 3 4, 5-mers: 20 40 60

# How many motifs were found?
nMotifs(gadem)

[1] 3

# get the number of occurrences
nOccurrences(gadem)

[1] 268 183 137

# view all sequences consensus
consensus(gadem)

[1] “nGsnGGGGsnGGrGssnGGGs” “nAAAAAnrArAn” “nCCCAGGsmn”

# Print motif
pwm <- getPWM(gadem)
pfm  <- new("pfm",mat=pwm[[1]],name="Novel Site 1")
plotMotifLogo(pfm)

# Number of instances of motif 1?
length(gadem@motifList[[1]]@alignList)

[1] 268

After rGADEM returns it’s results, the user can use MotIV package (Mercier and Gottardo 2014; Mahony, Auron, and Benos 2007; Mahony and Benos 2007; Mercier et al. 2011) to start the motif matching analysis (line 4 of listing below).

library(MotIV)
# motif matching analysis
analysis.jaspar <- motifMatch(pwm)

Ungapped Alignment
Scores read
Database read
Motif matches : 5

summary(analysis.jaspar)

Number of input motifs :  3 
Input motifs names :  GGGCCACCGGGGGGCCCGGGT GAAAAATGAGTT ACCCAGGCAG 
Number of matches per motif:  5 
Matches repartition : 
TFAP2A       EBF1 EWSR1-FLI1      FOXI1      Foxd3      Foxq1       IRF1 
     2          1          1          1          1          1          1 
 PLAG1       Pax4      RREB1      SOX10        SP1     Zfp423        Zfx 
     1          1          1          1          1          1          1 
Arguments used : 
  -metric name :  PCC 
  -alignment :  SWU 

plot(analysis.jaspar, ncol=2, top=5, rev=FALSE, main="", bysim=TRUE, cex=0.4)

#  visualize the quality of the results looking into the alignments
# E-value gives an estimation of the match.
alignment <- viewAlignments(analysis.jaspar)
print(alignment[[1]])

   SP1                     Pax4                            

seq “NGSNGGGGSNGGRGSSNGGGN” “NGSNGGGGSNGGRGSSNGGGN———” match “—GGGGGNGGGG——–” “GKRNNNNKNNNNNNNNNNNKNNWWWTWTTY” evalue “6.6989e-08” “3.3633e-05”
TFAP2A RREB1
seq “NGSNGGGGSNGGRGSSNGGGN” “NGSNGGGGSNGGRGSSNGGGN” match “——GCCNNNNGS——” “NGKKKGKGGGTGKTTTGGGG-” evalue “4.0753e-04” “8.4676e-04”
EWSR1-FLI1
seq “NGSNGGGGSNGGRGSSNGGGN” match “-GGAAGGAAGGAAGGAAGG–” evalue “1.553e-03”

Integration of DNA methylation & gene expression data

After finding differentially methylated CpG sites, one might ask whether nearby genes also undergo a change in expression either an increase or a decrease. DNA methylation at promoters of genes has been shown to be associated with silencing of the respective gene. The starburst plot is proposed to combine information from two volcano plots and is applicable for studies of DNA methylation and gene expression (Noushmehr et al. 2010). Even though, being desirable that both gene expression and DNA methylation data are from the same samples, the starburst plot can be applied as a meta-analysis tool, combining data from different samples (Siegmund 2011). We used the TCGAvisualize_starburst function to create a starburst plot. The \(log_{10}\) (FDR-corrected p-value) for DNA methylation is plotted on the x-axis, and for gene expression on the y-axis, for each gene. The horizontal black dashed line shows the FDR-adjusted p-value of \(10^{-2}\) for the expression data and the vertical ones shows the FDR-adjusted P value of \(10^{-2}\) for the DNA methylation data. While the argument met.p.cut and exp.p.cut control the black dashed lines, the arguments diffmean.cut and logFC.cut will be used to highlight the genes that respect these parameters (circled genes in Figure below). For the example below we set higher p.cuts trying to get the most significant list of pair gene/probes. But for the next sections we will use exp.p.cut = 0.01 and logFC.cut = 1 as the previous sections.

#------------------- Starburst plot ------------------------------
starburst <- TCGAvisualize_starburst(met,    # DNA methylation with results
                                     dataDEGs,    # DEG results
                                     group1 = "Glioblastoma Multiforme", 
                                     group2 = "Brain Lower Grade Glioma", 
                                     filename = "starburst.png",
                                     genome = "hg19",
                                     met.platform = "450K",
                                     met.p.cut = 0.05,
                                     exp.p.cut = 10^-2,
                                     diffmean.cut = 0.15,
                                     logFC.cut = 1,
                                     width = 15,height = 10,
                                     names = TRUE)

Due to the time constraints, with the reduced number of samples, genes and probes, the result showed previously might not be very descriptive. A result with a large dataset can be found in TCGAbiolinks vignette.

ChIP-seq analysis

ChIP-seq is used primarily to determine how transcription factors and other chromatin-associated proteins influence phenotype-affecting mechanisms. Determining how proteins interact with DNA to regulate gene expression is essential for fully understanding many biological processes and disease states. The aim is to explore significant overlap datasets for inferring co-regulation or transcription factor complex for further investigation. A summary of the association of each histone mark is shown in the table below.

Besides, ChIP-seq data exists in the ROADMAP database and can be obtained through the AnnotationHub package (T. D. Morgan M Carlson M and S., n.d.) or from Roadmap web portal. The table below shows the description of all the roadmap files that are available through AnnotationHub.

After obtaining the ChIP-seq data, we can then identify overlapping regions with the regions identified in the starburst plot. The narrowPeak files are the ones selected for this step. For a complete pipeline with Chip-seq data, Bioconductor provides excellent tutorials to work with ChIP-seq and we encourage our readers to review the following article (Aleksandra Pekowska 2015). The first step is shown in listing below is to download the chip-seq data. The function query received as argument the annotationHub database (ah) and a list of keywords to be used for searching the data, EpigenomeRoadmap is selecting the roadmap database, consolidated is selecting only the consolidate epigenomes, brain is selecting the brain samples, E068 is one of the epigenomes for the brain (keywords can be seen in the summary table) and narrowPeak is selecting the type of file. The data downloaded is a processed data from an integrative Analysis of 111 reference human epigenomes (Kundaje et al. 2015).

library(ChIPseeker)
library(pbapply)
library(ggplot2)

#------------------ Working with ChipSeq data ---------------
# Step 1: download histone marks for a brain and non-brain samples.
#------------------------------------------------------------
# loading annotation hub database
library(AnnotationHub)
ah = AnnotationHub()

# Searching for brain consolidated epigenomes in the roadmap database
bpChipEpi_brain <- query(ah , c("EpigenomeRoadMap", "narrowPeak", "chip", "consolidated","brain","E068"))
# Get chip-seq data
histone.marks <- pblapply(names(bpChipEpi_brain), function(x) {ah[[x]]})
names(histone.marks) <- names(bpChipEpi_brain) 
# OBS: histone.marks is available in TCGAWorkflowData package

The Chipseeker package (Yu, Wang, and He 2015) implements functions that use Chip-seq data to retrieve the nearest genes around the peak, to annotate genomic region of the peak, among others. Also, it provides several visualization functions to summarize the coverage of the peak, average profile and heatmap of peaks binding to TSS regions, genomic annotation, distance to TSS and overlap of peaks or genes.

After downloading the histone marks, it is useful to verify the average profile of peaks binding to hypomethylated and hypermethylated regions, which will help the user understand better the regions found. Listing below shows an example of code to plot the average profile.

To help the user better understand the regions found in the DMR analysis, we downloaded histone marks specific to brain tissue using the AnnotationHub package that can access the Roadmap database. Next, the Chipseeker was used to visualize how histone modifications are enriched onto hypomethylated and hypermethylated regions, (listing below). The enrichment heatmap and the average profile of peaks binding to those regions.

data(histoneMarks)
# Create a GR object based on the hypo/hypermethylated probes.
probes <- keepStandardChromosomes(rowRanges(met)[values(met)[,status.col] %in% c("Hypermethylated","Hypomethylated"),])
# Defining a window of 3kbp - 3kbp_probe_3kbp
ranges(probes) <- IRanges(start = c(start(ranges(probes)) - 3000), end = c(start(ranges(probes)) + 3000))

### Profile of ChIP peaks binding to TSS regions
# First of all, to calculate the profile of ChIP peaks binding to TSS regions, we should
# prepare the TSS regions, which are defined as the flanking sequence of the TSS sites.
# Then align the peaks that are mapping to these regions and generate the tagMatrix.
tagMatrixList <- pbapply::pblapply(histone.marks, function(x) {
  getTagMatrix(keepStandardChromosomes(x), windows = probes, weightCol = "score")
})
# change names retrieved with the following command: basename(bpChipEpi_brain$title)
names(tagMatrixList) <- c("H3K4me1","H3K4me3", "H3K9ac", "H3K9me3", "H3K27ac",  "H3K27me3", "H3K36me3")

To plot the enrichment heatmap use the function tagHeatmap

tagHeatmap(tagMatrixList, xlim=c(-3000, 3000),color = NULL)

To plot the average profile of peaks binding to those region use plotAvgProf:

p <- plotAvgProf(tagMatrixList, xlim = c(-3000,3000), xlab = "Genomic Region (5'->3', centered on CpG)")
# We are centreing in the CpG instead of the TSS. So we'll change the labels manually
p <- p + scale_x_continuous(breaks=c(-3000,-1500,0,1500,3000),labels=c(-3000,-1500,"CpG",1500,3000))
library(ggthemes)
p + theme_few() + scale_colour_few(name="Histone marks") +  guides(colour = guide_legend(override.aes = list(size=4)))

The hypomethylated and hypermethylated regions are enriched for H3K4me3, H3K9ac, H3K27ac, and H3K4me1 which indicates regions of enhancers, promoters and increased activation of genomic elements. However, these regions are associated neither with transcribed regions nor Polycomb repression as the H3K36me3 and H3K27me3 heatmaps do not show an enrichment nearby the position 0, and the average profile also does not show a peak at position 0.

Supervised and unsupervised analysis mode

ELMER is designed to identify differences between two sets of samples within a given dataset. But, before we start our analysis, it is better to understand the two analysis mode available in ELMER 2.0 (supervised and unsupervised mode) which affects how the samples are selected for each comparison.

In original version (L Yao et al. 2015), the first step - identification of differentially methylated CpG probes (DMCs) - was hard-coded to identify DMCs between non-cancer vs. cancer samples, and the subsequent step was unsupervised, identifying changes within any subset of tumors. In ELMER 2.0, we generalize these strategies so that they are applicable to any paired dataset, including disease vs. healthy tissue for any disease type.

In the Unsupervised mode, as in ELMER 1.0, it is assumed that at least one of our groups is heterogenous (i.e Breast Invasive Carcinoma primary solid tumors samples might belong to basal,Luminal A, Her2 molecular subtypes), for that reason for each comparison samples from two most extreme quintiles based on DNA methylation levels were used, generating both the \(M\) (methylated) and the \(U\) (unmethylated) groups. In the Supervised mode, it is assumed that our groups are homogenous, for that reason the \(U\) and \(M\) groups are defined strictly by sample group labels, and all samples in each group are used (Silva et al. 2018).

The next steps we will perform ELMER analysis using the supervised mode, which will use all samples from GBM and compare to all LGG samples in all analysis.

Identification of differentially methylated CpGs (DMCs)

For each distal probe, samples of each group (group 1 and group 2) are ranked by their DNA methylation beta values, those samples in the lower quintile (20% samples with the lowest methylation levels) of each group are used to identify if the probe is hypomethylated in group 1 compared to group 2, using an unpaired one-tailed t-test. The 20% is a parameter to the diff.meth function called minSubgroupFrac. For the (ungrouped) cancer case, this is set to 20% as in (L Yao et al. 2015), because we typically wanted to be able to detect a specific molecular subtype among the tumor samples; these subtypes often make up only a minority of samples, and 20% was chosen as a lower bound for the purposes of statistical power (high enough sample numbers to yield t-test p-values that could overcome multiple hypothesis corrections, yet low enough to be able to capture changes in individual molecular subtypes occurring in 20% or more of the cases.) This number can be set arbitrarily as an input to the diff.meth function and should be tuned based on sample sizes in individual studies.

As we are dealing with a low number of samples, and we are comparing two cancer types the minSubgroupFrac parameter is set to 100% to use all samples in this analysis.

#--------------------------------------
# STEP 3: Analysis                     |
#--------------------------------------
# Step 3.1: Get diff methylated probes |
#--------------------------------------
message("Get diff methylated probes")
Sig.probes <- get.diff.meth(data = mae, 
                            group.col = group.col,
                            group1 = group1,
                            group2 =  group2,
                            minSubgroupFrac = 1.0, # Use all samples
                            sig.dif = 0.2, # defualt is 0.3
                            diff.dir = direction, # Search for hypomethylated probes in group 1
                            cores = cores, 
                            dir.out = dir.out, 
                            pvalue = 0.1)

datatable(Sig.probes[1:10,], 
          options = list(scrollX = TRUE, keys = TRUE, pageLength = 5), 
          rownames = TRUE)

Identification of putative target gene(s) for differentially methylated distal probes

For each differentially methylated distal probe (DMC), the closest 10 upstream genes and the closest 10 downstream genes are tested for inverse correlation between methylation of the probe and expression of the gene (the number 10 can be changed using the numFlankingGenes parameter). To select these genes, the probe-gene distance is defined as the distance from the probe to the transcription start site specified by the ENSEMBL gene level annotations (Yates et al. 2015) accessed via the R/Bioconductor package biomaRt (Durinck et al. 2009,Durinck et al. (2005)). By choosing a constant number of genes to test for each probe, our goal is to avoid systematic false positives for probes in gene-rich regions. This is especially important given the highly non-uniform gene density of mammalian genomes. Thus, exactly 20 statistical tests were performed for each probe, as follows.

For each probe-gene pair, the samples (all samples from both groups) are divided into two groups: the M group, which consisted of the upper methylation quintile (the 20% of samples with the highest methylation at the enhancer probe), and the U group, which consists of the lowest methylation quintile (the 20% of samples with the lowest methylation.) The 20% ile cutoff is a configurable parameter minSubgroupFrac in the get.pair function. As with its usage in the diff.meth function, the default value of 20% is a balance, allowing for the identification of changes in a molecular subtype making up a minority (i.e. 20%) of cases, while also yielding enough statistical power to make strong predictions. For larger sample sizes or other experimental designs, this could be set even lower.

For each candidate probe-gene pair, the Mann-Whitney U test is used to test the null hypothesis that overall gene expression in group M is greater than or equal than that in group U. This non-parametric test was used in order to minimize the effects of expression outliers, which can occur across a very wide dynamic range. For each probe-gene pair tested, the raw p-value \(P_r\) is corrected for multiple hypothesis using a permutation approach as follows. The gene in the pair is held constant, and x random methylation probes are chosen to perform the same one-tailed U test, generating a set of x permutation p-values \(P_p\). We chose the x random probes only from among those that were “distal” (farther than \(2kb\) from an annotated transcription start site), in order to draw these null-model probes from the same set as the probe being tested (???). An empirical p-value \(P_e\) value was calculated using the following formula (which introduces a pseudo-count of 1):

Notice that in the Supervised mode, no additional filtering is necessary to ensure that the M and U group segregate by sample group labels. The two sample groups are segregated by definition, since these probes were selected for their differential methylation, with the same directionality, between the two groups.

In our example, we will reduce the number of closest genes evaluated to 5 (the 2 closest upstream genes and the closest 2 downstream genes). The mode is set to “supervised” which will set group 1 (GBM samples) as the \(U\) group, and the group 2 (LGG samples) as \(M\) group. Also, to make the example faster, the number of permutation is reduced to 5 (permu.size = 5), the raw p-value cut-off is loosened to 0.1 and the empirical p-value to 0.5, but it is recommended in a real analysis to use the default values.

#-------------------------------------------------------------
# Step 3.2: Identify significant probe-gene pairs            |
#-------------------------------------------------------------
# Collect nearby 20 genes for Sig.probes
message("Get nearby genes")
nearGenes <- GetNearGenes(data = mae,
                          numFlankingGenes = 4, # default is 20 genes
                          probes = Sig.probes$probe)

message("Get anti correlated probes-genes")
pair <- get.pair(data = mae,
                 group.col = group.col,
                 group1 = group1,
                 group2 =  group2,
                 nearGenes = nearGenes,
                 mode = "supervised",
                 minSubgroupFrac = 1, # % of samples to use in to create groups U/M
                 permu.dir = paste0(dir.out,"/permu"),
                 permu.size = 5, # Please set to 100000 to get significant results
                 raw.pvalue = 0.1,   # defualt is 0.001
                 Pe = 0.5, # Please set to 0.001 to get significant results
                 filter.probes = TRUE, # See preAssociationProbeFiltering function
                 filter.percentage = 0.05,
                 filter.portion = 0.3,
                 dir.out = dir.out,
                 diff.dir = direction,
                 cores = cores,
                 label = direction)

head(nearGenes)

datatable(pair[1:10,], 
          options = list(scrollX = TRUE, keys = TRUE, pageLength = 5), 
          rownames = TRUE)

Motif enrichment analysis

To identify enriched motifs and potential upstream regulatory TFs, HOCOMOCO v11 (Kulakovskiy et al. 2016) TF binding models were used as input for HOMER (Heinz et al. 2010) to find motif occurrences in a \(\pm 250bp\) region around each probe from EPIC and HM450 arrays. Transcription factor (TF) binding models are available at http://hocomoco.autosome.ru/downloads (using the HOMER specific format with threshold score levels corresponding to p-value \(>1^{-4}\)).

A motif enrichment analysis using Fisher’s exact test and Benjamini-Hochberg multiple hypothesis testing correction is performed using as correction background all distal probes.

A probe set was considered significantly enriched for a particular motif if the 95% confidence interval of the Odds Ratio was greater than \(1.1\) (specified by option lower.OR, \(1.1\) is default), and the motif occurred at least 10 times (specified by option min.incidence, \(10\) is default) in the probe set.

#-------------------------------------------------------------
# Step 3.3: Motif enrichment analysis on the selected probes |
#-------------------------------------------------------------
enriched.motif <- get.enriched.motif(data = mae,
                                     probes = pair$Probe, 
                                     dir.out = dir.out, 
                                     label = direction,
                                     pvalue = 1, # default is FDR < 0.05
                                     min.incidence = 10,
                                     lower.OR = 1.1)
# One of the output from the  previous function is a file with the motif, OR and Number of probes
# It will be used for plotting purposes
motif.enrichment <- read.csv(paste0(dir.out,"/getMotif.",direction, ".motif.enrichment.csv"))

head(enriched.motif[names(enriched.motif)[1]]) ## probes in the given set that have the first motif.

## $ZN524_HUMAN.H11MO.0.D
##  [1] "cg13787850" "cg14159539" "cg13857678" "cg13493766" "cg02602455"
##  [6] "cg14560364" "cg13562069" "cg26787801" "cg14262955" "cg14121282"
## [11] "cg13663116" "cg21210721" "cg12819393" "cg14186948" "cg14223856"
## [16] "cg00505502" "cg13753488" "cg10962166" "cg13595495" "cg14567085"
## [21] "cg26756625" "cg04681525" "cg14061378" "cg19264651" "cg00938819"
## [26] "cg14103872" "cg03916189" "cg13512951" "cg06752260"

datatable(motif.enrichment, 
          options = list(scrollX = TRUE, keys = TRUE, pageLength = 5), 
          filter = 'top',
          rownames = FALSE)

Identification of master regulator TFs

When a group of enhancers is coordinately altered in a specific sample subset, this is often the result of an altered upstream master regulator transcription factor in the gene regulatory network. ELMER tries to identify such transcription factors corresponding to each of the TF binding motifs enriched from the previous analysis step. For each enriched motif, ELMER takes the average DNA methylation of all distal probes (in significant probe-gene pairs) that contain that motif occurrence (within a \(\pm 250bp\) region) and compares this average DNA methylation to the expression of each gene annotated as a human TF.

In the Unsupervised mode, a statistical test is performed for each motif-TF pair, as follows. All samples are divided into two groups: the \(M\) group, which consists of the 20% of samples with the highest average methylation at all motif-adjacent probes, and the \(U\) group, which consisted of the 20% of samples with the lowest methylation. This step is performed by the get.TFs function, which takes minSubgroupFrac as an input parameter, again with a default of 20%. For each candidate motif-TF pair, the Mann-Whitney U test is used to test the null hypothesis that overall gene expression in group \(M\) is greater or equal than that in group \(U\). This non-parametric test was used in order to minimize the effects of expression outliers, which can occur across a very wide dynamic range. For each motif tested, this results in a raw p-value (\(P_r\)) for each of the human TFs.

All TFs are ranked by their \(-log_{10}(Pr)\) values, and those falling within the top 5% of this ranking were considered candidate upstream regulators. The best upstream TFs which are known to recognize to specific binding motif are automatically extracted as putative regulatory TFs, and rank ordered plots are created to visually inspect these relationships, as shown in the example below. Because the same motif can be recognized by many transcription factors of the same binding domain family, we define these relationships at both the family and subfamily classification level using the classifications from TFClass database (Wingender, Schoeps, and Dönitz 2013).

The Supervised mode uses the same approach as described for the identification of putative target gene(s) step. The \(U\) and \(M\) groups are one of the the label group of samples and the minSubgroupFrac parameter is set to 100% to use all samples from both groups in the statistical test.

#-------------------------------------------------------------
# Step 3.4: Identifying regulatory TFs                        |
#-------------------------------------------------------------
TF <- get.TFs(data = mae, 
              group.col = group.col,
              group1 = group1,
              group2 =  group2,
              mode = "supervised",
              enriched.motif = enriched.motif,
              dir.out = dir.out, 
              cores = cores,
              save.plots = FALSE,
              diff.dir = direction,
              label = direction)
    
# One of the output from the previous function is a file with the raking of TF,
# for each motif. It will be used for plotting purposes
TF.meth.cor <- get(load(paste0(dir.out,"/getTF.",direction,".TFs.with.motif.pvalue.rda")))

The output of this step is a data frame with the following columns:

• motif: enriched motif name.
• top_5percent_TFs: the top 5% TFs ranked
• potential.TFs.family: TF from the “top_5percent_TFs” that belongs to the same family as the TF of the motif,
  
• top.potential.TFs.family is the highest ranked TF belonging to the same family as the TF of the motif (same as the first TF from potential.TFs.family column).
• potential.TFs.subfamily and top.potential.TFs.subfamily are the same as potential.TFs.family and top.potential.TFs.family but considering the subfamily classification instead (you can check the classification at http://hocomoco11.autosome.ru/human/mono?full=true).

datatable(TF, 
          options = list(scrollX = TRUE, keys = TRUE, pageLength = 5), 
          rownames = FALSE)

datatable(TF.meth.cor[1:10,1:6], 
          options = list(scrollX = TRUE, keys = TRUE, pageLength = 5), 
          rownames = TRUE)

ELMER visualization functions

The pairs identified by ELMER can be visualized through a heatmap using the function heatmapPairs. The first heatmap shows the DNA methylation levels, while the second heatmap shows the gene expression levels of the target genes. The last heatmap shows the distance between the probe and the gene TSS.

heatmapPairs(data = mae, 
             group.col = group.col,
             group1 = group1, 
             group2 = group2, 
             annotation.col = c("gender"),
             pairs = pair,
             filename =  NULL)

When ELMER Identifies the enriched motifs for the distal enhancer probes which are significantly differentially methylated and linked to putative target gene, it will plot the Odds Ratio (x-axis) for the each motif found.

The list of enriched motifs for the hypo direction (probes hypomethylated in GBM group compared to the LGG group) with lower OR >= 1.1 is found in the Figure below.

motif.enrichment.plot(motif.enrichment = motif.enrichment, 
                      save = FALSE,
                      significant = list(lowerOR = 1.1)) # Filter motifs in the plot lowerOR > 1.3

After finding the enriched motifs, ELMER identifies regulatory transcription factors (TFs) whose expression is associated with DNA methylation at motifs. ELMER automatically creates a TF ranking plot for each enriched motif. This plot shows the TF ranking plots based on the association score \((-log(p-value))\) between TF expression and DNA methylation of the motif. This plot represents the rank of p-values of all TF, the top 3, the TF classified in the same family and sub-family according to TFClass database are highlighted. The dashed line represents the top 5% of all TF.

grid:TF.rank.plot(motif.pvalue=TF.meth.cor, motif=TF$motif[1], save=FALSE)

|=================================================|100% ~0 s remaining
|=================================================|100% Completed after 0 s

Also, for each motif, we can take a look at the three most relevant TFs. For example, Figure below shows the average DNA methylation level of sites with the first motif plotted against the expression of some top TFs associated with it. We can see that the GBM samples have a lower average methylation level of sites with the motif plotted and a higher average expression of the TFs.

png("TF.png",width = 800, height = 400)
scatter.plot(mae, 
             category = group.col, 
             save = FALSE, 
             lm_line = TRUE,
             byTF=list(TF=unlist(stringr::str_split(TF[1,"top_5percent_TFs"],";"))[1:4], 
                       probe=enriched.motif[[TF$motif[1]]]))
dev.off()

png 2