
==== Front
Brief Bioinform
Brief Bioinform
bib
Briefings in Bioinformatics
1467-5463
1477-4054
Oxford University Press

10.1093/bib/bbae164
bbae164
Problem Solving Protocol
AcademicSubjects/SCI01060
BEERS2: RNA-Seq simulation through high fidelity in silico modeling
https://orcid.org/0000-0002-6980-0079
Brooks Thomas G Institute for Translational Medicine and Therapeutics, University of Pennsylvania, PA, USA

Lahens Nicholas F Institute for Translational Medicine and Therapeutics, University of Pennsylvania, PA, USA

Mrčela Antonijo Institute for Translational Medicine and Therapeutics, University of Pennsylvania, PA, USA

Sarantopoulou Dimitra Institute for Translational Medicine and Therapeutics, University of Pennsylvania, PA, USA
Current address: National Institute on Aging, National Institutes of Health, Baltimore, MD, USA

Nayak Soumyashant Institute for Translational Medicine and Therapeutics, University of Pennsylvania, PA, USA
Current address: Statistics and Mathematics Unit, Indian Statistical Institute, Bengaluru, Karnataka, India

Naik Amruta Institute for Translational Medicine and Therapeutics, University of Pennsylvania, PA, USA
Children’s Hospital of Philadelphia, Philadelphia, PA, USA

Sengupta Shaon Institute for Translational Medicine and Therapeutics, University of Pennsylvania, PA, USA
Children’s Hospital of Philadelphia, Philadelphia, PA, USA
Department of Pediatrics, University of Pennsylvania Perelman School of Medicine, Philadelphia, Pennsylvania, USA

Choi Peter S Division of Cancer Pathobiology, Children’s Hospital of Philadelphia, Philadelphia, PA, USA
Department of Pathology & Laboratory Medicine, University of Pennsylvania Perelman School of Medicine, Philadelphia, PA, USA

Grant Gregory R Institute for Translational Medicine and Therapeutics, University of Pennsylvania, PA, USA
Department of Genetics, University of Pennsylvania, Philadelphia, PA, USA

Corresponding author: Thomas G. Brooks1, Institute for Translational Medicine and Therapeutics, University of Pennsylvania, PA, USA. Tel: +1 215-573-3736; Fax:+1 2156626300; E-mail: thobr@pennmedicine.upenn.edu
5 2024
11 4 2024
11 4 2024
25 3 bbae16412 7 2023
26 1 2024
26 3 2024
© The Author(s) 2024. Published by Oxford University Press.
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com

Abstract

Simulation of RNA-seq reads is critical in the assessment, comparison, benchmarking and development of bioinformatics tools. Yet the field of RNA-seq simulators has progressed little in the last decade. To address this need we have developed BEERS2, which combines a flexible and highly configurable design with detailed simulation of the entire library preparation and sequencing pipeline. BEERS2 takes input transcripts (typically fully length messenger RNA transcripts with polyA tails) from either customizable input or from CAMPAREE simulated RNA samples. It produces realistic reads of these transcripts as FASTQ, SAM or BAM formats with the SAM or BAM formats containing the true alignment to the reference genome. It also produces true transcript-level quantification values. BEERS2 combines a flexible and highly configurable design with detailed simulation of the entire library preparation and sequencing pipeline and is designed to include the effects of polyA selection and RiboZero for ribosomal depletion, hexamer priming sequence biases, GC-content biases in polymerase chain reaction (PCR) amplification, barcode read errors and errors during PCR amplification. These characteristics combine to make BEERS2 the most complete simulation of RNA-seq to date. Finally, we demonstrate the use of BEERS2 by measuring the effect of several settings on the popular Salmon pseudoalignment algorithm.

simulations
in silico
RNA-seq
short read sequencing
benchmark
pseudoalignment
Next Generation of RNA-Seq Simulators for Benchmarking Analyses R21-LM012763-01A1 National Center for Advancing Translational Sciences 10.13039/100006108 5UL1TR000003 NHLBI R01HL155934-01A1(SS) NHLBI-R01HL147472 DP2GM146251
==== Body
pmcINTRODUCTION

Since the introduction of RNA-Seq circa 2010, there has been a tremendous proliferation of bioinformatics tools published for every stage of data processing and analysis, for example over a dozen differential expression methods exist [1]. Such a landscape of tools naturally calls for unbiased benchmarking studies to compare and evaluate their performances and to develop best practice guidelines. Benchmarking studies are most informative when there are data for which the ground truth is known. In the case of real world omics data, the ground truth is impossible to determine without the use of the same tools being evaluated, and real-world data are inherently limited for evaluating metrics [2].

