1app.py

"""

Bioinformatics Protein Analysis Application

Optimized for Hugging Face Spaces & Docker Deployment

"""

import os
import io
import csv
import sys
import subprocess
import pandas as pd
import joblib
import numpy as np
import torch
import esm
import logging
from pathlib import Path
from datetime import datetime
from flask import Flask, request, render_template, send_file, session, redirect, url_for
from Bio.Seq import Seq

# ========== LOGGING SETUP ==========
logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s %(levelname)s %(message)s",
    handlers=[
        logging.FileHandler("app.log"),
        logging.StreamHandler(sys.stdout)
    ],
)
logger = logging.getLogger(__name__)

# ========== ENVIRONMENT & PATHS ==========
# Get base directory (works in Docker and local)
BASE_DIR = Path(os.getenv("BASE_DIR", "/app"))
DATA_DIR = BASE_DIR / "data"
TEMP_DIR = BASE_DIR / "temp"
MODEL_DIR = BASE_DIR / "models"
TEMPLATE_DIR = BASE_DIR / "templates"
STATIC_DIR = BASE_DIR / "static"

# Create directories if they don't exist
for dir_path in [DATA_DIR, TEMP_DIR, MODEL_DIR]:
    dir_path.mkdir(parents=True, exist_ok=True)

# Model paths - use defaults if not found
MODEL_PATH = MODEL_DIR / "best_model_LR.sav"
ENCODER_PATH = MODEL_DIR / "label_encoder.pkl"

# BLAST configuration
BLAST_PATH = os.getenv("BLAST_PATH", "/usr/bin/blastp")  # Linux default
BLAST_DB = DATA_DIR / "pathway_db"
MAPPING_FILE = DATA_DIR / "pathway_map.csv"

# Thresholds
BITSCORE_THRESHOLD = 80.0
CONFIDENCE_THRESHOLD = 0.7
RESULT_CSV = TEMP_DIR / "combined_results.csv"

# ========== MODEL LOADING ==========
model = None
encoder = None
esm_model = None
batch_converter = None

logger.info("===== Application Startup =====")
logger.info(f"BASE_DIR: {BASE_DIR}")
logger.info(f"DATA_DIR: {DATA_DIR}")
logger.info(f"TEMP_DIR: {TEMP_DIR}")

# Load Logistic Regression Model
try:
    if MODEL_PATH.exists():
        model = joblib.load(MODEL_PATH)
        logger.info("✅ LR model loaded successfully.")
    else:
        logger.warning(f"⚠️ Model not found at {MODEL_PATH}")
except Exception as e:
    logger.error(f"❌ Failed to load LR model: {e}")

# Load Label Encoder
try:
    if ENCODER_PATH.exists():
        encoder = joblib.load(ENCODER_PATH)
        logger.info("✅ Label encoder loaded successfully.")
    else:
        logger.warning(f"⚠️ Encoder not found at {ENCODER_PATH}")
except Exception as e:
    logger.error(f"❌ Failed to load label encoder: {e}")

# Load ESM2 Model (downloads if not cached)
try:
    logger.info("Loading ESM2 model (this may take a moment on first run)...")
    esm_model, alphabet = esm.pretrained.load_model_and_alphabet("esm2_t6_8M_UR50D")
    esm_model.eval()
    esm_model = esm_model.to("cpu")
    batch_converter = alphabet.get_batch_converter()
    logger.info("✅ ESM2 model loaded successfully on CPU.")
except Exception as e:
    logger.error(f"❌ Failed to load ESM2 model: {e}")

# ========== HELPER FUNCTIONS ==========

def get_best_frame_protein(nuc_seq):
    """Translate DNA to best reading frame protein"""
    nuc_seq = "".join(nuc_seq.split()).upper()
    trim_len = len(nuc_seq) - (len(nuc_seq) % 3)
    nuc_seq = nuc_seq[:trim_len]
    
    seq = Seq(nuc_seq)
    frames = [
        seq[i:].translate(to_stop=False) for i in range(3)
    ] + [
        seq.reverse_complement()[i:].translate(to_stop=False) for i in range(3)
    ]
    
    best_protein = ""
    for prot in frames:
        longest_orf = max(prot.split("*"), key=len)
        if len(longest_orf) > len(best_protein):
            best_protein = longest_orf
    
    return str(best_protein)


