Skip to content

Latest commit

 

History

History
426 lines (302 loc) · 19.3 KB

File metadata and controls

426 lines (302 loc) · 19.3 KB

SILOsim Simulator

Models bacterial evolution through substitution, gene gain/loss, and recombination events.

silosim simulator -h
usage:  simulator [-h] --ancestor ANCESTOR --years YEARS --taxa TAXA --lambda_rate LAMBDA_RATE --mu_rate MU_RATE
                  --substitution_model SUBSTITUTION_MODEL --model_parameters MODEL_PARAMETERS --mutation_rate
                  MUTATION_RATE [--use_weighted_mutation USE_WEIGHTED_MUTATION] --weighted_mutation_file
                  WEIGHTED_MUTATION_FILE [--use_gain_loss USE_GAIN_LOSS] --gain_rate GAIN_RATE --loss_rate
                  LOSS_RATE --bin BIN --position_coverage POSITION_COVERAGE --mge_data MGE_DATA --mge_fasta
                  MGE_FASTA --mge_entropy MGE_ENTROPY [--use_recombination USE_RECOMBINATION] --recombination_rate
                  RECOMBINATION_RATE --min_recombination_size MIN_RECOMBINATION_SIZE --mean_recombination_size
                  MEAN_RECOMBINATION_SIZE --nu NU [--export-per-event-genomes] [--output OUTPUT] [--seed SEED]
                  [--prefix PREFIX] [--conda_prefix CONDA_PREFIX] [--force]

options:
  -h, --help            show this help message and exit
  --ancestor ANCESTOR   Path to ancestral genome in FASTA format.
                        Generated by SILOSim Profiler as {genome}_concat.fasta.
  --years YEARS         Number of years to simulate evolution.
  --taxa TAXA           Number of taxa(offspring) at the end of the simulation.
  --lambda_rate LAMBDA_RATE       Speciation rate (lambda) for the birth-death tree simulation.
                        Rate at which new sub-lineages arise per existing lineage per year.
                        No default value. Must be a positive number.
  --mu_rate MU_RATE     Extinction rate (mu) for the birth-death tree simulation.
                        Rate at which sub-lineages are lost per lineage per year.
                        When mu = 0, reduces to a pure-birth (Yule) process.
                        No default value. Must be a non-negative number.
  --substitution_model SUBSTITUTION_MODEL
                        Substitution model to use for simulation. 
                        Default is GTR.
                        Options are JC69, K2P, K3P, GTR.
  --model_parameters MODEL_PARAMETERS
                        Substitution model parameters (comma-separated, no spaces).
                        
                        - JC69: No parameters required (leave empty)
                        
                        - K2P: Single float value
                          - Transition/transversion ratio (kappa) > 0
                        
                        - K3P: Three comma-separated floats (all > 0)
                          - Position 1: Transition rate (alpha)
                          - Position 2: Transversion rate beta (A↔C, G↔T)
                          - Position 3: Transversion rate gamma (A↔T, C↔G)
                        
                        - GTR: Six comma-separated floats for substitution rates
                          - Position 1: A↔G
                          - Position 2: A↔C
                          - Position 3: A↔T
                          - Position 4: G↔C
                          - Position 5: G↔T
                          - Position 6: C↔T
  --mutation_rate MUTATION_RATE
                        Mutation rate (substitutions per site per year) for simulation
  --use_weighted_mutation USE_WEIGHTED_MUTATION
                        Whether to use weighted mutation. Default is True.
                        To disable weighted mutation, set to False.
  --weighted_mutation_file WEIGHTED_MUTATION_FILE
                        Path to weighted mutation sites CSV file (position,entropy columns).
                         Generated by SILOSim Profiler as {genome}_entropy.csv.
                         Required when --use_weighted_mutation is True.
  --use_gain_loss USE_GAIN_LOSS
                        Whether to simulate gene gain/loss events. Default is True.
                        To disable gene gain/loss simulation, set to False.
  --gain_rate GAIN_RATE
                        Gene gain rate (λ) for Poisson process model of gain events.
                        Required when --use_gain_loss is True.
  --loss_rate LOSS_RATE
                        Gene loss rate (λ) for Poisson process model of loss events.
                        Required when --use_gain_loss is True.
  --bin BIN             Path to genome bin summary CSV file.
                        Generated by SILOSim Profiler as {genome}_bin_summary.csv.
                        Required when --use_gain_loss is True.
  --position_coverage POSITION_COVERAGE
                        Path to position coverage RDS file. 
                        Generated by SILOSim Profiler as {genome}_position_coverage.RDS. 
                        Required when --use_gain_loss is True.
  --mge_data MGE_DATA   Path to mobile genetic element (MGE) data RDS file. 
                        Generated by SILOSim Profiler as {genome}_mges.RDS. 
                        Required when --use_gain_loss is True.
  --mge_fasta MGE_FASTA
                        Path to MGE sequence database in FASTA format. 
                        Generated by SILOSim Profiler as {genome}_mges_seq.fasta. 
                        Required when --use_gain_loss is True.
  --mge_entropy MGE_ENTROPY
                        Path to directory containing entropy CSV files for all MGEs in the database. 
                        Generated by SILOSim Profiler as {genome}_mge_entropy/ directory. 
                        Required when --use_gain_loss is True.
  --use_recombination USE_RECOMBINATION
                        Whether to simulate recombination events. Default is True.
                        To disable recombination simulation, set to False.
  --recombination_rate RECOMBINATION_RATE
                        Recombination rate (event per genome per year) for simulation.
                        Required when --use_recombination is True.
  --min_recombination_size MIN_RECOMBINATION_SIZE
                        Minimum recombination size (bp) for simulation. Default is 1 bp. 
                        Required when --use_recombination is True.
  --mean_recombination_size MEAN_RECOMBINATION_SIZE
                        Mean recombination size (bp) for simulation.
                        Required when --use_recombination is True.
  --nu NU               Nu parameter (snps / total length of recombination in bp) for recombination simulation.
                        Required when --use_recombination is True.
  --export-per-event-genomes
                        Export a genome snapshot after each individual evolutionary event
                        (substitution, recombination, or gene gain/loss if applicable) is applied,
                        rather than only after all events in an interval have been accumulated.
                        By default, genomes are exported once per interval after all events are
                        applied. This option enables fine-grained inspection of sequential
                        evolutionary changes but may significantly increase output size and runtime.
  --output OUTPUT       Output directory for the simulation results. If not provided, defaults to silosim_simulator_output_<YYYY_MM_DD_HHMMSS> under the current working directory.
  --seed SEED           Random seed for reproducibility.
                                default is the current date in the format of YYYYMMDD
  --prefix PREFIX       Prefix for config file, output files, and analysis naming. If not provided, defaults to timestamp: YYYY_MM_DD_HHMMSS
  --conda_prefix CONDA_PREFIX
                        Directory for conda environments needed for this analysis. If not provided, defaults to <OUTPUT>/conda_envs_<YYYY_MM_DD_HHMMSS>
  --force               Bypass operating system compatibility checks and force execution of the pipeline.
                        This may cause instability or failures.

