RandomReadsMG

Script: randomreadsmg.sh Package: synth Class: RandomReadsMG.java

Generates synthetic reads from a set of FASTA assemblies, assigning each assembly its own coverage depth so that one run produces a mixed community with known abundance ratios. Read headers record the source file, contig, strand, position, insert size, and taxID (when present) for use in accuracy evaluation.

Overview

RandomReadsMG is the multi-genome version of RandomReads. Where randomreads.sh draws reads from one reference, randomreadsmg.sh takes many assemblies (one genome per file, or several contigs per file) and gives each file a coverage depth. Depths can be set explicitly, drawn from a distribution, or derived from a target read count.

Use it to:

Three read models are available: Illumina short reads (single-end or paired, the default), PacBio HiFi-style long reads, and ONT-style long reads. Illumina mode can add substitution errors, adapter contamination, and PCR duplicates. Long-read modes can add platform-style substitution, insertion, and deletion errors, and can embed an exact truth CIGAR in each header. PCR duplicates are available in every mode.

Each read header has the form f_(file number)_c_(contig)_s_(strand)_p_(position)_i_(insert)_r_(reference length), followed by d_(duplicate flag) when PCR duplicates are enabled, then tid_(taxID) when a taxID is found or name_(filename) otherwise, and optionally cigar_(CIGAR). All numbers are 0-based. Paired reads end with a pair number such as " 1:" or " 2:".

Basic Usage

randomreadsmg.sh in=<files> out=<file> [options]
randomreadsmg.sh *.fa out=reads.fq.gz
randomreadsmg.sh ecoli.fa=40 ruber.fa=0.1 phix.fa=10 out=reads.fq.gz

Input can be one file, a directory, a comma-delimited list, or bare file names with no "=" sign. Each input file is a separate genome. The first form gives every file a random depth. The second form sets the depth of each named file directly with <file>=depth.

With no mode or platform flag, reads are paired Illumina reads of 150 bp with an average insert of 300 bp. With only out= given, paired output is interleaved.

Parameters

Parameters are grouped as in the shell script usage text. The last group lists options the Java class accepts that the usage text does not mention.

File Parameters

in=<file,file>
Assembly input. Can be a single file, a directory of files, or a comma-delimited list. Unrecognized arguments with no "=" sign are also treated as input files if they exist and are readable. Default: none (required).
out=<file>
Synthetic read output destination. If paired reads are generated and only out= is given, output is interleaved. Default: null (no output written).
out2=<file>
Read 2 output, when twin files are wanted for paired reads. Requires out=. Default: null.

Depth and Processing Parameters

mindepth=1
Minimum average depth assigned to an assembly when depths are drawn randomly. Alias: mincov=.
maxdepth=256
Maximum average depth assigned to an assembly when depths are drawn randomly. Alias: maxcov=.
depth=
Sets mindepth and maxdepth to the same value, so every file gets this depth. Alias: cov=. Default: unset (mindepth=1, maxdepth=256).
reads=-1
If positive, set each file's depth from its size so that it yields approximately this many reads (read pairs when paired). Requires reading each input file twice: once to count bases, once to generate reads. Files with a custom depth (cov_x or x=) keep that depth. Custom depths set via cov_x or x= take precedence but are still modified by jitter when it is positive; only reads=-derived depths skip jitter.
readspercontig=-1
If positive, ignore depth and generate this many reads from each contig. The base budget is this count times the read length (times 2 when paired).
mode=min4
Distribution used to draw random depths: min4, exp, root, or uniform. Uniform is the same as linear. Bare words (e.g. exp) are also accepted. See Algorithm Details for the formulas.
cov_x=
Set a custom depth for the file named x. The name is matched against the file name with its path removed (cov_foo.fa=5 for foo.fa). x can instead be a taxID, matched when the file name starts with tid_x_ (cov_7=5 for tid_7_foo.fa). Alias: depth_x=.
<file>=x
Alternate way to set a custom depth. The file named on the left (which must exist) gets depth x. For example, ecoli.fa=40.
circular=f
Treat each contig as circular. Each contig's sequence is doubled before reads are drawn, so reads can span the origin. Depth is halved to compensate for the doubled sequence.
threads=
Maximum number of threads. Default: logical core count. Each input file uses one thread. The total is the smallest of threads, the number of files, and mcg. With exactly one input file and seed set to -1, this value also enables multithreaded processing of that file. See Technical Notes.
seed=-1
If non-negative, use this RNG seed. Output is then deterministic regardless of thread count: each file's reads are a function of (seed, seed2, file name). With seed=-1, output varies between runs.

