SpidermiR: An R/Bioconductor Package for Integrative Analysis with miRNA Data

Abstract

Gene Regulatory Networks (GRNs) control many biological systems, but how such network coordination is shaped is still unknown. GRNs can be subdivided into basic connections that describe how the network members interact e.g., co-expression, physical interaction, co-localization, genetic influence, pathways, and shared protein domains. The important regulatory mechanisms of these networks involve miRNAs. We developed an R/Bioconductor package, namely SpidermiR, which offers an easy access to both GRNs and miRNAs to the end user, and integrates this information with differentially expressed genes obtained from The Cancer Genome Atlas. Specifically, SpidermiR allows the users to: (i) query and download GRNs and miRNAs from validated and predicted repositories; (ii) integrate miRNAs with GRNs in order to obtain miRNA–gene–gene and miRNA–protein–protein interactions, and to analyze miRNA GRNs in order to identify miRNA–gene communities; and (iii) graphically visualize the results of the analyses. These analyses can be performed through a single interface and without the need for any downloads. The full data sets are then rapidly integrated and processed locally.

1. Introduction

Gene regulatory networks (GRNs) play a crucial role in many key biological processes, such as cell differentiation, metabolism, cell cycle, and signal transduction. For instance, when a pathological process is ongoing, the dynamics of GRNs are most likely altered. Thus, the differences observed in GRNs in healthy and pathological conditions may reveal the mechanisms behind disease onset and progression.

A wide number of computational tools have been developed for analyzing and displaying GRN aided by the ever-increasing public availability of (big) data.

GeneMANIA [1] is one of the most popular web-based tools for studying validated biological networks. It enables different types of protein-protein and gene–gene interactions and networks to be predicted and visualized. Initially released in 2010, it currently indexes more than 2000 association networks, and maps more than 500,000,000 interactions from nine different organisms, extracted from a large number of publicly available sources.

MicroRNAs (miRNAs), which are small non-coding RNAs, are emerging as key drivers in the regulation of GRNs [2,3]. This is evidenced by the growing number of papers on this topic (from approximately 200 articles published in 2010 to 600 in 2015). MicroRNAs are particularly appealing as non-invasive tools for disease diagnosis, prognosis, and also therapy, since they have been found to be stably expressed in biofluids (such as serum and plasma) in several pathologies [4,5].

Some studies have provided insights into miRNA–gene regulation, showing that miRNAs tend to target highly connected genes in cellular networks, thus identifying miRNA–gene communities [6,7]. Since community structure conveys important information on specific functions [8], several community detection algorithms have been implemented to identify community structural properties, e.g., centrality measures, path length properties, and vertex and edge attributes.

Degree centrality (i.e., in the graphical representation of GRN, the number of neighboring vertices to which a node is directly connected) has been widely used for the analysis of biological networks [9]. It enables high-density nodes to be detected and ranked, e.g., proteins, genes, and miRNAs.

Several online resources provide collections of multiple databases that enable biologists to easily identify the biological functions and regulatory relationships between a group of known/putative miRNAs and validated/predicted target proteins, genes, and also associations with diseases or drugs (see, e.g., miRTAR [10], miRwalk [11], DIANA [12], Miranda [13], PicTar [14], TargetScan [15], miRandola [16], miR2disease [17], and Pharmaco-miR [18]). However, these resources need to be imported into other software tools for processing, tabulation, analysis, and graphing.

To overcome these limitations, more advanced tools have been recently made available to the scientific community such as collections of predicted and validated miRNA–target interactions (e.g., miRNATap [19], multiMiR [20], and MAGIA² [21]). Although each of these tools supports the user with particular features, only MAGIA² is able to automatically integrate the information on miRNA interaction networks with GRNs but with some limitations. For instance, there is no information on validated protein–protein or gene–gene interactions, or miRNA associations with diseases and drugs. MAGIA² only provides information on the interaction between transcription factors and the genes regulated by those transcription factors, omitting any other information regarding gene–gene interaction.

