batchProtParam.py

#!/usr/bin/env python3
"""
Batch ProtParam Analysis Script (Biopython Version)
---------------------------------------------------

This script calculates protein properties (ProtParam-style) for all FASTA
files in a folder. It is designed for batch processing of many proteins
(100s–100,000s) without using a web browser.

WHAT THIS SCRIPT CALCULATES (per protein sequence):

- Amino acid counts (20 standard amino acids)
- Amino acid percentages
- Molecular weight
- Aromaticity
- Theoretical isoelectric point (pI)
- Secondary structure fraction (helix, turn, sheet)
- GRAVY (hydrophobicity score)
- Instability index
- Flexibility (mean, min, max, stdev)

The results are written as CSV files that open directly in Excel.


---------------------------------------------------
HOW TO RUN THE SCRIPT
---------------------------------------------------

This guide assumes Windows + PowerShell.

Step 1) Put files in the right place

Create a folder with this structure:

    project_folder/
        batchProtParam.py
        fastas/
            proteins1.fasta
            proteins2.fasta

All FASTA files you want to process should be inside the "fastas" folder.

Step 2) Open PowerShell in project_folder

In File Explorer:
    - Open project_folder
    - Click the address bar
    - Type: powershell
    - Press Enter

Step 3) Run one of these commands

Option A (default): One CSV per FASTA file

    python batchProtParam.py --in_dir .\fastas --out_dir .\results --output_mode per_fasta

What you get:
    - A "results" folder
    - One CSV output file for each FASTA input file

Option B: One CSV for all FASTA files together

    python batchProtParam.py --in_dir .\fastas --out_dir .\results_combined --output_mode all_fastas

What you get:
    - A "results" folder
    - One combined CSV file containing all sequences from all FASTA files

Optional: choose the combined CSV filename

    python batchProtParam.py --in_dir .\fastas --out_dir .\results_combined --output_mode all_fastas --all_fastas_name my_results.csv


---------------------------------------------------
AMBIGUOUS AMINO ACIDS (IMPORTANT)
---------------------------------------------------

Some sequences contain letters that are NOT one of the 20 standard amino acids.

Examples:
    X = unknown amino acid
    B = D or N
    Z = E or Q
    U = selenocysteine
    O = pyrrolysine

You must choose how to handle these using --ambiguous:

Default (recommended for most datasets):

    --ambiguous drop

    Removes non-standard letters before calculation.

Strict mode (best for publication-quality work):

    --ambiguous fail

    Skips any sequence containing non-standard letters.

Advanced mode:

    --ambiguous keep

    Keeps all letters exactly as-is.
    May cause errors depending on sequence content.

If unsure, use the default (drop).

"""

from __future__ import annotations

import argparse
import csv
import statistics
from pathlib import Path
from typing import Dict, Iterable, List, Optional, Tuple

from Bio import SeqIO
from Bio.SeqUtils.ProtParam import ProteinAnalysis

# Common FASTA extensions (protein FASTA often uses .faa/.fasta/.fa/.fas)
FASTA_EXTS = {".fa", ".fasta", ".faa", ".fsa", ".fas", ".fna"}

AA20 = list("ACDEFGHIKLMNPQRSTVWY")

FIELDNAMES = (
    ["seq_id", "length_aa"]
    + [f"count_{aa}" for aa in AA20]
    + [f"pct_{aa}" for aa in AA20]
    + [
        "molecular_weight",
        "aromaticity",
        "theoretical_pi",
        "ss_helix",
        "ss_turn",
        "ss_sheet",
        "gravy",
        "instability_index",
        "flex_mean",
        "flex_min",
        "flex_max",
        "flex_stdev",
        "source_fasta",
        "warnings",
        "status",
        "error_type",
    ]
)


def iter_fasta_files(in_dir: Path) -> Iterable[Path]:
    for p in in_dir.iterdir():
        if p.is_file() and p.suffix.lower() in FASTA_EXTS:
            yield p


def clean_seq(seq: str, mode: str) -> Tuple[str, Optional[str]]:
    """
    Returns (cleaned_seq, warning_message).

    mode:
      - "drop": remove any non-standard residues
      - "fail": if any non-standard residues exist, return empty seq with warning
      - "keep": keep them (may cause ProteinAnalysis methods to error)
    """
    s = str(seq).upper().replace("*", "")  # remove stop symbol if present
    allowed = set(AA20)
    nonstd = sorted(set([aa for aa in s if aa not in allowed]))

    if mode == "drop":
        cleaned = "".join([aa for aa in s if aa in allowed])
        warn = f"dropped_nonstandard={''.join(nonstd)}" if nonstd else None
        return cleaned, warn

    if mode == "fail":
        if nonstd:
            return "", f"nonstandard_residues={''.join(nonstd)}"
        return s, None

    if mode == "keep":
        warn = f"kept_nonstandard={''.join(nonstd)}" if nonstd else None
        return s, warn

    raise ValueError(f"Unknown mode: {mode}")


