Differential expression reveals novel breast-cancer associated miRNA

We performed a differential expression analysis comparing clinically relevant subgroups of breast cancer. We measured differential expression of a specific miRNA on a pair of sample subpopulations (e.g. ER positive vs ER negative). Fold-change was defined as the ratio (log₂) between median expression of both sets. We applied Wilcoxon Rank-Sum (WRS) 1-tailed tests and resulting p-values were corrected across miRNAs using false discovery rates (FDR). Figs 1–3 showcases our differential expression analysis results for ER status. In Fig 1 scatter plot, we observe that the AQuN normalized dataset yields more significant results (lower Q-values) for most miRNAs (482/655). Fig 2 volcano plot illustrates that the increase in significance is not necessarily correlated with effect size (i.e. fold change), and that we gain confidence on lower effect sizes as anticipated by our power analysis (more details in the 4^th paragraph of the Discussion section). In Fig 3 cumulative distribution function (CDF) plot we depict the overall trend of increased statistical significance, contrasted by even lower statistical significance that would be obtained from performing the differential expression analysis on each dataset separately (shown as dashed lines). In addition, we present the CDF plots that would be obtained by (individually) applying four commonly used normalization methods (shown as dotted lines). Evaluated normalization methods include:

Mean ratio: scales each sample by $$ M\left(i,j\right)=\frac{M\left(i,j\right)}{AvgM(:,j)}$$.

Median subtraction: subtracts the median of each sample, then sets the minimum of each sample to the (global) minimum across samples. I.e.:

M(i,j) = M(i,j)−Medial M(:, j). min M(:,j)←min M(:,:)

Vanilla quantile: MATLAB’s implementation of Quantile Normalization also known as Quantile Standardization [27].

ComBat [28]: empirical Bayes batch effect mitigation employing a design matrix that includes dataset batching along with clinical labels and status of Tumor grade, Subtype, ER, PR, HER2 and TP53. We apply the QR decomposition [29] to mitigate any co-linearity in the design matrix.

Fig 1

Differential miRNA expression between ER positive and negative.

A scatter plot of differential expression p-values (-log₁₀, Wilcoxon Rank-sum) for the unnormalized (x) vs normalized (y) joint dataset. Title contains sample size details and dataset distribution.

Fig 2

Differential miRNA expression between ER positive and negative.

Volcano plot showing the fold change and corresponding Wilcoxon Rank-sum FDR corrected Q value ratio between the normalized and unnormalized datasets. Dashed arrow connects the unnormalized (gray circles) and normalized (red circles) results on a particular miRNA. High absolute values in X axis correspond to substantial difference in median expression between ER negative over ER positive samples (for a particular miRNA). High values in Y axis correspond to miRNAs that present substantial difference *after* normalization but not before. Low values in Y axis correspond to miRNAs that present substantial difference *before* normalization but not after. Vertical dashed lines represent a Fold change threshold of 2x (log₂(2) = 1) and horizontal dashed lines represent a Q-value threshold of 0.05 (-log₁₀(0.05)≅1.3).

Fig 3

Differential miRNA expression between ER positive and negative.

A CDF plot showing many more substantially differentially expressed miRNAs after normalization (red line) than before normalization (blue line), and substantially more than would be expected at random (compared to 20 random permutations of labels, dashed black lines). Also shown are dashed colored lines corresponding to each appropriate single-dataset Q values exemplifying the advantage of a joint-dataset analysis. Note that at Q = 10⁻¹⁸ we can find 10 miRNAs under AQuN but none under other normalization approaches or per dataset analyses.