def detect_and_translate(seq):
    """Detect if DNA or protein and handle accordingly"""
    seq = seq.upper()
    
    # Check if DNA
    if all(base in "ATGCN" for base in seq):
        logger.info("Detected DNA sequence → translating to protein.")
        if len(seq) < 150:
            logger.warning("DNA sequence too short for translation")
            return None
        
        prot_seq = get_best_frame_protein(seq)
        if len(prot_seq) < 50:
            logger.warning("Translated protein too short")
            return None
        return prot_seq
    else:
        logger.info("Detected protein sequence → using directly.")
        if len(seq) < 50:
            logger.warning("Protein sequence too short")
            return None
        return seq


def esm2_320_embed(sequence):
    """Generate ESM2 embeddings for sequence"""
    if esm_model is None or batch_converter is None:
        raise RuntimeError("ESM model not available.")
    
    batch_labels, batch_strs, batch_tokens = batch_converter([("seq1", sequence)])
    
    with torch.no_grad():
        results = esm_model(batch_tokens, repr_layers=[6], return_contacts=False)
    
    token_representations = results["representations"][6]
    seq_repr = token_representations[0, 1:-1].detach().cpu().numpy()
    
    return seq_repr.mean(axis=0)


def run_blast_and_get_dataframe(temp_fasta, blast_output):
    """Execute BLAST and return results as DataFrame"""
    try:
        # Check if BLAST is available
        if not Path(BLAST_PATH).exists():
            logger.warning(f"⚠️ BLAST not found at {BLAST_PATH}, skipping BLAST search")
            return pd.DataFrame(
                columns=["qseqid", "sseqid", "pident", "length", "evalue", "bitscore"]
            )
        
        # Check if database exists
        if not BLAST_DB.exists():
            logger.warning(f"⚠️ BLAST database not found at {BLAST_DB}, skipping BLAST search")
            return pd.DataFrame(
                columns=["qseqid", "sseqid", "pident", "length", "evalue", "bitscore"]
            )
        
        cmd = [
            str(BLAST_PATH),
            "-query", str(temp_fasta),
            "-db", str(BLAST_DB),
            "-out", str(blast_output),
            "-outfmt", "6 qseqid sseqid pident length evalue bitscore"
        ]
        
        logger.info(f"🧬 Running BLAST...")
        proc = subprocess.run(cmd, capture_output=True, text=True, timeout=300)
        
        if proc.returncode != 0:
            logger.error(f"❌ BLAST error: {proc.stderr}")
            return pd.DataFrame(
                columns=["qseqid", "sseqid", "pident", "length", "evalue", "bitscore"]
            )
        
        if not Path(blast_output).exists() or Path(blast_output).stat().st_size == 0:
            logger.warning("⚠️ No BLAST results generated.")
            return pd.DataFrame(
                columns=["qseqid", "sseqid", "pident", "length", "evalue", "bitscore"]
            )
        
        cols = ["qseqid", "sseqid", "pident", "length", "evalue", "bitscore"]
        df = pd.read_csv(blast_output, sep="\t", names=cols)
        df["bitscore"] = df["bitscore"].astype(float)
        logger.info(f"✅ BLAST completed successfully with {len(df)} hits")
        
        return df
    
    except subprocess.TimeoutExpired:
        logger.error("❌ BLAST timed out (>300s)")
        return pd.DataFrame(
            columns=["qseqid", "sseqid", "pident", "length", "evalue", "bitscore"]
        )
    except Exception as e:
        logger.error(f"❌ Unexpected error in BLAST: {e}")
        return pd.DataFrame(
            columns=["qseqid", "sseqid", "pident", "length", "evalue", "bitscore"]
        )


# ========== FLASK APPLICATION ==========

