Count-level TCGA isomiR expression data ( n = 11089 samples) were downloaded from the IsoMiRmap tool (Loher et al., 2021) official website (https://cm.jefferson.edu/isomirmap/). Reads exclusively mapped to the miRNA space were selected for further processing. Transcript-level mRNA-seq data for n = 10530 TCGA samples were downloaded from the UCSC Xena portal (Goldman et al., 2020) (https://xenabrowser.net/), TOIL RSEM expected_count dataset built on GENCODE 23 human reference genome was used.

Normalization of both miRNA-seq and mRNA-seq TCGA count-level data was performed using the median of ratios algorithm implemented as the “estimateSizeFactors” function in the DESeq2 R package (Love et al., 2014); “fpm” and “fpkm” functions were used to generate the final normalized expression tables for miRNA-seq and mRNA-seq data, respectively (i.e., mRNA-seq data was additionally normalized by the transcript lengths). Finally, we applied log 2 x + 1 transformation for both miRNA-seq and mRNA-seq data.

PAM50 molecular subtyping results of primary breast cancer samples and normal tissues annotation were extracted from the UCSC Xena portal. We set a minimum of 100 sample thresholds for all subtypes to eliminate possible sample size-related biases in the downstream analysis. As a result, four sample groups were considered: luminal A primary tumors ( n = 561 samples), luminal B primary tumors ( n = 210 ), basal-like primary tumors ( n = 183 ), and adjacent normal tissues ( n = 104 ). Sample identifiers and corresponding group labels are available in Supplementary Table S1.

Since 3′-end isomiR variations do not alter the isomiR targetome much (seed region not affected), we summed up the expression of isomiRs originating from the same miRNA, which have identical 5′-end sequences (5′-isomiRs). We used our previous notation for 5′-isomiRs (Nersisyan et al., 2022b): a number after the “|” symbol stands for 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.

For further processing, we selected 139 5′-isomiRs with the highest median expression—this was a minimum number of isomiRs that covered 99% of the total median isomiR expression. The used threshold corresponded to approximately 100 DESeq2-normalized FPM. A less strict thresholding was used for the mRNA-seq data: transcripts with zero reads in more than half of the analyzed samples were discarded. The final list of analyzed 5′-isomiRs and their median expression levels is available in Supplementary Table S2.

Two tools were used to predict 5′-isomiR targets: RNA22 (Miranda et al., 2006; Loher and Rigoutsos, 2012) and TargetScan 7.2 (Agarwal et al., 2015). Precomputed RNA22 predictions for the canonical miRNAs were downloaded from the official tool’s website (https://cm.jefferson.edu/rna22-full-sets-of-predictions/; ENSEMBL 96 and miRBase 22 version). For the remaining 30 non-canonical isomiRs, a program that allows one to submit target prediction batch requests was used (https://cm.jefferson.edu/rna22/Interactive/remoteRNA22v2.zip, ENSEMBL 96 transcript sequences). The default target predictions settings were used. Since human reference transcriptome versions differed for mRNA-seq reads mapping and isomiR target prediction, we selected transcripts with identical sequences in both versions for further analysis.

The set of predicted isomiR-mRNA interactions was classified into five types based on seed region binding motifs: 8mer, 7mer-m8, 7mer-A1, 6mer (classical seed binding types), and other (including unlimited G:U wobbles and up to one mismatch or unpaired nucleotide in a seed region). Custom TargetScan predictions were done as previously described (Nersisyan et al., 2022a): publicly available Perl scripts and 3′-UTR sequences were downloaded from the official website (https://www.targetscan.org/cgi-bin/targetscan/data_download.vert72.cgi). Assessment of binding with imperfect seed complementarity (3′-compensatory sites) was not possible with TargetScan custom prediction mode. The results of RNA22 and TargetScan predictions are available in figshare repository (https://doi.org/10.6084/m9.figshare.21579729.v1). MiRTarBase 9 database (Huang et al., 2022) was used to retrieve experimentally verified miRNA-mRNA interactions (luciferase reporter assays).

For a given 5′-isomiR, we composed a set of all mRNAs which do not contain the 6mer seed region binding site in mRNA sequences (including an unlimited number of G:U wobbles). The computational approach for solving this problem was based on the composition of the hash table, which mapped all possible 6mers to their positions in all mRNA sequences. The set of negative control interactions is available in Supplementary Data 1.

We calculated Spearman’s correlation coefficients and corresponding p -values between all 5′-isomiR and mRNA expression levels (separately for the four groups of BC samples). Benjamini–Hochberg procedure was used to adjust p -values.

For a given 5′-isomiR, we denote the number of predicted targets as n t a r and the number of predicted targets supported by a significant negative correlation of expression levels ( r < − 0.3 ,   F D R < 0.05 ) as m t a r ( m t a r ≤ n t a r ). The analogous terms were used for the isomiR non-targets: n n o n t a r (number of non-targets predicted) and m n o n t a r (number of anti-correlated non-targets: r < − 0.3 ,   F D R < 0.05 ). For the downstream text, we denote m t a r as isomiR targeting activity (ITA) and m t a r − n t a r m n o n t a r n n o n t a r as adjusted ITA. The statistical significance of an isomiR activity was assessed with the one-sided Fisher’s exact test applied to the following contingency table: m t a r n t a r − m t a r m n o n t a r n n o n t a r − m n o n t a r

Benjamini–Hochberg procedure was used to adjust p -values.

We used DESeq2 1.36 (Love et al., 2014) to identify differentially expressed isomiRs for each pair of the considered sample groups (luminal A, luminal B, and basal-like BCs, normal samples). An isomiR was considered as differentially expressed if the fold change was lower than 0.5 or higher than two and the adjusted p -value was lower than 0.05.

We used SciPy stats 1.9 (Virtanen et al., 2020) for conducting statistical analysis (Spearman’s correlation, Fisher’s exact test). Pandas 1.3 (McKinney, 2010) and NumPy 1.23 (Harris et al., 2020) were used for miscellaneous computations. Plots were constructed with Seaborn 0.11 (Waskom, 2021) and ggVennDiagram (Gao et al., 2021).

The study workflow is outlined in Figure 1. We used miRNA-seq and mRNA-seq data of TCGA-BRCA primary tumor samples corresponding to three molecular subtypes (luminal A, luminal B, basal-like) and adjacent normal tissue samples. We identified 139 highly expressed 5′-isomiRs and predicted their target mRNAs using two sequence-based target prediction tools: RNA22 and TargetScan. The tools are based on different ideas. Namely, TargetScan favors evolutionary conserved binding sites, relies on the classical seed region binding types (8mer, 7mer-m8, 7mer-A1, 6mer), and predicts targets only within 3′-UTR. In contrast, RNA22 is a pattern-based method, which does not consider target site conservation, predicts targets within 5′-UTR, CDS, and 3′-UTR, and is significantly more tolerant to mismatches in a seed region. Aside from target predictions, we also composed sets of non-targets (i.e., negative controls) for each 5′-isomiR. Non-targets were defined as mRNAs that do not contain regions complementary to an isomiR’s seed region (including possible G:U wobbles). Finally, we calculated Spearman’s correlation coefficient between expression levels of all 5′-isomiRs and their targets and non-targets.

Before analyzing the joint expression profiles of isomiRs and their putative targets, we compared the sets of predicted isomiR-target interactions between RNA22 and TargetScan tools. Only 9.9% of RNA22-predicted target sites corresponded to the classical seed region binding (4.0% for 6mer, 3.2% for 7mer-m8, 1.6% for 7mer-A1, and 1.2% for 8mer), while the vast majority of putative targets had non-classical seed binding (90.1%, Figure 2A). High numbers of target sites were found in 3′-UTR and CDS (57.3% and 34.3%, respectively), while 5′-UTR binding was reported less frequently (8.3%, Figure 2B). TargetScan predictions included only classical seed region binding types within 3′-UTRs, with a prevalence of 6mer sites (53%, Figure 2C). RNA22 predicted 1.2 times more isomiR-mRNA interactions than TargetScan, and the tools shared 26.7% of the putative interactions (Figure 2D). At the same time, RNA22 missed the majority of classical seed region binding sites in 3′-UTRs compared to TargetScan (Figures 2E–H). For the downstream analysis, we used the union of RNA22-and TargetScan-predicted isomiR-mRNA interactions. With such a target prediction strategy, we were able to cover 89.1% out of 1,608 experimentally validated interactions between canonical miRNAs and their target mRNAs. In particular, 49.6% of validated interactions were predicted by both RNA22 and TargetScan, 6.5%—only by RNA22, and 33.0%—only by TargetScan.

isomiRs tend to stronger anti-correlate with the putative targets rather than with the non-targets.

To systematically compare correlation profiles of 5′-isomiRs with their targets and non-targets, we introduced two measures of isomiR activity. The first one, denoted as isomiR targeting activity (ITA), was calculated as a number of predicted targets, which were additionally supported by the significant negative correlation. However, a negative correlation between expression levels of a 5′-isomiR and its target gene could be explained by non-random side effects, such as opposite-sided regulation of both molecules by a transcription factor. Therefore, to estimate the number of such false positive interactions (i.e., negatively correlated 5′-isomiRs and their target genes with no direct effect), we used the non-target mRNA sets. Specifically, background ITA was calculated as the fraction of non-targets that were significantly anti-correlated with the 5′-isomiR, multiplied by the number of predicted isomiR targets. Finally, adjusted ITA was calculated by subtracting background ITA from unadjusted ITA. In other words, adjusted ITA reflects the difference between the number of negatively correlated potential isomiR targets and the expected number of negatively correlated non-targets.

Note that near-zero or negative adjusted ITA does not imply the absence of targeting. For example, consider the synthetic case where the only target of a 5′-isomiR is a hub TF with multiple downstream target genes, which are not direct targets of the isomiR. In this case, unadjusted ITA equals one, while background activity is much greater since the isomiR is anti-correlated with TF targets by transitivity. Thus, the adjusted ITA of the isomiR is negative despite having one functional target.

A comparison of adjusted ITA values calculated using 5′-isomiR and mRNA expression levels for the union of RNA22 and TargetScan predictions (including all seed binding types and mRNA regions) is shown in Figure 3. As can be seen, the distributions were essentially skewed toward the positive direction for each seed region binding type and each sample group. Thus, for the majority of 5′-isomiRs, ITA was higher than the background correlations level. Notably, adjusted ITA values were dramatically higher in the normal tissues compared to the cancerous ones. Thus, we observed global downregulation of isomiR targeting in breast cancer. The tables with adjusted ITA scores (union of RNA22 and TargetScan) for three BC subtypes and normal samples are available in Supplementary Table S3-S6.

In consistence with the previous section, RNA22 and TargetScan had different contributions to the correlation analysis. Namely, in the case of RNA22 predictions, non-classical seed binding types (referred to as “Other” in Figure 3) had significantly higher adjusted ITA scores compared to the four classical seed types (Supplementary Figure S1). Moreover, in the case of luminal A and luminal B BC subtypes, the whole distribution of adjusted ITA values (except outliers) was concentrated near zero, which was not the case for the “Other” binding type. Near-equal adjusted ITA scores were found for the 3′-UTR and CDS binding, while adjusted ITA in 5′-UTRs was systematically lower. In contrast to the RNA22 case, an essential mass of the TargetScan-based adjusted ITA distribution was in the positive zone for each classical seed binding type (Supplementary Figure S2). Thus, combining two target prediction algorithms ultimately allowed us to make valid target predictions both for classical and non-classical seed binding sites. It is also worth noting that adjusted ITA scores calculated separately using RNA22 and TargetScan were highly correlated: Spearman’s r = 0.70 , p = 5.25 × 10 − 22 for Luminal A, r = 0.51 , p = 1.73 × 10 − 10 for Luminal B, r = 0.88 , p = 3.54 × 10 − 47 for Basal-like, and r = 0.85 , p = 1.31 × 10 − 39 for Normal samples.

We were not able to identify any further subgroups within the predicted targets, which would have significantly different adjusted ITA scores. The tested groupings included partitioning of RNA22-predicted interactions by ten quantiles of binding energy, separating RNA22-predicted non-classical target sites by a number of G:U wobbles and mismatches, and partitioning of TargetScan-predicted interactions by ten quantiles of weighted context++ score.

Comparison of the numbers of anti-correlated isomiR-mRNA pairs between the sets of predicted targets and non-targets with Fisher’s exact test allowed us to identify 5′-isomiRs with a statistically significant activity over the background. The mutual arrangement of the significant isomiR sets for four sample groups (luminal A, luminal B, basal-like, normal) is presented in Figure 4A. As can be seen, isomiR activity showed clear subtype-specific patterns. Consistently with the previous section, the maximum number of significantly active isomiRs was detected in normal breast tissue samples: 103 out of 139 considered isomiRs. Normal tissues were followed by the basal-like (56 isomiRs), luminal A (34 isomiRs), and luminal B (12 isomiRs) BCs. Only eight isomiRs were significantly active in all BC subtypes: hsa-miR-17-3p|0, hsa-miR-30b-5p|0, hsa-miR-93-5p|0, hsa-miR-101-3p|-1, hsa-miR-106b-5p|0, hsa-miR-182-5p|0, hsa-miR-200b-3p|0, and hsa-miR-210-3p|0.

We then wondered whether the results of differential isomiR activity are in line with expression levels of these isomiRs in four considered sample groups. First, we noticed that the unsupervised clustering of samples based on the isomiR expression profiles was similar enough to the four considered sample groups (importantly, normal samples clustered separately from the cancerous ones, Supplementary Figure S3). Then, we conducted differential isomiR expression analysis for each pair of sample groups (Figure 4B). In contrast to the higher number of active isomiRs in normal tissues, we saw a higher number of isomiR upregulation events in cancer than downregulation events (this held for each BC subtype). The other interesting observation could be made while comparing luminal A and luminal B subtypes: only four isomiRs passed the differential expression significance thresholds (fold change and adjusted p -value), while the Jaccard index for the sets of the active isomiRs (9 isomiRs in the intersection divided by 37 isomiRs in the union) was equal to just 0.24. The detailed analysis of intersections between the sets of active/not active and upregulated/downregulated isomiRs for each pair of sample groups also did not allow us to find any strong associations between an isomiR’s expression level and significance of adjusted ITA score (Supplementary Figure S4). Finally, it is worth noting that we did not find a correlation between adjusted ITA scores and isomiR median expression levels neither in cancerous nor in normal samples (Spearman’s r < 0.18 ). Thus, the analyzed data told us that the variation of isomiR expression levels and correlations between isomiRs and putative targets constitute two independent dimensions.