Required Input

All required input files are generated by SILOSim Profiler. While custom files following the same format are accepted, using SILOSim Profiler output is strongly recommended to ensure compatibility and correctness.

  1. Ancestor Genome(--ancestor):

    Path to the concatenated ancestral genome in FASTA format for simulation.

    Generated by SILOSim Profiler as profiler/{genome}_concat.fasta.

  2. Simulation Length(--years):

    Number of years to simulate evolution.

  3. Number of Taxa(--taxa):

    Number of taxa at the end of the simulation.

  4. Lambda(--lambda_rate):

    Speciation rate for the birth-death tree simulation. This is the rate at which new sub-lineages arise per existing lineage per year. Higher values produce more branching events and a bushier tree.

    No default value. Must be a positive number.

  5. Mu(--mu_rate):

    Extinction rate for the birth-death tree simulation. This is the rate at which existing sub-lineages are lost per lineage per year, reflecting lineage extinction. When mu = 0, the model reduces to a pure-birth (Yule) process where all lineages persist.

    No default value. Must be a non-negative number.

    Note on lambda and mu: The simulated tree represents the evolutionary history of a strain lineage, capturing diversification regardless of whether it occurs within a single host or across multiple hosts through transmission. R0 = lambda / mu (undefined when mu = 0) must exceed 1 for the lineage to persist, but R0 alone does not fully characterize the simulation. Two parameter sets with the same R0 but different absolute magnitudes produce different trees: higher absolute rates generate shorter branch lengths (less sequence divergence between tips), while lower absolute rates generate longer branch lengths (more sequence divergence). Both R0 and the absolute magnitudes of lambda and mu should be chosen to produce branch lengths compatible with the desired level of sequence divergence in downstream analysis.

  6. Substitution Model(--substitution_model):

    Substitution model to use for simulation. Options are JC69, K2P, K3P, JC, GTR.

    Default is GTR.

  7. Substitution Model Parameters(--model_parameters):

    Substitution model parameters (comma-separated, no spaces) based on the selected substitution model.

     - JC69: No parameters required (leave empty)
    
     - K2P: Single float value
             - Transition/transversion ratio (kappa) > 0
    
     - K3P: Three comma-separated floats (all > 0)
             - Position 1: Transition rate (alpha)
             - Position 2: Transversion rate beta (A↔C, G↔T)
             - Position 3: Transversion rate gamma (A↔T, C↔G)
    
     - GTR: Six comma-separated floats for substitution rates
             - Position 1: A↔G
             - Position 2: A↔C
             - Position 3: A↔T
             - Position 4: G↔C
             - Position 5: G↔T
             - Position 6: C↔T
    

    For exampple, to set GTR model with rates A↔G=0.289, A↔C=0.081, A↔T=0.082, G↔C=0.046, G↔T=0.069, C↔T=0.433, use:

    --model_parameters 0.289,0.081,0.082,0.046,0.069,0.433

  8. Mutation Rate(--mutation_rate):

    Substitutions per site per year for simulation.

  9. Use Weighted Mutation(--use_weighted_mutation):

    Whether to use weighted mutation. Default is True.

    To disable weighted mutation, set to False.

  10. Weighted Mutation File(--weighted_mutation_file):

    Path to weighted mutation sites CSV file containing position and entropy columns.

    Generated by SILOSim Profiler as profiler/{genome}_entropy.csv.

    Required when --use_weighted_mutation is set to True.

  11. Use Gene Gain/Loss(--use_gain_loss):

    Whether to simulate gene gain/loss events.

    Default is True.To disable gene gain/loss simulation, set to False.

  12. Gain Rate(--gain_rate):

    Gene gain rate (λ) for Poisson process model of gain events.

    Required when --use_gain_loss is set to True.

  13. Loss Rate(--loss_rate):

    Gene loss rate (μ) for Poisson process model of loss events.

    Required when --use_gain_loss is set to True.

  14. Genome Bin Summary(--bin):

    Path to genome bin summary CSV file.

    Generated by SILOSim Profiler as profiler/{genome}_bin_summary.csv.

    Required when --use_gain_loss is set to True.

  15. Position Coverage(--position_coverage):

    Path to position coverage RDS file.

    Generated by SILOSim Profiler as profiler/{genome}_position_coverage.RDS.

    Required when --use_gain_loss is set to True.

  16. MGE Data(--mge_data):

    Path to mobile genetic element (MGE) data RDS file.

    Generated by SILOSim Profiler as profiler/{genome}_mges.RDS.

    Required when --use_gain_loss is set to True.

  17. MGE Fasta(--mge_fasta):

    Path to mobile genetic element (MGE) FASTA file.

    Generated by SILOSim Profiler as profiler/{genome}_mges.fasta.

    Required when --use_gain_loss is set to True.

  18. MGE Entropy(--mge_entropy):

    Path to directory containing entropy CSV files for all MGEs in the database.

    Generated by SILOSim Profiler as profiler/{genome}_mge_entropy/ directory.

    Required when --use_gain_loss is set to True.

  19. Use Recombination(--use_recombination):

    Whether to simulate recombination events.

    Default is True.To disable recombination simulation, set to False.

  20. Recombination Rate(--recombination_rate):

    Recombination rate (event per genome per year) for simulation.

    Required when --use_recombination is set to True.

  21. Minimum Recombination Size(--min_recombination_size):

    Minimum recombination size (bp) for simulation. Default is 1 bp.

    Required when --use_recombination is set to True.

  22. Mean Recombination Size(--mean_recombination_size): Mean recombination size (bp) for simulation.

    Required when --use_recombination is set to True.

  23. Nu Parameter(--nu): Nu parameter (snps / total length of recombination in bp) for recombination simulation.

    Required when --use_recombination is set to True.

