Metadata-Version: 2.4
Name: raep
Version: 1.0.1
Summary: Random Forest Enzyme Prediction
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.7
Description-Content-Type: text/markdown
Requires-Dist: numpy
Requires-Dist: pandas
Requires-Dist: scikit-learn
Requires-Dist: xgboost
Requires-Dist: joblib

# RAEP: Rapid Enzyme/Non-Enzyme Prediction

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

**RAEP (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**: Achieves fast and accurate enzyme/non-enzyme classification using optimized feature extraction and the XGBoost model.
- **Multi-Mode Support**: Supports single-sequence prediction, multi-sequence batch prediction, and FASTA file batch prediction.
- **Rich Feature Set**: Utilizes multi-physicochemical property pseudo-amino acid composition (Pseudo-AAC), CTD features, and windowed amino acid composition.
- **User-Friendly**: Offers a concise Python API that is easy to integrate into existing projects.
- **Multi-Process Optimization**: Employs multi-process parallel processing in the feature extraction step to improve processing efficiency for large-scale datasets.

## 📦 Dependencies

- `joblib`: Used for model saving/loading and parallel processing.
- `numpy`: Numerical computation.
- `pandas`: Data processing.
- `scikit-learn`: Machine learning utilities and evaluation metrics.
- `xgboost`: Implementation of the gradient boosting tree algorithm.

  **Requirements**: Python 3.7 or higher.

## 📥 Installation

### Install from PyPI (Recommended)

```bash
pip install raep
```

### Quick start with CLI

```bash
raep --input /your_fasta/file.fasta --output /path_to_your_result/result.json
```

## 💻 Basic Usage

### Import and Initialization

```python
from raep import RAEP

# Default initialization (uses built-in model)
predictor = RAEP()

# Initialization with a custom model path (optional)
# predictor = RAEP(model_path="path/to/your/model.pkl")
```

### Single Sequence Prediction

```python
# Predict a single protein sequence
sequence = "MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHGKKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTPAVHASLDKFLASVSTVLTSKYR"

prediction = predictor.predict(sequence)
probability = predictor.predict_proba(sequence)
non_enzyme_prob = 1.0 - enzyme_prob

print(f"Prediction result: {'Enzyme' if prediction == 1 else 'Non-Enzyme'}")
print(f"Prediction probabilities: Non-Enzyme={non_enzyme_prob:.4f}, Enzyme={enzyme_prob:.4f}")
```

### Batch Prediction from FASTA Files

```python
# Batch predict from a FASTA file
fasta_path = "test_sequences.fasta"
results = predictor.predict_fasta(fasta_path)

print(f"Prediction results ({len(results)} sequences):")
for protein_id, result in results.items():
    pred = result['prediction']
    prob = result['probability']
    print(f"{protein_id}: {'Enzyme' if pred == 1 else 'Non-Enzyme'} (Enzyme probability: {prob:.4f})")
```

## 📚 API Reference

### `RAEP` Class

#### Initialization

`RAEP(model_path=None)`

- **Purpose**: Instantiates the RAEP predictor, automatically loads the model and initializes feature extraction parameters (e.g., `LAG=10`, `W=0.05`), ensuring consistency in subsequent prediction workflows.
- **Parameters**:
  - `model_path`: Optional. Path to a custom model file. If not provided, the built-in `enzyme_xgb_model.pkl` model will be used.

#### Methods

**`predict(sequence)`**

- **Purpose**: Predicts whether a single protein sequence is an enzyme.
- **Parameters**:
  - `sequence` (String): The protein sequence to be predicted.
- **Returns**:
  - `prediction` (Int): 0 = Non-enzyme, 1 = Enzyme.

**`predict_proba(sequence)`**

- **Purpose**: Predicts whether a single protein sequence is an enzyme.
- **Parameters**:
  - `sequence` (String): The protein sequence to be predicted.
- **Returns**:
  - `probability` (float): Probability of the sequence being an enzyme.

**`predict_fasta(fasta_path)`**

- **Purpose**: Performs batch prediction for protein sequences from a FASTA file.
- **Parameters**:
  - `fasta_path` (String): Path to the FASTA file.
- **Returns**: Dictionary where keys are protein IDs (from FASTA headers) and values are dictionaries with `prediction` (0 for non-enzyme, 1 for enzyme) and `probability` (float, enzyme probability).

## 📝 Notes

1. **Sequence Format Requirements**: Input sequences should only contain single-letter codes (uppercase) for the 20 standard amino acids.
2. **Sequence Length**: The tool automatically processes sequences of different lengths, but excessively short sequences (e.g., < 10 amino acids) may affect prediction accuracy.
3. **Multi-Process Processing**: The feature extraction process uses multi-processing acceleration by default, which automatically adjusts based on the number of CPU cores in the system.
4. **Model File**: Ensure the model file exists and is accessible, especially when using a custom model path.

## 🌟 Real-World Project Example

RAEP 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 how to use RAEP for 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=["PredLabel", "Prob_Enzyme"])
    total_seqs = len(df)
    if total_seqs == 0:
        continue
    
    is_enzyme = (df["PredLabel"] == "Enzyme")
    enzyme_count = is_enzyme.sum()
    ratio = enzyme_count / total_seqs
    
    if enzyme_count > 0:
        mean_prob = df.loc[is_enzyme, "Prob_Enzyme"].mean()
        sampled_probs = df.loc[is_enzyme, "Prob_Enzyme"].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
