Detection of miRNA regulatory effect on triple negative breast cancer transcriptome

Functional miRNAs in breast cancer

The analyzed dataset includes 129 tumour and healthy breast tissues for which both miRNA and mRNA expression profiles are available (see Material and Methods). In this study, we used sequence-based target sets obtained from TargetScan version 6.2 [8, 25], a widely used algorithm which takes into account sequences that match the seed region of each miRNA and evaluate their conservation in several species. The confidence of target predictions is calculated as described in Methods. We include in our study 394 miRNAs for which both expression data and TargetScan predictions are available.

We identified 136 miRNAs as potential regulators with FDR < 0.1 (results are reported in Additional File 1, Table S1). Among the top significant miRNAs, we found several ones that are known to function as oncomirs, such as the members of miR 17-92 cluster and its paralog cluster miR-106b-25, miR-15 and miR-16 [26].

To validate these results, we applied our algorithm to an independent dataset of 521 healthy and cancerous breast tissue samples from the TCGA project (see Material and Methods). This study includes 260 miRNA for which both expression data and TargetScan predictions are available. Of these, 142 miRNAs are identified as regulators with FDR < 0.1. The intersection between results obtained in both datasets contains 44 miRNAs (see Figure 2a and Table 1) which corresponds to a statistically significant overlap (hypergeometric p − value < 0.05).

For a limited subset of 7 commonly identified miRNAs, the enrichment score obtained for the Maire dataset has opposite sign compared to that obtained for the TCGA dataset. In these cases, predicted targets are enriched among genes whose expression is positively correlated with the expression of the miRNA. The observation of miRNAs positively correlated with their predicted targets is in agreement with analogous integrated analysis of miRNA-mRNA correlation in tissue samples [27]. An interesting hypothesis suggests that this effect can be explained by common transcriptional regulation conferring robustness to gene expression program and ensuring tissue identity. Consistently, architectural features of the mammalian miRNA regulatory network reveal that the coordinated transcriptional regulation of a miRNA and its targets is an abundant motif in gene networks [28–30]. In our analysis, we consider both positive and negative correlation of predicted targets as a good evidence to infer miRNA regulation. We investigated whether the correlation sign of the miRNA expression with that of its targets is associated to a different proportion of correlated targets or a different correlation strength. For each significant miRNA we extracted the leading-edge subset of genes, corresponding to those genes in the set S _m that appear in the ranked list before the point where the running sum achieves the ES (see Figure 3). This can be interpreted as the core of a gene set that accounts for the enrichment signal. We plotted the size of the leading-edge subset and the maximum correlation value as function of the correlation sign. Despite the size of the leading-edge subset and its maximum correlation value do not differ significantly according to the sign of the correlation, results show a trend toward stronger correlations of negatively correlated leading targets (Additional File 2, Figure S1).

Triple negative breast cancer specific study

A group of 20 miRNAs is identified specifically in the TNBC study. Remarkably, some of these are already associated with aggressive breast cancer: miR-29a, miR-29c and miR-148a have been shown to be downregulated and associated to aberrant hypermethylation in basal-like cell line [32], miR-27a is involved in endothelial differentiation of breast cancer in a basal-like cell line [33] and in the MDA-MB-231 basal-like cell line [34] and miR-139-5p is described as a regulator of breast cancer cell motility and invasion [35].

Identification of biological processes targeted by miRNAs

To assess the biological relevance of miRNAs identified as potential regulators, we investigated whether the genes that account for the enrichment signal of a given miRNA participate in the same cellular process or signalling pathway. The assumption that some miRNAs downregulate a group of genes participating in the same pathway is supported by multiple experimental studies [36–38]. Based on this hypothesis, for each significant miRNA we tested the leading-edge subset of genes for functional enrichment using curated annotation databases, such as Gene Ontology [39], KEGG [40], BioCarta [41], Reactome [42] and ACSN [43].

For miRNAs identified in both Maire and TCGA datasets, the corresponding subsets of leading-edge targets were analyzed separately. Interestingly, these subsets of genes display highly significant overlap and similar enriched categories, supporting the relevance of miRNA regulatory role in breast cancer. We report in Additional File 3, Table S2 the complete list of enriched categories with Bonferroni corrected p-values below 10−2. The most enriched categories include breast cancer related processes like cell substrate adherens junction, regulation of cell migration, nuclear pore complex and integrin pathway.

Computational framework

The ranking scheme and the enrichment score definition are described as follows.