def protparam_metrics(seq_id: str, seq: str) -> Dict[str, object]:
    """
    Compute ProtParam-like metrics using Biopython ProteinAnalysis.
    Returns a dict with keys matching a subset of FIELDNAMES.
    """
    pa = ProteinAnalysis(seq)

    counts = pa.count_amino_acids()         # dict AA -> count
    # Biopython: get_amino_acids_percent() is deprecated; use attribute when available.
    perc = getattr(pa, "amino_acids_percent", None)
    if perc is None:
        perc = pa.get_amino_acids_percent()  # fallback for older Biopython

    ss_h, ss_t, ss_s = pa.secondary_structure_fraction()
    flex = pa.flexibility()                # list[float]

    # Flexibility summary stats (for CSV friendliness)
    if flex:
        flex_mean = statistics.mean(flex)
        flex_min = min(flex)
        flex_max = max(flex)
        flex_stdev = statistics.pstdev(flex) if len(flex) > 1 else 0.0
    else:
        flex_mean = flex_min = flex_max = float("nan")
        flex_stdev = float("nan")

    row: Dict[str, object] = {
        "seq_id": seq_id,
        "length_aa": len(seq),
        "molecular_weight": pa.molecular_weight(),
        "aromaticity": pa.aromaticity(),
        "theoretical_pi": pa.isoelectric_point(),
        "ss_helix": ss_h,
        "ss_turn": ss_t,
        "ss_sheet": ss_s,
        "gravy": pa.gravy(),
        "instability_index": pa.instability_index(),
        "flex_mean": flex_mean,
        "flex_min": flex_min,
        "flex_max": flex_max,
        "flex_stdev": flex_stdev,
    }

    # Force all 20 AAs to exist in output columns
    for aa in AA20:
        row[f"count_{aa}"] = int(counts.get(aa, 0))
        row[f"pct_{aa}"] = float(perc.get(aa, 0.0)) * 100.0  # percent

    return row


def write_csv(path: Path, rows: List[Dict[str, object]]) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    with path.open("w", encoding="utf-8", newline="") as f:
        writer = csv.DictWriter(f, fieldnames=FIELDNAMES)
        writer.writeheader()
        for r in rows:
            # Ensure every row has all fields (DictWriter will fill missing with "")
            writer.writerow(r)


def main() -> None:
    ap = argparse.ArgumentParser(description="Batch ProtParam analysis for FASTA files (directory input).")
    ap.add_argument("--in_dir", required=True, help="Directory containing FASTA files.")
    ap.add_argument("--out_dir", required=True, help="Directory to write CSV outputs.")
    ap.add_argument(
        "--output_mode",
        choices=["per_fasta", "all_fastas"],
        default="per_fasta",
        help="Output behavior: one CSV per FASTA ('per_fasta') or one CSV for all FASTAs ('all_fastas').",
    )
    ap.add_argument(
        "--all_fastas_name",
        default="protparam_all.csv",
        help="Filename to use when --output_mode all_fastas (default: protparam_all.csv).",
    )
    ap.add_argument(
        "--ambiguous",
        choices=["drop", "fail", "keep"],
        default="drop",
        help="How to handle non-standard amino acids (B, Z, X, U, O, etc.). Default: drop.",
    )
    ap.add_argument(
        "--format",
        default="fasta",
        help="SeqIO format (default: fasta). Change only if you know inputs are another format.",
    )
    args = ap.parse_args()

    in_dir = Path(args.in_dir).expanduser().resolve()
    out_dir = Path(args.out_dir).expanduser().resolve()

    if not in_dir.exists() or not in_dir.is_dir():
        raise SystemExit(f"Input directory not found or not a directory: {in_dir}")

    # Create output directory up front
    out_dir.mkdir(parents=True, exist_ok=True)

    fasta_files = sorted(iter_fasta_files(in_dir))
    if not fasta_files:
        raise SystemExit(f"No FASTA files found in {in_dir} (extensions: {sorted(FASTA_EXTS)})")

    all_rows: List[Dict[str, object]] = []

    for fasta_path in fasta_files:
        rows: List[Dict[str, object]] = []
        base = fasta_path.stem
        out_path = out_dir / f"{base}.protparam.csv"

        for rec in SeqIO.parse(str(fasta_path), args.format):
            seq_id = rec.id
            raw_seq = str(rec.seq)

            cleaned, warn = clean_seq(raw_seq, args.ambiguous)
            warnings = warn or ""

            if not cleaned:
                # skipped due to ambiguous mode=fail or empty after dropping
                row = {k: "" for k in FIELDNAMES}
                row.update(
                    {
                        "seq_id": seq_id,
                        "length_aa": len(cleaned),
                        "source_fasta": fasta_path.name,
                        "warnings": warnings,
                        "status": "skipped",
                        "error_type": "",
                    }
                )
                rows.append(row)
                continue

            try:
                m = protparam_metrics(seq_id, cleaned)
                m.update(
                    {
                        "source_fasta": fasta_path.name,
                        "warnings": warnings,
                        "status": "ok",
                        "error_type": "",
                    }
                )
                rows.append(m)
            except Exception as e:
                # Keep going even if one sequence fails
                row = {k: "" for k in FIELDNAMES}
                row.update(
                    {
                        "seq_id": seq_id,
                        "length_aa": len(cleaned),
                        "source_fasta": fasta_path.name,
                        "warnings": (warnings + (";" if warnings else "") + f"error={type(e).__name__}"),
                        "status": "error",
                        "error_type": type(e).__name__,
                    }
                )
                rows.append(row)

        all_rows.extend(rows)

        if args.output_mode == "per_fasta":
            write_csv(out_path, rows)
            print(f"Wrote: {out_path}  (n={len(rows)})")

    if args.output_mode == "all_fastas":
        combined_path = out_dir / args.all_fastas_name
        write_csv(combined_path, all_rows)
        print(f"Wrote: {combined_path}  (n={len(all_rows)})")


if __name__ == "__main__":
    main()