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 package.

Installation

# 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).

# 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 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:

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

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

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

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

# 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.

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.

Contact

For questions, issues, or contributions, please contact:

• Peter Simmonds (peter.simmonds@utu.fi)
• Terry Jones (terence.jones@charite.de)