In this work, we describe a new software tool, namely “SpidermiR”, which automatically integrates with GRNs the information on miRNA interaction networks, and the association with drugs and diseases. It enables users to query, download, analyze, and visualize miRNA data with respect to specific GRNs, including different types of gene–gene (protein–protein) interactions (e.g., co-expression, genetic interactions, pathways, physical interactions, and shared protein domains).

SpidermiR has already been integrated with the most comprehensive repository of human cancer molecular and clinical data, The Cancer Genome Atlas (TCGA) [22,23,24], collecting 33 different tumor types by sampling across hundreds of cases per tumor type.

Most studies on TCGA data have focused on differentially expressed genes (DEGs) among different groups of tissue samples, and these analyses have not yet been able to fully reveal the relationships among genes or their dependency with miRNAs. Thus, the identification of miRNA-gene communities enriched by DEG, as is possible with SpidermiR, could simplify the problem of data interpretation and clarify the role of miRNAs in the onset and development of a specific cancer.

SpidermiR is available as open source and open development software in the bioconductor platform and thus is easily visible and accessible to an active user community, thereby promoting the reproducibility and the transparency of results.

2. Results and Discussion

In this section, we present the key features implemented in the SpidermiR package, and highlight the advantages of this software compared to other packages currently available to the scientific community with similar purposes.

We then present two case studies that will help in clarifying the utility of SpidermiR for users. The two case studies are shown in addition to the executable R code in the online documentation [25,26].

2.1. Key Features

Figure 1 highlights the key features of SpidermiR.

SpidermiR allows users to: (i) query and download GRNs from GeneMANIA and miRNAs from validated and predicted repositories, and harmonize annotations for miRNAs, genes, and proteins (query/download/annotation harmonization); (ii) integrate miRNA data with GRNs in order to obtain miRNA-regulated networks (miRNA–gene–gene and miRNA–protein–protein interactions) (enrichment), analyze miRNA GRNs to select specific interactions (interaction selection) and identify miRNA–gene communities (community detection); and (iii) graphically visualize and quantitatively summarize the miRNA-GRNs (graphics/metrics).

These analyses can be performed without having to navigate and access different web-based databases, without the need to download data, and by integrating and locally processing the full data sets in a short time.

SpidermiR is implemented as an R package, licensed under the General Public License (GPLv3). It is freely available, through the Bioconductor repository, at [25,26].

2.2. Benchmarking

SpidermiR offers easy access to both GRNs and miRNAs, and integrates this information with differentially expressed genes obtained from TCGA. It has several advantages with respect to other tools for integrative network analysis with miRNA data: the main difference is that in the competing tools, the role of miRNAs in the gene network has scarcely been considered and assessed.

As specified in the introduction, GeneMANIA is a web portal enabling users to build and visualize a composite gene–gene (both predicted and validated) or protein–protein (validated) functional interaction network. Currently, there are no functions that allow users to download these data in R or in Bioconductor in order to analyze them and integrate the networks with other information.

miRNAtap is a software package integrated in Bioconductor, solely focused on four predicted miRNA–gene databases (Miranda, PicTar, TargetScan, miRandola), and it integrates miRNA target prediction from different sources, aggregating them with various methods. Compared to miRNAtap, SpidermiR provides more features such as miRNA pharmaco/disease association, network data, analysis function, and miRNA–gene databases.

MultimiR is an R script that provides users with both predicted data from eight external databases (DIANA, EIMMo, MicroCosm, miRanda, miRDB, PicTar, PITA, and TargetScan), and validated miRNA–gene interactions from miRecords, miRTarBase, and TarBase. It also provides the users with both disease and drug associations with miRNAs from several databases (miR2Disease, Pharmaco-miR and PhenomiR). Although MultimiR has a greater number of predicted databases it does not allow the user to integrate the GRN and TCGA data. Consequently, although several tools have recently been developed to improve miRNA–gene identification almost none of them address the complex problem of reconstructing miRNA–gene target–genes in the different network types.