app = Flask(
    __name__,
    template_folder=str(TEMPLATE_DIR),
    static_folder=str(STATIC_DIR)
)
app.secret_key = os.getenv("SECRET_KEY", "default-secret-key-change-in-production")
app.config['ENV'] = os.getenv("FLASK_ENV", "production")
app.config['DEBUG'] = os.getenv("FLASK_DEBUG", "False") == "True"

latest_predictions = []


@app.after_request
def add_cache_control(response):
    """Prevent caching of responses"""
    response.headers["Cache-Control"] = "no-cache, no-store, must-revalidate"
    response.headers["Pragma"] = "no-cache"
    response.headers["Expires"] = "0"
    return response


@app.route("/", methods=["GET", "POST"])
def index():
    """Main index page - handle sequence input and analysis"""
    global latest_predictions
    latest_predictions = []
    
    if request.method == "POST":
        sequences = []
        headers = []
        
        # Parse uploaded FASTA file
        uploaded_file = request.files.get("fasta_file")
        if uploaded_file and uploaded_file.filename != "":
            seq = ""
            header = ""
            for line in uploaded_file:
                line = line.decode().strip()
                if not line:
                    continue
                if line.startswith(">"):
                    if seq:
                        sequences.append(seq)
                        headers.append(header)
                        seq = ""
                    header = line[1:]
                else:
                    seq += line
            if seq:
                sequences.append(seq)
                headers.append(header)
        
        # Parse textarea input
        sequence_text = request.form.get("sequence_text", "").strip()
        if sequence_text:
            seq = ""
            header = ""
            for line in sequence_text.splitlines():
                line = line.strip()
                if not line:
                    continue
                if line.startswith(">"):
                    if seq:
                        sequences.append(seq)
                        headers.append(header)
                        seq = ""
                    header = line[1:]
                else:
                    seq += line
            if seq:
                sequences.append(seq)
                headers.append(header or "sequence_from_text")
        
        if not sequences:
            return render_template(
                "index.html",
                error="⚠️ No valid sequences provided.",
                current_year=datetime.now().year
            )
        
        # Detect and translate sequences
        translated_sequences = [detect_and_translate(s) for s in sequences]
        filtered_sequences = [s for s in translated_sequences if s is not None]
        filtered_headers = [
            h for s, h in zip(translated_sequences, headers) if s is not None
        ]
        
        if not filtered_sequences:
            return render_template(
                "index.html",
                error="⚠️ All sequences too short or invalid.",
                current_year=datetime.now().year
            )
        
        # Gene Family Prediction using ESM2 embeddings
        gene_preds = []
        try:
            logger.info(f"Generating embeddings for {len(filtered_sequences)} sequences...")
            features = np.array([esm2_320_embed(seq) for seq in filtered_sequences])
            
            if model and encoder:
                probs = model.predict_proba(features)
                max_probs = probs.max(axis=1)
                pred_indices = probs.argmax(axis=1)
                
                for idx, prob in zip(pred_indices, max_probs):
                    if prob >= CONFIDENCE_THRESHOLD:
                        gene_preds.append(encoder.inverse_transform([idx])[0])
                    else:
                        gene_preds.append("Unknown")
            else:
                logger.warning("Model or encoder not available, marking all as Unknown")
                gene_preds = ["Unknown"] * len(filtered_sequences)
        
        except Exception as e:
            logger.error(f"Gene family prediction failed: {e}")
            gene_preds = ["Unknown"] * len(filtered_sequences)
        
        # Save temporary FASTA for BLAST
        temp_fasta = TEMP_DIR / "temp_input.fasta"
        with open(temp_fasta, "w") as f:
            for h, s in zip(filtered_headers, filtered_sequences):
                f.write(f">{h.replace(' ', '_')}\n{s}\n")
        
        # Run BLAST
        blast_output = TEMP_DIR / "blast_results.txt"
        blast_df = run_blast_and_get_dataframe(temp_fasta, blast_output)
        
        # Load pathway mapping file
        map_df = pd.DataFrame()
        if MAPPING_FILE.exists():
            try:
                map_df = pd.read_csv(MAPPING_FILE, dtype=str)
                logger.info(f"Loaded {len(map_df)} pathway mappings")
            except Exception as e:
                logger.error(f"Failed to read mapping file: {e}")
        
        # Build results
        combined_results = []
        
        if not blast_df.empty:
            top_hits_idx = blast_df.groupby("qseqid")["bitscore"].idxmax()
            top_hits = blast_df.loc[top_hits_idx]
        else:
            top_hits = pd.DataFrame(
                columns=["qseqid", "sseqid", "pident", "length", "evalue", "bitscore"]
            )
        
        for i, orig_header in enumerate(filtered_headers):
            safe_h = orig_header.replace(" ", "_")
            top_hit_row = top_hits[top_hits["qseqid"] == safe_h]
            
            if not top_hit_row.empty:
                row = top_hit_row.iloc[0]
                top_hit = row["sseqid"]
                bitscore = float(row["bitscore"])
                pathways = "No Pathways Found"
                
                if (
                    not map_df.empty
                    and "Entry" in map_df.columns
                    and "Pathways" in map_df.columns
                ):
                    match = map_df.loc[map_df["Entry"] == top_hit, "Pathways"]
                    if not match.empty and bitscore >= BITSCORE_THRESHOLD:
                        pathways = match.values[0]
            else:
                top_hit = "No Hit"
                bitscore = 0.0
                pathways = "No Pathways Found"
            
            gene_family = gene_preds[i] if i < len(gene_preds) else "Unknown"
            combined_results.append([orig_header, top_hit, bitscore, pathways, gene_family])
        
        # Save results to CSV
        result_df = pd.DataFrame(
            combined_results,
            columns=["Query/Header", "Top Hit", "Bitscore", "Pathways", "Predicted Gene Family"]
        )
        result_df.to_csv(RESULT_CSV, index=False)
        latest_predictions = combined_results
        
        # Store in session
        session['predictions'] = combined_results
        session['result_summary'] = {
            'total_sequences': len(filtered_sequences),
            'successful_predictions': len([p for p in gene_preds if p != "Unknown"]),
            'total_pathways_found': len([r for r in combined_results if r[3] != "No Pathways Found"])
        }
        
        return redirect(url_for('results'))
    
    return render_template("index.html", predictions=[], current_year=datetime.now().year)


