isomiRTar: a comprehensive portal of pan-cancer 5′-isomiR targeting

TCGA data collection and processing

RNA-seq and miRNA-seq read count tables for primary solid tumor samples were downloaded from the GDC portal (https://portal.gdc.cancer.gov/). Thirty-one TCGA projects were selected for the downstream analysis (Table S1). First, library size normalization of RNA-seq data was independently conducted for each cancer type using edgeR TMM algorithm implementation (Robinson, McCarthy & Smyth, 2010), and the default low-expressed gene removal procedure was applied. After normalizing for transcript length, TMM-normalized fragments per kilobase of transcript per million mapped reads (TMM-FPKM) tables were generated. Finally, TMM-FPKM tables were ommllog2-transformed.

Processing of miRNA-seq data was started from the 5′-isomiR annotation from *.mirbase21.isoforms.quantification.txt files. For that, pri-miRNA genomic coordinates and local coordinates of canonical mature miRNA sequences were extracted from miRBase v21 (Kozomara, Birgaoanu & Griffiths-Jones, 2019). The standard 5′-isomiR nomenclature was utilized: a number after the “—” character denotes a shift from the canonical 5′-end in the 5′-3′ direction. For example, hsa-miR-192-5p—+1 differs from the canonical hsa-miR-192-5p miRNA by the absence of the first nucleotide on its 5′-end. As in the case of mRNA expression data, the edgeR package was used to filter out low expressed isomiRs and generate ommllog2-transformed TMM-normalized reads per million mapped reads (TMM-RPM) tables. It is well known that the overwhelming majority of miRNA sequencing reads correspond to the dozens of the most highly expressed isomiRs (Shao et al., 2010; Nersisyan et al., 2020a). Given that we marked a minimal set of 5′-isomiRs accounting for 95% of sequencing reads in a TCGA project as highly expressed (Table S1). A typical read count distribution is presented in Fig. S1 for the TCGA-COAD project.

Sequence- and correlation-based analysis of 5′-isomiR targeting

IsomiR target prediction was conducted with two commonly used tools: miRDB v6.0 (Chen & Wang, 2020) and TargetScan v7.2 (Agarwal et al., 2015).

The custom prediction mode of the miRDB web portal (http://mirdb.org) was programmatically called for each isomiR sequence using the Requests module of Python. The default 80 threshold value was set for miRDB target scores to select reliable interactions. Due to the technical limitations in miRDB, target prediction was performed only for isomiRs with a sequence length between 17 and 29 nucleotides.

Publicly available TargetScan Perl scripts and UTR data were obtained from the official web portal (https://www.targetscan.org/vert_72/). Dr. George Bell kindly provided Perl scripts for calculating cumulative weighted context++ scores (CWCS). TargetScan predictions almost always covered miRDB predictions (see Results). Because of that, and since there is no conventional CWCS cut-off value, we equalized the numbers of miRDB/TargetScan predicted targets for each 5′-isomiR. Specifically, for each 5′-isomiR we sorted the predictions according to the CWCS (lowest to highest) and then selected top n entries of the list, where n is the number of targets of this isomiR predicted by miRDB. Finally, filtered miRDB and TargetScan predictions were merged.

For a pair of 5′-isomiR and predicted target, we calculated Spearman’s correlation coefficient r between corresponding expression levels in each cancer type. The following thresholds were used to define the significance of interaction: r <  − 0.3, adjusted p-value <0.05. The same standard criteria for considering miRNA-target correlation significant was previously used in many works (Hu et al., 2016; Paquette et al., 2019; Nersisyan et al., 2020b). Multiple testing correction was performed differently depending on the analysis purpose:

1. Analysis of a cancer-specific isomiR-target regulatory network. In this case, the Benjamini–Hochberg procedure was applied to the list of p-values corresponding to all putative isomiR-target interactions (referred to as FDR_global in the isomiRTar interface).

2. Analysis of specific 5′-isomiR targets or specific target gene regulators in a fixed cancer type (out of other interactions context). In this case, the Benjamini–Hochberg procedure was applied only to the narrow list of corresponding p-values (referred to as FDR_local in the isomiRTar interface).

Functional annotation analysis

Functional annotation analysis of gene sets was performed with DAVID 2021 (Sherman et al., 2022). The default 0.05 cut-off value for FDR was used to select significantly enriched KEGG pathways.

Network and statistical analyses, web portal implementation

Networks were constructed and analyzed using the NetworkX v2.8 Python module (Hagberg, Schult & Swart, 2008). Maximal cliques of isomiR-isomiR co-expression networks were found with the find_cliques function.

SciPy v1.5 Python module (Virtanen et al., 2020) was used for the statistical analysis. Pandas v1.3 (McKinney, 2010), NumPy v1.19 (Harris et al., 2020), Seaborn v0.11 (Waskom, 2021), and Matplotlib v3.5 (Hunter, 2007) Python modules were extensively used for the miscellaneous computations and visualization.

Implementation of the isomiRTar website was based on Flask v2.0 Python framework and PostgreSQL database. All isomiRTar source codes have been made available on GitHub (https://github.com/s-a-nersisyan/isomiRTar).

Pan-cancer expression landscape of 5′-isomiRs

We compiled the comprehensive list of 1022 5′-isomiRs expressed in 31 cancer types using TCGA miRNA sequencing data. Out of them, 168 5′-isomiRs were marked as highly expressed in at least one cancer type, which included 31 non-canonical isomiRs. Approximately 38% of highly expressed isomiRs (63 entries) were expressed in 10 or more cancer types, including 11 non-canonical 5′-isomiRs: hsa-miR-10a-5p—+1, hsa-miR-10b-5p—+1, hsa-miR-22-3p—+1, hsa-miR-29a-3p—-1, hsa-miR-30a-3p—+1, hsa-miR-30e-5p—+1, hsa-miR-101-3p—-1, hsa-miR-142-3p—+1, hsa-miR-183-5p—+1, hsa-miR-192-5p—+1, hsa-miR-203a-3p—+1. Finally, 18 miRNAs (including almost all representatives of let-7 and miR-30 families) were expressed at high levels in all 31 TCGA projects (Fig. 1A). The rest of 5′-isomiRs (105 entries, 62%) were expressed in a cancer-specific manner (less than 10 cancers, Fig. 1B).

In most cases, the high expression of a non-canonical 5′-isomiR was tied with the high expression of its corresponding canonical form. There were only two exceptions: the canonical forms of hsa-miR-302a-5p—+3 and hsa-miR-519a-5p—-1 (both were highly expressed in TCGA-TGCT) were not detectable in the analyzed data, suggesting “incorrect” miRBase annotation or systematical sequencing biases. Interestingly, several miRNAs had three highly expressed 5′-isomiRs: hsa-miR-10a-5p (—-1, —0, —+1), hsa-miR-199a-3p (—-1, —0, —+1), hsa-miR-199b-3p (—-1, —0, —+1), hsa-miR-375 (—-1, —0, —+1), hsa-miR-514a-3p (—-3, —-2, —0) and hsa-miR-9-5p (—0, —+1, —+2).

Large-scale sequence-based prediction of 5′-isomiR target genes

We applied two state-of-the-art tools (miRDB and TargetScan) to predict target genes for all considered 5′-isomiRs. The tools are based on different ideas: the miRDB machine learning algorithm was trained on RNA-seq data of HeLa cells transfected with 25 miRNA mimics to predict miRNA targets based on binding site sequence and structure (Liu & Wang, 2019). TargetScan algorithm also scores miRNA bindings sites according to their cross-species conservation (Agarwal et al., 2015). TargetScan systematically yielded a significantly higher number of isomiR targets than miRDB, and the majority of miRDB targets were included in TargetScan predictions: median 95.5% with 95% confidence interval (CI) 85.7–100.0%. To harmonize these data, we equalized the numbers of predicted targets for each 5′-isomiR by selecting the respective number of the best-ranked targets according to TargetScan (see Materials and Methods for details). After such a harmonization, the predictions of miRDB and TargetScan had a median of 26.5% common targets, 95% CI [0.0–48.9]% (Fig. 2A). From here onwards, we consider the union of predictions by miRDB and TargetScan for each 5′-isomiR, unless indicated otherwise.

The number of sequence-based predicted targets significantly varied between 5′-isomiRs. The median number of targets was 209, and 95%-CI was wide: 6-1082 (Fig. 2B). The canonical form of hsa-miR-3613-3p had the highest number of putative target genes (2551), while five isomiRs had no targets which passed the defined predictions significance thresholds: hsa-miR-1307-5p—0, hsa-miR-1307-3p—+2, hsa-miR-1307-3p—0, hsa-miR-1307-3p—+1 and hsa-miR-371a-3p—+2.

We used the sequence-based predictions to assess the impact of one nucleotide difference in a miRNA 5′-end on its targetome. For each pair of 5′-isomiRs that differ in one nucleotide shift from the 5′-end (adjacent 5′-isomiRs), we calculated the Jaccard index (JI) for the sets of putative targets—the size of the intersection divided by the size of the union of the predicted target sets (Gösgens et al., 2021). JI = 1 indicates the complete equality of the target sets, JI = 0 indicates that the target sets are disjoint and intermediate values between 0 and 1 indicate the degree of the overlap. For example, the Jaccard index for two 5′-isomiRs with 20 common and 60 targets each equals 20/(60 + 60 − 20) = 0.2. Surprisingly, the median Jaccard index was equal to 0.06, 95% CI [0.00–0.31], indicating the average poor intersection between targetomes of adjacent 5′-isomiRs (Fig. 2C). Thus, different 5′-isomiRs originating from the same hairpin can generally be considered as independent functional units. Representative pairs of adjacent 5′-isomiRs are shown in Fig. 2: hsa-miR-10b-3p—0, hsa-miR-10b-3p—+1 (no common targets, Fig. 2D), hsa-miR-101-3p—0, hsa-miR-101-3p—+1 (median overlap with JI = 0.06, Fig. 2E) and hsa-miR-590-3p—-1, hsa-miR-590-3p—0 (maximum overlap with JI = 0.39, Fig. 2F). Next, we searched for nucleotide motifs associated with high values of JI. For that, 6-mer seed regions of the adjacent 5′-isomiRs were merged into 7-mer sequences (merged seed regions). Then, duplicate merged seed regions originating from closely related mature miRNAs (e.g., hsa-miR-19a-3p and hsa-miR-19b-3p) were removed. It turned out that 49 out of 50 merged seed regions corresponding to 50 adjacent 5′-isomiRs with the highest JI values started with uracil (χ² test p = 1.29 × 10⁻³⁰). Interestingly, we did not find any motifs for the seed sequences of adjacent 5′-isomiRs with non-overlapping targetomes.

Pan-cancer correlation analysis of 5′-isomiRs and their target genes

The conducted sequence-based target prediction allowed us to compare the putative non-cancer-specific targetomes of 5′-isomiRs. To determine which predicted targets could be downregulated by the corresponding 5′-isomiRs on transcriptome level, we selected isomiR-gene interactions with a significant negative correlation of corresponding expression levels in a set of TCGA samples (the analysis was separately performed for each cancer type). As a result, we could construct the 5′-isomiR-gene interaction networks for each cancer from TCGA. For a given 5′-isomiR and cancer type, we denote the number of significantly anti-correlated 5′-isomiR targets (i.e., the out-degree) as 5′-isomiR targeting activity (ITA).

ITA values substantially varied between cancer types, but generally, the median ITA did not exceed 10 (Fig. 3A), which was an order of magnitude lower than the median number of sequence-based predicted targets. Nevertheless, there were several isomiRs with broad functional activity (ITA >100) in almost all cancers. The number of sequence-based predicted targets was the primary determinant of ITA: the corresponding Spearman’s correlation coefficients were all positive and significantly different from zero for all cancers (see the vertical axis labels in Fig. 3A). Interestingly, these coefficients strongly varied between cancer types from 0.29 (TCGA-CHOL) to 0.79 (TCGA-TGCT) and were significantly correlated with the median ITA values (Spearman’s correlation = 0.8, p = 4.71 × 10⁻⁸).

Expression levels of 5′-isomiRs also positively and statistically significantly correlated with ITA for all cancer types (p < 10⁻³), though correlation coefficients were not high (between 0.2 and 0.3, Fig. 3B). We hypothesized that some low-expressed 5′-isomiRs had high ITA because of a co-expression with other active molecules. Given this, we first identified 5′-isomiR co-expression modules in each cancer type. Modules were defined as maximal cliques in corresponding co-expression graphs, where two isomiRs were connected with an undirected edge if absolute Spearman’s correlation between their expression levels exceeded 0.5. Then, a representative 5′-isomiR with the highest median expression level was selected from each module. Limiting the analysis to these representative isomiRs systematically increased the correlation values between ITA and average expression in all except two cancer types (Fig. 3B). Notably, there was a complete absence of correlation between ITA and expression in the set of highly expressed 5′-isomiRs in each cancer, indicating that high expression is a necessary but insufficient condition for the downregulation of multiple target genes. The same holds for 18 miRNAs with high expression in all 31 cancer types: ITA values varied from less than ten anti-correlated targets (hsa-miR-100-5p, hsa-miR-99b-5p, hsa-miR-10b-5p) to dozens and hundreds (Fig. S2). Indeed, low values of ITA do not imply the absence of isomiR activity, e.g., in the case of very target-specific regulation.

Overview of isomiRTar portal

To make the 5′-isomiR pan-cancer expression profiles and the results of targeting analysis easily accessible, we developed the isomiRTar web portal (https://isomirtar.hse.ru). The portal’s main page contains the forms for the five query modes: by cancer type, by 5′-isomiR, by gene, by 5′-isomiR, gene, and cancer triple.

The main feature of the “query by cancer” page is an interaction network between 5′-isomiRs and their target genes in a specified cancer type. As the network is highly interactive, a user can retrieve information about a 5′-isomiR, target gene, or their interaction by simply hovering over corresponding nodes or edges (Fig. 4A). Aside from the graphical representation, the page contains summary statistics for 5′-isomiRs in a selected cancer: median expression, number of predicted targets, number of predicted targets supported by the significant negative correlation. The second query type allows one to explore a pan-cancer expression profile for all 5′-isomiRs of a specified miRNA (Fig. 4B).

The analysis of targets for a given 5′-isomiR or 5′-isomiR regulators of a given mRNA can be performed in two modes: cancer-specific or pan-cancer. In the first case, the portal page contains a table with the results of correlation analysis between a selected 5′-isomiR and its predicted targets or predicted 5′-isomiR regulators and a selected target gene (Fig. 4C). By clicking to interaction, a user can see the scatterplot of corresponding expression values (Fig. 4D).

In the second case (pan-cancer mode), isomiR-target interactions are sorted according to the number of cancers, where we observed the significant anti-correlation between corresponding expression levels. Thus, “universal” (not cancer-specific) isomiR-gene interactions can be inferred from isomiRTar. To analyze the landscape of such “universal” interactions, we selected abundantly expressed 5′-isomiRs (highly expressed in 10 or more cancer types) and incident interactions, which were supported by significant negative correlations in at least half of cancers where the corresponding 5′-isomiR was highly expressed. The resulting network is illustrated in Fig. 4E. Aside from crucial and extensively studied interactions such as the downregulation of ZEB1 and ZEB2 by miR-200 (Hill, Browne & Tulchinsky, 2013; Nersisyan et al., 2021b), the network had several regulatory interactions mediated by non-canonical 5′-isomiRs: hsa-miR-29a-3p—-1, hsa-miR-101-3p—-1 and hsa-miR-192-5p—+1. The regulation of centrosomal protein 55 (CEP55) by hsa-miR-101-3p—-1 was the most pan-cancer conserved: a significant negative correlation was observed in 15 cancer types. Since CEP55 is a well-studied oncogene (Jeffery et al., 2016), this 5′-isomiR of hsa-miR-101-3p acted as a potential tumor suppressor, which agrees with the tumor suppressor role of the canonical hsa-miR-101-3p form (Wang et al., 2018).

The interaction network of 5′-isomiRs and their targets in colorectal cancer

We selected colorectal cancer (TCGA-COAD) to illustrate the detailed analysis of 5′-isomiR-gene interaction network (Fig. 5). Out of 55 highly expressed isomiRs 39 molecules (71%) had at least one anti-correlated target gene. While most nodes corresponded to canonical miRNAs, four non-canonical 5′-isomiRs had high out-degrees: hsa-miR-142-3p—+1 (ITA = 34), hsa-miR-22-3p—+1 (ITA = 31), hsa-miR-203a-3p—+1 (ITA = 17) and hsa-miR-192-5p—+1 (ITA = 12). Consistent with the sequence-based target prediction, targetomes of the adjacent 5 ′-isomiRs had a little or no overlap: seven common targets of hsa-miR-203a-3p—0 and hsa-miR-203a-3p—+1 (JI = 0.17), two common targets of hsa-miR-143-3p—-1 and hsa-miR-143-3p—0 (JI = 0.29), one common target for hsa-miR-192-5p—0 and hsa-miR-192-5p—+1 (JI = 0.04), while the rest of adjacent 5′-isomiR pairs did not share target genes.

We selected two pairs of adjacent 5′-isomiR to analyze the cancer-related targeting: hsa-miR-203a-3p (—0, —+1) and hsa-miR-192-5p (—0, —+1). Interestingly, all these isomiRs showed clear tumor-suppressive patterns. Both canonical and non-canonical isomiRs targeted well-known drivers of tumor progression, including master regulators of epithelial-mesenchymal transition (EMT), such as ZEB1 and ZEB2 (see ‘Discussion’ for details).

Next, we used sets of anti-correlated target genes as an input for the functional enrichment analysis. The significant findings were regulatory interactions mediated by hsa-miR-148a-3p—0: multiple pro-cancer pathways were enriched in the set of its target genes. This included focal adhesion (KEGG hsa04510), ECM-receptor interaction (hsa04512), TGF-β signaling pathway (hsa04350), and PI3K-Akt signaling pathway (hsa04151). Furthermore, the mentioned pathways shared many hsa-miR-148a-3p—0 targets: COL4A1, ITGA5, ITGA11 and LAMB2 were present in three pathways, FLT1 and ROCK1 were present in two pathways, while CSF1, FBN1, LTBP1, and TGFB2 were pathway-specific. Thus, we showed that the canonical miRNA hsa-miR-148a-3p potentially suppresses colorectal tumors.