Metadata-Version: 2.4
Name: pyasara
Version: 0.4.0
Summary: Pythonic interface for YASARA molecular modeling software
Author: Qiaole HE; Enzymaster; PyYASARA Contributors
License: BSD-3-Clause
Project-URL: Homepage, https://github.com/your-org/pyasara
Project-URL: Repository, https://github.com/your-org/pyasara
Keywords: molecular-modeling,bioinformatics,yasara,protein-structure,docking,md
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: BSD License
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.8
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Classifier: Topic :: Scientific/Engineering :: Chemistry
Requires-Python: >=3.8
Description-Content-Type: text/markdown

# pyasara

Pythonic interface for YASARA molecular modeling software.

## What is pyasara?

pyasara wraps YASARA's molecular modeling capabilities into a clean Python API
via a socket-based control protocol. No manual path configuration, no console
management — just `import pyasara` and start modeling.

## Requirements

- **YASARA** v25+ installed at `/opt/yasara` (or custom path via `set_path()`)
- **Python** 3.8+
- **Xvfb** (optional) for molecular visualization on headless servers

## Architecture

```
┌─────────────────┐     socket (MSG_EXECUTE/MSG_RESULT)     ┌──────────┐
│  pyasara        │ ◄─────────────────────────────────────► │  YASARA  │
│ (Python client) │                                         │(process) │
└─────────────────┘                                         └──────────┘
```

YASARA starts as a subprocess and connects back to pyasara via a TCP socket
(`-pym` protocol). All commands are sent as YASARA macro strings; results
are returned as pickled Python objects.

## Features

| Module | Capability |
|---|---|
| `homology` | Homology modeling (PSI-BLAST → multi-template → quality report) |
| `structure` | Structure preparation (CleanAll, pH, minimization, SMILES build, mutagenesis) |
| `docking` | Molecular docking (VINA protocol) |
| `md` | Molecular dynamics (force field, solvation, NVT/NPT, production → .xtc) |
| `analysis` | Structural analysis (atoms, mass, charge, surface, volume, SS, distances, angles, dihedrals, H-bonds, contacts, binding energy, Ramachandran, PCA, interaction fingerprint, pocket detection) |
| `scanning` | Mutagenesis scanning (saturation, double-mutant, parallel scanning, ΔΔG classification, substitution typing, library design & enrichment, hotspot recommendation, conservation analysis) |
| `visualization` | Molecular visualization (styles, colors, surfaces, PNG export, 3D mutation heatmap) |
| `alignment` | Structural alignment (AlignObj, SupAtom, RMSD matrix) |
| `commands` | Low-level YASARA macro wrappers (load/save PDB, SDF, Mol2, YOB, SCE) |
| `workflow` | Pipeline orchestration (enzyme engineering, virtual screening, mutation characterisation, model validation, design→MD, conformational analysis) |

## Quick Start

```python
from pyasara import run_enzyme_pipeline

result = run_enzyme_pipeline(
    structure="enzyme.pdb",
    ligand="ligand.mol2",
    scan_positions=["res 108", "res 112"],
    n_workers=4,
)
print(f"Best hotspot: {result.hotspot_rankings[0]}")
print(f"Library designs: {len(result.library_design)} positions")
print(f"Heatmap image: {result.heatmap_image}")
# Generated files: enzyme_library.csv, enzyme_hotspots.csv, enzyme_heatmap.png
```

## Examples (Function API)

### Energy Units

YASARA uses **kJ/mol** as default. To switch to kcal/mol:

```python
from pyasara.commands import _cmd
_cmd("EnergyUnit kcal/mol")
```

1 kcal = 4.184 kJ.

### Build Molecule from SMILES

```python
from pyasara import build_smiles

obj_num = build_smiles("CCO")                     # ethanol → object 1
obj_num = build_smiles("C1CCCCC1", sort=True)     # cyclohexane
```

### Structure Preparation

```python
from pyasara import prepare_structure

result = prepare_structure("raw.pdb", ph=7.0)
```

### File Format Conversion