In Figs 4 and 5 we demonstrate the impact of normalization on single miRNAs (hsa-miR-190b and hsa-miR-18a-5p, accordingly) across samples and on their differential expression in the context of ER status. This is done by detailing expression values for each sample in the joint dataset prior to (top row) and following (bottom row) AQuN normalization. We present the medians of each clinical group (dashed horizontal lines) and a breakdown of how samples of both clinical groups are distributed when sorted by value and when compared to a uniform null model. This provides a qualitative view of the effect normalization has on both individual samples and datasets in the context of the investigated differential expression. Previous studies [30] have shown hsa-miR-190b to be linked to ER status and further suggested its use as a potential biomarker. Similarly, hsa-miR-18a-5p is an oncogene and prognostic biomarker [31]. As we have shown in the volcano plot in Fig 2, hsa-miR-190b would not have been identified as differentially expressed in ER positive vs negative samples prior to normalization. Similar plots for the top 40 differentially expressed miRNA (post-normalization) are available compressed in S1 Data.

Fig 4

Visualizing expression of hsa-miR-190b across datasets and samples and in regard to estrogen receptor (ER) positive (pos) vs. negative (neg) differential expression.

(A, D) Expression values (log₂) of each sample before quantile normalization. Samples are ranked by ER status label, then by dataset and finally by ascending expression value. (A, B)-Unnormalized joint dataset. (C, D)-Normalized joint dataset. (B, C) Actual vs expected (via a uniform null model) rank distribution of ER negative (neg) vs positive (pos). Diagonal straight lines bounding a polygon represent a null uniform distribution of positive and negative samples (when ranked by expression value). The colored surface area represents the magnitude of deviation from a uniform distribution. The boundary of the surface is calculated by the cumulative number of ER negative (x axis) vs ER positive (y axis) samples in the ranked (descending) expression vector. Top-illustrating the rank distribution per-dataset (without normalization). Bottom-comparing the joint-dataset distributions when ranking before or after normalization.

Fig 5

Visualizing expression of hsa-miR-18a across datasets and samples and in regard to estrogen receptor (ER) positive (pos) vs. negative (neg) differential expression.

Caption description matches the one provided in Fig 4.

When inspecting the differential expression results of all normalization methods, the unnormalized data and each dataset separately, there are 33 unique miRNAs that are only shown as significantly (Q value<0.05) differentially expressed in ER positive vs ER negative as identified by our normalization method (S2 Fig, S2 Text, list in S6 Table). Contrastingly, other approaches yield far fewer significantly differentially expressed miRNAs. Of the 33 miRNAs uniquely detected by our method, we present four in Table 1 that have fold change greater than 0.15 (absolute log₂ > 0.15, which translates to > 10% change between median expression of ER positive and ER negative tumors).

Table 1

Top differentially expressed miRNA sorted by fold change.

miRNA	Fold Change (log2)	Q-value
hsa-miR-601	-0.18	0.048
hsa-miR-424-3p	-0.17	0.0003
hsa-miR-936	-0.15	0.027
hsa-miR-193b-5p	0.19	0.0002

We apply AQuN normalization on the joint dataset and not detected by other approaches. Fold change is defined as $$ {\log }_2\frac{MedianERpositive}{MedianERnegative}$$.

To study any breast cancer related functional significance of these top differentially expressed miRNAs we performed miRNA gain-of-function studies in the MCF-7 breast cancer cell line. Here, cell viability was measured as an endpoint after overexpression of the miRNAs. Indeed, one of the miRNAs, hsa-miR-193b-5p, showed a significant reduction in cell viability compared to miRNA negative controls (Fig 6). Furthermore, we looked into data from another functional experiment previously published [32] in the HER2 positive breast cancer cell line KPL4 and here we found that hsa-miR-193b-5p induced apoptosis (as measured by the levels of cleaved PARP), and downregulated the levels of HER2 and phosphorylated ERK upon overexpression. Altogether, these results suggest that miR-193b-5p may exert a tumor-suppressor function in breast cancer, both in an ER+ and a HER2+ context. Interestingly, the other miRNA originating from the same precursor, hsa-miR-193b-3p has been previously shown to directly target ESR1 mRNA and is thus a direct regulator of the ER [33].

Fig 6

Functional experiment results.

