Simulating Illumina metagenomic data with InSilicoSeq

1With the release of a growing number of bioinformatics tools, it has become challenging to know which tool performs best or is best suited for a particular experiment. The simulation of genomics and metagenomics data holds a prominent role both in the planning of an experiment and the development of new methods. On the contrary to real data, simulated data can be produced with controlled parameters, such as—in the case of metagenomics—the abundance of the species present in a sample. Fixing such parameters allows for benchmarking and testing of new tools in controlled conditions, as well as provides researchers with mock data for testing new tools or pipelines. Such an environment for testing is especially important for fast-growing sub-fields, such as metagenomics.

2Additionally, simulated data has proven to be very useful in the classroom, where teachers often need mock datasets that are small enough to be analysed quickly and yield meaningful and clear results that are easy to interpret for the students.

3Surprisingly, only a few such simulation software exist for metagenomics, and the existing solutions are often difficult or inconvenient to use as well as poorly maintained. Here, we describe InSilicoSeq, a software that simulates realistic Illumina reads from (meta)genomes. InSilicoSeq is multi-threaded, well-documented and easily installed via Python’s package manager pip. InSilicoSeq aims at making the benchmarking and testing of (meta)genomics software easier.

4InSilicoSeq is written in Python, can accurately model PHRED scores, supports substitution, insertion and deletion errors, as well as insert size distribution and GC bias. InSilicoSeq implements Kernel Density Estimation (KDE) to model base quality and insert size. Briefly, KDE is a non-parametric class of estimators that generally produces a smoother estimation of a distribution than histograms. Simulation of insertions, deletions, substitutions and GC bias is made empirically and calculated from aligned reads.

5In the current release, InSilicoSeq comes with pre-built error models for MiSeq, HiSeq and NovaSeq instruments, but we provide a command to generate error models from any bam file containing aligned reads. The provided error models are calculated from aligned reads in the bam format, generated from three datasets: PRJEB20178 for the MiSeq instrument and public datasets from Illumina Basespace for HiSeq and NovaSeq instruments. The three datasets were assembled with megahit and the reads were mapped back to the assembly using bowtie2 with default parameters.

6InSilicoSeq being designed for metagenomics, it will generate reads from multiple genomes according to a log-normal abundance distribution per default. Other distributions are built-in, as well as the possibility to provide the software with the exact abundance for each input genome.

7Existing software for simulating metagenomics include MetaSim, NeSSM, BEAR, FASTQSim, GemSim, Grinder, pIRS and FunctionSim (; https://cals.arizona.edu/∼anling/software/FunctionSIM.htm).

8We attempted to install and run all the aforementioned software as well as ART, a popular single genome simulator; of the nine tested software, only ART, Grinder and pIRS could be installed and run without issues. This is symptomatic of software development in several areas of science, including biology and was one of the main drivers behind the development of InSilicoSeq (Refer to Supplementary Material for more information on usability issues of the other simulators).

9InSilicoSeq can simulate half a million reads in under 10 min (Supplementary Fig. S1) using 4 CPUs and less than 1 G of RAM, and produces more realistic datasets than the other tested simulators. While pIRS ran under 1 min, Grinder took on average more than 13 h to generate half a million reads.

10Figure 1 shows the per-base quality distribution of forward reads simulated with InSilcioSeq, ART and pIRS compared to real data. InSilicoSeq and ART model very closely the base quality of the MiSeq dataset, while pIRS reports all the bases with a PHRED score of 40. For a figure including Grinder as well as reverse reads, refer to Supplementary Figure S2.

11One difficult part of generating realistic Illumina data is generating low-quality sequences. InSilicoSeq models this by clustering the sequences by mean quality before modelling the per-base quality distribution. Supplementary Figure S2 shows that our approach out-performs ART, grinder and pIRS: InSilicoSeq is the only simulator to produce low-quality sequences with a mean quality below 20.

12We developed a simulator that is free, open-source, well-tested, easy to install, has sufficient documentation and consists of a unified command (iss). InSilicoSeq produces realistic Illumina data with errors models based on recent Illumina machines and chemistry. New models can be produced from bam files in less than an hour, making it easy to keep them up to date. InSilicoSeq produces more realistic data than existing metagenomics simulation methods and is useful for planning experiments and benchmarking new methods.

13This work was supported by the Swedish Research Council, grant number 2015-03443_VR.

14Conflict of Interest: none declared.