seq2squiggle 0.3.2

End-to-end simulation of nanopore sequencing signals with feed-forward transformers

seq2squiggle

seq2squiggle is a deep learning-based tool for generating artifical nanopore signals from DNA sequence data.

Please cite the following publication if you use seq2squiggle in your work:

• Beslic D, Kucklick M, Engelmann S, Fuchs S, Renard BY, Körber N. End-to-end simulation of nanopore sequencing signals with feed-forward transformers. Bioinformatics. 2024; btae744. doi:10.1093/bioinformatics/btae744

Installation

Dependencies

seq2squiggle requires Python >= 3.10.

We recommend to run seq2squiggle in a separate conda / mamba environment. This keeps the tool and its dependencies isolated from your other Python environments.

conda create -n seq2squiggle-env python=3.10
conda activate seq2squiggle-env 

Install with pip

pip install seq2squiggle 

Install with conda or mamba

mamba install seq2squiggle -c bioconda

Install from source

git clone https://github.com/ZKI-PH-ImageAnalysis/seq2squiggle.git
cd seq2squiggle
pip install . 

Download training data and model weights

seq2squiggle requires compatible pretrained model weights to make predictions, which can be specified using the --model command-line parameter.

If you do not provide a model file, seq2squiggle will automatically attempt to download a compatible model file to ensure predictions can be made.

Predict signals from FASTA/Q file

seq2squiggle simulates artificial signals based on an input FASTA or FASTQ file. It supports two main simulation modes.

Reference mode (default)

Input: FASTA file (e.g. genome or transcriptome reference)

Goal: Generate synthetic reads and simulate corresponding signals.

How it works

• Synthetic reads are sampled from the reference sequence.
• Each read is used to simulate a raw signal.
• Use -n, -c, or -r to control number, coverage, or average read length.
• Use --distr to specify the read length distribution.

Example:

# Generate 10,000 reads & signals from a fasta reference file
seq2squiggle predict example.fasta -o example.blow5 -n 10000

# Generate 10,000 reads & signals using R9.4.1 chemistry
seq2squiggle predict example.fasta -o example.blow5 -n 10000 --profile dna_r9_min

# Export as pod5
seq2squiggle predict example.fastq -o example.pod5 -n 10000

# Generate reads with a coverage of 30 and an average read length of 5,000:
seq2squiggle predict example.fasta -o example.blow5 -c 30 -r 5000

Read mode (--read-input)

Input: FASTA or FASTQ file containing pre-existing reads (e.g. basecalled reads)

Goal: Simulates signals directly from input reads - no additional fragmentation

There are two behaviors within Read Mode:

One signal per input read

seq2squiggle predict example.fastq -o example.blow5 --read-input

This generates exactly one signal per input read (1:1 mapping).

Multiple signals based on input reads (-n)

seq2squiggle predict example.fastq -o example.blow5 --read-input -n 1000

This samples n total signals by selecting input reads without splitting them. It's not simulating n signals per read, but rather simulating n signals across the input reads.

This hybrid approach is useful when:

• You want to simulate larger datasets from real reads without generating synthetic sequences.
• You want to preserve real read characteristics while exploring noise or signal variability.

Notes

Output format

By default, the output is in SLOW5/BLOW5 format. Exporting to the new POD5 format is also supported, though BLOW5 is preferred for its stability. You will need to specify the path to the model through the configuration file.

Performance tips

For optimal performance, running seq2squiggle on a GPU is recommended, especially to speed up inference for a high number of samples. However, the tool also works on CPU-only systems, though at a slower inference speed.

Noise options

Signal generation can be customized with:

• Noise Sampler: predicts per-k-mer signal amplitude variation
• Duration Sampler: predicts per-k-mer dwell time (event length) variation

You can disable either or both to use static distributions instead.

Examples using different noise options

# Default configuration (noise sampler and duration sampler enabled)
seq2squiggle predict example.fasta -o example.blow5

# Using the noise sampler with increased noise standard deviation and the duration sampler
seq2squiggle predict example.fasta -o example.blow5 --noise-std 1.5

# Using a static normal distribution for the amplitude noise and the duration sampler
seq2squiggle predict example.fasta -o example.blow5 --noise-std 1.0 --noise-sampling False