```python
from pyasara import convert_format

# Load SDF and save as Mol2
convert_format("ligand.sdf", "ligand.mol2")

# Load Mol2 and save as PDB
convert_format("ligand.mol2", "ligand.pdb")

# Load PDB and save as YOB (YASARA binary)
convert_format("protein.pdb", "protein.yob")
```

### Molecular Visualization

```python
from pyasara import visualize

# Requires Xvfb on headless servers:
#   Xvfb :99 -screen 0 1920x1080x24 &
#   DISPLAY=:99 python your_script.py

visualize("protein.pdb", output="render.png", style="ribbon")
```

### Visualization on Headless Servers

`SavePNG` requires YASARA GUI mode and a display. On headless servers, use Xvfb:

```bash
sudo apt install xvfb
Xvfb :99 -screen 0 1920x1080x24 &
export DISPLAY=:99
python your_script.py
```

### Conservation Analysis

```python
from pyasara import compute_conservation

msa = [
    "MALWMRLLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGFFYTPKT",
    "MALWMRLLPLLALLALWAPDPAAAFVNQHLCGSHLVEALYLVCGERGFFYTPKT",
    "MALWIRLLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGFFYSPKS",
]
scores = compute_conservation(msa)  # {0: 1.0, 1: 0.78, ...}
```

### Homology Modeling

```python
from pyasara import homology_model

result = homology_model(
    sequence="MKWVTFISLLFLFSSAYS...",
    template="template.pdb",
    speed="fast",
)
print(f"Model: {result.path}  ZScore: {result.score:.2f}")
```

**Parameters:**

| Parameter | Type | Default | Description |
|-----------|------|---------|-------------|
| `sequence` | str | "" | Target protein sequence (FASTA format) |
| `template` | str | "" | Path to a single template PDB file (use OR `templates`, not both) |
| `templates` | List[str] | None | List of multiple template PDB files for multi-template modeling |
| `speed` | str | "" | Set "fast" for accelerated modeling |
| `sequence_file` | str | "" | Path to FASTA file containing the sequence |
| `output` | str | "model.pdb" | Output model filename |
| `output_dir` | str | "" | Output directory for model and report files |
| `options` | HomologyModelingOptions | None | Advanced parameters |

**Advanced options via `HomologyModelingOptions`:**

```python
from pyasara import homology_model, HomologyModelingOptions

options = HomologyModelingOptions(
    psi_blast_rounds=5,
    evalue=0.1,
    max_templates=3,
    loop_samples=10,
)
result = homology_model(
    sequence="MKWVTFISLLFLFSSAYS...",
    template="template.pdb",
    options=options,
)
```

| Option | Type | Default | Description |
|--------|------|---------|-------------|
| `psi_blast_rounds` | int | 3 | Number of PSI-BLAST iterations |
| `evalue` | float | 0.5 | E-value cutoff for sequence alignments |
| `max_templates` | int | 5 | Maximum number of templates to use |
| `max_alignments` | int | 5 | Maximum number of sequence alignments |
| `oligostate` | int | 10 | Oligomeric state for modeling |
| `term_extension` | int | 10 | Terminal extension length |
| `loop_samples` | int | 5 | Number of loop refinement samples |
| `loop_len_max` | int | 0 | Maximum loop length to refine |
| `animation` | str | "normal" | Animation mode - "normal" or "fast" |
| `struct_profile` | bool | False | Use structural profile |
| `fix_model_res` | bool | False | Fix model residues |
| `report` | bool | True | Generate HTML report |

**`HomologyModelResult` fields:**

| Field | Type | Description |
|-------|------|-------------|
| `path` | str | Path to the output PDB model file |
| `score` | float | Overall quality Z-score from YASARA (typically -2.0 to 0.0, higher is better) |
| `template` | str | Template PDB file used for modeling |
| `atom_count` | int | Number of atoms in the model |
| `sequence_identity` | float | Sequence identity percentage (aligned region only) |
| `coverage` | float | Target sequence coverage percentage |
| `log_file` | str | Path to the JSON report file |

**Quality Z-score interpretation:**
- \>= 0: Excellent quality
- -1.0 to 0: Good quality
- -2.0 to -1.0: Satisfactory quality
- < -2.0: Poor quality, significant errors expected

