Metadata-Version: 2.4
Name: SSE_for_Py
Version: 0.1.0rc9
Author-email: Peter Simmonds <peter.simmonds@utu.fi>, Terry Jones <terence.jones@charite.de>
License-Expression: LicenseRef-Academic-NonCommercial
Project-URL: Homepage, https://github.com/psimmond/SSE_for_Py
Project-URL: Repository, https://github.com/psimmond/SSE_for_Py
Project-URL: Issues, https://github.com/psimmond/SSE_for_Py/issues
Keywords: RNA,structure,MFE,bioinformatics,sequence-analysis,scrambling
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
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
Classifier: Programming Language :: C
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: prseq>=0.0.38
Dynamic: license-file

# PRE-RELEASE

**The code in this repository is under active development and testing
prior to an official release. Until that time, the project should be
considered unstable and liable to change without warning.**

# pysse — sequence scrambling and mutation tools

`pysse` provides tools for **scrambling** and **mutating** RNA/DNA genome
sequences while holding chosen properties of those sequences constant —
for example, preserving the encoded protein while randomizing codon usage,
or nudging a sequence's dinucleotide composition (e.g., CpG content) toward
a target value. These operations are used to build null-model or control
sequences for studying compositional selection pressures in viral and
other genomes.

The scrambling and mutation algorithms were originally written in
PowerBasic. That code was translated into C, and Python bindings were
added on top of the C extension, all written by Claude, Anthropic's AI
coding agent. The C code was checked (via the Python bindings) for
equivalence with the original PowerBasic implementation (see
[Testing](#testing) below).

## Methods

The following methods are implemented identically in the PowerBasic
original and the C code (for which Python bindings are provided). They
are described here in a language-independent way; see
[Implementations](#implementations) below for how to call them from
either language.

All methods operate on a nucleotide sequence (A, C, G, T). Methods whose
description mentions "codons" require the sequence length to be a
multiple of 3 and treat it as a coding sequence read in one frame,
starting at the first base.

- **NDR** (Nucleotide Distribution Randomization) — repeatedly swaps pairs
  of individual nucleotides, choosing swaps that preserve the sequence's
  n-tuple (e.g., trinucleotide) frequency distribution.
- **CDLR** (Codon Distribution at Leucine and aRginine) — swaps
  nucleotides between synonymous codons so as to preserve dinucleotide
  frequencies while leaving the encoded protein sequence unchanged.
- **NOR** (Nucleotide Only Randomization) — shuffles all nucleotides in
  the sequence with no constraints, other than leaving gap (`-`) and
  ambiguous (`N`) characters in their original positions. This destroys
  both the reading frame and the encoded protein.
- **COR** (Codon Only Randomization) — shuffles whole codons to new
  positions with no constraints, preserving the overall codon composition
  but not the encoded protein sequence or codon-to-position mapping.
- **CLR** (Codon Like Randomization) — shuffles synonymous codons (those
  encoding the same amino acid) among the positions that encode that
  amino acid, preserving the protein sequence while randomizing synonymous
  codon usage.
- **CLS** (Codon Like Swap) — like CLR, but performs the randomization as
  a series of pairwise swaps between synonymous codons, rather than a
  full shuffle.
- **CLM** (Codon Like Maximal) — for each codon position, independently
  draws a new, random synonymous codon, maximally randomizing synonymous
  codon usage while preserving the protein sequence.

All of the codon-aware methods (CDLR, COR, CLR, CLS) accept a genetic code
number (1–23, following the NCBI genetic code tables) so that non-standard
codon tables can be used.

In addition to scrambling, a **mutation** method nudges a sequence toward
a target frequency for a chosen dinucleotide (e.g., CpG or ApT), expressed
as a ratio relative to the frequency expected from the sequence's base
composition. It can preserve the encoded protein, optionally preserve
overall mononucleotide frequencies while it mutates, and can be told to
leave specific other (named) dinucleotides unchanged.

## Implementations

### Python

#### Installation

`pysse` ships as a compiled Python package (the C code is built into a
Python extension module), with prebuilt wheels for common platforms.

Using pip:

```bash
pip install pysse
```

Using [uv](https://docs.astral.sh/uv/):

```bash
uv add pysse
```

#### Command-line tools

Installing `pysse` provides the following commands:

- **`scramble`** — scramble sequences in a FASTA file (or from stdin, or a
  single sequence given directly) using any of the methods described
  above.

  ```bash
  # Scramble a FASTA file using the CDLR method, with genetic code 1
  scramble cdlr input.fasta --genetic-code 1 --output scrambled.fasta

  # Scramble a single sequence using NOR, with a given random seed
  scramble nor --sequence ACGTACGTACGT --seed 42

  # Read from stdin, scramble with NDR
  scramble ndr --n-tuple 3 < input.fasta > scrambled.fasta
  ```

- **`mutate`** — mutate sequences in a FASTA file (or from stdin, or a
  single sequence) toward a target dinucleotide frequency ratio.

  ```bash
  # Reduce/increase CpG toward a 2x target ratio (the default)
  mutate input.fasta --output mutated.fasta

  # Target ApT at a different ratio, with a fixed seed
  mutate input.fasta --dinuc-mutate ApT --diverge 1.5 --seed 123
  ```

- **`encode-alignment`** — convert a FASTA alignment into the internal,
  digit-encoded alignment format used elsewhere in the package.

  ```bash
  encode-alignment alignment.fasta --output alignment.enc
  ```

- **`sse-benchmark`** — run performance benchmarks over the scrambling and
  mutation methods, optionally writing results as HTML, Markdown, or CSV.

  ```bash
  sse-benchmark --markdown results.md
  ```

All of the FASTA-processing commands (`scramble`, `mutate`) share a common
set of options, including `--seed` (for reproducibility), `--output`,
`--wrap` (output line wrapping), `--repeat` (apply the operation more than
once per input sequence), and `--checksum` (print an Adler-32 checksum of
the output, used when validating against the PowerBasic reference — see
[Testing](#testing)).

#### Library functions

The high-level, recommended entry points are `SSE_for_Py.scramble.scramble`
(a generator that dispatches to any scrambling method) and `SSE_for_Py.mutate.mutate`:

```python
from SSE_for_Py.scramble import scramble

# Scramble a single sequence directly
for seq_id, scrambled in scramble("cdlr", sequence="ATGCGTAAATAG", genetic_code=1):
    print(scrambled)

# Scramble every sequence in a FASTA file
for seq_id, scrambled in scramble("ndr", "input.fasta", n_tuple=3):
    print(f"{seq_id}: {scrambled}")
```

```python
from SSE_for_Py.mutate import mutate

result = mutate("ATGCGTAAATAG", dinuc_mutate="CpG", dinuc_diverge=2.0)
```

Each scrambling method is also available directly as its own function,
if you don't want to go through the `scramble` dispatcher:

```python
from SSE_for_Py.scramble_ndr import scramble_ndr
from SSE_for_Py.scramble_cdlr import scramble_cdlr
from SSE_for_Py.scramble_nor import scramble_nor
from SSE_for_Py.scramble_cor import scramble_cor
from SSE_for_Py.scramble_clr import scramble_clr
from SSE_for_Py.scramble_cls import scramble_cls
from SSE_for_Py.scramble_clm import scramble_clm

scramble_ndr("ACGTACGTACGT", s_method=2, no_swaps=10, n_tuple=3)
scramble_cdlr("ATGCGTAAATAG", genetic_code=1)
scramble_nor("ATGCGTAAATAG")
scramble_cor("ATGCGTAAATAG", genetic_code=1)
scramble_clr("ATGCGTAAATAG", genetic_code=1)
scramble_cls("ATGCGTAAATAG", genetic_code=1)
scramble_clm("ATGCGTAAATAG")
```

Invalid arguments (e.g., a sequence whose length isn't a multiple of 3 for
a codon-aware method) raise `ValueError`. Failures in the underlying C
code, or a scramble that could not make enough progress to be considered
meaningful (see the `min_swap_fraction` argument on `scramble_ndr` and
`scramble_cdlr`), raise `SSE_for_Py.exceptions.ScrambleError`.

### C library

The algorithms themselves live in a small, dependency-free C library
under `src/c/`, with one source file per method (`scramble_ndr.c`,
`scramble_cdlr.c`, `scramble_nor.c`, `scramble_cor.c`, `scramble_clr.c`,
`scramble_cls.c`, `scramble_clm.c`), plus `mutate_sequence.c` for the
mutation method and `lib.c`/`lib.h` for shared utilities (codon table
setup and translation, RNG seeding).

All of that is compiled into a single Python extension module,
`SSE_for_Py._combined`, built from `combined_module.c` plus the files above.
Everything is built as one extension, rather than one per method, so
that every method shares a single C-level RNG state — important for
reproducibility when a script calls more than one scrambling method with
the same seed.

The C code has no dependencies beyond the C standard library and `libm`,
and can also be compiled and used standalone, outside of Python — see the
`benchmark-bin` target in the `Makefile` for an example of linking it into
a plain C binary.

## Testing

Clone the repository and install the development dependencies with
[uv](https://docs.astral.sh/uv/):

```bash
uv sync
```

Run the test suite:

```bash
make test
```

which is equivalent to `uv run pytest`. To run the tests across every
supported Python version (currently 3.10 through 3.14), use:

```bash
make nox
```

### Validating against the PowerBasic reference

Because the C code is a translation of the original PowerBasic
implementation, the C extension can optionally be compiled with a
simplified, deterministic random number generator instead of the
platform's real one. Building with this test RNG (`SSE_TEST_RNG=1`)
makes the C code draw from the exact same pseudo-random integer
sequence as the PowerBasic code, given the same seed. This allows the
two implementations' outputs be compared directly rather than only
statistically.

That comparison is done by running the same input through both
implementations and comparing an Adler-32 checksum of their output
(the `--checksum` option on the `scramble` and `mutate` commands
prints this).  A checksum match on a given input and seed is strong
evidence that the C translation faithfully reproduces the PowerBasic
original's logic for that case. This technique was used throughout
development to verify the C translation against the PowerBasic
reference implementation, and is exercised by the test suite via:

```bash
make nox-with-rng-checksum
```

which builds the extension with both `SSE_TEST_RNG=1` (the shared
deterministic RNG) and `PREVENT_REVERSIONS=1` (an extra guard against a
specific class of double-free bug) before running the tests.

### Benchmarking

```bash
make benchmark
```

builds the C benchmark binary and the Python extension, then runs
`sse-benchmark` to compare performance across methods and implementations.

### benchmark results

# SSE benchmark results

- **Sequence lengths:** 100, 1000, 10000
- **Sequences per length:** 100
- **Length SD:** 0.0
- **Methods:** cdlr, clm, clr, cls, cor, ndr, nor, translate
- **Repeats per sequence:** 3
- **Warmup calls per sequence:** 0
- **Seed:** 1
- **Generated:** 2026-08-24 16:08:09 UTC

## Nucleotides per second

### Sequence length 100

| Method | C nt/s | Python nt/s | C / Python |
| --- | ---: | ---: | ---: |
| cdlr | 3,169,342 | 3,139,684 | 1.0x |
| clm | 197,419,355 | 164,881,017 | 1.2x |
| clr | 8,913,487 | 8,502,869 | 1.0x |
| cls | 9,161,677 | 8,994,179 | 1.0x |
| cor | 53,310,105 | 45,940,360 | 1.2x |
| ndr | 3,846,154 | 3,892,160 | 1.0x |
| nor | 25,478,768 | 24,327,750 | 1.0x |
| translate | 209,589,041 | 225,625,582 | 0.9x |

### Sequence length 1000

| Method | C nt/s | Python nt/s | C / Python |
| --- | ---: | ---: | ---: |
| cdlr | 6,378,915 | 6,357,758 | 1.0x |
| clm | 229,115,854 | 210,912,714 | 1.1x |
| clr | 8,733,548 | 8,519,392 | 1.0x |
| cls | 10,308,288 | 9,736,246 | 1.1x |
| cor | 41,582,515 | 40,226,251 | 1.0x |
| ndr | 10,796,638 | 10,718,490 | 1.0x |
| nor | 18,030,230 | 18,298,231 | 1.0x |
| translate | 136,202,990 | 128,342,080 | 1.1x |

### Sequence length 10000

| Method | C nt/s | Python nt/s | C / Python |
| --- | ---: | ---: | ---: |
| cdlr | 6,449,089 | 6,439,097 | 1.0x |
| clm | 231,563,513 | 210,873,606 | 1.1x |
| clr | 8,285,014 | 8,196,008 | 1.0x |
| cls | 10,424,397 | 10,263,846 | 1.0x |
| cor | 35,344,421 | 34,134,994 | 1.0x |
| ndr | 11,555,925 | 11,540,557 | 1.0x |
| nor | 14,413,142 | 13,958,460 | 1.0x |
| translate | 26,553,512 | 26,086,465 | 1.0x |

<details><summary>Call/timing details</summary>

#### Sequence length 100

| Method | Language | Calls | Skipped | Discarded | Total nt | Total seconds | Min (s) | Max (s) | StdDev (s) |
| --- | --- | ---: | ---: | ---: | ---: | ---: | ---: | ---: | ---: |
| cdlr | C | 300 | 0 | 0 | 30,600 | 0.0097 | 0.000025 | 0.000101 | 0.000005 |
| cdlr | Python | 300 | 0 | 0 | 30,600 | 0.0097 | 0.000025 | 0.000112 | 0.000007 |
| clm | C | 300 | 0 | 0 | 30,600 | 0.0002 | 0.000000 | 0.000001 | 0.000000 |
| clm | Python | 300 | 0 | 0 | 30,600 | 0.0002 | 0.000001 | 0.000010 | 0.000001 |
| clr | C | 300 | 0 | 0 | 30,600 | 0.0034 | 0.000010 | 0.000023 | 0.000001 |
| clr | Python | 300 | 0 | 0 | 30,600 | 0.0036 | 0.000011 | 0.000022 | 0.000002 |
| cls | C | 300 | 0 | 0 | 30,600 | 0.0033 | 0.000009 | 0.000016 | 0.000001 |
| cls | Python | 300 | 0 | 0 | 30,600 | 0.0034 | 0.000010 | 0.000019 | 0.000001 |
| cor | C | 300 | 0 | 0 | 30,600 | 0.0006 | 0.000001 | 0.000008 | 0.000001 |
| cor | Python | 300 | 0 | 0 | 30,600 | 0.0007 | 0.000002 | 0.000012 | 0.000001 |
| ndr | C | 300 | 0 | 0 | 30,600 | 0.0080 | 0.000017 | 0.000059 | 0.000006 |
| ndr | Python | 300 | 0 | 0 | 30,600 | 0.0079 | 0.000016 | 0.000070 | 0.000006 |
| nor | C | 300 | 0 | 0 | 30,600 | 0.0012 | 0.000003 | 0.000012 | 0.000001 |
| nor | Python | 300 | 0 | 0 | 30,600 | 0.0013 | 0.000004 | 0.000008 | 0.000001 |
| translate | C | 300 | 0 | 0 | 30,600 | 0.0001 | 0.000000 | 0.000009 | 0.000001 |
| translate | Python | 300 | 0 | 0 | 30,600 | 0.0001 | 0.000000 | 0.000003 | 0.000000 |

#### Sequence length 1000

| Method | Language | Calls | Skipped | Discarded | Total nt | Total seconds | Min (s) | Max (s) | StdDev (s) |
| --- | --- | ---: | ---: | ---: | ---: | ---: | ---: | ---: | ---: |
| cdlr | C | 300 | 0 | 0 | 300,600 | 0.0471 | 0.000140 | 0.000271 | 0.000014 |
| cdlr | Python | 300 | 0 | 0 | 300,600 | 0.0473 | 0.000138 | 0.000226 | 0.000012 |
| clm | C | 300 | 0 | 0 | 300,600 | 0.0013 | 0.000003 | 0.000015 | 0.000002 |
| clm | Python | 300 | 0 | 0 | 300,600 | 0.0014 | 0.000005 | 0.000013 | 0.000001 |
| clr | C | 300 | 0 | 0 | 300,600 | 0.0344 | 0.000106 | 0.000213 | 0.000010 |
| clr | Python | 300 | 0 | 0 | 300,600 | 0.0353 | 0.000107 | 0.000174 | 0.000008 |
| cls | C | 300 | 0 | 0 | 300,600 | 0.0292 | 0.000089 | 0.000158 | 0.000008 |
| cls | Python | 300 | 0 | 0 | 300,600 | 0.0309 | 0.000091 | 0.000135 | 0.000007 |
| cor | C | 300 | 0 | 0 | 300,600 | 0.0072 | 0.000021 | 0.000038 | 0.000003 |
| cor | Python | 300 | 0 | 0 | 300,600 | 0.0075 | 0.000022 | 0.000046 | 0.000004 |
| ndr | C | 300 | 0 | 0 | 300,600 | 0.0278 | 0.000082 | 0.000187 | 0.000010 |
| ndr | Python | 300 | 0 | 0 | 300,600 | 0.0280 | 0.000082 | 0.000229 | 0.000013 |
| nor | C | 300 | 0 | 0 | 300,600 | 0.0167 | 0.000050 | 0.000090 | 0.000005 |
| nor | Python | 300 | 0 | 0 | 300,600 | 0.0164 | 0.000050 | 0.000081 | 0.000005 |
| translate | C | 300 | 0 | 0 | 300,600 | 0.0022 | 0.000006 | 0.000014 | 0.000001 |
| translate | Python | 300 | 0 | 0 | 300,600 | 0.0023 | 0.000007 | 0.000018 | 0.000002 |

#### Sequence length 10000

| Method | Language | Calls | Skipped | Discarded | Total nt | Total seconds | Min (s) | Max (s) | StdDev (s) |
| --- | --- | ---: | ---: | ---: | ---: | ---: | ---: | ---: | ---: |
| cdlr | C | 300 | 0 | 0 | 3,000,600 | 0.4653 | 0.001424 | 0.001747 | 0.000067 |
| cdlr | Python | 300 | 0 | 0 | 3,000,600 | 0.4660 | 0.001408 | 0.002351 | 0.000094 |
| clm | C | 300 | 0 | 0 | 3,000,600 | 0.0130 | 0.000038 | 0.000059 | 0.000004 |
| clm | Python | 300 | 0 | 0 | 3,000,600 | 0.0142 | 0.000042 | 0.000064 | 0.000004 |
| clr | C | 300 | 0 | 0 | 3,000,600 | 0.3622 | 0.001108 | 0.001386 | 0.000055 |
| clr | Python | 300 | 0 | 0 | 3,000,600 | 0.3661 | 0.001108 | 0.001386 | 0.000056 |
| cls | C | 300 | 0 | 0 | 3,000,600 | 0.2878 | 0.000873 | 0.001066 | 0.000041 |
| cls | Python | 300 | 0 | 0 | 3,000,600 | 0.2923 | 0.000885 | 0.001074 | 0.000044 |
| cor | C | 300 | 0 | 0 | 3,000,600 | 0.0849 | 0.000263 | 0.000405 | 0.000020 |
| cor | Python | 300 | 0 | 0 | 3,000,600 | 0.0879 | 0.000265 | 0.000389 | 0.000019 |
| ndr | C | 300 | 0 | 0 | 3,000,600 | 0.2597 | 0.000783 | 0.001266 | 0.000054 |
| ndr | Python | 300 | 0 | 0 | 3,000,600 | 0.2600 | 0.000789 | 0.001362 | 0.000065 |
| nor | C | 300 | 0 | 0 | 3,000,600 | 0.2082 | 0.000644 | 0.000829 | 0.000032 |
| nor | Python | 300 | 0 | 0 | 3,000,600 | 0.2150 | 0.000646 | 0.000852 | 0.000039 |
| translate | C | 300 | 0 | 0 | 3,000,600 | 0.1130 | 0.000346 | 0.000446 | 0.000020 |
| translate | Python | 300 | 0 | 0 | 3,000,600 | 0.1150 | 0.000351 | 0.000488 | 0.000023 |

</details>
