The sufficient minimal set of miRNA seed types.

MOTIVATION: Pairing between the target sequence and the 6-8 nt long seed sequence of the miRNA presents the most important feature for miRNA target site prediction. Novel high-throughput technologies such as Argonaute HITS-CLIP afford meanwhile a detailed study of miRNA:mRNA duplices. These interaction maps enable a first discrimination between functional and non-functional target sites in a bulky fashion. Prediction algorithms apply different seed paradigms to identify miRNA target sites. Therefore, a quantitative assessment of miRNA target site prediction is of major interest.
RESULTS: We identified a set of canonical seed types based on a transcriptome wide analysis of experimentally verified functional target sites. We confirmed the specificity of long seeds but we found that the majority of functional target sites are formed by less specific seeds of only 6 nt indicating a crucial role of this type. A substantial fraction of genuine target sites arenon-conserved. Moreover, the majority of functional sites remain uncovered by common prediction methods.

1 INTRODUCTION

The relation between miRNAs and their targets in higher eukaryotes is part of the highly complex gene regulation network. To unravel the functional specific interactions, the available information on the interaction of the short RNAs as presented by the RISC complex and their mRNA counterparts is insufficient to reliably predict all functional pairs modulating translation and mRNA decay (Bagga ; Guo ; Lim ).

The basic prerequisite for miRNA targeting in metazoans is a short perfect match complemented by imperfect matches in close vicinity. This region is called the seed sequence and is considered to be a 6–8 nt long substring within the first 8 nt at the 5′-end of the miRNA (Lewis ). It is regarded to be the most important feature for target recognition by miRNAs in mammalians (Bartel, 2009; Nielsen ).

Naturally, merely seeking for short sequence matches yields a plethora of putative target sites with a large fraction of false positives. To dodge a priori the majority of false positives in silico miRNA target site prediction approaches concentrate on the subset of target sites equipped with long perfect seed matches. In addition, several miRNA targeting determinants beyond the seed have been proposed (Grimson ; Hausser ; Kertesz ) to extract authentic target sites from the set of seed matches. A common strategy to increase specificity is to require conservation of the seed match. But there is evidence that non-conserved miRNA targeting is even more widespread (Baek ; Farh ).

To date the effect of different types of seed matches has been assessed by means of signal-to-noise ratio (Lewis , 2005), degree of mRNA (Grimson , Nielsen ) or protein repression (Baek , Selbach ). Based on that, a set of canonical seed types that differ in abundance and intensity of the regulatory effect has been defined (Bartel, 2009). Recent experimental approaches allow for the identification of Argonaute (Ago)-miRNA:mRNA ternary complexes using an in vivo cross-linking protocol and subsequent high-throughput sequencing (Chi ; Hafner ). Chi et al. (2009) analyzed miRNA:mRNA interactions in Mus musculus neocortex tissue samples and published an interaction map containing a set of verified target sites in the transcriptome of the murine brain.

Complementing previous studies, we determined seed types using functional target sites of the interaction map. We identified a minimal and sufficient set of six seed types. The precise mapping of Ago footprints allowed us to distinguish between miRNA:target and higher resolved miRNA:target site interaction. We quantified the impact of individual seed types on recall and specificity. Additional target site conservation analyses revealed short seeds to be less conserved than long seeds.

2 MATERIALS AND METHODS

Chi et al. (2009) provided a miRNA:mRNA interaction map that contains the absolute chromosomal positions of sites full complementary to miRNA seeds (murine genome assembly of 2006). These sites are located almost at the center of an average Ago-mRNA footprint. This is a defined region of mRNA complexed with Ago determined by Ago-mRNA clusters, where Ago bound within 62 nt of cluster peaks ≥ 95% of the time. For each chromosomal coordinate, we determined the longest protein-coding mature mRNA transcript and its corresponding relative position by means of the NCBI reference sequence database (Pruitt ). Sites that were located within an intron (4%) or upstream of the 3′UTR (45%) were removed. Ago HITS-CLIP included 20 miRNAs, whereas 18 of which are broadly conserved [according to (Friedman )]. We proceeded with conserved miRNAs. All mRNA and miRNA data were downloaded from UCSC (Karolchik ) and miRBase (Griffiths-Jones, 2010) on October 2010.

