Metadata-Version: 2.5
Name: MFED_for_Py
Version: 0.1.0rc3
Project-URL: Homepage, https://github.com/psimmond/MFED_for_Py
Project-URL: Repository, https://github.com/psimmond/MFED_for_Py
Project-URL: Issues, https://github.com/psimmond/MFED_for_Py/issues
Author-email: Peter Simmonds <peter.simmonds@utu.fi>, Terry Jones <terence.jones@charite.de>
License-Expression: LicenseRef-Academic-NonCommercial
License-File: LICENSE
Keywords: MFE,RNA,bioinformatics,sequence-analysis,structure
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
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: Topic :: Scientific/Engineering :: Bio-Informatics
Requires-Python: >=3.10
Requires-Dist: numpy>=1.24.0
Requires-Dist: openpyxl>=3.1.0
Requires-Dist: prseq>=0.0.36
Requires-Dist: rich>=13.0.0
Requires-Dist: sse-for-py>=0.1.0rc12
Requires-Dist: viennarna>=2.7.2
Description-Content-Type: text/markdown

# MFED_for_Py - Minimum Free Energy Difference Analysis

This package provides tools for detecting structured RNA elements by comparing
the minimum free energy (MFE) of native RNA/DNA sequence fragments to MFE
values of scrambled controls, using sequence-scrambling utilities from the
companion [SSE_for_Py](https://github.com/psimmond/SSE_for_Py) package.

## Installation

```bash
# Clone the repository
git clone <repository-url>
cd MFED_for_Py

# Install with uv
uv sync
```

## Overview

The `MFED_for_Py` package includes tools for:

1. **MFED Scanning** - Detecting structured RNA elements by comparing native
   sequences to scrambled controls
2. **Sequence Scrambling** - Preserving nucleotide/n-tuple/codon
   frequencies while randomizing sequences (provided by `SSE_for_Py`)
3. **MFE Calculation** - Computing minimum free energy for RNA secondary
   structures (via ViennaRNA)

## Main Tool: mfed-scan

The `mfed-scan` command is the primary tool for detecting structured RNA
elements using the MFED (Minimum Free Energy Difference) method.

### How It Works

1. **Fragment the input sequences** into overlapping windows of specified length
2. **Calculate MFE** for each native fragment using ViennaRNA
3. **Generate scrambled controls** using one of the `SSE_for_Py` scrambling
   methods (see below)
4. **Calculate MFE** for each scrambled sequence
5. **Compute statistics** (mean, standard deviation, Z-score) comparing native to scrambled sequences
6. **Output results** in tab-delimited or Excel format

A **negative Z-score** indicates the native sequence is more structured than expected by chance (lower MFE = more stable structure).

### Basic Usage

`--scramble-method` is required and must be given in uppercase; all examples
below use `NDR` (Nucleotide Distribution Randomization, dinucleotide
context).

```bash
# Basic usage with default parameters
mfed-scan input.fasta --scramble-method NDR

# Save output to a file
mfed-scan input.fasta --scramble-method NDR --output results.dat

# Save output as Excel
mfed-scan input.fasta --scramble-method NDR --output results.xlsx

# With custom fragment parameters
mfed-scan input.fasta --scramble-method NDR \
  --frag-length 200 \
  --frag-step 30 \
  --scrambles 10 \
  --output results.dat

# With multiprocessing and verbose output
mfed-scan input.fasta --scramble-method NDR \
  --processes 8 \
  --verbose \
  --output results.dat
```

### Command-Line Options

#### Required Arguments

- `fasta` - Input FASTA file containing one or more sequences
- `--scramble-method`, `--sm {CDLR,CLM,CLR,CLS,COR,NDR,NDS,NOR,NQR,NQS,NTR,NTS}` -
  Sequence scrambling method to use as the control, given in uppercase (see
  [Scrambling Methods](#scrambling-methods) below)

#### Fragment Parameters

- `--frag-length LENGTH` - Length of fragments to analyze (default: 200)
- `--frag-step STEP` - Step size between fragment start positions (default: 30)
  - Example: With length=200 and step=30, fragments will be at positions 1-200, 31-230, 61-260, etc.
- `--offset N` - Number of leading bases to skip in each input sequence before fragmenting begins (default: 0)

#### Scrambling Parameters

- `--scrambles N` - Number of scrambled sequences to generate per fragment (default: 10)
- `--max-retries N` - Maximum attempts per swap to find a valid pick before
  giving up on that swap. Used by the `CDLR` method and the NDR-family
  methods (`NDR`, `NTR`, `NQR`, `NDS`, `NTS`, `NQS`) (default: -1, meaning
  use the built-in default)
- `--min-swap-fraction F` - Minimum required fraction (0.0-1.0) of swap
  attempts that must succeed, else a `ScrambleError` is raised. Used by the
  `CDLR` method and the NDR-family methods (default: 0.0)
- `--genetic-code N` - Genetic code (1-23). Used by the `CDLR`, `CLS`, and
  `COR` methods (default: 1)
- `--failed-scramble-action {ignore,warn,exit}` - How to handle a fragment
  whose scrambling fails: drop it silently, drop it and warn, or abort
  immediately (default: warn)

NDR-family options (used only when `--scramble-method` is one of `NDR`,
`NTR`, `NQR`, `NDS`, `NTS`, `NQS`):

- `--no-swaps N` - Number of swap operations to perform (default: `frag_length / 2`)

#### Performance Options

- `--processes N` - Number of parallel processes to use (default: CPU cores - 2)
- `--cache-size SIZE` - Total LRU cache size for MFE calculations across all
  workers, e.g. `100M`, `1G`, `500K` (default: `600M`); divided among worker
  processes
- `--verbose` - Display detailed progress and timing information with colored output
- `--scramble-only` - Stop after the scrambling phase, skipping MFE calculation

#### Output Options

- `--output PATH` - Output file path (default: stdout)
  - If ends with `.xlsx`: writes Excel format with formatting
  - Otherwise: writes tab-delimited text
- `--fasta-output PATH` - Optional FASTA file to write all native and scrambled sequences
- `--seed N` - Random seed for reproducible results (default: random)

### Output Format

#### Tab-Delimited / Excel Output

The main output contains these columns:

| Column | Description |
|--------|-------------|
| Sequence | Name of the input sequence |
| Fragment | Fragment number (1, 2, 3, ...) |
| Start | Start position in original sequence (1-based) |
| End | End position in original sequence (1-based, inclusive) |
| Native | MFE of the native sequence fragment |
| Scr_1, Scr_2, ... | MFE of each scrambled sequence |
| Mean | Mean MFE of all scrambled sequences |
| StdDev | Standard deviation of scrambled MFE values |
| Z-score | (Native - Mean) / StdDev |

**Interpreting Z-scores:**
- **Negative Z-score**: Native sequence is more structured (lower/more negative MFE) than scrambled controls
- **Positive Z-score**: Native sequence is less structured than scrambled controls
- **Z-score near 0**: Native sequence has similar structure to random sequences

#### FASTA Output (Optional)

When `--fasta-output` is specified, a FASTA file is created containing all native and scrambled sequences with descriptive IDs:

```
>SequenceName_Frag=1_Start=1_End=200_Native
ACGTACGT...
>SequenceName_Frag=1_Start=1_End=200_Scramble=1
ACGTTACG...
>SequenceName_Frag=1_Start=1_End=200_Scramble=2
ATCGACGT...
```

### Examples

#### Example 1: Quick scan with minimal scrambles

```bash
mfed-scan viral_genomes.fasta --scramble-method NDR \
  --frag-length 150 \
  --frag-step 50 \
  --scrambles 5 \
  --seed 42 \
  --output quick_scan.dat
```

#### Example 2: High-resolution scan with many controls

```bash
mfed-scan candidate_sequences.fasta --scramble-method NQR \
  --frag-length 200 \
  --frag-step 10 \
  --scrambles 100 \
  --processes 8 \
  --verbose \
  --seed 42 \
  --output detailed_results.xlsx \
  --fasta-output all_sequences.fasta
```

#### Example 3: Sequential (nearest-neighbor) dinucleotide-preserving scrambles

```bash
mfed-scan sequences.fasta --scramble-method NDS \
  --frag-length 200 \
  --frag-step 30 \
  --scrambles 10 \
  --output results.dat
```

#### Example 4: Reproducible analysis with parallel processing

```bash
# Run 1
mfed-scan input.fasta --scramble-method NDR --seed 12345 --processes 10 --output run1.dat

# Run 2 (will produce identical results)
mfed-scan input.fasta --scramble-method NDR --seed 12345 --processes 10 --output run2.dat
```

### Performance Considerations

The tool uses **two-phase multiprocessing** for optimal performance:

1. **Phase 1**: Scramble all fragments in parallel
2. **Phase 2**: Calculate MFE for all sequences in parallel
3. **Phase 3**: Aggregate results and write output

**Performance factors:**

- **Fragment length**: Longer fragments take more time to fold (MFE calculation is O(n^3))
- **Number of scrambles**: More scrambles = better statistics but longer runtime
- **Step size**: Smaller steps = more fragments = longer runtime but higher resolution
- **Number of processes**: More processes = faster runtime (up to CPU limit)
- **Sequence length**: Total runtime scales with number of fragments generated

**Recommended settings for different use cases:**

- **Quick exploratory scan**: `--frag-length 150 --frag-step 50 --scrambles 10 --processes 4`
- **Standard analysis**: `--frag-length 200 --frag-step 30 --scrambles 100 --processes 8`
- **High-resolution scan**: `--frag-length 200 --frag-step 10 --scrambles 100 --processes 16`

### Verbose Output Example

With `--verbose`, the tool displays colored progress information:

```
───────────────────────── MFED Scanner ─────────────────────────
Processes: 8
Fragment length: 200
Fragment step: 30
Scrambles per fragment: 10
Random seed: 42

Phase 0: Reading sequences and extracting fragments...
  ✓ Found 5 sequences, 150 fragments, 1650 total jobs (0.01s)

Phase 1: Scrambling 1650 fragments using 8 processes...
  ✓ Scrambling complete (2.34s)

Phase 2: Calculating MFE for 1650 sequences using 8 processes...
  ✓ MFE calculation complete (12.56s)

Phase 3: Aggregating results and writing output...
  ✓ Output written (0.23s)

───────────────────────── Complete ─────────────────────────
```

## Scrambling Methods

`--scramble-method` selects the control-generation method, provided by the
`SSE_for_Py` package. Method names must be given in **uppercase**.

| Method | Name | Description |
|--------|------|--------------|
| `NDR` | Nucleotide Di-tuple Randomization | Random swaps preserving dinucleotide (n_tuple=2) context |
| `NTR` | Nucleotide Tri-tuple Randomization | Random swaps preserving trinucleotide (n_tuple=3) context |
| `NQR` | Nucleotide Quad-tuple Randomization | Random swaps preserving tetranucleotide (n_tuple=4) context |
| `NDS` | Nucleotide Di-tuple Sequential | Sequential (nearest-neighbor) swaps preserving dinucleotide context |
| `NTS` | Nucleotide Tri-tuple Sequential | Sequential (nearest-neighbor) swaps preserving trinucleotide context |
| `NQS` | Nucleotide Quad-tuple Sequential | Sequential (nearest-neighbor) swaps preserving tetranucleotide context |
| `NOR` | Nucleotide Only Randomization | Randomly shuffles nucleotides with no constraints, preserving only gaps (`-`) and ambiguous bases (`N`) in place |
| `CDLR` | Codon Distribution at Leucine and aRginine | Randomizes codon choice at Leucine/Arginine positions subject to genetic-code constraints |
| `CLM` | Codon Like Maximal | Maximally scrambles a coding sequence by randomly choosing synonymous codons at each position |
| `CLR` | Codon Like Randomization | Scrambles codons encoding the same amino acid, preserving the amino acid sequence |
| `CLS` | Codon Like Swap | Swaps pairs of codons encoding the same amino acid, preserving the amino acid sequence |
| `COR` | Codon Only Randomization | Randomly shuffles whole codons without constraints, preserving codon structure |

The six `N**` methods (`NDR`, `NTR`, `NQR`, `NDS`, `NTS`, `NQS`) are all
presets of the same underlying nucleotide-swap algorithm, each fixing a
tuple size (Di/Tri/Quad, i.e. 2/3/4) and a swap style (Random/Sequential),
and all accept `--no-swaps`.

### NDR-Family Algorithm Details

The `N**` methods implement n-tuple preserving sequence scrambling:

1. **Scan sequence** to categorize positions by their flanking nucleotide context
2. **Encode positions** using a 4-ary encoding of flanking bases (A=0, C=1, G=2, T=3)
3. **Randomly (`*R`) or sequentially (`*S`) select** positions with identical flanking context
4. **Swap nucleotides** at these positions, preserving the n-tuple frequencies
5. **Validate swaps** to avoid overlapping positions

These methods preserve:
- Mononucleotide frequencies (always)
- Dinucleotide frequencies (`ND*` methods, and higher-order methods)
- Higher-order n-tuple frequencies (`NT*` and `NQ*` methods)
- Overall sequence composition

### MFE Calculation

Uses the ViennaRNA package to:
1. Predict the minimum free energy secondary structure
2. Return the energy value in kcal/mol
3. Lower (more negative) values indicate more stable structures

## Dependencies

- Python >= 3.10
- ViennaRNA - RNA secondary structure prediction
- NumPy - Numerical operations
- SSE_for_Py - Sequence scrambling utilities
- prseq - Sequence processing utilities
- openpyxl - Excel file output
- rich - Colored terminal output

## License

See [LICENSE](LICENSE).

## Contact

For questions, issues, or contributions, please contact:
- Peter Simmonds (peter.simmonds@utu.fi)
- Terry Jones (terence.jones@charite.de)