@app.route("/results")
def results():
    """Display analysis results"""
    predictions = session.get("predictions", [])
    summary = session.get("result_summary", {})
    
    if not predictions:
        return redirect(url_for("index"))
    
    return render_template(
        "results.html",
        predictions=predictions,
        summary=summary,
        current_year=datetime.now().year
    )


@app.route("/download_csv")
def download_csv():
    """Download results as CSV"""
    global latest_predictions
    
    if not latest_predictions:
        return "⚠️ No predictions to download."
    
    output = io.StringIO()
    writer = csv.writer(output)
    writer.writerow(["Query/Header", "Top Hit", "Bitscore", "Pathways", "Predicted Gene Family"])
    writer.writerows(latest_predictions)
    output.seek(0)
    
    return send_file(
        io.BytesIO(output.getvalue().encode()),
        mimetype="text/csv",
        as_attachment=True,
        download_name="bioinformatics_results.csv"
    )


@app.route("/health")
def health():
    """Health check endpoint"""
    return {
        "status": "healthy",
        "esm_model": esm_model is not None,
        "lr_model": model is not None,
        "encoder": encoder is not None
    }


# ========== ERROR HANDLERS ==========

@app.errorhandler(404)
def not_found(error):
    return render_template("error.html", error="Page not found"), 404


@app.errorhandler(500)
def internal_error(error):
    logger.error(f"Internal error: {error}")
    return render_template("error.html", error="Internal server error"), 500


# ========== MAIN ==========

if __name__ == "__main__":
    # Get port from environment or default to 7860 (HF Spaces)
    port = int(os.getenv("PORT", 7860))
    host = os.getenv("HOST", "0.0.0.0")
    debug = os.getenv("FLASK_DEBUG", "False") == "True"
    
    logger.info(f"Starting Flask app on {host}:{port} (debug={debug})")
    
    app.run(
        host=host,
        port=port,
        debug=debug,
        threaded=True,
        use_reloader=False
    )