Metadata-Version: 2.4
Name: neat-genreads
Version: 4.6.1
Summary: NGS Simulation toolkit
License: BSD 3-Clause License
License-File: LICENSE.md
Author: Joshua Allen
Author-email: jallen17@illinois.edu
Requires-Python: >=3.10,<4.0
Classifier: License :: Other/Proprietary License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: 3.14
Requires-Dist: PyYAML (>=6.0,<7.0)
Requires-Dist: biopython (==1.85)
Requires-Dist: frozendict (>=2.4.6,<3.0.0)
Requires-Dist: matplotlib (>=3.10,<4.0)
Requires-Dist: numpy (>=2.2,<3.0)
Requires-Dist: pip (>=25.2,<26.0)
Requires-Dist: pkginfo (>=1.8.3,<2.0.0)
Requires-Dist: pysam (>=0.23,<0.24)
Requires-Dist: scipy (>=1.14.1,<2.0.0)
Description-Content-Type: text/markdown

# The NEAT Project v4.5

Welcome to the NEAT project, the NExt-generation sequencing Analysis Toolkit, version 4.5. See the [ChangeLog](ChangeLog.md) for the full version history.

## NEAT v4.5

NEAT 4.5 is the current release. It adds the `neat compare-vcfs` subcommand, removes the deprecated `cleanup_splits`/`reuse_splits` config keys, and ships a number of performance and correctness fixes accumulated across the 4.4.x line. See the [ChangeLog](ChangeLog.md) for details.

We consider NEAT 4.x to be a stable release line and will continue to provide bug fixes and support. We will consider new features and pull requests. Please include justification for major changes. See [contribute](CONTRIBUTING.md) for more information. If you'd like to use some of our code in your own, please review the [license](LICENSE.md) first.

We've deprecated NEAT's command-line interface options for the most part, opting to simplify things with configuration files. If you require the CLI for legacy purposes, NEAT 3.4 was our last release to be fully supported via command-line interface. Please convert your CLI commands to the corresponding configuration file for future runs.

### Statement of Need

Developing and validating bioinformatics pipelines depends on access to genomic data with known ground truth. As a result, many research groups rely on simulated reads, and it can be useful to vary the parameters of the sequencing process itself. NEAT addresses this need as an open-source Python package that can integrate seamlessly with existing bioinformatics workflows—its simulations account for a wide range of sequencing parameters (e.g., coverage, fragment length, sequencing error models, mutational frequencies, ploidy, etc.) and allow users to customize their sequencing data.

NEAT is a fine-grained read simulator that simulates real-looking data using models learned from specific datasets. It was originally designed to simulate short reads and is adaptable to different machines, with custom error models and the capability to handle single-base substitutions, indel errors, and other types of mutations. Unlike simulators that rely solely on fixed error profiles, NEAT can learn empirical mutation and sequencing models from real datasets and use these models to generate realistic sequencing data, providing outputs in several common file formats (e.g., FASTQ, BAM, and VCF). There are several supporting utilities for generating models used for simulation and for comparing the outputs of alignment and variant calling to the golden BAM and golden VCF produced by NEAT.

To cite this work, please use both of the following:

1. > Stephens, Z. D., Hudson, M. E., Mainzer, L. S., Taschuk, M., Weber, M. R., & Iyer, R. K. (2016). Simulating Next-Generation Sequencing Datasets from Empirical Mutation and Sequencing Models. _PLOS ONE_, _11_(11), e0167047. https://doi.org/10.1371/journal.pone.0167047

2. > Allen, J. M., Gandhi, K. R., Alhazmy, R., Wasnik, Y., & Fliege, C. E. (2026). Enhancing next-generation sequencing simulation: Updates to NEAT. _Journal of Open Source Software_, _11_(121), 9056. https://doi.org/10.21105/joss.09056