Multi-sample Simulation Parameters

These options generate several correlated samples from one genome set, for testing differential-coverage binners. Each sample is one invocation.

depthseed=-1
Seed for per-file base depth assignment. -1 uses seed. Invocations that share a depthseed, but have different seeds, draw the same base depth for each genome, so the samples are correlated.
jitter=0.0
Multiply each file's depth by a random factor exp(u), where u is uniform in [-jitter, +jitter]. The variation is symmetric in log space: jitter=0.1 gives roughly +/-10% per-sample variation around the base depth. With jitter=0 and a shared depthseed, samples have identical depths.
seed2=-1
Optional separate seed for the jitter stream. -1 uses seed. Use it to reproduce one specific jitter pattern independently of the generation seed. It is also mixed into the generation stream, so samples that share a depthseed but differ in seed or seed2 produce independent reads.
zeroprob=0
Probability that each genome is absent (depth 0) from a sample. The decision is keyed on depthseed, so samples that share a depthseed share their presence/absence pattern, and independent depthseeds draw independently. Applies only to randomly chosen depths, not to cov_x= or reads=.
zeroprobbias=0
Biases absence probability by input-file rank. 0 gives every file the flat zeroprob. 1 ramps linearly so the last file is about twice as likely to be absent as average, and the first about never absent. The mean stays at zeroprob. Files are ranked in argument order.

Example: four correlated samples with about 10% depth variation, drawn from the same community:

for S in 1 2 3 4; do
  randomreadsmg.sh config=genomes.txt depthseed=777 seed=$S jitter=0.1 \
    out=sample$S.fq.gz
done

Artifact Parameters

pcr=0.0
PCR duplicate rate, 0 to 1. A duplicate reuses the start position, strand, and insert of the previous read on the same contig, and its header carries d_1. Originals carry d_0. Alias: pcrrate=.
randomkmer=f
Bias read start sites with random k-mer priming. Giving kprime a numeric value also turns this on. Aliases: randompriming=, hexamer=, kmerprime=.
kprime=6
Length k of the random primer k-mer. Giving a numeric value sets randomkmer=true. Aliases: kmerprime=, hexamer=.
kpower=0.5
Raises the linear primer distribution to this power. Must be greater than 0. Higher powers increase priming bias. Aliases: kmerpower=, kexp=, kmerexp=.
minkprob=0.1
Minimum primer k-mer probability. Alias: minkmerprob=.

Platform Parameters

illumina
Use Illumina length and error model. This is the default. Bare word, or platform=ILLUMINA (alias ILL).
pacbio
Use PacBio HiFi length and error model. Bare word, or platform=PACBIO (alias PB).
ont
Use ONT length and error model. Bare word, or platform=ONT (alias NANOPORE).
paired=true
Generate paired reads. Applies only to Illumina mode. PacBio and ONT modes ignore this flag and write single reads.
length=150
Read length for Illumina mode. Aliases: len=, readlen=, readlength=. In long-read mode it also sets the mean read length (see Technical Notes). Default: 150.
avginsert=300
Average insert size for paired reads. Each insert is drawn from a Gaussian with mean avginsert and standard deviation 0.25 times avginsert. Alias: insert=.

Long-read Length Parameters

minlen=1000
Minimum read length for PacBio and ONT modes. Alias: minlength=. Contigs shorter than minlen+50 are skipped.
meanlen=15000
Mean read length for PacBio and ONT modes. Also sets readlen, the same field that length= sets (see Technical Notes). Alias: avglen=.
maxlen=100000
Maximum read length for PacBio and ONT modes. Alias: maxlength=.
tailfactor=0.2
ONT only. Probability that a read length is drawn from the heavy exponential tail rather than the log-normal core. Alias: ontLongTailFactor=.
pbsigma=0.5
PacBio only. Standard deviation of the log-normal read length distribution, in log space. Alias: pacbiosigma=.

Error Parameters (All Platforms)

adderrors=f
Add model-specific errors. For Illumina, substitutions are drawn from per-base quality scores. For PacBio and ONT, substitution, insertion, and deletion errors are applied using the rates in the long-read error parameters. Without this flag, long reads carry no model errors. Alias: addErrors=.
subrate=0.0
Add substitutions at this rate per base, independent of platform model. Alias: snprate=.
insrate=0.0
Add length-1 insertions at this rate, independent of platform model.
delrate=0.0
Add length-1 deletions at this rate, independent of platform model.
indelrate=
Sets insrate and delrate each to half of this value.

Illumina-specific Parameters