Breast cancer cell lines were transfected with miRNA mimics (20nM) and assayed for functional effects 72 hours after transfection. A) Cell viability measured in MCF7 breast cancer cells. B) Apoptosis measured by levels of cleaved PARP (cPARP), HER2 and phosphorylated ERK (pERK) protein levels measured in KPL4 cells. The dashed lines indicate cut-off points that were considered significant (see Materials and Methods). Asterisks denote significant effects. Original data from b) are taken from [32].

Further investigation of the three other top differentially expressed miRNAs shows prior evidence linking them to cancer. For example, hsa-miR-601 is a known prognostic marker and potential tumor-suppressor in breast cancer [34] and hsa-miR-936 was identified as a potential tumor-suppressor miRNA in ovarian cancer [35]. While these findings do not directly validate our findings in ER differential expression, they support the potential association of these miRNA through related mechanisms of cancer pathogenesis.

Joint analysis with mRNA data

A similar pipeline to the one described in subsection “Dataset pre-processing and coverage” was used to parse the Oslo2 cohort mRNA expression data, using Limma.

We wanted to assess the effect of AQuN normalization on the results of enrichment analysis as performed using both mRNA and miRNA data. To this end we first formed a ranked list of transcripts as follows. For each miRNA, μ, we ranked all mRNAs according to the (ascending) Spearman correlation between the miRNA expression pattern across the entire dataset and the mRNA expression pattern across the entire dataset (paired on matching samples). We denote the resulting ranked gene list, with μ as a pivot, as $$ {\mathcal{G}}_{\mu }$$.

Effect on gene target enrichment

For the first analysis we investigated the impact of AQuN normalization on correlations between miRNA and the expression levels of their expected mRNA targets in the Oslo2 dataset. We expect stronger negative correlation after normalization to direct gene targets. To validate this hypothesis, we applied a non-parametric, rank-based analysis using the MiTEA [36, 37] approach. MiTEA is used to evaluate the statistical association between $$ {\mathcal{G}}_{\mu }$$ and $$ {\mathcal{C}}_{\nu }$$, where $$ {\mathcal{C}}_{\nu }$$ is a ranked list of genes wherein the ranking is based on the affinity of the gene as a target candidate for the miRNA ν, taken from TargetScan [38]. A short overview of MiTEA’s algorithm is available in the Materials and Methods section.

We declare a matching if MiTEA returns a significant (≤0.001) P-value when ν = μ. To recapitulate, a matching occurs if the top of two lists of genes overlap to a high degree: the prominent predicted gene targets (by TargetScan) of miRNA μ and the list of genes ranked according to their sample-wise anti-correlation with their matched expression levels of miRNA μ. When applying this procedure on a non-normalized miRNA expression we find no matchings. When applying the same procedure on AQuN normalized data we find 6 matchings as detailed in Table 2. For each matched miRNA we also provide supporting evidence of several studies describing its role in breast cancer. We included an extension of this analysis across other datasets and normalization approaches in S2 Table.

Table 2

Resulting MiTEA matchings on normalized miRNA expression.

miRNA	P-value	Q-value	Corroborating studies
hsa-miR-29b	1.28E-08	1.73E-06	[39–41]
hsa-miR-106b	1.96E-06	1.11E-04	[42–44]
hsa-miR-200b	1.06E-04	5.54E-03	[45–47]
hsa-miR-30d	4.38E-04	1.19E-02	[48, 49]
hsa-miR-96	9.02E-05	1.53E-02	[50, 51]
hsa-miR-182	4.58E-04	4.43E-02	[52, 53]

P and Q values are color coded by magnitude where from green (more significant results) to red (less significant results). None of these statistically significant associations between pivot miRNAs and their targets is observed when using the raw, un-normalized data. Nor is any other matching miRNA target enrichment observed in the unnormalized data.

We show one such analysis in detail for hsa-miR-29b in Figs 7 and 8. Here we follow MiTEA’s approach to obtain a statistical assessment of target enrichment for μ = ν = hsa-miR-29b and $$ B=\left\{1,\dots ,\left|{\mathcal{C}}_{\nu }\right|\right\}$$ binary vectors $$ \mathcal{B}\left(\mu ,\nu ,B\right)$$. We present the results on various Bs and the optimal B* for both unnormalized and normalized miRNA expression.