### Batch Homology Modeling

```python
from pyasara import homology_model_batch

result = homology_model_batch(
    sequence_file="sequences.fasta",
    template="template.pdb",
    output_dir="./models",
    speed="fast",
    n_workers=4,  # parallel modeling
)
print(f"Success: {result.success_count}, Failed: {result.failure_count}")
for r in result.results:
    print(f"  {r.path} - ZScore: {r.score:.2f}")
for name, error in result.failed:
    print(f"  FAILED: {name} - {error}")
```

**Parameters:**

| Parameter | Type | Default | Description |
|-----------|------|---------|-------------|
| `sequence_file` | str | - | **Required** Path to FASTA file |
| `template` | str | "" | Single template PDB file |
| `templates` | List[str] | None | Multiple template PDB files |
| `output_dir` | str | "" | Output directory |
| `speed` | str | "" | Set "fast" for accelerated modeling |
| `n_workers` | int | 1 | Parallel workers (1 = serial) |
| `options` | HomologyModelingOptions | None | Advanced options |

### Structural Alignment

```python
from pyasara import align_structures, compute_rmsd

# Sequence-based alignment → Cα RMSD (standard in structural biology)
result = align_structures(reference="native.pdb", mobile="model.pdb")
print(f"RMSD: {result.rmsd:.3f} Å  Aligned: {result.aligned_atoms} residues")

# Structural superposition → all-atom RMSD (more fine-grained)
rmsd = compute_rmsd("model1.pdb", "model2.pdb")
matrix = compute_rmsd_matrix(["model1.pdb", "model2.pdb", "model3.pdb"])
```

**`align_structures` vs `compute_rmsd`:**

| | `align_structures()` | `compute_rmsd()` |
|---|---|---|
| Method | `AlignObj` (sequence alignment) | `SupAtom` (structural superposition) |
| RMSD type | Cα RMSD only | All-atom RMSD (matched atoms) |
| Input matching | Sequence alignment (handles diff. sequences) | `Match=Yes` (by atom/residue/name) |
| Returns | `AlignmentResult` (RMSD + aligned count) | `float` (RMSD value only) |
| Use case | Comparing models to reference structures | Comparing similar structures/conformations |
| Compatible with | PyMOL, BioPython, ProFit (standard Cα RMSD) | More sensitive to side-chain differences |

### Molecular Dynamics

```python
from pyasara import md_simulate

result = md_simulate(
    structure="protein.pdb",
    duration=1000.0, # ps unit
)
print(f"Trajectory: {result.trajectory}")  # Returns .xtc format for GROMACS compatibility
```

**Parameters:**

| Parameter | Type | Default | Description |
|-----------|------|---------|-------------|
| `structure` | str | - | **Required** Path to input PDB file |
| `duration` | float | 1000.0 | Simulation duration in picoseconds (ps) |
| `temperature` | str | "298K" | Simulation temperature |
| `ph` | float | 7.4 | pH value for protonation |
| `ions` | str | "0.9" | Ion concentration (NaCl) |
| `force_field` | str | "AMBER14" | Force field to use |
| `cell_shape` | str | "Cube" | Simulation cell shape |
| `boundary` | str | "periodic" | Boundary conditions |
| `cutoff` | int | 8 | Cutoff distance for interactions |
| `use_gpu` | bool | True | Use GPU acceleration |
| `cpu_threads` | int | 0 | CPU threads (0 = auto-detect) |
| `speed` | str | "normal" | Simulation speed mode |
| `save_interval` | int | 100000 | Trajectory save interval |

**`MDSimulationResult` fields:**

| Field | Type | Description |
|-------|------|-------------|
| `trajectory` | str | Path to trajectory file (.xtc) |
| `structure` | str | Final structure PDB file |
| `duration` | float | Simulation duration (ps) |
| `temperature` | float | Simulation temperature (K) |
| `energy` | float | Final system energy (kJ/mol) |
| `snapshots` | int | Number of trajectory frames |

### Pocket Detection

