diff --git a/README.md b/README.md index 7cd74da..3285900 100644 --- a/README.md +++ b/README.md @@ -1,12 +1,12 @@ # Bygul: Amplicon & Metagenomics Read Simulator -**Bygul** is a Python 3 tool designed for simulating sequencing reads in wastewater surveillance and other metagenomic applications. It allows users to simulate complex multi-sample datasets with customizable proportions using industry-standard backends like `wgsim` and `mason`. +**Bygul** is a Python 3 tool designed for simulating sequencing reads in wastewater surveillance and other metagenomic applications. It allows users to simulate complex multi-sample datasets with customizable proportions using industry-standard backends like `wgsim`, `mason`, and ART Illumina. --- ## 🏗 Installation -Bygul requires **Python 3**. Since it relies on external simulators (`wgsim` and `mason`), we recommend using Conda to manage dependencies.For more info on wgsim and mason simulator please check their documentations. +Bygul requires **Python 3**. Since it relies on external simulators (`wgsim`, `mason`, and ART), we recommend using Conda to manage dependencies.For more info on wgsim, mason simulator, and ART please check their documentations. ### Option 1: Via Conda (Recommended) ```bash @@ -17,7 +17,7 @@ conda create -n bygul bioconda::bygul ```bash pip install bygul ``` -*Note: Some binary dependencies (wgsim/mason) may need to be installed manually or built from source if using this method.* +*Note: Some binary dependencies (wgsim/mason/ART) may need to be installed manually or built from source if using this method.* ### Option 3: Local Build from Source ```bash @@ -43,9 +43,10 @@ bygul simulate-proportions [SAMPLE1.fasta,SAMPLE2.fasta] --primers [primer.bed] bygul simulate-proportions sample1.fasta,sample2.fasta --primers primer.bed --reference reference.fasta --outdir results/ --maxmismatch 2 ``` * **Switching Simulators:** - Use `mason` instead of the default `wgsim`. + Use `mason` or `art` instead of the default `wgsim`. ```bash bygul simulate-proportions sample1.fasta,sample2.fasta --primers primer.bed --reference reference.fasta --simulator mason + bygul simulate-proportions sample1.fasta,sample2.fasta --primers primer.bed --reference reference.fasta --simulator art ``` * **Custom Error Rates & Lengths:** Pass simulator-specific parameters (e.g. indel fraction `-R`) directly. @@ -76,11 +77,13 @@ Bygul acts as a wrapper. While most flags are passed directly to the underlying - `--readcnt`: Number of reads per amplicon. - `--wgsim_insert_size`: Insert size for wgsim. - `--wgsim_read_length` / `--wgsim_error_rate`. +- `--art_read_length` / `--art_insert_size` / `--art_insert_sd` / `--art_seq_system`: paired-end ART Illumina defaults. To see all available backend flags, run: ```bash wgsim --help mason_simulator --help +art_illumina --help ``` Please note that some dependencies are not available through pypi. You need to install them using conda or build from source. @@ -106,9 +109,10 @@ Simulate reads from different samples without defining proportions (will be assi bygul simulate-proportions sample.fasta,sample2.fasta --primers primer.bed --reference reference.fasta --outdir results/ --maxmismatch 2 ``` Simulate reads with user-defined proportions and specifing read simulator. -bygul uses wgsim as a simulator but you can change it to mason. +bygul uses wgsim as a simulator but you can change it to mason or ART. ``` bygul simulate-proportions sample.fasta,sample2.fasta --primers primer.bed --reference reference.fasta --proportions 0.2,0.8 --simulator mason +bygul simulate-proportions sample.fasta,sample2.fasta --primers primer.bed --reference reference.fasta --proportions 0.2,0.8 --simulator art --art_read_length 150 ``` Simulate reads with user-defined proportions and number of reads per amplicon. ``` @@ -126,7 +130,7 @@ It is recommended to define the number of reads per amplicon to be greater than ### 🧬 Input BED File Format The pipeline expects a tab-delimited BED file where the first six columns represent standard genomic coordinates (`chrom`, `chromStart`, `chromEnd`, `name`, `poolName`, `strand`). Crucially, the fourth column (`name`) must follow a strict naming convention to prevent downstream parsing failures in variant-calling tools: **`[Scheme-Name]_[AmpliconNumber]_[Direction]_[OptionalSuffix]`** (e.g., `SARS-CoV-2_3_LEFT` or `SARS-CoV-2_3_LEFT_alt`). To ensure structural boundaries are parsed correctly, the prefix must not contain underscores, and any optional trailing modifiers must be restricted to standard alternative tags (`_alt`, `_ALT1`) or tracking indexes (`_0`, `_1`). Multi-level pool formatting, such as `SARS-CoV-2_400_1_LEFT_1`, is malformed and will fail validation.. The maximum number of mismatches allowed for each primer sequence is 1 SNP. To change this number, you may use the `--maxmismatches` flag. #### Complete set of available parameters -To learn more about how to adjust other parameters for the simulator please read the documentation for wgsim and mason simulator. Users can pass any simulator parameter directly in their command. The only parameters set through bygul are `--readcnt` and `--wgsim_insert_size` for amplicon sequencing mode. +To learn more about how to adjust other parameters for the simulator please read the documentation for wgsim, mason simulator, and ART. Users can pass any simulator parameter directly in their command. The only parameters set through bygul are `--readcnt`, `--wgsim_insert_size`, `--wgsim_read_length`, `--wgsim_error_rate`, `--art_read_length`, `--art_insert_size`, `--art_insert_sd`, and `--art_seq_system`. #### Simulated reads output Simulated reads from all samples are located in `provided_output_path/reads.fastq` #### Information about amplicon dropouts @@ -142,7 +146,8 @@ Simulate reads from different samples without defining proportions (will be assi bygul simulate-proportions sample.fasta,sample2.fasta --outdir results/ --simulation_mode metagenomics ``` -Specify proportions for each sample and add other simulator specific parameters. To access simulator parameters, please read wgsim and mason documentation. +Specify proportions for each sample and add other simulator specific parameters. To access simulator parameters, please read wgsim, mason, and ART documentation. ``` bygul simulate-proportions sample.fasta,sample2.fasta --proportions 0.5,0.5 --outdir results/ --simulation_mode metagenomics --simulator mason --illumina-read-length 200 +bygul simulate-proportions sample.fasta,sample2.fasta --proportions 0.5,0.5 --outdir results/ --simulation_mode metagenomics --simulator art --art_read_length 150 ``` diff --git a/bygul/_cli.py b/bygul/_cli.py index 3ef8c90..ec32940 100644 --- a/bygul/_cli.py +++ b/bygul/_cli.py @@ -51,8 +51,8 @@ def cli(): @click.option( "--simulator", default="wgsim", - type=click.Choice(["wgsim", "mason"], case_sensitive=False), - help="Select the simulator to use (wgsim or mason)", + type=click.Choice(["wgsim", "mason", "art"], case_sensitive=False), + help="Select the simulator to use (wgsim, mason, or art)", ) @click.option( "--simulation_mode", @@ -69,17 +69,33 @@ def cli(): @click.option( "--wgsim_insert_size", default=150, help="Outer distance for simulation using wgsim in amplicon" - "simulation mode." + " simulation mode." ) @click.option( "--wgsim_read_length", default=150, help="Read length for simulation using wgsim in amplicon" - "simulation mode." + " simulation mode." ) @click.option( "--wgsim_error_rate", default=0.0001, help="Error rate for simulation using wgsim in amplicon" - "simulation mode." + " simulation mode." +) +@click.option( + "--art_read_length", default=150, + help="Read length for simulation using ART Illumina." +) +@click.option( + "--art_insert_size", default=300, + help="Mean fragment size for paired-end simulation using ART Illumina." +) +@click.option( + "--art_insert_sd", default=10, + help="Fragment size standard deviation for ART Illumina." +) +@click.option( + "--art_seq_system", default="MSv3", + help="ART Illumina sequencing system/profile." ) @click.option("--readcnt", default=500, help="Number of reads per amplicon") @click.option( @@ -97,6 +113,10 @@ def simulate_proportions( wgsim_insert_size, wgsim_read_length, wgsim_error_rate, + art_read_length, + art_insert_size, + art_insert_sd, + art_seq_system, outdir, readcnt, maxmismatch, @@ -226,6 +246,10 @@ def simulate_proportions( wgsim_insert_size, wgsim_read_length, wgsim_error_rate, + art_read_length, + art_insert_size, + art_insert_sd, + art_seq_system, extra_flags=extra_simulator_flags ) # Merging logic remains the same @@ -253,6 +277,10 @@ def simulate_proportions( os.path.join(outdir, name, "reads"), cnt, simulator, + art_read_length, + art_insert_size, + art_insert_sd, + art_seq_system, extra_flags=extra_simulator_flags ) read_path1 = os.path.join( diff --git a/bygul/utils.py b/bygul/utils.py index 533f16f..a964568 100644 --- a/bygul/utils.py +++ b/bygul/utils.py @@ -419,6 +419,10 @@ def run_simulation_on_fasta( wgsim_insert_size, wgsim_read_length, wgsim_error_rate, + art_read_length=150, + art_insert_size=200, + art_insert_sd=10, + art_seq_system="MSv3", extra_flags=None ): """Runs simulator on a single FASTA file with the given parameters.""" @@ -470,7 +474,7 @@ def run_simulation_on_fasta( if extra_flags: command.extend(extra_flags) - else: + elif simulator == "mason": # Adjust Mason command command = [ "mason_simulator", @@ -485,6 +489,37 @@ def run_simulation_on_fasta( ] if extra_flags: command.extend(extra_flags) + + elif simulator == "art": + art_output_prefix = os.path.join( + output_dir, f"{output_prefix}_contig{contig_idx + 1}_" + ) + output1 = f"{art_output_prefix}1.fq" + output2 = f"{art_output_prefix}2.fq" + command = [ + "art_illumina", + "-ss", + str(art_seq_system), + "-p", + "-na", + "-i", + fasta_file, + "-l", + str(art_read_length), + "-c", + str(int(reads_per_contig)), + "-m", + str(art_insert_size), + "-s", + str(art_insert_sd), + "-o", + art_output_prefix, + ] + if extra_flags: + command.extend(extra_flags) + + else: + raise ValueError(f"Unsupported simulator: {simulator}") # Run the simulator command and capture any errors try: subprocess.run( @@ -502,6 +537,10 @@ def run_simulation_on_fasta_single_genome( output_dir, read_cnt, simulator, + art_read_length=150, + art_insert_size=300, + art_insert_sd=10, + art_seq_system="MSv3", extra_flags=None ): """Runs simulator on a single FASTA file with the given parameters.""" @@ -524,7 +563,7 @@ def run_simulation_on_fasta_single_genome( ] if extra_flags: command.extend(extra_flags) - else: + elif simulator == "mason": # Adjust Mason command command = [ "mason_simulator", @@ -539,10 +578,40 @@ def run_simulation_on_fasta_single_genome( ] if extra_flags: command.extend(extra_flags) + elif simulator == "art": + art_output_prefix = os.path.join(output_dir, "reads_") + art_output1 = f"{art_output_prefix}1.fq" + art_output2 = f"{art_output_prefix}2.fq" + command = [ + "art_illumina", + "-ss", + str(art_seq_system), + "-p", + "-na", + "-i", + fasta_file, + "-l", + str(art_read_length), + "-c", + str(int(read_cnt)), + "-m", + str(art_insert_size), + "-s", + str(art_insert_sd), + "-o", + art_output_prefix, + ] + if extra_flags: + command.extend(extra_flags) + else: + raise ValueError(f"Unsupported simulator: {simulator}") # Run the simulator command and capture any errors try: subprocess.run( command, check=True, capture_output=True, text=True) + if simulator == "art": + shutil.move(art_output1, output1) + shutil.move(art_output2, output2) except subprocess.CalledProcessError as e: print(f"An error occurred while running the command: {e}") diff --git a/ci/conda_requirements.txt b/ci/conda_requirements.txt index 39e2c4b..64adc94 100644 --- a/ci/conda_requirements.txt +++ b/ci/conda_requirements.txt @@ -3,7 +3,8 @@ regex numpy pandas click +art wgsim pytest minimap2 -freyja \ No newline at end of file +freyja diff --git a/environment.yml b/environment.yml index e83b6dd..b8f080e 100644 --- a/environment.yml +++ b/environment.yml @@ -20,6 +20,7 @@ dependencies: - libuuid=2.38.1=h0b41bf4_0 - libxcrypt=4.4.36=hd590300_1 - libzlib=1.2.13=h4ab18f5_6 + - art=2016.06.05 - mason=2.0.9=hdcf5f25_3 - ncurses=6.5=he02047a_1 - openssl=3.3.2=hb9d3cd8_0