Fig 7

Impact of normalization on the correlation between hsa-miR-29b expression and its in-silico predicted targets according to TargetScan.

(A) AQuN normalized vs Unnormalized (B) miRNA showing normalized is more negatively correlated to the prominent hsa-miR-29b targets in TargetScan as evident in stronger enrichment values.

Fig 8

Impact of normalization on the correlation between hsa-miR-29b expression and its in-silico predicted targets according to TargetScan.

Scatter plot of spearman correlation on normalized miRNA or unnormalized miRNA expression. If the target mRNA appears in TargetScan it is highlighted in orange. The marginal distributions are shown parallel to the axes and corresponding Kolmogorov-Smirnov test p-values display an overall lowered correlation for TargetScan candidates on normalized data.

Effect on gene ontology (GO) enrichment

We applied GOrilla [36] to identify gene ontology enrichment in $$ {\mathcal{G}}_{\mu }$$ on both unnormalized miRNA expression and on normalized miRNA expression. Given a ranked list $$ {\mathcal{G}}_{\mu }$$, GOrilla produces a binary vector $$ \mathcal{B}\left({\mathcal{G}}_{\mu },\omega \right)$$ for each gene ontology term, ω, in which a gene is labeled as binary ‘1’ if it belongs to ω. Next, GOrilla computes mHG p-values, correcting them across GO terms. Fig 9 is a scatterplot comparing between our results on unnormalized and normalized hsa-miR-29b lists. The findings from this analysis are in line with previous studies that have linked the miR-29 family with tumor growth and metastasis [40, 54, 55].

Fig 9

GOrilla enrichment analysis comparison of hsa-miR-29b correlation with gene expression before and after miRNA normalization with AQuN.

Showing scatter of GO term Q-values before and after AQuN. Red dots depicted with “≥P98” are above the 98^th percentile of Normalized–Unnormalized Q-values (-log₁₀) and green dots are for Unnormalized–Normalized. Right side panel shows a list of GO terms in the ≥P98 group.

Ethics statement

The Stavanger cohort was approved by REC Region West, approval number 2010/2014. By this approval, none of the patients were required to provide written informed consent to participate. All insights in a patient’s journal were monitored electronically, and all except the treating physician were required to state the reason why they needed to read that patient’s journal. This log was always open for the patient to view.

Overview

We used miRNA expression data from three previously published breast cancer datasets along with a newly released, fourth, miRNA dataset. These datasets were acquired from fresh-frozen material with different minimal number of tumor cells criteria, using different technologies and experimental protocols as overviewed in Table 3. In addition, we utilized mRNA expression to further investigate the effect of normalization using one of the cohorts. We examine miRNA normalization also in the context of jointly analyzing these measurements. Below we elaborate our considerations in the selections made during the normalization process and our means of providing evidence for validating these results.

Table 3

Technical details of platforms used for expression measurements for the four different cohorts.

Color code	Dataset	Manufacturer	Technology	Version	Accession number	Number of samples
	DBCG[25]–miRNA	Agilent	Human miRNA Microarray Kit	(V2 G4470B) design id 019118	GSE46934	149
	Oslo2[15]–miRNA	Agilent	Human miRNA Microarray Kit (V2)	v14 Rev.2 design id 029297	GSE81000	425
	Oslo2[15]–mRNA	Agilent	SurePrint G3 Human GE 8x60K Microarray	(Probe Name Version) 028004	GSE80999	381
	Micma[11]–miRNA	Agilent	Human miRNA Microarray Kit	(V2 G4470B) design id 019118	GSE19536	101
	Stavanger–miRNA	Exiqon	miRCURY LNA Array	v.11.0		109

Datasets are color coded consistently throughout the paper. miRNA expression colors are highlighted compared to mRNA measurements.

Dataset pre-processing and coverage