```python
from pyasara import find_binding_pocket

# Ligand-based: residues within cutoff distance of ligand
pocket = find_binding_pocket(
    "complex.pdb",
    ligand_selection="Lig",  # LIG Lig UNK 1EK HETATM name for your ligand 
    cutoff=5.0,              # distance threshold in Å (default: 6.0)
)
print(f"Pocket residues: {len(pocket.residues)} at {pocket.center}")

# Surface-based: largest solvent-exposed concave region
pocket_surf = find_binding_pocket("protein.pdb", method="surface")
```

**Parameters:**

| Parameter | Type | Default | Description |
|-----------|------|---------|-------------|
| `path` | str | - | **Required** Structure file |
| `ligand_selection` | str | "" | YASARA ligand object name (required for `method="ligand"`) |
| `cutoff` | float | 6.0 | Distance (Å) for ligand-based detection |
| `method` | str | "auto" | `"ligand"`, `"surface"`, or `"auto"` |

### Molecular Docking

```python
from pyasara import dock

result = dock(
    receptor="protein.pdb",
    ligand="ligand.mol2",
    center=(10.0, 20.0, 30.0),
    size=(20.0, 20.0, 20.0),
)
print(f"Best pose: {result.best_pose.energy:.2f} kJ/mol")
```

**Parameters:**

| Parameter | Type | Default | Description |
|-----------|------|---------|-------------|
| `receptor` | str | - | **Required** Path to receptor PDB file |
| `ligand` | str | - | **Required** Path to ligand file (PDB/SDF/MOL2) |
| `center` | tuple | (0,0,0) | Docking grid center (x, y, z) in Å |
| `size` | tuple | (20,20,20) | Docking grid dimensions (x, y, z) in Å |
| `exhaustiveness` | int | 8 | Vina search exhaustiveness |
| `options` | DockingOptions | None | Advanced options |

**DockingOptions:**

| Option | Type | Default | Description |
|--------|------|---------|-------------|
| `num_modes` | int | 9 | Maximum number of binding modes to return |
| `work_dir` | str | "" | Working directory for intermediate files |
| `keep_files` | bool | False | Keep intermediate files after completion |

**`DockingResult` fields:**

| Field | Type | Description |
|-------|------|-------------|
| `best_pose` | DockingPose | Lowest energy binding pose |
| `all_poses` | List[DockingPose] | All binding poses |
| `receptor` | str | Receptor file path |
| `ligand` | str | Ligand file path |
| `method` | str | Docking method ("VINA") |
| `log_file` | str | Output PDBQT log file |

**`DockingPose` fields:**

| Field | Type | Description |
|-------|------|-------------|
| `energy` | float | Binding energy (kJ/mol) |
| `rmsd_lb` | float | RMSD from best pose (lower bound) |
| `rmsd_ub` | float | RMSD from best pose (upper bound) |

### Batch Docking

Virtual screening across multiple ligands:

```python
from pyasara import batch_dock, DockingOptions

opts = DockingOptions(num_modes=10, keep_files=True)
result = batch_dock(
    receptor="protein.pdb",
    ligands=["lig1.mol2", "lig2.sdf", "lig3.pdb"],
    center=(10.0, 20.0, 30.0),
    size=(20.0, 20.0, 20.0),
    exhaustiveness=8,
    n_workers=4,  # parallel docking
    options=opts,
)
print(f"Top hits:")
for h in result.hits[:5]:
    print(f"  {h.ligand_name}: {h.best_energy:+.2f} kJ/mol")
```

**Parameters:**

| Parameter | Type | Default | Description |
|-----------|------|---------|-------------|
| `receptor` | str | - | **Required** Path to receptor PDB file |
| `ligands` | List[str] | - | **Required** List of ligand file paths |
| `center` | tuple | (0,0,0) | Docking grid center (x, y, z) in Å |
| `size` | tuple | (20,20,20) | Docking grid dimensions (x, y, z) in Å |
| `exhaustiveness` | int | 8 | Vina search exhaustiveness |
| `n_workers` | int | 1 | Parallel workers (1 = serial) |
| `options` | DockingOptions | None | Advanced options |
| `verbose` | bool | True | Print progress |