Simulation is a popular solution; however, the further one goes downstream in the analysis pipeline of RNA-Seq, the more realistic the simulation must be. Benchmarking genome alignment with simulated data only requires simulating realistic reads because genome aligners only work on one read at a time. However, benchmarking the downstream tasks of quantification and normalization requires simulating entire realistic samples of reads because these steps combine the information across reads, and modern methods like pseudoalignment [3, 4] can combine traditional alignment and quantification into a single step. Therefore, there is a strong need for sophisticated RNA-Seq simulation, and current simulators are too simplistic [5].

Developing read-level simulated data from scratch is labor-intensive and depending on the specific application may not be necessary. Indeed, much of the benchmarking of RNA-seq tasks in the literature [6–9] involves directly simulating the final normalized spreadsheet of gene, or transcript, level quantifications without ever simulating reads. However, this spreadsheet level simulation cannot fully capture the difficulties of handling real reads [10] and read-level simulation is necessary to cover the full spectrum of RNA-Seq analysis methods. Both alignment and quantification, for example, require read-level data.

When embarking on a benchmarking study with simulated read-level data, it is necessary to budget a fair amount of time, effort and thought into generating the simulated data. Otherwise, the results are only as meaningful as the data are realistic and it is very easy to oversimplify when simulating biology. In this work we have provided powerful tools to achieve reads and samples of reads that are as realistic as possible. But these are not push-button applications, they are aids to science that must be employed with expert knowledge.

We built the original BEERS RNA-Seq read-level simulator in 2011 [11] specifically to benchmark alignment. However, in the years that followed, we realized that BEERS was being used to benchmark all stages of the processing and analysis pipeline, which was far beyond its intended purpose [12–14]. To this end, we present BEERS2, in which we have endeavored to model every step of the library prep pipeline in silico.

In addition to earlier versions of BEERS, there are existing RNA-seq read-level simulators. Polyester [15] supports configurable fragment GC bias and differential expression. The ASimulatoR [16] project adds alternative-splice configuration to Polyester. The RSEM simulator [17] allows for modeling position bias such as from polyA selection. The Flux Simulator [18] can simulate GC bias in polymerase chain reaction (PCR) amplification, hexamer priming sequence bias and a 3′ bias due to poly-dT primers but not due to polyA selection. SimBA [19] adds genomic mutations onto the Flux Simulator. CuReSim [20] and RandomReads [21] generate reads with configurable errors but no simulated biases. PBSIM [22] can simulate long reads from Pacific Biosciences instruments. The rlsim [23] package allows simulating with hexamer priming bias and PCR amplification with GC and length biases. However, none of these can model all four of the PolyA selection, RiboZero, GC and primer biases. In addition, they lack some more advanced features: errors in PCR that propagate realistically to multiple descendant molecules, read errors in molecule barcodes leading to demultiplexing difficulties and diploid genomes with allele-specific effects.

Simulating RNA-Seq data for general purpose involves first simulating ‘expression’ and then simulating ‘sequencing’. Given the complexity of both of these steps, we have separated them into two different applications, reflecting the separation of ‘expression’ and ‘measurements’ models of RNA-seq [24]. The expression simulator is called CAMPAREE [25] and it generates full length RNA molecules from a diploid genome inferred separately for each sample. CAMPAREE takes real RNA-Seq data as input, from which it infers variants and expression levels and generates diploid genomes. These are then used to simulate a sample of RNA-seq molecules with a known set of underlying properties.

BEERS2 takes the set of ‘molecules’ output by CAMPAREE as input, but it can just as easily take as input any properly formatted set of molecules, which are, typically, full length RNA, but can be RNA molecules of any sort, see Figure 1A. BEERS2 transforms these input molecules into RNA-Seq reads, by mimicking one of the various library prep protocols followed by sequencing. Briefly, this involves ribosomal depletion, fragmentation, various PCR steps, labeling and then sequencing by synthesis, which involves flowcell hybridization and bridge amplification, and base assessment, see Figure 1B. Errors and other features of sequencing are incurred at the various steps that they occur during library prep. The system is modular, so that modeling new library prep protocols or implementing changes to library preps requires only local changes to the modules for the steps that have been affected. For example, polyA selection can be easily swapped out for a different ribosomal depletion protocol without changing other steps in library prep simulation. Further, this modular nature also allows for isolating specific steps to investigate the effects of their specific biases on the analysis under question.

Figure 1 BEERS2 Overview. (A) BEERS2 is a customizable pipeline which processes input transcripts (typically fully length mRNA transcripts with polyA tails) from either customizable input or from CAMPAREE simulated RNA samples. It produces realistic reads of these transcripts as FASTQ, SAM or BAM formats with the SAM or BAM formats containing the true alignment to the reference genome. It also produces true transcript-level quantification values. (B) The default pipeline simulates the Illumina Stranded mRNA library preparation process, followed by sequencing. Each step influences one or more aspects of the output sample through configurable parameters.

METHODS