where r _mj is the correlation between the expression of gene j and the expression of miRNA m, < s _m · |r _m | ^α > is the average of prediction confidence-weighted correlations for the set S _m . The parameter α controls the contribution of the correlation r _mj such that correlation values can contribute linearly (α = 1) or non-linearly in the running sum. Similar to the original GSEA algorithm, in this study, we set α = 1. According to this equation, the running sum is incremented by value (s _mj · |r _mj | ^α − < s _m · |r _m | ^α >) when encountering a gene in S _m and decreased by < s _m · |r _m | ^α > when not.

The statistical significance of the enrichment score ES(S _m ) is assessed by a permutation based p-value. The enrichment score of the randomly shuffled list G _m is computed for N _p =1000 permutation rounds and the empirical null distribution of ES(S _m ) is generated. An empirical p-value PV(ES(S _m )) is estimated for the positive and the negative region of the distribution by the proportion of permutations that result in larger ES(S _m ) than originally observed (or lower ES(S _m ) for the negative region). Once PV(ES(S _m )) for a miRNA m is below a fixed threshold, the regulatory activity of miRNA m on the transcriptome is inferred. The PV threshold is set according to False Dicovery Rate (FDR) obtained by the Benjamini-Hochberg method [44]. In our study, the PV threshold was set according to FDR < 0.1.

Tissue collection

Healthy samples from mammary plastic surgery and tumor samples were obtained from patients treated at the hospital of Institut Curie (Biological Resource Center, Paris, France) as described previously [22, 23]. Experiments were performed in agreement with the Bioethic Law No. 2004-800 and the Ethic Charter from the French National Institute of Cancer (INCa), and after approval of the ethics committee of our Institution. From the initial dataset, we retained only the samples for which both mRNA and miRNA data were available, numbering to 129, including 11 healthy breast tissue samples, 37 TNBC, 28 ER-/HER2+, 24 Luminal A and 29 Luminal B breast tumour samples as characterized by immunohistochemical staining.

MiRNA expression data

Samples were hybridized on the Agilent miRNA microarray kit (V3). One hundred ng of total RNAs was hybridized to the microarrays according to the manufacturer's instructions. Hybridized microarrays were scanned with a DNA microarray scanner (Agilent G2565BA) and features were extracted using the Agilent Feature Extraction image analysis tool with default protocols and settings. Data were first transformed using the reverse hyperbolic sine function and quantile normalized. The data were then corrected for a hybridization batch effect using a linear model including the hybridization series as a fixed effect. Next, probes with negative intensity values in 95% or more of the arrays were discarded, leaving 516 miRNAs for analysis. When technical replicates for a sample were present, they were subsequently averaged. Two samples displayed an outlier behavior evident from Principal Component Analysis (data not shown) and were discarded from this analysis. The dataset can be downloaded here at the following address: http://bioinfo-out.curie.fr/projects/targetrunningsum.

Gene expression data

For the protein-coding transcriptome, the data from Affymetrix U133plus2 arrays was processed as described [22]. In summary, we used the brainarray HGU133Plus2-Hs-Entrez version 13 custom chipset definition file [45], data were then normalized with GC-RMA [46], technical batch artefacts were corrected using a linear model, and genes with noise-level expression in 95% or more arrays were filtered out leaving 11543 probesets each corresponding to a unique Entrez Identifier.

MiRNA and gene expression data from TCGA

We conducted our study on the publicly available data of tumour and healthy breast tissues from TCGA described in [24]. We selected 500 tumors and 21 tumor-adjacent normal breast tissue samples for which both mRNA and miRNA data were available. Of the 500 tumors, 82 were assigned to the TNBC subtype according to ER, PR and HER2 negative status. In this dataset mRNA expression levels were determined by Agilent custom 244K whole genome microarrays and miRNA abundance was measured by Illumina sequencing technology. Level 3 released data contain quality controlled and processed data done by Broad Institute's TCGA workgroup with expression call for genes per samples. Gene level expression data were normalized by using Robust Multi-array Average (RMA) and expression values were gene centered. MiRNA expression was given as read counts normalized to relative read frequency in each sample. Detected miRNAs were defined as having more than 10 reads in at least 10% of the samples, leaving 332 miRNAs for analysis.

Target predictions by TargetScan

In the TargetScan algorithm, the prediction score of a miRNA binding site depends on the level of conservation and sequence context criteria such as the distance of the target from the 3'UTR end and the AU composition of the flanking area. For each miRNA, we take as prediction confidence weight of its targets the total context scores generated by the algorithm for 3'UTRs aggregating the binding site scores. By construction, the total context scores range between -1 and 0. When total context score for multiple 3'UTR isoforms are predicted, we take the total context score of the longest 3'UTR isoform.