**`ScreeningResult` fields:**

| Field | Type | Description |
|-------|------|-------------|
| `receptor` | str | Receptor file path |
| `hits` | List[ScreeningHit] | All hits sorted by energy |

**`ScreeningHit` fields:**

| Field | Type | Description |
|-------|------|-------------|
| `ligand_name` | str | Ligand file name |
| `ligand_path` | str | Full path to ligand file |
| `best_energy` | float | Lowest binding energy (kJ/mol) |
| `all_energies` | List[float] | Energy of all poses |
| `n_poses` | int | Number of poses found |
| `log_file` | str | Log file or error message |

### Structural Analysis

```python
from pyasara import analyze_structure

analysis = analyze_structure("protein.pdb")
print(f"Atoms: {analysis.atom_count}")
print(f"MW: {analysis.molecular_weight:.1f} Da")
print(f"Charge: {analysis.charge:.1f}")
print(f"Helix: {analysis.secondary_structure['helix']:.0%}")
print(f"H-bonds: {len(analysis.hydrogen_bonds)}")
```

**Quick measurements:**

```python
from pyasara import (
    compute_distance, compute_angle, compute_dihedral,
    compute_contacts, compute_surface, compute_hbonds,
    compute_binding_energy, check_structure, compute_ramachandran,
)

# Atom indices are 1-based (YASARA convention)
dist = compute_distance("protein.pdb", 1, 10)       # atom 1 to atom 10 distance
ang  = compute_angle("protein.pdb", 1, 2, 3)        # angle at atom 2 (atom 1-2-3)
dih  = compute_dihedral("protein.pdb", 1, 2, 3, 4)  # dihedral atom 1-2-3-4

# Contacts and surface
contacts = compute_contacts("protein.pdb", "Protein", "Ligand", cutoff=5.0)
surface  = compute_surface("protein.pdb", surface_type="accessible")
hbonds   = compute_hbonds("protein.pdb", "Protein", "Protein")

# Binding energy (requires non-periodic boundary: Boundary Type=Wall)
be = compute_binding_energy("complex.pdb", selection="Ligand")
print(f"Binding energy: {be.total:.2f} {be.unit}")

# Structure check
check = check_structure("protein.pdb", selection="all", type="Bonds")
print(f"Problems found: {check.problems}")

# Ramachandran plot
rama = compute_ramachandran("protein.pdb")
for r in rama[:5]:
    print(f"Residue {r['residue']}: φ={r['phi']:.1f}° ψ={r['psi']:.1f}°")
```

**Parameters:**

| Function | Parameters | Type | Default | Description |
|----------|------------|------|---------|-------------|
| `compute_distance` | path, atom1, atom2 | str, int, int | - | Distance between two atoms (Å) |
| `compute_angle` | path, atom1, atom2, atom3 | str, int, int, int | - | Angle at atom2 (degrees) |
| `compute_dihedral` | path, atom1-4 | str, int, int, int, int | - | Dihedral atom1-2-3-4 (degrees) |
| `compute_contacts` | path, selection1, selection2, cutoff | str, str, str, float | "Protein", "Protein", 5.0 | Residue contacts within cutoff |
| `compute_surface` | path, selection, surface_type | str, str, str | "Protein", "accessible" | Surface area; types: "accessible", "extended" |
| `compute_hbonds` | path, selection1, selection2 | str, str, str | "Protein", "Protein" | H-bond pairs |
| `compute_binding_energy` | path, selection | str, str | "Ligand" | Binding energy (kJ/mol); requires `Boundary Type=Wall` |
| `check_structure` | path, selection, type | str, str, str | "all", "Bonds" | Check types: `Bonds`, `Angles`, `Omega`, `BFactor`, `Occupancy`, `Chirality` |
| `compute_ramachandran` | path, selection | str, str | "Protein" | Phi/psi dihedrals per residue |

**Result types:**
- `BindingEnergyResult`: `total`, `vdw`, `elec`, `solv`, `unit`
- `StructureCheckResult`: `problems` (count), `details` (list of issues)
- Ramachandran returns `List[Dict]` with keys: `residue`, `phi`, `psi`