BEERS2 is developed in the Python v3.11 language and uses the Snakemake [26] pipeline library v7.24.0 for parallelization and reproducibility. Since BEERS2 simulates every individual molecule through the entire library preparation and sequencing process, realistic BEERS2 executions require significant compute resources. Our Snakemake pipeline can be configured to work on many common cluster compute environments, but to further ease the process of installing the necessary software and configuring it to work on a particular cluster, we prepared a stack template which can be easily deployed in any AWS account. The deployed stack will contain both the software and hardware necessary to execute the BEERS2 pipeline, and researchers can start the execution immediately after providing the required input for BEERS2 (Figure S1). We estimate a $10–15 cost on AWS per sample, primarily for compute, and a duration of approximately 4.5 h. Runtime costs and durations depend substantially on configuration parameters, but the runtime is approximately constant for running multiple samples as these are done in parallel.

The modular design allows mixing and matching library pipeline steps through the configuration file. For example, a DNA-seq pipeline could be simulated by reusing the shared aspects of the RNA-seq pipeline (such as fragmentation, size selection and PCR amplification) and dropping the distinct steps (such as polyA-tail selection). This allows BEERS2 to keep up with ever-changing protocols.

BEERS2 includes modules simulating the individual steps of the two standard RNA library preps, Illumina Stranded messenger RNA (mRNA) and Illumina Stranded Total RNA, followed by simulated sequencing by an Illumina sequencing machine. The individual steps are modeled as closely as possible to the real process. First, input molecules are selected for having a polyA tail, where the selection criteria consist of configurable acceptance probabilities as a function of tail length. This step also introduces an optional 3′ bias by modeling truncated molecules by breaking bonds between bases with uniform probability. Alternatively, the RiboZero selection module computes sequence-similarity to a population of rRNA reference oligos and degrades molecules at sites with probability dependent upon the sequence-similarity.

Next, the fragmentation step breaks molecules by default uniformly, meaning that every bond between base pairs is equally likely to break, or with a configurable positional and fragment length bias, see Supplemental Methods. Afterwards, the hexamer priming and complementary DNA (cDNA) synthesis is simulated. Hexamers are modeled as first selecting a random number of primers to bind (according to a binomial distribution with n equal to the number of potential binding sites), and then all possible binding locations are weighted according to a configurable hexamer bias model to simulate the known hexamer sequence non-uniformity [27]. Finally, the specified number of binding hexamers are then placed on random locations according to the weights and the 5′ most hexamer is used to synthesize the cDNA. Therefore, at synthesis of both cDNA strands, some fraction of end sequence may be lost if no hexamer binds at the very end of the fragment.

After cDNA is generated, the molecules are size selected with a configurable window retention probability as a function of molecule length. Then, adapters are ligated to either end. Finally, PCR amplification is modeled. Since most resulting molecules do not end up being sequenced, resulting from the process of diluting and loading samples onto an Illumina flow cell, we model the subsequent down sampling at the same time as PCR amplification. This enables efficient modeling with many rounds of PCR amplification without computing molecules that do not end up being retained. Appropriate retention rates for this step can be chosen from the known PCR cycle count with the desired PCR dupe rate, see Supplemental Methods.

Moreover, PCR amplification can have configurable copy errors (insertions, deletions and single base errors) during each round. Errors in early rounds may therefore propagate to multiple sequenced PCR duplicates. PCR amplification also allows a configurable GC bias, with PCR amplification failing with a rate depending upon the overall GC content of the molecule.

After PCR amplification with down-sampling, molecules are run through a simulated sequencer. Molecules are bridge amplified, including single base copy errors in the amplification. This process is approximated by storing the total counts of each of A/C/G/T for each position in the molecule. Then, sequencing by synthesis is simulated by estimating fluorescence from each base pair, including phasing, followed by base-calling and quality scores derived from the process used by Illumina’s Bustard algorithm as described in [28].

Mimicking real samples

BEERS2 is designed to be parameterizable with each step having configurable effects on the output. Any simulated dataset output by BEERS2 is the result of a large number of steps, meaning that the simulation parameters do not specify what the output dataset looks like but rather they specify what happens to the input dataset. This allows BEERS2 to flexibly generate a wide range of scenarios including impossible scenarios like an idealized ‘unbiased’ sample if perfect sequencing were available. However, it means that it takes significant effort if the goal is to mimic a specific real sample. To do this, we recommend running BEERS2 on a selection of parameter values and optimizing it by comparison with a real sample. To make this feasible, we additionally recommend restricting this optimization procedure to a small input sample.

As a demonstration, we selected a single real sample with polyA-selected RNA-seq data. We chose several criteria on which we wanted to mimic the real sample. First, the non-uniformity of coverage data, measured by the coefficient of variation of coverage plots on short (under 3 kb) transcripts. Secondly, the 3′ bias, measured by the slope of the regression line of the log of coverage depth versus the position within the transcript on long (over 3 kb) transcripts. Thirdly, the GC content of reads. And lastly, the read fragment length distributions. To measure these, we selected the 100 most highly expressed genes that have a single annotated transcript, so that position within the gene is unambiguous. This gives a distribution of values, one for each transcript or (for fragment size and GC content) each read. Although BEERS2’s GC content bias uses the GC content of the entire fragment and not just the two paired-end reads, for this process we measure the GC content of the reads alone as it is more straightforward to compare this with the GC content of a real sample, which only contains the reads.

