Reducing ligation bias of small RNAs in libraries for next generation sequencing

Abstract

Background: The use of nucleic acid-modifying enzymes has driven the rapid advancement in molecular biology. Understanding their function is important for modifying or improving their activity. However, functional analysis usually relies upon low-throughput experiments. Here we present a method for functional analysis of nucleic acid-modifying enzymes using next generation sequencing.

Findings: We demonstrate that sequencing data of libraries generated by RNA ligases can reveal novel secondary structure preferences of these enzymes, which are used in small RNA cloning and library preparation for NGS. Using this knowledge we demonstrate that the cloning bias in small RNA libraries is RNA ligase-dependent. We developed a high definition (HD) protocol that reduces the RNA ligase-dependent cloning bias. The HD protocol doubled read coverage, is quantitative and found previously unidentified microRNAs. In addition, we show that microRNAs in miRBase are those preferred by the adapters of the main sequencing platform.

Conclusions: Sequencing bias of small RNAs partially influenced which microRNAs have been studied in depth; therefore most previous small RNA profiling experiments should be re-evaluated. New microRNAs are likely to be found, which were selected against by existing adapters. Preference of currently used adapters towards known microRNAs suggests that the annotation of all existing small RNAs, including miRNAs, siRNAs and piRNAs, has been biased.

Figure 1

Scheme depicting the experimental approach and HD adapters. a Data were generated to analyse the sequence preferences of T4 Rnl1 and T4 Rnl2 using a degenerate RNA library (N21 RNA). b HD adapters include degenerate tags at the end of the adapters that allow the formation of stable secondary structures for more sequences and reduce RNA ligase-dependent sequence bias. Panel (c) shows the structure of miR-29b with the Illumina adapters (top) and some of the structures formed by HD adapters (bottom). We found 1,031 distinct structures originating from 12,479 tag combinations.

Figure 2

Sequencing cDNA generated from N21 RNA libraries. a Number of reads for the 100 most abundant sequences in the N21 libraries, prepared with Illumina (red) or HD adapters (blue). b-d Frequencies of predicted nucleotide base-pairing per position for N21 insert (b), N21 insert and 3’ adapter (c) and 5’ adapter, insert and 3’ adapter (d). In (c) and (d) vertical dotted line indicates ligation point. Red line denotes data obtained with Illumina protocol, blue line with HD protocol and grey line randomly generated sets of 21nt sequences. Bars indicate minimum and maximum values in all replicates. Horizontal bars at bottom indicate sequence region: green, insert; red, 3’ adapter; blue, 5’ adapter. For insert folding frequencies obtained with random sequences are more closely matched by HD data (R² = 0.83) than by Illumina data (R² = 0.60). e Comparison of T4 Rnl2 ligase activity on substrates with ss flaps of differing nucleotide lengths upstream or downstream of ligation site. In vitro ligation assay of RNA-DNA duplexes with either a nick (0NT) or ss flaps up- or downstream from the ligation site was carried out at 25°C for 30 min. Substrates with ss flaps >2nt in length upstream from the ligation site are inefficiently ligated. The diagram illustrates the position of the flaps, the fluorescein reporter group (star) and the backbone oligonucleotide (black). If ligation occurs the size of the nucleic acid attached to the fluorescein increases as visualised by 15% PAGE.

Figure 3

cDNA library preparation protocols distort miRNA research. a Comparison of change in miRNA level between wild-type and Dicer KO DLD cells obtained in Illumina (x axis) and HD samples (y axis). R² = 0.62. b Number of known miRNAs found in DLD cells at different thresholds using Illumina or HD adapters. Regardless of chosen threshold, HD adapters identify more miRNAs. c Absolute quantification of eight known miRNAs (let-7i, miR-10a, miR-19b, miR-21, miR-25, miR-29b, miR-93, miR-375) obtained by Northern blot compared with number of times these miRNAs were sequenced using Illumina or HD adapters in DLD cell line. Data obtained with HD adapters correlates better with absolute quantifications (R² = 0.70) than Illumina data (R² = 0.12). d Number of PubMed citations and number of reads per experiment (data obtained from miRbase v17) of miRNAs conserved between mouse and human. MiRNAs with higher number of reads tend to be more extensively studied [R² = 0.58, p-value < 10⁽⁻¹⁵⁾]. e-f Distributions of minimum free energy (MFE) of known human miRNAs concatenated with 5’ and 3’ adapter sequences. Using Illumina adapter sequences sRNA cloning kit V1.5 the set of miRNAs found by Illumina has lower average MFE than the set of miRNAs found by 454 (Wilcoxon test p = 0.01). We found the same result using the 3' adapter from sRNA cloning kit V1.0 (data not shown). e Conversely, using 454 adapter sequences average MFE is lower for set of miRNAs found by 454 (p = 0.07). f Analogous results for concatenation of miRNA only with 3’ adapter display a similar trend (see Additional file 5: Figure S7).