Optional Input

  1. Output Directory(--output):

    Output directory for the simulation results.

    If not provided, defaults to silosim_simulator_output_<YYYY_MM_DD_HHMMSS> under the current working directory.

  2. Random Seed(--seed):

    Random seed for reproducibility.

    Default is the current date in the format of YYYYMMDD.

  3. Prefix(--prefix):

    Prefix for config file, output files, and analysis naming.

    If not provided, defaults to timestamp: YYYY_MM_DD_HHMMSS.

  4. Conda Environment Directory(--conda_prefix):

    Directory for conda environments needed for this analysis.

    If not provided, defaults to <OUTPUT>/conda_envs_<YYYY_MM_DD_HHMMSS>.

  5. Export Per-Event Genomes(--export-per-event-genomes):

    Export a genome snapshot after each individual evolutionary event is applied, rather than only after all events within an interval have been accumulated.

    Events include substitutions, recombination (if --use_recombination is True), and gene gain/loss (if --use_gain_loss is True).

    By default, genomes are exported once per interval after all events are applied. This option enables fine-grained inspection of sequential evolutionary changes but may significantly increase output size and runtime.

Output

  1. Simulated Phylogeny:

    Dated tree (Newick format): <OUTPUT>/treesim/simu_tree.nwk

  2. Simulated Genomes FASTA file: <OUTPUT>/simulator/output/simulated_genomes.fasta

    When --export-per-event-genomes is enabled, two additional files are exported:

    • <OUTPUT>/simulator/output/per_event_genomes.fasta: Contains only the per-event genome snapshots, each capturing the genome state immediately after an individual evolutionary event is applied. Headers follow the format >{node}_{event_type}_{ID} (e.g., >N3_mutation_1, >N5_gain_1, >N3_recomb_int_3).

    • <OUTPUT>/simulator/output/combined_simulated_genomes.fasta: Concatenation of simulated_genomes.fasta and per_event_genomes.fasta.

  3. Event Tables:

  • Mutation Events:

    • File: <OUTPUT>/simulator/output/mutations.csv

    • Columns in the mutations.csv file:

      • node: Name of the node or taxa receiving the mutation

      • alignment_position: Position in the concatenated genome where the mutation occurred.

      • profile_position: Position in the profile of the concatenated genome where the mutation occurred. The value of profile_position might be different from alignment_position since the gain/loss events alter the genome size.

      • original_base: Reference base at the position before mutation

      • mutated_base: Altered base at the position after mutation

  • Gain Events:

    • File: <OUTPUT>/simulator/output/gain.csv
    • Columns in the gain.csv file:
      • node: Name of the node or taxa receiving the MGE
      • insert_bin: ID of the bin that gained MGE
      • insertion_profile_position: Position in the profile of the concatenated genome where the MGE was inserted
      • insert_mge_id: ID of the MGE that was gained
      • insert_mge_length: Length of the MGE that was gained
  • Loss Events:

    • File: <OUTPUT>/simulator/output/loss.csv
    • Columns in the loss.csv file:
      • node: Name of the node or taxa losing the MGE
      • loss_bin: ID of the bin(MGE) that was lost
      • loss_mge_id: ID of the MGE that was lost
      • loss_mge_length: Length of the MGE that was lost
  • Recombination Events:

    • File: <OUTPUT>/simulator/output/recombination.csv
    • Columns in the recombination.csv file:
      • donor: Name of the node or taxa giving the recombination

      • recipient: Name of the node or taxa receiving the recombination

      • start: Start position in the concatenated genome where the recombination occurred

      • end: End position in the concatenated genome where the recombination occurred

      • size: Size of the recombination event in base pairs

      • snp_introduced: Number of SNPs introduced by the recombination event

      • snp_positions: list of positions in the concatenated genome where SNPs were introduced by the recombination event, seperated by semicolon (;)

  • Summary Statistics:

    • File: <OUTPUT>/simulator/output/stats.txt
    • Summary of simulation statistics including:
      • rm
      • Nu
      • Total Mutations
      • Total Recombination Events
      • Total Gain Events
      • Total Loss Events
  • SNP Matrix:

    • Text Matrix: <OUTPUT>/simulator/output/snp_matrix.txt

    • Visualization of SNP Matrix: <OUTPUT>/simulator/output/snp_matrix_visualization.pdf

  • Visualization:

    • File: <OUTPUT>/simulator/output/simulation_visualization.png

    • Generated when both gain/loss simulation and weighted mutation are enabled

    • Panels (top to bottom):

      1. Site coverage: Coverage across alignment positions identifying MGE regions
      2. Gain/loss events: Blue indicates gains, red indicates losses
      3. Mutation probability: Site-specific mutation probabilities from Genome Profiler
      4. Simulated mutations: Simulated mutations across simulated genomes

Reference

This tool is adapted from CoreSimul (https://github.com/lbobay/CoreSimul, DOI: 10.1186/s12859-020-03619-x) Please cite the original publication when using this tool.