Based on the set of conserved miRNA sequences and mRNA 3′ UTR sequences, we determined all sites complementary to a minimum of six contiguous nucleotides beginning at either position one, two or three relative to the 5′-end of the miRNA. Seed matches were classified functional or non-functional by means of their distance to Ago HITS-CLIP sites. To account for all seed start positions, each seed match located within a distance of two nucleotides to an Ago HITS-CLIP site was tagged functional. Since the Ago HITS-CLIP sites were located almost at the center of an average Ago-mRNA footprint, matches located within a distance of 3–31 nt could also be functional. Since the locations of the footprints were not available, an unambiguous classification was not feasible. To avoid false positives, these sites remained unclassified. All seed matches located beyond the footprint (distance > 31) were classified as non-functional. Further, two miRNAs whose target sites were not significantly enriched in the footprints were removed from the dataset (Supplementary Table S2). Finally, we got 7342 functional, 64 689 non-functional and 1755 unclassified seed sites. Verifying a required minimum target site length of 6 nt, we determined all 5mer matches. The frequency of seed matches within a footprint (distance ≤ 31) and beyond of it was calculated for each seed match length. Additionally, to support our results we prepared the data of the PAR-CLIP experiment in a quite similar fashion (Hafner ) (Supplementary Material).

We defined the background set Ω based on the functional and non-functional sites. A seed match s ∈ Ω was distinguished by its start position relative to the miRNA 5′UTR (1 = α, 2 = β, 3 = γ) and its length. The outcome of this were 20 match types S for a length k and a start position type p. The distributions of all seed match types were disjoint that is each seed match was graded by the longest possible type. To reduce unnecessary complexity of the seed match type set, we merged iteratively non-significant seed match types with their superset. Due to the hierarchical structure of Ω, we were able to apply a separate-and-conquer algorithm (Supplementary Algorithm S1). First we divided the target sites by their seed match start position. Thus, we got three supersets composed of seed matches of a minimum length of 6 nt containing all seed types: S+α,6, S+β,6, S+γ,6. These sets were separated into 6mers having a mismatch at their subsequent position (Sα,6, Sβ,6, Sγ,6) and seed matches having a minimum length of 7 nt S+. We tested the null hypothesis stating that the distribution of functional and non-functional target sites is independent of a mismatch at the 3′ most subsequent position of a seed match. Thus, if the proportions of functional to non-functional target sites between the S and the S+ seed types were not significantly varying (P > 0.05), the separation terminated otherwise the procedure was continued for the next seed type length. A P-value was calculated by means of a two-tailed Fisher's exact test. The α seed site separation terminated after three steps, the β seed matches contained two significant subsets and γ yielded no significant subsets. We termed the found significant seed types based on their start position and their length: S = ‘kmerp’. For standardization, we renamed the endmost subsets: S+α,8 = 8merα, S+β,7 = 7merβ, S+γ,6 = 6merγ.

To estimate the significance of our seed type set, we compared the distribution of the functional sites with a randomized pool of functional seed matches. We drew without replacement a subset of 7803 instances of the multinomial distribution Ω from functional and non-functional seed matches. A P-value was calculated by means of a χ2 test of independence.

The impact of the seed types to miRNA target site prediction was evaluated in terms of recall and specificity. The recall estimates how many of the functional target sites OP are covered by a certain seed type S and the specificity computes the fraction of correctly excluded non-functional target sites ON.

The recall and specificity of each miRNA target prediction algorithm was determined in terms of pure seed finding. Their seed type selection was assigned as described in the related literature. Due to ambiguous seed type assignments based on the first position of the target sequence, the specificity and recall values for TargetScan were computed by executing predictions on our mRNA set.

To estimate the miRNA seed type usage, we calculated the relative frequencies of a seed type for a certain miRNA. These values were normalized by the mean μ and the SD σ:

The conservation of each seed site was determined as described in (Betel ). We used the software package PHAST (Siepel ). The algorithm PhastCons is based on a phylogenetic hidden Markov model, which is fitted to the input sequence by maximum likelihood. Each nucleotide gets a score, which measures the evolutionary conservation across 17 vertebrates. For each seed match, the absolute chromosomal coordinates were determined and a conservation score was calculated. Only if the score of each nucleotide within a functional seed match exceeded the threshold of 0.57 (Betel ), the site was tagged conserved in mammals. The background conservation of a seed type was computed by calculating the fraction of conserved nucleotides of a non-redundant set of UTRs holding a specific seed type. For all statistical computations, the R programming language was applied (R Development Core Team, 2010).

3 RESULTS

3.1 Canonical seed types of miRNA target recognition