In fact, the association miRNA–gene is not sufficient to study miRNA activity. miRNA activity should be studied as part of a compound system where the single elements can interact among them on multiple levels. To understand such a complex system, it is crucial to study the multiple relationships between miRNA–gene target-genes and their roles in various diseases.

MAGIA² is an updated extension and evolution of MAGIA [27], a web tool for the integrated analysis of miRNAs and gene expression data aimed at the identification of GRNs. MAGIA² reports validated and predicted miRNA–gene interactions from several miRNA databases supporting multiple organisms (human, mouse, rat, and drosophila) and using a large list of prediction algorithms. Although MAGIA² is the most advanced tool currently available for integrative miRNA–gene analyses, it is not able to integrate the information on validated protein–protein or gene–gene interactions, such as those reported by GeneMANIA, or those related to diseases and drugs as those implemented in MultimiR. MAGIA² can estimate some of the relationships between transcription factor (TF)-gene interactions and miRNA–gene interactions, but the predicted regulatory circuits are limited to gene-TF interactions. Conversely, SpidermiR performs an integrative analysis on more comprehensive interaction types, and it provides all the features offered by the other competing tools, also including extracellular and circulating miRNAs and validated gene–gene interactions.

Compared to MAGIA², SpidermiR, has other advantages: (1) to identify regulatory circuits, users do not need to locally download the gene/miRNA expression data for each analysis; (2) since it is open software available in Bioconductor, users can integrate other analyses that are not already implemented in the tool; and (3) users can integrate miRNA–gene interaction into different network data, such as co-expression, physical interactions, genetic interactions, shared protein domains, co-localization and pathways, and are able to select miRNA–GRNs communities.

In addition, compared to online web portal, packages in the statistical software R provide more user-controlled features with no limitations on the number of miRNA/genes to be used.

The particular features (and sub features) of GeneMANIA, mirnaTap, multiMiR, MAGIA², and SpidermiR tools are summarized in

Table 1. Each column represents a software tool, and each row represents a feature. The checked cells indicate features that are implemented in the tool. Available platform abbreviations are defined as: R (only R script), B (integrated in Bioconductor package), and W (available as a web portal). * only gene-transcriptional factor, ** co-expression, genetic interaction, co-localization, and pathway, *** physical interactions, shared protein domain.

Features	Sub-Features	SpidermiR	GeneMANIA	miRNATap	multiMir	Magia2
Availability	Platform	B	W	B	R	W
	Functions for data query/download/annotation harmonization	●		●	●	●
	Functions for GRNs enrichment	●				●
	Functions for interaction selection	●				●
	Functions for community detection	●
	Functions for miRNA–GRNs graphics	●				●
	Functions for computation of miRNA–GRNs metrics	●				●
	Expression level	●				●
Interaction Type	Predicted miRNA–gene	●		●	●	●
	Validated miRNA–gene	●			●	●
	Disease–miRNA	●			●
	miRNA–gene–drug	●			●
	Extracellular/circulating miRNA	●
	Predicted gene–gene	●	●			● *
	Validated gene–gene	● **	● **
	Validated protein–protein	● ***	● ***
	Validated gene–gene–miRNA	●
	Validated protein–protein–miRNA	●
	Predicted gene–gene–miRNA	●				● *

.

2.3. Case Study No. 1. Prostate Cancer: Role of miRNAs in Shared Protein Domains

In this case study, we illustrate an application of SpidermiR aimed at identifying miRNAs that control shared protein domains in aggressive prostate cancer (PC). Combining gene expression data on aggressive PC samples from TCGA and shared protein domains from GeneMANIA, SpidermiR enabled us to detect a community of shared protein domains enriched by DEGs with direct gene-gene interactions, regulated by specific miRNAs. These miRNAs regulators could be interesting as therapeutic tools in aggressive PC.

Domains are generally considered as the conserved structural and functional units of proteins, and previous studies have shown that about 80% of protein pairs sharing the same domain also share the same function and have an a high degree of sequence similarity [28,29]. In the literature, there are several examples of proteins with similar domains, which are involved in cancer development: for instance, the protein family of Src, containing among other proteins Fyn and Yes, possesses all SH domains that are mutated in breast cancer, colon cancer, head and neck carcinoma, and non-small cell lung cancer [30].

