Abstract

In order to make light of cancer development, it is crucial to understand which genes play a role in the mechanisms linked to this disease and moreover which role that is. Commonly biological processes such as proliferation and apoptosis have been linked to cancer progression. Based on expression data we perform functional enrichment analysis, infer gene regulatory networks and upstream regulator analysis to score the importance of well-known biological processes with respect to the studied cancer. We then use these scores to predict two specific roles: genes that act as tumor suppressor genes (TSGs) and genes that act as oncogenes (OCGs). This methodology not only allows us to identify genes with dual role (TSG in one cancer type and OCG in another) but also to elucidate the underlying biological processes.

Introduction

Cancer development is influenced by mutations in two distinctly different categories of genes, known as tumor suppressor genes (TSG) and oncogenes (OCG). The occurrence of mutations in genes of the first category leads to faster cell proliferation while mutations in genes of second category increases or changes their function. We propose MoonlightR a new approach to define TSGs and OCGs based on functional enrichment analysis, infer gene regulatory networks and upstream regulator analysis to score the importance of well-known biological processes with respect to the studied cancer.

Moonlight’s pipeline

Moonlight’s pipeline is shown below:

Moonlight’s proposed workflow

The proposed pipeline consists of following eight steps:

1. getDataTCGA & getDataGEO for Data collection: expression levels of genes in all samples obtained with IlluminaHiSeq RNASeqV2 in 18 normal tissues (NT) and 18 cancer tissues (CT) according to TCGA criteria, and GEO data set matched to one of the 18 given TCGA cancer types as described in following Table TCGA / GEO.
2. DPA Differential Phenotype Analysis (DEA) to identify genes or probes that are different significantly with two phenotypes such as normal and tumor, or normal and stageI, normal and molecular subtype.
3. FEA Functional Enrichment Analysis (EA), using Fisher’s test, to identify gene sets (with biological functions linked to cancer¹) significantly enriched by RG.
4. GRN Gene regulatory network inferred between each single DEG (sDEG) and all genes by means of mutual information, obtaining for each DEG a list of regulated genes (RG).
5. URA Upstream Regulator Analysis for DEGs in each enriched gene set, we applied z-score being the ratio between the sum of all predicted effects for all the gene involved in the specific function and the square-root of the number of all genes.
6. PRA Pattern recognition analysis identifies candidate TCGs (down) and OCGs (up). We either use user defined biological processes or random forests.
7. We applied the above procedure to multiple cancer types to obtain cancer-specific lists of TCGs and OCGs. We compared the lists for each cancer: if a sDEG was TSG in a cancer and OCG in another we defined it as dual-role TSG-OCG. Otherwise if we found a sDEG defined as OCG or TSG only in one tissue we defined it tissue specific biomarker.
8. We use the COSMIC database to define a list of gold standard TSG and OCGs to assess the accuracy of the proposed method.

¹ For the devel version of MoonlightR we use a short extract of ten biological functions from QIAGEN’S Ingenuity Pathway Analysis (IPA). We are still working to integrate the package.

Installation

To install use the code below.

source("https://bioconductor.org/biocLite.R")
biocLite("MoonlightR")

Download: Get TCGA data

You can search TCGA data using the getDataTCGA function.

dataFilt <- getDataTCGA(cancerType = "LUAD", 
                          dataType = "Gene expression",
                          directory = "data",
                          nSample = 4)

dataFilt <- getDataTCGA(cancerType = "BRCA", 
                          dataType = "Methylation", 
                          directory = "data",nSample = 4)

Download: Get GEO data

You can search GEO data using the getDataGEO function.

GEO_TCGAtab: a 18x12 matrix that provides the GEO data set we matched to one of the 18 given TCGA cancer types

knitr::kable(GEO_TCGAtab, digits = 2, 
             caption = "Table with GEO data set matched to one 
             of the 18 given TCGA cancer types ",
             row.names = TRUE)

dataFilt <- getDataGEO(GEOobject = "GSE20347",platform = "GPL571")

dataFilt <- getDataGEO(TCGAtumor = "ESCA")

Analysis: To analyze TCGA data

dataDEGs <- DPA(dataFilt = dataFilt,
                dataType = "Gene expression")

data(GEO_TCGAtab)
DataAnalysisGEO<- "../GEO_dataset/"
i<-5

cancer <- GEO_TCGAtab$Cancer[i]
cancerGEO <- GEO_TCGAtab$Dataset[i]
cancerPLT <-GEO_TCGAtab$Platform[i]
fileCancerGEO <- paste0(cancer,"_GEO_",cancerGEO,"_",cancerPLT, ".RData")

dataFilt <- getDataGEO(TCGAtumor = cancer)

GEOdegs <- DPA(dataConsortium = "GEO",
               gset = dataFilt ,
               colDescription = "title",
               samplesType  = c(GEO_TCGAtab$GEO_Normal[i],
                                GEO_TCGAtab$GEO_Tumor[i]),
               fdr.cut = 0.01,
               logFC.cut = 1,
               gsetFile = paste0(DataAnalysisGEO,fileCancerGEO))

library(TCGAbiolinks)
TCGAVisualize_volcano(DEGsmatrix$logFC, DEGsmatrix$FDR,
                      filename = "DEGs_volcano.png",
                      x.cut = 7,
                      y.cut = 10^-5,
                      names = rownames(DEGsmatrix),
                      color = c("black","red","dodgerblue3"),
                      names.size = 2,
                      xlab = " Gene expression fold change (Log2)",
                      legend = "State",
                      title = "Volcano plot (Normal NT vs Tumor TP)",
                      width = 10)

data(DEGsmatrix)

BPselected <- "apoptosis"

BPannotations <- DiseaseList[[match(BPselected,names(DiseaseList))]]$ID
dataLPA <- LPA(dataDEGs = DEGsmatrix[1:5,],
               BP =  BPselected,
               BPlist = BPannotations)
DiseaseListNew <- dataLPA
names(DiseaseListNew) <- BPselected

data(DEGsmatrix)
dataFEA <- FEA(DEGsmatrix = DEGsmatrix)

plotFEA(dataFEA = dataFEA, additionalFilename = "_exampleVignette", height = 20, width = 10)