We compared the distributions of these values in the real data with their equivalents in the BEERS2-generated samples according to the earth mover’s distance (EMD), normalized by the standard deviation of the probability distributions. Then we selected a variety of relevant BEERS2 parameters and ran BEERS2 50 times on samples of just those 100 genes with varying settings of those parameters. This was then repeated four more times for a total of 250 BEERS2 runs, with each batch of 50 using the top-performing parameters of the previous batch as initial values before perturbing them to obtain new parameters.

This optimization reduced the total EMD from 2.01 (using the initial, hand-selected values) to 0.94 in the best-performing configuration, see Figure S2 and Table S1. Despite that improvement, we found only moderate alignment of BEERS2 to the real dataset. This highlighted that BEERS2, like other models of RNA-seq [29], has yet to capture the entire complexity of the data. One specific area for future improvement that this identified is that CAMPAREE generates transcripts that match the annotation, but many of these genes in the real dataset showed abrupt drops in coverage long before the annotated end of the 3′ untranslated region. This means that our measure of 3′ bias was often negative in the real dataset but nearly always positive in BEERS2-generated data. To minimize this, we trimmed the last 1 kb off of all genes before computing the regression, but still observed some level of negative 3′ biases.

This demonstrates that while BEERS2 is highly flexible, situations where exactly mimicking a specific real dataset is required could call for alternative methods such as non-parametric methods [30] that more directly reuse real samples.

RESULTS

BEERS2 generates compelling RNA-seq reads under varying sequencing conditions, see Figure 2. We provide two example applications of BEERS2 to demonstrate its use in providing actionable recommendations in RNA-seq analysis through benchmarking. First, we focus on evaluating the performance of a single pseudoaligner called Salmon [4] in order to provide quick and simple applications. To properly perform a full-scale benchmarking comparison of all tools requires a full-length article. Salmon quantifies expression at the transcript level using an RNA-Seq pseudo-alignment approach. Our previous benchmarking study found that Salmon was a top transcript-level quantifier [31], which is the motivation to evaluate its performance further.

Figure 2 Coverage Tracks of Real and Simulated Data. UCSC Genome Browser tracks of the coverage on the gene Lrrc58 from, on top, three real mouse lung RNA-seq datasets with polyA selection, RiboZero selection and no selection step. Below, six simulated datasets under varying configurations: ideal data (base), hexamer priming sequence bias (seq_bias), GC bias in PCR amplification (GC_bias), RiboZero sequence-specific bias (RiboZero), 3′ bias from polyA selection (3prime_bias) and all three biases of a typical polyA run (GC_3prime_seq_bias). The last dataset represents the most realistic scenario. Deviations between BEERS2 coverage and real sample coverage indicate that the full causes of coverage non-uniformity are still incompletely modeled and understood. Lrrc58 gene model included at the top of the figure.

Application 1: Benchmark of bias correction in Salmon

Salmon provides options to correct for three known bias factors in RNA-seq data: position of the fragment relative to the start and end of the transcript (Pos); base sequence near the start and end of the fragment (Seq) and GC content of the fragment (GC). Biases related to these factors can be simulated in BEERS2 by including Poly-A selection with 3′ bias (for Pos bias), primer sequence bias (for Seq bias) and PCR amplification GC bias (for GC bias). This provides an ideal test-case for BEERS2’s bias support by quantifying how Salmon’s correction options affect its accuracy in quantifying transcripts and genes under varying levels of bias. Pseudoalignment-based methods such as Salmon require read-level RNA-Seq benchmarking data such as provided by BEERS2.

CAMPAREE [25] was used to simulate eight mouse liver samples based on eight real samples from a previous study [32]. Simulated sequencing of these samples was performed with BEERS2 under seven conditions with varying levels of bias:

With no bias factors,

With medium or high 3′ bias from polyA selection,

With medium or high GC bias from PCR amplification,

With primer sequence bias or

With high levels of all three biases (All Bias).

Salmon was then run on the FASTQ files generated by BEERS2 for each of the eight samples under each of the eight combinations of the possible bias correction settings (with/without Pos correction, with/without Seq correction, with/without GC correction), see Figure S3. Salmon-quantified TPM values were then compared with the true TPM values and Salmon read counts were compared with the true read counts. Concordance was assessed by the mean absolute relative difference (MARD) of these values, which is defined as the mean (across all transcripts) of ARD = |true – inferred|/|true + inferred|, or in other words the difference of quantified and true values divided by their sum. Lower values of MARD indicate higher mean concordance of quantified and true values. The denominator also acts as a normalization factor putting low and high expressed genes on a comparable scale.

