Metadata-Version: 2.4
Name: sxlaep
Version: 1.0.5
Summary: sxLaep: Rapid Enzyme/Non-Enzyme Prediction using sequence features and XGBoost
Author-email: DHY <dhy.scut@outlook.com>
License: MIT
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Intended Audience :: Science/Research
Requires-Python: >=3.9
Description-Content-Type: text/markdown
Requires-Dist: numpy
Requires-Dist: pandas
Requires-Dist: scikit-learn
Requires-Dist: xgboost
Requires-Dist: joblib

# sxlaep: Rapid Enzyme/Non-Enzyme Prediction

[![PyPI version](https://img.shields.io/pypi/v/sxlaep.svg)](https://pypi.org/project/sxlaep/)
[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT)
[![Python](https://img.shields.io/badge/python-3.9-blue.svg)](https://www.python.org/)

**sxlaep (Rapid Enzyme/Non-Enzyme Prediction)** is an efficient enzyme/non-enzyme prediction tool for protein sequences. It is built on multi-physicochemical property features and the XGBoost machine learning algorithm.

## Features

- **Efficient Prediction**: Fast and accurate enzyme/non-enzyme classification using optimized feature extraction and a pre-trained XGBoost model.
- **Multi-Mode Support**: Single-sequence prediction, batch prediction, and FASTA file batch prediction.
- **Rich Feature Set**: Multi-property pseudo-amino acid composition (Pseudo-AAC), CTD features, and windowed amino acid composition.
- **Simple API**: Clean Python API with standalone functions — no class instantiation needed.
- **Command-Line Interface**: Predict directly from the terminal with a single command.
- **Parallel Processing**: Multi-process feature extraction for large-scale datasets.

## Dependencies

- `numpy`, `pandas`, `scikit-learn`, `xgboost`, `joblib`
- **Requirements**: Python 3.9 or higher.

## Installation

```bash
pip install sxlaep
```

## Quick Start

### Command Line

```bash
# Shorthand: use the bundled pre-trained model
sxlaep --input sequences.fasta --output predictions.csv

# Or with explicit model path
sxlaep predict --model enzyme_xgb_model.ubj --fasta sequences.fasta --output predictions.csv
```

### Python API

```python
from pathlib import Path
import sxlaep
from sxlaep.model import load_model, predict_sequences

# Load the bundled model
ubj = Path(sxlaep.__file__).resolve().parent / "enzyme_xgb_model.ubj"
model = load_model(ubj)

# Single sequence
df = predict_sequences(model, ["MKVLWALIFLLKSAF"])
print(df)  # columns: pred_label, enzyme_probability
```

FASTA to CSV in one call:

```python
from sxlaep import run_prediction_pipeline

df = run_prediction_pipeline(
    model_path="enzyme_xgb_model.ubj",
    fasta_path="sequences.fasta",
    output_csv="predictions.csv",
)
print(df[["sequence_id", "pred_label", "enzyme_probability"]].head())
```

## API Reference

### Model Loading

**`load_model(model_path)`**

Load an XGBoost classifier from native UBJSON (`.ubj`), JSON, or legacy joblib pickle (`.pkl`/`.joblib`).

```python
from sxlaep.model import load_model

model = load_model("enzyme_xgb_model.ubj")
```

### Prediction

**`predict_sequences(model, sequences, config=None, n_jobs=1)`**

Predict enzyme probabilities for a list of protein sequences. Returns a DataFrame with columns `pred_label` (0/1) and `enzyme_probability` (float).

**`predict_fasta(model_path, fasta_path, output_csv, config=None, n_jobs=1)`**

Load a model, predict all sequences in a FASTA file, and write results to CSV. Returns a DataFrame with columns `sequence_id`, `description`, `pred_label`, `enzyme_probability`.

**`run_prediction_pipeline(model_path, fasta_path, output_csv, feature_config=None, n_jobs=1)`**

One-shot convenience: `predict_fasta` with a slightly different signature.

### Feature Extraction

**`extract_features(sequence, config=None)`** — Extract the full feature vector for one sequence.

**`extract_feature_matrix(sequences, config=None, n_jobs=1)`** — Extract a feature matrix from multiple sequences.

**`sanitize_sequence(sequence)`** — Clean a sequence to only the 20 standard amino acids.

### FASTA Utilities

**`read_fasta(path)`** — Return a list of sequence strings from a FASTA file.

**`read_fasta_records(path)`** — Return a list of `FastaRecord` objects with `identifier`, `description`, and `sequence` fields.

### Configuration

**`FeatureConfig(lag=10, weight=0.05, n_segments=3, add_length=True, properties=...)`**

Configure feature extraction parameters. The bundled model uses defaults; customize only when loading a custom-trained model.

## CLI Reference

### Shorthand (bundled model)

```bash
sxlaep --input proteins.fasta --output predictions.csv
```

| Argument | Required | Description |
|----------|----------|-------------|
| `--input`, `-i` | Yes | Input FASTA path |
| `--output`, `-o` | No | Output CSV path (default: `sxlaep_predictions.csv`) |

### `sxlaep predict` (custom model)

```bash
sxlaep predict --model path/to/model.ubj --fasta sequences.fasta --output results.csv
```

| Argument | Required | Description |
|----------|----------|-------------|
| `--model` | Yes | Model path (`.ubj` / `.pkl` / `.joblib`) |
| `--fasta` | Yes | Input FASTA path |
| `--output` | No | Output CSV path (default: `results/predictions.csv`) |
| `--lag` | No | Pseudo-AAC lag (default: 10) |
| `--weight` | No | Pseudo-AAC weight (default: 0.05) |
| `--segments` | No | Window-AAC segments (default: 3) |
| `--add-length` / `--no-add-length` | No | Append sequence length (default: enabled) |
| `--properties` | No | Physicochemical properties: `hydro`, `polar`, `charge` (default: all three) |
| `--n-jobs` | No | Parallel workers (default: 1) |

### Output CSV Format

| Column | Description |
|--------|-------------|
| `sequence_id` | FASTA record ID (first token of header) |
| `description` | Full FASTA header |
| `pred_label` | 0 = non-enzyme, 1 = enzyme |
| `enzyme_probability` | Estimated enzyme class probability |

## Notes

1. **Sequence Format**: Input sequences should contain single-letter codes for the 20 standard amino acids. Non-standard characters are automatically stripped.
2. **Sequence Length**: Excessively short sequences (< 10 amino acids) may yield less reliable predictions.
3. **Model Format**: The pre-trained model ships as `enzyme_xgb_model.ubj` (XGBoost native UBJSON format), which loads without pickle version warnings.

## Real-World Project Example

sxlaep has been successfully applied in a large-scale marine metagenomics research project. The [Marine_Enzyme_Analysis_Pipeline_EN.ipynb](https://github.com/labxscut/sxRaep/blob/main/notebooks/Marine_Enzyme_Analysis_Pipeline_EN.ipynb) notebook demonstrates comprehensive analysis of enzyme sequences from 16,240 Marine Metagenome-Assembled Genomes (MAGs).

**Project Highlights:**
- **Dataset Scale**: 33,479,647 predicted protein sequences
- **Enzyme Identification**: 17,492,382 enzyme sequences identified
- **Analysis Scope**: Metabolic scaling law, ecological strategies, and biochemical property characterization

### 1. Parse Prediction Results and Extract Ecological Data

```python
import os
import glob
import numpy as np
import pandas as pd
from tqdm.notebook import tqdm

PRED_DIR = "/path/to/prediction_results"
csv_files = glob.glob(os.path.join(PRED_DIR, "*.csv"))
print(f"Found {len(csv_files)} MAG prediction result files.")

mag_ratios = []
mag_total_seqs = []
mag_mean_probs = []
all_probabilities = []

for f in tqdm(csv_files, desc="Parsing CSVs"):
    df = pd.read_csv(f, usecols=["pred_label", "enzyme_probability"])
    total_seqs = len(df)
    if total_seqs == 0:
        continue

    enzyme_count = (df["pred_label"] == 1).sum()
    ratio = enzyme_count / total_seqs

    if enzyme_count > 0:
        mean_prob = df.loc[df["pred_label"] == 1, "enzyme_probability"].mean()
        sampled_probs = df.loc[df["pred_label"] == 1, "enzyme_probability"].sample(frac=0.01, random_state=42).tolist()
        all_probabilities.extend(sampled_probs)
    else:
        mean_prob = 0.0

    mag_ratios.append(ratio * 100)
    mag_total_seqs.append(total_seqs)
    mag_mean_probs.append(mean_prob)

mag_ratios = np.array(mag_ratios)
mag_total_seqs = np.array(mag_total_seqs)
mag_enzyme_counts = mag_total_seqs * (mag_ratios / 100.0)
mag_mean_probs = np.array(mag_mean_probs)
all_probabilities = np.array(all_probabilities)
```

### 2. Ecological Strategy Classification

```python
categories = []
for r_val in mag_ratios:
    if r_val < 40:
        categories.append("Streamlined\n(<40%)")
    elif r_val > 60:
        categories.append("Generalist\n(>60%)")
    else:
        categories.append("Intermediate\n(40-60%)")
```

![MAG Enzyme Ratio Distribution](https://raw.githubusercontent.com/CirinMok/Picture_Raep/main/mag_enzyme_ratio_distribution.png)

*Figure: Distribution of enzyme ratios across 16,240 marine MAGs, revealing three distinct ecological strategies.*

### 3. Metabolic Scaling Law Analysis

```python
from scipy import stats

r, p_value = stats.pearsonr(mag_total_seqs, mag_enzyme_counts)
p_str = "P < 1e-10" if p_value < 1e-10 else f"P = {p_value:.2e}"

z = np.polyfit(mag_total_seqs, mag_enzyme_counts, 1)
p = np.poly1d(z)

print(f"Pearson R: {r:.3f}, {p_str}")
```

### 4. Prediction Probability Analysis

```python
import seaborn as sns
import matplotlib.pyplot as plt

sns.histplot(all_probabilities, bins=50, kde=True, color="#3C5488")
plt.axvline(np.median(all_probabilities), color='red', linestyle='dashed',
            label=f'Median: {np.median(all_probabilities):.3f}')
plt.xlabel('XGBoost Probability')
plt.title('Prediction Probability Distribution (Sampled)')
plt.legend()
plt.show()
```

### 5. Biochemical Properties Calculation (Reservoir Sampling)

```python
from Bio.SeqUtils.ProtParam import ProteinAnalysis

faa_files = glob.glob(os.path.join(PRED_DIR, "*_enzymes.faa"))
SAMPLE_RATE = 100000 / 17500000.0

sampled_seqs = []
lengths = []

def clean_seq(seq):
    return ''.join([aa for aa in seq if aa in "ACDEFGHIKLMNPQRSTVWY"])

for faa in tqdm(faa_files, desc="Sampling FASTA"):
    with open(faa, 'r') as f:
        seq_chunks = []
        for line in f:
            if line.startswith(">"):
                if seq_chunks:
                    seq = "".join(seq_chunks)
                    lengths.append(len(seq))
                    if random.random() < SAMPLE_RATE:
                        sampled_seqs.append(seq)
                seq_chunks = []
            else:
                seq_chunks.append(line.strip())

pis = []
mws = []

for seq in tqdm(sampled_seqs, desc="Computing ProtParam"):
    clean_s = clean_seq(seq)
    if not clean_s:
        continue
    pa = ProteinAnalysis(clean_s)
    try:
        pis.append(pa.isoelectric_point())
        mws.append(pa.molecular_weight() / 1000.0)  # kDa
    except Exception:
        pass
```

---

## Troubleshooting

- **Import error**: Ensure the package is installed (`pip show sxlaep`).
- **Model loading failure**: Verify the model file path exists.
- **Prediction errors**: Check that input sequences contain standard amino acid characters.
- **Performance**: For very large datasets, consider batch processing to avoid memory overflow.

## Getting Help

- **Author**: DHY
- **Email**: <dhy.scut@outlook.com>

## License

MIT License.

---

*Version: 1.0.5 | Last updated: 2026*
