miRDRN—miRNA disease regulatory network: a tool for exploring disease and tissue-specific microRNA regulatory networks

Data integration

miRDRM integrated data from several existing database on disease-miRNA association, miRNA-target gene association, gene ontology, biological pathway, and PPI. For disease-specific cases disease-associated miRNAs and targets were obtained from human microRNA and disease associations database (HMDD) (Li et al., 2014) (v2.0, http://www.cuilab.cn/hmdd). For non-disease specific cases, miRNA/siRNA and targets were obtained from TarBase (Vlachos et al., 2015) (v7.0, http://carolina.imis.athena-innovation.gr/diana_tools/web/index.php?r=tarbasev8%2Findex/). In disease-specific cases, the optional filter requiring miRNA-target association be assay validated used TarBase data; the filter excludes miRNA-target pairs appearing in HMDD but not in TarBase (if and when this happens). Gene-tissue associations were taken from NCBI-Entrez (NCBI Resource Coordinators et al., 2018); RSP associations with known pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2017) (http://www.genome.jp/kegg); gene-tumor associations from TAG (Chen et al., 2013) (http://www.binfo.ncku.edu.tw/TAG); OMIM IDs of target genes and genes on regulatory networks from NCBI-OMIM (NCBI Resource Coordinators et al., 2018) (https://www.ncbi.nlm.nih.gov/omim); gene IDs and transcription factor and/or receptor from NCBI-GeneBank (NCBI Resource Coordinators et al., 2018) (https://www.ncbi.nlm.nih.gov/genbank); data on human PPI from biological general repository for interaction datasets (BioGRID) (Chatr-aryamontri et al., 2013) (https://thebiogrid.org/); and gene ontology, gene annotation and gene product information from GO (Ashburner et al., 2000) (http://geneontology.org/). (Table 1).

Construction of miRNA-associated target-specific regulatory sub-pathways

Consider a linked sequence (M, T, G₁, G₂) (Fig. 1), where M is a miRNA, T is its regulatory target gene, G₁ is a gene whose encoded protein (p₁) interacts (according to PPI data) with the protein (p_T) encoded by T, and G₂ is a gene whose encoded protein (p₂) interacts with p₁. In what follows, when there is little risk of misunderstanding, the same symbol will be used to represent a gene or the protein it encodes. We call the sequence (T, G₁, G₂) a target-specific RSP, or simply a RSP, and (M, T, G₁, G₂) a miRNA-specific RSP (MRSP). Given a target gene T, we use PPI data from BioGRID to collect all RSPs by extending from T two levels of interaction.

Jaccard score of a regulatory sub-pathway

Jaccard similarity coefficients (Ng, Liu & Lee, 2009) were used to score the RSPs, based on the assumption that there is a tendency for two directly interacting proteins to participate in the same set of BP or share the same set of molecular functions (MF). Given two sets S1 and S2 (in the current application, a set will be either a list of BP or a list of MF, both according to GO), the Jaccard (similarity) coefficient (JC) of S1 and S2 is defined as, $$\mathrm{J}\mathrm{C}(S1,S2)=\frac{\left|S1\cap S2\right|}{\left|S1\cup S2\right|}$$

Where ∪ is the union (of two sets), ∩ is the intersection, and |Z| is the cardinality of Z. JC, which ranges from zero to one, is a quantitative measure of the similarity between two sets. For example, when S1 = {a, b, c} and S2 = {b, c, d}, JC (S1, S2) = 2/4 = 0.5.

Let (T, G₁, G₂) be an RSP as defined in the previous section and denote by [G] the set of BP (or pathways) (Kanehisa et al., 2017; Ashburner et al., 2000) that involve the gene G. We define the Jaccard score, or JS, of RSP as, $$\mathrm{J}{\mathrm{S}}_{\mathrm{X}}(T,G_1,G_2)=\frac{1}{2}\left(\mathrm{J}\mathrm{C}\left({\left[T\right]}_{\mathrm{X}},{\left[G_1\right]}_{\mathrm{X}}\right)+\mathrm{J}\mathrm{C}\left({\left[G_1\right]}_{\mathrm{X}},{\left[G_2\right]}_{\mathrm{X}}\right)\right)$$

Where X may be BP or MF. If the pair [T] and [G₁] do not share a common term, then the corresponding JC has a zero value; similarly for the pair [G₁] and [G₂]. In either case the RSP is considered to be not viable and discarded. In other words, miRDRN excludes any RSP with zero JC score.

p-Value of a regulatory sub-pathway

A p-value for an RSP (T, G₁, G₂) was assigned as follows. Let the total number of BP (or MF, as the case may be) terms be N, and the number of terms in [T], [G₁], [G₂], [T] ∉ [G₁], [G₁] ∉ [G₂] be x, y, z, n₁, and n₂, respectively, then the p-values, P₁ and P₂, for (T, G₁) and (G₁, G₂) are, respectively $$P_1={\displaystyle \frac{C_{n_1}^N{C}_{x-n_1}^{N-n_1}{C}_{y-n_1}^{N-x}}{C_x^N{C}_y^N}}$$and $$P_2={\displaystyle \frac{C_{n_2}^N{C}_{y-n_2}^{N-n_2}{C}_{z-n_2}^{N-y}}{C_y^N{C}_z^N}}$$

The p-value for the RSP was set to be the greater of P₁ and P₂.

Assembly and storage of target-specific regulatory sub-pathways

A union set of miRNA-associated target genes were collected from HMDD and TarBase and for every target a complete set of RSPs, with BP- and MF-type JC scores and p-values assigned, was assembled. The entire set of RSPs for all targets was stored in miRDRN.

Construction of disease-specific miRNA regulatory network

A user-initiated construction of a disease-specific miRNA regulatory network (RRN) proceeds as follows. Step 1. Select a disease. Step 2. Collect from HMDD all miRNAs (M’s) and target genes associated with the disease. Step 3. For each M and each of its targets retrieve from miRDRN storage all target-specific RSPs, thus forming a set of MRSPs. The union of the sets of MRSPs over all M’s is the set of disease-specific MRSPs for the selected disease. Step 4. Construct the disease-specific RRN from the set of disease-specific MRSPs by linking all unlinked pairs of genes/proteins if they have interaction according to BioGRID (Fig. 2).

miRNA disease regulatory network (miRDRN)—A database and web service platform

We built miRDRN (http://mirdrn.ncu.edu.tw/mirdrn/), a web-based service that allows the user to construct a disease and tissue-specific, p-valued, miRNA-protein regulatory network, or RRN. The current version of miRDRN contains 6,973,875 p-valued target-specific RSPs constructed through 389 miRNA-regulated genes from 207 diseases-associated miRNAs associated with 116 diseases (Table 2).

Comparison of miRDRN with other miRNA-related databases

A number of databases and/or web service platforms on miRNA-related topics are publicly available (Table 3). Aside from HMDD and TarBase on which miRDRN was built (Table 1), PhenomiR (Ruepp, Kowarsch & Theis, 2012) is a database on disease-miRNA association, miRwayDB (Das, Saha & Chakravorty, 2018) is a database on disease-miRNA-target and target-KEGG term association, and miRPathDB (Backes et al., 2017) is a database on miRNA-pathway association. New and unique as a database, miRDRN stores the 6,973,875 p-valued target-specific RSPs it has assembled (Table 2). As a web service platform miRDRN is a tool that facilitates the construction and visualization of disease-specific RRNs using these RSPs in combination with resources from HMDD, TarBase, and several other databases (Table 1).

Brief description of usage of miRDRN

miRNA disease regulatory network is reasonably user friendly; its many features are easily discovered by user exploration. Here, we give a brief description of its main features.

User may use miRDRN to explore a single disease, or the comorbidity of a disease-pair. In the course of either type of study, all relevant miRNAs, genes, and RSPs are made accessible to the user in tabulated form, and RRNs in the form of interactive maps, both of which may be downloaded by the user. Often a map is too large for practical visualization, and in such a case the user may use options such as setting a p-value cut-off, or requiring a specific gene to be present in the map, or both, to obtain a partial RRN.

The entrance interface of miRDRN (http://mirdrn.ncu.edu.tw/mirdrn/) asks the user to select “Single Search” to explore a single disease (or miRNA/siRNA) or “Comorbidity Search” to explore the comorbidity of a disease-pair (Fig. 3). The user is then asked to specify the disease or disease-pair to be explored and tissue/tumor types, and p-value threshold for RSP evaluation, and to click on (or not) several optional filters, respectively, on targets and on RSPs. The filter on miRNA targets allows the user to admit only targets positively validated by the seven direct experimental methods: HITS-CLIP, PAR-CLIP, IMPACT-Seq, CLASH, Luciferase Reporter Assay, 3LIFE, and Genetic Testing (Vlachos et al., 2015); filters on RSP allow the user to select only those RSPs with some or all of the proteins to be cancer related (Fig. 4). The user may then click on “Query” to start the computation. Tabulated results of disease-associated miRNAs and their target genes (Fig. 5), a multi-page list of all RSPs (Fig. 6) and, in the case of Comorbidity Search, a list of all comorbid genes (Fig. 7) will then automatically appear. After the first, automatic iteration, the user may reduce the size of the RSP-list by using the “Gene filter” and “Show top … sub-pathways” options (Fig. 6). The next interface (Fig. 8), in ready mode on first appearance, waits for the user to select one of three network layouts: “Tree,” “Circle,” or “Radial.” After “Go” is clicked on, the platform displays an interactive map showing the RRN built from RSPs selected by user-specified options (Fig. 8). When the mouse is placed on a node (a miRNA or a gene) on the map a small pop-up window opens to show the name of the node/gene and the number of other nodes it is linked to, and annotation on the node from GO, OMIM, KEGG, and GeneBank databases.

Case 1. A single disease study of colorectal neoplasm

Here, we demonstrate a single disease application of miRDRN. On the query interface select “colorectal neoplasms” (or colorectal cancer (CRC)), tissue type “colorectal tumor,” pathway ranking by “Jaccard index (MF),” and p-value < “0.001,” and no option on target genes or RSP. The result yielded 33 miRNAs and 37 target genes (if the option “miRNA target experimental validation is positive and direct" on the query interface was selected then there would be 20 miRNAs and 23 target genes (Fig. 5)), and 45,565 RSPs involving 3,079 genes (reduced to 2,111 RSPs and 1,650 genes when target is restricted to being “positive and direct” and cancer related (Fig. 6)) (Table 4).

By default, the interface “target-specific RSPs” (Fig. 6) lists all the constructed RSPs, namely all 45,565 of them in the present case and, if requested, would present a map including all the RSPs which, however, would be impractical to visualize, not to say interact with. On the same interface are two options for displaying/using a smaller RSP set: “Gene filter,” where the user can restrict the set to only those RSPs containing a specified gene; and “Show top … sub-pathways,” where the user can ask for only the top-N RSPs having the smallest p-values be listed and used for network construction. The interface “Disease specific RRN” then allows the user to choose one among the layouts “Tree,” “Circle,” and “Radial.” Here, a tree-map, with several disconnected parts, built from the top-70 RSPs is shown (Fig. 9).

The largest connected sub-RRN, or “Network-1” (Fig. 10), is composed of six miRNAs targeting four genes connected to 52 other genes (Table 5). Of the 56 genes in Network-1, 22 have known CRC connections (CORECG database, http://lms.snu.edu.in/corecg) (Agarwal et al., 2016), and 26 others have references linking them either directly or indirectly to CRC (Kang et al., 2016; Sirvent et al., 2010; Jeong et al., 2018; Xie et al., 2012; Wu, Wu & Jiang, 2017; Xiang, Wang & Xiang, 2014; Masuda & Yamada, 2017; Ali et al., 2014; Yun et al., 2018; Vázquez-Cedeira & Lazo, 2012; Zhang et al., 2009; Csukasi et al., 2018; Rey et al., 2016; Goyal et al., 2009; Sabir et al., 2015; Bjerrum et al., 2014; Hanna et al., 2018; Record et al., 2010; Zhou et al., 2017; Gong et al., 2018; Yuan et al., 2018; Kim et al., 2015; Li et al., 2013; Guo et al., 2012; Peng et al., 2014; Alonso et al., 2017; Qi & Ding, 2018; Jin et al., 2015) (Table 6). Among these, TNIK (Masuda & Yamada, 2017) and TNK2 (Qi & Ding, 2018) have been used as drug targets for CRC treatment. We consider the remaining eight genes—PRKACA, MAP3K12, LRRK1, RIOK2, OXSR1, CDK17, EIF2AK1, TSSK4—to be potential novel CRC-related genes. Noticeably, Network-1 has two parts, one 28 nodes (five miRNAs targeting three genes) and the other 34 nodes (one miRNA targeting one gene), connected by a single link, or PPI. The three types of genes, known CRC-related, reference-supported, and potential CRC-related, are more or less proportionately distributed in these two parts.

Table 5:

Statistics and gene information in the network-1, the largest connected sub-network of the CRC-specific miRNA regulatory network.

		Number	Item set
Network-1	miRNAs	6	hsa-mir-199a, hsa-mir-34a, hsa-mir-199b, hsa-mir-139, hsa-mir-497, hsa-mir-106a
	Target genes	4	AXL, IGF1R, RAP1B, TGFBR2
	Gene set (including target genes)	56	AXL, CSK, TNK2, LCK, PRKACA, FGR, MAPK15, IGF1R, MERTK, ERBB2, PTK2, EGFR, JAK2, JAK1, PRKCD, TEC, EPHB2, PHKG2, ROR1, FES, MAP3K12, RAP1B, MST4, PAK1, LRRK1, MAP2K3, CDK11B, ACVR1, TGFBR2, RIOK2, TGFBR1, MAP3K7, NEK8, NUAK2, OXSR1, CDK1, ACVRL1, MKNK2, STK35, CDK17, EIF2AK4, DAPK2, EIF2AK1, TSSK4, ZAK, MAP2K6, SIK3, VRK2, PINK1, TAOK2, TNIK, MAPK6, PRKACB, WNK1, PAK6, PKMYT1

DOI: 10.7717/peerj.7309/table-5

Table 6:

Known, literature supported, and potential novel CRC-related genes.

		Number	Item set
Gene set (Network-1)	Known CRC genes	22	AXL, LCK, FGR, IGF1R, MERTK, ERBB2, PTK2, EGFR, JAK2, JAK1, EPHB2, FES, PAK1, MAP2K3, ACVR1, TGFBR2, TGFBR1, CDK1, EIF2AK4, DAPK2, MAP2K6, PAK6
	Reference supported (Kang et al., 2016; Sirvent et al., 2010; Jeong et al., 2018; Xie et al., 2012; Wu, Wu & Jiang, 2017; Xiang, Wang & Xiang, 2014; Masuda & Yamada, 2017; Ali et al., 2014; Yun et al., 2018; Vázquez-Cedeira & Lazo, 2012; Zhang et al., 2009; Csukasi et al., 2018; Rey et al., 2016; Goyal et al., 2009; Sabir et al., 2015; Bjerrum et al., 2014; Hanna et al., 2018; Record et al., 2010; Zhou et al., 2017; Gong et al., 2018; Yuan et al., 2018; Kim et al., 2015; Li et al., 2013; Guo et al., 2012; Peng et al., 2014; Alonso et al., 2017; Qi & Ding, 2018; Jin et al., 2015)	26	CSK, TNK2*, MAPK15, PRKCD, TEC, PHKG2, ROR1, RAP1B, MST4, CDK11B, MAP3K7, NEK8, NUAK2, ACVRL1, MKNK2, STK35, ZAK, SIK3, VRK2, PINK1, TAOK2, TNIK*, MAPK6, PRKACB, PKMYT1, WNK1
	Potential novel cancer-related gene	8	PRKACA, MAP3K12, LRRK1, RIOK2, OXSR1, CDK17, EIF2AK1, TSSK4

DOI: 10.7717/peerj.7309/table-6

Note:

* Known target genes used in treatment of CRC.

The “Gene filter” option (Fig. 6) allows the user to focus on a specific gene in RRN construction. As example, TNK2, a key drug target for the treatment of metastatic CRC (Qi & Ding, 2018), was selected as the filter, together with the “Show top 70 RSPs” option. The result was a nine-node sub-RRN: the target gene AXL regulated by three miRNAs—hsa-mir-199b, hsa-mir-34a, hsa-mir-199a—and linked (by PPI) to TNK2, itself linked to four other genes AXL(OCG), MAGI3, HSP90AB2P, MERTK(OCG), KAT8 (Fig. 11).