In this work, we defined a set of canonical seed types by analyzing the seed matches of experimentally verified functional target sites in the 3′UTR. The Ago HITS-CLIP miRNA:mRNA interaction map (murine assembly of 2006) (Chi ) lists 15 665 chromosomal positions of perfectly matching seed sites of length 6–8 nt belonging to 20 miRNAs frequently bound in Ago complexes. We mapped these sites to annotated protein-coding mRNA transcripts and retained sites located within the 3′UTR. For each miRNA, we scanned the 3′UTRs of the transcript set for all sites complementary to a miRNA subsequence beginning at either position one (α-position), two (β-position) or three (γ-position) relative to the miRNA 5′-end. We required a minimum length of 6 nt. Seed matches of length five [as reported by (Brennecke )] were not significantly enriched in average Ago footprints (Supplementary Table S1). We classified these sites by means of their distance to an Ago HITS-CLIP site and retained miRNAs significantly enriched in footprints. This resulted in 2369 murine genes containing 7070 Ago HITS-CLIP sites of 16 broadly conserved miRNAs.

Each contiguous seed match was defined by its start position type and its length. The dataset was composed of eight α-, seven β- and five γ- seed match types (Supplementary Table S3). Following the law of Occam's razor, the simplest seed type setting for target prediction should usually be the correct one. To reduce unnecessary complexity of the seed type set, we identified unique seed types differing significantly from their superset in terms of functional and non-functional site distribution. For the murine and the human dataset, we achieved six different, disjunct types of seeds: three 6mers either beginning at the first nucleotide (6merα), the second nucleotide (6merβ) or the third nucleotide (6merγ), two 7mers either starting at position one (7merα) or position two (7merβ) and one 8mer beginning at the first nucleotide (8merα) (Supplementary Fig. S1A). These canonical seed types terminated within the first 8 nt of the miRNA in 97% of cases. This underscores the importance of the octamer at the miRNA 5′-end. The significance of this seed type set was evaluated by a sampling approach. The log odds ratio of long seed types is above zero, pointing to a better discrimination between functional and non-functional sites (Table 1).

Table 1.

Determined canonical seed types

	Functional		Non-functional
Seed type	Frequency	%	Frequency	%	LORa	P-value
6merα	1793	24	20 746	32	−0.12	1.20E−028
6merβ	1382	19	13 500	21	−0.04	2.57E−004
6merγ	1755	24	17 954	28	−0.06	2.26E−009
7merα	760	10	5036	8	0.12	2.03E−013
7merβ	959	13	5250	8	0.21	1.34E−042
8merα	693	9	2203	3	0.44	7.60E−132

aLog odds ratio based on sampling.

In a previous work, (Bartel, 2009) defined seeds of miRNA target recognition. Comparing this previous definition with our canonical set of seed types, we recover this set and extend it by additional seed types starting at the α-position (Fig. 1). The miRNA seed match starting at position two and requiring a length of at least 6 nt was described as miRNA core seed (Bartel, 2009; Friedman ; Grimson ). In our canonical set, it is covered by seed types 6merβ, 7merα, 7merβ and 8merα.

Fig. 1.

Definition of seed types. The seed types were termed by the start position relative to the 5′-end of the miRNA and the length of the consecutive seed match. The defined set of canonical seed types can be surjectively projected to the seed type set of (Bartel, 2009). Equivalent definitions could be found for 6merγ, 6merβ and 7merβ. In the case of miRNAs having a seed sequence beginning with an uracile, 7merα complies with 7mer-A1 and 8merα is equal to 8mer. Otherwise 6merβ equates 7mer-A1 and 7mer β complies with 8mer. If the first position within the target sequence is not an adenine, 8merα equates 7mer-m8 and 7merα is equal to 6mer. Additionally, our set considered 6mer matches that are complementary to the first position of a miRNA seed (6merα). Common target site prediction tools focus on seeds of length seven and eight to increase precision.

3.2 Majority of functional sites are based on 6mer seeds

We examined recall and specificity affected by the individual seed types (Fig. 2). Focusing on the relative contribution of each seed type to functional sites, 6mer seeds make up the highest fraction of true target sites (recall: 0.67). On the other hand, 6mer types involve many false positives leading in sum to a low specificity (0.19) and precision (0.09, Supplementary Table S6). In terms of in silico target site classification, the usage of a short seed type causes an inverse prediction [Matthews correlation coefficient (MCC) < 0, Supplementary Table S6], suggesting the avoidance of such a type. In this case, reversing the classification would yield a result superior to an average random prediction.

Fig. 2.

Accuracy evaluation. (A and C) The impact of each seed type on miRNA target site prediction was determined by means of recall and specificity. The effect of the (default) seed type selection is shown for several prediction algorithms. These values present respectively the minimum specificity and the maximum recall of the tools. The dashed line shows an average random prediction. (B and D) Removing non-conserved target sites increases the precision, but lowers the recall. Note that panels (C and D) do not reflect the ranking of predictions based on the algorithms' scoring schemes.