First, we consider the TPM values, see Figure 3A. When no biases are simulated, read depth is still not uniform primarily due to edge effects of fragmentation and size selection, see Figure S3A. Not surprisingly, applying the Pos correction improves accuracy only slightly. Applying either GC or Seq correction shows no improvement, and in fact slightly impaired accuracy unless Pos correction is also applied simultaneously.

Figure 3 Salmon Transcript-Level Accuracy with Bias Correction. We generated BEERS2 samples with varying levels of 3′ bias from PolyA-selection, GC bias from PCR amplification and start and end sequence bias from primer bias for each of eight samples. Then we ran Salmon with each combination of bias correction options for fragment start/end sequence (Seq), position within transcript (Pos) and fragment GC content (GC). MARD computes the mean value of the difference of true and estimated values over their sum. Lower MARD indicates more accurate estimates. MARD values vary between samples considerably due to the relative difficulty of quantification under different bias conditions and read counts. (A) MARD of TPM values. (B) MARD of read count values. The boxes show quartiles.

As expected, when GC bias alone was present, GC correction helped performance and with 3′ bias alone, positional correction resulted in substantial improvement. Meanwhile with primer bias alone, the Seq correction improved performance but primarily only when Pos correction was also enabled.

When all three biases were included, GC and Pos correction showed substantial improvements in quantification, while Seq impaired it slightly.

Next, we considered quantified read count values instead of TPM values, see Figure 3B. Effects were largely similar but smaller than in TPM quantification with almost no changes in read count quantification on the run with all biases enabled. Read counts are overall more accurate than TPM and so may not have as much room to improve. Surprisingly, Pos correction slightly impaired instead of improved read count values when no bias or only primer bias were included, in opposite direction to the TPM values.

Comparing TPM and count values with truth by Spearman correlation instead of MARD yielded a similar picture, see Figure S4, except that Pos correction impaired, rather than improved, quantifications when no 3′ bias was present.

Since transcript-level counts were mostly unaffected by bias correction, we conclude that Salmon’s correction factors are largely affecting TPM by influencing the estimated effective length of the transcript rather than by influencing where reads counts are assigned in the expectation maximization step. One line of future developments therefore could be the use of bias corrections to help disambiguate reads between isoforms. For example, when there are two isoforms with different 3′ transcription stop sites, then accounting for the 3′ bias could help to disambiguate reads that occur far from the 3′ end, thereby improving transcript read count estimates.

Lastly, we assessed gene-level TPM values, see Figure 4. These show modest but consistent improvement from inclusion of all three salmon bias corrections. The only exception was for the medium GC bias sample, which showed a small performance hit when measuring by MARD but not Spearman correlation. We have therefore confirmed that improved transcript effective length estimates also improve the gene-level TPM quantification. We did not assess performance of gene-level count data since this is recommended to be done with additional tools, such as tximport [33], and is best evaluated in the context of differential expression. Furthermore, differential expression analysis can be complicated by differential splicing when isoforms of one gene have different effective lengths, which impacts both gene and transcript level quantified counts.

Figure 4 Salmon Gene-Level Accuracy with Bias Correction. Gene-level TPM values were obtained by summing transcript-level TPMs. Comparison of true gene-level TPM values with Salmon-quantified gene-level TPM values by (A) MARD and (B) Spearman rank correlation (corr). Lower MARD and higher correlations are indicative of better performance. Here, only Salmon results with no bias corrections are compared with Salmon results with all bias corrections (for GC content, positional and sequence biases) applied. In all cases (except MARD in GC_bias_med), applying bias correction yields an increase in accuracy.

Application 2: Inclusion of pre-mRNA in Salmon index

Salmon indices are generated using a FASTA file containing the target transcripts and their sequences which should be quantified. These can be prepared by the user as they wish, but one common source is to use Ensembl’s cDNA files, which do not include pre-mRNA. The official Salmon Getting Started guide uses this file for Arabidopsis thaliana as an example. Since real RNA-seq samples invariably include intronic expression, we investigated whether Salmon’s performance would improve when including pre-mRNA in the index files. One application of pre-mRNA quantification is in single-cell RNA-seq, particularly for RNA velocity, which has been previously assessed [30, 34, 35], while we focus instead on the impacts on mature mRNA quantification in bulk RNA-seq.

CAMPAREE generates pre-mRNA signal at realistic levels, so we used the same BEERS2-generated data from Application 1 above to assess the accuracy of Salmon with and without pre-mRNA in the index. For each gene, we added one pre-mRNA isoform to the Salmon index which contained the full genomic stretch from the first to the last base of all isoforms of the gene.