# Using the noise sampler and a static normal distribution for event durations
seq2squiggle predict example.fasta -o example.blow5 --duration-sampling False --dwell-std 4.0

# Using the noise sampler with ideal event lengths (each k-mer event will have a length of 10):
seq2squiggle predict example.fasta -o example.blow5 --duration-sampling False --dwell-mean 10.0 --dwell-std 0.0

# Using a static normal distribution for amplitude noise and ideal event lengths:
seq2squiggle predict example.fasta -o example.blow5 --duration-sampling False --dwell-mean 10.0 --dwell-std 0.0 --noise-sampling False --noise-std 1.0

# Generating reads with no amplitude noise and ideal event lengths:
seq2squiggle predict example.fasta -o example.blow5 --duration-sampling False --dwell-mean 10.0 --dwell-std 0.0 --noise-sampling False --noise-std 0.0 

Train a new model

seq2squiggle can be trained using segmented sequence and signal data obtained from a signal-to-reference alignment. These alignments can be generated using uncalled4 or f5c, which produce an eventalign table mapping raw signals to their corresponding reference sequences.

1. Generate training data with uncalled4 or f5c

You need event-level signal alignments as training data. These can be generated using either uncalled4 or f5c.

Using uncalled4

Run the following command to align reads and extract event-level signals:

uncalled4 align input.ref input.slow5 \
    --bam-in input.bam \
    --eventalign-out events.tsv \
    --eventalign-flags print-read-names,signal-index,samples \
    --pore-model dna_r10.4.1_400bps_9mer \
    --flowcell FLO-MIN114 \
    --kit SQK-LSK114  

Using f5c

Alternatively, you can use f5c to perform the alignment:

f5c eventalign \
    -r input.fastq \
    -b input.bam \
    -g input.ref \
    --slow5 input.slow5 \
    -o events.tsv \
    --samples \
    --signal-index \
    --print-read-names \
    --collapse-events \
    --pore r10

2. Filter and transform into pA values (optional)

Since uncalled4 version 4.1.0, it no longer outputs scaled pA values directly. Instead, it outputs the normalized pore model values, which have a mean of 0 and a standard deviation of 1. seq2squiggle requires the signals to be in pA, so they must be converted back.

To denormalize signals, you first need to compute the average mean and average standard deviation of the pA values from your data. This can be done using sigtk:

sigtk pa -n input.slow5 | \
    cut -f3 | \
    sed 's/,/\n/g' | \  
    datamash mean 1 sstdev 1 | \
    awk '{print "Mean pa_mean:", $1, "Mean pa_std:", $2}' > sigtk.txt

The sigtk.txt will look like this

Mean pa_mean: 80.646167196571
Mean pa_std: 17.826207455846

Once you have the sigtk.txt file, you can run the standardize-events.py script to denormalize the signal values:

./src/seq2squiggle/standardize-events.py events.tsv events_norm.tsv --sigtk sigtk.txt

3. Preprocess the event table to .npy chunks

Before training, the eventalign table needs to be preprocessed into .npy chunks. It is recommended to first split the data into training and validation sets. This can be done based on chromosome regions (e.g., using chromosome 1 for training and chromosome 21 for validation) or by applying a random split with a specified fraction.

Once the data is split, preprocess each subset:

seq2squiggle preprocess events_norm_train.tsv train_dir --max-chunks -1 --config my_config.yml
seq2squiggle preprocess events_norm_valid.tsv valid_dir --max-chunks -1 --config my_config.yml

4. Train

Once preprocessing is complete, you can start training:

seq2squiggle train train_dir valid_dir --config my_config.yml --model last.ckpt

For a detailed explanation of the configuration parameters, see CONFIG.md

Acknowledgement

The model is based on xcmyz's implementation of FastSpeech. Some code snippets for preprocessing DNA-signal chunks have been taken from bonito. We also incorporated code snippets from Casanovo for different functionalities, including downloading weights, logging, and the design of the main function. Additionally, we used parameter profiles from squigulator for various chemistries to set digitisation, sample-rate, range, median_before, and other signal parameters. These profiles are detailed in squigulator's documentation.