RandomReadsMG

Basic Usage

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

Generate synthetic reads from all FASTA assemblies with random coverage levels.

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

Generate synthetic reads with custom coverage levels and realistic errors for specific genomes.

randomreadsmg.sh in=M.ruber.fa.gz reads=20k adderrors pacbio out=pacbio_reads.fq

Generate exactly 20,000 PacBio reads with errors (demonstrates 'k' suffix support for read counts).

Parameters

Parameters are organized by their function in the synthetic read generation process. Each assembly can be assigned random or custom coverage levels, and reads contain taxonomic information when filenames follow the 'tid_x_' convention.

Examples

# Generate Illumina paired-end reads from multiple assemblies
randomreadsmg.sh *.fa out=synthetic_reads.fq.gz

# Generate single-end reads with custom depth and errors
randomreadsmg.sh input.fa depth=50 paired=f adderrors out=single_reads.fq

# Generate exactly 20,000 PacBio reads with errors
randomreadsmg.sh in=M.ruber.fa.gz reads=20k adderrors pacbio out=pacbio_reads.fq

# Generate ONT reads with realistic error profile
randomreadsmg.sh genomes.fa ont meanlen=20000 adderrors=t out=ont_reads.fq.gz

# Custom coverage with specific genomes at different abundances
randomreadsmg.sh ecoli.fa=100 staph.fa=10 env_sample.fa=1 out=custom_coverage.fq.gz

# Add realistic coverage variation and PCR duplicates
randomreadsmg.sh genomes.fa sinewave numwaves=6 waveamp=0.8 pcr=0.05 out=realistic_reads.fq.gz

# Generate reads from taxonomically labeled genomes
randomreadsmg.sh tid_511145_ecoli_k12.fa tid_1280_staph_aureus.fa out=tax_reads.fq.gz

# Custom coverage by taxonomy ID
randomreadsmg.sh tid_*.fa cov_511145=50 cov_1280=25 out=custom_tax_reads.fq.gz

# Process large dataset with multiple threads
randomreadsmg.sh large_genomes/ threads=16 seed=12345 out=parallel_reads.fq.gz

# Single complex genome with multi-threading
randomreadsmg.sh complex_genome.fa threads=8 sinewave numwaves=10 out=complex_reads.fq.gz

Algorithm Details

Coverage Distribution Models

RandomReadsMG implements four coverage distribution modes to model realistic metagenomic abundance patterns:

• MIN4 (default): Takes the minimum of 4 random square-root values, creating a heavily skewed distribution toward low-abundance organisms typical in environmental samples.
• EXP: Uses exponential distribution with lambda=1/√(min_depth×max_depth), producing natural long-tailed abundance patterns.
• ROOT: Square-root distribution providing moderate skew toward lower abundances.
• LINEAR: Uniform distribution across the depth range for even coverage testing.

Platform-Specific Error Modeling

Each sequencing platform uses distinct error characteristics:

Illumina Error Model

Quality-based error introduction using PHRED scores. Generates random quality scores within the specified range (qavg ± qrange), then introduces substitutions probabilistically based on quality-to-error probability conversion. Adapter contamination is modeled when insert sizes are shorter than read length.

PacBio Error Model

Uses log-normal read length distribution with configurable sigma parameter (default 0.5). Error rates: substitutions 0.00015, insertions/deletions 0.000055/0.000045, homopolymer bonus 0.000015. Length generation applies bias correction: μ = log(mean) - 0.5σ².

ONT Error Model

Mixed distribution approach combining log-normal core with exponential heavy tail (controlled by tailfactor=0.2). Higher error rates: substitutions 0.0025, insertions/deletions 0.0055/0.0045, homopolymer bonus 0.02. Models the bimodal length patterns typical of ONT sequencing.

Sine Wave Coverage Modeling

The sinewave flag enables within-contig coverage variation modeling with realistic spatial bias patterns:

• Multiple sine waves: Combines numwaves overlapping sine waves with random periods between minperiod and maxperiod.
• Origin bias: Models replication origin effects with linear coverage decrease controlled by oribias parameter.
• Coverage bounds: minprob prevents coverage from dropping below a threshold, avoiding unrealistic assembly gaps.
• Wave amplitude: waveamp controls the maximum variation strength, with values near 0.7-0.8 producing realistic patterns.

Multi-threaded Architecture

ProcessThread workers use atomic coordination for work distribution. Each thread maintains independent random number generators seeded deterministically (seed + thread_id) for reproducible results. Thread-local statistics are accumulated atomically at completion. Single-file multi-threading splits depth targets across threads while maintaining total coverage.

Taxonomic Integration

Automatically parses taxonomy IDs from filenames following 'tid_x_' convention or read headers. Integrates with TaxTree for proper species classification. Read headers encode source information: f_[file]_c_[contig]_s_[strand]_p_[position]_i_[insert]_d_[duplicate]_tid_[taxid]. Enables benchmarking of metagenomic analysis tools with ground-truth data.

PCR Duplicate Modeling

Caches read start positions, insert sizes, and strand orientations to generate realistic PCR duplicates at the specified rate. Duplicate reads maintain identical genomic coordinates while receiving distinct read numbers. Single-end and paired-end duplicates are handled appropriately for each platform.

Notes

• Reads headers contain the TaxID of the originating genome if the filename starts with 'tid_x_' where x is a positive integer.
• Custom coverage can be set per file using cov_filename=depth or filename=depth syntax.
• The tool supports automatic file format detection and compression (gzip, bzip2).
• Memory usage scales with the number of concurrent genomes (maxConcurrentGenomes parameter).
• For deterministic output, use seed parameter with threads=1 or ensure consistent thread counts.
• Adapter contamination only occurs when insert size is shorter than read length in paired mode.
• Long-read error rates can be overridden for any platform, including Illumina.
• Coverage variation (variance parameter) is applied unless sine wave modeling is enabled.

Support

For questions and support:

• Email: bbushnell@lbl.gov
• Documentation: bbmap.org