app.py

"""
RNA-seq Processing App

A Streamlit-based GUI for processing RNA-seq data using Salmon.
"""

# Raise file descriptor limit before any heavy imports
# macOS app bundles launched via Finder inherit a 256 soft limit from launchd,
# which is too low for Streamlit + the scientific Python stack.
import resource
try:
    _soft, _hard = resource.getrlimit(resource.RLIMIT_NOFILE)
    resource.setrlimit(resource.RLIMIT_NOFILE, (min(_hard, 8192), _hard))
except (ValueError, resource.error):
    pass

import streamlit as st
import pandas as pd
import plotly.express as px
import plotly.graph_objects as go
import numpy as np
import re
from pathlib import Path
import os
from datetime import datetime
import subprocess
import json
import time
import io
import threading
from PIL import Image


def _start_idle_watchdog(timeout=60):
    """Shut down Streamlit if no browser is connected for `timeout` seconds."""
    def watchdog():
        idle_since = None
        while True:
            time.sleep(10)
            try:
                # Check for active WebSocket connections on port 8501
                result = subprocess.run(
                    ["lsof", "-i", ":8501", "-sTCP:ESTABLISHED"],
                    capture_output=True, text=True
                )
                has_clients = bool(result.stdout.strip())
            except Exception:
                has_clients = True  # assume connected if check fails

            if has_clients:
                idle_since = None
            else:
                if idle_since is None:
                    idle_since = time.time()
                elif time.time() - idle_since > timeout:
                    os._exit(0)

    t = threading.Thread(target=watchdog, daemon=True)
    t.start()

# Start the idle watchdog once per process (not on every Streamlit rerun)
if not getattr(_start_idle_watchdog, "_started", False):
    _start_idle_watchdog(timeout=60)
    _start_idle_watchdog._started = True


from pipeline import (
    SampleInfo,
    PipelineConfig,
    check_salmon_installed,
    get_salmon_version,
    download_reference_index,
    build_index_from_fasta,
    validate_fasta_is_transcriptome,
    run_pipeline,
    merge_counts,
    get_mapping_stats,
    validate_fastq
)

# App configuration
st.set_page_config(
    page_title="GenomifySeq",
    page_icon="🧬",
    layout="wide"
)

# Initialize session state
if "pipeline_running" not in st.session_state:
    st.session_state.pipeline_running = False
if "results" not in st.session_state:
    st.session_state.results = None
if "log_messages" not in st.session_state:
    st.session_state.log_messages = []
if "detected_samples" not in st.session_state:
    st.session_state.detected_samples = None
if "metadata_df" not in st.session_state:
    st.session_state.metadata_df = None
if "analysis_results" not in st.session_state:
    st.session_state.analysis_results = None
if "de_results" not in st.session_state:
    st.session_state.de_results = None
if "gsea_results" not in st.session_state:
    st.session_state.gsea_results = None
if "gmt_file_path" not in st.session_state:
    st.session_state.gmt_file_path = None
if "volcano_overlay" not in st.session_state:
    st.session_state.volcano_overlay = None
if "geneset_overlay_genes" not in st.session_state:
    st.session_state.geneset_overlay_genes = None
if "geneset_overlay_name" not in st.session_state:
    st.session_state.geneset_overlay_name = None
if "volcano_overlay_r" not in st.session_state:
    st.session_state.volcano_overlay_r = None

# Paths - save user data to Documents folder (not inside app bundle)
APP_DIR = Path(__file__).parent
USER_DATA_DIR = Path.home() / "Documents" / "GenomifySeq"
INDEX_DIR = USER_DATA_DIR / "indices"
RESULTS_DIR = USER_DATA_DIR / "results"
TX2GENE_DIR = USER_DATA_DIR / "tx2gene"
SESSIONS_DIR = USER_DATA_DIR / "sessions"

# Ensure directories exist
for d in [INDEX_DIR, RESULTS_DIR, TX2GENE_DIR, SESSIONS_DIR]:
    d.mkdir(parents=True, exist_ok=True)


def detect_species():
    """Auto-detect species from the selected index path name."""
    index_path = st.session_state.get("index_path")
    if index_path and "mouse" in str(index_path).lower():
        return "mouse"
    return "human"


def add_log(message: str):
    """Add a message to the log."""
    timestamp = datetime.now().strftime("%H:%M:%S")
    st.session_state.log_messages.append(f"[{timestamp}] {message}")


def png_to_tiff(png_path: Path) -> bytes:
    """Convert a PNG file to TIFF format and return as bytes."""
    img = Image.open(png_path)
    # Convert to RGB if necessary (TIFF doesn't support all PNG modes well)
    if img.mode in ('RGBA', 'LA', 'P'):
        # Create white background for transparent images
        background = Image.new('RGB', img.size, (255, 255, 255))
        if img.mode == 'P':
            img = img.convert('RGBA')
        background.paste(img, mask=img.split()[-1] if img.mode == 'RGBA' else None)
        img = background
    elif img.mode != 'RGB':
        img = img.convert('RGB')

    buffer = io.BytesIO()
    img.save(buffer, format='TIFF', compression='tiff_lzw', dpi=(300, 300))
    buffer.seek(0)
    return buffer.getvalue()


def plotly_to_tiff(fig: go.Figure, width: int = 1200, height: int = 900) -> bytes:
    """Convert a Plotly figure to TIFF format and return as bytes."""
    # First convert to PNG using Plotly's to_image
    try:
        png_bytes = fig.to_image(format="png", width=width, height=height, scale=2)
        img = Image.open(io.BytesIO(png_bytes))

        # Convert to RGB
        if img.mode != 'RGB':
            img = img.convert('RGB')

        buffer = io.BytesIO()
        img.save(buffer, format='TIFF', compression='tiff_lzw', dpi=(300, 300))
        buffer.seek(0)
        return buffer.getvalue()
    except Exception as e:
        # Fallback: kaleido may not be installed
        raise RuntimeError(f"Could not export to TIFF. Install kaleido: pip install kaleido. Error: {e}")