Each miRNA dataset is read from a single-channel image analysis output file acquired from their corresponding GEO repositories (referenced in Table 3) and preprocessed in R using the Limma [60] package. We note that while Stavanger (Exiqon) data contains a pooled-reference second channel, this measurement is not utilized in our analysis (further discussed in S1 Text). Initially, control probes are removed, and the data is corrected by background intensity normalization [61]. Same-probe replicate samples are replaced by their median value. Probe ids are mapped to their corresponding miRbase v22 accession using miRBaseConverter [62]. Missing or deleted accession IDs are discarded. Multiple probes that map to the same miRNAs are again replaced by their median value. Next, we apply arrayQualityMetrics [63] (resulting Quality Control reports are available compressed in S2 and S3 Data) and filter out samples that are marked as outliers by all three outlier detection criteria (L₁-Distance between arrays, Boxplot, MA plot). We thereby filtered out 6, 30, 12 and 2 outliers from DBCG, Oslo2, Micma and Stavanger, respectively. Next, we apply minimum subtraction to avoid log scaling issues with negative numbers where applicable. The joint dataset table is then compiled by applying a “full outer-join” relational operation on the miRbase accession IDs as key. The resulting miRNA cross-dataset table is visualized in Fig 12 (and available in S3 Table as raw data and S4 Table as normalized data with corresponding clinical labels in S5 Table).

Fig 12

Overview of the miRNA coverage in the dataset.

Each row represents one miRNA. Each entry represents the intensity (log₁₀) in a specific sample. Dashed vertical lines separate between samples from the four datasets. Dashed horizontal lines separate between groups of miRNAs by their dataset availability. Blank (white) entries correspond to miRNAs that are missing from a dataset.

Batch effects in joint data

We tested for rank-order consistency of miRNA expression in pairs of datasets (Fig 13). For each miRNA we take the median of its expression, or similarly, intra-sample percentile, across all samples belonging to the same dataset. We display the resulting values for each pair of datasets in a scatterplot matrix considering the miRNAs (n = 655) present in all four cohorts. This analysis shows that the Stavanger data appears to behave differently, presumably due to its fundamentally different measurement technology (Exiqon LNA—Locked Nucleic Acids vs Agilent Microarray).

Fig 13

Scatterplot matrix of quantile normalized data showing miRNA expression reproducibility across dataset pairs.

Each subplot depicts a pair of datasets. In the upper-diagonal-subplots, each point corresponds to a single miRNA’s median (across samples) rank (intra-sample) in each dataset. Similarly, the bottom-diagonal shows median log2 expressions in place of ranks. A second degree polynomial curve is fitted and prediction intervals at confidence level 0.8 are plotted as dashed lines. Spearman correlation is given for each subplot. Figures at the diagonal show percentile plotted against expression and a circle represents the dataset colorcode as related to other figures in the paper.

We further visualize the batch-clustering behavior of the unnormalized joint dataset in Fig 14. On the left panel (A) we present hierarchical clustering of the data. Edges of sub-trees in the dendrogram are color-coded by the dataset when all leaves in the subtree belong to samples from the same original dataset. We observe a visual clustering of colors, especially evident for yellow (Stavanger) being clustered as an outgroup. In the middle panel (B) we show a silhouette plot, depicting the clustering consistency according to dataset. We observe a substantial portion of samples that are well assigned to their cluster with large silhouette values, and only a small portion are mis-assigned, again showcasing how batch effects dominate sample pairwise-distance pattern behavior. Finally, on the right panel (C) we present a visualization of the sample-wise pairwise Euclidean distance matrix with dashed lines separating between samples of the same dataset. The block-diagonal structure that evidently results from coloring according to distances corresponds well to the dashed lines separating samples from different datasets. This analysis demonstrates the prevalence of batch effects in the joint datasets.

Fig 14

Visualizing batch effects in the combined cross-tech miRNA dataset considering the unnormalized data.

(A) Dendrogram with edges colored by dataset. Note that the tree root is outside the displayed axis range. (B) Silhouette plot [67] showing that most samples cluster according to the dataset they originate from. (C) Pairwise Euclidean distances showing a block structure that agrees with the sample dataset of origin.

