Many tools have been developed to profile microRNA (miRNA) expression from small RNA-seq data. These tools must contend with several issues: the small size of miRNAs, the small number of unique miRNAs, the fact that similar miRNAs can be transcribed from multiple loci, and the presence of miRNA isoforms known as isomiRs. Methods failing to address these issues can return misleading information. We propose a novel quantification method designed to address these concerns.
We present miR-MaGiC, a novel miRNA quantification method, implemented as a cross-platform tool in Java. miR-MaGiC performs stringent mapping to a core region of each miRNA and defines a meaningful set of target miRNA sequences by collapsing the miRNA space to “functional groups”. We hypothesize that these two features, mapping stringency and collapsing, provide more optimal quantification to a more meaningful unit (i.e., miRNA family). We test miR-MaGiC and several published methods on 210 small RNA-seq libraries, evaluating each method’s ability to accurately reflect global miRNA expression profiles. We define accuracy as total counts close to the total number of input reads originating from miRNAs. We find that miR-MaGiC, which incorporates both stringency and collapsing, provides the most accurate counts.
MicroRNAs (miRNAs) are endogenous small (~ 23 nt) RNA molecules that contribute to post-transcriptional regulation of target messenger RNAs (mRNAs) in plants and animals [1, 2]. In recent years, many tools have been developed to estimate miRNA expression from small RNA-seq data. These include CAP-miRSeq , Chimira , CPSS , iSRAP , miRanalyzer , the miRDeep2 quantifier , miRExpress , miRge , miRNAKey , mirTools , Oasis , omiRAs , Shortran , and sRNAbench . Table 1 summarizes these methods. In a typical workflow, the read counts form the foundation for downstream analyses such as differential expression and co-expression analysis. Therefore, accurate expression quantification is essential for the validity of downstream results.
The effectiveness of quantification methods may be affected by three issues particular to miRNAs. One issue involves mapping accuracy. The small size of miRNA molecules leads to short sequencing reads after adapter removal. Short reads are less likely to be aligned uniquely to the genome ; this issue could be compounded by individual genetic variation at the endogenous locus producing the read . The second issue involves challenges of functional interpretation. Identical or near-identical miRNAs are often transcribed from multiple genomic loci [19, 20]. So as not to introduce count bias, quantification methods must deal with reads that map ambiguously to multiple loci or miRNA sequences. In addition, there are many fewer unique miRNA molecules than large RNAs. Normalization methods such as total read count or quantile normalization are less robust with fewer features and highly skewed distributions. Therefore, the handling of multi-mapped reads can have a larger impact on normalized counts for miRNAs compared to larger RNAs. Third, isomiRs—miRNA variants that can be expressed in a cell type specific manner—present a challenge for mapping and functional interpretation. Research suggests that the three main classes of isomiRs (5′ isomiRs, 3′ isomiRs, and polymorphic isomiRs) may have differing functional consequences [21, 22]. The question of whether isomiRs should be counted and, if so, which ones should be merged with their parent miRNA for expression analysis, is nontrivial and should be addressed by quantification methods.
Methods that fail to adequately address these issues can return misleading quantification results. We examined the accuracy of several published methods as well as a novel quantification pipeline that incorporates stringent mapping and collapsing of the miRNA space into meaningful functional units.
We designed a quantification method with the following objectives: (1) perform highly stringent mapping to a core region of miRNA sequences, minimizing the number of ambiguous mappings, and (2) perform collapsing to associate reads with functional classes of miRNAs instead of individual annotated miRNAs. Functional classes of miRNAs, subsequently referred to as “functional groups”, are defined by the user to be groups of miRNAs that are considered equivalent in the context of the study goals. For instance, if the study aims to address binding of target mRNAs, families of highly similar miRNAs that bind the same targets can be considered equivalent. This consideration allows reads to be counted at most once per functional group; counts are then returned at the group level. We implemented a pipeline, miR-MaGiC, that incorporates these features. For details of the software and workflow, see Additional file 1: Additional material and Figure S1.
We tested miR-MaGiC and several publicly available methods on 210 mouse brain small RNA-seq libraries. This dataset was chosen due to the large number of samples and high sequencing depth, making it a valuable test case for comparing methods, while the variability in proportion of miRNA reads between libraries provided an interesting testing scenario. We ran 7 quantification schemes for each library: iSRAP , the miRDeep2 quantifier , miRge , a modified version of miRge, and three collapsing conditions for miR-MaGiC. Our modified version of miRge removed its final round of alignments to mature miRNAs, a highly permissive alignment step that allowed up to two mismatches per read; we suspected that this step may introduce noise to the counts. See Additional file 1: Table S1 and Additional material.
To evaluate the methods, we reasoned that methods which correctly handle the issues particular to miRNA quantification should return total counts that reflect the number of reads originating from miRNAs in the input library. We estimated the number of miRNA reads in each library as the number of adapter-clipped reads between 19 and 23 nucleotides in length; 95% of miRNA loci and 91% of unique mature miRNAs in miRBase fall in this length range. The libraries each had between 50% and 72% of reads in this range. We examined how well each method reflected this estimated number of input miRNA reads in terms of total output read count, calculating the mean squared error between the estimated number of input miRNA reads and the output total counts. A lower score would indicate more accurate counts and therefore less distortion and bias introduced during normalization by the method-dependent total count.
Due to different implementation choices, the methods systematically return different levels of total absolute counts (Fig. 1). The miRDeep2 quantifier returns the highest counts because it first matches mature miRNAs to precursors in a many-to-many mapping, then counts every instance of a read matching one of these mature miRNA/precursor pairs. As expected, miR-MaGiC returns reduced total counts when functional group collapsing is performed, as opposed to no collapsing. Because the read counts for miRNAs are right skewed (Additional file 1: Figure S2), double counting in any of the highly expressed miRNAs can dramatically change the total read count. See Additional file 1: Additional material and Figure S3 for a case study of miRNAs that are treated differently by different methods.
Comparing miR-MaGiC to published software, miR-MaGiC with collapsing by functional group showed the best accuracy (Fig. 2). The least accurate method is the miRDeep2 quantifier, probably due to double counting reads that map to multiple precursors. The closest method to miR-MaGiC is miRge, which also incorporates collapsing but uses permissive mapping. As expected, miR-MaGiC with no functional group collapsing is less accurate than with collapsing. When we modified the miRge code to remove the final round of highly permissive alignments, performance improved dramatically and the method gained a slight advantage over miR-MaGiC with collapsing. One possible explanation for why the published version of miRge is less accurate than the more stringent modified version is that the permissive alignment step allows some non-miRNA reads to be mapped to miRNAs.
We have proposed a quantification method, miR-MaGiC, that addresses several issues particular to miRNAs, including their small size, low complexity, family structure, and isoforms. miR-MaGiC uses stringent mapping to reduce noise associated with the small size and low complexity of miRNAs, while allowing for uncertainty at the endpoints of reads and miRNAs. Final counts are returned at the group level instead of the individual miRNA level. Recommended group tables are provided for common species on the miR-MaGiC web page, https://github.com/KechrisLab/miR-MaGiC.
We tested miR-MaGiC as well as three published methods on a set of 210 small RNA-seq libraries. We evaluated the faithfulness of the final total counts to the original number of miRNA reads per library. Importantly, we found that methods which specifically address the above issues produced the greatest accuracy in overall counts. The novelty of miR-MaGiC is the combination of stringent mapping to a core region of each miRNA and collapsing by functional group.
To evaluate this combination of features we tested miR-MaGiC with and without collapsing, observing that collapsing in fact improves accuracy. Regarding mapping stringency, the published version of miRge, which performs collapsing, performed poorly according to our accuracy metric, but we suspected this may be due to over-permissiveness of one of its alignment steps. Once we modified this detail, miRge emerged as comparable to miR-MaGiC, with a slight advantage in accuracy. In summary, when methods use one feature but not the other (i.e., miR-MaGiC_noCollapse and miRge in Fig. 2), or neither feature (i.e., iSRAP and miRDeep2 in Fig. 2) there is a notable drop in accuracy.
Our analysis of miRge indicated that more noise than signal is introduced if methods try to capture isomiRs simply by allowing more mismatches. miR-MaGiC uses stringent mapping to reduce noise associated with the small size and low complexity of miRNAs. This decision effectively causes 5′ and 3′ isomiRs to be merged with their parent miRNA while discarding polymorphic isomiRs. 3′ isomiRs are the most common class of isomiR and are thought to be largely functionally redundant, while 5′ and polymorphic isomiRs are less common but can affect target binding [21, 22]. Therefore, miR-MaGiC merges most functionally redundant miRNA isoforms with their parent miRNA while also possibly including 5′ isoforms that may affect function. This decision has the effect of including the largest class of isomiRs which are currently believed to be largely functionally redundant while excluding polymorphic isomiRs which may have distinct functions.
In this work, we examined accurate quantification of miRNA expression based on sequencing. Several issues particular to miRNAs can affect the accuracy of quantification methods based on small RNA-seq. These issues include the small size of miRNAs, the low complexity of the overall repertoire of miRNAs, the fact that highly similar miRNAs can be processed from different genomic loci, and the presence of isomiRs. Furthermore, it is important that quantification be performed at an appropriate level of granularity to be functionally meaningful. Implementation choices at the quantification step can have a significant impact on common downstream steps such as normalization and interpretation of expression results. When counts are split over multiple features, the multiple testing burden is increased and statistical power is reduced. In addition, the relatively low complexity of the miRNA repertoire means that a handful of highly expressed miRNAs can have an impact on the library size used for normalization.
Our work demonstrates the importance of identifying the most meaningful unit of information when studying miRNA expression. We find that results are most accurate when we associate each read with one meaningful unit such as a miRNA family. To accomplish this, our proposed method, miR-MaGiC, looks for a stringent match to one or more members of the family and then ignores which member(s) it matched and reports results for the family. The mapping is stringent in one sense, but also flexible at the ends of each miRNA, as these can be affected by isomiRs or artifacts in the reads. The most meaningful level of granularity for a particular study may vary. We therefore recommend that investigators understand the implementation details of various quantification methods and choose a method that will return the most meaningful expression profile for their study.
Materials and methods
Known miRNAs and creation of individualized miRNA sequences
We used the mouse miRNA database in miRBase version 21 . See Additional file 1 for details.
Defining functional groups of miRNAs
Our pipeline, miR-MaGiC, counts mappings of reads to functional groups of miRNAs instead of individual miRNAs. We evaluated three different groupings of miRNAs. The first was no collapsing by functional group. The second combined miRNAs with the same miRBase accession number (“MIMAT” number) before an underscore. The final grouping combined miRNAs with the same core number, letter (if applicable), and 3p/5p identifier. See Additional file 1 for details.
Test with publicly available software packages
We chose publicly available methods to include in our comparison based on several criteria: (1) ability to be run in batch jobs on a Linux cluster, (2) success of installation and execution on our Linux environment, and (3) methods representing a variety of quantification strategies. These criteria led to choosing iSRAP , the miRDeep2 quantifier , and miRge . 210 mouse whole brain small RNA-seq libraries were analyzed. Run details are in Additional file 1: Table S1 and Additional material.
Our analysis demonstrates that for short sequences from a low-complexity repertoire, a high level of mapping stringency is important for minimizing noise. However, a limitation of this high stringency is that errors in reads or individual variation in miRNAs could lead to incorrectly missed read mappings, i.e., an increase in false negative mappings. Another limitation is that miR-MaGiC only generates counts and does not perform analyses such as normalization and differential expression, in contrast to other small RNA-seq analysis tools that perform multiple analyses in a pipeline fashion. Nonetheless, the resulting miR-MaGiC quantification is easily plugged into other downstream analyses.
mean squared error
reads per million
Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116:281–97.
Sun Z, Evans J, Bhagwate A, Middha S, Bockol M, Yan H, et al. CAP-miRSeq: a comprehensive analysis pipeline for microRNA sequencing data. BMC Genomics. 2014;15:423. https://doi.org/10.1186/1471-2164-15-423.
Hackenberg M, Rodríguez-Ezpeleta N, Aransay AM. miRanalyzer: an update on the detection and analysis of microRNAs in high-throughput sequencing experiments. Nucleic Acids Res. 2011;39:W132–8. https://doi.org/10.1093/nar/gkr247.
Friedländer MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40:37–52. https://doi.org/10.1093/nar/gkr688.
Wang W-C, Lin F-M, Chang W-C, Lin K-Y, Huang H-D, Lin N-S. miRExpress: analyzing high-throughput sequencing data for profiling microRNA expression. BMC Bioinform. 2009;10:328. https://doi.org/10.1186/1471-2105-10-328.
Wu J, Liu Q, Wang X, Zheng J, Wang T, You M, et al. mirTools 20 for non-coding RNA discovery, profiling, and functional annotation based on high-throughput sequencing. RNA Biol. 2013;10:1087–92. https://doi.org/10.4161/rna.25193.
Müller S, Rycak L, Winter P, Kahl G, Koch I, Rotter B. omiRas: a Web server for differential expression analysis of miRNAs derived from small RNA-Seq data. Bioinformatics. 2013;29:2651–2. https://doi.org/10.1093/bioinformatics/btt457.
Barturen G, Rueda A, Hamberg M, Alganza A, Lebron R, Kotsyfakis M, Shi B-J, Koppers-Lalic D, Hackenberg M. sRNAbench: profiling of small RNAs and its sequence variants in single or multi-species high-throughput experiments. Methods Next-Generation Seq. 2014;1:21–31.
Li W, Freudenberg J, Miramontes P. Diminishing return for increased mappability with longer sequencing reads: implications of the k-mer distributions in the human genome. BMC Bioinform. 2014;15:2. https://doi.org/10.1186/1471-2105-15-2.
PHR carried out the analysis and took the lead with writing the manuscript. BV, WS, and PDR provided feedback on the development of the analysis and interpretation of results. RR and RD generated the genetic data for the study. PHR, LS and KK conceived the analysis and LS and KK were in charge of overall direction and planning. All authors read and approved the final manuscript.
We thank Spencer Mahaffey, University of Colorado Skaggs School of Pharmacy, for his assistance in making the sequencing data available over the Web.
The authors declare they have no competing interests.
This work was supported by the National Institute on Alcohol Abuse and Alcoholism of the National Institutes of Health (Grant Numbers R01AA021131, R01AA016957). WS acknowledges support from National Library of Medicine Institutional Training Grant (Grant Number T15LM009451).
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Authors and Affiliations
Department of Biostatistics and Informatics, Colorado School of Public Health, Aurora, CO, 80045, USA
Pamela H. Russell, Pratyaydipta D. Rudra & Katerina Kechris
Computational Bioscience Program, University of Colorado, Aurora, CO, 80045, USA
Department of Molecular, Cellular, and Developmental Biology, University of Colorado, Boulder, CO, 80309, USA
Department of Pharmaceutical Sciences, University of Colorado Skaggs School of Pharmacy and Pharmaceutical Sciences, Aurora, CO, 80045, USA
Richard Radcliffe & Laura Saba
Center for Genes, Environment and Health, National Jewish Health, Denver, CO, 80206, USA
Additional materials, methods, figures, and tables.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.