Flux Simulator for RNA-seq Simulation
High-throughput sequencing of cDNA libraries constructed from cellular RNA complements (RNA-Seq) naturally provides a digital quantitative measurement for every expressed RNA molecule. Nature, impact and mutual interference of biases in different experimental setups are, however, still poorly understood- mostly due to the lack of data from intermediate protocol steps. We analysed multiple RNA-Seq experiments, involving different sample preparation protocols and sequencing platforms: we broke them down into their common–and currently indispensable–technical components (reverse transcription, fragmentation, adapter ligation, PCR amplification, gel segregation and sequencing), investigating how such different steps influence abundance and distribution of the sequenced reads. For each of those steps, we developed universally applicable models, which can be parameterised by empirical attributes of any experimental protocol. Our models are implemented in a computer simulation pipeline called the Flux Simulator, and we show that read distributions generated by different combinations of these models reproduce well corresponding evidence obtained from the corresponding experimental setups. We further demonstrate that our in silico RNA-Seq provides insights about hidden precursors that determine the final configuration of reads along gene bodies; enhancing or compensatory effects that explain apparently controversial observations can be observed. Moreover, our simulations identify hitherto unreported sources of systematic bias from RNA hydrolysis, a fragmentation technique currently employed by most RNA-Seq protocols.
From the discussion section:
If we accept that our bioinformatics models capture the main origins of the experimental biases, our simulation enables us to investigate intermediate stages of sequencing protocolsusually hidden layers of RNA-Seq (Figures 24). Specifically, we give computational and experimental evidence as to why insert size distributions obtained by hydrolysis differ substantially between transcripts of different lengths (Figure 2). In the light of the uniform random fragmentation model we developed, the dependence of the RNA molecules geometry on their length can be interpreted as shorter molecules being more linearised when hydrolysed, whereas longer RNA polymersin spite of strongly denaturing conditionsstill tend to form higher order structures. Therefore, size filtering alters the way transcripts are represented in the library as a function of the length of the original RNA molecule.
In addition, our models show why fragments obtained by sequence-independent fragmentation processes, as for instance cDNA nebulisation or RNA hydrolysis, are not uniformly distributed along the fragmented molecule, but occur more frequently at rather specific points: the ends of nebulised fragments accumulate at the midpoints of recursively split molecules (Figure 6D), whereas fragment density obtained by RNA hydrolysis propagates from a transcripts ends towards its centre in patterns produced by characteristic Weibull distribution of the obtained insert sizes (Figure 3). Onto these patterns one has to superimpose sequence-specific biases (Figures 4 and ?and5).5). If heterogeneous transcripts are considered together, however, the recognition of these biases on large scale is complicated (Figure 6).
As for the computational analysis of RNA-Seq experiments, we consider our simulation-based studies as a serious motivation to debunk the widespread belief that all biases should affect the interpretation of data negatively: in fact, well-understood biases of systematic nature are valuable as additional sources of information. Therefore, we are convinced that the critical evaluation of experiments mimicked in silico will have an increasing impact on design and evaluation of bioinformatics approaches to RNA-Seq.
