pIRS

What it simulates

pIRS (profile-based Illumina paired-end reads simulator) simulates Illumina reads from a reference FASTA. It cuts fragments with a normal length distribution, then adds errors drawn from profiles measured on real sequencing data.

The inputs are a FASTA file, a read length, a coverage and an insert length. The outputs are two FASTQ files plus logs of every read, the insert lengths and the error positions.

Three empirical profiles shape the reads:

• A base-calling profile gives substitution errors and quality scores for each cycle. The bundled profile is humNew.PE100, built from human data with 100 bp reads.
• An indel profile gives sequencing-error insertions and deletions. The bundled one comes from PhiX.
• A GC-bias profile lowers the coverage of fragments with extreme GC content.

When to use it

Use pIRS when error and quality patterns must look like real Illumina data. Typical questions:

• Does my trimming or error-correction step work on position-dependent errors?
• How does GC bias change coverage in my target?
• Does my pipeline cope with realistic quality strings?

wgsim and DWGSIM use a flat error rate and plant variants for you. Use them when you need a truth set of SNPs. InSilicoSeq fits better when the input is several genomes at different abundances. pIRS has a separate pirs diploid command to make a variant copy of a genome, and that command is not part of the example.

Install it yourself

The Bioconda package is pirs 2.0.2. The recipe builds for linux/amd64 only, so on an ARM machine it runs under emulation.

conda create -n pirs -c conda-forge -c bioconda pirs=2.0.2
conda activate pirs

The minimal command takes a FASTA file and an output prefix:

pirs simulate -l 100 -x 5 -m 180 -v 18 -S 42 --threads=1 -o sim reference.fa

• -l sets the read length. The bundled profile covers up to 100 bp.
• -x sets the average coverage.
• -m sets the mean insert length. The default is 180.
• -v sets the insert standard deviation. The default is 10% of the mean.
• -S sets the seed.
• --threads=1 keeps the run reproducible. With more threads each thread uses its own generator.
• -o sets the output prefix.
• -e sets a substitution error rate. Without it the profile's own rate applies.
• --no-gc-bias, --no-indels and --no-subst-errors switch the three models off.
• -d marks a diploid input and halves the coverage, so give the two sequences from pirs diploid.

Reading the output

The reference is one 20,000 bp sequence called synthetic_ref. The run asks for 5x coverage with 100 bp reads. That gives 20,000 times 5 divided by 200 bases per pair, or 500 read pairs.

sim_100_180_1.fq and sim_100_180_2.fq each hold 500 records and are 109,892 bytes. The names in the file encode the mean insert, so sim_100_180 means reads of 100 bp and an insert of 180 bp. Read names look like @read_180_1/1. Quality characters use the Phred+33 offset and vary along each read, from H (Phred 39) down to 5 and below.

pirs_log.txt (4,577 bytes in the example) reports what was loaded and what happened:

• The base-calling profile has 4 reference bases, 200 cycles and 41 quality scores.
• The GC-bias profile humNew.gcdep_200.dat was chosen for this insert length. Its relative abundance is 0.78 at 20% GC, 0.97 at 30% GC and 0.06 at 80% GC.
• 500 of 500 read pairs were simulated, with 360 substitution errors, an average rate of 0.360%.
• There were 3 deletions and no insertions.
• 21.14% of fragments were affected by GC bias.

sim_100_180.insert_len.distr shows the realized insert lengths. The target was 180 with a standard deviation of 18. The actual mean was 180.384 and the actual standard deviation was 18.5168.

sim_100_180.read.info (74,531 bytes) has one line per read. The columns are read ID, reference file, contig, position, orientation, insert size, masked end length, substitutions, insertions and deletions. The line read_180_1/2 ... 18912 - 160 0 92,A->C; says read 2 of the first pair lies on the minus strand at position 18912, comes from a 160 bp insert and carries the substitution 92,A->C. This file is the truth for scoring a mapper.

sim_100_180.error_rate.distr lists the error frequency and rate for each of the 200 cycles. The first 10 cycles of the example have no errors.

The log prints pIRS version 2.0.0 although the package is 2.0.2.

Your own runs

These runs used the recipe's image and the example reference. Compare each with the baseline above.

• Higher coverage. Set -x 20. The run makes 2,000 pairs and the FASTQ files grow from 109,892 to 440,893 bytes. The average substitution error rate is 0.384% and the realized insert mean is 179.468.

• Longer inserts. Set -m 400 -v 40. The realized insert mean is 399.468. The log still loads humNew.gcdep_200.dat, so the GC profile did not change at this insert length. The FASTQ names become sim_100_400_1.fq.

• Your own error rate. Add -e 0.01. The average substitution error rate in the log becomes 1.247% against the 1.0% requested. The documentation warns that the realized rate can be off by 20% or more.

Pitfalls

• The bundled base-calling profile covers reads up to 100 bp. The recipe limits -l to 100.
• The realized error rate differs from -e. Read the final rate in pirs_log.txt and adjust.
• Several threads make a run non-reproducible. Use --threads=1 with -S.
• Coverage is for a haploid genome. For a diploid input pass -d, or you get twice the reads you wanted.
• Output names carry the read length and insert, as in sim_100_180_1.fq. Match them with sim_*_1.fq.

Further reading

• pIRS README, with all options: https://github.com/galaxy001/pirs
• Hu X et al. (2012). pIRS: profile-based Illumina pair-end reads simulator. Bioinformatics 28:1533-1535. https://doi.org/10.1093/bioinformatics/bts187