Add art

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 <a href="https://github.com/lh3/wgsim">wgsim</a> and <a href="https://github.com/seqan/seqan/blob/main/apps/mason2/README.mason_simulator">mason simulator</a> 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 <a href="https://github.com/lh3/wgsim">wgsim</a>, <a href="https://github.com/seqan/seqan/blob/main/apps/mason2/README.mason_simulator">mason simulator</a>, and <a href="https://www.niehs.nih.gov/research/resources/software/biostatistics/art">ART</a> 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
  ```

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(

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}")
 

ci/conda_requirements.txt

@@ -3,7 +3,8 @@ regex
 numpy
 pandas
 click
+art
 wgsim
 pytest
 minimap2
-freyja
+freyja

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