We started with a network of shared protein domain interactions, consisting of 16,502 nodes (proteins) and 1,041,003 edges (interactions)—using SpidermiRquery, SpidermiRdownload, and SpidermiRprepare—and we enriched this network with miRNA-gene interactions, including only miRNAs already found to be deregulated in PC according to the literature, using SpidermiRanalyze_mirna_gene_complnet. This second network consisted of 13,338 nodes, (74 miRNAs and 13,264 genes) and 100,016 edges.

We identified 421 DEGs from TCGA data between aggressive PC (Gleason Sore ≥ 7; n = 452) versus normal samples (n = 52), and we integrated DEGs into the second network, thus obtaining a third network of only 80 nodes (13 miRNAs and 67 genes) and 86 edges, containing DEGs which interact directly, using SpidermiRanalyze_DEnetworkTCGA [25,26].

This network was processed to find communities based on density metrics and we then focused on the community with the highest number of elements, resulting in a subnetwork of 47 nodes (5 miRNAs and 42 genes) and 59 edges (fourth network), using SpidermiRanalyze_Community_detection and SpidermiRanalyze_Community_detection_net. Three nodes were miRNAs with direct interactions with three genes with a high centrality (hubs), which are directly connected to all the other genes in the community. The modulation of these three miRNAs could help regulate the community, thus pointing to these miRNAs as potential therapeutic tools in aggressive PC. A graphical visualization of these communities is presented in Figure 2.

miR-17-3p has already been found to play a role in PC growth and invasion [31], and has been proposed as a biomarker secreted in the urine of PC patients [5]. miR-29a and b have been proposed as tumor suppressors in PC [32,33,34].

The HOX gene family includes various transcription factors that help to determine cell and tissue identity in the early development of PC [35], and increase greatly during the development of PC [36].

FAT2 belongs to those Cadherin superfamily members that are homologous to Drosophila Fat, which functions as a regulator of Drosophila wing development [37]. FAT2 has been described as a regulator of invasiveness in breast cancer [38]. The literature has not yet clarified the role of FAT genes in PC, nor in other types of human cancers, where they could act as a controller of the invasion ability of the tumor cells by regulating the FAT Hippo signaling cascades [39].

Prostatic calcification is commonly observed and has been associated with PC [40]. Studies on tumor endothelial cells isolated from mice that spontaneously developed PC revealed that the calcification process in prostate vessels is sustained by the chondrogenic differentiation of endothelial cells. This population of cells, which is differentiated from chondrocytes, expresses increased levels of cartilage-specific markers, such as the type II collagen (COL2A1) gene, thus encoding the extracellular matrix protein type II collagen [41]. This gene was found in our analysis as one of the potential hubs regulated by miR-29a/b.

In addition several members of the HOX family (e.g., HOXB6, B7) that belong to this subnetwork, are involved in PC onset [35] or development [36].

2.4. Case Study No. 2. Breast Cancer: The Role of miRNAs in Regulating High Degree Centrality Proteins in Physical Interactions

In this case study, we illustrate an application of SpidermiR aimed at identifying a network of key miRNAs involved in breast cancer (BC), which directly interact with proteins in physical interaction (PI). PIs are referred to as protein–protein interactions, and are crucial for assessing the structural and functional architecture of the cell in terms of how direct associations between molecules constitute protein complexes, signal transduction pathways and other cellular machinery [42]. There is a positive correlation between protein connectivity, described by PI parameter, and regulatory miRNAs, as demonstrated by Liang et al. [43]. Moreover, to assess protein connectivity, degree centrality parameter is a key measure. Proteins with a high degree centrality are usually more essential for the survival and growth of the organism than proteins with low degree centrality values [9,43]. However, to our knowledge, there is no comprehensive study on the relationship between proteins with high degree centrality and their regulating miRNAs with high degree centrality in this case study, we assessed the hypothesis that miRNAs with a high degree centrality may preferably regulate proteins with high degree centrality.

