cRegulome: an R package for accessing microRNA and transcription factor-gene expression correlations in cancer

Data sources

The database file was built based on two data sources: Cistrome Cancer and miRCancerdb databases (Mei et al., 2017; Ahmed et al., 2018b). Based on an integrative analysis of TCGA and public ChIP-seq data, Cistrome Cancer provided the expression correlations of transcription factors (n = 305) and their target genes in different cancer types (n = 29). In addition, Cistrome Cancer calculated the regulatory potential of transcription factors with target or non-target genes. Similarly, miRCancerdb used TCGA data and TargetScan annotations to determine the expression correlation between microRNAs (n = 684) and target or non-target genes in cancer types (n = 34).

Database file

The cRegulome can obtain and query a pre-built SQLite database file of the Cistrome Cancer and miRCancerdb databases. The details of this build were described in cRegulomedb (https://github.com/MahShaaban/cRegulomedb). In addition, this repository contains the scripts and documentation to pull, format and deposit the data at an on-line repository. Briefly, the SQLite database consists of four tables cor_mir and cor_tf for correlation values; and targets_mir and targets_tf for microRNA miRBase ID and transcription factor symbols to gene mapping. Two indices were created to facilitate the database search using the miRBase IDs and transcription factor symbols. The database file can be downloaded using the function get_db.

Database query

To query the database using cRegulome, we provided two main functions with get_mir and get_tf for microRNA and transcription factor correlations, respectively. Users need to provide the proper IDs for microRNA, transcription factor symbols and/or TCGA study identifiers. microRNAs are referred to by the official miRBase IDs, transcription factors by their corresponding official gene symbols and TCGA studies with their common identifiers. In either cases, the output of calling these functions is a tidy data.frame of four columns; mirna_base or tf, feature, cor and study. They correspond to each of the miRBase IDs or transcription factor coding gene symbols, target gene symbols, correlation values and the TCGA study identifiers.

Objects

Two S3 (simple storage service) objects are defined in cRegulome to store and dispatch methods on the correlation data: cmicroRNA and cTF for microRNA and transcription factors, respectively. The structure of these objects is very similar. Basically, as all S3 objects, it’s a list of four items; microRNA or TF for the regulome element, features for the gene hits, studies for the TCGA studies and finally corr is a named list of data.frames-one for each queried study. Each of these data.frames has the regulome element (microRNAs or transcription factors) in columns and features/genes in rows. To construct these objects, users need to call a construct function with the corresponding names on the data.frame output form get_*. The reverse is possible by calling the function cor_tidy on the object to get back the tidy data.frame.

Methods

cRegulome provides S3 methods to interact with and visualize the correlations data in the cmicroRNA and cTF objects. Table 1 provides an overview of these functions. These methods dispatch directly on the objects and could be customized and manipulated in the same way as their generics. Generally, the functions fall in one of two categories. One is for obtaining and querying the database file and another for managing and visualizing the output data.

Identifying common regulators of PEBP1 and autophagy/EMT genes

To illustrate the use of the cRegulome package in answering relevant biological questions, we studied candidate transcriptional regulators of RKIP, an anti-cancer gene product known also as phosphatidylethanolamine binding protein (PEBP1) (Ahmed et al., 2018a). In this study, we used the weighted-gene co-expression network analysis (WGCNA) to find the epithelial to mesenchymal transition (EMT) and autophagy gene products that potentially interact with PEBP1 during the development of prostate cancer. These interactions could be either a direct physical binding or an indirect relation. Gene products that share common upstream regulators are expected to show high co-expression at mRNA level. Using the cRegulome package, we identified two transcription factors and three microRNAs that commonly regulate PEBP1 and one or more of the eight gene products of interest.

Database query and target filtering

We first used the cRegulome to identify transcription factors and microRNAs that target PEBP1 in prostate cancer study (PRAD). Two transcription factors were also found to target one or more of the eight genes of interest and three microRNAs were found to target PEBP1 only (Table 2). The transcription factors were RCC Excision Repair 6, Chromatin Remodeling Factor (ERCC6) and Vascular Endothelial Zinc Finger 1 (VEZF1). The three microRNAs were hsa-mir-23c, hsa-mir-378c and hsa-mir-761. The full list of targets and their co-expression correlations were retained to investigate their global regulatory roles.

Visualization of the query output

cRegulome provide a set of functions to visualize the query outputs. Figure 2 is a dot plot of the filtered query. The expression correlations of each regulator are shown as a dot. The size of the dot is the magnitude of the correlation. Two points are immediately obvious in this figure. First, both negative and positive regulations are present for this set of target. Second, the target TBC1D5 has the strongest regulation with most regulators. The same regulatory connection between ERCC6 and TBC1D5 appear in multiple other cancer types (Table S1).

Another way to visualize the identified regulations is to use network graphs. In Fig. 3, we used genes, transcription factors and microRNAs as graph nodes. The edges between the nodes were mapped to the co-expression correlations. In more denser graphs, clustering and node influence can be identified. In this case, hsa-mir-378c and hsa-mir-23c, having the highest number of edges, are expected to exert a higher influence on the network. This functionality can serve as both visualization method and an interface to other network analysis workflows.

In addition to identifying regulatory relations of a set of genes of interest, cRegulome can be used to learn about the global role of certain regulators in different types of cancer. The two identified transcription factors, and to lesser extent the microRNAs, appeared to have both positive and negative regulatory roles in prostate cancer (Fig. 4). When this pattern was investigated in other types of cancer, both transcription factors tended towards the positive regulation (Fig. S1). In Fig. 5, we showed that while the two transcription factors share a significant number of target genes, the microRNAs do not. Although the shared number of targets between the two factors vary, it remains significant in several other cancer types.

Code for reproducing the examples in the article

In this section, we discussed the steps of a typical cRegulome database query, data filtering and result visualization. The full code for reproducing this analysis was attached to this article (Additional File 1) and is available online at https://github.com/BCMSLab/cregart.

First, we start by installing the package, loading the required libraries and downloading the database file in an R session.

 
 
# Load  required  libraries 
library(RSQLite) 
library(cRegulome) 
library(ggplot2) 
 
# Installing  cRegulome 
install.packages(’cRegulome’)                       # stable version 
devtools::install_github(’ropensci/cRegulome’)   # development  version 
    

 
 
# Downlaod  database  file 
if(!file.exists(’cRegulome.db’)) { 
   get_db(test = FALSE, 
           destfile = ’./cRegulome.db’) 
} 

 
 
# Alternatively: the  database  file  can be  downloaded  using wget 
wget  https://s3-eu-west-1.amazonaws.com/pfigshare-u-files/9537385/ 
     cRegulome.db.gz 
gunzip  cRegulome.db.gz    

Second, we queried the database for the expression correlation of two transcription factors (ERCC6 VEZF1) and three microRNAs (hsa-mir-23c, hsa-mir-378c and hsa-mir-761). and limited the output to the nine genes of interest.

 
 
# Identify  target  genes 
targets  <- c(’PEBP1’, ’PIK3C3’, ’PIK3CB’, ’TBC1D25’, ’TBC1D5 
    ’, ’TOLLIP’, ’WDR45’, ’WIPI1’, ’TGFBR1’) 
# Query  mir  with  targets = False 
conn  <- dbConnect(SQLite(), ’./cRegulome.db’) 
 
mir  <- get_mir(conn, 
                  mir = c(’hsa-mir-23c’, ’hsa-mir-378c’, ’hsa- 
                      mir-761’), 
                  study = ’PRAD’, 
                  targets = targets) 
 
tf <- get_tf(conn, 
                tf = c(’ERCC6’, ’VEZF1’), 
                study = ’PRAD’, 
                targets_only = TRUE, 
                targets = targets) 
 
dbDisconnect(conn)    

The code to generate Figure 2:

 
 
ob <- cTF(tf) 
 
p1 <- cor_plot(ob) + 
   theme(legend.position = ’top’, 
          legend.direction = ’vertical’) + 
   scale_size_continuous(breaks = c(.5, .65)) + 
   labs(x = ’’ ’’, y = ’’ ’’) 
 
ob <- cmicroRNA(mir) 
 
p2 <- cor_plot(ob) + 
   theme(legend.position = ’top’, 
          legend.direction = ’vertical’) + 
   scale_size_continuous(breaks = c(.1, .3)) + 
   labs(x = ’’ ’’, y = ’’ ’’)    

The code to generate Figure 3:

 
 
ob <- cTF(tf) 
 
p1 <- cor_igraph(ob, directed = TRUE) 
 
ob <- cmicroRNA(mir) 
 
p2 <- cor_igraph(ob, directed = TRUE)    

Next, we constructed a query to find all targets of the five regulators in prostate cancer.

 
 
# Query  the  database 
conn  <- dbConnect(SQLite(), ’./cRegulome.db’) 
 
all_tf <- get_tf(conn, 
                     tf = c(’ERCC6’, ’VEZF1’), 
                     study = ’PRAD’, 
                     targets_only = TRUE) 
 
all_mir  <- get_mir(conn, 
                       mir = c(’hsa-mir-23c’, ’hsa-mir-378c’, ’ 
                           hsa-mir-761’), 
                       study = ’PRAD’, 
                       targets_only = TRUE) 
 
dbDisconnect(conn) 
    

The code to generate Figure 4:

 
 
ob <- cTF(alltf) 
 
p1 <- cor_joy(ob) + 
   labs(x = ’Transcription  factor’, 
         y = ’Expression  Correlation’) 
 
ob <- cmicroRNA(allmir) 
 
p2 <- cor_joy(ob)+ 
   labs(x = ’microRNA’, 
         y = ’Expression  Correlation’)    

The code to generate Figure 5:

 
 
ob <- cTF(alltf) 
 
cor_upset(ob) 
 
ob <- cmicroRNA(allmir) 
 
cor_upset(ob)