Barely one-third of all genuine target sites are covered by seeds of length 7 and 8. Among these seed types, 7merβ holds the highest recall (0.13) and 8merα shows the best specificity (0.97). The combined set of 7- and 8mer matches achieves a specificity of 0.8 (precision: 0.19). miRNAs perform fine-tuning of gene expression, in particular 6mer seed matches are associated with low repressive effects (Friedman ). As most of the functional sites are formed by short seed sites, one can infer marginal reduction of the mRNA level to be the predominant effect of global miRNA-mediated regulation.

Computing recall and specificity in terms of miRNA:mRNA interactions resulted in a growth of both measurements for each seed type (Supplementary Fig. S2, Table S6 and S7). Here, only the presence of a site on a mRNA matters, whereas in terms of miRNA target site determination the location of a seed match relative to an Ago footprint is important. Multiple matches of one miRNA on a target mRNA are combined into one miRNA:mRNA interaction. Consequently, multiple false positive seed matches may be united to one true positive miRNA:3′UTR interaction. Conversely, multiple true negative target sites may be merged to one false positive interaction. Obviously, recall benefits but specificity suffers from these facts.

Moreover, we wondered if the seed type distributions differ between miRNAs. The relative frequencies of the seed types were computed for each miRNA. A z-score indicates miRNAs holding a frequency over or below the mean frequency given a specific seed type (Fig. 3). Interestingly, 6mer seed types and long seed types are grouped to clusters, respectively, demonstrating that a miRNA either binds to long sites or to short sites but not to both. Further, two main miRNA cluster appeared. The larger group contains miRNAs binding primarily to 6mer-based functional sites. Seven of the 16 miRNAs carry out stronger repression by pairing to rather long seed matches.

Fig. 3.

Heatmap showing the seed type distribution for each miRNA. The colors affected by the row z-score indicate the bias of miRNAs to prefer targets holding a specific seed type. A red/blue coloration implies a higher/lower usage of a seed type compared to other miRNAs.

The importance of short seed types gains further support by the observation that 37% of the 3′UTRs contain exclusively seed matches of length six in their Ago footprints (Supplementary Table S4 lists the numbers of 3′UTRs containing seed matches of exclusively one type.) Interestingly, the sequences of this subset of 3′UTRs are significantly shorter than that of the superset (t-test, P = 4.53E−06). Stark et al. (2005) studied the impact of miRNA regulation on 3′UTR evolution and found that short 3′UTRs indicate avoidance of miRNA regulation. This goes well with our observation of short 3′UTRs regulated by less effective 6mer matches.

3.3 Non-conserved targeting relies on short seeds

We used a strategy established by (Betel ) to identify seed sites conserved across mammals (Fig. 4 and Supplementary Fig. S1B). The majority of functional target sites is conserved (60%). All seed types have a higher fraction of conserved sites than one would expect by chance, given the conservation of their 3′UTRs (Supplementary Table S5). The 6mer sites reveal an almost equal partitioning in conserved and non-conserved sites. A clear discrepancy between the numbers of conserved and non-conserved sites emerges for 7- and 8mer seeds. Particularly, 8merα seed matches exhibit a significant tendency to be conserved. The number of conserved sites in this case is more than three times as high as the number of non-conserved sites. In terms of 7mer seeds, about two-thirds of the seed matches are conserved, whereas 7merα exceeds 7merβ.

Fig. 4.

Observed and expected fraction of conserved seed matches for each seed type illustrated for functional target sites.

In summary, the mean probability to be conserved is about 55% for a 6mer seed. In contrast, 7mer and 8mer seeds have a probability up to 77% to be conserved. Further, a total of 75% of the functional non-conserved sites are covered by 6mer seeds. Therefore, non-conserved or species-specific targeting relies to a large extent on target sites containing short seeds.

Keeping only the conserved sites from the set of seed matches lifts specificity of all seed types (Fig. 2, Supplementary Table S9). In particular, the 6mer seeds show a significant increase of specificity leading to a classification better than an average random prediction (MCC > 0, Supplementary Table S6).

3.4 Target prediction focuses on 7- and 8mer seed matches

