Metadata-Version: 2.5
Name: MFED_for_Py
Version: 0.1.0rc2
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
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; all examples below use `ndr`
(Nucleotide Distribution Randomization).

```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,nor}` - Sequence
  scrambling method to use as the control (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` and `ndr` methods (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` and `ndr` 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-specific options (used only when `--scramble-method ndr`):

- `--s-method {2,3}` - NDR sub-method (default: 2)
  - `2` = random swaps
  - `3` = sequential (nearest-neighbor) swaps
- `--n-tuple N` - Size of tuple to preserve during scrambling (default: 3)
  - `2` = preserve dinucleotide frequencies
  - `3` = preserve trinucleotide context
  - `4` = preserve tetranucleotide context
- `--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 ndr \
  --frag-length 200 \
  --frag-step 10 \
  --scrambles 100 \
  --n-tuple 4 \
  --processes 8 \
  --verbose \
  --seed 42 \
  --output detailed_results.xlsx \
  --fasta-output all_sequences.fasta
```

#### Example 3: Dinucleotide-preserving scrambles

```bash
mfed-scan sequences.fasta --scramble-method ndr \
  --frag-length 200 \
  --frag-step 30 \
  --scrambles 10 \
  --n-tuple 2 \
  --s-method 2 \
  --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 | Name | Description |
|--------|------|--------------|
| `ndr` | Nucleotide Distribution Randomization | Swaps nucleotides sharing the same flanking n-tuple context, preserving n-tuple frequencies (`--s-method`, `--n-tuple`, `--no-swaps`) |
| `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 |

### NDR Algorithm Details

The default `ndr` method implements 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 select** positions with identical flanking context
4. **Swap nucleotides** at these positions, preserving the n-tuple frequencies
5. **Validate swaps** to avoid overlapping positions

This method preserves:
- Mononucleotide frequencies (always)
- Dinucleotide frequencies (when n_tuple >= 2)
- Higher-order n-tuple frequencies (when n_tuple > 2)
- 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)
