Identification of neoplasm-specific signatures of miRNA interactions by employing a systems biology approach

Acquisition of miRNA-neoplasm associations

The miRNA and neoplasm associations were retrieved from four comprehensive publicly available databases: (a) miR2Disease (Jiang et al., 2009), a manually curated repository of miRNAs deregulated in diverse human diseases, (b) miRCancer (Xie et al., 2013) contains miRNA–cancer association data extracted through literature mining, (c) dbDEMC (Yang et al., 2017), a database of differentially regulated miRNAs in cancers, and (d) HMDD v3.2 (Huang et al., 2019b), a repository for experimentally verified human miRNA-disease relationships. In-house R scripts were used to extract the miRNA-cancer pairwise associations from those databases; duplicate entries were removed.

Nomenclature

In order to maintain a consistent nomenclature for the neoplasm names from the four databases, the disease terms were annotated according to the International Classification of Diseases and Related Health Problems (ICD), revision 11 (ICD-11) (Pickett & Anderson, 2018). ICD is maintained by the World Health Organization (WHO) and represents the global, comprehensive reference for the collection, processing and hierarchical classification of diseases, which can be used accurately and effectively to investigate the relationships among different diseases. NCBI’s Medical Subject Headings (MeSH) (Lipscomb, 2000), a comprehensive, hierarchically-organized vocabulary for subject indexing, cataloguing and searching biomedical and life sciences-related information, was used for assigning the neoplastic diseases to broader groups.

Official gene symbols were assigned to gene names according to the HUGO Gene Nomenclature Committee (HGNC) (Tweedie et al., 2021); the official miRNA names in miRBase (Kozomara, Birgaoanu & Griffiths-Jones, 2019) were used in this study.

Prediction of miRNA target genes

The corresponding gene targets of the miRNAs under investigation were identified. To this end, both the experimentally supported and predicted mRNA targets of the neoplasm-associated miRNAs were retrieved with the usage of the following state-of-the-art software tools: (i) microT_CDS has implemented the target prediction algorithm DIANA-microT-CDS (Paraskevopoulou et al., 2013); (ii) TargetScan searches for the presence of conserved sites that match the seed region of each miRNA (Agarwal et al., 2015); (iii) miRDB predicts functional miRNA gene targets by machine learning modeling of target gene expression data (Chen & Wang, 2020; Liu & Wang, 2019); (iv) PicTar relies on an algorithm for the identification of common microRNA targets (Krek et al., 2005); (v) miRanda, which is available online as part of the miRanda-mirSVR tool, depends on evolutionarily conserved miRNA-targets (Enright et al., 2003). To enhance the robustness of the prediction accuracy, only those target genes predicted by more than three tools were included in the study.

Integration of miRNA-miRNA interactions

The miRNA-miRNA pairwise functional relationships were assembled and inferred from two diverse sources, subsequently followed by merging and removal of duplicate miRNA-miRNA interactions.

GOSemSim (Yu et al., 2010), utilizes information-theoretic measurement to find similarities between genes based on Gene Ontology (GO) (Ashburner et al., 2000) terms. Thus, GOSemSim assigns a function to genes according to shared GO annotations. A miRNA-miRNA matrix was created based on the degree of functional similarity of their target genes; the higher the degree of gene relatedness, the higher the association score of the respective miRNAs.

In addition, miRNA-miRNA pairwise relationships were obtained from MIRGOFs (Yang et al., 2018), which measures functional similarities for miRNAs based on the GO annotation of their corresponding target genes by adopting two novel features: (i) it utilizes a new GO semantic similarity metric which considers both common ancestors and descendants of GO terms; (ii) it computes similarity between GO sets in an asymmetric manner, and weights each GO term by its statistical significance.

Only those miRNA pairwise interactions, from GoSemSim and MIRGOFs, with weight score ≥ 0.6 (where “1” is the highest) were selected to ensure robustness.

Construction of miRNA-neoplasm networks

Firstly, miRNA-neoplasm networks were constructed for each group of neoplasms, where the links were assigned two weights, one from MIRGOFs (Yang et al., 2018) and the second one from the calculated miRNA matrix based on GOSemSim.

The topological features of the generated networks were investigated through Cytoscape (v.3.9.1) (Shannon et al., 2003), a software platform for network manipulation and visualization. Arena3Dweb, an interactive, web-based tool for the visualization of multilayer networks in tree-dimensional space was also used for network depiction (Karatzas et al., 2021).

Identification of mutually exclusive cancer-relevant miRNA modules

To achieve the maximum possible mutual exclusivity, first, unique miRNAs were detected in each neoplasm. Second, the igraph package in R was used to calculate the node degree distribution (Barabasi, Gulbahce & Loscalzo, 2011) and clustering coefficient (i.e., an estimate of the tendency of the nodes in a graph to cluster together) (Watts & Strogatz, 1998) in each miRNA-neoplasm network (a total of thirty graphs). Then, the degree distribution was applied and a simultaneous comparison was performed between the 30 networks in order to find modules of miRNAs that are included in one neoplasm and not in the others. An in-house R script was used in order to handle the complexity of the datasets. Every dataset was included in a hash and the mutually exclusive modules were detected.

Neoplasm classification

All neoplasms were classified into thirty major groups using first the ICD-11 standard hierarchy for medical terms and then the MeSH vocabulary for assigning those terms to broader classes (Table 1). For example, the ICD terms ‘adenocarcinoma of lung’, ‘follicular carcinoma of thyroid gland’, and ‘adenoid cystic carcinoma’ were classified under the broader MeSH term ‘adenocarcinoma’. The ambiguous terms were removed from the subsequent stages of the analysis.