def format_time(seconds: float) -> str:
    """Format seconds into human-readable time string."""
    if seconds < 60:
        return f"{int(seconds)}s"
    elif seconds < 3600:
        mins = int(seconds // 60)
        secs = int(seconds % 60)
        return f"{mins}m {secs}s"
    else:
        hours = int(seconds // 3600)
        mins = int((seconds % 3600) // 60)
        return f"{hours}h {mins}m"


def serialize_samples(detected_samples: dict) -> dict:
    """Serialize detected_samples dict for JSON storage (convert Path to str)."""
    if not detected_samples:
        return None
    result = {
        "read_type": detected_samples.get("read_type"),
        "samples": []
    }
    for s in detected_samples.get("samples", []):
        sample_data = {
            "name": s["name"],
            "r1": str(s["r1"]) if s.get("r1") else None,
            "r2": str(s["r2"]) if s.get("r2") else None,
        }
        result["samples"].append(sample_data)
    return result


def deserialize_samples(data: dict) -> dict:
    """Deserialize samples dict from JSON (convert str back to Path)."""
    if not data:
        return None
    result = {
        "read_type": data.get("read_type"),
        "samples": []
    }
    for s in data.get("samples", []):
        sample_data = {
            "name": s["name"],
            "r1": Path(s["r1"]) if s.get("r1") else None,
            "r2": Path(s["r2"]) if s.get("r2") else None,
        }
        result["samples"].append(sample_data)
    return result


def save_session(name: str) -> Path:
    """Save current session to JSON + CSV files."""
    session_dir = SESSIONS_DIR / f"{name}_{datetime.now().strftime('%Y%m%d_%H%M%S')}"
    session_dir.mkdir(parents=True, exist_ok=True)

    # Save metadata_df as CSV
    if st.session_state.metadata_df is not None:
        st.session_state.metadata_df.to_csv(session_dir / "metadata.csv", index=False)

    # Save everything else as JSON
    session_data = {
        "data_dir_path": st.session_state.get("data_dir_path"),
        "detected_samples": serialize_samples(st.session_state.get("detected_samples")),
        "index_path": str(st.session_state.get("index_path", "")) if st.session_state.get("index_path") else None,
        "results": st.session_state.get("results"),
        "analysis_results": st.session_state.get("analysis_results"),
        "de_results": st.session_state.get("de_results"),
        "gsea_results": st.session_state.get("gsea_results"),
        "gmt_file_path": st.session_state.get("gmt_file_path"),
    }
    (session_dir / "session.json").write_text(json.dumps(session_data, indent=2, default=str))
    return session_dir


def load_session(session_dir: Path):
    """Restore session from saved files."""
    # Load JSON
    session_data = json.loads((session_dir / "session.json").read_text())

    # Restore to session state
    st.session_state.data_dir_path = session_data.get("data_dir_path")
    st.session_state.detected_samples = deserialize_samples(session_data.get("detected_samples"))
    if session_data.get("index_path"):
        st.session_state.index_path = Path(session_data["index_path"])
    st.session_state.results = session_data.get("results")
    st.session_state.analysis_results = session_data.get("analysis_results")
    st.session_state.de_results = session_data.get("de_results")
    st.session_state.gsea_results = session_data.get("gsea_results")
    st.session_state.gmt_file_path = session_data.get("gmt_file_path")

    # Load metadata CSV - ensure all columns are strings for data_editor compatibility
    metadata_file = session_dir / "metadata.csv"
    if metadata_file.exists():
        metadata_df = pd.read_csv(metadata_file, dtype=str)
        # Replace NaN with empty strings
        metadata_df = metadata_df.fillna("")
        st.session_state.metadata_df = metadata_df


def create_interactive_volcano(
    df: pd.DataFrame,
    lfc_threshold: float = 1.0,
    padj_threshold: float = 0.05,
    search_genes: list = None,
    label_genes: list = None,
    highlight_genes: list = None,
    title: str = "Volcano Plot"
) -> go.Figure:
    """
    Create an interactive volcano plot with Plotly.

    Parameters:
    -----------
    df : DataFrame with columns 'gene', 'log2FoldChange', 'padj' (and optionally 'pvalue', 'baseMean')
    lfc_threshold : Log2 fold change threshold for significance
    padj_threshold : Adjusted p-value threshold for significance
    search_genes : List of genes to highlight from search
    label_genes : List of genes to label on the plot
    highlight_genes : List of genes to highlight (e.g., from gene set)
    title : Plot title

    Returns:
    --------
    Plotly Figure object
    """
    # Make a copy to avoid modifying original
    plot_df = df.copy()

    # Calculate -log10(padj) for y-axis
    plot_df["neg_log10_padj"] = -np.log10(plot_df["padj"].clip(lower=1e-300))

    # Determine significance category
    def get_significance(row):
        if pd.isna(row["padj"]) or pd.isna(row["log2FoldChange"]):
            return "NS"
        if row["padj"] < padj_threshold and row["log2FoldChange"] > lfc_threshold:
            return "Up"
        elif row["padj"] < padj_threshold and row["log2FoldChange"] < -lfc_threshold:
            return "Down"
        else:
            return "NS"

    plot_df["category"] = plot_df.apply(get_significance, axis=1)

    # Color mapping (blue=up, red=down)
    color_map = {"Up": "#3498db", "Down": "#e74c3c", "NS": "#95a5a6"}

    # Create base figure
    fig = go.Figure()

    # Add scatter traces for each category (for legend)
    for category, color in color_map.items():
        mask = plot_df["category"] == category
        subset = plot_df[mask]

        if len(subset) == 0:
            continue

        # Check if these genes should be highlighted
        if search_genes:
            search_mask = subset["gene"].isin(search_genes)
            # Non-searched genes (smaller, more transparent)
            non_search = subset[~search_mask]
            if len(non_search) > 0:
                fig.add_trace(go.Scatter(
                    x=non_search["log2FoldChange"],
                    y=non_search["neg_log10_padj"],
                    mode="markers",
                    marker=dict(color=color, size=6, opacity=0.3),
                    name=f"{category} ({len(non_search)})",
                    text=non_search["gene"],
                    customdata=np.stack([
                        non_search["log2FoldChange"].round(3),
                        non_search["padj"].apply(lambda x: f"{x:.2e}"),
                        non_search.get("baseMean", pd.Series(["N/A"] * len(non_search))),
                    ], axis=-1) if "baseMean" in non_search.columns else np.stack([
                        non_search["log2FoldChange"].round(3),
                        non_search["padj"].apply(lambda x: f"{x:.2e}"),
                    ], axis=-1),
                    hovertemplate="<b>%{text}</b><br>" +
                                  "log2FC: %{customdata[0]}<br>" +
                                  "padj: %{customdata[1]}<br>" +
                                  "<extra></extra>",
                    showlegend=True,
                    legendgroup=category,
                ))
            # Searched genes (larger, fully opaque)
            search_subset = subset[search_mask]
            if len(search_subset) > 0:
                fig.add_trace(go.Scatter(
                    x=search_subset["log2FoldChange"],
                    y=search_subset["neg_log10_padj"],
                    mode="markers",
                    marker=dict(color=color, size=12, opacity=1.0,
                               line=dict(color="black", width=2)),
                    name=f"Search match ({len(search_subset)})",
                    text=search_subset["gene"],
                    customdata=np.stack([
                        search_subset["log2FoldChange"].round(3),
                        search_subset["padj"].apply(lambda x: f"{x:.2e}"),
                    ], axis=-1),
                    hovertemplate="<b>%{text}</b><br>" +
                                  "log2FC: %{customdata[0]}<br>" +
                                  "padj: %{customdata[1]}<br>" +
                                  "<extra></extra>",
                    showlegend=True,
                    legendgroup="search",
                ))
        else:
            # Normal display without search
            fig.add_trace(go.Scatter(
                x=subset["log2FoldChange"],
                y=subset["neg_log10_padj"],
                mode="markers",
                marker=dict(color=color, size=6, opacity=0.7),
                name=f"{category} ({len(subset)})",
                text=subset["gene"],
                customdata=np.stack([
                    subset["log2FoldChange"].round(3),
                    subset["padj"].apply(lambda x: f"{x:.2e}"),
                ], axis=-1),
                hovertemplate="<b>%{text}</b><br>" +
                              "log2FC: %{customdata[0]}<br>" +
                              "padj: %{customdata[1]}<br>" +
                              "<extra></extra>",
                showlegend=True,
                legendgroup=category,
            ))

    # Highlight significant genes from gene set overlay (matching R behavior)
    if highlight_genes:
        highlight_mask = (plot_df["gene"].isin(highlight_genes)) & (plot_df["category"] != "NS")
        highlight_subset = plot_df[highlight_mask]
        if len(highlight_subset) > 0:
            # Add black outline circles around significant genes in the set
            fig.add_trace(go.Scatter(
                x=highlight_subset["log2FoldChange"],
                y=highlight_subset["neg_log10_padj"],
                mode="markers",
                marker=dict(
                    color="rgba(0,0,0,0)",  # Transparent fill
                    size=7,  # Slightly larger than points (size 6) for visibility
                    line=dict(color="black", width=1.5),
                    symbol="circle"
                ),
                name=f"Gene set sig. ({len(highlight_subset)})",
                text=highlight_subset["gene"],
                customdata=np.stack([
                    highlight_subset["log2FoldChange"].round(3),
                    highlight_subset["padj"].apply(lambda x: f"{x:.2e}"),
                ], axis=-1),
                hovertemplate="<b>%{text}</b><br>" +
                              "log2FC: %{customdata[0]}<br>" +
                              "padj: %{customdata[1]}<br>" +
                              "<extra></extra>",
                showlegend=True,
            ))

            # Add labels for top 10 significant genes from the gene set (sorted by padj)
            top_highlight = highlight_subset.nsmallest(10, "padj")
            if len(top_highlight) > 0:
                fig.add_trace(go.Scatter(
                    x=top_highlight["log2FoldChange"],
                    y=top_highlight["neg_log10_padj"],
                    mode="text",
                    text=top_highlight["gene"],
                    textposition="top center",
                    textfont=dict(size=10, color="black"),
                    name=f"Labels ({len(top_highlight)})",
                    hoverinfo="skip",
                    showlegend=False,
                ))

    # Add gene labels
    if label_genes:
        label_mask = plot_df["gene"].isin(label_genes)
        label_subset = plot_df[label_mask]
        if len(label_subset) > 0:
            # Add labeled points with larger markers
            fig.add_trace(go.Scatter(
                x=label_subset["log2FoldChange"],
                y=label_subset["neg_log10_padj"],
                mode="markers+text",
                marker=dict(
                    color="purple",
                    size=10,
                    opacity=1.0,
                    line=dict(color="white", width=2),
                ),
                text=label_subset["gene"],
                textposition="top center",
                textfont=dict(size=11, color="black"),
                name=f"Labeled ({len(label_subset)})",
                customdata=np.stack([
                    label_subset["log2FoldChange"].round(3),
                    label_subset["padj"].apply(lambda x: f"{x:.2e}"),
                ], axis=-1),
                hovertemplate="<b>%{text}</b><br>" +
                              "log2FC: %{customdata[0]}<br>" +
                              "padj: %{customdata[1]}<br>" +
                              "<extra></extra>",
                showlegend=True,
            ))

    # Add threshold lines
    neg_log10_padj_threshold = -np.log10(padj_threshold)

    # Horizontal line for p-value threshold
    fig.add_hline(y=neg_log10_padj_threshold, line_dash="dash",
                  line_color="gray", opacity=0.7,
                  annotation_text=f"padj = {padj_threshold}",
                  annotation_position="bottom right")

    # Vertical lines for fold change thresholds
    fig.add_vline(x=lfc_threshold, line_dash="dash",
                  line_color="gray", opacity=0.7)
    fig.add_vline(x=-lfc_threshold, line_dash="dash",
                  line_color="gray", opacity=0.7,
                  annotation_text=f"|log2FC| = {lfc_threshold}",
                  annotation_position="bottom left")

    # Update layout
    fig.update_layout(
        title=dict(text=title, x=0.5, xanchor="center"),
        xaxis_title="log2 Fold Change",
        yaxis_title="-log10(adjusted p-value)",
        hovermode="closest",
        legend=dict(
            yanchor="top",
            y=0.99,
            xanchor="left",
            x=1.02,
            bgcolor="rgba(255,255,255,0.8)"
        ),
        plot_bgcolor="white",
        width=800,
        height=600,
    )

    # Add gridlines
    fig.update_xaxes(showgrid=True, gridwidth=1, gridcolor="lightgray",
                     zeroline=True, zerolinewidth=1, zerolinecolor="gray")
    fig.update_yaxes(showgrid=True, gridwidth=1, gridcolor="lightgray",
                     zeroline=True, zerolinewidth=1, zerolinecolor="gray")

    return fig


def open_folder_dialog(initial_dir: str = None) -> str:
    """Open native macOS folder picker dialog using AppleScript."""
    initial = initial_dir or str(Path.home())

    script = f'''
    set initialPath to POSIX file "{initial}" as alias
    try
        set selectedFolder to choose folder with prompt "Select FASTQ Directory" default location initialPath
        return POSIX path of selectedFolder
    on error
        return ""
    end try
    '''

    try:
        result = subprocess.run(
            ['osascript', '-e', script],
            capture_output=True,
            text=True,
            timeout=120
        )
        return result.stdout.strip()
    except Exception:
        return ""


def open_file_dialog(initial_dir: str = None) -> str:
    """Open native macOS file picker dialog using AppleScript."""
    initial = initial_dir or str(Path.home())

    script = f'''
    set initialPath to POSIX file "{initial}" as alias
    try
        set selectedFile to choose file with prompt "Select FASTA File" default location initialPath of type {{"fa", "fasta", "gz", "public.data"}}
        return POSIX path of selectedFile
    on error
        return ""
    end try
    '''

    try:
        result = subprocess.run(
            ['osascript', '-e', script],
            capture_output=True,
            text=True,
            timeout=120
        )
        return result.stdout.strip()
    except Exception:
        return ""


def open_gmt_dialog(initial_dir: str = None) -> str:
    """Open native macOS file picker for GMT files."""
    initial = initial_dir or str(Path.home())

    script = f'''
    set initialPath to POSIX file "{initial}" as alias
    try
        set selectedFile to choose file with prompt "Select GMT Gene Set File" default location initialPath
        return POSIX path of selectedFile
    on error
        return ""
    end try
    '''

    try:
        result = subprocess.run(
            ['osascript', '-e', script],
            capture_output=True,
            text=True,
            timeout=120
        )
        return result.stdout.strip()
    except Exception:
        return ""


def get_basename(filename: str) -> str:
    """Strip FASTQ extensions to get base name."""
    name = filename
    # Remove .gz first if present
    if name.endswith('.gz'):
        name = name[:-3]
    # Remove .fastq or .fq
    if name.endswith('.fastq'):
        name = name[:-6]
    elif name.endswith('.fq'):
        name = name[:-3]
    return name


def detect_read_pair(filename: str) -> tuple[str, int]:
    """
    Detect if a file is Read 1 or Read 2 based on END of filename.

    Returns (sample_base_name, read_number) where read_number is 1, 2, or 0 (unknown).

    Supports patterns at end of basename:
    - _1, _2 (e.g., sample_1.fastq.gz, sample_2.fastq.gz)
    - _R1, _R2 (e.g., sample_R1.fastq.gz, sample_R2.fastq.gz)
    - .R1, .R2 (e.g., sample.R1.fastq.gz)
    - _read1, _read2 (e.g., sample_read1.fastq.gz)
    """
    basename = get_basename(filename)

    # Patterns to check at END of basename (order matters - check more specific first)
    # Each tuple: (regex pattern, read number, replacement to get sample name)
    end_patterns = [
        (r'[_.]read1$', 1),
        (r'[_.]read2$', 2),
        (r'[_.]R1$', 1),
        (r'[_.]R2$', 2),
        (r'[_.]1$', 1),
        (r'[_.]2$', 2),
    ]

    for pattern, read_num in end_patterns:
        match = re.search(pattern, basename, re.IGNORECASE)
        if match:
            sample_name = basename[:match.start()]
            return (sample_name, read_num)

    # No paired pattern found
    return (basename, 0)


def scan_fastq_directory(directory: Path) -> dict:
    """
    Scan a directory for FASTQ files and detect paired-end samples.

    Returns dict with 'samples' list and 'read_type' (paired/single).
    """
    fastq_extensions = {'.fastq', '.fq', '.fastq.gz', '.fq.gz'}

    # Find all FASTQ files
    fastq_files = []
    for f in directory.iterdir():
        if f.is_file():
            # Check extension (handle .fastq.gz)
            if f.suffix == '.gz':
                ext = ''.join(f.suffixes[-2:]) if len(f.suffixes) >= 2 else f.suffix
            else:
                ext = f.suffix
            if ext.lower() in fastq_extensions:
                fastq_files.append(f)

    if not fastq_files:
        return {"samples": [], "read_type": None, "error": "No FASTQ files found"}

    # Detect read pairs based on END of filename
    file_info = []
    for f in fastq_files:
        sample_name, read_num = detect_read_pair(f.name)
        file_info.append({
            "file": f,
            "sample_name": sample_name,
            "read_num": read_num
        })

    # Group by sample name
    samples_dict = {}
    for info in file_info:
        name = info["sample_name"]
        if name not in samples_dict:
            samples_dict[name] = {"r1": None, "r2": None, "unpaired": []}

        if info["read_num"] == 1:
            samples_dict[name]["r1"] = info["file"]
        elif info["read_num"] == 2:
            samples_dict[name]["r2"] = info["file"]
        else:
            samples_dict[name]["unpaired"].append(info["file"])

    # Determine if we have valid paired-end samples
    paired_samples = []
    single_files = []

    for name, files in samples_dict.items():
        if files["r1"] and files["r2"]:
            # Valid pair
            paired_samples.append({
                "name": name,
                "r1": files["r1"],
                "r2": files["r2"]
            })
        else:
            # Collect as single-end
            if files["r1"]:
                single_files.append(files["r1"])
            if files["r2"]:
                single_files.append(files["r2"])
            single_files.extend(files["unpaired"])

    # Return paired if we have pairs, otherwise single
    if paired_samples:
        paired_samples.sort(key=lambda x: x["name"])
        return {"samples": paired_samples, "read_type": "paired"}
    elif single_files:
        single_files.sort(key=lambda x: x.name)
        samples = []
        for f in single_files:
            sample_name = get_basename(f.name)
            samples.append({
                "name": sample_name,
                "r1": f,
                "r2": None
            })
        return {"samples": samples, "read_type": "single"}

    return {"samples": [], "read_type": None, "error": "Could not determine read type"}


def main():
    st.title("🧬 GenomifySeq")
    st.markdown("*End-to-end RNA-seq analysis: quantification, differential expression, and gene set enrichment*")

    # Sidebar - System Status
    with st.sidebar:
        st.header("System Status")

        # Check Salmon
        salmon_ok = check_salmon_installed()
        if salmon_ok:
            st.success(f"Salmon: {get_salmon_version()}")
        else:
            st.error("Salmon not found!")
            st.markdown("""
            Please run the setup script first:
            ```bash
            ./setup.sh
            ```
            Then launch with:
            ```bash
            ./run_app.sh
            ```
            """)
            return

        st.divider()

        # Available indices
        st.subheader("Available Indices")
        existing_indices = [d for d in INDEX_DIR.iterdir() if d.is_dir()] if INDEX_DIR.exists() else []
        if existing_indices:
            for idx in existing_indices:
                st.text(f"• {idx.name}")
        else:
            st.text("No indices downloaded yet")

        st.divider()

        # Settings
        st.subheader("Settings")
        threads = st.slider("CPU threads", 1, os.cpu_count() or 4, min(4, os.cpu_count() or 4))

        st.divider()

        # About / Methods section for publications
        st.subheader("About / Methods")
        with st.expander("View methods for publication"):
            st.markdown("""
**RNA-seq Quantification and Analysis Methods**

*Copy and adapt this text for your Methods section:*

---

**Transcript Quantification**

Raw RNA-seq reads were quantified using Salmon v1.10 (Patro et al., 2017) with the `--validateMappings` flag for improved accuracy. Reads were mapped to the [human/mouse] GENCODE reference transcriptome (version [X]). Salmon employs a quasi-mapping approach that provides fast and accurate transcript-level abundance estimates.

**Gene-Level Summarization**

Transcript-level counts were summarized to gene-level using tximport (Soneson et al., 2015) with the `lengthScaledTPM` method, which accounts for changes in average transcript length across samples.

**Normalization and Quality Control**

Count data were normalized using the median-of-ratios method implemented in DESeq2 (Love et al., 2014). Principal component analysis (PCA) was performed on variance-stabilized transformed (VST) counts to assess sample clustering and identify potential outliers.

**Differential Expression Analysis**

Differential gene expression was assessed using DESeq2 with default parameters. Genes were considered significantly differentially expressed with an adjusted p-value (Benjamini-Hochberg) < 0.05 and |log2 fold change| > 1.0.

**Gene Set Enrichment Analysis**

Gene set enrichment analysis was performed using fgsea (Korotkevich et al., 2021) with genes pre-ranked by the Wald statistic from DESeq2. Gene sets were obtained from MSigDB (Liberzon et al., 2015).

---

**References**

- Patro R, et al. (2017) Salmon provides fast and bias-aware quantification of transcript expression. *Nature Methods* 14:417-419.
- Soneson C, et al. (2015) Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. *F1000Research* 4:1521.
- Love MI, et al. (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. *Genome Biology* 15:550.
- Korotkevich G, et al. (2021) Fast gene set enrichment analysis. *bioRxiv* doi:10.1101/060012.
- Liberzon A, et al. (2015) The Molecular Signatures Database Hallmark Gene Set Collection. *Cell Systems* 1:417-425.
            """)

            # Add copy button for the methods text
            methods_text = """Transcript Quantification: Raw RNA-seq reads were quantified using Salmon v1.10 (Patro et al., 2017) with the --validateMappings flag. Reads were mapped to the GENCODE reference transcriptome.

Gene-Level Summarization: Transcript-level counts were summarized to gene-level using tximport (Soneson et al., 2015) with the lengthScaledTPM method.

Normalization and Quality Control: Count data were normalized using DESeq2 (Love et al., 2014). PCA was performed on variance-stabilized transformed counts.

Differential Expression: Differential expression was assessed using DESeq2. Genes with adjusted p-value < 0.05 and |log2FC| > 1.0 were considered significant.

Gene Set Enrichment: GSEA was performed using fgsea (Korotkevich et al., 2021) with genes ranked by DESeq2 Wald statistic."""

            st.download_button(
                "Download Methods Text",
                methods_text,
                file_name="rnaseq_methods.txt",
                mime="text/plain"
            )

        st.divider()

        # Uninstall
        with st.expander("Uninstall"):
            st.markdown(
                "To fully remove GenomifySeq, delete these two locations:\n\n"
                "1. **The app** — drag *GenomifySeq.app* from Applications to Trash\n"
                "2. **Dependencies (~2 GB)** — delete the hidden folder `~/.genomifyseq`\n"
                "3. **Your data** *(optional)* — delete `~/Documents/GenomifySeq`"
            )
            if st.button("Reveal dependency folder in Finder", use_container_width=True):
                dep_dir = Path.home() / ".genomifyseq"
                dep_dir.mkdir(exist_ok=True)
                subprocess.Popen(["open", str(dep_dir)])
            if st.button("Reveal data folder in Finder", use_container_width=True):
                USER_DATA_DIR.mkdir(parents=True, exist_ok=True)
                subprocess.Popen(["open", str(USER_DATA_DIR)])

        st.divider()

        if st.button("Quit App", type="secondary", use_container_width=True):
            st.info("Shutting down...")
            os._exit(0)

    # Main content - Tabs
    tab1, tab2, tab3, tab4, tab5 = st.tabs(["📁 Select Data", "📚 Reference Index", "▶️ Run Pipeline", "📊 Results", "🔬 Analysis"])

    # Tab 1: Select Data (Local Directory)
    with tab1:
        st.header("Select FASTQ Data Directory")

        st.markdown("""
        Enter the path to a folder containing your FASTQ files. The app will automatically detect:
        - **Paired-end** samples (files with `_R1`/`_R2` or `_1`/`_2` patterns)
        - **Single-end** samples

        Supported formats: `.fastq`, `.fq`, `.fastq.gz`, `.fq.gz`
        """)

        # Directory input with Browse button
        col_input, col_browse = st.columns([4, 1])

        with col_input:
            data_dir = st.text_input(
                "FASTQ directory path:",
                value=st.session_state.get("data_dir_path", ""),
                help="Enter the full path to your FASTQ files folder",
                placeholder="Click 'Browse' or paste a path..."
            )

        with col_browse:
            st.markdown("<br>", unsafe_allow_html=True)  # Align button with input
            if st.button("📂 Browse"):
                selected = open_folder_dialog(st.session_state.get("data_dir_path"))
                if selected:
                    st.session_state.data_dir_path = selected
                    st.rerun()

        col1, col2 = st.columns([1, 3])
        with col1:
            scan_button = st.button("🔍 Scan Directory", type="primary", disabled=not data_dir)

        if scan_button and data_dir:
            data_path = Path(data_dir)

            if not data_path.exists():
                st.error(f"Directory not found: {data_dir}")
            elif not data_path.is_dir():
                st.error(f"Not a directory: {data_dir}")
            else:
                # Warn if files are on cloud storage (Google Drive, iCloud, Dropbox, OneDrive)
                cloud_keywords = ["CloudStorage", "Google Drive", "GoogleDrive", "iCloud", "Dropbox", "OneDrive"]
                if any(kw in data_dir for kw in cloud_keywords):
                    st.warning(
                        "These files appear to be on cloud storage (e.g. Google Drive). "
                        "Large FASTQ files streamed from cloud storage often fail mid-processing. "
                        "**Copy your FASTQ files to a local folder first** for reliable results."
                    )

                st.session_state.data_dir_path = data_dir
                with st.spinner("Scanning for FASTQ files..."):
                    result = scan_fastq_directory(data_path)
                    st.session_state.detected_samples = result

        # Show detected samples
        if st.session_state.detected_samples:
            result = st.session_state.detected_samples

            if "error" in result:
                st.error(result["error"])
            elif result["samples"]:
                read_type = result["read_type"]
                samples = result["samples"]

                st.success(f"Found {len(samples)} {read_type}-end samples")

                # Display samples in a table
                st.subheader("Detected Samples")

                if read_type == "paired":
                    sample_data = []
                    for s in samples:
                        sample_data.append({
                            "Sample Name": s["name"],
                            "Read 1": s["r1"].name,
                            "Read 2": s["r2"].name
                        })
                    st.dataframe(pd.DataFrame(sample_data), width="stretch")
                else:
                    sample_data = []
                    for s in samples:
                        sample_data.append({
                            "Sample Name": s["name"],
                            "File": s["r1"].name
                        })
                    st.dataframe(pd.DataFrame(sample_data), width="stretch")

                # Store for pipeline
                st.session_state.samples_ready = True

    # Tab 2: Reference Index
    with tab2:
        st.header("Reference Transcriptome Index")

        st.markdown("""
        Select a reference transcriptome for quantification. You can either:
        1. **Download pre-built** - Human or Mouse from GENCODE
        2. **Use existing** - Select a previously downloaded index
        3. **Build from FASTA** - Provide your own transcriptome file path
        """)

        index_source = st.radio(
            "Index source:",
            ["Download pre-built", "Use existing", "Build from FASTA"],
            horizontal=True
        )

        if index_source == "Download pre-built":
            organism = st.selectbox(
                "Select organism:",
                ["human_gencode_v36", "human_gencode_v44", "mouse_gencode_vM25", "mouse_gencode_vM33"],
                format_func=lambda x: {
                    "human_gencode_v36": "Human (GENCODE v36 / GRCh38.p13)",
                    "human_gencode_v44": "Human (GENCODE v44 / GRCh38.p14)",
                    "mouse_gencode_vM25": "Mouse (GENCODE vM25 / GRCm38)",
                    "mouse_gencode_vM33": "Mouse (GENCODE vM33 / GRCm39)"
                }.get(x, x)
            )

            index_path = INDEX_DIR / organism

            if (index_path / "complete_ref_lens.bin").exists():
                st.success(f"Index already downloaded: {index_path}")
                st.session_state.index_path = index_path
            else:
                st.info("📦 This will download ~500MB and build the index (~5-10 minutes)")
                if st.button("Download Index", type="primary"):
                    # Create visual elements for progress
                    progress_bar = st.progress(0, text="Starting download...")
                    status_text = st.empty()

                    # Animated messages for each phase
                    download_messages = [
                        "📥 Downloading transcriptome data...",
                        "📥 Still downloading... (this file is ~500MB)",
                        "📥 Download in progress... please wait",
                        "📥 Fetching sequence data from GENCODE...",
                        "📥 Almost there... downloading remaining data",
                    ]
                    index_messages = [
                        "🔧 Building Salmon index...",
                        "🔧 Indexing transcripts...",
                        "🔧 Constructing k-mer index...",
                        "🔧 Processing sequences...",
                        "🔧 Finalizing index structure...",
                    ]

                    # Shared state for thread communication
                    # Progress breakdown: Download 0-50%, Index 50-90%, Mapping 90-100%
                    download_state = {
                        "phase": "download",  # download, index, mapping, done
                        "download_percent": 0,  # 0-100 for download phase
                        "finished": False,
                        "result": None,
                        "error": None,
                        "log_messages": []  # Thread-safe log storage
                    }

                    def update_progress(msg):
                        # Check for download progress updates
                        if msg.startswith("DOWNLOAD_PROGRESS:"):
                            percent = int(msg.split(":")[1])
                            download_state["download_percent"] = percent
                            return  # Don't log these frequent updates

                        # Store log in thread-safe dict (not session_state)
                        timestamp = datetime.now().strftime("%H:%M:%S")
                        download_state["log_messages"].append(f"[{timestamp}] {msg}")
                        if "Downloading" in msg:
                            download_state["phase"] = "download"
                        elif "Building" in msg:
                            download_state["phase"] = "index"
                        elif "Generating" in msg:
                            download_state["phase"] = "mapping"
                        elif "complete" in msg.lower():
                            download_state["phase"] = "done"

                    def run_download():
                        try:
                            result = download_reference_index(
                                organism,
                                INDEX_DIR,
                                progress_callback=update_progress
                            )
                            download_state["result"] = result
                        except Exception as e:
                            download_state["error"] = str(e)
                        finally:
                            download_state["finished"] = True

                    # Start download in background thread
                    download_thread = threading.Thread(target=run_download, daemon=True)
                    download_thread.start()

                    # Track progress based on actual phase completion
                    # Download: 0-50%, Index: 50-90%, Mapping: 90-100%
                    tick = 0
                    index_start_time = None
                    INDEX_ESTIMATED_SECONDS = 180  # ~3 minutes for index build

                    while not download_state["finished"]:
                        phase = download_state["phase"]

                        if phase == "download":
                            # Real progress from download (0-50%)
                            progress = int(download_state["download_percent"] * 0.5)
                            msg_idx = (tick // 10) % len(download_messages)
                            msg = download_messages[msg_idx]
                            label = "Downloading transcriptome"

                        elif phase == "index":
                            # Index phase: 50-90% based on elapsed time
                            if index_start_time is None:
                                index_start_time = time.time()
                            elapsed = time.time() - index_start_time
                            index_progress = min(elapsed / INDEX_ESTIMATED_SECONDS, 0.95)
                            progress = 50 + int(index_progress * 40)  # 50-90%
                            msg_idx = (tick // 10) % len(index_messages)
                            msg = index_messages[msg_idx]
                            label = "Building Salmon index"

                        elif phase == "mapping":
                            # Mapping is quick: 90-99%
                            progress = 95
                            msg = "📋 Creating transcript-to-gene mappings..."
                            label = "Generating gene mappings"

                        else:
                            progress = 99
                            msg = "✅ Finishing up..."
                            label = "Completing"

                        progress = min(progress, 99)
                        progress_bar.progress(progress, text=f"{label}... {progress}%")
                        status_text.markdown(f"*{msg}*")
                        tick += 1
                        time.sleep(0.5)

                    # Sync log messages back to session state
                    for log_msg in download_state["log_messages"]:
                        st.session_state.log_messages.append(log_msg)

                    # Handle result
                    if download_state["error"]:
                        st.error(f"Failed to download index: {download_state['error']}")
                    else:
                        st.session_state.index_path = download_state["result"]
                        progress_bar.progress(100, text="Complete! 100%")
                        status_text.markdown("✅ *All done!*")
                        st.success("✅ Index downloaded and built successfully!")
                        st.rerun()

        elif index_source == "Use existing":
            existing_indices = [d for d in INDEX_DIR.iterdir() if d.is_dir()] if INDEX_DIR.exists() else []

            if existing_indices:
                selected = st.selectbox(
                    "Select index:",
                    existing_indices,
                    format_func=lambda x: x.name
                )
                st.session_state.index_path = selected
                st.success(f"Using index: {selected}")
            else:
                st.warning("No existing indices found. Please download or build one.")

        else:  # Build from FASTA
            st.markdown("Select a transcriptome FASTA file.")

            col_input, col_browse = st.columns([4, 1])

            with col_input:
                fasta_path_str = st.text_input(
                    "Transcriptome FASTA path:",
                    value=st.session_state.get("fasta_path", ""),
                    help="Full path to .fa, .fasta, or .fa.gz file",
                    placeholder="Click 'Browse' or paste a path..."
                )

            with col_browse:
                st.markdown("<br>", unsafe_allow_html=True)
                if st.button("📂 Browse", key="browse_fasta"):
                    selected = open_file_dialog()
                    if selected:
                        st.session_state.fasta_path = selected
                        st.rerun()

            index_name = st.text_input("Index name:", "custom_index")

            if fasta_path_str and index_name:
                fasta_path = Path(fasta_path_str)

                fasta_error = (
                    f"File not found: {fasta_path_str}" if not fasta_path.exists()
                    else validate_fasta_is_transcriptome(fasta_path)
                )
                if fasta_error:
                    st.error(fasta_error)
                elif st.button("Build Index", type="primary"):
                    progress_placeholder = st.empty()

                    def update_progress(msg):
                        progress_placeholder.text(msg)
                        add_log(msg)

                    try:
                        result_path = build_index_from_fasta(
                            fasta_path,
                            INDEX_DIR,
                            index_name,
                            threads=threads,
                            progress_callback=update_progress
                        )
                        st.session_state.index_path = result_path
                        st.success("Index built successfully!")
                        st.rerun()
                    except Exception as e:
                        st.error(f"Failed to build index: {e}")

    # Tab 3: Run Pipeline
    with tab3:
        st.header("Run Quantification Pipeline")

        # Check prerequisites
        has_samples = (
            st.session_state.get("detected_samples")
            and st.session_state.detected_samples.get("samples")
        )
        has_index = st.session_state.get("index_path")

        col1, col2 = st.columns(2)
        with col1:
            if has_samples:
                n_samples = len(st.session_state.detected_samples["samples"])
                read_type = st.session_state.detected_samples["read_type"]
                st.success(f"✓ {n_samples} {read_type}-end samples ready")
            else:
                st.warning("⚠ No samples detected - go to Select Data tab")

        with col2:
            if has_index:
                st.success(f"✓ Index: {Path(st.session_state.index_path).name}")
            else:
                st.warning("⚠ No index selected - go to Reference Index tab")

        if has_samples and has_index:
            st.divider()

            # Output folder naming
            st.markdown("**Output Settings**")
            # Initialize run name in session state only once
            if "run_name" not in st.session_state:
                st.session_state.run_name = datetime.now().strftime("%Y%m%d_%H%M%S")

            run_name_input = st.text_input(
                "Run name",
                key="run_name",
                help="Name for this analysis run. Results will be saved to a folder with this name."
            )
            # Sanitize run name (remove invalid characters)
            run_name = "".join(c for c in run_name_input if c.isalnum() or c in "_-").strip()
            if not run_name:
                run_name = datetime.now().strftime("%Y%m%d_%H%M%S")

            st.caption(f"📁 Results will be saved to: `~/Documents/GenomifySeq/results/{run_name}/`")

            st.divider()

            # Preprocessing options (optional - add significant time)
            st.markdown("**Preprocessing Options**")
            st.caption("These steps are not required but recommended for publication-quality analysis.")

            col1, col2 = st.columns(2)
            with col1:
                run_fastqc = st.checkbox("Run FASTQC", value=False)
                st.caption(
                    "*Generates comprehensive quality metrics including per-base sequence quality, "
                    "GC content distribution, sequence duplication levels, and adapter contamination detection. "
                    "Essential for identifying problematic samples before downstream analysis.*"
                )
            with col2:
                run_trim_galore = st.checkbox("Run Trim Galore", value=False)
                st.caption(
                    "*Performs adapter trimming and quality filtering using Cutadapt. Removes Illumina adapter "
                    "sequences and low-quality bases from read ends, improving mapping rates and reducing "
                    "technical artifacts in quantification.*"
                )

            st.divider()

            # Time estimate based on real benchmarks:
            # - Salmon only: ~1.5 min/sample @ 14 CPUs
            # - With FASTQC + Trim Galore: ~15 min/sample @ 14 CPUs
            n_samples = len(st.session_state.detected_samples["samples"])

            # Base Salmon time scales with threads
            salmon_time = n_samples * 1.5 * (14 / threads)

            # Preprocessing time (scales less with threads - I/O bound)
            preprocessing_time = 0
            if run_fastqc:
                preprocessing_time += n_samples * 3 * (14 / max(threads, 4))  # ~3 min/sample @ 14 CPUs
            if run_trim_galore:
                preprocessing_time += n_samples * 10 * (14 / max(threads, 4))  # ~10 min/sample @ 14 CPUs

            est_minutes = round(salmon_time + preprocessing_time)
            if est_minutes < 1:
                est_minutes = 1
            st.info(f"⏱️ Estimated processing time: **~{est_minutes} minutes** for {n_samples} samples at {threads} threads")

            with st.expander("View quantification command"):
                st.code(
                    f"salmon quant \\\n"
                    f"  --index <index_dir> \\\n"
                    f"  --libType A \\\n"
                    f"  --mates1 <R1.fastq.gz> \\\n"
                    f"  --mates2 <R2.fastq.gz> \\\n"
                    f"  --output <output_dir> \\\n"
                    f"  --threads {threads} \\\n"
                    f"  --validateMappings \\\n"
                    f"  --gcBias \\\n"
                    f"  --seqBias",
                    language="bash"
                )

            # Run button
            if st.button("🚀 Start Quantification", type="primary", disabled=st.session_state.pipeline_running):
                st.session_state.pipeline_running = True
                st.session_state.log_messages = []

                # Create run directory with user-specified name
                run_dir = RESULTS_DIR / run_name
                run_dir.mkdir(parents=True, exist_ok=True)

                # Build sample objects from detected samples
                detected = st.session_state.detected_samples
                samples = []

                for s in detected["samples"]:
                    samples.append(SampleInfo(
                        name=s["name"],
                        read1=s["r1"],
                        read2=s.get("r2")
                    ))

                # Enhanced progress display
                st.markdown("---")
                progress_header = st.empty()
                progress_cols = st.columns([3, 2, 2, 2])
                progress_bar_container = progress_cols[0].empty()
                sample_counter = progress_cols[1].empty()
                elapsed_display = progress_cols[2].empty()
                eta_display = progress_cols[3].empty()
                current_sample_display = st.empty()
                fun_status = st.empty()
                stop_hint = st.empty()
                log_display = st.empty()

                # Fun status messages to rotate through
                fun_messages = [
                    "Aligning reads to transcriptome...",
                    "Counting molecules...",
                    "Crunching numbers...",
                    "Teaching salmon to swim upstream...",
                    "Herding RNA molecules...",
                    "Convincing genes to express themselves...",
                    "Calibrating the quantification engine...",
                    "Mapping reads at light speed...",
                    "Decoding the transcriptome...",
                    "Making science happen...",
                    "Wrangling transcripts...",
                    "Sequencing at warp speed...",
                    "Converting biology to numbers...",
                    "Feeding the salmon algorithm...",
                    "Discovering gene expression patterns...",
                ]

                # Spinning animation frames
                spinner_frames = ["🧬", "🔬", "🧪", "📊", "🧫", "💊", "🔭", "⚗️"]

                total_samples = len(samples)

                # Shared state for thread communication (not session_state - not thread safe)
                pipeline_state = {
                    "current_sample": 0,
                    "completed_samples": 0,
                    "start_time": time.time(),
                    "sample_start_time": time.time(),
                    "sample_durations": [],
                    "current_sample_name": "",
                    "last_log": "",
                    "log_messages": [],  # Thread-safe log storage
                    "finished": False,
                    "results": None,
                    "error": None,
                }

                # Configure pipeline
                config = PipelineConfig(
                    index_path=Path(st.session_state.index_path),
                    output_dir=run_dir / "quant",
                    threads=threads,
                    run_fastqc=run_fastqc,
                    run_trim_galore=run_trim_galore
                )

                def pipeline_progress(msg):
                    """Callback from pipeline - updates shared state (thread-safe)."""
                    now = time.time()
                    timestamp = datetime.now().strftime("%H:%M:%S")
                    pipeline_state["last_log"] = msg

                    # When starting a new sample (first message for each sample)
                    if "Processing sample" in msg:
                        pipeline_state["current_sample"] += 1
                        pipeline_state["sample_start_time"] = now
                        # Extract sample name from "Processing sample X/Y: sample_name"
                        if ": " in msg:
                            sample_name = msg.split(": ", 1)[1]
                        else:
                            sample_name = f"Sample {pipeline_state['current_sample']}"
                        pipeline_state["current_sample_name"] = sample_name
                        pipeline_state["log_messages"].append(f"[{timestamp}] Analyzing {sample_name}...")

                    # Log preprocessing and quantification steps
                    elif "FASTQC" in msg or "Trimming" in msg or "Quantifying" in msg:
                        pipeline_state["log_messages"].append(f"[{timestamp}] {msg}")

                    # When a sample completes
                    elif "Completed" in msg:
                        pipeline_state["completed_samples"] += 1
                        sample_duration = now - pipeline_state["sample_start_time"]
                        pipeline_state["sample_durations"].append(sample_duration)
                        sample_name = pipeline_state["current_sample_name"]
                        pipeline_state["log_messages"].append(f"[{timestamp}] Completed {sample_name} ({format_time(sample_duration)})")

                def run_pipeline_thread():
                    """Run the pipeline in a background thread."""
                    try:
                        results = run_pipeline(
                            samples, config,
                            progress_callback=pipeline_progress
                        )
                        pipeline_state["results"] = results
                    except Exception as e:
                        pipeline_state["error"] = str(e)
                    finally:
                        pipeline_state["finished"] = True

                # Calculate expected time per sample based on options
                # Base: ~90s for Salmon only, add ~180s for FASTQC, ~600s for Trim Galore
                base_seconds_per_sample = 90  # Salmon quantification
                if run_fastqc:
                    base_seconds_per_sample += 180  # ~3 min for FASTQC
                if run_trim_galore:
                    base_seconds_per_sample += 600  # ~10 min for Trim Galore

                # Initialize displays
                progress_header.markdown("### Pipeline Progress")
                progress_bar_container.progress(0, text="0%")
                sample_counter.metric("Samples", f"0 / {total_samples}")
                elapsed_display.metric("Elapsed", "0s")
                # Initial ETA based on selected options
                initial_eta = base_seconds_per_sample * total_samples
                eta_display.metric("ETA", f"~{format_time(initial_eta)}")
                current_sample_display.info("Starting pipeline...")
                fun_status.markdown("🧬 *Warming up the quantification engine...*")
                stop_hint.caption("💡 *To stop the analysis, refresh the page or close this tab.*")

                # Add initial log message to thread-safe storage
                timestamp = datetime.now().strftime("%H:%M:%S")
                pipeline_state["log_messages"].append(f"[{timestamp}] Starting pipeline ({len(samples)} samples)...")

                # Prevent system sleep during processing (macOS)
                caffeinate_proc = None
                try:
                    caffeinate_proc = subprocess.Popen(
                        ["caffeinate", "-i", "-s"],
                        stdout=subprocess.DEVNULL,
                        stderr=subprocess.DEVNULL
                    )
                except Exception:
                    pass  # caffeinate not available (non-macOS)

                # Start pipeline in background thread
                pipeline_thread = threading.Thread(target=run_pipeline_thread, daemon=True)
                pipeline_thread.start()

                # Main UI update loop - runs every 0.5 seconds
                msg_idx = 0
                while not pipeline_state["finished"]:
                    now = time.time()
                    elapsed = now - pipeline_state["start_time"]

                    # Update elapsed time (real-time!)
                    elapsed_display.metric("Elapsed", format_time(elapsed))

                    # Update progress bar and sample counter with smooth interpolation
                    current = pipeline_state["current_sample"]
                    completed = pipeline_state["completed_samples"]

                    if current > 0:
                        # Base progress from completed samples
                        base_progress = completed / total_samples

                        # Interpolate progress within current sample
                        sample_elapsed = now - pipeline_state["sample_start_time"]
                        if pipeline_state["sample_durations"]:
                            # Use average of completed samples as estimate
                            avg_duration = sum(pipeline_state["sample_durations"]) / len(pipeline_state["sample_durations"])
                        else:
                            # First sample: use estimate based on selected options
                            avg_duration = float(base_seconds_per_sample)

                        # Calculate fraction of current sample completed (cap at 0.95 to avoid overshooting)
                        sample_fraction = min(sample_elapsed / avg_duration, 0.95)

                        # Total progress = completed samples + fraction of current sample
                        progress_pct = base_progress + (sample_fraction / total_samples)
                        progress_pct = min(progress_pct, 0.99)  # Never show 100% until truly done

                        progress_bar_container.progress(progress_pct, text=f"{int(progress_pct * 100)}%")
                        sample_counter.metric("Samples", f"{current} / {total_samples}")

                    # Update current sample display
                    sample_name = pipeline_state["current_sample_name"]
                    if sample_name:
                        current_sample_display.info(f"Processing: **{sample_name}**")

                    # Calculate and display ETA (show estimate immediately)
                    if current > 0:
                        sample_elapsed = now - pipeline_state["sample_start_time"]

                        if pipeline_state["sample_durations"]:
                            # Use actual completed sample times
                            avg_time = sum(pipeline_state["sample_durations"]) / len(pipeline_state["sample_durations"])
                        else:
                            # First sample: use estimate based on selected options
                            avg_time = float(base_seconds_per_sample)

                        # Calculate remaining time:
                        # - Remaining time in current sample
                        # - Plus all remaining samples after current
                        remaining_samples_after_current = total_samples - current
                        current_sample_remaining = max(0, avg_time - sample_elapsed)
                        eta_seconds = current_sample_remaining + (remaining_samples_after_current * avg_time)
                        eta_display.metric("ETA", f"~{format_time(eta_seconds)}")

                    # Rotate fun message and spinner every 5 seconds (10 iterations)
                    fun_msg_idx = msg_idx // 10  # Change message every 5 seconds
                    spinner = spinner_frames[msg_idx % len(spinner_frames)]
                    fun_msg = fun_messages[fun_msg_idx % len(fun_messages)]
                    fun_status.markdown(f"{spinner} *{fun_msg}*")
                    msg_idx += 1

                    # Update log display from thread-safe storage
                    log_display.text("\n".join(pipeline_state["log_messages"][-5:]))

                    # Sleep briefly before next update
                    time.sleep(0.5)

                # Pipeline finished - process results
                results = pipeline_state["results"]
                stop_hint.empty()  # Clear the stop hint

                # Sync log messages back to session state (now safe - back in main thread)
                st.session_state.log_messages = pipeline_state["log_messages"].copy()

                if pipeline_state["error"]:
                    st.error(f"Pipeline failed: {pipeline_state['error']}")
                    add_log(f"ERROR: {pipeline_state['error']}")
                else:
                    # Normal completion
                    current_sample_display.info("Merging count matrices...")

                    # Merge counts
                    if results and (results.get("success") or any(s["success"] for s in results.get("samples", []))):
                        add_log("Merging count matrices...")
                        sample_dirs = [Path(s["output_dir"]) for s in results["samples"] if s["success"]]

                        if sample_dirs:
                            counts_file = run_dir / "counts_raw.csv"
                            tpm_file = run_dir / "counts_tpm.csv"
                            merge_counts(sample_dirs, counts_file, "NumReads")
                            merge_counts(sample_dirs, tpm_file, "TPM")
                            results["counts_file"] = str(counts_file)
                            results["tpm_file"] = str(tpm_file)

                    # Final progress update
                    total_elapsed = time.time() - pipeline_state["start_time"]
                    progress_bar_container.progress(1.0, text="100%")
                    elapsed_display.metric("Elapsed", format_time(total_elapsed))
                    eta_display.metric("ETA", "Complete!")
                    sample_counter.metric("Samples", f"{total_samples} / {total_samples}")
                    fun_status.markdown("🎉 *All done! Your RNA has been thoroughly sequenced.*")

                    st.session_state.results = results
                    st.session_state.results["run_dir"] = str(run_dir)
                    st.session_state.results["index_path"] = str(st.session_state.index_path)

                    # Save index path to run directory for future sessions
                    (run_dir / "index_path.txt").write_text(str(st.session_state.index_path))

                    successful = sum(1 for s in results["samples"] if s["success"])

                    if successful == len(samples):
                        current_sample_display.success("Pipeline completed successfully!")
                        add_log("Pipeline completed successfully!")
                        st.success("✓ Quantification complete! Go to Results tab to view.")
                    elif successful > 0:
                        current_sample_display.warning(f"Completed with {len(samples) - successful} failures")
                        st.warning(f"Completed with {len(samples) - successful} failures. Check Results tab.")
                    else:
                        current_sample_display.error("Pipeline failed")
                        st.error(f"Pipeline failed: {results.get('errors', ['Unknown error'])}")

                st.session_state.pipeline_running = False

                # Stop preventing sleep
                if caffeinate_proc:
                    caffeinate_proc.terminate()

        # Show log
        if st.session_state.log_messages:
            with st.expander("View Log", expanded=False):
                st.code("\n".join(st.session_state.log_messages))

    # Tab 4: Results
    with tab4:
        st.header("Results")

        # Find previous runs
        previous_runs = []
        if RESULTS_DIR.exists():
            for run_folder in sorted(RESULTS_DIR.iterdir(), reverse=True):
                if run_folder.is_dir() and (run_folder / "quant").exists():
                    # Try to get sample count
                    quant_dir = run_folder / "quant"
                    sample_count = len([d for d in quant_dir.iterdir() if d.is_dir()]) if quant_dir.exists() else 0
                    previous_runs.append({
                        "path": run_folder,
                        "name": run_folder.name,
                        "samples": sample_count
                    })

        # Run selector
        if previous_runs or st.session_state.results:
            st.subheader("Select Run")

            run_options = []
            if st.session_state.results:
                run_options.append("Current Session")
            run_options.extend([f"{r['name']} ({r['samples']} samples)" for r in previous_runs])

            selected_run = st.selectbox(
                "Load results from:",
                options=run_options,
                index=0
            )

            # Determine which results to show
            if selected_run == "Current Session":
                results = st.session_state.results
                run_dir = Path(results["run_dir"]) if results else None
            else:
                # Find the selected previous run
                run_idx = run_options.index(selected_run) - (1 if st.session_state.results else 0)
                run_dir = previous_runs[run_idx]["path"]

                # Load results from the previous run
                quant_dir = run_dir / "quant"
                sample_dirs = [d for d in quant_dir.iterdir() if d.is_dir()] if quant_dir.exists() else []

                # Recover index path from saved file or auto-detect from transcript IDs
                saved_index = None
                index_path_file = run_dir / "index_path.txt"
                if index_path_file.exists():
                    saved_index = index_path_file.read_text().strip()
                else:
                    # Auto-detect species from quant.sf transcript IDs
                    for sd in sample_dirs:
                        qf = sd / "quant.sf"
                        if qf.exists() and qf.stat().st_size > 0:
                            with open(qf) as f:
                                f.readline()  # skip header
                                first_tid = f.readline().split('\t')[0]
                            is_mouse = first_tid.startswith("ENSMUST")
                            # Find matching index directory
                            if INDEX_DIR.exists():
                                for idx_dir in sorted(INDEX_DIR.iterdir(), reverse=True):
                                    if idx_dir.is_dir() and (idx_dir / "tx2gene.tsv").exists():
                                        if is_mouse and "mouse" in idx_dir.name:
                                            saved_index = str(idx_dir)
                                            break
                                        elif not is_mouse and "human" in idx_dir.name:
                                            saved_index = str(idx_dir)
                                            break
                            break

                results = {
                    "run_dir": str(run_dir),
                    "index_path": saved_index,
                    "samples": [],
                    "success": True
                }

                for sample_dir in sample_dirs:
                    quant_file = sample_dir / "quant.sf"
                    results["samples"].append({
                        "name": sample_dir.name,
                        "output_dir": str(sample_dir),
                        "quant_file": str(quant_file),
                        "success": quant_file.exists()
                    })

                # Check for count files
                counts_file = run_dir / "counts_raw.csv"
                tpm_file = run_dir / "counts_tpm.csv"
                if counts_file.exists():
                    results["counts_file"] = str(counts_file)
                if tpm_file.exists():
                    results["tpm_file"] = str(tpm_file)

                # Update session state so Analysis tab can use this run
                st.session_state.results = results

            st.divider()
        else:
            results = None
            run_dir = None

        if results is None:
            st.info("No results yet. Run the pipeline first, or previous runs will appear here.")
        else:
            run_dir = Path(results["run_dir"])

            # Summary
            st.subheader("Summary")
            col1, col2, col3 = st.columns(3)

            successful = sum(1 for s in results["samples"] if s["success"])
            failed = len(results["samples"]) - successful

            col1.metric("Samples Processed", len(results["samples"]))
            col2.metric("Successful", successful)
            col3.metric("Failed", failed)

            st.caption(f"Results saved to: `{run_dir}`")

            # Mapping statistics
            st.subheader("Mapping Statistics")

            stats_data = []
            for sample in results["samples"]:
                if sample["success"]:
                    sample_dir = Path(sample["output_dir"])
                    stats = get_mapping_stats(sample_dir)
                    if stats:
                        stats_data.append({
                            "Sample": sample["name"],
                            "Total Reads": f"{stats.get('num_processed', 0):,}",
                            "Mapped Reads": f"{stats.get('num_mapped', 0):,}",
                            "Mapping Rate (%)": round(stats.get("mapping_rate", 0), 2)
                        })

            if stats_data:
                stats_df = pd.DataFrame(stats_data)
                st.dataframe(stats_df, width="stretch")

                # Mapping rate plot
                plot_df = pd.DataFrame([{
                    "Sample": d["Sample"],
                    "Mapping Rate (%)": d["Mapping Rate (%)"]
                } for d in stats_data])

                fig_mapping = px.bar(
                    plot_df,
                    x="Sample",
                    y="Mapping Rate (%)",
                    title="Mapping Rate by Sample",
                    color="Mapping Rate (%)",
                    color_continuous_scale="RdYlGn",
                    range_color=[0, 100]
                )
                fig_mapping.update_layout(showlegend=False)
                st.plotly_chart(fig_mapping, width="stretch")

                # Export mapping rate chart
                exp_col1, exp_col2, exp_col3 = st.columns([1, 1, 2])
                with exp_col1:
                    try:
                        png_bytes = fig_mapping.to_image(format="png", width=1000, height=600, scale=2)
                        st.download_button(
                            "📥 Chart (PNG)",
                            png_bytes,
                            file_name="mapping_rate.png",
                            mime="image/png",
                            key="mapping_png"
                        )
                    except Exception:
                        pass
                with exp_col2:
                    try:
                        tiff_bytes = plotly_to_tiff(fig_mapping, width=1000, height=600)
                        st.download_button(
                            "📥 Chart (TIFF)",
                            tiff_bytes,
                            file_name="mapping_rate.tiff",
                            mime="image/tiff",
                            key="mapping_tiff"
                        )
                    except Exception:
                        pass

            # Count matrix preview
            st.subheader("Count Matrix Preview")

            if "counts_file" in results and Path(results["counts_file"]).exists():
                counts_df = pd.read_csv(results["counts_file"], index_col=0)

                # Show top genes
                st.markdown(f"**Top 20 transcripts by total counts** ({len(counts_df):,} total)")
                top_genes = counts_df.sum(axis=1).nlargest(20)
                st.dataframe(counts_df.loc[top_genes.index].round(2), width="stretch")

                # Expression heatmap
                if len(counts_df.columns) > 1 and len(top_genes) > 1:
                    st.subheader("Expression Heatmap (Top 20 transcripts)")

                    heatmap_data = np.log2(counts_df.loc[top_genes.index] + 1)

                    fig_heatmap = px.imshow(
                        heatmap_data,
                        aspect="auto",
                        color_continuous_scale="Viridis",
                        labels={"color": "log2(count+1)"}
                    )
                    st.plotly_chart(fig_heatmap, width="stretch")

                    # Export heatmap
                    hm_col1, hm_col2, hm_col3 = st.columns([1, 1, 2])
                    with hm_col1:
                        try:
                            png_bytes = fig_heatmap.to_image(format="png", width=1000, height=800, scale=2)
                            st.download_button(
                                "📥 Heatmap (PNG)",
                                png_bytes,
                                file_name="expression_heatmap.png",
                                mime="image/png",
                                key="heatmap_png"
                            )
                        except Exception:
                            pass
                    with hm_col2:
                        try:
                            tiff_bytes = plotly_to_tiff(fig_heatmap, width=1000, height=800)
                            st.download_button(
                                "📥 Heatmap (TIFF)",
                                tiff_bytes,
                                file_name="expression_heatmap.tiff",
                                mime="image/tiff",
                                key="heatmap_tiff"
                            )
                        except Exception:
                            pass

            # Download section
            st.subheader("Download Results")

            col1, col2 = st.columns(2)

            with col1:
                if "counts_file" in results:
                    counts_path = Path(results["counts_file"])
                    if counts_path.exists():
                        with open(counts_path, "rb") as f:
                            st.download_button(
                                "📥 Download Raw Counts (CSV)",
                                f,
                                file_name="counts_raw.csv",
                                mime="text/csv"
                            )

            with col2:
                if "tpm_file" in results:
                    tpm_path = Path(results["tpm_file"])
                    if tpm_path.exists():
                        with open(tpm_path, "rb") as f:
                            st.download_button(
                                "📥 Download TPM Counts (CSV)",
                                f,
                                file_name="counts_tpm.csv",
                                mime="text/csv"
                            )

            # Show failed samples if any
            if failed > 0:
                st.subheader("Failed Samples")
                for sample in results["samples"]:
                    if not sample["success"]:
                        st.error(f"**{sample['name']}**: {sample.get('error', 'Unknown error')}")

    # Tab 5: Analysis
    with tab5:
        st.header("Differential Expression Analysis")

        st.markdown("""
        Perform PCA and differential expression analysis on your quantified samples.

        **Steps:**
        1. Fill in the sample metadata table below
        2. Select analysis options
        3. Run analysis
        """)

        with st.expander("View analysis code (analysis.R)"):
            r_script = Path(__file__).parent / "analysis.R"
            if r_script.exists():
                st.code(r_script.read_text(), language="r")
            else:
                st.warning("analysis.R not found.")

        # Check if we have results to analyze
        if st.session_state.results is None:
            st.warning("⚠️ No quantification results found. Please run the pipeline first (Run Pipeline tab).")
        else:
            results = st.session_state.results
            run_dir = Path(results.get("run_dir", ""))
            quant_dir = run_dir / "quant"

            # Get successful samples
            successful_samples = [s for s in results.get("samples", []) if s.get("success")]

            if not successful_samples:
                st.error("No successfully quantified samples found.")
            else:
                st.success(f"Found {len(successful_samples)} samples ready for analysis")

                # Create/edit metadata table
                st.subheader("1. Sample Metadata")

                st.markdown("""
                Fill in the metadata for each sample. The **condition** column is used for grouping samples in the PCA plot.
                You can add additional columns by typing new column names.
                """)

                # Initialize metadata dataframe if not exists or samples changed
                sample_names = [s["name"] for s in successful_samples]

                # Only initialize if metadata_df doesn't exist or has different samples
                needs_init = False
                if st.session_state.metadata_df is None:
                    needs_init = True
                elif list(st.session_state.metadata_df["sample"]) != sample_names:
                    # Only reset if the actual sample list changed (not just set comparison)
                    needs_init = True

                if needs_init:
                    # Create new metadata with sample names
                    st.session_state.metadata_df = pd.DataFrame({
                        "sample": sample_names,
                        "condition": [""] * len(sample_names),
                        "replicate": [""] * len(sample_names)
                    })
                    # Clear the editor key to force refresh
                    if "metadata_editor" in st.session_state:
                        del st.session_state["metadata_editor"]

                # Option 1: Upload metadata CSV
                st.markdown("**Option 1: Upload metadata CSV**")
                uploaded_csv = st.file_uploader(
                    "Upload CSV with 'sample' and 'condition' columns",
                    type="csv",
                    help="CSV must have 'sample' and 'condition' columns. Additional columns will be preserved.",
                    key="metadata_csv_upload"
                )

                if uploaded_csv:
                    try:
                        uploaded_df = pd.read_csv(uploaded_csv)

                        # Validate required columns exist
                        if "sample" not in uploaded_df.columns or "condition" not in uploaded_df.columns:
                            st.error("CSV must have 'sample' and 'condition' columns")
                        else:
                            # Match to detected samples - merge with existing sample list
                            # Keep all samples from quantification, fill in metadata from CSV where available
                            merged_df = pd.DataFrame({"sample": sample_names})

                            # Get all columns from uploaded CSV except 'sample'
                            for col in uploaded_df.columns:
                                if col != "sample":
                                    merged_df[col] = ""

                            # Fill in values from uploaded CSV where sample names match
                            for idx, row in merged_df.iterrows():
                                sample = row["sample"]
                                match = uploaded_df[uploaded_df["sample"] == sample]
                                if len(match) > 0:
                                    for col in uploaded_df.columns:
                                        if col != "sample":
                                            merged_df.at[idx, col] = match.iloc[0][col]

                            # Ensure condition and replicate columns exist
                            if "condition" not in merged_df.columns:
                                merged_df["condition"] = ""
                            if "replicate" not in merged_df.columns:
                                merged_df["replicate"] = ""

                            st.session_state.metadata_df = merged_df

                            # Count matches
                            matched = sum(1 for s in sample_names if s in uploaded_df["sample"].values)
                            st.success(f"Metadata loaded! Matched {matched} of {len(sample_names)} samples.")

                            if matched < len(sample_names):
                                unmatched = [s for s in sample_names if s not in uploaded_df["sample"].values]
                                st.warning(f"Unmatched samples (please fill manually): {', '.join(unmatched[:5])}" +
                                          (f"... and {len(unmatched)-5} more" if len(unmatched) > 5 else ""))

                    except Exception as e:
                        st.error(f"Failed to read CSV: {e}")

                st.markdown("**Option 2: Fill in manually**")
                st.caption("Fill in the table below, then click 'Save Metadata' to confirm your changes.")

                # Build column config dynamically based on current columns
                column_config = {
                    "sample": st.column_config.TextColumn(
                        "Sample",
                        help="Sample name (from quantification)",
                        disabled=True
                    ),
                    "condition": st.column_config.TextColumn(
                        "Condition",
                        help="Experimental condition/group (e.g., WT, KO, Treatment)",
                        required=True
                    ),
                    "replicate": st.column_config.TextColumn(
                        "Replicate",
                        help="Replicate number or ID"
                    )
                }

                # Add config for any additional columns
                for col in st.session_state.metadata_df.columns:
                    if col not in column_config:
                        column_config[col] = st.column_config.TextColumn(col)

                # Use a form to prevent rerun while editing
                with st.form("metadata_form"):
                    edited_df = st.data_editor(
                        st.session_state.metadata_df,
                        num_rows="fixed",
                        use_container_width=True,
                        column_config=column_config,
                        hide_index=True
                    )

                    col_save, col_add = st.columns([1, 2])
                    with col_save:
                        save_clicked = st.form_submit_button("💾 Save Metadata", type="primary")

                    if save_clicked:
                        st.session_state.metadata_df = edited_df
                        st.success("Metadata saved!")

                # Add new column option (outside form)
                with st.expander("Add custom column"):
                    new_col_name = st.text_input("New column name:", key="new_col_name")
                    if st.button("Add Column") and new_col_name:
                        if new_col_name not in st.session_state.metadata_df.columns:
                            st.session_state.metadata_df[new_col_name] = ""
                            st.rerun()
                        else:
                            st.warning(f"Column '{new_col_name}' already exists")

                st.divider()

                # Analysis options
                st.subheader("2. Analysis Options")

                # Get available columns for coloring (exclude 'sample')
                color_options = [c for c in st.session_state.metadata_df.columns if c != "sample"]
                color_by = st.selectbox(
                    "Color PCA by:",
                    options=color_options,
                    index=0 if color_options else 0,
                    help="Select which metadata column to use for coloring the PCA plot"
                )

                # Auto-detect species from index path (used by analysis.R fallback)
                species = detect_species()

                st.divider()

                # Run analysis button
                st.subheader("3. Run Analysis")

                # Check if metadata is filled
                conditions = st.session_state.metadata_df["condition"].tolist()
                metadata_complete = all(c and c.strip() for c in conditions)

                if not metadata_complete:
                    st.warning("⚠️ Please fill in the 'condition' column for all samples before running analysis.")

                if st.button("🔬 Run PCA Analysis", type="primary", disabled=not metadata_complete):
                    with st.spinner("Running analysis..."):
                        try:
                            # Create analysis output directory
                            analysis_dir = run_dir / "analysis"
                            analysis_dir.mkdir(exist_ok=True)

                            # Save metadata to CSV
                            metadata_file = analysis_dir / "metadata.csv"
                            st.session_state.metadata_df.to_csv(metadata_file, index=False)

                            # Find tx2gene file from the index used for this run
                            tx2gene_arg = "NULL"
                            run_index_path = results.get("index_path") or st.session_state.get("index_path")
                            if run_index_path:
                                tx2gene_in_index = Path(run_index_path) / "tx2gene.tsv"
                                if tx2gene_in_index.exists():
                                    tx2gene_arg = str(tx2gene_in_index)

                            # Fallback: check data/tx2gene cache
                            if tx2gene_arg == "NULL":
                                tx2gene_cache = TX2GENE_DIR / f"{species}_tx2gene.tsv"
                                if tx2gene_cache.exists():
                                    tx2gene_arg = str(tx2gene_cache)

                            # Run R analysis script
                            add_log("Starting R analysis...")
                            add_log(f"Quant directory: {quant_dir}")
                            add_log(f"Metadata file: {metadata_file}")

                            cmd = [
                                "Rscript",
                                str(APP_DIR / "analysis.R"),
                                "pca",
                                str(quant_dir),
                                str(metadata_file),
                                str(analysis_dir),
                                tx2gene_arg,
                                species,
                                color_by
                            ]

                            result = subprocess.run(
                                cmd,
                                capture_output=True,
                                text=True,
                                cwd=str(APP_DIR)
                            )

                            if result.returncode != 0:
                                st.error("Analysis failed!")
                                st.code(result.stderr)
                                add_log(f"ERROR: {result.stderr}")
                            else:
                                add_log("Analysis completed successfully!")

                                # Store results
                                st.session_state.analysis_results = {
                                    "analysis_dir": str(analysis_dir),
                                    "pca_plot": str(analysis_dir / "pca_plot.png"),
                                    "pca_data": str(analysis_dir / "pca_data.csv"),
                                    "variance_explained": str(analysis_dir / "variance_explained.csv"),
                                    "normalized_counts": str(analysis_dir / "normalized_counts.csv")
                                }

                                st.success("✓ Analysis complete!")
                                st.rerun()

                        except Exception as e:
                            st.error(f"Analysis failed: {e}")
                            add_log(f"ERROR: {e}")

                # Display results
                if st.session_state.analysis_results:
                    st.divider()
                    st.subheader("4. Results")

                    analysis = st.session_state.analysis_results
                    analysis_dir = Path(analysis["analysis_dir"])

                    # Show PCA plot
                    pca_plot_path = Path(analysis["pca_plot"])
                    if pca_plot_path.exists():
                        st.image(str(pca_plot_path), caption="PCA Plot", width="stretch")

                    # Show variance explained
                    var_file = Path(analysis["variance_explained"])
                    if var_file.exists():
                        var_df = pd.read_csv(var_file)
                        st.markdown("**Variance Explained by Principal Components:**")

                        # Bar chart of variance
                        fig = px.bar(
                            var_df.head(10),
                            x="PC",
                            y="variance_explained",
                            title="Variance Explained by PC",
                            labels={"variance_explained": "Variance Explained (%)"}
                        )
                        st.plotly_chart(fig, width="stretch")

                    # Show PCA data table
                    pca_data_file = Path(analysis["pca_data"])
                    if pca_data_file.exists():
                        with st.expander("View PCA Data"):
                            pca_df = pd.read_csv(pca_data_file)
                            st.dataframe(pca_df, width="stretch")

                    # Download buttons
                    st.subheader("Download Results")

                    col1, col2, col3, col4 = st.columns(4)

                    with col1:
                        if pca_plot_path.exists():
                            with open(pca_plot_path, "rb") as f:
                                st.download_button(
                                    "📥 PCA (PNG)",
                                    f,
                                    file_name="pca_plot.png",
                                    mime="image/png"
                                )

                    with col2:
                        pdf_path = analysis_dir / "pca_plot.pdf"
                        if pdf_path.exists():
                            with open(pdf_path, "rb") as f:
                                st.download_button(
                                    "📥 PCA (PDF)",
                                    f,
                                    file_name="pca_plot.pdf",
                                    mime="application/pdf"
                                )

                    with col3:
                        if pca_plot_path.exists():
                            try:
                                tiff_data = png_to_tiff(pca_plot_path)
                                st.download_button(
                                    "📥 PCA (TIFF)",
                                    tiff_data,
                                    file_name="pca_plot.tiff",
                                    mime="image/tiff"
                                )
                            except Exception:
                                pass  # Skip if conversion fails

                    with col4:
                        norm_counts_file = Path(analysis["normalized_counts"])
                        if norm_counts_file.exists():
                            with open(norm_counts_file, "rb") as f:
                                st.download_button(
                                    "📥 Counts (CSV)",
                                    f,
                                    file_name="normalized_counts.csv",
                                    mime="text/csv"
                                )

                # --- Differential Expression Analysis ---
                st.divider()
                st.subheader("5. Differential Expression Analysis")

                st.markdown("""
                Compare two conditions to find differentially expressed genes.
                Select a reference (control) and test condition from your metadata.
                """)

                # Get unique conditions from metadata
                if st.session_state.metadata_df is not None and "condition" in st.session_state.metadata_df.columns:
                    conditions_list = sorted(st.session_state.metadata_df["condition"].dropna().unique().tolist())
                    conditions_list = [c for c in conditions_list if c.strip()]
                else:
                    conditions_list = []

                if len(conditions_list) < 2:
                    st.warning("⚠️ Need at least 2 different conditions in the metadata to run DE analysis.")
                else:
                    col1, col2 = st.columns(2)

                    with col1:
                        ref_condition = st.selectbox(
                            "Reference condition (control):",
                            options=conditions_list,
                            index=0,
                            help="The baseline/control condition (denominator in fold change)"
                        )

                    with col2:
                        # Filter out the reference from test options
                        test_options = [c for c in conditions_list if c != ref_condition]
                        test_condition = st.selectbox(
                            "Test condition:",
                            options=test_options,
                            index=0,
                            help="The condition to compare against reference (numerator in fold change)"
                        )

                    # Thresholds
                    col1, col2 = st.columns(2)
                    with col1:
                        lfc_threshold = st.number_input(
                            "Log2 fold change threshold:",
                            min_value=0.0,
                            value=1.0,
                            step=0.5,
                            help="Genes with |LFC| above this are considered significant"
                        )
                    with col2:
                        padj_threshold = st.number_input(
                            "Adjusted p-value threshold:",
                            min_value=0.001,
                            max_value=0.1,
                            value=0.05,
                            step=0.01,
                            format="%.3f",
                            help="Genes with padj below this are considered significant"
                        )

                    # Show sample counts per condition
                    ref_n = sum(st.session_state.metadata_df["condition"] == ref_condition)
                    test_n = sum(st.session_state.metadata_df["condition"] == test_condition)
                    st.info(f"Comparing **{test_condition}** (n={test_n}) vs **{ref_condition}** (n={ref_n})")

                    if st.button("🔬 Run Differential Expression", type="primary"):
                        with st.spinner("Running DESeq2 analysis..."):
                            try:
                                analysis_dir = Path(results["run_dir"]) / "analysis"
                                de_dir = analysis_dir / "de"
                                de_dir.mkdir(parents=True, exist_ok=True)

                                # Save metadata
                                metadata_file = analysis_dir / "metadata.csv"
                                st.session_state.metadata_df.to_csv(metadata_file, index=False)

                                # Find tx2gene from the index used for this run
                                tx2gene_arg = "NULL"
                                run_index_path = results.get("index_path") or st.session_state.get("index_path")
                                if run_index_path:
                                    tx2gene_in_index = Path(run_index_path) / "tx2gene.tsv"
                                    if tx2gene_in_index.exists():
                                        tx2gene_arg = str(tx2gene_in_index)
                                if tx2gene_arg == "NULL":
                                    tx2gene_cache = TX2GENE_DIR / f"{species}_tx2gene.tsv"
                                    if tx2gene_cache.exists():
                                        tx2gene_arg = str(tx2gene_cache)

                                add_log(f"Running DE: {test_condition} vs {ref_condition}")

                                cmd = [
                                    "Rscript",
                                    str(APP_DIR / "analysis.R"),
                                    "de",
                                    str(quant_dir),
                                    str(metadata_file),
                                    str(de_dir),
                                    tx2gene_arg,
                                    "condition",
                                    ref_condition,
                                    test_condition,
                                    str(lfc_threshold),
                                    str(padj_threshold)
                                ]

                                result = subprocess.run(
                                    cmd,
                                    capture_output=True,
                                    text=True,
                                    cwd=str(APP_DIR)
                                )

                                if result.returncode != 0:
                                    st.error("DE analysis failed!")
                                    st.code(result.stderr)
                                    add_log(f"ERROR: {result.stderr}")
                                else:
                                    add_log("DE analysis completed!")
                                    comparison = f"{test_condition}_vs_{ref_condition}"

                                    st.session_state.de_results = {
                                        "de_dir": str(de_dir),
                                        "comparison": comparison,
                                        "volcano": str(de_dir / f"{comparison}_volcano.png"),
                                        "volcano_pdf": str(de_dir / f"{comparison}_volcano.pdf"),
                                        "all_genes": str(de_dir / f"{comparison}_all_genes.csv"),
                                        "sig_genes": str(de_dir / f"{comparison}_significant_genes.csv"),
                                        "summary": str(de_dir / f"{comparison}_summary.csv"),
                                    }
                                    st.success("✓ DE analysis complete!")
                                    st.rerun()

                            except Exception as e:
                                st.error(f"DE analysis failed: {e}")
                                add_log(f"ERROR: {e}")

                    # Display DE results
                    if st.session_state.de_results:
                        de = st.session_state.de_results
                        de_dir = Path(de["de_dir"])

                        st.divider()
                        st.subheader("6. DE Results")

                        # Summary stats
                        summary_file = Path(de["summary"])
                        if summary_file.exists():
                            summary_df = pd.read_csv(summary_file)
                            col1, col2, col3 = st.columns(3)
                            col1.metric("Total Genes", f"{summary_df['total_genes'].iloc[0]:,}")
                            col2.metric("Upregulated", summary_df["upregulated"].iloc[0])
                            col3.metric("Downregulated", summary_df["downregulated"].iloc[0])

                        # Interactive Volcano plot
                        all_genes_file = Path(de["all_genes"])
                        volcano_path = Path(de["volcano"])

                        if all_genes_file.exists():
                            all_genes_df = pd.read_csv(all_genes_file)

                            st.markdown("**Interactive Volcano Plot**")
                            st.caption("Hover over points to see gene info. Use the controls below to search and label genes.")

                            # Search and label controls
                            vol_col1, vol_col2 = st.columns(2)

                            with vol_col1:
                                # Search genes
                                search_input = st.text_input(
                                    "Search genes:",
                                    placeholder="Enter gene names separated by commas (e.g., TP53, BRCA1, MYC)",
                                    help="Search for specific genes to highlight them on the plot",
                                    key="volcano_search"
                                )
                                search_genes = []
                                if search_input:
                                    search_genes = [g.strip().upper() for g in search_input.split(",") if g.strip()]
                                    # Match case-insensitively
                                    gene_list_upper = all_genes_df["gene"].str.upper().tolist()
                                    search_genes = [g for g in search_genes if g in gene_list_upper]
                                    # Get actual gene names with original case
                                    search_genes = all_genes_df[all_genes_df["gene"].str.upper().isin(search_genes)]["gene"].tolist()

                                    if search_genes:
                                        st.success(f"Found {len(search_genes)} gene(s): {', '.join(search_genes[:10])}" +
                                                  (f"... and {len(search_genes)-10} more" if len(search_genes) > 10 else ""))
                                    elif search_input:
                                        st.warning("No matching genes found")

                            with vol_col2:
                                # Label genes - get significant genes for suggestions
                                sig_file = Path(de["sig_genes"])
                                sig_gene_list = []
                                if sig_file.exists():
                                    sig_df_temp = pd.read_csv(sig_file)
                                    sig_gene_list = sig_df_temp["gene"].tolist()

                                label_genes = st.multiselect(
                                    "Significant genes:",
                                    options=sig_gene_list if sig_gene_list else all_genes_df["gene"].head(100).tolist(),
                                    default=[],
                                    help="Select significant genes to label on the volcano plot",
                                    key="volcano_labels"
                                )

                            # Custom genes to label (text input for genes not in dropdown)
                            custom_labels = st.text_input(
                                "Additional genes to label (comma-separated):",
                                placeholder="e.g., GAPDH, ACTB",
                                help="Enter additional gene names to label, if not in the dropdown above",
                                key="volcano_custom_labels"
                            )
                            if custom_labels:
                                custom_list = [g.strip() for g in custom_labels.split(",") if g.strip()]
                                # Match to actual genes in data
                                custom_matched = all_genes_df[all_genes_df["gene"].str.upper().isin([g.upper() for g in custom_list])]["gene"].tolist()
                                label_genes = list(set(label_genes + custom_matched))

                            # Create interactive volcano plot
                            comparison_name = de.get("comparison", "DE Analysis")
                            fig = create_interactive_volcano(
                                all_genes_df,
                                lfc_threshold=lfc_threshold,
                                padj_threshold=padj_threshold,
                                search_genes=search_genes if search_genes else None,
                                label_genes=label_genes if label_genes else None,
                                title=f"Volcano Plot: {comparison_name.replace('_', ' ')}"
                            )

                            # Display the plot
                            st.plotly_chart(fig, width="stretch", key="volcano_main")

                            # Show static version option
                            with st.expander("View static volcano plot (R-generated)"):
                                if volcano_path.exists():
                                    st.image(str(volcano_path), caption="Static Volcano Plot", width="stretch")

                        elif volcano_path.exists():
                            # Fallback to static image if no CSV
                            st.image(str(volcano_path), caption="Volcano Plot", width="stretch")

                        # Significant genes table
                        sig_file = Path(de["sig_genes"])
                        if sig_file.exists():
                            sig_df = pd.read_csv(sig_file)
                            if len(sig_df) > 0:
                                with st.expander(f"View Significant Genes ({len(sig_df)})"):
                                    display_df = sig_df[["gene", "log2FoldChange", "padj", "significance"]].copy()
                                    display_df["log2FoldChange"] = display_df["log2FoldChange"].round(3)
                                    display_df["padj"] = display_df["padj"].apply(lambda x: f"{x:.2e}")
                                    st.dataframe(display_df, width="stretch")

                        # Download buttons
                        st.markdown("**Download DE Results:**")
                        col1, col2, col3, col4, col5 = st.columns(5)

                        with col1:
                            if volcano_path.exists():
                                with open(volcano_path, "rb") as f:
                                    st.download_button(
                                        "📥 Volcano (PNG)",
                                        f,
                                        file_name="volcano_plot.png",
                                        mime="image/png"
                                    )

                        with col2:
                            volcano_pdf = Path(de["volcano_pdf"])
                            if volcano_pdf.exists():
                                with open(volcano_pdf, "rb") as f:
                                    st.download_button(
                                        "📥 Volcano (PDF)",
                                        f,
                                        file_name="volcano_plot.pdf",
                                        mime="application/pdf"
                                    )

                        with col3:
                            if volcano_path.exists():
                                try:
                                    tiff_data = png_to_tiff(volcano_path)
                                    st.download_button(
                                        "📥 Volcano (TIFF)",
                                        tiff_data,
                                        file_name="volcano_plot.tiff",
                                        mime="image/tiff"
                                    )
                                except Exception:
                                    pass

                        with col4:
                            if all_genes_file.exists():
                                with open(all_genes_file, "rb") as f:
                                    st.download_button(
                                        "📥 All Genes",
                                        f,
                                        file_name="all_genes.csv",
                                        mime="text/csv"
                                    )

                        with col5:
                            if sig_file.exists():
                                with open(sig_file, "rb") as f:
                                    st.download_button(
                                        "📥 Sig Genes",
                                        f,
                                        file_name="significant_genes.csv",
                                        mime="text/csv"
                                    )

                        # --- Gene Set Overlay on Volcano ---
                        st.divider()
                        st.subheader("7. Gene Set Overlay on Volcano Plot")

                        st.markdown("""
                        Highlight genes from a specific gene set on the volcano plot.
                        Select a GMT file (e.g., from [MSigDB](https://www.gsea-msigdb.org/gsea/msigdb/)) and choose a gene set to overlay.
                        """)

                        # GMT file selection
                        col_gmt_input, col_gmt_browse = st.columns([4, 1])

                        with col_gmt_input:
                            gmt_path_input = st.text_input(
                                "GMT file path:",
                                value=st.session_state.get("gmt_file_path") or "",
                                help="Path to a .gmt gene set file (e.g., c2.all.v2026.1.Hs.symbols.gmt)",
                                placeholder="Click 'Browse' or paste path to .gmt file...",
                                key="gmt_path_input"
                            )

                        with col_gmt_browse:
                            st.markdown("<br>", unsafe_allow_html=True)
                            if st.button("📂 Browse", key="browse_gmt"):
                                selected_gmt = open_gmt_dialog()
                                if selected_gmt:
                                    st.session_state.gmt_file_path = selected_gmt
                                    st.rerun()

                        if gmt_path_input:
                            st.session_state.gmt_file_path = gmt_path_input

                        # Load pathway names for dropdown
                        gmt_file = st.session_state.get("gmt_file_path")
                        de_dir_path = Path(de["de_dir"])
                        gsea_dir = de_dir_path / "gsea"

                        # Check if pathway names are already cached
                        pathway_names_file = gsea_dir / "pathway_names.txt"
                        pathway_names = []

                        if pathway_names_file.exists():
                            pathway_names = pathway_names_file.read_text().strip().split("\n")
                            pathway_names = [p for p in pathway_names if p.strip()]

                        if gmt_file and Path(gmt_file).exists() and not pathway_names:
                            # Load pathway names from GMT file directly
                            try:
                                with open(gmt_file, "r") as f:
                                    pathway_names = sorted([line.split("\t")[0] for line in f if line.strip()])
                            except Exception:
                                pathway_names = []

                        if pathway_names:
                            selected_geneset = st.selectbox(
                                "Select gene set to highlight:",
                                options=pathway_names,
                                index=0,
                                help="Choose a gene set from the GMT file to overlay on the volcano plot"
                            )

                            # Load gene set from GMT file and generate R plot
                            geneset_genes = []
                            if st.button("🔬 Highlight Gene Set on Volcano", type="primary"):
                                try:
                                    # Load gene set
                                    with open(gmt_file, "r") as f:
                                        for line in f:
                                            parts = line.strip().split("\t")
                                            if parts[0] == selected_geneset:
                                                geneset_genes = parts[2:]
                                                break

                                    if geneset_genes:
                                        st.session_state.geneset_overlay_genes = geneset_genes
                                        st.session_state.geneset_overlay_name = selected_geneset

                                        # Also generate R version for publication-quality download
                                        with st.spinner("Generating publication-quality plot..."):
                                            overlay_output = de_dir_path / f"volcano_overlay_{selected_geneset[:50]}.png"
                                            all_genes_csv = Path(de["all_genes"])

                                            cmd = [
                                                "Rscript",
                                                str(APP_DIR / "analysis.R"),
                                                "volcano",
                                                str(all_genes_csv),
                                                gmt_file,
                                                selected_geneset,
                                                str(overlay_output),
                                                str(lfc_threshold),
                                                str(padj_threshold),
                                                selected_geneset
                                            ]

                                            result = subprocess.run(
                                                cmd,
                                                capture_output=True,
                                                text=True,
                                                cwd=str(APP_DIR)
                                            )

                                            if result.returncode == 0:
                                                st.session_state.volcano_overlay_r = str(overlay_output)
                                            else:
                                                st.session_state.volcano_overlay_r = None

                                        st.success(f"Loaded {len(geneset_genes)} genes from {selected_geneset}")
                                    else:
                                        st.warning(f"Gene set '{selected_geneset}' not found or empty")

                                except Exception as e:
                                    st.error(f"Failed to load gene set: {e}")

                            # Display interactive overlay if gene set is loaded
                            if st.session_state.get("geneset_overlay_genes"):
                                overlay_genes = st.session_state.geneset_overlay_genes
                                overlay_name = st.session_state.get("geneset_overlay_name", "Gene Set")

                                # Load DE data for overlay plot
                                all_genes_csv = Path(de["all_genes"])
                                if all_genes_csv.exists():
                                    overlay_df = pd.read_csv(all_genes_csv)

                                    # Match genes (case-insensitive)
                                    overlay_genes_upper = [g.upper() for g in overlay_genes]
                                    matched_genes = overlay_df[overlay_df["gene"].str.upper().isin(overlay_genes_upper)]["gene"].tolist()

                                    st.info(f"Matched {len(matched_genes)} of {len(overlay_genes)} genes from **{overlay_name}**")

                                    # Create interactive volcano with gene set highlight
                                    comparison_name = de.get("comparison", "DE Analysis")
                                    fig_overlay = create_interactive_volcano(
                                        overlay_df,
                                        lfc_threshold=lfc_threshold,
                                        padj_threshold=padj_threshold,
                                        highlight_genes=matched_genes,
                                        title=f"Volcano Plot with {overlay_name[:40]}{'...' if len(overlay_name) > 40 else ''}"
                                    )

                                    st.plotly_chart(fig_overlay, width="stretch", key="volcano_overlay")

                                    # Show which genes from the set are significant
                                    sig_in_set = overlay_df[
                                        (overlay_df["gene"].str.upper().isin(overlay_genes_upper)) &
                                        (overlay_df["padj"] < padj_threshold) &
                                        (abs(overlay_df["log2FoldChange"]) > lfc_threshold)
                                    ]
                                    if len(sig_in_set) > 0:
                                        with st.expander(f"Significant genes in {overlay_name} ({len(sig_in_set)})"):
                                            display_set = sig_in_set[["gene", "log2FoldChange", "padj"]].copy()
                                            display_set["log2FoldChange"] = display_set["log2FoldChange"].round(3)
                                            display_set["padj"] = display_set["padj"].apply(lambda x: f"{x:.2e}")
                                            st.dataframe(display_set, width="stretch")

                                    # Clear overlay button
                                    if st.button("Clear overlay"):
                                        st.session_state.geneset_overlay_genes = None
                                        st.session_state.geneset_overlay_name = None
                                        st.session_state.volcano_overlay_r = None
                                        st.rerun()

                                    # R-generated publication quality version
                                    if st.session_state.get("volcano_overlay_r"):
                                        r_overlay_path = Path(st.session_state.volcano_overlay_r)
                                        if r_overlay_path.exists():
                                            st.divider()
                                            st.markdown("**View static volcano plot (R-generated):**")
                                            st.image(str(r_overlay_path), use_container_width=True)

                                            st.markdown("**Export Static Plot:**")
                                            r_col1, r_col2, r_col3 = st.columns(3)
                                            with r_col1:
                                                with open(r_overlay_path, "rb") as f:
                                                    st.download_button(
                                                        "📥 PNG",
                                                        f.read(),
                                                        file_name="volcano_overlay_publication.png",
                                                        mime="image/png",
                                                        key="overlay_r_png"
                                                    )
                                            with r_col2:
                                                r_pdf_path = r_overlay_path.with_suffix(".pdf")
                                                if r_pdf_path.exists():
                                                    with open(r_pdf_path, "rb") as f:
                                                        st.download_button(
                                                            "📥 PDF",
                                                            f.read(),
                                                            file_name="volcano_overlay_publication.pdf",
                                                            mime="application/pdf",
                                                            key="overlay_r_pdf"
                                                        )
                                            with r_col3:
                                                r_tiff_path = r_overlay_path.with_suffix(".tiff")
                                                if r_tiff_path.exists():
                                                    with open(r_tiff_path, "rb") as f:
                                                        st.download_button(
                                                            "📥 TIFF",
                                                            f.read(),
                                                            file_name="volcano_overlay_publication.tiff",
                                                            mime="image/tiff",
                                                            key="overlay_r_tiff"
                                                        )

                        elif gmt_file and not Path(gmt_file).exists():
                            st.warning(f"GMT file not found: {gmt_file}")
                        elif not gmt_file:
                            st.info("Select a GMT file above to enable gene set overlay.")

                        # --- GSEA Analysis ---
                        st.divider()
                        st.subheader("8. Gene Set Enrichment Analysis (GSEA)")

                        st.markdown("""
                        Run GSEA using **fgsea** to identify enriched gene sets from the differential expression results.
                        Uses the same GMT file selected above.
                        """)

                        if gmt_file and Path(gmt_file).exists():
                            all_genes_csv = Path(de["all_genes"])

                            if not all_genes_csv.exists():
                                st.warning("DE results file not found. Please run differential expression first.")
                            else:
                                if st.button("🔬 Run GSEA", type="primary"):
                                    with st.spinner("Running GSEA with fgsea..."):
                                        try:
                                            gsea_dir.mkdir(parents=True, exist_ok=True)

                                            cmd = [
                                                "Rscript",
                                                str(APP_DIR / "analysis.R"),
                                                "gsea",
                                                str(all_genes_csv),
                                                gmt_file,
                                                str(gsea_dir)
                                            ]

                                            result = subprocess.run(
                                                cmd,
                                                capture_output=True,
                                                text=True,
                                                cwd=str(APP_DIR)
                                            )

                                            if result.returncode != 0:
                                                st.error("GSEA failed!")
                                                st.code(result.stderr)
                                            else:
                                                st.session_state.gsea_results = {
                                                    "gsea_dir": str(gsea_dir),
                                                    "barplot": str(gsea_dir / "gsea_barplot.png"),
                                                    "all_results": str(gsea_dir / "gsea_all_results.csv"),
                                                    "top_pathways": str(gsea_dir / "gsea_top_pathways.csv"),
                                                }
                                                st.success("GSEA complete!")
                                                st.rerun()

                                        except Exception as e:
                                            st.error(f"GSEA failed: {e}")

                                # Display GSEA results
                                if st.session_state.gsea_results:
                                    gsea = st.session_state.gsea_results

                                    # GSEA barplot
                                    barplot_path = Path(gsea["barplot"])
                                    if barplot_path.exists():
                                        st.image(str(barplot_path), caption="GSEA - Top Enriched Pathways", width="stretch")

                                    # Top pathways table
                                    top_file = Path(gsea["top_pathways"])
                                    if top_file.exists():
                                        top_df = pd.read_csv(top_file)
                                        if len(top_df) > 0:
                                            # Summary metrics
                                            all_res_file = Path(gsea["all_results"])
                                            if all_res_file.exists():
                                                all_gsea_df = pd.read_csv(all_res_file)
                                                n_sig = sum(all_gsea_df["padj"] < 0.05)
                                                n_total = len(all_gsea_df)

                                                col1, col2, col3 = st.columns(3)
                                                col1.metric("Total Gene Sets Tested", f"{n_total:,}")
                                                col2.metric("Significant (padj < 0.05)", n_sig)
                                                col3.metric("Top Pathways Shown", len(top_df))

                                            with st.expander(f"View Top Pathways ({len(top_df)})"):
                                                display_cols = ["pathway", "NES", "padj", "size"]
                                                display_cols = [c for c in display_cols if c in top_df.columns]
                                                disp = top_df[display_cols].copy()
                                                if "NES" in disp.columns:
                                                    disp["NES"] = disp["NES"].round(3)
                                                if "padj" in disp.columns:
                                                    disp["padj"] = disp["padj"].apply(lambda x: f"{x:.2e}")
                                                st.dataframe(disp, width="stretch")

                                    # Download buttons
                                    st.markdown("**Download GSEA Results:**")
                                    col1, col2, col3, col4 = st.columns(4)

                                    with col1:
                                        if barplot_path.exists():
                                            with open(barplot_path, "rb") as f:
                                                st.download_button(
                                                    "📥 GSEA (PNG)",
                                                    f,
                                                    file_name="gsea_barplot.png",
                                                    mime="image/png",
                                                    key="dl_gsea_png"
                                                )

                                    with col2:
                                        gsea_pdf = Path(gsea["gsea_dir"]) / "gsea_barplot.pdf"
                                        if gsea_pdf.exists():
                                            with open(gsea_pdf, "rb") as f:
                                                st.download_button(
                                                    "📥 GSEA (PDF)",
                                                    f,
                                                    file_name="gsea_barplot.pdf",
                                                    mime="application/pdf",
                                                    key="dl_gsea_pdf"
                                                )

                                    with col3:
                                        if barplot_path.exists():
                                            try:
                                                tiff_data = png_to_tiff(barplot_path)
                                                st.download_button(
                                                    "📥 GSEA (TIFF)",
                                                    tiff_data,
                                                    file_name="gsea_barplot.tiff",
                                                    mime="image/tiff",
                                                    key="dl_gsea_tiff"
                                                )
                                            except Exception:
                                                pass

                                    with col4:
                                        all_res_path = Path(gsea["all_results"])
                                        if all_res_path.exists():
                                            with open(all_res_path, "rb") as f:
                                                st.download_button(
                                                    "📥 Results (CSV)",
                                                    f,
                                                    file_name="gsea_all_results.csv",
                                                    mime="text/csv",
                                                    key="dl_gsea_csv"
                                                )
                        else:
                            st.info("Select a GMT file in the Gene Set Overlay section above to enable GSEA.")


if __name__ == "__main__":
    main()