We downloaded PI data in Homo sapiens from GeneMANIA, and we obtained a first network with 15,407 nodes (all proteins) and 189,030 edges, using SpidermiRquery, SpidermiRdownload, and SpidermiRprepare. We integrated miRNAs already found as deregulated in BC in the literature within this network, thus obtaining a second network of miRNA-protein PI including 830 nodes (62 miRNAs and 768 proteins) and 1001 edges, by using SpidermiRanalyze_mirna_network [25,26].

Next, the degree centrality of the proteins in the first network and of miRNAs in the second network were computed by using the function SpidermiRanalyze_degree_centrality. Proteins were clustered according to their degree centrality (50 proteins/cluster), identifying more than 300 clusters of proteins, and we then assessed the relationships between miRNAs with higher centrality and their target proteins with higher centrality. Statistical results showed that proteins with higher centrality are effectively targets of miRNAs with higher centrality (Pearson correlation coefficient = 0.7) (see Figure 3). We then focused on the cluster of proteins with the highest degree centrality (cluster No. 1): we found eight proteins were interacting with seven miRNAs with a high degree of centrality (see

Table 2. The 10 top miRNAs with the highest degree of centrality (d.c.) and their target proteins with the highest degree of centrality.

miRNA	Target Protein
miR-17-5p (d.c.296)	APP (d.c. 2008) HSP90AA1 (d.c. 773) MYC (d.c. 570)
miR-125b (d.c.55)	TP53 (d.c. 630)
miR-146a (d.c.38)	TRAF6 (d.c. 517) BRCA1 (d.c. 416)
miR-30a-3p (d.c.28)	EP300 YWHAE (d.c. 390)
let-7a (d.c. 26)	MYC (d.c. 570)
miR-429 (d.c.26)	MYC EP300 (d.c. 492)
miR-145 (d.c.23)	MYC

).

Figure 4 shows this network, describing the key miRNAs involved in BC and directly interacting with proteins in PI (only the direct interactions among BIs of the network are shown). For some of the identified miRNAs, the regulation of genes with a high degree of centrality has already been demonstrated. For instance, let-7a [44] and miR-145 [45] are known regulators of the high degree of centrality protein MYC, and miR-125b is the controller of the high degree of centrality protein p53 [46]. miR-17-5p, the miRNA with the highest degree centrality in the network, has V-Myc Avian Myelocytomatosis Viral Oncogene Homolog (MYC), Amyloid β Precursor Protein (APP), and HSP90AA1 as target mRNAs, which also present a high degree of centrality.

3. Materials and Methods

The SpidermiR package consists of functions that can be grouped into three main levels: (i) data; (ii) analysis; and (iii) visualization. In the following sections, we briefly describe the pipelines that can be used in case studies through these functions.

3.1. Data

The SpidermiR data function includes three main function categories: (i) query; (ii) download; and (iii) annotation harmonization.

(i)

“Query” enables users to query: (1) recent and archived data and to identify the elements to download; and (2) species and network type via GeneMANIA database.

(ii)

“Download” enables users to download: (1) gene–gene networks as previously queried; (2) miRNA validated data targets using miRTar and miRWalk databases [10,11]; (3) miRNA predicted data targets using DIANA, Miranda, PicTar, and TargetScan databases [12,13,14,15]; (4) extracellular/circulating miRNAs using the miRandola database [16]; (5) the associations among miRNAs, genes, and diseases, using the miR2Disease database [17], and among miRNAs, genes and drugs, using the Pharmaco-miR database [18].

(iii)

“Harmonization” enables users to process the data for downstream analyses and it prepares a matrix of gene networks by mapping Ensembl Gene ID to gene symbols. Gene symbols are needed to integrate miRNAdata.

3.2. Analyses

The analysis functions are designed to process network data through standard and novel computational methods. Once the network data have been prepared with Gene symbols ID, the downstream analyses can be divided into: (i) enrichment; (ii) interaction selection; and (iii) community detection.