We compared accuracy by comparing true and estimated CPM values. Inclusion of pre-mRNA in the index slightly improved Salmon’s performance in quantifying mature mRNA among genes whose expression was at least 10% pre-mRNA (MARD 0.24 without pre-mRNA in the index versus 0.22 with; Spearman correlation of 0.935 without pre-mRNA and 0.938 with). In the genes with less than 10% pre-mRNA, including pre-mRNA in the index resulted in slightly worse performance (MARD of 0.256 without and 0.276 with pre-mRNA; Spearman correlation of 0.904 without and 0.892 with pre-mRNA).

Lastly, running Salmon with bootstrapped outputs to estimate technical variance, we noted that the relative quantification range (difference of max to min bootstrap values, divided by the mean bootstrap value) correlated modestly with ARD (Spearman correlation of 0.31 in all highly expressed transcripts and 0.39 in genes with at least 10% pre-mRNA expression). However, genes whose quantification regressed after inclusion of pre-mRNA often showed low technical variance as estimated by the Salmon bootstrap despite having large errors compared with the true simulated counts. This indicates that the difficulties in quantification are poorly captured by the bootstrap process.

We then reran all these comparisons on BEERS2 data with all biases enabled. Results were similar except that in all comparisons, inclusion of pre-mRNA made slight regressions in quantification accuracy, even in mRNA transcripts of genes with at least 10% pre-mRNA (MARD of 0.668 without and 0.700 with, and Spearman correlation of 0.799 without and 0.781 with).

We inspected genes whose quantification was most sensitive to the inclusion of pre-mRNA in the index and identified that these genes primarily displayed one of the following challenging characteristics: a large number of pseudo-genes (e.g. Rpl36a and other ribosomal proteins); annotated retained intron events (e.g. Slc35b2); overlap with an exon or intron of another gene on the same strand (e.g. Eno1b); or large overlap between an existing transcript and the included pre-mRNA (e.g. Mrpl41). However, these characteristics were common both in genes that saw improvements and those that saw regressions in quantification, indicating that they are markers of variability and not systematically affected by inclusion of the pre-mRNA index.

We conclude that inclusion of pre-mRNA in Salmon indexes has a relatively small effect and may hurt mature mRNA quantification. This leads to the recommendation to avoid pre-mRNA inclusion unless directly interested in pre-mRNA quantifications.

DISCUSSION

We demonstrate that the detailed simulation of molecule-level sequencing through BEERS2 provides valuable insights to RNA-seq analysis. In particular, it allows us to separate out individual aspects of the library preparation and sequencing to assess the impact of these steps on quantification.

In combination, CAMPAREE and BEERS2 allow simulations with substantially more detailed and configurable values. CAMPAREE includes sample-dependent features such as indel and single-nucleotide variants that are then consistent across all reads from a location, as well as a realistic distribution of reads including allele-specific expression. However, BEERS2 also remains flexible enough to accept custom input from sources other than CAMPAREE or to combine CAMPAREE-generated data and custom spike-in data. This makes BEERS2 highly adaptable for future benchmarking designs. BEERS2 allows configuration of every step of a realistic library preparation, allowing the effect of specific but known variables to be quantified.

When used with CAMPAREE, BEERS2 also supports sample-dependent genomic variants and allele-specific diploid expression. If the user prepares appropriate input transcript samples with differential expression already present, then BEERS2 also supports differential expression. All these factors contribute to creating more realistic simulated samples. Furthermore, BEERS2’s modular design makes it flexible enough to keep up with the fast-evolving technology of high-throughput sequencing. One future development could be the support of single-cell RNA-seq pipelines. These could reuse several of the sequencing steps implemented in BEERS2 (such as PCR amplification) but would also require implementing new steps such as for applying cell-specific barcodes. Likely the most complex part would be preparation of the input sample, analogous to the bulk RNA-seq samples generated by CAMPAREE. These would require realistic cell types, realistic populations of cells of those types and realistic transcript expression within each cell type, which remains a topic of active development even in spreadsheet-level simulations of small conditional RNA-seq [36].

One limitation of the BEERS2 model is that values are configured according to how each individual step behaves. However, real data are always the result of many steps done together. Therefore, configuration values do not directly correspond to measured values and so users may have difficulty determining parameters that give realistic values for their use-case. Moreover, BEERS can only account for effects that we can model mechanistically in the pipeline and the work of exhaustively understanding RNA-seq effects is ongoing and incomplete, meaning that our generated datasets never perfectly match real data. Incorporation of these effects, as they are characterized, will be an ongoing effort.

Key Points

BEERS2 provides detailed simulation of individual reads in RNA sequencing, including every step of a realistic library preparation and short-read sequencing experiment.

BEERS2 produces FASTQ and SAM/BAM files with a known origin of every read and known abundances of transcripts.

These results include known, configurable bias effects that create non-ideal data, including: 3′ bias from polyA-selection; GC bias during PCR amplification; start sequence biases from hexamer priming; sequence depletion from ribosomal RNA removal; errors in sequencing of barcodes; and copy errors during PCR.

BEERS2 is highly adaptable with a modular design that allows re-use of common sequencing components when simulating new library preparation techniques