### Interaction Fingerprint

```python
from pyasara import compute_interaction_fingerprint, fingerprint_to_dataframe

fp = compute_interaction_fingerprint(
    "complex.pdb",
    receptor="Protein",  # YASARA 对象名（根据分子类型自动分配）
    ligand="Ligand",     # 如配体命名为其他名称需修改
    cutoff=5.0,
)
for f in fp[:5]:
    print(f"{f.residue_name}{f.residue}: "
          f"dist={f.min_distance:.2f}Å  "
          f"hbond={f.hbond}  hydrophobic={f.hydrophobic}  "
          f"pi_pi={f.pi_pi}  salt_bridge={f.salt_bridge}")

df = fingerprint_to_dataframe(fp)
```

**Parameters:**

| Parameter | Type | Default | Description |
|-----------|------|---------|-------------|
| `path` | str | - | **Required** Structure file with receptor + ligand |
| `receptor` | str | "Protein" | Receptor YASARA selection |
| `ligand` | str | "Ligand" | Ligand YASARA selection |
| `cutoff` | float | 5.0 | Max distance (Å) for contact detection |

**`ResidueInteraction` fields:**

| Field | Type | Description |
|-------|------|-------------|
| `residue` | str | Residue identifier (e.g. "108") |
| `residue_name` | str | Amino acid code (e.g. "ALA", "GLU") |
| `chain` | str | Chain identifier |
| `min_distance` | float | Min distance to ligand (Å) |
| `hbond` | bool | H-bond present |
| `hbond_type` | str | H-bond donor/acceptor type |
| `hydrophobic` | bool | Hydrophobic interaction |
| `pi_pi` | bool | π-π stacking |
| `salt_bridge` | bool | Salt bridge |
| `cation_pi` | bool | Cation-π interaction |

### Structure Mutation

```python
from pyasara import build_mutant, classify_ddg, classify_substitution

# Single mutation — auto-named from prefix
result = build_mutant(
    "protein.pdb",
    mutations=[("res 108", "ALA")],
)
print(f"Mutant saved to: {result.output_file}")  # → mutant_mutated.pdb

# Custom output path
result = build_mutant(
    "protein.pdb",
    mutations=[("res 108", "ALA")],
    output_path="K108A.pdb",
)

# Custom prefix (auto → {prefix}_mutated.pdb)
result = build_mutant(
    "protein.pdb",
    mutations=[("res 108", "ALA")],
    prefix="K108A",
)

# Multiple mutations (double mutant, etc.)
result = build_mutant(
    "protein.pdb",
    mutations=[("res 108", "ALA"), ("res 112", "VAL")],
    isomer="L",  # L- or D-amino acids
)

# Functional helpers (pure functions, no YASARA needed)
print(classify_ddg(-3.0))        # "stabilizing"   (≤ -2 kJ/mol)
print(classify_substitution("LEU", "ILE"))  # "conservative"
```

**`build_mutant` parameters:**

| Parameter | Type | Default | Description |
|-----------|------|---------|-------------|
| `structure_path` | str | - | **Required** Input PDB/YOB/SCE file |
| `mutations` | list | - | **Required** List of `(selection, target_aa)` pairs |
| `output_path` | str | `""` | Full output path (auto-generated from `prefix` if empty) |
| `prefix` | str | `"mutant"` | Prefix for auto-generated filename → `{prefix}_mutated.pdb` |
| `isomer` | str | `"L"` | Chirality: `"L"` or `"D"` |


## Builder Pattern

For complex workflows requiring step-by-step control, pyasara provides Builder classes:

### HomologyModeler

```python
from pyasara import HomologyModeler

model = (
    HomologyModeler()
    .set_sequence("MKWVTFISLLFLFSSAYS...")
    .set_template("template.pdb")
    .set_speed("fast")
    .set_evalue(0.1)
    .set_max_templates(3)
    .build()
)
print(f"Model: {model.path}  ZScore: {model.score:.2f}")
```