(i)

“Enrichment” enables users to: (1) enrich the networks with some further biological information. For example, for each network users can integrate miRNA databases (validated or predicted) in order to find miRNA–gene target interactions in the downloaded gene network; (2) retrieve the information on miRNA–gene and gene-pharmaco from the Pharmaco-miR database; (3) retrieve the extracellular/circulating miRNA database in order to find miRNA–gene target interactions in the downloaded gene network; (4) enrich a chosen network with DEGs. Users can simply choose the type of tumor, platform, and the ID samples from the TCGA portal and then obtain the directed interactions of DEG among them [23,24].

In the enrichment step, SpidermiR combines interactions found in all validated databases, and it combines only interactions commonly found in at least two predicted databases.

(ii)

“Interaction Selection”. In this step, users can play with the obtained network. For example, user can: (1) find sub-networks including all direct interactions involving at least one of the biomarkers of interest (BIs)—this is carried out on the basis of a set of BIs, genes, miRNA, or both; (2) search for sub-networks including all direct interactions involving only BIs; (3) can search for sub-networks including all direct and indirect interactions involving at least one of the BIs; (4) find the number of direct neighbors of a BI and select those BIs with a number of direct neighbours higher than a given cut-off value.

(iii)

“Community detection”. In this step, users can analyze the network to detect communities using algorithms developed in the study by Csardi et al. [47], and characterize them in terms of the number of community elements (both genes and miRNAs). On the basis of a community to which some BIs belong, the community can be characterized as a network of elements (both genes and miRNAs), and users can find out whether or not a set of BIs is included within such a community.

3.3. Visualization

The visualization functions enable users to display the results of the analyses through a graphical representation of networks (i.e., with vertices, nodes and edges) and plots (for other results).

For example, users can see a 3D representation of the network [48] in different colors for miRNAs, genes, and drugs, and manage the network directly by moving/shifting the nodes and the edges, according to the analysis and visualization needs. Users can highlight specific BIs within a network with a different color.

Some metrics are computed and plotted for the network, such as the number of direct neighbors of BI (i.e., the degree centrality), the cumulative frequency distribution of degree centrality of communities, and the adjacency matrix of the community, representing the degree of connections among the nodes. A summary of the networks is also provided showing the number of edges, nodes, and miRNAs.

4. Conclusions

The huge amount of data on validated biological networks currently made available in public repositories offers an excellent opportunity to interpret molecular mechanisms and molecular dynamics, as well as to promote discoveries of novel diagnoses, prognoses, treatments, and monitoring protocols to improve clinical outcomes for patients.

However, the tools available to exploit these data are not exhaustive, and do not offer the user a workflow to download, integrate, and/or analyze the data with miRNAs in an environment that can also provide access to other statistical analysis methods, such as those provided by the increasingly popular Bioconductor repository.

In this paper we have presented SpidermiR, a software tool freely available on the Bioconductor platform, which provides functions to query, download, and process biological network data.

With the various functions, users can incorporate additional biological information retrieved from public databases in the selected network, such as miRNA data (validated, predicted, and extracellular circulating miRNA), disease annotations, drug associations, and TCGA data.

In addition, using interactive and community detection users can apply standard and novel methods for the analysis of miRNA-gene-gene interactions, and this is especially useful in the field of genomics and epigenomics research, also in identifying diagnostic and prognostic biomarkers for specific tumors and simulating candidate miRNAs as therapeutic agents. The visualization functions enable users to display the results generated by the analyses through a graphical representation of networks.

These functions provide an easy, time-saving computational instrument for complex investigations, without having to navigate through different web and data portals and integrating different computational platforms.

Our two case studies highlight the typical uses of our package, such as protein-protein interactions (shared_protein_domains and physical_interactions), gene–miRNA interactions, gene-gene–miRNA, miRNA disease associations, TCGA data, community detection, and degree centrality. Our studies also demonstrate how the package can generate candidate biomarkers, and network communities and their miRNA regulators, which could subsequently be tested experimentally in the laboratory.