illuminanames=f
Generate Illumina-format headers of the form instrument:number:number:flowcell:lane:tile:x:y, then read number, filter flag, control bits, and barcode. The mate's header is derived from the first read's header with the read number changed to 2. Alias: illuminaHeaders=.
qavg=25
Average quality score, used for generating Illumina errors. Aliases: qscore=, avgq=, avgqual=.
qrange=0
Quality score range. Each base's quality is drawn uniformly from qavg-qrange to qavg+qrange.
qflat=f
Use one quality value for all bases in a read, to improve compression when qrange>0. Aliases: flat=, flatqual=.
addadapters=f
Add adapter sequence to paired reads whose insert is shorter than the read length. Bases beyond the adapter are set to G at quality 8 or lower.
adapter1=
Custom R1 adapter, as observed in R1. Default: the standard Illumina R1 adapter. Aliases: fragadapter=, fragadapter1=, adapter=.
adapter2=
Custom R2 adapter, as observed in R2. Default: the standard Illumina R2 adapter. Alias: fragadapter2=.
barcode=
Barcode written into Illumina headers. Default: a random 10+10 base barcode.
machine=
Machine name written into Illumina headers. Default: "LH000" followed by a random number.

Long-read Error Parameters

These apply in addition to subrate, insrate, and delrate, and only when adderrors=t. They can be overridden for any platform.

srate=-1
Substitution rate per base. Default: 0.0025 for ONT, 0.00015 for PacBio.
irate=-1
Insertion rate per base. Default: 0.0055 for ONT, 0.000055 for PacBio.
drate=-1
Deletion rate per base. Default: 0.0045 for ONT, 0.000045 for PacBio.
hrate=-1
Homopolymer error boost. The indel chance rises by this much for each additional base of the same identity in a run. Default: 0.02 for ONT, 0.000015 for PacBio.
addcigar=f
Append an exact extended CIGAR (=, X, I, D, M) to PacBio and ONT read headers. The CIGAR is in forward-reference order. Requires a long-read platform, and cannot be combined with circular=t. Alias: cigar=. Default false, which keeps the legacy header format.

Coverage Variation Parameters

sinewave=f
Enable coverage variation within contigs, using a sum of sine waves with an origin-of-replication trend. Alias: wavecov=. Setting waves to a positive number also enables this mode.
waves=4
Number of sine waves to combine. More waves create more irregular peaks and valleys. Aliases: sinewaves=, numwaves=. Any positive value also enables sinewave mode.
waveamp=0.70
Maximum variation in coverage from the sine waves. Higher values (0 to 1) give larger differences between high and low coverage regions. Alias: maxamplitude=.
oribias=0.25
Strength of the origin-of-replication bias: the maximum linear decrease in coverage from the start to the end of a contig.
minprob=0.10
Minimum coverage probability, as a fraction of the target. Keeps coverage from dropping so low that gaps appear in the assembly. Alias: minwaveprob=.
minperiod=2k
Minimum sine wave period, in bp. Accepts K and M suffixes. Default: 2000.
maxperiod=80k
Maximum sine wave period, in bp. Accepts K and M suffixes. Default: 80000.
variance=0.0
Variation in coverage between contigs within an assembly, as plus or minus this factor. Alias: depthvariance=. In the current source this is realized as a random linear decline in coverage along each contig, not as a per-contig depth ratio. It is ignored when sinewave mode is on. See Technical Notes.

Java Parameters

-Xmx
Set Java heap memory limit. Format: -Xmx20g (20 GB) or -Xmx200m (200 MB). Overrides autodetection. The maximum is typically 85% of physical memory.
-eoom
Exit when an out-of-memory exception occurs. Requires Java 8u92+.
-da
Disable Java assertions.

Additional Options Accepted by the Source

The following keys are parsed by RandomReadsMG.java but do not appear in the shell script usage text. Defaults are from the source.

platform=ILLUMINA
Platform by name: ILLUMINA (or ILL), ONT (or NANOPORE), PACBIO (or PB). Same effect as the bare-word flags.
ref=
Alias for in=. Input files are the same.
verbose=f
Print extra progress messages, such as stream start and end.
loud=t
Print each file's depth to stderr, as "File n, name ...: depth=...". Turned off automatically in single-file multithreaded mode.
tree=f
Load the taxonomy tree (TaxTree) for taxID handling. Alias: usetree=. Default false.
sfthreads=1
Threads used for a single input file. Aliases: singleFileThreads=, sft=. Only used when seed is negative. Raises the global thread count to at least this value.
mcg=2048
Maximum number of input files processed concurrently. Aliases: concurrency=, maxConcurrentGenomes=.
lane=
Lane number written into Illumina headers. Default: random.
readlen=, avglen=, readlength=, len=
Aliases for length=. Each sets both the read length and the long-read mean length.
depthmode=
Alias for mode=.
interleaved=, int=
Aliases for paired=.
jitterseed=, depthjitter=, seed2variance=, absentprob=, zerobias=
Aliases for seed2=, jitter=, jitter=, zeroprob=, and zeroprobbias= respectively.
cigar=
Alias for addcigar=.
pacbiosigma=, ontLongTailFactor=, pcrrate=
Aliases for pbsigma=, tailfactor=, and pcr= respectively.