# Classify MAGs based on enzyme ratio
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: Streamlined (<40%), Intermediate (40-60%), and Generalist (>60%) genomes.*

### 3. Metabolic Scaling Law Analysis

```python
from scipy import stats

# Calculate Pearson correlation between genome size and enzyme count
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}"

# Linear regression
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

# Plot probability distribution
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()
```

![Prediction Probability Distribution](https://raw.githubusercontent.com/CirinMok/Picture_Raep/main/probability_distribution.png)

*Figure: Distribution of enzyme prediction probabilities, demonstrating the model's confidence in classification.*

![MAG Confidence Distribution](https://raw.githubusercontent.com/CirinMok/Picture_Raep/main/mag_confidence_distribution.png)

*Figure: Distribution of mean prediction confidence across MAGs.*

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

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

faa_files = glob.glob(os.path.join(PRED_DIR, "*_enzymes.faa"))
SAMPLE_RATE = 100000 / 17500000.0  # Sample ~100k sequences

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
```

![Sequence Length Distribution](https://raw.githubusercontent.com/CirinMok/Picture_Raep/main/length_distribution.png)

*Figure: Length distribution of predicted enzyme sequences from marine metagenomes.*

![Molecular Weight Distribution](https://raw.githubusercontent.com/CirinMok/Picture_Raep/main/mw_distribution.png)

*Figure: Molecular weight distribution of predicted enzyme sequences (kDa).*

![Isoelectric Point Distribution](https://raw.githubusercontent.com/CirinMok/Picture_Raep/main/pi_distribution.png)

*Figure: Isoelectric point (pI) distribution of predicted enzyme sequences.*

### 6. Visualization of Biochemical Properties

```python
fig, axes = plt.subplots(2, 2, figsize=(15, 10))
sns.despine()

# Length Distribution
lengths_arr = np.array(lengths)
p99_len = np.percentile(lengths_arr, 99)
sns.histplot(lengths_arr[lengths_arr <= p99_len], bins=80, kde=True, color="#00A087", ax=axes[0, 0])
axes[0, 0].axvline(np.median(lengths_arr), color='red', linestyle='dashed', 
                   label=f'Median: {np.median(lengths_arr):.0f}')
axes[0, 0].set_title('Enzyme Sequence Length Distribution')
axes[0, 0].legend()

# pI Distribution
sns.histplot(pis, bins=50, kde=True, color="#E64B35", ax=axes[1, 0])
axes[1, 0].set_title('Isoelectric Point (pI) Distribution')

# Molecular Weight Distribution
p99_mw = np.percentile(mws, 99)
sns.histplot([m for m in mws if m <= p99_mw], bins=50, kde=True, color="#8491B4", ax=axes[1, 1])
axes[1, 1].set_title('Molecular Weight Distribution')
axes[1, 1].set_xlabel('Molecular Weight (kDa)')

plt.tight_layout()
plt.show()
```

![Amino Acid Composition](https://raw.githubusercontent.com/CirinMok/Picture_Raep/main/aa_composition.png)

*Figure: Amino acid composition analysis of marine enzyme sequences.*

---

The complete analysis pipeline includes:
1. **Batch prediction** using `predict_fasta()` for large-scale sequence analysis
2. **Ecological strategy classification** based on enzyme ratios (Streamlined/Intermediate/Generalist)
3. **Metabolic scaling law analysis** (genome size vs enzyme count correlation)
4. **Prediction confidence evaluation** for quality assessment
5. **Biochemical property calculation** (pI, molecular weight, amino acid composition)
6. **Publication-quality figure generation** following Science/Nature journal standards

## 🔧 Troubleshooting

- **Failed to import the RAEP package**: Ensure the package is correctly installed in the current Python environment (`pip show raep`).
- **Model loading failure**: Verify that the model file path is correct and the file exists at the specified location.
- **Prediction errors**: Check if the input sequence format is valid and contains only standard amino acid characters.
- **Performance issues**: For extremely large datasets, consider processing in batches to avoid memory overflow.

## 🤝 Getting Help

If you encounter any problems, please contact the author:

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

## 📄 License

This project is licensed under the MIT License. See the `LICENSE` file for details.

## 🙏 Acknowledgments

The development of this project is supported by several open-source tools, especially machine learning libraries such as **XGBoost** and **scikit-learn**.

***

*Version: 1.0.1 | Last updated: 2026*