Construction of miRNA-miRNA graphs per neoplasm

The miRNA-neoplasm bipartite graphs were constructed by relying entirely on experimentally known miRNA-cancer associations assembled and integrated from four diverse resources.

The miRNA-miRNA interactions in the corresponding thirty monopartite networks were inferred from two different sources, directly and indirectly. In the former case, the pairwise miRNA associations were obtained from miRGOFS which includes miRNA-miRNA similarities computed with respect to the functional relationships of their target genes according to shared GO terms. In the latter case, miRNA targets were first detected through the consensus output of more than three miRNA-gene interactions databases, which include both experimentally supported and predicted entries. Then, GOSemSim was deployed, where the similarities for each pair of miRNAs were computed based on the functional annotation of their target genes according to GO; the more semantically related two GO terms, the higher the degree of association of the corresponding genes. The miRNA-miRNA pairs from both sources were combined and the duplicates were removed; a total of 2604 unique miRNAs were detected across the thirty neoplasms (Table S1).

The miRNAs hsa-miR-101-3p, hsa-miR-145-5p, hsa-miR-223-3p were found in all 30 types of cancers (Table S1). Furthermore, the well-studied oncomir miR-21, as well as the suppressomir miR-34a, were detected in 29 neoplasms, except testicular cancers. Fifteen unique miRNAs were found in ten neoplasms, namely, brain (hsa-miR-291-5p), colorectal (hsa-miR-103a-1-5p), head and neck (hsa-miR-6836-3p), liver (hsa-miR-199a-1-5p, hsa-miR-199a-2-5p), lung (hsa-miR-886-5p, hsa-miR-886-3p, hsa-mir-320a-3p), ovarian (hsa-miR-526a-5p), skin (hsa-miR-3689c), stomach (hsa-miR-7156-3p), thyroid (hsa-miR-5686, hsa-miR-199a-3p, hsa-miR-199b-3p), urinary bladder (hsa-miR-498-3p).

All thirty networks are very dense with an average estimated clustering coefficient ∼0.89 (where ‘1’ is the highest value). This indicates that the nodes in the given graphs exhibit a high tendency to cluster together. Despite the stringent cutoff values, the number of significant links is quite large, and as such, no visually discernible modules could be detected in any of the networks. Based on topological analyses, all graphs have a random degree distribution, where the node degrees (i.e., the number of connections of a given node to neighboring nodes) follow a Poisson distribution, unlike the biological networks which are scale-free and are dominated by few ‘hubs’ (i.e., highly connected nodes) (Barabasi, 2009; Barabasi, Gulbahce & Loscalzo, 2011; Goh et al., 2007; Kontou et al., 2016). This is probably attributed to the “promiscuous” nature of miRNAs to potentially target numerous genes, resulting in this way to a high number of corresponding miRNA connections.

In order to infer less dense and more robust networks, only those links with weights above the 90th percentile per network were selected. In this way, the highest scoring miRNA-miRNA interactions were used for network construction and subsequent module detection. After pruning, there were 1,893 unique miRNAs across neoplasms, with the number of miRNAs per neoplasm ranging from 39 (endocrine gland neoplasm) to 1,821 (colorectal neoplasms). Consequently, the number of the corresponding pairwise miRNA-miRNA associations was reduced significantly (Table 1 and Table S2). Noteworthy, this is observed in a more pronounced way in the more populated graphs, since the number of potential miRNA connections increases exponentially with an increase in the number of member miRNAs. As it is shown in Fig. 2, the endocrine gland neoplasm miRNA network becomes less dense with a reduced number of nodes (miRNAs) and connections (miRNA relations) when the upper 10th percentile of miRNA-miRNA links was used for graph construction.

Detection of neoplasm-specific miRNA modules

The goal of this study was to detect distinct subgraphs within the neoplasm graphs which cannot be contained by another subgraph. These subgraphs, which represent local clusters or modules, are considered semiautonomous and therefore are of high biological importance (Pelaez & Carthew, 2012). The top scoring mutually exclusive modules per neoplastic disease were selected based on the criteria that they exhibited significant mutual exclusivity and their constituent miRNAs were more highly connected to each other than to other miRNAs in the entire graph.

Modules of miRNA interactions that displayed a statistically significant level of mutual exclusivity were detected in seventeen neoplasms (Table 2). A single miRNA-miRNA pair (i.e., hsa-miR-1296-3p and hsa-miR-6748-5p) was detected in bladder cancer, which was excluded from the analysis. In the other twelve neoplasms under study, no distinguishable modules could be detected. The seventeen modules were comprised of three to twenty-three miRNAs (Table 2). Of note, only a few miRNAs were found in up to five cancer modules (hsa-miR-4273, hsa-miR-4422, hsa-miR-4720-3p and hsa-miR-5002-3p); this is because mutual exclusivity was a prerequisite. Among them were modules consisting of a relatively large number of miRNAs, such as the esophageal (23), colorectal neoplasms (20), and breast neoplasms (17), and also modules with a moderate number of component miRNAs including the liver (8) and lung neoplasms (7). Cliques of only three miRNAs were detected in the sarcoma and squamous cell carcinomas (Fig. 3). In addition, PubMed (https://pubmed.ncbi.nlm.nih.gov/) was searched with each of the neoplasms and the core miRNA components of the corresponding modules (Table 2). Nonetheless, in spite of extensive searches, any reported combinations of any two miRNAs that compose the neoplasm-specific modules could not be detected (Table 2).