Examples

Mock Community with Default Depths

randomreadsmg.sh *.fa out=reads.fq.gz

Each FASTA file receives a random depth between 1x and 256x (MIN4 distribution). Output is interleaved paired Illumina reads of 150 bp, with 300 bp average inserts, written to a gzipped FASTQ.

Explicit Depths per Genome

randomreadsmg.sh ecoli.fa=40 ruber.fa=0.1 phix.fa=10 out=reads.fq.gz

Sets ecoli.fa to 40x, ruber.fa to 0.1x, and phix.fa to 10x. The three files must exist.

Single-end Reads from a Directory

randomreadsmg.sh in=genomes/ out=reads.fq.gz paired=f length=100 mindepth=5 maxdepth=50 mode=exp seed=42

Reads every FASTA file in genomes/. Depths are drawn from an exponential distribution between 5x and 50x. Single-end 100 bp reads, with a fixed seed so the output is reproducible.

Target Read Count per Genome

randomreadsmg.sh in=genomes/ out=reads.fq.gz reads=1000000 seed=42

Scales each file's depth so that it yields about 1,000,000 read pairs, since paired is on by default. Each file is read twice.

PacBio HiFi Long Reads

randomreadsmg.sh in=genomes/ out=hifi.fq.gz pacbio depth=20 meanlen=15000 minlen=1000 maxlen=50000 seed=7

Generates HiFi-style reads with a log-normal length distribution (mean 15 kb), at 20x for each genome. The HiFi error rates apply only when adderrors=t is added.

ONT Reads with Errors and Truth CIGAR

randomreadsmg.sh in=genomes/ out=ont.fq.gz ont depth=10 meanlen=8000 adderrors=t addcigar=t seed=3

Generates ONT-style reads with platform error rates applied and an exact truth CIGAR in each header (cigar_ field). The CIGAR can be used to grade alignments at the indel level.

Correlated Multi-sample Simulation

for S in 1 2 3 4; do
  randomreadsmg.sh config=genomes.txt depthseed=777 seed=$S jitter=0.1 \
    out=sample$S.fq.gz
done

Four samples share one base depth per genome (depthseed=777), and each gets about 10% independent per-sample depth variation (jitter=0.1). Each sample has different reads, because its seed differs.

Sparse Community

randomreadsmg.sh in=genomes/ out=sparse.fq.gz zeroprob=0.3 zeroprobbias=1 depthseed=11 seed=1

About 30% of genomes are absent from the sample. Absence is more likely for genomes later in the file order, so the same genomes tend to be missing across samples that share a depthseed.

Algorithm Details

The following descriptions are taken from RandomReadsMG.java. Line numbers refer to that file.

Per-file Depth Assignment

Depth Distributions

Contig Processing

Illumina Read Generation

Long-read Generation

PCR Duplicates

Within-contig Coverage Variation

Header Format

Threading and Determinism

Memory

RandomReadsMG does not build a genome index. Each file is streamed through a Streamer in batches of contigs (lines 649-655, 1339-1368). Memory use is set by the Java heap (see -Xmx), by the number of concurrent files, and by the largest contig, which is copied when circular=t. The Java source contains no explicit memory estimate.

Technical Notes

Discrepancies Between the Usage Text and the Source

These differences were found when reading RandomReadsMG.java. The source is treated as correct.

Parameters Handled by Shared Classes

out=, out2=, overwrite=, append=, and extension overrides are passed to the shared Parser class (line 358). Its handling was not read for this page. The usage text documents out= and out2=. The field defaults are overwrite=true and append=false (lines 1911-1913), but Parser values are copied over them (lines 140-147).

Platform Defaults

When a long-read platform is set, the substitution, insertion, deletion, and homopolymer rates are set from the platform defaults unless given explicitly (lines 449-466). Shared.FAKE_QUAL is set to 20 for ONT and 35 for PacBio (lines 456, 462).

Constraints

Performance

Not Verified for This Page

Support

For questions and support: