shUMI: simulating high-throughput UMI-based single-genome sequencing datasets

shUMI simulates high-throughput UMI-tagged molecule amplification and sequencing experiments from source genomes, through molecular tagging and amplification, to raw sequencing reads.

The current implementation models a staged assay workflow:

1. Input molecule population generation or loading
2. Single- or double-UMI tagging, with transcription errors, indels, and recombination
3. PCR amplification with stocastic copying (efficiency), substitutions, indels, recombination
4. Sequencing by PacBio CCS, ONT, or a lightweight generic read simulator

The recommended entry point is a YAML configuration file with the shumi full command. Each stage can also be run directly as a subcommand for modular workflows.

Installation

For a conda-based install with compatible compiled dependencies, the following installation commands will configure a new environment and install shUMI:

conda create -n shumi-env \
  -c conda-forge -c bioconda \
  --strict-channel-priority \
  python=3.10 pip \
  biopython h5py matplotlib numpy scipy pysam pyyaml tqdm \
  "samtools>=1.16" seqkit pbsim3 pbccs

conda activate shumi-env
git clone git@github.com:niaid/shumi.git
cd shumi
python -m pip install --no-deps -e .
shumi --version

For manual installation, just make sure the following requirements are available on PATH, and then install shUMI via pip.

Requirements

• Python 3.10 or newer
• Python packages listed in pyproject.toml
• pbsim3, providing a pbsim executable and model files
• samtools (≥1.16)
• seqkit
• PacBio pbccs

Quick Start

Run a complete configured simulation:

shumi full \
  --config config/shumi-minimal.yaml \
  --run-id example \
  --out-dir shumi_out

Inspect available commands:

shumi --help
shumi full --help
shumi populate --help
shumi umi --help
shumi pcr --help
shumi pbccs --help
shumi ont --help
shumi genericseq --help

A complete reference of options available in each subcommand is included at config/shumi-reference.yaml.

Command Overview

shumi full : Run the end-to-end workflow: population generation/loading, UMI tagging and transcription, PCR, sequencing, visualizations, and template summary exports.

shumi populate : Create or load an input molecule population from a single sequence, an aligned population FASTA, or a random sequence.

shumi umi : Attach one or two UMI tags to input molecules. The default mode is umi; use --mode 2umi with --second-umi-tag for a second UMI-tagging stage.

shumi pcr : Amplify templates while modeling amplification efficiency, substitution errors, insertions, deletions, recombination, molecule sampling, and optional orientation shuffling.

shumi pbccs : Simulate PacBio CCS reads from PCR molecules using PBSIM, samtools, ccs, and seqkit.

shumi ont : Simulate ONT reads from PCR molecules using PBSIM and seqkit.

shumi genericseq : Generate simple substitution-only FASTQ/FASTA reads from PCR molecules.

Input Populations

The simulated population represents the assayed amplicon region, not the entire genome sequence. Population inputs are configured in the populate section:

• input_seq: a FASTA file containing one base sequence. Use with clonal: true for one haplotype, or with dirichlet and nhaps for a diversified population.
• input_pop: a FASTA file containing an existing population or alignment. Relative abundances are inferred from size=N annotations when present, or by counting duplicate sequences otherwise.
• random: generate a random sequence of the requested length.

Input FASTA sequences are normalized to DNA alphabet at the input boundary: U is converted to T, and ambiguity codes are rejected for molecule input except for alignment gaps in input_pop.

UMI Tag Schemas

UMI tag schemas support standard bases, IUPAC ambiguity symbols, grouped homopolymers such as (N)3, and explicit UMI blocks in square brackets. The bracketed region defines the UMI sequence that is tracked as the molecule tag. IUPAC ambiguous bases outside the brackets are treated as degenerate primer bases: they are sampled independently when the primer is attached, but they are not included in the UMI identity.

Use ambiguous symbols inside brackets when those positions are part of the UMI:

umi:
  mode: umi
  umi_tag: GGCCACGTACGTACGT[NNNNNNNN]AAGCGCTTAAGCGCTT
  fp: AATTGGCCAAGGTTCC

Keep ambiguous primer bases outside the brackets when they are part of the primer itself:

umi:
  mode: 2umi
  umi_tag: GTTTGGCACCTCGATGTCG[NNNYRNNNYRNNNYRNNN]GATCTCGGTGGTCGCCGTATCATT
  second_umi_tag: CAAGCAGAAGACGGCATACGAGAT[NNNYRNNNYRNNNYRNNN]AGRGTTYGATYMTGGCTCAG

Two-UMI mode uses a first tag and a second tag:

umi:
  mode: 2umi
  umi_tag: CAGGTACCTTTAAGACCAATGAC[NNNNNNNN]GTATAGTATAGTGAGTCGTCGGACGGAGCGG
  second_umi_tag: GATCTACACTCTTTCCCTACACGAC[NNNNNNNN]TAGGCATCTCCTATGGCAGGAAGAAG

First-tag error parameters use the tag_* prefix. Second-tag error parameters use the ss_tag_* prefix:

umi:
  tag_error_rate: 1.0e-4
  tag_recomb_rate: 0.0
  ss_tag_error_rate: 4.95e-6
  ss_tag_recomb_rate: 0.0

Configuration Files

Config files are YAML documents with sections matching pipeline stages.

global:
  seed: 271828
  threads: 4
  out_dir: shumi_out
  run_id: demo

populate:
  input_seq: sequences/sars-cov-2-spike-amplicon.fasta
  clonal: true

umi:
  mode: umi
  n_templates: 500
  umi_tag: GGCCACGTACGTACGT[NNNNNNNN]AAGCGCTTAAGCGCTT
  fp: AATTGGCCAAGGTTCC

pcr:
  cycles: 12
  sampled_molecules: 1000
  copy_prob: 0.95
  error_rate: 4.40e-5
  del_rate: 5.0e-6
  ins_rate: 5.0e-6
  recomb_rate: 1.0e-3
  max_molecules: 1000000

genericseq:
  reads: 1000
  error_rate: 0.01

Output-path interpolation is supported for string fields:

pcr:
  input_templates: "{umi.output_templates}"
  output_molecules: "{out_dir}/{run_id}.pcr_molecules.fasta"

genericseq:
  input_molecules: "{pcr.output_molecules}"
  output_fastq: "{out_dir}/{run_id}.genericseq.fastq.gz"

Stage Examples

Generate a clonal input molecule population:

shumi populate \
  --input-seq sequences/sars-cov-2-spike-amplicon.fasta \
  --clonal \
  --output shumi_out/demo.molecules.fasta \
  --output-abundances shumi_out/demo.abundances.tsv

Generate UMI-tagged templates:

shumi umi \
  --input-molecules shumi_out/demo.molecules.fasta \
  --abundances shumi_out/demo.abundances.tsv \
  --umi-tag 'GGCCACGTACGTACGT[NNNNNNNN]AAGCGCTTAAGCGCTT' \
  --fp AATTGGCCAAGGTTCC \
  --n-templates 100 \
  --output-templates shumi_out/demo.templates.fasta

Run PCR and sample output molecules:

shumi pcr \
  --input-templates shumi_out/demo.templates.fasta \
  --output-molecules shumi_out/demo.pcr_molecules.fasta \
  --cycles 6 \
  --copy-prob 0.95 \
  --sampled-molecules 2000

Generate generic sequencing reads:

shumi genericseq \
  --input-molecules shumi_out/demo.pcr_molecules.fasta \
  --output-fastq shumi_out/demo.genericseq.fastq.gz \
  --reads 100 \
  --error-rate 0.01

Generate PacBio CCS reads:

shumi pbccs \
  --input-molecules shumi_out/demo.pcr_molecules.fasta \
  --output-fastq shumi_out/demo.ccs.fastq.gz \
  --reads 1000 \
  --movie-time 24 \
  --overhead 1.5

pbccs samples input molecules without replacement before running the PacBio simulation, so reads * overhead must not exceed the number of available input molecules.

Generate ONT reads:

shumi ont \
  --input-molecules shumi_out/demo.pcr_molecules.fasta \
  --output-fastq shumi_out/demo.ont.fastq.gz \
  --reads 1000 \
  --model QSHMM-ONT-HQ.model \
  --mean-accuracy 0.92

Outputs

Output names depend on out_dir, run_id, and any explicit paths in the config or command line. A full run typically produces:

• input molecule population FASTA and abundance TSV
• UMI-tagged template FASTA and mutation TSV
• optional two-UMI template FASTA and mutation TSV
• sampled post-PCR molecule FASTA
• optional PCR history HDF5 file
• sequencing FASTQ/FASTA outputs for enabled sequencing stages
• sequencing read-to-molecule mapping files
• PCR and sequencing summary visualizations
• a template summary CSV with PCR and read representation

PCR molecule names include lineage metadata:

mol_1_template=<template_id>_copies=<n>_muts=<n>_recomb=<bool>_revcomp=<bool>

PCR history HDF5 files store enough structure to reconstruct sampled PCR molecules and inspect their lineages:

pcr_history.h5
  templates/
    ids
    sequences

  molecules/
    mol_<n>/
      attrs: parent_template_id, copy_count, total_mutations
      segments
      events/
        seg_<n>/
          substitutions
          deletions
          insertions

Reproducibility

Use --seed or global.seed to make stochastic stages reproducible. Use --threads or global.threads to control multiprocessing where supported.

Citation

Citation TBD.