Supplementary Material

BEERS2_12-17-2023_supplemental_bbae164

ACKNOWLEDGEMENTS

We would like to thank Cris Lawrence for her programming expertise, Kaitlyn Forrest for processing samples for sequencing and Jonathan Schug for insightful discussions about sequencing.

FUNDING

This work was supported by R21-LM012763-01A1: ‘The Next Generation of RNA-Seq Simulators for Benchmarking Analyses’ (PI: GR. Grant) and by the National Center for Advancing Translational Sciences Grant (5UL1TR000003), NHLBI R01HL155934-01A1(SS) and NHLBI-R01HL147472 (SS) and DP2GM146251 (PSC).

DATA AVAILABILITY

Input data used for the example applications have been previously published on the Gene Expression Omnibus repository under accession GSE98562. The generated simulated data are available for use at http://bioinf.itmat.upenn.edu/BEERS2/paper/index.html.

CODE AVAILABILITY

The code for the BEERS2 simulator is available under an open-source GPLv3 license at https://github.com/itmat/BEERS2. Code used for generating and analyzing the sample application of BEERS2 is available at https://github.com/tgbrooks/BEERS2_salmon_benchmark. The pipeline to optimize BEERS2 parameters to mimic a real sample is available at https://github.com/tgbrooks/BEERS2_optimize.

Author Biographies

Thomas G. Brooks is a research associate in bioinformatics at the Institute for Translational Medicine and Therapeutics at the University of Pennsylvania.

Nicholas F. Lahens is a bioinformatics researcher at the Institute for Translational Medicine and Therapeutics, University of Pennsylvania.

Antonijo Mrčela is a post-doc in bioinformatics at the Institute for Translational Medicine and Therapeutics, University of Pennsylvania.

Dimitra Sarantopoulou is a bioinformatician at the National Institute on Aging.

Soumyashant Nayak is an assistant professor at the Theoretical Statistics and Mathematics Unit of the Indian Statistical Institute.

Amruta Naik is a senior research investigator at the Gene Therapy Program, University of Pennsylvania.

Shaon Sengupta is an attending neonatologist in the Division of Neonatology at Children's Hospital of Philadelphia.

Peter S. Choi is an assistant professor of Pathology and Laboratory Medicine and an investigator with the Center for Childhood Cancer Research at Children's Hospital of Philadelphia.

Gregory R. Grant is a research associate professor of Genetics at the University of Pennsylvania and is the director of bioinformatics at the Institute for Translational Medicine and Therapeutics, University of Pennsylvania.
==== Refs
References

