miR-MaGiC improves quantification accuracy for small RNA-seq

Results

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 [6], the miRDeep2 quantifier [8], miRge [10], 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.

Conclusions

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.

Discussion

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