### DockingRunner

```python
from pyasara import DockingRunner, DockingOptions

opts = DockingOptions(num_modes=20, exhaustiveness=16, keep_files=True)
result = (
    DockingRunner()
    .set_receptor("protein.pdb")
    .set_ligand("ligand.mol2")
    .set_center(10.0, 20.0, 30.0)
    .set_size(20.0, 20.0, 20.0)
    .set_exhaustiveness(8)
    .set_options(opts)
    .run()
)
print(f"Best pose: {result.best_pose.energy:.2f} kJ/mol")
```

### StructureBuilder

```python
from pyasara import StructureBuilder

result = (
    StructureBuilder()
    .load("raw.pdb")
    .set_ph(7.0)
    .clean()
    .minimize(steps=1000)
    .save("prepared.pdb")
)
```

**Multiple mutations:** Call `.mutate()` chained multiple times:
```python
result = (
    StructureBuilder()
    .load("protein.pdb")
    .mutate("res 108", "ALA")
    .mutate("res 112", "VAL")
    .save("double_mutant.pdb")
)
```

### MDSimulator

```python
from pyasara import MDSimulator

traj = (
    MDSimulator()
    .set_structure("protein.pdb")
    .set_duration(1000.0)
    .set_temperature("298K")
    .set_ph(7.4)
    .set_ions("0.9")
    .set_force_field("AMBER14")
    .set_cell_shape("Cube")
    .set_use_gpu(True)
    .run()
)
print(f"Trajectory: {traj.trajectory}")
```

### MutagenesisScanner

```python
from pyasara import MutagenesisScanner, HotspotRecommender, visualize_scanning_results

scanner = MutagenesisScanner("protein.pdb")

# Single-point scan
results = scanner.run_scan(
    positions=["res 108", "res 112"],
    residues=["ALA", "VAL", "LEU"],
)
```

#### ΔΔG Classification

```python
from pyasara import classify_ddg, classify_substitution

# Single mutation
results = scanner.run_scan(positions=["res 108"], residues=["ALA", "VAL"])
for r in MutagenesisScanner.to_dicts(results):
    print(f"{r['position']} {r['wildtype']}→{r['mutant']}: "
          f"ΔΔG={r['delta_energy']:+.2f}  {r['stability_class']}")

# Double mutant (epistasis calculation)
double = scanner.run_double_mutant_scan(
    positions=["res 108", "res 112"],
    residues=["ALA", "VAL", "SER"],  # amino acids to substitute at both positions
)
for d in double:
    print(f"Epistasis: {d.epistasis:+.2f}  ({d.epistasis_class})")
```

**`SingleMutationResult` fields:** `position`, `wildtype`, `mutant`, `delta_energy`, `stability_class`, `substitution_type`

**`DoubleMutationResult` fields:** All SingleMutationResult fields plus `epistasis` (ε = ΔΔG_double − ΔΔG_single1 − ΔΔG_single2), `epistasis_class` ("synergistic"/"additive"/"antagonistic")

#### Library Design

```python
results = scanner.run_scan(...)
designs = scanner.design_library(results, cutoff=-1.0)

for d in designs:
    print(f"Position {d.position} ({d.wildtype}):")
    print(f"  Best mutants: {[(r.mutant, round(r.delta_energy, 2)) for r in d.best_mutants(3)]}")

csv_path = MutagenesisScanner.library_to_csv(designs, "library.csv")
enrichment = MutagenesisScanner.compute_library_enrichment(designs)
```

#### Parallel Scanning

```python
# Single-point — 8 workers
results = scanner.run_scan_parallel(n_workers=8)

# Double-mutant — 4 workers
results2 = scanner.run_double_mutant_scan_parallel(
    positions=["res 108", "res 112"],
    n_workers=4,
)

import pandas as pd
df = pd.DataFrame(MutagenesisScanner.to_dicts(results2))
```

**Note:** Each worker reloads the PDB independently and runs in a separate YASARA instance.

#### 3D Mutation Heatmap

```python
visualize_scanning_results(
    "protein.pdb",
    results,
    output="heatmap.png",
    ligand="Ligand",
)
# Blue = stabilising, Red = destabilising (white = neutral)
```