1. Froussios K , SchurchNJ, MackinnonK, et al. How well do RNA-Seq differential gene expression tools perform in a complex eukaryote? A case study in Arabidopsis thaliana. Bioinformatics 2019;35 (18 ):3372–7.30726870
2. Sandve GK , GreiffV. Access to ground truth at unconstrained size makes simulated data as indispensable as experimental data for bioinformatics methods development and benchmarking. Bioinformatics 2022;38 (21 ):4994–6.36073940
3. Bray NL , PimentelH, MelstedP, PachterL. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol 2016;34 (5 ):525–7.27043002
4. Patro R , DuggalG, LoveMI, et al. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods 2017;14 (4 ):417–9.28263959
5. Srivastava A , MalikL, SarkarH, et al. Alignment and mapping methodology influence transcript abundance estimation. Genome Biol 2020;21 :239.32894187
6. Baik B , YoonS, NamD. Benchmarking RNA-seq differential expression analysis methods using spike-in and simulation data. PloS One 2020;15 :e0232271.32353015
7. Li X , BrockGN, RouchkaEC, et al. A comparison of per sample global scaling and per gene normalization methods for differential expression analysis of RNA-seq data. PloS One 2017;12 :e0176185.28459823
8. Assefa AT , De PaepeK, EveraertC, et al. Differential gene expression analysis tools exhibit substandard performance for long non-coding RNA-sequencing data. Genome Biol 2018;19 . 10.1186/s13059-018-1466-5.
9. Rigaill G , BalzergueS, BrunaudV, et al. Synthetic data sets for the identification of key ingredients for RNA-seq differential analysis. Brief Bioinform 2018;19 (1 ):bbw092–76.
10. Benidt S , NettletonD. SimSeq: a nonparametric approach to simulation of RNA-sequence datasets. Bioinformatics 2015;31 (13 ):2131–40.25725090
11. Grant GR , FarkasMH, PizarroAD, et al. Comparative analysis of RNA-Seq alignment algorithms and the RNA-Seq unified mapper (RUM). Bioinformatics 2011;27 (18 ):2518–28.21775302
12. Kalari KR , NairAA, BhavsarJD, et al. MAP-RSeq: Mayo Analysis Pipeline for RNA sequencing. BMC Bioinformatics 2014;15 :224.24972667
13. Slaff B , RadensCM, JewellP, et al. MOCCASIN: a method for correcting for known and unknown confounders in RNA splicing analysis. Nat Commun 2021;12 :3353 PubMed PMID: 34099673; PMCID: PMC8184769.34099673
14. Li HD , FunkCC, PriceND. iREAD: a tool for intron retention detection from RNA-seq data. BMC Genomics 2020;21 . 10.1186/s12864-020-6541-0 PubMed PMID: 32028886; PMCID: PMC7006120.
15. Frazee AC , JaffeAE, LangmeadB, LeekJT. Polyester: simulating RNA-seq datasets with differential transcript expression. Bioinformatics 2015;31 (17 ):2778–84.25926345
16. Manz Q , TsoyO, FennA, et al. ASimulatoR: splice-aware RNA-Seq data simulation. Bioinformatics 2021;37 (18 ):3008–10.33647976
17. Li B , DeweyCN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 2011;12 :323.21816040
18. Griebel T , ZacherB, RibecaP, et al. Modelling and simulating generic RNA-Seq experiments with the flux simulator. Nucleic Acids Res 2012;40 (20 ):10073–83.22962361
19. Audoux J , SalsonM, GrossetCF, et al. SimBA: a methodology and tools for evaluating the performance of RNA-Seq bioinformatic pipelines. BMC Bioinformatics 2017;18 :428 PubMed PMID: 28969586; PMCID: PMC5623974.28969586
20. Caboche S , AudebertC, LemoineY, HotD. Comparison of mapping algorithms used in high-throughput sequencing: application to ion torrent data. BMC Genomics 2014;15 (1 ):264.24708189
21. Bushnell B. BBMap short read aligner, and other bioinformatic tools. [Access date 3/8/2023]. Available from: https://sourceforge.net/projects/bbmap/.
22. Ono Y , AsaiK, HamadaM. PBSIM: PacBio reads simulator—toward accurate genome assembly. Bioinformatics 2012;29 (1 ):119–21.23129296
23. Botond Sipos GS . Tim Massingham, Nick Goldman. Realistic simulations reveal extensive sample-specificity of RNA-seq biases arXiv. 2013; doi: https://arxiv.org/abs/1308.3172.
24. Sarkar A , StephensM. Separating measurement and expression models clarifies confusion in single-cell RNA sequencing analysis. Nat Genet 2021;53 (6 ):770–7.34031584
25. Lahens NF , BrooksTG, SarantopoulouD, et al. CAMPAREE: a robust and configurable RNA expression simulator. BMC Genomics 2021;22 (1 ):692.34563123
26. Molder F , JablonskiKP, LetcherB, et al. Sustainable data analysis with Snakemake. F1000Res 2021;10 :33.34035898
27. Hansen KD , BrennerSE, DudoitS. Biases in Illumina transcriptome sequencing caused by random hexamer priming. Nucleic Acids Res 2010;38 :e131 PubMed PMID: 20395217; PMCID: PMC2896536.20395217
28. Kao WC , StevensK, SongYS. BayesCall: a model-based base-calling algorithm for high-throughput short-read sequencing. Genome Res 2009;19 (10 ):1884–95. 10.1101/gr.095299.109.19661376
29. Love MI , HogeneschJB, IrizarryRA. Modeling of RNA-seq fragment sequence bias reduces systematic errors in transcript abundance estimation. Nat Biotechnol 2016;34 (12 ):1287–91.27669167
30. Kaminow B , YunusovD, DobinA. STARsolo: accurate, fast and versatile mapping/quantification of single-cell and single-nucleus RNA-seq data. bioRxiv. 2021 2021.05.05.442755. 10.1101/2021.05.05.442755.
31. Sarantopoulou D , BrooksTG, NayakS, et al. Comparative evaluation of full-length isoform quantification from RNA-Seq. BMC Bioinformatics 2021;22 (1 ):266.34034652
32. Lahens NF , RicciottiE, SmirnovaO, et al. A comparison of Illumina and ion torrent sequencing platforms in the context of differential gene expression. BMC Genomics 2017;18 :602 PubMed PMID: 28797240; PMCID: PMC5553782.28797240
33. Soneson C , LoveM, RobinsonM. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences [version 2; peer review: 2 approved]. F1000Research 2016;4 (1521 ):1521.
34. Soneson C , SrivastavaA, PatroR, StadlerMB. Preprocessing choices affect RNA velocity results for droplet scRNA-seq data. PLoS Comput Biol 2021;17 :e1008585 PubMed PMID: 33428615; PMCID: PMC7822509.33428615
35. He D , ZakeriM, SarkarH, et al. Alevin-fry unlocks rapid, accurate and memory-frugal quantification of single-cell RNA-seq data. Nat Methods 2022;19 (3 ):316–22.35277707
36. Cao Y , YangP, YangJYH. A benchmark study of simulation methods for single-cell RNA sequencing data. Nat Commun 2021;12 (1 ):6911.34824223