[![DOI](https://joss.theoj.org/papers/10.21105/joss.09056/status.svg)](https://doi.org/10.21105/joss.09056)

## Table of Contents

* [The NEAT Project v4.5](#the-neat-project-v45)
* [NEAT v4.5](#neat-v45)
* [Table of Contents](#table-of-contents)
  * [Prerequisites](#prerequisites)
  * [Installation](#installation)
  * [Usage](#usage)
  * [Functionality](#functionality)
    * [Estimated runtimes](#estimated-runtimes)
  * [Examples](#examples)
    * [Whole genome simulation](#whole-genome-simulation)
    * [Targeted region simulation](#targeted-region-simulation)
    * [Insert specific variants](#insert-specific-variants)
    * [Single-end reads](#single-end-reads)
    * [Large single-end reads](#large-single-end-reads)
    * [Parallelizing simulation](#parallelizing-simulation)
  * [Utilities](#utilities)
    * [`neat model-fraglen`](#neat-model-fraglen)
    * [`neat gen-mut-model`](#neat-gen-mut-model)
    * [`neat model-seq-err`](#neat-model-seq-err)
    * [`neat model-qual-score`](#neat-model-qual-score)
    * [`neat model-gc-bias`](#neat-model-gc-bias)
    * [`neat compare-vcfs`](#neat-compare-vcfs)
    * [`neat bacterial-wrapper`](#neat-bacterial-wrapper)
  * [Validating NEAT Outputs](#validating-neat-outputs)
  * [Tests](#tests)
    * [Guide to run locally](#guide-to-run-locally)
    * [Note on Sensitive Patient Data](#note-on-sensitive-patient-data)

## Prerequisites

The most up-to-date requirements are found in `environment.yml`.

* Some version of Anaconda to set up the environment
* `python` >= 3.10
* `poetry` >= 2.0
* `biopython` == 1.85
* `pkginfo`
* `matplotlib`
* `numpy` >= 2.2
* `pyyaml`
* `scipy`
* `pysam` >= 0.23
* `frozendict`

NEAT assumes a Linux environment and was tested primarily on Ubuntu Linux. It should work on most Linux systems. If you 
use another operating system, please install WSL or a similar tool to create a Linux environment to operate NEAT from.
For setting up NEAT, you will need `Anaconda` (or `miniconda`). The method described here installs NEAT as a base package 
in the active `conda` environment, so whenever you want to run NEAT, you can first activate the environment, then run
from any place on your system. If you desire VCF files, please install also `bcftools`. For your convenience, we have
added `bcftools` to the environment file, as it is available from `conda`. You may remove this line if you do not want or 
need VCF files with the variants NEAT added.

## Installation
To install NEAT, you must create a virtual environment using a tool such as `conda`. 

First, clone the environment and move to the NEAT directory:

```bash
$ git clone https://github.com/ncsa/NEAT.git
$ cd NEAT
```

A quick form of installation uses `bioconda`. You must run these commands inside the NEAT project directory.

```bash
(base) $ conda create -n neat -c conda-forge -c bioconda neat
(base) $ conda activate neat
(neat) $ neat --help # tests that NEAT has installed correctly
```

Alternatively, instead of the `bioconda` method, you can use the `poetry` module in build a wheel file, which can then be `pip` installed.

Once `conda` is installed, the following command can be run for easy setup. 
In the NEAT repository, at the base level is the `environment.yml` file you will need. Change directories into the NEAT repository
and run:

```bash
(base) $ conda env create -f environment.yml
(base) $ conda activate neat
(neat) $ poetry install
(neat) $ neat --help # tests that NEAT has installed correctly
```

Assuming you have installed `conda`, run `source activate` or `conda activate`.

NEAT is tested on Linux. macOS may work but is not regularly tested; on Windows, use WSL.

Alternatively, if you wish to work with NEAT in the development-only environment, you can use `poetry install` within
the NEAT repo, after creating the `conda` environment:

```bash
$ conda env create -f environment.yml
$ conda activate neat
$ poetry install
```

Any updates to the NEAT code will be instantly reflected with a `poetry install` version.

Notes: If any packages are struggling to resolve, check the channels and try to manually `pip` install the package to see 
if that helps (but note that NEAT is not tested on the `pip` versions). If `poetry` hangs for you, try the following fix 
(from https://github.com/python-poetry/poetry/issues/8623):

```bash
export PYTHON_KEYRING_BACKEND=keyring.backends.null.Keyring
```
Then, re-run `poetry install`.

Test your install by running:
```bash
$ neat --help
```

You can also try running it using the Python command directly:
```bash
$ python -m neat --help
```

## Usage
NEAT's core functionality is invoked using the `read-simulator` command. Here's the simplest invocation of read-simulator using default parameters. This command produces a single-ended FASTQ file with reads of length 151, ploidy 4, coverage 15x, default sequencing substitution, and default mutation rate models.

Contents of `neat_config.yml`:
```yml
reference: /path/to/my/genome.fa
read_len: 151
ploidy: 4
coverage: 15
```

```bash
neat read-simulator -c neat_config.yml -o simulated_data
```

The `--output` (`-o`) option sets the folder to place output data. If the folder does not exist, Python will attempt to create it. To specify a common filename prefix, you can additionally add the `--prefix` (`-p`), e.g., `-p ecoli_20x` will result in output files `ecoli_20x.fastq.gz`, `ecoil_20x.vcf.gz`, etc.

A config file is required. The config is a `.yml` file specifying the input parameters. The following is a brief
description of the potential inputs in the config file. See `NEAT/config_template/template_neat_config.yml` for a template config file to copy and use for your runs. It provides a detailed description of all the parameters that NEAT uses.

To run the simulator in multithreaded mode, set the `threads` value in the config to something greater than 1.

`reference`: Full path to a FASTA file to generate reads from.  

`read_len`: The length of the reads for the FASTQ (if using). _Integer value, default 151._

`coverage`: Desired average coverage depth. _Numeric value (fractional allowed, e.g. `0.5` for low-pass simulation), must be > 0, default = 10._

`ploidy`: Desired value for ploidy (# of copies of each chromosome in the organism, where if ploidy > 2, "heterozygous"  mutates floor(ploidy / 2) chromosomes). _Default is 2._

`paired_ended`: If paired-ended reads are desired, set this to `True`. Setting this to `True` requires either entering values for `fragment_mean` and `fragment_st_dev` or entering the path to a valid `fragment_model`.    

`fragment_mean`: Use with paired-ended reads, setting a fragment length mean manually.    

`fragment_st_dev`: Use with paired-ended reads, setting the standard deviation of the fragment length dataset.

The following values can be set to `True` or omitted to use defaults. If `True`, NEAT will produce the file type.

The default is given:

`produce_bam`: `False`    
`produce_vcf`: `False`    
`produce_fastq`: `True`    

More parameters are below:

| Parameter           | Description                                                                                                                                                                                                   |
|---------------------|---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| `error_model`       | Full path to an error model or quality score model generated by NEAT. Leave empty to use default model (default model based on human, sequenced by Illumina).                                                 |
| `mutation_model`    | Full path to a mutation model generated by NEAT. Leave empty to use a default model (default model based on human data sequenced by Illumina).                                                                |
| `fragment_model`    | Full path to fragment length model generated by NEAT. Leave empty to use default model (default model based on human data sequenced by Illumina).                                                             |
| `gc_model`          | Full path to GC-bias model generated by NEAT. Leave empty for no GC bias.                                                                                                                                     |
| `avg_seq_error`     | Average sequencing error rate for the sequencing machine. Use to increase or decrease the rate of errors in the reads. Float between 0 and 0.3. Default is set by the error model.                            |
| `rescale_qualities` | Rescale the quality scores to reflect the `avg_seq_error` rate above. Set `True` to activate if you notice issues with the sequencing error rates in your dataset.                                            |
| `quality_offset`    | ASCII offset for quality score encoding. Default `33` (Sanger/Illumina 1.8+ Phred+33). Only change this if your data uses a non-standard encoding (valid range 33–64).                                       |
| `include_vcf`       | Full path to list of variants in VCF format to include in the simulation. These will be inserted as they appear in the input VCF into the final VCF, and the corresponding FASTQ and BAM files, if requested. |
| `target_bed`        | Full path to list of regions in BED format to target. All areas outside these regions will have coverage of 0.                                                                                                |
| `discard_bed`       | Full path to a list of regions to discard, in BED format.                                                                                                                                                     |
| `mutation_rate`     | Desired rate of mutation for the dataset. Float between 0.0 and 0.3 (default is determined by the mutation model).                                                                                            |
| `mutation_bed`      | Full path to a list of regions with a column describing the mutation rate of that region, as a float with values between 0 and 0.3. The mutation rate must be in the third column as, e.g., `mut_rate`=0.00.  |
| `rng_seed`          | Manually enter a seed for the random number generator. Used for repeating runs. Must be an integer.                                                                                                           |
| `min_mutations`     | Set the minimum number of mutations that NEAT should add, per contig. Default is 0. We recommend setting this to at least one for small chromosomes, so NEAT will produce at least one mutation per contig.   |
| `overwrite_output`  | If `true`, existing output files with the same name will be overwritten. Default `false` (NEAT exits with an error if output files already exist).                                                            |
| `no_coverage_bias`  | If `true`, disables statistical coverage models and forces uniform coverage across the genome. Intended for debugging. Default `false`.                                                                        |
| `threads`           | Number of threads to use. More than 1 will use multi-threading to speed up processing. With `threads > 1`, NEAT splits each contig into chunks; with `threads == 1`, one chunk per contig is used.            |
| `parallel_block_size` | Per-chunk size in bases when `threads > 1`. Default `0` (auto-tune from total genome length and thread count, targeting ~8 chunks per thread). Set to a positive integer to override. Ignored when `threads == 1`. |

The command line options for NEAT are as follows:

Universal options can be applied to any subfunction. The commands should come before the function name (e.g., neat --log-level DEBUG read-simulator ...), except -h or --help, which can appear anywhere in the command.
| Universal Options   | Description                          |
|---------------------|--------------------------------------|
| -h, --help          | Displays usage information           |
| --no-log            | Turn off log file creation           |
| --log-name LOG_NAME | Custom name for log file, can be a full path (default is current working directory with a name starting with a timestamp)|
| --log-level VALUE   | VALUE must be one of [DEBUG, INFO, WARN, WARNING, ERROR] - sets level of log to display |
| --log-detail VALUE   | VALUE must be one of [LOW, MEDIUM, HIGH] - how much info to write for each log record |
| --silent-mode       | Writes logs, but suppresses stdout messages |

`read-simulator` command line options
| Option              | Description                         |
|---------------------|-------------------------------------|
| -c VALUE, --config VALUE | The VALUE should be the name of the config file to use for this run |
| -o OUTPUT_DIR, --output_dir OUTPUT_DIR | The path to the directory to write the output files |
| -p PREFIX, --prefix String | The prefix for file names |

## Functionality

![Diagram describing the way that genReads simulates datasets](docs/NEATNEAT.png "Diagram describing the way that gen_reads simulates datasets")

NEAT produces simulated sequencing datasets. It creates FASTQ files with reads sampled from a provided reference genome, using sequencing error rates and mutation rates learned from real sequencing data. The strength of NEAT lies in the ability for the user to customize many sequencing parameters, produce 'golden,' true-positive datasets. We are working on expanding the functionality even further to model more species, generate larger variants, model tumor/normal data, and more!

Features:

- Simulate single-end and paired-end reads 
- Custom read length
- Can introduce random mutations and/or mutations from a VCF file
  - Supported mutation types include SNPs, indels (of any length), inversions, translocations, duplications
  - Can emulate multi-ploid heterozygosity for SNPs and small indels
- Can simulate targeted sequencing via BED input specifying regions to sample from
- Can accurately simulate large, single-end reads with high indel error rates (PacBio-like) given a model
- Specify simple fragment length model with mean and standard deviation or an empirically learned fragment distribution
- Simulates quality scores using either the default model or empirically learned quality scores using `neat gen_mut_model`
- Introduces GC-bias using an empirically learned model from `neat model-gc-bias`
- Introduces sequencing substitution errors using either the default model or empirically learned in `utilities`
- Output a VCF file with the 'golden' set of true positive variants. These can be compared to bioinformatics workflow output (includes coverage and allele balance information)
- Output a BAM file with the 'golden' set of aligned reads. These indicate where each read originated and how it should be aligned with the reference
- Create paired tumour/normal datasets using characteristics learned from real tumour data

### Output ordering

The BAM file NEAT produces is coordinate-sorted (by construction at write time — no separate sort pass is run, which used to allocate ~1 GB of sort buffer at stitch time).

The FASTQ files are written in the natural fragment-sampling order. This is *roughly* random — fragments are drawn from batched random positions — but is not a strict uniform shuffle. If your downstream tooling assumes a real-sequencer-style shuffle, pipe the FASTQ through [`seqkit shuffle`](https://bioinf.shenwei.me/seqkit/usage/#shuffle):

```sh
# Single-end
seqkit shuffle reads.fastq.gz -o reads.shuffled.fastq.gz

# Paired-end — use a shared seed so the two files stay mate-aligned
seqkit shuffle -s 42 reads_r1.fastq.gz -o reads_r1.shuffled.fastq.gz
seqkit shuffle -s 42 reads_r2.fastq.gz -o reads_r2.shuffled.fastq.gz
```

### Estimated runtimes

> **Note:** The tables below are from the original NEAT 4.4 (v4.4.0) benchmark. NEAT v4.4.4 is roughly **13× faster on multi-threaded ecoli** and uses **~3× less memory** thanks to the performance work in v4.4.3 and v4.4.4. The relative shape of the tables (size scaling, contig vs. size mode tradeoffs) remains accurate, but absolute runtimes should be divided by ~10–15× for v4.4.4+. See `ChangeLog.md` (v4.4.3 and v4.4.4 entries) for a detailed before/after table on ecoli and a c_elegans scale-test.

To give users a sense of how long `neat read-simulator` runs may take, we benchmarked NEAT 4.4 on several reference genomes. All runs were paired-end, with read length of 150 bp, coverage of 10, fragment mean of 300 bp, and fragment standard deviation of 50 bp. Runtimes are reported as the average across three unique runs (`Avg. time (ms)`) and the corresponding runtime in minutes. Cells marked with N/A indicate that NEAT was not able to run to completion.

Benchmarks were run on a System76 Meerkat with a 13th Gen Intel Core i3-1315U (8 logical cores, up to 4.50 GHz) and 16 GiB RAM, using a 512 GB SSD and Ubuntu 24.04.3 LTS (Linux kernel 6.14). Actual runtimes will vary depending on your hardware.

#### NEAT Single-Threaded (contig mode)

The inputs in single-threaded, contig-based mode most closely replicate the behavior of legacy versions of NEAT.

The configuration used:

- `threads: 1`  (NEAT processes one contig per chunk in single-thread mode)

| Organism        | File size (bytes) | Avg. runtime (ms) | Avg. runtime (min) |
|-----------------|-------------------|-------------------|--------------------|
| H1N1            | 13,599            | 702               | 0.0117             |
| Herpes virus    | 174,029           | 4,909             | 0.0818             |
| Pneumonia       | 2,137,444         | 64,424            | 1.0737             |
| *E. coli*       | 1,258,807         | 146,843           | 2.4474             |
| *S. cerevisiae* | 12,310,392        | 339,842           | 5.6640             |
| Honeybee        | 228,091,137       | N/A               | N/A                |
| Rice            | 394,543,607       | N/A               | N/A                |
| *Miscanthus*    | 2,718,242,062     | N/A               | N/A                |

For small genomes, single-threaded performance is already fast (on the order of seconds). Larger genomes could not be fully benchmarked under this configuration.

#### NEAT Multi-Threaded (size mode)

Here we enabled NEAT’s multi-threaded mode, which splits contigs into size-based blocks and processes them concurrently. The configuration used:

- `parallel_block_size: 500000`
- `produce_bam: false`
- `threads: 7`

| Organism        | File size (bytes) | Avg. runtime (ms) | Avg. runtime (min) |
|-----------------|-------------------|-------------------|--------------------|
| H1N1            | 13,599            | 199               | 0.0033             |
| Herpes virus    | 174,029           | 5,267             | 0.0878             |
| Pneumonia       | 2,137,444         | 31,113            | 0.5186             |
| *E. coli*       | 1,258,807         | 58,927            | 0.9821             |
| *S. cerevisiae* | 12,310,392        | 139,148           | 2.3191             |
| Honeybee        | 228,091,137       | 3,040,336         | 50.6723            |
| Rice            | 394,543,607       | 4,335,126         | 72.2521            |
| *Miscanthus*    | 2,718,242,062     | 24,876,744        | 414.6              |

For mid-sized genomes (e.g., *E. coli* and *S. cerevisiae*), enabling parallelization reduced runtimes by roughly a factor of two to three compared to the base configuration. For larger genomes (honeybee and rice), the parallel configuration may make multi-hour simulations feasible.

Access to high-performance computing systems may improve runtimes.

## Examples

The following commands are examples for common types of data to be generated. The simulation uses a reference genome in fasta format to generate reads of 126 bases with default 10X coverage. Outputs paired fastq files, a BAM file and a VCF file. The random variants inserted into the sequence will be present in the VCF and the reads will show their proper alignment in the BAM. Unless specified, the simulator will also insert some "sequencing error"—random variants in some reads that represents false positive results from sequencing.

### Whole genome simulation
Simulate whole genome dataset with random variants inserted according to the default model:

Contents of neat_config.yml:
```yml
reference: hg19.fa
read_len: 126
produce_bam: True
produce_vcf: True
paired_ended: True
fragment_mean: 300
fragment_st_dev: 30
```
```bash
neat read-simulator                 \
        -c neat_config.yml          \
        -o /home/me/simulated_reads/
```

### Targeted region simulation
Simulate a targeted region of a genome (e.g., an exome) with a targeted bed:

Contents of neat_config.yml:
```yml
reference: hg19.fa
read_len: 126
produce_bam: True
produce_vcf: True
paired_ended: True
fragment_mean: 300
fragment_st_dev: 30
targed_bed: hg19_exome.bed
```
```bash
neat read-simulator                 \
        -c neat_config              \
        -o /home/me/simulated_reads/
```
### Parallelizing simulation
Split the contig into blocks rather than reading by contig. Even in single-threaded mode this is likely to be faster. The chunk size auto-tunes from total genome length and thread count, targeting roughly 8 chunks per thread for good load balancing — on small bacterial genomes you get ~500 kb chunks (similar to NEAT's old hardcoded default), on human-scale genomes you get ~6 Mb chunks (a few hundred total instead of thousands). Specify `parallel_block_size` explicitly if you want to override.

Also, we demonstrate the situation where you do not want any logs produced:

`neat_config.yml`:
```yml
reference: giant_bacterium.fa
read_len: 150
produce_bam: True
paired_ended: True
fragment_mean: 350
fragment_st_dev: 50
threads: 12
# parallel_block_size omitted: auto-tuned from genome length and thread count.
# Set explicitly (e.g. parallel_block_size: 1000000) only if you have a reason.
```
Then run with the command:
```
neat read-simulator                 \
        --no-log                    \
        -c neat_config.yml          \
        -o /home/me/simulated_reads/
```

### Multi-node deployment on HPC clusters

NEAT runs on a single node using Python's `multiprocessing`. To use multiple nodes on a supercomputer, dispatch one NEAT job per contig (or contig group) as a job-array element and concatenate the outputs afterwards. Each array task gets its own node and uses all available cores on it. NEAT itself doesn't need to know about the cluster.

SLURM example with one task per human chromosome (24 tasks):

```bash
#!/bin/bash
#SBATCH --job-name=neat-array
#SBATCH --array=1-24
#SBATCH --cpus-per-task=64
#SBATCH --mem=64G
#SBATCH --time=04:00:00

CHROM=$(awk "NR==$SLURM_ARRAY_TASK_ID" chroms.txt)   # e.g. chr1, chr2, ...
samtools faidx hg38.fa "$CHROM" > "ref_${CHROM}.fa"   # one chromosome per task

cat > "config_${CHROM}.yml" <<YAML
reference: ref_${CHROM}.fa
read_len: 150
coverage: 30
paired_ended: true
fragment_mean: 400
fragment_st_dev: 50
produce_bam: true
produce_fastq: true
produce_vcf: true
threads: ${SLURM_CPUS_PER_TASK}
YAML

neat read-simulator -c "config_${CHROM}.yml" \
                    -o "out/${CHROM}/" \
                    -p "${CHROM}"
```

After the array completes, concatenate the per-chromosome outputs:

```bash
# BAMs are already coordinate-sorted within each chromosome; cat in chromosome order
# produces a globally-sorted whole-genome BAM with no further      2_golden.bam ... out/chrY/chrY_golden.bam
samtools index all.bam

# FASTQs: gzip-concat preserves a valid gzip stream
cat out/*/chr*_r1.fastq.gz > all_r1.fastq.gz
cat out/*/chr*_r2.fastq.gz > all_r2.fastq.gz

# VCFs: use bcftools concat for proper merging
bcftools concat -o all.vcf.gz -Oz out/*/chr*_golden.vcf.gz
bcftools index all.vcf.gz
```

This gives you whole-genome simulation in roughly `single-chromosome-time / nodes` wall time. For a 30× human genome with 24 nodes (each 64 cores), a per-chromosome run takes ~10 min and the array finishes in ~10 min wall, plus a few minutes for the final concat.

### Insert specific variants
Simulate a whole genome dataset with only the variants in the provided VCF file by setting `include_vcf` to the VCF path and `mutation_rate: 0`:

Contents of neat_config.yml:
```yml
reference: hg19.fa
read_len: 126
produce_bam: True
produce_vcf: True
paired_ended: True
fragment_mean: 300
fragment_st_dev: 30
include_vcf: NA12878.vcf
mutation_rate: 0
```
```bash
neat read-simulator                 \
        -c neat_config.yml          \
        -o /home/me/simulated_reads/
```

### Single-end reads
Simulate single-end reads by omitting paired-ended options:

Contents of neat_config.yml:
```yml
reference: hg18.fa
read_len: 126
produce_bam: True
produce_vcf: True
```
```bash
neat read-simulator                 \
        -c neat_config.yml          \
        -o /home/me/simulated_read/ \
        -p 126_frags
```

### Large single-end reads
Simulate PacBio-like reads by providing an error model. Note: full long-read support is planned for a separate future release and is not yet fully validated in NEAT 4.4.

Contents of neat_config.yml:
```yml
reference: hg19.fa
read_len: 5000
error_model: errorModel_pacbio.pickle.gz
avg_seq_error: 0.1
```
```bash
neat read-simulator                 \
        -c neat_config.yml          \
        -o /home/me/simulated_reads
```

## Utilities

Several scripts are distributed with `gen_reads` that are used to generate the models used for simulation.

### `neat model-fraglen`

Computes empirical fragment length distribution from sample paired-end data. NEAT uses the template length (tlen) attribute calculated from paired-ended alignments to generate summary statistics for fragment lengths, which can be input into NEAT.

```bash
    neat model-fraglen    \
        -i input.bam      \
        -o /path/to/prefix
```

and creates `fraglen.pickle.gz` model in working directory.

### `neat gen-mut-model`

Takes reference genome and VCF file to generate mutation models:

```bash
neat gen-mut-model reference.fa input_variants.vcf   \
        -o /home/me/models
```

Trinucleotides are identified in the reference genome and the variant file. The mutation model uses trinucleotide context and selects mutation sites and alternate alleles with a transition matrix. Frequencies of each trinucleotide transition are calculated and output as a pickle file. Mutations are simulated to reflect the same context-dependent biases as the training data. As of NEAT 4.4, we have only made minor optimizations to improve the speed, and the underlying statistical models are similar to those described in the original NEAT manuscript.

| Option          | Description                                                                   |
|-----------------|-------------------------------------------------------------------------------|
| -o <str>        | Path to output file and prefix; defaults to "neat_sim" in current working dir |
| --bed           | Flag that indicates you are using a bed-restricted VCF and FASTA (see below)  |
| --save-trinuc   | Save trinucleotide counts for reference                                       |
| --human-sample  | Use to skip unnumbered scaffolds in human references                          |
| --skip-common   | Do not save common snps or high mutation areas                                |

### `neat model-seq-err`

Generates sequencing error model for NEAT.

The sequencing error model uses real FASTQ/BAM data, summarizing position-specific substitution, insertion, and deletion frequencies. These are stored as discrete transition matrices that are sampled during read generation to reproduce empirical error profiles. In this way, it is similar to a simple version of `gen-mut-model`.

This script needs revision to improve the quality-score model and to include code to learn sequencing errors from pileup data.

Note that `model-seq-err` does not allow for SAM inputs. If you would like to use
a BAM/SAM file, please use `samtools` to convert to a FASTQ, then use the FASTQ as input.

Note additionally that the file must either be unzipped or bgzipped. If your file is currently gzipped, you can
use `samtools`' built-in `bgzip` utility to perform this conversion.

```bash
gzip -d my_file.fastq.gz
bgzip my_file.fastq
```

The blocked-zip format allows for indexing of the file.

For quality scores, note that using a single number will check quality scores for every number. As this could 
potentially slow down model creation, binned quality scores are advisable.

Soon, we will take a `samtools mpileup` output as input and have some plotting features.

```bash
neat model-seq-err                                    \
        -i input_read.fq (.gz)                        \
        -o /path/to/prefix                            \
        -q quality score offset [33]                  \
        -Q maximum quality score [2, 11, 24, 37]      \
        -m maximum number of reads to process [all]   \
```

Please note that `-i2` can be used in place of `-i` to produce paired data.

### `neat model-qual-score`

Typical usage:

```bash
neat model-qual-score \
    -i input_reads.fastq(.gz)            \
    -q 33                                \
    -Q 42                                \
    -m 1000000                           \
    --markov                             \
    -o /path/to/models                   \
    -p my_qual_model
```

Similarly, use `-i2` to produce a model for paired-ended data. `-q` denotes the quality score offset, while `-Q` is the maximum quality score.

`-m` denotes the maximum number of reads to process. Use a large number or input -1 to use all reads. `--markov` fits a quality model from the input data using a Markov chain process instead of the baseline quality score model (optional).

Finally, `-o` is the output directory for the model file and `-p` is the prefix for the output model, such that the file will be written as `<prefix>.p.gz` inside the output folder.

#### Binned quality scoring for modern Illumina instruments

Modern Illumina instruments (NovaSeq 6000, NovaSeq X, NextSeq 2000) compress
Phred quality scores into a small discrete set of bins rather than emitting a
continuous range. To train a model that faithfully reproduces this behaviour,
use either the named `--quality-preset` flag or an explicit bin list via `-Q`.

**Named presets** (recommended):

```bash
# NovaSeq 6000 / NovaSeq X — Q2, Q12, Q23, Q37
neat model-qual-score -i reads.fastq.gz --quality-preset novaseq \
    -o /path/to/models -p novaseq_model

# NextSeq 2000 — Q2, Q12, Q26, Q37
neat model-qual-score -i reads.fastq.gz --quality-preset nextseq2000 \
    -o /path/to/models -p nextseq2000_model

# NextSeq 500 / MiniSeq — Q2, Q12, Q23, Q27, Q37
neat model-qual-score -i reads.fastq.gz --quality-preset nextseq500 \
    -o /path/to/models -p nextseq500_model
```

`--quality-preset` implies `--markov`; you do not need to pass both.

**Explicit bin list** — if your instrument uses non-standard bins, pass them
directly with `-Q`:

```bash
neat model-qual-score -i reads.fastq.gz -Q 2 12 23 37 --markov \
    -o /path/to/models -p custom_binned_model
```

When bins are specified, both the Markov and traditional quality models will
constrain simulation output to those exact Phred values.

### `neat model-gc-bias`

Computes GC-bias model from a BAM file and reference genome. It calculates the relative weight of fragments based on their GC content.

```bash
neat model-gc-bias    \
    -i input.bam      \
    -r reference.fa   \
    -o /path/to/prefix
```

and creates `gc_model.pickle.gz` model in the working directory.

### `neat compare-vcfs`

Compares a downstream variant caller's VCF against a NEAT-simulated truth VCF,
attributing each false negative to the simulator's own configuration. The
intended workflow is: run NEAT to produce a synthetic FASTQ and golden VCF →
run your aligner + caller on the FASTQ → compare the caller's VCF against the
golden with `neat compare-vcfs`.

Variant equivalence (multi-allelic normalization, haplotype-level matching) is
delegated to `hap.py`. The NEAT-specific value-add is the false-negative
attribution: each FN is tagged with one or more reasons drawn from the run
config recorded in `simulation_summary.json`.

```bash
neat compare-vcfs golden.vcf called.vcf            \
        --neat-run-dir /path/to/neat/output/dir    \
        --output-dir /path/to/comparison/output    \
        --reference reference.fa                   \
        [--target-bed target.bed]                  \
        [--happy-bin /abs/path/to/hap.py]          \
        [--plot]                                   \
        [--chrom-aliases aliases.tsv]
```

Outputs (in `--output-dir`):

| File | Purpose |
|------|---------|
| `comparison_summary.json`   | Machine-readable: counts (TP/FN/FP), precision/recall/F1, per-reason FN counts |
| `comparison_summary.txt`    | Human-readable rollup of the same |
| `FN_with_reasons.vcf`       | hap.py's false-negative records, each annotated with a `NEAT_REASON` INFO tag |
| `happy.vcf.gz` (+ siblings) | Raw hap.py output preserved for inspection |
| `fn_attribution.png`        | Optional — only written when `--plot` is set |

**Chromosome-name conventions:** NEAT does NOT silently normalize chrom names
across the reference, golden VCF, called VCF, and BED files. If your
`mutation_bed` uses `1`, `2`, … but the reference uses `chr1`, `chr2`, …,
`compare-vcfs` detects the mismatch up front and writes a warning into
`comparison_summary.json` suggesting the rename. Pass `--chrom-aliases
aliases.tsv` (two-column TSV: `bed_name<TAB>canonical_name`) to apply the
rename at load time. This also handles common mitochondrial variants
(`M`/`MT`/`chrM`/`chrMT`).

**False-negative reason categories:**

| Tag | Meaning |
|-----|---------|
| `outside_simulated_contigs` | FN's chromosome wasn't part of the NEAT run at all |
| `outside_mutation_bed`      | `mutation_bed` was set in the config and the FN falls outside its regions |
| `outside_target_bed`        | `target_bed` was set and the FN falls outside its regions |
| `unknown`                   | None of the above — NEAT has no specific explanation |

**`simulation_summary.json` prerequisite:** Every `neat read-simulator` run now
emits `simulation_summary.json` into its output dir as part of the standard
artifacts; `compare-vcfs` reads this file from `--neat-run-dir` to drive
attribution. No extra step required when running NEAT yourself; if you're
working with pre-NEAT-4.5 output, re-run the simulator to produce one.

**`hap.py` install:** `hap.py` is an external dependency. The cleanest install
path is a dedicated conda env:

```bash
conda create -n hap_py_env -c bioconda -c conda-forge hap.py -y
```

then point `--happy-bin` at `/path/to/hap_py_env/bin/hap.py` (or put that
directory on `$PATH`). Without `hap.py` available, `neat compare-vcfs` exits
with a clear install hint.

### `neat bacterial-wrapper`

Runs NEAT's read simulator twice on the same reference — once with a standard
mutation model ("Regular") and once with a higher mutation rate model
("Wrapped") representing a bacterium under selection pressure — then stitches
the two output sets together for downstream comparison.

```bash
neat bacterial-wrapper reference.fa bacteria_name \
        -c config.yml                             \
        -o /path/to/output/dir
```

Outputs are written into `<output_dir>/Regular/` and `<output_dir>/Wrapped/`
subdirectories, then stitched into the parent output directory. The config file
follows the same format as `neat read-simulator`.

## Validating NEAT Outputs

NEAT does not ship its own validation utilities. Use the standard bioinformatics
tools below, which handle gzipped inputs, are actively maintained, and are
typically already present in any conda/bioconda environment.

### FASTQ validation

**FastQC** — format validation plus quality metrics (GC content, adapter
contamination, per-base quality scores):

```bash
fastq read1.fastq.gz            # single-end
fastqc read1.fastq.gz read2.fastq.gz   # paired-end
```

**fastp** — lightweight format check without trimming or filtering:

```bash
fastp --in1 read1.fastq.gz --in2 read2.fastq.gz \
      --disable_adapter_trimming --disable_quality_filtering \
      --disable_length_filtering --thread 4
```

### BAM validation

**samtools quickcheck** — fast EOF and header sanity check; exits non-zero on
failure, making it CI-friendly:

```bash
samtools quickcheck output.bam
```

**samtools flagstat** — alignment statistics; a non-zero exit or obviously wrong
counts (e.g., 0 mapped reads when coverage > 0) indicate a problem:

```bash
samtools flagstat output.bam
```

**Picard ValidateSamFile** — comprehensive structural validation (CIGAR
consistency, mate-pair pairing, flag conflicts):

```bash
picard ValidateSamFile -I output.bam -MODE SUMMARY
```

### VCF validation

**bcftools stats** — reports site counts, ts/tv ratio, and indel length
distribution; useful for a sanity check against expected mutation rates:

```bash
bcftools stats output_golden.vcf.gz | grep "^SN"
```

## Tests

We provide unit tests (e.g., mutation and sequencing error models) and basic integration tests for the CLI.

### Guide to run locally
```bash
conda env create -f environment.yml
conda activate neat
poetry install --with dev
pytest -q tests
```

Please see `CONTRIBUTING.md` for more information and further instructions.

### Note on Sensitive Patient Data
ICGC's "Access Controlled Data" documentation can be found at <a href = https://docs.icgc.org/portal/access/ target="_blank">https://docs.icgc.org/portal/access/</a>. To have access to controlled germline data, a DACO must be submitted. Open tier data can be obtained without a DACO, but germline alleles that do not match the reference genome are masked and replaced with the reference allele. Controlled data includes unmasked germline alleles.

### rneat
We would like to introduce a sister tool, `rneat`, a Rust implementation of NEAT based largely on the original NEAT 2.1 codebase. It shares NEAT's core feature set while adding a native cancer-simulation workflow — tumor/normal mixtures with configurable purity, an origin-tagged truth VCF, complex structural variants (CNV/BND/INV), and per-tissue cancer models. It also keeps a low, flat memory footprint, produces byte-identical output for a given seed regardless of thread count, and ships as a single self-contained binary (no Python environment to manage). Check our repo here: <a href=https://github.com/ncsa/rusty-neat/ target="_blank">https://github.com/ncsa/rusty-neat</a>, where you can navigate to releases to download the latest binary, or find it on <a href=https://zenodo.org/records/20577688>Zenodo</a>. Check out `rneat` and let us know what you think!