Adjusted quantile normalization (AQuN)

In this section we describe our quantile-normalization-based strategy for analyzing combined cross-technology miRNA datasets.

Let M be a batch collected, joint dataset. M∈R^n×m where M(i,j) is the log measured intensity value of miRNA i in sample j. We note that the j-th column of M, denoted M(:,j), corresponds to the j-th sample and that the i-th row, M(i,:) corresponds to the i-th miRNA in the joint dataset.

Define β(M(:,j)) = k as the experiment batch id during which sample j was collected.

We note the following distinction between missing values in M:

Let MFP(i,j) = 1 if miRNA i is missing from platform β(M(:,j)) and MFP(i,j) = 0 otherwise (indicates if i is missing in the platform j was measured in).

	Adjusted Quantile Normalization (M):
1.	$$	Jitter M to break rank ties.
2.	Let P(i,j) = the percentile of $$ within $$.	nans are ignored in percentile computation. Note: P(i,j)∈[0,100]
3.	$$	Transforms values to the cross-sample-median of the corresponding per-sample-quantile.
4.	Q(i,j) = nan if MFP(i,j) = 1

A description of this process in words is that it replaces all present expression values with the corresponding median value of all samples within the same percentile. The underlying assumption is that a measured expression is volatile due to technical differences and measurement noise, however, (sample-based) percentiles are assumed to be stable up to the biological differences between samples. In addition—platform coverage differences are addressed.

The overall impact of applying AQuN to the distribution of expression values and to quantified batch effects as measured by the silhouette coefficient is further presented in S1 Fig.

Functional experiments

Functional experiments were performed as previously described [32, 33] with the breast cancer cell lines MCF7 and KPL-4. The lysate microarray data measuring apoptosis in the form of cleaved PARP (cPARP), HER2 and phosphorylated ERK (pERK) protein levels after 72 hours were previously published (data taken from S2 Table of referenced paper and provided as S7 Table herein) [32]. Values ±2 × standard deviation (SD) were considered as significant, which corresponded to a threshold of |1.96|. For the cell viability data, MCF7 cells were transfected with the Dharmacon miRIDIAN microRNA mimic library v.10.1 (20 nM) and incubated for 72 hours. The cell viability was measured with CellTiter-Glo assay (Promega Corp, Madison, WI, USA) according to manufacturer’s protocol. The experiments were done with two biological replicates. The data were normalized by a Loess method [64] and log2-transformed. Values ±2 × SD, were considered as significant, which corresponded to a threshold of |0.2|. In both experiments the average of two different miRNA mimic controls from two replicates was used as negative controls (miRIDIAN microRNA Mimic Negative Control #1 from Dharmacon and pre-miR negative control #2 from Ambion). The transfection efficiency of miRNA mimics has been determined previously [33].

MiTEA algorithm overview

Briefly, for each prefix $$ {\mathrm{\Pi }}_{\mathrm{B}}\left({\mathcal{C}}_{\nu }\right)$$ of B most-prominent candidate targets in $$ {\mathcal{C}}_{\nu }$$, MiTEA produces a binary vector, $$ \mathcal{B}\left(\mu ,\nu ,B\right)$$, such that, g_i, the i-th gene in $$ {\mathcal{G}}_{\mu }$$ is labeled “1” if and only if it is in the candidate prefix, i.e. $$ g_i\in {\mathrm{\Pi }}_{\mathrm{B}}\left({\mathcal{C}}_{\nu }\right)$$. MiTEA then computes an approximate minimum hypergeometric (mHG [36, 57]) P-value to quantify whether the B proposed targets are enriched at the top of the $$ {\mathcal{G}}_{\mu }$$ list or not. Finally–MITEA applies an FDR correction (using the Benjamini-Hochberg procedure [65]) across evaluated νs and reports the set of miRNAs associated with the ranked target list $$ {\mathcal{G}}_{\mu }$$ and their associated Q-values.