#### Hotspot Recommendation

```python
from pyasara import compute_interaction_fingerprint

fp = compute_interaction_fingerprint("complex.pdb", ligand="Ligand")
recommender = HotspotRecommender("complex.pdb")
rankings = recommender.recommend(scan_results=results, fingerprint=fp, top_n=10)

df = HotspotRecommender.rankings_to_dataframe(rankings)
print(df[["rank", "position", "combined_score"]])
```

### Visualizer

```python
from pyasara import Visualizer

viz = (
    Visualizer()
    .load("protein.pdb")
    .style("ribbon")
    .color_by("secondary_structure")
    .export_image("render.png")
)
```

## Recommended Workflows

pyasara's modules can be combined into powerful end‑to‑end workflows.

### Enzyme Engineering

Design and prioritise enzyme mutants via docking, fingerprint, and saturation scanning.

```python
result = run_enzyme_pipeline("enzyme.pdb", ligand_selection="substrate.mol2")
```

Steps: prepare → dock → fingerprint → parallel scan → library design → 3D heatmap.

### Virtual Screening

Dock a library of ligands against a target and rank hits by binding energy.

```python
result = run_virtual_screening_pipeline(
    "target.pdb",
    ligands=["lig1.mol2", "lig2.mol2", "lig3.sdf"],
    detect_pocket=True,
    n_workers=4,
)
df = result.screening.to_dataframe()
print(df.head(10))
```

### Mutation Characterisation

Compare interaction fingerprints before and after a point mutation.

```python
result = run_mutation_characterization_pipeline(
    "complex.pdb",
    mutation_selection="res 108",
    mutation_target="ALA",
    ligand_selection="Ligand",
)
comp = result.fingerprint_comparison
print(f"Gained: {comp.gained}, Lost: {comp.lost}")
```

### Model Validation

Assess a homology model against its template by RMSD and Ramachandran quality.

```python
result = run_model_validation_pipeline(
    "model.pdb",
    template_path="template.pdb",
    ligand_selection="Ligand",
)
print(f"RMSD: {result.alignment_rmsd:.3f} Å")
print(f"Ramachandran outliers: {result.ramachandran_outliers}")
```

### Design → MD Validation

Scan hotspots, build the best mutant, and run short MD to check stability.

```python
result = run_design_md_pipeline(
    "enzyme.pdb",
    scan_positions=["res 108", "res 112"],
    duration=5000.0,
)
print(f"Top mutant: {result.top_mutant_path}")
```

### Conformational Analysis

Generate a trajectory via MD and extract principal motions with PCA.

```python
result = run_conformational_analysis_pipeline(
    "protein.pdb",
    duration=10000.0,
)
pca = result.pca_result
print(f"PCA modes: {pca.modes}, top eigenvalue: {pca.eigenvalues[0]:.1f}")
```

## Feature Details

### PCA / NMA

```python
from pyasara import compute_pca

result = compute_pca("trajectory.xtc", selection="Protein")
print(f"Modes: {result.modes}, Top eigenvalue: {result.eigenvalues[0]:.1f}")
```

**Parameters:**

| Parameter | Type | Default | Description |
|-----------|------|---------|-------------|
| `path` | str | - | Trajectory file (.xtc) |
| `selection` | str | "Protein" | YASARA selection |
| `mode` | str | "Atom" | PCA mode type |

**Note:** Requires trajectory data and YASARA Dynamics license.

### YASARA Command Format Notes

YASARA's raw macro parser accepts both comma-separated and keyword-arg formats:

- **Comma-separated**: `ColorAtom all,Element`
- **Keyword-arg**: `BuildSMILES String=CCO,Sort=Yes`
- **Returns lists**: Most commands return `[value]` (list), not bare values

## Tests

```bash
# Unit tests (requires YASARA installed)
pytest tests/test_basic.py -v

# Integration tests (E2E with real YASARA socket)
pytest tests/test_integration.py -v

# All tests
pytest tests/ -v
```

## License

BSD-3-Clause