We reviewed frequently used approaches for target prediction in mammals with regard to the implemented seed types (Table 2). The TargetScan algorithm (Grimson ) seeks mainly for seeds of length seven and eight via seed types 7mer-A1, 7mer-m8 and 8mer. The 7mer-A1 sites may be of type 6merβ in the event the miRNA sequence starts with a nucleotide different to uracile. However, the majority of mammalian miRNAs begins with an U (Lewis ). Both PicTar (Krek ) and EIMMO (Gaidatzis ) require stringent seed pairing between 7 nt starting at either the α or the β-position. A novel approach called TargetSpy seed (Sturm ) restricts the set of seed matches to predictions containing a perfect 7mer.

Table 2.

Default miRNA seed type selection of prediction algorithms

	Seed type
Algorithm	6merα	6merβ	6merγ	7merα	7merβ	8merα
PITAa		✓		✓	✓	✓
TargetScan		✓b		✓c	✓	✓
PicTar				✓	✓	✓
EIMMO				✓	✓	✓
TargetSpy S.				✓	✓	✓
PACMITa					✓	✓

aConfigurable seed length, default seed types ensure high precision.

bIf miRNA seed sequence starts with an adenine, guanine, cytosine.

cIf miRNA seed sequence starts with an uracile.

Some algorithms allow for custom-defined seed searching: PITA (Kertesz ) seeks by default for sites of length six, seven and eight that start at position two of the miRNA. The standard setting of PACMIT (Marn and Vaníček, 2011) is even more restrictive by considering merely sites matching to miRNA positions two to eight. Both tools enable adjusting of the site length by the user. RNAhybrid (Krüger and Rehmsmeier, 2006) as well as IntaRNA (Busch ) are more flexible by providing a couple of additional parameters to customize the seed search, e.g. setting the start position. The latter is a general approach to predict RNA:RNA interactions. Both do not suggest default seed search parameters.

The impact of the (default) seed type selection of prediction algorithms on recall and specificity was evaluated (Fig. 2, Supplementary Tables S8 and S9).

Prediction methods implement scoring schemes to value target site characteristics beside the seed. In contrast to common evaluations of miRNA target site prediction algorithms, this is not an assessment of a subset of top scored instances but of all predictions. Therefore, the denoted specificity values represent the minima while the recall values show the maxima for the (default) seed choice, respectively. Subsets composed of top scored predictions would achieve significantly higher specificity values.

Obviously, all prediction models exhibit a considerable constraint regarding their ability of finding potential target sites. PITA holds the highest recall of 52% (specificity: 60%) owing to the exhaustive search for 6merβ seed matches, whereas PACMIT has the lowest recall of 23% (specificity: 88%) restricted to find less than a quarter of all functional seed sites. Additional filtering by removing conserved sites increases the specificity but consequently lowers the recall. Here, PACMIT could only find 16% of all functional sites (specificity: 73%). A higher recall but a lower specificity can be observed for the prediction of miRNA:mRNA interactions (Supplementary Fig. S2). Concluding, due to the significant gain of precision, tool developers recommend to use long seeds. Our study quantified the loss of recall accompanied by this proceeding.

4 CONCLUSION

In this study, we present an analysis of the most important feature for miRNA target recognition, the so-called miRNA seed, using a large-scale dataset of functional target sites. Based on the Ago HITS-CLIP and PAR-CLIP miRNA:mRNA interaction maps, we analyzed seeds properties and their influences on miRNA target site prediction methods. Due to the definite specification of Ago binding sites, we were able to classify miRNA recognition elements contained in the mRNA 3′UTR as either functional or non-functional. We defined a minimal set of seed types that is sufficient for accurate miRNA target site predictions. The final data pool allows for enhanced analysis of miRNA target prediction algorithms compared to earlier studies that were restricted by experimental constraints (Alexiou ; Selbach ). We found that most conserved miRNAs interact predominantly with target sites endowed with short seed matches; 67% of functional sites are based on 6mer seeds. In contrast, common prediction algorithms focus mainly on seeds of length seven or eight. At present, prediction algorithms have to accept severe deficiencies of recall to ensure high specificity that is naturally considered to be more important. Moreover, the preferential search for long seeds lifts the proportion of conserved sites. But we found that a substantial fraction (40%) of all functional target sites is not conserved. Target sites including 6mer seeds are enriched among these.

Concluding, the problem of recall can be easily translated to the problem of precision. However, this strongly intensifies the need for features beyond seed pairing that realistically describe miRNA targeting, in particular non-conserved target sites. It may also raise the basic question for the potential of seed-based approaches in discriminating between functional and non-functional sites.

Funding: Helmholtz Association within the Helmholtz Alliance on Systems Biology (project CoReNe); Federal Ministry of Education and Research (BMBF) within project SysMBo (Förderkennzeichen: 0315494A).

Conflict of Interest: none declared.