#!/usr/bin/env python3
"""
Phylogeny阶段：对基因家族成员进行多序列比对和系统发育分析

输入：Identify阶段的final/{prefix}.final_members.pep.fa
输出：
  - final/{prefix}.pep.msa.fa (多序列比对)
  - final/{prefix}.tree.nwk (系统发育树)
  - final/{prefix}.phylogeny_report.md (分析报告)
  - final/{prefix}.phylogeny.manifest.json (机器清单)
"""

import os
import argparse
import json
import re
import shlex
import subprocess
import sys
from datetime import datetime
from pathlib import Path
from typing import Dict, List, Any, Optional
import hashlib
import time
import shutil  # 需要添加这个导入，用于复制文件

PROJECT_ROOT = Path(__file__).resolve().parents[2]
if str(PROJECT_ROOT) not in sys.path:
    sys.path.insert(0, str(PROJECT_ROOT))

from gfa_core.gene_family_scope import read_scope, sanitize_token

CUSTOM_PHYLOGENY_ROUTES = {"custom_phylogeny", "comparative_phylogeny"}


def route_output_suffix(analysis_route: str) -> str:
    route = str(analysis_route or "").strip()
    if route in CUSTOM_PHYLOGENY_ROUTES:
        return sanitize_token(route)
    return ""


from gfa_core.gene_family_stage_context import apply_stage_context, build_stage_context, user_output_prefix_for, write_legacy_manifest_copy

# Karpathy: 添加conda环境bin目录到PATH
CONDA_PREFIX = Path(sys.executable).resolve().parent.parent
CONDA_BIN = str(CONDA_PREFIX / "bin")
if CONDA_BIN not in os.environ.get("PATH", "").split(os.pathsep):
    os.environ["PATH"] = CONDA_BIN + os.pathsep + os.environ.get("PATH", "")


def resolve_iqtree_executable() -> str:
    """Resolve the IQ-TREE executable used by this Python environment."""
    candidates = [
        Path(sys.executable).resolve().parent / "iqtree3",
        Path(CONDA_BIN) / "iqtree3",
        shutil.which("iqtree3"),
        shutil.which("iqtree2"),
        shutil.which("iqtree"),
    ]
    for candidate in candidates:
        if not candidate:
            continue
        path = Path(candidate)
        if path.exists() and os.access(path, os.X_OK):
            return str(path)
    return "iqtree3"


def get_tool_versions() -> Dict[str, str]:
    """获取MAFFT和IQ-TREE版本信息"""
    versions = {}
        
    # MAFFT版本 - 修复版本获取
    try:
        # 尝试不同的方法来获取MAFFT版本
        # MAFFT的--version可能输出到stderr
        result = subprocess.run(
            ["mafft", "--version"], 
            capture_output=True, 
            text=True, 
            check=False,
            timeout=5
        )
        
        # 检查stdout和stderr
        output = result.stdout.strip() or result.stderr.strip()
        
        if output:
            # 提取版本信息
            # MAFFT输出格式可能是:
            # 1. "MAFFT v7.505 (2023/Jul/25)"
            # 2. "v7.505" 
            # 3. 或其他格式
            lines = output.split('\n')
            for line in lines:
                line = line.strip()
                if line:
                    # 查找包含版本信息的行
                    if "MAFFT v" in line:
                        versions["mafft"] = line
                        break
                    elif line.startswith("v") and "." in line:
                        versions["mafft"] = f"MAFFT {line}"
                        break
            
            # 如果没找到特定格式，使用第一行
            if "mafft" not in versions and lines:
                versions["mafft"] = lines[0]
        else:
            versions["mafft"] = "unknown (no output)"
            
    except (FileNotFoundError, subprocess.CalledProcessError):
        versions["mafft"] = "not_found"
    except subprocess.TimeoutExpired:
        versions["mafft"] = "timeout"
    except Exception as e:
        versions["mafft"] = f"error: {str(e)[:50]}"

    # IQ-TREE版本（与实际执行命令保持一致）
    try:
        iqtree_exe = resolve_iqtree_executable()
        result = subprocess.run(
            [iqtree_exe, "--version"],
            capture_output=True, 
            text=True, 
            check=False
        )
        if result.returncode == 0:
            # 提取版本信息
            lines = result.stdout.strip().split('\n')
            for line in lines:
                if "IQ-TREE version" in line or "IQ-TREE" in line:
                    versions["iqtree"] = line.strip()
                    break
            if "iqtree" not in versions:
                versions["iqtree"] = lines[0] if lines else "unknown"
        else:
            versions["iqtree"] = f"error_{result.returncode}"
    except (FileNotFoundError, subprocess.CalledProcessError):
        versions["iqtree"] = "not_found"
    
    return versions


def compute_file_hash(filepath: Path) -> str:
    """计算文件的SHA256哈希值"""
    sha256_hash = hashlib.sha256()
    try:
        with open(filepath, "rb") as f:
            for byte_block in iter(lambda: f.read(4096), b""):
                sha256_hash.update(byte_block)
        return sha256_hash.hexdigest()
    except Exception:
        return ""




def infer_effective_prefix(workspace: Path, requested_prefix: str = "") -> str:
    """Return a stable non-hidden prefix for phylogeny outputs."""
    if requested_prefix:
        return requested_prefix

    final_dir = workspace / "final"

    for pep in sorted(final_dir.glob("Multi.*.final_members.pep.fa")):
        suffix = ".final_members.pep.fa"
        if pep.name.endswith(suffix):
            return pep.name[:-len(suffix)]

    for manifest_path in sorted(final_dir.glob("*.identify.manifest.json")):
        if manifest_path.name.startswith("."):
            continue
        try:
            data = json.loads(manifest_path.read_text(encoding="utf-8"))
            prefix = data.get("prefix")
            if prefix:
                return prefix
            gene_family = data.get("gene_family")
            if gene_family:
                return f"Multi.{gene_family}"
        except Exception:
            continue

    legacy_manifest = final_dir / ".identify.manifest.json"
    if legacy_manifest.exists():
        try:
            data = json.loads(legacy_manifest.read_text(encoding="utf-8"))
            prefix = data.get("prefix")
            if prefix:
                return prefix
            gene_family = data.get("gene_family")
            if gene_family:
                return f"Multi.{gene_family}"
        except Exception:
            pass

    raise FileNotFoundError("Unable to infer prefix: no Multi.*.final_members.pep.fa or identify manifest found")

def validate_identify_output(workspace: Path, prefix: str) -> Dict[str, Any]:
    """
    验证Identify阶段的输出文件并获取必要信息
    
    Args:
        workspace: 工作空间目录（workspace/<job_id>）
        prefix: 分析前缀，如 Ft.SPL
    
    Returns:
        包含输入文件路径和Identify manifest数据的字典
    
    Raises:
        FileNotFoundError: Identify输出不存在
        ValueError: 文件无效或格式错误
    """
    print(f"正在验证Identify阶段输出: {workspace}")
    
    # 1. 查找Identify的manifest
    identify_manifest_path = workspace / "final" / f"{prefix}.identify.manifest.json"
    single_species_protein: Optional[Path] = None
    if not identify_manifest_path.exists() and "." in prefix:
        species, family = prefix.split(".", 1)
        if species not in {"Target", "Multi"} and family:
            candidate_manifest = workspace / "final" / f"Target.{family}.identify.manifest.json"
            species_candidates = [
                workspace / "final" / f"Target.{species}.{family}.identified_members.pep.fa",
                workspace / "final" / f"Target.{species}.{family}.final_members.pep.fa",
                workspace / "identify" / species / "final" / f"{species}.final_gfmembers.pep.fa",
            ]
            single_species_protein = next((p for p in species_candidates if p.exists() and p.stat().st_size > 0), None)
            if candidate_manifest.exists() and single_species_protein:
                identify_manifest_path = candidate_manifest
    if not identify_manifest_path.exists():
        raise FileNotFoundError(f"Identify manifest未找到: {identify_manifest_path}")
    
    print(f"  找到Identify manifest: {identify_manifest_path}")
    
    # 2. 读取Identify manifest
    with open(identify_manifest_path, 'r') as f:
        identify_manifest = json.load(f)
    
    # 3. 验证Identify状态
    if identify_manifest.get("status") != "success":
        raise ValueError(f"Identify阶段未成功完成，状态: {identify_manifest.get('status')}")
    
    # 4. 查找蛋白序列文件
    # artifacts可能是列表（旧格式）或字典（新格式）
    artifacts = identify_manifest.get("artifacts", {})
    protein_artifact = None
    
    if isinstance(artifacts, list):
        # 列表格式：查找artifact为"final_proteins"的项
        for artifact in artifacts:
            if isinstance(artifact, dict) and artifact.get("artifact") == "final_proteins":
                protein_artifact = artifact
                break
    elif isinstance(artifacts, dict):
        # 字典格式：直接查找"final_proteins"键
        protein_artifact = artifacts.get("final_proteins")
    
    if not protein_artifact:
        # 尝试其他可能的键名
        for key in ["identified_members", "final_proteins", "final_members", "proteins"]:
            if key in artifacts:
                protein_artifact = artifacts[key] if isinstance(artifacts, dict) else None
                break
    
    if single_species_protein:
        protein_path = single_species_protein
        protein_artifact = {
            "full_path": str(single_species_protein),
            "artifact": "species_final_members",
        }
    elif not protein_artifact:
        raise ValueError("在Identify manifest中找不到final_proteins artifact")

    # 获取文件路径
    if not single_species_protein:
        protein_path = manifest_artifact_path(workspace, identify_manifest_path, protein_artifact)
        if not protein_path:
            raw_full_path = Path(protein_artifact.get("full_path", ""))
            protein_path = raw_full_path if raw_full_path.exists() and inside_current_job(workspace, raw_full_path) else Path("")
        if not protein_path.exists():
            # 尝试相对路径
            rel_path = protein_artifact.get("path", "")
            if rel_path:
                protein_path = resolve_job_path(workspace, rel_path, identify_manifest_path) or (workspace / "final" / Path(rel_path).name)
    
    if not protein_path.exists():
        raise FileNotFoundError(f"蛋白序列文件未找到: {protein_path}")
    
    # 检查文件是否为空
    if protein_path.stat().st_size == 0:
        raise ValueError(f"蛋白序列文件为空: {protein_path}")
    
    print(f"  找到蛋白序列文件: {protein_path}")
    print(f"  文件大小: {protein_path.stat().st_size} 字节")
    
    return {
        "protein_fasta": protein_path,
        "identify_manifest": identify_manifest,
        "identify_metrics": identify_manifest.get("metrics", {}),
        "identify_parameters": identify_manifest.get("parameters", {})
    }


def read_fasta_records(fasta_file: Path) -> List[Dict[str, str]]:
    records: List[Dict[str, str]] = []
    current_header = ""
    current_seq: List[str] = []
    with fasta_file.open("r", encoding="utf-8", errors="ignore") as handle:
        for raw in handle:
            line = raw.strip()
            if not line:
                continue
            if line.startswith(">"):
                if current_header:
                    records.append({
                        "id": current_header.split()[0],
                        "header": current_header,
                        "seq": "".join(current_seq),
                    })
                current_header = line[1:].strip()
                current_seq = []
            else:
                current_seq.append(line)
    if current_header:
        records.append({
            "id": current_header.split()[0],
            "header": current_header,
            "seq": "".join(current_seq),
        })
    return records


def write_fasta_records(records: List[Dict[str, str]], output_file: Path) -> None:
    output_file.parent.mkdir(parents=True, exist_ok=True)
    with output_file.open("w", encoding="utf-8") as handle:
        for record in records:
            handle.write(f">{record['header']}\n")
            seq = record["seq"].replace("*", "")
            for idx in range(0, len(seq), 80):
                handle.write(seq[idx:idx + 80] + "\n")


def display_path(path: Path | str) -> str:
    """Return a stable user-facing absolute path without expanding symlinks."""
    return str(Path(path).absolute())


def current_input_dir(workspace: Path) -> Path:
    return workspace.parent.parent / "input" / workspace.name


def inside_current_job(workspace: Path, path: Path) -> bool:
    input_dir = current_input_dir(workspace)
    try:
        resolved = path.resolve()
        return resolved.is_relative_to(workspace.resolve()) or resolved.is_relative_to(input_dir.resolve())
    except OSError:
        return False


def rebase_current_job_path(workspace: Path, path: Path) -> Optional[Path]:
    """Map copied legacy workspace/input paths back onto the current job roots."""
    parts = path.parts
    for root_name, current_root in (("workspace", workspace), ("input", current_input_dir(workspace))):
        for idx, part in enumerate(parts):
            if part != root_name:
                continue
            if idx + 1 >= len(parts) or parts[idx + 1] != workspace.name:
                continue
            suffix = Path(*parts[idx + 2:]) if idx + 2 < len(parts) else Path()
            candidate = current_root / suffix
            if candidate.exists():
                return candidate
    input_dir = current_input_dir(workspace)
    by_name_candidates = [
        input_dir / path.name,
        workspace / "final" / path.name,
        workspace / "phylogeny" / "analysis" / path.name,
    ]
    for candidate in by_name_candidates:
        if candidate.exists():
            return candidate
    return None


def resolve_job_path(workspace: Path, raw_path: str | Path, manifest_path: Optional[Path] = None) -> Optional[Path]:
    """Resolve a copied job path without trusting stale absolute full_path values."""
    raw_text = str(raw_path or "").strip()
    if not raw_text:
        return None
    path = Path(raw_text)
    if path.is_absolute():
        if inside_current_job(workspace, path):
            return path
        rebased = rebase_current_job_path(workspace, path)
        if rebased:
            return rebased
        return None

    candidates: List[Path] = [workspace / path, current_input_dir(workspace) / path]
    if manifest_path:
        candidates.extend([manifest_path.parent / path, manifest_path.parent / path.name])
    if str(path).startswith("final/") and manifest_path:
        candidates.append(manifest_path.parent / path.name)
    for candidate in candidates:
        if candidate.exists():
            return candidate
    return None


def manifest_artifact_path(workspace: Path, manifest_path: Path, artifact: Dict[str, Any]) -> Optional[Path]:
    if not isinstance(artifact, dict):
        return None
    for key in ("path", "relpath_or_hint", "file", "full_path"):
        resolved = resolve_job_path(workspace, artifact.get(key) or "", manifest_path)
        if resolved:
            return resolved
    return None


def rebase_job_path_text(workspace: Path, value: str) -> str:
    """Rebase old job paths embedded in manifest strings or command history."""
    if not value:
        return value
    if value.startswith("/"):
        direct = rebase_current_job_path(workspace, Path(value))
        if direct:
            return display_path(direct)
    if "workspace/" + workspace.name not in value and "input/" + workspace.name not in value:
        return value

    def replace_match(match: re.Match[str]) -> str:
        raw_path = match.group(0)
        trailing = ""
        while raw_path and raw_path[-1] in {",", ";", ")", "]", "}"}:
            trailing = raw_path[-1] + trailing
            raw_path = raw_path[:-1]
        rebased = resolve_job_path(workspace, raw_path)
        return (display_path(rebased) if rebased else raw_path) + trailing

    return re.sub(r"/[^\s\"']*(?:workspace|input)/" + re.escape(workspace.name) + r"/[^\s\"']*", replace_match, value)


def rebase_manifest_payload_paths(workspace: Path, value: Any) -> Any:
    if isinstance(value, dict):
        return {key: rebase_manifest_payload_paths(workspace, item) for key, item in value.items()}
    if isinstance(value, list):
        return [rebase_manifest_payload_paths(workspace, item) for item in value]
    if isinstance(value, str):
        return rebase_job_path_text(workspace, value)
    return value


def first_existing_scope_file(workspace: Path, scope_data: Dict[str, Any], species: str, file_types: List[str]) -> Optional[Path]:
    species_files = scope_data.get("files", {}).get(species, {})
    for file_type in file_types:
        for record in species_files.get(file_type, []) or []:
            path_text = record.get("path") or ""
            if path_text:
                path = resolve_job_path(workspace, path_text) or Path(path_text)
                if path.exists() and path.stat().st_size > 0:
                    return path
    return None


def infer_input_dir(workspace: Path, scope_file: Optional[Path]) -> Path:
    """Resolve the job input directory without requiring the caller to pass it."""
    candidates: List[Path] = []
    if scope_file:
        # Common generated scope path:
        # workspace/<job_id>/final/<prefix>.<scope_id>.scope.json
        try:
            scope_parent = scope_file.parent
            if scope_parent.name == "final":
                job_workspace = scope_parent.parent
                candidates.append(job_workspace.parent.parent / "input" / job_workspace.name)
        except Exception:
            pass
    candidates.extend([
        workspace.parent.parent / "input" / workspace.name,
        workspace / "input",
    ])
    for candidate in candidates:
        if candidate.exists():
            return candidate
    return candidates[0]


def find_species_family_pep(
    workspace: Path,
    input_dir: Path,
    species: str,
    family: str,
    scope_data: Dict[str, Any],
) -> Optional[Path]:
    """Find a known family protein FASTA for reference/comparative species."""
    family_tokens = [family] if family else []
    prefix_tokens = [f"{species}.{token}" for token in family_tokens]
    candidates: List[Path] = []
    for prefix_token in prefix_tokens:
        candidates.extend([
            workspace / "final" / f"{prefix_token}.pep.fa",
            workspace / "final" / f"{prefix_token}.pep.fasta",
            workspace / "final" / f"{prefix_token}.faa",
            workspace / "final" / f"{prefix_token}.sequence_extract.pep.fa",
            workspace / "final" / f"{prefix_token}.sequence_extract.pep.fasta",
            workspace / "final" / f"{prefix_token}.ncbi_fetch.protein.fa",
            workspace / "final" / f"{prefix_token}.ncbi_fetch.protein.fasta",
            input_dir / f"{prefix_token}.pep.fa",
            input_dir / f"{prefix_token}.pep.fasta",
            input_dir / f"{prefix_token}.faa",
            input_dir / f"{prefix_token}.sequence_extract.pep.fa",
            input_dir / f"{prefix_token}.sequence_extract.pep.fasta",
            input_dir / f"{prefix_token}.ncbi_fetch.protein.fa",
            input_dir / f"{prefix_token}.ncbi_fetch.protein.fasta",
            input_dir / species / f"{prefix_token}.pep.fa",
            input_dir / species / f"{prefix_token}.pep.fasta",
            input_dir / species / f"{prefix_token}.faa",
            input_dir / species / f"{prefix_token}.sequence_extract.pep.fa",
            input_dir / species / f"{prefix_token}.sequence_extract.pep.fasta",
            input_dir / species / f"{prefix_token}.ncbi_fetch.protein.fa",
            input_dir / species / f"{prefix_token}.ncbi_fetch.protein.fasta",
        ])

    for candidate in candidates:
        if candidate.exists() and candidate.stat().st_size > 0:
            return candidate
    scoped = first_existing_scope_file(workspace, scope_data, species, ["family_pep"])
    if scoped:
        return scoped
    return None


def find_refined_species_members(
    workspace: Path,
    input_dir: Path,
    species: str,
    family: str,
    scope_data: Dict[str, Any],
) -> Optional[Path]:
    """Resolve the target species member FASTA for scoped phylogeny.

    The subfamily evidence route is pre-final and must prefer curated_members.
    Regular gene-family routes remain anchored on final_members.
    """
    family = str(family or "").strip()
    candidates: List[Path] = []
    scope_prefix = str(scope_data.get("prefix") or "").strip()
    analysis_route = str(scope_data.get("analysisRoute") or scope_data.get("analysis_route") or "").strip()
    stems = ["curated_members"] if analysis_route == "subfamily_evidence" else ["final_members"]

    for stem in stems:
        if family:
            candidates.extend([
                workspace / "final" / f"Target.{species}.{family}.{stem}.pep.fa",
                input_dir / "final" / f"Target.{species}.{family}.{stem}.pep.fa",
            ])
        # In custom multi-species scopes the scope prefix can be an aggregate or a
        # stale user-facing file. Prefer the explicit per-species Target artifact
        # above; keep the scope-prefix fallback only after that.
        if scope_prefix and scope_prefix.split(".", 1)[0] == species:
            candidates.append(workspace / "final" / f"{scope_prefix}.{stem}.pep.fa")
            candidates.append(input_dir / "final" / f"{scope_prefix}.{stem}.pep.fa")
        if family:
            candidates.extend([
                workspace / "final" / f"{species}.{family}.{stem}.pep.fa",
                input_dir / "final" / f"{species}.{family}.{stem}.pep.fa",
            ])
        candidates.extend([
            workspace / "final" / f"{species}.{stem}.pep.fa",
            input_dir / "final" / f"{species}.{stem}.pep.fa",
        ])

    seen: set[Path] = set()
    for candidate in candidates:
        try:
            resolved = candidate.absolute()
        except Exception:
            resolved = candidate
        if resolved in seen:
            continue
        seen.add(resolved)
        if candidate.exists() and candidate.stat().st_size > 0:
            return candidate
    return None


def prepare_phylogeny_scope_fasta(
    workspace: Path,
    prefix: str,
    identify_info: Dict[str, Any],
    scope_file: Optional[Path],
    fallback_fasta: Path,
    analysis_route: str = "",
) -> tuple[Path, Dict[str, Any]]:
    if not scope_file:
        records = read_fasta_records(fallback_fasta)
        return fallback_fasta, {
            "enabled": False,
            "input_fasta": display_path(fallback_fasta),
            "sequence_count": len(records),
        }

    scope_data = read_scope(scope_file)
    if analysis_route and not (scope_data.get("analysisRoute") or scope_data.get("analysis_route")):
        scope_data["analysisRoute"] = analysis_route
        scope_data["analysis_route"] = analysis_route
    scope_id = scope_data.get("scope_id", "")
    selected_species = scope_data.get("selected_species") or scope_data.get("target_species") or []
    if not selected_species:
        selected_species = identify_info.get("identify_manifest", {}).get("target_species", [])

    identify_manifest = identify_info.get("identify_manifest", {})
    artifacts = identify_manifest.get("artifacts", {}) or {}
    identify_species = identify_manifest.get("target_species", []) or []
    reference_species = scope_data.get("reference_species") or []
    comparative_species = scope_data.get("comparative_species") or []
    family = str(scope_data.get("family") or prefix.split(".")[-1] or "").strip()
    input_dir = infer_input_dir(workspace, scope_file)
    route = str(analysis_route or scope_data.get("analysisRoute") or scope_data.get("analysis_route") or "").strip()
    target_species_set = {str(item).strip() for item in scope_data.get("target_species") or [] if str(item).strip()}

    records: List[Dict[str, str]] = []
    seen_ids: set[str] = set()
    species_counts: Dict[str, int] = {}
    sources: Dict[str, str] = {}
    missing_species: List[str] = []
    missing_details: Dict[str, List[str]] = {}

    def add_records(species: str, fasta: Path, source: str) -> None:
        added = 0
        for record in read_fasta_records(fasta):
            record_id = record["id"]
            if record_id in seen_ids:
                continue
            seen_ids.add(record_id)
            records.append(record)
            added += 1
        species_counts[species] = species_counts.get(species, 0) + added
        sources[species] = display_path(fasta)
        print(f"  Scope FASTA add {species}: {added} sequences from {source}: {fasta}")

    for species in selected_species:
        refined_fasta = find_refined_species_members(
            workspace=workspace,
            input_dir=input_dir,
            species=species,
            family=family,
            scope_data=scope_data,
        )
        if refined_fasta:
            source_label = "curated_members" if route == "subfamily_evidence" else "final_members"
            add_records(species, refined_fasta, source_label)
            continue

        if route == "subfamily_evidence" and species in target_species_set:
            missing_species.append(species)
            missing_details[species] = [
                str(workspace / "final" / f"Target.{species}.{family}.curated_members.pep.fa"),
                str(input_dir / "final" / f"Target.{species}.{family}.curated_members.pep.fa"),
                str(workspace / "final" / f"{species}.{family}.curated_members.pep.fa"),
                "请先完成 identify 成员初筛确认，生成 curated_members 后再运行亚家族证据进化树。",
            ]
            continue

        artifact = artifacts.get(f"{species}_final_members") if isinstance(artifacts, dict) else None
        artifact_path = manifest_artifact_path(workspace, workspace / "final" / f"{prefix}.identify.manifest.json", artifact) if artifact else None
        if artifact_path and artifact_path.exists() and artifact_path.stat().st_size > 0:
            add_records(species, artifact_path, "identify")
            continue

        if species in identify_species:
            fallback = workspace / "identify" / species / "final" / f"{species}.final_gfmembers.pep.fa"
            if fallback.exists() and fallback.stat().st_size > 0:
                add_records(species, fallback, "identify_fallback")
                continue

        if species in reference_species or species in comparative_species:
            family_file = find_species_family_pep(
                workspace=workspace,
                input_dir=input_dir,
                species=species,
                family=family,
                scope_data=scope_data,
            )
            if family_file:
                add_records(species, family_file, "family_pep")
                continue

        missing_species.append(species)
        missing_details[species] = [
            f"identify artifact: {species}_final_members",
            str(workspace / "identify" / species / "final" / f"{species}.final_gfmembers.pep.fa"),
            str(input_dir / f"{species}.{family}.pep.fa"),
            str(input_dir / f"{species}.{family}.pep.fasta"),
            str(input_dir / species / f"{species}.{family}.pep.fasta"),
            str(workspace / "final" / f"{species}.{family}.pep.fa"),
            str(workspace / "final" / f"{species}.{family}.sequence_extract.pep.fa"),
            str(workspace / "final" / f"{species}.{family}.ncbi_fetch.protein.fa"),
        ]

    if missing_species:
        raise RuntimeError(
            "Phylogeny scope 缺少以下物种的家族蛋白序列: "
            + ",".join(missing_species)
            + "。请先运行对应物种 identify，或提供 Species.Family.pep.fa。"
            + " 已检查路径: "
            + json.dumps(missing_details, ensure_ascii=False)
        )
    if not records:
        raise RuntimeError("Phylogeny scope 没有可用于建树的蛋白序列。")

    output_fasta = workspace / "phylogeny" / "analysis" / f"{prefix}.{sanitize_token(scope_id)}.input.pep.fa"
    write_fasta_records(records, output_fasta)
    return output_fasta, {
        "enabled": True,
        "scope_file": display_path(scope_file),
        "scope_id": scope_id,
        "scope_mode": scope_data.get("scope_mode", ""),
        "family": scope_data.get("family", ""),
        "target_species": scope_data.get("target_species") or [],
        "reference_species": reference_species,
        "comparative_species": comparative_species,
        "selected_species": selected_species,
        "input_fasta": display_path(output_fasta),
        "sequence_count": len(records),
        "species_counts": species_counts,
        "sources": sources,
        "input_dir": display_path(input_dir) if input_dir.exists() else str(input_dir),
    }


def prepare_directories(workspace: Path) -> Dict[str, Path]:
    """创建Phylogeny阶段的标准目录结构"""
    paths = {
        "analysis_dir": workspace / "phylogeny" / "analysis",
        "final_dir": workspace / "final",
        "logs_dir": workspace / "logs",
    }
    
    print("创建目录结构:")
    for name, path in paths.items():
        path.mkdir(parents=True, exist_ok=True)
        print(f"  {name}: {path}")
    
    return paths


def force_option_value(args: list[str], option: str, value: str) -> list[str]:
    normalized: list[str] = []
    skip_next = False
    prefix = f"{option}="
    for arg in args:
        if skip_next:
            skip_next = False
            continue
        if arg == option:
            skip_next = True
            continue
        if arg.startswith(prefix):
            continue
        normalized.append(arg)
    normalized.extend([option, value])
    return normalized


def run_mafft_alignment(
    input_fasta: Path,
    output_msa: Path,
    cpu: int = 8,
    mafft_params: Optional[str] = None
) -> Dict[str, Any]:
    """
    运行MAFFT进行多序列比对
    
    Args:
        input_fasta: 输入FASTA文件
        output_msa: 输出MSA文件
        cpu: CPU核心数
        mafft_params: 自定义MAFFT参数
    
    Returns:
        包含命令、返回码、输出和指标的字典
    """
    print(f"运行MAFFT比对: {input_fasta} -> {output_msa}")
    
    if mafft_params:
        params = force_option_value(shlex.split(mafft_params), "--thread", str(cpu))
        cmd = ["mafft"] + params
        print(f"  使用自定义参数: {' '.join(params)}")
    else:
        cmd = ["mafft", "--auto", "--thread", str(cpu)]
        print(f"  使用默认参数: --auto --thread {cpu}")

    cmd.append(str(input_fasta))
    
    # 记录开始时间
    start_time = time.time()
    
    # 执行比对
    try:
        with open(output_msa, "w") as f_out:
            result = subprocess.run(
                cmd, 
                stdout=f_out, 
                stderr=subprocess.PIPE, 
                text=True,
                check=False
            )
        
        elapsed = time.time() - start_time
        
        if result.returncode == 0:
            print(f"  MAFFT成功完成，耗时: {elapsed:.1f}秒")
            
            # 计算MSA统计
            metrics = compute_msa_metrics(output_msa)
            metrics["alignment_time_sec"] = elapsed
            
            return {
                "command": " ".join(cmd),
                "returncode": result.returncode,
                "stderr": result.stderr,
                "metrics": metrics,
                "output_file": output_msa,
                "success": True
            }
        else:
            print(f"  MAFFT失败，返回码: {result.returncode}")
            if result.stderr:
                print(f"  错误信息: {result.stderr[:200]}")
            
            return {
                "command": " ".join(cmd),
                "returncode": result.returncode,
                "stderr": result.stderr,
                "metrics": {},
                "output_file": None,
                "success": False
            }
            
    except Exception as e:
        print(f"  MAFFT执行异常: {e}")
        return {
            "command": " ".join(cmd),
            "returncode": -1,
            "stderr": str(e),
            "metrics": {},
            "output_file": None,
            "success": False
        }


def compute_msa_metrics(msa_file: Path) -> Dict[str, Any]:
    """计算MSA的统计指标"""
    try:
        sequences = []
        current_seq = ""
        current_id = ""
        
        with open(msa_file, 'r') as f:
            for line in f:
                line = line.strip()
                if line.startswith('>'):
                    if current_id:
                        sequences.append(current_seq)
                    current_id = line[1:].split()[0] if line[1:] else ""
                    current_seq = ""
                elif line:
                    current_seq += line
            
            if current_id:
                sequences.append(current_seq)
        
        if not sequences:
            return {
                "sequences_count": 0,
                "msa_length": 0,
                "gap_percentage": 0.0,
                "min_length": 0,
                "max_length": 0
            }
        
        # 检查序列长度一致性
        lengths = [len(seq) for seq in sequences]
        if len(set(lengths)) > 1:
            print(f"警告: MSA序列长度不一致: {lengths}")
        
        # 计算gap百分比
        seq_count = len(sequences)
        msa_length = max(lengths) if lengths else 0
        total_chars = seq_count * msa_length
        
        if total_chars > 0:
            gap_chars = sum(seq.count('-') + seq.count('.') for seq in sequences)
            gap_percentage = (gap_chars / total_chars * 100)
        else:
            gap_percentage = 0.0
        
        return {
            "sequences_count": seq_count,
            "msa_length": msa_length,
            "gap_percentage": round(gap_percentage, 2),
            "min_length": min(lengths) if lengths else 0,
            "max_length": max(lengths) if lengths else 0
        }
        
    except Exception as e:
        print(f"计算MSA指标时出错: {e}")
        return {
            "sequences_count": 0,
            "msa_length": 0,
            "gap_percentage": 0.0,
            "min_length": 0,
            "max_length": 0
        }


def count_tree_tips(tree_file: Path) -> int:
    try:
        from Bio import Phylo

        tree = Phylo.read(str(tree_file), "newick")
        return len(tree.get_terminals())
    except Exception as exc:
        print(f"计算tree tip数量失败: {exc}")
        return 0

def run_iqtree_tree(
    msa_file: Path,
    output_dir: Path,
    prefix: str,
    cpu: int = 8,
    iqtree_params: Optional[str] = None
) -> Dict[str, Any]:
    """Run IQ-TREE and return the generated tree file."""
    print(f"Run IQ-TREE: {msa_file}")
    output_dir.mkdir(parents=True, exist_ok=True)
    base_name = output_dir / f"{prefix}.iqtree"

    # IQ-TREE reuses checkpoints for the same -pre prefix. In this pipeline the
    # same scope prefix can be rerun after upstream member IDs or sequences are
    # corrected, so stale treefiles must not be accepted as current output.
    for stale in output_dir.glob(f"{prefix}.iqtree*"):
        if stale.is_file():
            try:
                stale.unlink()
            except OSError as exc:
                print(f"  Warning: failed to remove stale IQ-TREE output {stale.name}: {exc}")

    if iqtree_params:
        params = force_option_value(shlex.split(iqtree_params), "-nt", str(cpu))
        cmd = [resolve_iqtree_executable(), "-s", str(msa_file), "-pre", str(base_name)] + params
        print(f"  Custom IQ-TREE params: {' '.join(params)}")
    else:
        cmd = [resolve_iqtree_executable(), "-s", str(msa_file), "-pre", str(base_name), "-m", "MFP", "-bb", "1000", "-nt", str(cpu)]
        print(f"  Default IQ-TREE params: -m MFP -bb 1000 -nt {cpu}")

    start_time = time.time()
    try:
        result = subprocess.run(cmd, capture_output=True, text=True, check=False)
        elapsed = time.time() - start_time
        stdout = result.stdout or ""
        stderr = result.stderr or ""

        tree_candidates = [
            base_name.with_suffix(".treefile"),
            output_dir / f"{prefix}.iqtree.treefile",
            base_name.with_suffix(".contree"),
            output_dir / f"{prefix}.iqtree.contree",
        ]
        tree_file = next((p for p in tree_candidates if p.exists() and p.stat().st_size > 0), None)
        best_model = extract_best_model(stdout + "\n" + stderr)

        checkpoint_finished = "previous run successfully finished" in (stdout + stderr)
        success = tree_file is not None and (result.returncode == 0 or checkpoint_finished)
        if success:
            print(f"  IQ-TREE completed, elapsed: {elapsed:.1f}s")
        else:
            print(f"  IQ-TREE completed, elapsed: {elapsed:.1f}s")
            if stderr:
                print(f"  Error: {stderr[:200]}")

        return {
            "command": " ".join(cmd),
            "returncode": result.returncode,
            "stdout": stdout,
            "stderr": stderr,
            "tree_file": tree_file,
            "best_model": best_model,
            "tree_building_time_sec": elapsed,
            "success": success
        }
    except Exception as e:
        print(f"  IQ-TREE exception: {e}")
        return {
            "command": " ".join(cmd),
            "returncode": -1,
            "stdout": "",
            "stderr": str(e),
            "tree_file": None,
            "best_model": "",
            "tree_building_time_sec": 0,
            "success": False
        }

def run_msa_plotting(
    msa_file: Path,
    output_prefix: str,
    method: str = "mafft",
    rscript_path: Optional[Path] = None
) -> Dict[str, Any]:
    """
    运行R脚本绘制多序列比对图
    
    Args:
        msa_file: 比对好的MSA文件
        output_prefix: 输出文件前缀
        method: 比对方法（用于报告）
        rscript_path: R脚本路径，如果None则自动查找
    
    Returns:
        包含绘图结果的字典
    """
    print(f"📊 绘制MSA比对图: {msa_file}")
    
    # 自动查找R脚本
    if rscript_path is None:
        # 尝试多个可能的路径
        script_candidates = [
            Path(__file__).parent.parent / "scripts" / "msa_figure_from_aligned.R",
            Path(__file__).parent.parent.parent / "scripts" / "msa_figure_from_aligned.R",
            Path.cwd() / "scripts" / "msa_figure_from_aligned.R",
        ]
        
        for candidate in script_candidates:
            if candidate.exists():
                rscript_path = candidate
                break
        
        if rscript_path is None:
            # 如果专用脚本不存在，尝试原始脚本
            if method == "mafft":
                rscript_path = Path(__file__).parent.parent / "scripts" / "msa_figure.mafft.R"
            else:
                rscript_path = Path(__file__).parent.parent / "scripts" / "msa_figure.R"
        # ========== 修复点1: 提前定义 output_dir ==========
    output_dir = None
    output_prefix_path = None

    # 确保输出目录存在
    if rscript_path and rscript_path.exists():
        output_prefix_path = Path(output_prefix)
        output_dir = output_prefix_path.parent
        output_dir.mkdir(parents=True, exist_ok=True)
    # ========== 修复结束 ==========
    
    if not rscript_path or not rscript_path.exists():
        return {
            "success": False,
            "error": f"R脚本未找到: {rscript_path}",
            "output_files": []
        }
    
    print(f"  使用R脚本: {rscript_path}")
    
    # 构建R命令
    rscript_exe = Path(shutil.which("Rscript") or "Rscript")
    if not rscript_exe.exists():
        rscript_exe = Path("Rscript")

    cmd = [
        str(rscript_exe),
        str(rscript_path),
        display_path(msa_file),
       # str(output_prefix),
        str(output_prefix_path.name),  # 只传递文件名，不含路径
        "--method", method
    ]
    
    # 记录开始时间
    start_time = time.time()
    
    try:
        result = subprocess.run(
            cmd,
            capture_output=True,
            text=True,
            cwd=output_dir,  # 关键：在输出目录中运行R脚本
            timeout=600,  # 10分钟超时
            check=False
        )
        
        elapsed = time.time() - start_time
        
        if result.returncode == 0:
            print(f"  ✅ MSA绘图成功，耗时: {elapsed:.1f}秒")
            
            # Collect generated image files without constructing invalid glob patterns.
            output_dir = output_prefix_path.parent
            output_stem = output_prefix_path.name
            output_files = []
            for f in output_dir.iterdir():
                if f.is_file() and f.suffix.lower() in {".png", ".pdf"} and f.name.startswith(output_stem):
                    output_files.append(str(f))
            output_files = sorted(set(output_files))
            
            print(f"  生成 {len(output_files)} 个图片文件在: {output_dir}")
            for i, f in enumerate(output_files[:3], 1):
                print(f"    {i}. {Path(f).name}")


            return {
                "success": True,
                "command": " ".join(cmd),
                "returncode": result.returncode,
                "stdout": result.stdout[-500:],  # 只保留最后500字符
                "stderr": result.stderr[-500:] if result.stderr else "",
                "output_files": output_files,
                "elapsed_sec": elapsed
            }
        else:
            print(f"  ⚠️  MSA绘图失败，返回码: {result.returncode}")
            if result.stderr:
                print(f"  错误信息: {result.stderr[:200]}")
            
            return {
                "success": False,
                "command": " ".join(cmd),
                "returncode": result.returncode,
                "stdout": result.stdout[-500:],
                "stderr": result.stderr[-500:] if result.stderr else "",
                "output_files": [],
                "elapsed_sec": elapsed
            }
            
    except subprocess.TimeoutExpired:
        print(f"  ❌ MSA绘图超时 (>10分钟)")
        return {
            "success": False,
            "error": "timeout",
            "output_files": []
        }
    except Exception as e:
        print(f"  ❌ MSA绘图异常: {e}")
        return {
            "success": False,
            "error": str(e),
            "output_files": []
        }



def run_msa_conserved_region_plotting(
    workspace: Path,
    msa_file: Path,
    prefix: str,
    window_len: int = 60,
    min_nongap: float = 0.6,
    scope_file: str = "",
    analysis_route: str = "",
    member_set_id: str = "",
) -> Dict[str, Any]:
    """Generate an additional automatically selected conserved-region MSA plot."""
    script = PROJECT_ROOT / "cli" / "run_msa_region_plot.py"
    python_exe = Path(sys.executable)
    cmd = [
        str(python_exe),
        str(script),
        "--workspace", display_path(workspace),
        "--prefix", prefix,
        "--msa-file", display_path(msa_file),
        "--auto-conserved",
        "--window-len", str(window_len),
        "--min-nongap", str(min_nongap),
        "--force",
    ]
    if scope_file:
        cmd.extend(["--scope-file", str(scope_file)])
    if analysis_route:
        cmd.extend(["--analysis-route", str(analysis_route)])
    if member_set_id:
        cmd.extend(["--member-set-id", str(member_set_id)])
    start_time = time.time()
    try:
        result = subprocess.run(
            cmd,
            capture_output=True,
            text=True,
            cwd=str(PROJECT_ROOT),
            timeout=900,
            check=False,
        )
        elapsed = time.time() - start_time
        final_dir = workspace / "final"
        output_files = []
        if final_dir.exists():
            for f in final_dir.iterdir():
                if f.is_file() and f.name.startswith(f"{prefix}.msa_conserved_region_"):
                    output_files.append(str(f))
        output_files = sorted(set(output_files))
        success = result.returncode == 0 and bool(output_files)
        return {
            "success": success,
            "command": " ".join(cmd),
            "returncode": result.returncode,
            "stdout": result.stdout[-1000:],
            "stderr": result.stderr[-1000:] if result.stderr else "",
            "output_files": output_files,
            "elapsed_sec": elapsed,
        }
    except subprocess.TimeoutExpired:
        return {"success": False, "error": "timeout", "output_files": []}
    except Exception as exc:
        return {"success": False, "error": str(exc), "output_files": []}



def infer_species_family_from_output_prefix(prefix: str) -> tuple[str, str]:
    """Infer Species/Family from output labels such as Lc.SPL.*."""
    parts = [p for p in str(prefix or "").split(".") if p]
    if len(parts) >= 2:
        return parts[0], parts[1]
    return "", ""


def find_cross_scope_subfamily_file(final_dir: Path, prefix: str, input_dir: Path) -> Optional[Path]:
    """Resolve confirmed labels produced for a target+reference scope.

    Prediction tables are not auto-detected here. A caller that wants to plot
    unconfirmed prediction labels must pass them explicitly through --subfamily.
    """
    species, family = infer_species_family_from_output_prefix(prefix)
    candidates: list[Path] = []
    if species and family:
        stem = f"{species}.{family}.target_plus_reference"
        candidates.extend([
            final_dir / f"{stem}.confirmed_subfamily.tsv",
            input_dir / f"{stem}.confirmed_subfamily.tsv",
        ])
        for pattern in (
            f"{species}.{family}*target_plus_reference*confirmed_subfamily.tsv",
            f"*target-{species}*target_plus_reference*confirmed_subfamily.tsv",
        ):
            candidates.extend(sorted(final_dir.glob(pattern)))
    for candidate in candidates:
        if candidate.exists() and candidate.stat().st_size > 0:
            return candidate
    return None


def run_tree_plotting(
    tree_file: Path,
    output_prefix: Optional[str] = None,
    rscript_path: Optional[Path] = None,
    subfamily_file: Optional[Path] = None,
    disable_subfamily_auto: bool = False,
    tree_layout: str = "both",
    label_offset: float = 0.4,
    tip_label_size: float = 3.0,
    branch_color: str = "black",
    branch_width: float = 0.5,
    show_bootstrap: bool = True,
    bootstrap_size: float = 2.0,
    tip_point_size: float = 2.0,
    target_shape: str = "triangle",
    reference_shape: str = "circle",
    gene_shape: str = "circle",
    marker_mode: str = "auto",
    circular_label_offset: float = 0.8,
    circular_label_size: float = 1.8,
    subfamily_arc_enabled: bool = True,
    subfamily_arc_gap: float = 1.2,
    subfamily_arc_label_gap: float = 0.6,
    subfamily_arc_label_size: float = 0.0,
) -> Dict[str, Any]:
    """Run the R tree plotting script and return the concrete files it generated."""
    print(f"Plot tree: {tree_file}")
    tree_layout = (tree_layout or "both").lower()
    if tree_layout not in {"rectangular", "circular", "both"}:
        return {"success": False, "error": f"Invalid tree layout: {tree_layout}", "output_files": []}

    if output_prefix is None:
        output_prefix = tree_file.stem

    prefix = tree_file.stem.replace('.tree', '').replace('.nwk', '')

    if subfamily_file is None and not disable_subfamily_auto:
        final_dir = tree_file.parent
        workspace_dir = final_dir.parent
        job_id = workspace_dir.name
        project_root = workspace_dir.parent.parent
        input_dir = project_root / "input" / job_id
        potential_confirmed_final = final_dir / f"{prefix}.confirmed_subfamily.tsv"
        potential_confirmed_input = input_dir / f"{prefix}.confirmed_subfamily.tsv"
        potential_cross_scope = find_cross_scope_subfamily_file(final_dir, prefix, input_dir)

        print(f"  [DEBUG] job_id: {job_id}")
        print(f"  [DEBUG] project_root: {project_root}")
        print(f"  [DEBUG] input_dir: {input_dir}")
        print(f"  [DEBUG] lookup confirmed final: {potential_confirmed_final} (exists: {potential_confirmed_final.exists()})")
        print(f"  [DEBUG] lookup confirmed input: {potential_confirmed_input} (exists: {potential_confirmed_input.exists()})")
        print(f"  [DEBUG] lookup cross-scope confirmed subfamily: {potential_cross_scope} (exists: {bool(potential_cross_scope and potential_cross_scope.exists())})")

        if potential_confirmed_final.exists():
            subfamily_file = potential_confirmed_final
            print(f"  Using confirmed subfamily from final: {subfamily_file}")
        elif potential_confirmed_input.exists():
            subfamily_file = potential_confirmed_input
            print(f"  Using confirmed subfamily from input: {subfamily_file}")
        elif potential_cross_scope and potential_cross_scope.exists():
            subfamily_file = potential_cross_scope
            print(f"  Using cross-scope confirmed subfamily: {subfamily_file}")
        else:
            print(f"  Confirmed subfamily file not found: {prefix}.confirmed_subfamily.tsv")

    if subfamily_file and subfamily_file.exists():
        script_name = "tree_plot.R"
        use_subfamily = True
    else:
        script_name = "plot_tree.R"
        use_subfamily = False

    if rscript_path is None:
        script_candidates = [
            Path(__file__).parent.parent / "scripts" / script_name,
            Path(__file__).parent.parent.parent / "scripts" / script_name,
            Path.cwd() / "scripts" / script_name,
        ]
        for candidate in script_candidates:
            if candidate.exists():
                rscript_path = candidate
                break

    if not rscript_path or not rscript_path.exists():
        return {"success": False, "error": f"R script not found: {rscript_path}", "output_files": []}

    print(f"  R script: {rscript_path}")
    print(f"  Mode: {'subfamily' if use_subfamily else 'normal'}")
    print(f"  Layout: {tree_layout}")

    tree_dir = tree_file.parent
    style_args_base = [
        str(label_offset),
        str(tip_label_size),
        str(branch_color or "black"),
        str(branch_width),
        "true" if show_bootstrap else "false",
        str(bootstrap_size),
        str(tip_point_size),
    ]
    if use_subfamily:
        subfamily_plot_args = [
            *style_args_base,
            str(target_shape or "triangle"),
            str(reference_shape or "circle"),
            str(gene_shape or target_shape or "triangle"),
            str(circular_label_offset),
            str(circular_label_size),
            str(marker_mode or "auto"),
            "true" if subfamily_arc_enabled else "false",
            str(subfamily_arc_gap),
            str(subfamily_arc_label_gap),
            str(subfamily_arc_label_size if subfamily_arc_label_size and subfamily_arc_label_size > 0 else ""),
        ]
        cmd = ["Rscript", str(rscript_path), str(tree_file.name), display_path(subfamily_file), str(output_prefix), tree_layout, "", "", *subfamily_plot_args]
    else:
        normal_plot_args = [
            *style_args_base,
            str(gene_shape or target_shape or "triangle"),
            str(circular_label_offset),
            str(circular_label_size),
            str(marker_mode or "auto"),
        ]
        cmd = ["Rscript", str(rscript_path), str(tree_file.name), str(tree_dir), tree_layout, str(output_prefix or ""), *normal_plot_args]

    start_time = time.time()
    try:
        result = subprocess.run(cmd, capture_output=True, text=True, cwd=tree_dir, timeout=600, check=False)
        elapsed = time.time() - start_time

        if result.returncode == 0:
            print(f"  Tree plotting completed, elapsed: {elapsed:.1f}s")
            output_files = []
            if use_subfamily:
                expected_names = [
                    f"{output_prefix}.png",
                    f"{output_prefix}.pdf",
                    f"{output_prefix}.circular_tree.png",
                    f"{output_prefix}.circular_tree.pdf",
                ]
                for name in expected_names:
                    f = tree_dir / name
                    if f.exists():
                        output_files.append(str(f))
                        print(f"    输出: {f.name}")
            else:
                normal_stem = str(output_prefix or tree_file.stem)
                expected_names = [
                    f"{normal_stem}.png",
                    f"{normal_stem}.pdf",
                    f"{normal_stem}.circular_tree.png",
                    f"{normal_stem}.circular_tree.pdf",
                ]
                for name in expected_names:
                    f = tree_dir / name
                    if f.exists():
                        output_files.append(str(f))
                        print(f"    输出: {f.name}")
                output_dir = tree_dir / f"{tree_file.stem}_output"
                if output_dir.exists():
                    extra_files = [str(f) for f in output_dir.glob("*") if f.is_file()]
                    output_files.extend(extra_files)
                    print(f"  Generated {len(extra_files)} files in: {output_dir}")
            output_files = sorted(set(output_files))
            return {
                "success": True,
                "command": " ".join(cmd),
                "returncode": result.returncode,
                "stdout": result.stdout[-500:],
                "stderr": result.stderr[-500:] if result.stderr else "",
                "output_files": output_files,
                "output_dir": str(tree_dir),
                "elapsed_sec": elapsed,
                "mode": "subfamily" if use_subfamily else "normal",
                "subfamily_file": str(subfamily_file) if subfamily_file else None,
                "plot_params": {
                    "label_offset": label_offset,
                    "tip_label_size": tip_label_size,
                    "branch_color": branch_color,
                    "branch_width": branch_width,
                    "show_bootstrap": show_bootstrap,
                    "bootstrap_size": bootstrap_size,
                    "tip_point_size": tip_point_size,
                    "target_shape": target_shape,
                    "reference_shape": reference_shape,
                    "gene_shape": gene_shape,
                    "marker_mode": marker_mode,
                    "circular_label_offset": circular_label_offset,
                    "circular_label_size": circular_label_size,
                    "subfamily_arc_enabled": subfamily_arc_enabled,
                    "subfamily_arc_gap": subfamily_arc_gap,
                    "subfamily_arc_label_gap": subfamily_arc_label_gap,
                    "subfamily_arc_label_size": subfamily_arc_label_size,
                }
            }

        print(f"  Tree plotting failed, return code: {result.returncode}")
        if result.stderr:
                print(f"  Error: {stderr[:200]}")
        return {
            "success": False,
            "command": " ".join(cmd),
            "returncode": result.returncode,
            "stdout": result.stdout[-500:],
            "stderr": result.stderr[-500:] if result.stderr else "",
            "output_files": [],
            "elapsed_sec": elapsed
        }
    except subprocess.TimeoutExpired:
        print("  ? ??????? (>10??)")
        return {"success": False, "error": "timeout", "output_files": []}
    except Exception as e:
        print(f"  ? ???????: {e}")
        return {"success": False, "error": str(e), "output_files": []}

def extract_best_model(iqtree_output: str) -> str:
    """从IQ-TREE输出中提取最佳模型"""
    if not iqtree_output:
        return "unknown"
    
    lines = iqtree_output.split('\n')
    for line in lines:
        if "Best-fit model:" in line:
            return line.split("Best-fit model:")[1].strip()
        elif "Model of substitution:" in line:
            return line.split("Model of substitution:")[1].strip()
        elif "Best model:" in line:
            return line.split("Best model:")[1].strip()
    
    return "unknown"


def write_report(
    workspace: Path,
    prefix: str,
    metrics: Dict[str, Any],
    parameters: Dict[str, Any],
    tool_versions: Dict[str, str],
    mafft_result: Dict[str, Any],
    iqtree_result: Dict[str, Any],
    msa_plot_result: Optional[Dict[str, Any]],  # 新增
    tree_plot_result: Optional[Dict[str, Any]],  # 新增
    artifacts: Dict[str, Dict]
) -> Path:
    """写入人可读的分析报告"""
    report_path = workspace / "final" / f"{prefix}.phylogeny_report.md"
    
    print(f"生成分析报告: {report_path}")
    
    with open(report_path, "w") as f:
        f.write(f"# Phylogeny分析报告 - {prefix}\n\n")
        f.write(f"**生成时间**: {datetime.now().isoformat()}\n\n")
        
        f.write("## 1. 分析概述\n")
        f.write(f"- **分析阶段**: Phylogeny\n")
        f.write(f"- **分析状态**: {'成功' if iqtree_result.get('success', False) else '部分成功'}\n")
        f.write(f"- **序列数量**: {metrics.get('sequences_count', 'N/A')}\n")
        f.write(f"- **比对长度**: {metrics.get('msa_length', 'N/A')} 氨基酸\n")
        f.write(f"- **Gap比例**: {metrics.get('gap_percentage', 'N/A')}%\n")
        f.write(f"- **最佳模型**: {iqtree_result.get('best_model', 'N/A')}\n\n")
        
        f.write("## 2. 工具版本\n")
        for tool, version in tool_versions.items():
            f.write(f"- **{tool}**: {version}\n")
        f.write("\n")
        
        f.write("## 3. 执行参数\n")
        for param, value in parameters.items():
            f.write(f"- **{param}**: {value}\n")
        f.write("\n")
        
        f.write("## 4. 执行命令\n")
        f.write("### MAFFT比对\n")
        f.write(f"```bash\n{mafft_result.get('command', 'N/A')}\n```\n\n")
        
        if mafft_result.get('returncode', -1) != 0:
            f.write(f"**返回码**: {mafft_result.get('returncode', 'N/A')}\n")
            if mafft_result.get('stderr'):
                f.write(f"**错误信息**:\n```\n{mafft_result['stderr'][:500]}\n```\n\n")
        
        f.write("### IQ-TREE建树\n")
        f.write(f"```bash\n{iqtree_result.get('command', 'N/A')}\n```\n\n")
        
        if iqtree_result.get('returncode', -1) != 0:
            f.write(f"**返回码**: {iqtree_result.get('returncode', 'N/A')}\n")
            if iqtree_result.get('stderr'):
                f.write(f"**错误信息**:\n```\n{iqtree_result['stderr'][:500]}\n```\n\n")
        ##开始 #  
                # MSA绘图结果
        f.write("### 多序列比对图\n")
        if msa_plot_result and msa_plot_result.get("success", False):
            f.write("✅ **MSA比对图绘制成功**\n")
            f.write(f"- 耗时: {msa_plot_result.get('elapsed_sec', 0):.1f}秒\n")

            # 列出生成的图片
            msa_images = [k for k, v in artifacts.items() if isinstance(v, dict) and v.get("artifact") == "msa_plot"]
            if msa_images:
                f.write("- 生成图片:\n")
                for img_key in msa_images[:5]:  # 最多显示5个
                    img_info = artifacts[img_key]
                    f.write(f"  - `{img_info.get('path', '')}` ({img_info.get('size', 0)}字节)\n")
        else:
            f.write("❌ **MSA比对图绘制失败**\n")
            if msa_plot_result:
                f.write(f"- 错误: {msa_plot_result.get('error', '未知错误')}\n")

        f.write("\n")

        # 进化树绘图结果
        f.write("### 系统发育树图\n")
        if tree_plot_result and tree_plot_result.get("success", False):
            f.write("✅ **进化树图绘制成功**\n")
            f.write(f"- 耗时: {tree_plot_result.get('elapsed_sec', 0):.1f}秒\n")

            # 列出生成的图片
            tree_images = [k for k, v in artifacts.items() if isinstance(v, dict) and v.get("artifact") == "tree_plot"]
            if tree_images:
                f.write("- 生成图片:\n")
                for img_key in tree_images[:5]:  # 最多显示5个
                    img_info = artifacts[img_key]
                    f.write(f"  - `{img_info.get('path', '')}` ({img_info.get('size', 0)}字节)\n")
        else:
            f.write("❌ **进化树图绘制失败**\n")
            if tree_plot_result:
                f.write(f"- 错误: {tree_plot_result.get('error', '未知错误')}\n")

        f.write("\n")
        ###结束


        f.write("## 5. 产出文件\n")
        for artifact_name, artifact_info in artifacts.items():
            if isinstance(artifact_info, dict):
                f.write(f"- **{artifact_name}**: {artifact_info.get('path', 'N/A')}\n")
                if artifact_info.get('size'):
                    f.write(f"  - 大小: {artifact_info['size']} 字节\n")
                if artifact_info.get('sha256'):
                    f.write(f"  - SHA256: {artifact_info['sha256'][:16]}...\n")
        f.write("\n")
        
        f.write("## 6. 详细统计\n")
        f.write("```json\n")
        f.write(json.dumps(metrics, indent=2, ensure_ascii=False))
        f.write("\n```\n")
    
    return report_path


def write_manifest(
    workspace: Path,
    job_id: str,
    run_id: str,
    prefix: str,
    start_time: datetime,
    end_time: datetime,
    artifacts: Dict[str, Dict],
    metrics: Dict[str, Any],
    parameters: Dict[str, Any],
    tool_versions: Dict[str, str],
    cmd_history: List[str],
    identify_inputs: List[Dict],
    status: str = "success",
    scope_info: Optional[Dict[str, Any]] = None,
    output_prefix: Optional[str] = None,
    analysis_route: str = "",
    member_set_id: str = "",
    subfamily_file: str = "",
) -> Path:
    """
    写入符合schema v1的manifest文件
    
    注意：artifacts必须是字典格式，不是列表
    """
    elapsed_sec = (end_time - start_time).total_seconds()
    output_prefix = output_prefix or prefix
    scope_id = str((scope_info or {}).get("scope_id") or "").strip()
    
    manifest_data = {
        # 基础标识字段（必需）
        "manifest_version": "1.0",
        "job_id": job_id,
        "stage": "phylogeny",
        "run_id": run_id,
        "prefix": prefix,
        "output_prefix": output_prefix,
        
        # 状态和时间字段（必需）
        "timestamp": end_time.isoformat(),
        "status": status,
        "start_time": start_time.isoformat(),
        "end_time": end_time.isoformat(),
        "elapsed_sec": round(elapsed_sec, 2),
        
        # 执行信息（必需）
        "cmd": rebase_manifest_payload_paths(workspace, cmd_history),
        "tool_versions": tool_versions,
        
        # 目录路径（必需）
        "workspace_root": display_path(workspace.parent),
        "workspace": display_path(workspace / "final"),
        "input_dir": display_path(workspace.parent.parent / "input" / job_id),
        
        # 输入文件（必需）
        "inputs": rebase_manifest_payload_paths(workspace, identify_inputs),
        
        # 输出产物（必需）- 必须是字典格式
        "artifacts": artifacts,
        
        # 参数和指标（必需）
        "parameters": rebase_manifest_payload_paths(workspace, parameters),
        "metrics": metrics,
        "scope_id": scope_id,
        "scope": rebase_manifest_payload_paths(workspace, scope_info or {}),
        
        # 溯源信息（可选但推荐）
        "provenance": {
            "derived_from": f"{prefix}.identify.manifest.json",
            "stage_sequence": ["identify", "phylogeny"]
        }
    }
    context = build_stage_context(
        stage="phylogeny",
        prefix=prefix,
        scope_data=scope_info or {},
        scope_file=(scope_info or {}).get("scope_file", ""),
        output_prefix=output_prefix,
        subfamily_file=subfamily_file,
        analysis_route=analysis_route,
        member_set_id=member_set_id,
    )
    apply_stage_context(manifest_data, context)
    
    if scope_id:
        manifest_path = workspace / "final" / f"{output_prefix}.phylogeny.manifest.json"
        manifest_data["scope_manifest"] = display_path(manifest_path)
    else:
        manifest_path = workspace / "final" / f"{prefix}.phylogeny.manifest.json"
    
    print(f"生成Manifest文件: {manifest_path}")
    
    # 确保JSON严格性：无注释，确保ASCII
    with open(manifest_path, "w", encoding='utf-8') as f:
        json.dump(manifest_data, f, indent=2, ensure_ascii=False)
    
    return manifest_path


def main():
    """主函数：执行Phylogeny分析流程"""
    parser = argparse.ArgumentParser(
        description="Phylogeny阶段：多序列比对和系统发育分析",
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog="""
示例：
  python run_phylogeny.py --prefix Ft.SPL --workspace /path/to/workspace --cpu 4
  
输入要求：
  - Identify阶段必须已成功运行
  - final/{prefix}.final_members.pep.fa 必须存在
  
输出产物：
  - final/{prefix}.pep.msa.fa (多序列比对)
  - final/{prefix}.tree.nwk (系统发育树)
  - final/{prefix}.phylogeny_report.md (分析报告)
  - final/{prefix}.phylogeny.manifest.json (机器清单)
        """
    )
    
    parser.add_argument("--prefix", required=True, 
                       help="分析前缀，格式：物种.家族（如 Ft.SPL）")
    parser.add_argument("--workspace", required=True,
                       help="工作空间目录（workspace/<job_id>/)")
    parser.add_argument("--cpu", type=int, default=8,
                       help="CPU核心数（默认：8）")
    parser.add_argument("--proteins", "--protein-fasta",
                       dest="proteins",
                       help="Protein FASTA override; defaults to Identify manifest final_members protein FASTA")
    parser.add_argument("--scope-file",
                       help="Scope JSON; when provided, build phylogeny input FASTA from selected species")
    parser.add_argument("--analysis-route", default="", help="Explicit analysis route identity for manifest/status")
    parser.add_argument("--member-set-id", default="", help="Explicit member-set identity for manifest/status")
    parser.add_argument("--subfamily", default="", help="Optional subfamily table for tree plotting")
    parser.add_argument("--mafft-params",
                       help="MAFFT参数覆盖（默认：--auto --thread {cpu}）")
    parser.add_argument("--iqtree-params",
                       help="IQ-TREE参数覆盖（默认：-m MFP -bb 1000 -nt {cpu}；-nt 会被 --cpu 覆盖）")
    parser.add_argument("--tree-layout", choices=["rectangular", "circular", "both"], default="both",
                       help="进化树绘图布局：rectangular、circular 或 both（默认：both）")
    parser.add_argument("--tree-subfamily-arc-enabled", choices=["true", "false"], default="true",
                       help="圈图是否显示亚家族外圈弧线（默认：true）")
    parser.add_argument("--tree-subfamily-arc-gap", type=float, default=1.2,
                       help="圈图亚家族弧线相对基因ID外边界的额外距离（默认：1.2）")
    parser.add_argument("--tree-subfamily-arc-label-gap", type=float, default=0.6,
                       help="圈图亚家族名称相对弧线的额外距离（默认：0.6）")
    parser.add_argument("--tree-subfamily-arc-label-size", type=float, default=0.0,
                       help="圈图亚家族名称字号；0 表示自动（默认：0）")
    parser.add_argument("--log-file",
                       help="日志文件路径（默认：workspace/logs/{timestamp}_phylogeny.log）")
    
    args = parser.parse_args()
    
    # 记录开始时间
    start_time = datetime.now()
    run_id = f"phylogeny_{start_time.strftime('%Y%m%d_%H%M%S_%f')[:-3]}"
    
    # 设置工作空间
    workspace = Path(args.workspace)
    job_id = workspace.name  # workspace/<job_id>中的job_id
    effective_prefix = infer_effective_prefix(workspace, args.prefix)
    
    print("=" * 60)
    print(f"🚀 启动Phylogeny分析")
    print("=" * 60)
    print(f"分析前缀: {args.prefix}")
    print(f"工作空间: {workspace}")
    print(f"任务ID: {job_id}")
    print(f"运行ID: {run_id}")
    print(f"CPU核心: {args.cpu}")
    print(f"树图布局: {args.tree_layout}")
    
    try:
        # 1. 准备目录结构
        print("\n📁 准备目录结构...")
        paths = prepare_directories(workspace)
        
        # 2. 获取工具版本
        print("\n🔧 检查工具版本...")
        tool_versions = get_tool_versions()
        print(f"   MAFFT: {tool_versions.get('mafft', 'unknown')}")
        print(f"   IQ-TREE: {tool_versions.get('iqtree', 'unknown')}")
        
        # 3. 验证Identify输出
        print("\n🔍 验证Identify阶段输出...")
        identify_info = validate_identify_output(workspace, effective_prefix)
        protein_fasta = identify_info["protein_fasta"]
        scope_path = Path(args.scope_file).absolute() if args.scope_file else None
        if scope_path and not scope_path.exists():
            raise FileNotFoundError(f"Scope file not found: {scope_path}")
        scope_info = {
            "enabled": False,
            "input_fasta": display_path(protein_fasta),
            "sequence_count": len(read_fasta_records(protein_fasta)),
        }
        if args.proteins:
            protein_override = resolve_job_path(workspace, args.proteins) or Path(args.proteins)
            if not protein_override.exists():
                raise FileNotFoundError(f"Specified protein FASTA not found: {protein_override}")
            if protein_override.stat().st_size == 0:
                raise ValueError(f"Specified protein FASTA is empty: {protein_override}")
            protein_fasta = protein_override
            print(f"  使用指定蛋白序列文件: {protein_fasta}")
            scope_info = {
                "enabled": bool(scope_path),
                "scope_file": str(scope_path) if scope_path else "",
                "input_fasta": display_path(protein_fasta),
                "sequence_count": len(read_fasta_records(protein_fasta)),
                "note": "--proteins override was used; scope selected species were not used to build FASTA",
            }
        elif scope_path:
            print("\n🧭 准备Phylogeny scope输入FASTA...")
            protein_fasta, scope_info = prepare_phylogeny_scope_fasta(
                workspace=workspace,
                prefix=effective_prefix,
                identify_info=identify_info,
                scope_file=scope_path,
                fallback_fasta=protein_fasta,
                analysis_route=args.analysis_route,
            )
            print(f"  Scope输入FASTA: {protein_fasta}")
            print(f"  Scope序列数: {scope_info.get('sequence_count')}")
        output_prefix_token = effective_prefix
        if scope_info.get("scope_id"):
            output_prefix_token = user_output_prefix_for(
                "phylogeny",
                effective_prefix,
                scope_info,
                analysis_route=args.analysis_route,
                member_set_id=args.member_set_id,
            )
        print(f"  输出前缀: {output_prefix_token}")
        
        # 准备Identify的inputs信息
        identify_inputs = []
        if "inputs" in identify_info["identify_manifest"]:
            identify_inputs = identify_info["identify_manifest"]["inputs"]
        
        # 4. 运行MAFFT比对
        print("\n🔄 运行MAFFT多序列比对...")
        msa_intermediate = paths["analysis_dir"] / f"{output_prefix_token}.pep.mafft.msa.fa"
        mafft_result = run_mafft_alignment(
            input_fasta=protein_fasta,
            output_msa=msa_intermediate,
            cpu=args.cpu,
            mafft_params=args.mafft_params
        )
        
        if not mafft_result.get("success", False):
            print("❌ MAFFT比对失败")
            raise RuntimeError(f"MAFFT失败: {mafft_result.get('stderr', '未知错误')}")
        
        msa_metrics = mafft_result.get("metrics", {})
        print(f"  比对成功: {msa_metrics.get('sequences_count', '未知')} 条序列")
        
        # 5. 复制MSA到final目录
        msa_final = paths["final_dir"] / f"{output_prefix_token}.pep.msa.fa"
        msa_final.write_bytes(msa_intermediate.read_bytes())
        print(f"  MSA保存到: {msa_final}")
        
        # 6. 运行IQ-TREE建树
        print("\n🌳 运行IQ-TREE构建系统发育树...")
        iqtree_result = run_iqtree_tree(
            msa_file= paths["final_dir"] / f"{output_prefix_token}.pep.msa.fa",
            output_dir=paths["analysis_dir"],
            prefix=output_prefix_token,
            cpu=args.cpu,
            iqtree_params=args.iqtree_params
        )
        
        # 7. 复制树文件到final目录（即使部分成功也尝试）
        tree_final = None
        if iqtree_result.get("tree_file") and iqtree_result["tree_file"].exists():
            tree_final = paths["final_dir"] / f"{output_prefix_token}.tree.nwk"
            tree_final.write_bytes(iqtree_result["tree_file"].read_bytes())
            print(f"  树文件保存到: {tree_final}")
        elif iqtree_result.get("success", False):
            # 尝试查找其他可能的树文件
            tree_files = list(paths["analysis_dir"].glob(f"{output_prefix_token}.iqtree.*"))
            tree_files = [f for f in tree_files if f.suffix in ['.treefile', '.nwk', '.tre']]
            if tree_files:
                tree_final = paths["final_dir"] / f"{output_prefix_token}.tree.nwk"
                tree_final.write_bytes(tree_files[0].read_bytes())
                print(f"  树文件保存到: {tree_final}")
            else:
                print("⚠️  警告: 未找到树文件，但IQ-TREE报告成功")


          # ========== 8. 绘制MSA比对图 ==========
        msa_plot_result = None
        msa_conserved_region_result = None
        tree_plot_result = None
        
        # 确保MSA文件路径正确
        msa_final = paths["final_dir"] / f"{output_prefix_token}.pep.msa.fa"
        
        if msa_final.exists():
            print("\n" + "="*60)
            print("📸 开始绘制MSA比对图...")
            msa_plot_result = run_msa_plotting(
                msa_file=msa_final,
                output_prefix=paths["final_dir"] / f"{output_prefix_token}.msa_plot",
                method="mafft"
            )
            
            if msa_plot_result and msa_plot_result.get("success", False):
                print(f"✅ MSA比对图绘制完成")
            else:
                print(f"⚠️  MSA比对图绘制失败")

            print("\n📌 自动筛选保守MSA区域并绘图...")
            msa_conserved_region_result = run_msa_conserved_region_plotting(
                workspace=workspace,
                msa_file=msa_final,
                prefix=output_prefix_token,
                window_len=60,
                min_nongap=0.6,
                scope_file=scope_info.get("scope_file", ""),
                analysis_route=args.analysis_route,
                member_set_id=args.member_set_id,
            )
            if msa_conserved_region_result and msa_conserved_region_result.get("success", False):
                print("✅ 保守区MSA图绘制完成")
            else:
                print("⚠️  保守区MSA图绘制失败")
        
        # ========== 9. 绘制进化树图 ==========
        if tree_final and tree_final.exists():
            print("\n" + "="*60)
            print("🌳 开始绘制进化树图...")
            subfamily_arg = str(args.subfamily or "").strip()
            subfamily_disabled = subfamily_arg.lower() in {"none", "skip", "disabled"}
            subfamily_override = None if (not subfamily_arg or subfamily_disabled) else (resolve_job_path(workspace, subfamily_arg) or Path(subfamily_arg)).absolute()
            tree_plot_result = run_tree_plotting(
                tree_file=tree_final,
                output_prefix=f"{output_prefix_token}.tree",
                tree_layout=args.tree_layout,
                subfamily_file=subfamily_override,
                disable_subfamily_auto=subfamily_disabled,
                subfamily_arc_enabled=str(args.tree_subfamily_arc_enabled).lower() != "false",
                subfamily_arc_gap=args.tree_subfamily_arc_gap,
                subfamily_arc_label_gap=args.tree_subfamily_arc_label_gap,
                subfamily_arc_label_size=args.tree_subfamily_arc_label_size,
            )
            
            if tree_plot_result and tree_plot_result.get("success", False):
                print(f"✅ 进化树图绘制完成")
            else:
                print(f"⚠️  进化树图绘制失败")
       


       # 8. 确定分析状态
        status = "success"
        if not iqtree_result.get("success", False):
            if mafft_result.get("success", False):
                status = "partial_success"
                print("⚠️  分析部分成功: MAFFT成功但IQ-TREE失败")
            else:
                status = "failed"
        
        # 9. 准备artifacts信息（字典格式）
        artifacts = {}
        
        if msa_final.exists():
            artifacts["msa"] = {
                "artifact": "msa",
                "path": f"final/{output_prefix_token}.pep.msa.fa",
                "full_path": display_path(msa_final),
                "type": "fasta",
                "size": msa_final.stat().st_size,
                "sha256": compute_file_hash(msa_final),
                "mtime": msa_final.stat().st_mtime
            }
        
        if tree_final and tree_final.exists():
            artifacts["tree"] = {
                "artifact": "tree",
                "path": f"final/{output_prefix_token}.tree.nwk",
                "full_path": display_path(tree_final),
                "type": "newick",
                "size": tree_final.stat().st_size,
                "sha256": compute_file_hash(tree_final),
                "mtime": tree_final.stat().st_mtime
            }
       
                # 添加MSA绘图文件到artifacts
        if msa_plot_result and msa_plot_result.get("success", False):
            
            # for i, img_file in enumerate(msa_plot_result.get("output_files", []), 1):
            #     img_path = Path(img_file)
            #     if img_path.exists():
            png_added = False
            for img_file in msa_plot_result.get("output_files", []):
                img_path = Path(img_file)
                if not img_path.exists():
                    continue

                # 只处理 PNG 图片用于 Gallery 显示
                if img_path.suffix.lower() == ".png":
                    artifact_key = f"msa_plot_{img_path.stem}"
                    artifacts[artifact_key] = {
                        "artifact": "msa_plot",
                        "path": str(img_path.relative_to(workspace)),
                        "full_path": display_path(img_path),
                        "type": "image",  # 只有 PNG 才是 image
                        "format": "png",
                        "size": img_path.stat().st_size,
                        "sha256": compute_file_hash(img_path),
                        "caption": "MSA Alignment"
                    }
                    png_added = True
                    print(f"[DEBUG] Added MSA PNG artifact: {img_path.name}")
                
                # PDF 作为 document 类型，不进入 Gallery
                elif img_path.suffix.lower() == ".pdf":
                    artifact_key = f"msa_plot_{img_path.stem}_pdf"
                    artifacts[artifact_key] = {
                        "artifact": "msa_plot_pdf",
                        "path": str(img_path.relative_to(workspace)),
                        "full_path": display_path(img_path),
                        "type": "document",  # 不是 image！
                        "format": "pdf",
                        "size": img_path.stat().st_size,
                        "sha256": compute_file_hash(img_path)
                    }
                    print(f"[DEBUG] Added MSA PDF artifact: {img_path.name}")
            
            if not png_added:
                print(f"[WARN] No MSA PNG file found in output_files")
                    # artifact_key = f"msa_plot_{i}"
                    # if img_path.suffix == ".png":
                    #     artifact_key += "_png"
                    # elif img_path.suffix == ".pdf":
                    #     artifact_key += "_pdf"
                    
                    # artifacts[artifact_key] = {
                    #     "artifact": "msa_plot",
                    #     "path": str(img_path.relative_to(workspace)),
                    #     "full_path": str(img_path),
                    #     "type": "image",
                    #     "format": img_path.suffix[1:],  # 去掉点号
                    #     "size": img_path.stat().st_size,
                    #     "sha256": compute_file_hash(img_path)
                    # }
        
        if msa_conserved_region_result and msa_conserved_region_result.get("success", False):
            for output_file in msa_conserved_region_result.get("output_files", []):
                conserved_path = Path(output_file)
                if not conserved_path.exists():
                    continue
                suffix = conserved_path.suffix.lower()
                if suffix == ".png":
                    key = f"msa_conserved_region_plot_{conserved_path.stem}"
                    artifacts[key] = {
                        "artifact": "msa_conserved_region_plot",
                        "path": str(conserved_path.relative_to(workspace)),
                        "full_path": display_path(conserved_path),
                        "type": "image",
                        "format": "png",
                        "size": conserved_path.stat().st_size,
                        "sha256": compute_file_hash(conserved_path),
                        "caption": "Automatically selected conserved MSA region",
                    }
                elif suffix == ".pdf":
                    key = f"msa_conserved_region_plot_{conserved_path.stem}_pdf"
                    artifacts[key] = {
                        "artifact": "msa_conserved_region_plot_pdf",
                        "path": str(conserved_path.relative_to(workspace)),
                        "full_path": display_path(conserved_path),
                        "type": "document",
                        "format": "pdf",
                        "size": conserved_path.stat().st_size,
                        "sha256": compute_file_hash(conserved_path),
                    }
                elif suffix in {".tsv", ".fa", ".json"}:
                    key = f"msa_conserved_region_{conserved_path.stem}_{suffix[1:]}"
                    artifacts[key] = {
                        "artifact": "msa_conserved_region_file",
                        "path": str(conserved_path.relative_to(workspace)),
                        "full_path": display_path(conserved_path),
                        "type": "file",
                        "format": suffix[1:],
                        "size": conserved_path.stat().st_size,
                        "sha256": compute_file_hash(conserved_path),
                    }

        # Add tree plot files returned by run_tree_plotting only.
        if tree_plot_result and tree_plot_result.get("success", False):
            for output_file in tree_plot_result.get("output_files", []):
                plot_path = Path(output_file)
                if not plot_path.exists():
                    continue
                if plot_path.suffix.lower() == ".png":
                    key = f"tree_plot_{plot_path.stem}"
                    artifacts[key] = {
                        "artifact": "tree_plot",
                        "path": str(plot_path.relative_to(workspace)),
                        "full_path": display_path(plot_path),
                        "type": "image",
                        "format": "png",
                        "size": plot_path.stat().st_size,
                        "sha256": compute_file_hash(plot_path),
                        "caption": "Tree Plot"
                    }
                elif plot_path.suffix.lower() == ".pdf":
                    key = f"tree_plot_{plot_path.stem}_pdf"
                    artifacts[key] = {
                        "artifact": "tree_plot_pdf",
                        "path": str(plot_path.relative_to(workspace)),
                        "full_path": display_path(plot_path),
                        "type": "document",
                        "format": "pdf",
                        "size": plot_path.stat().st_size,
                        "sha256": compute_file_hash(plot_path)
                    }
            if tree_plot_result.get("output_dir"):
                output_dir = Path(tree_plot_result["output_dir"])
                if output_dir.exists():
                    artifacts["tree_plot_output_dir"] = {
                        "artifact": "tree_plot_directory",
                        "path": str(output_dir.relative_to(workspace)),
                        "full_path": display_path(output_dir),
                        "type": "directory",
                        "file_count": len(list(output_dir.glob("*")))
                    }

        # 10. ??metrics
        combined_metrics = {}
        if msa_metrics:
            combined_metrics.update(msa_metrics)
        
        if iqtree_result.get("best_model"):
            combined_metrics["tree_model"] = iqtree_result["best_model"]
        
        if mafft_result.get("metrics", {}).get("alignment_time_sec"):
            combined_metrics["alignment_time_sec"] = mafft_result["metrics"]["alignment_time_sec"]
        
        if iqtree_result.get("tree_building_time_sec"):
            combined_metrics["tree_building_time_sec"] = iqtree_result["tree_building_time_sec"]
        if tree_final and tree_final.exists():
            combined_metrics["tree_tip_count"] = count_tree_tips(tree_final)
        
        # 添加绘图相关metrics
        if msa_plot_result:
            combined_metrics["msa_plotting_success"] = msa_plot_result.get("success", False)
            if msa_plot_result.get("elapsed_sec"):
                combined_metrics["msa_plotting_time_sec"] = msa_plot_result["elapsed_sec"]

        if msa_conserved_region_result:
            combined_metrics["msa_conserved_region_success"] = msa_conserved_region_result.get("success", False)
            if msa_conserved_region_result.get("elapsed_sec"):
                combined_metrics["msa_conserved_region_time_sec"] = msa_conserved_region_result["elapsed_sec"]
        
        if tree_plot_result:
            combined_metrics["tree_plotting_success"] = tree_plot_result.get("success", False)
            if tree_plot_result.get("elapsed_sec"):
                combined_metrics["tree_plotting_time_sec"] = tree_plot_result["elapsed_sec"]
            if tree_plot_result.get("subfamily_file"):
                combined_metrics["tree_subfamily_file"] = tree_plot_result["subfamily_file"]
                combined_metrics["tree_plot_mode"] = tree_plot_result.get("mode", "")
        combined_metrics["scope_sequence_count"] = scope_info.get("sequence_count")
        if scope_info.get("species_counts"):
            combined_metrics["scope_species_counts"] = scope_info.get("species_counts")


        # 11. 准备parameters
        parameters = {
            "cpu": args.cpu,
            "prefix": effective_prefix,
            "output_prefix": output_prefix_token,
            "mafft_params": args.mafft_params or "--auto --thread {cpu}",
            "iqtree_params": args.iqtree_params or f"-m MFP -bb 1000 -nt {args.cpu}",
            "tree_layout": args.tree_layout,
            "tree_subfamily_arc_enabled": str(args.tree_subfamily_arc_enabled).lower() != "false",
            "tree_subfamily_arc_gap": args.tree_subfamily_arc_gap,
            "tree_subfamily_arc_label_gap": args.tree_subfamily_arc_label_gap,
            "tree_subfamily_arc_label_size": args.tree_subfamily_arc_label_size,
            "mafft_version": tool_versions.get("mafft", "unknown"),
            "iqtree_version": tool_versions.get("iqtree", "unknown"),
            "stage": "phylogeny",
            "job_id": job_id,
            "run_id": run_id,
            "manifest_version": "1.0",
            "status": status,
            "scope_id": scope_info.get("scope_id", ""),
            "scope_file": scope_info.get("scope_file", ""),
            "scope_input_fasta": scope_info.get("input_fasta", ""),
            "tree_subfamily_file": tree_plot_result.get("subfamily_file") if tree_plot_result else "",
            "tree_subfamily_disabled": bool(locals().get("subfamily_disabled", False)),
            "tree_plot_mode": tree_plot_result.get("mode") if tree_plot_result else "",
        }
        
        # 12. 记录命令历史
        cmd_history = []
        if mafft_result.get("command"):
            cmd_history.append(mafft_result["command"])
        if iqtree_result.get("command"):
            cmd_history.append(iqtree_result["command"])
        
        # 13. 写入报告
        print("\n📝 生成分析报告...")
        report_path = write_report(
            workspace=workspace,
            prefix=output_prefix_token,
            metrics=combined_metrics,
            parameters=parameters,
            tool_versions=tool_versions,
            mafft_result=mafft_result,
            iqtree_result=iqtree_result,
            msa_plot_result=msa_plot_result,  # 新增
            tree_plot_result=tree_plot_result,  # 新增
            artifacts=artifacts
        )
        
        # 添加报告到artifacts
        if report_path.exists():
            artifacts["report"] = {
                "artifact": "report",
                "path": f"final/{output_prefix_token}.phylogeny_report.md",
                "full_path": display_path(report_path),
                "type": "markdown",
                "size": report_path.stat().st_size,
                "sha256": compute_file_hash(report_path)
            }
        
        # 14. 写入manifest
        print("\n📋 生成Manifest文件...")
        end_time = datetime.now()
        manifest_path = write_manifest(
            workspace=workspace,
            job_id=job_id,
            run_id=run_id,
            prefix=effective_prefix,
            start_time=start_time,
            end_time=end_time,
            artifacts=artifacts,
            metrics=combined_metrics,
            parameters=parameters,
            tool_versions=tool_versions,
            cmd_history=cmd_history,
            identify_inputs=identify_inputs,
            status=status,
            scope_info=scope_info,
            output_prefix=output_prefix_token,
            analysis_route=args.analysis_route,
            member_set_id=args.member_set_id,
            subfamily_file=tree_plot_result.get("subfamily_file") if tree_plot_result else "",
        )
        
        # 添加manifest自身到artifacts
        if manifest_path.exists():
            artifacts["manifest"] = {
                "artifact": "manifest",
                "path": str(manifest_path.relative_to(workspace)),
                "full_path": display_path(manifest_path),
                "type": "json",
                "size": manifest_path.stat().st_size,
                "sha256": compute_file_hash(manifest_path)
            }
        
        # 15. 更新manifest包含完整的artifacts
        if manifest_path.exists():
            with open(manifest_path, 'r') as f:
                manifest_data = json.load(f)
            manifest_data["artifacts"] = artifacts
            with open(manifest_path, 'w') as f:
                json.dump(manifest_data, f, indent=2, ensure_ascii=False)
            if scope_info.get("scope_id"):
                legacy_manifest = write_legacy_manifest_copy(manifest_path, effective_prefix, str(scope_info.get("scope_id") or ""), "phylogeny")
                if legacy_manifest:
                    print(f"   Legacy scoped manifest: {legacy_manifest}")
        
        # 16. 输出成功信息
        elapsed = (end_time - start_time).total_seconds()
        print("\n" + "=" * 60)
        print(f"✅ Phylogeny分析完成!")
        print(f"   总耗时: {elapsed:.1f} 秒")
        print(f"   状态: {status}")
        print(f"   Manifest: {manifest_path}")
        print(f"   报告: {report_path}")
        
        if msa_final.exists():
            print(f"   MSA: {msa_final}")
        
        if tree_final and tree_final.exists():
            print(f"   系统发育树: {tree_final}")

         # 添加绘图文件信息
        if msa_plot_result and msa_plot_result.get("success", False):
            output_files = msa_plot_result.get("output_files", [])
            print(f"   MSA比对图: {len(output_files)}个文件")
            for i, f in enumerate(output_files[:3], 1):  # 最多显示3个
                print(f"      {i}. {Path(f).name}")
        
        if tree_plot_result and tree_plot_result.get("success", False):
            output_files = tree_plot_result.get("output_files", [])
            print(f"   进化树图: {len(output_files)}个文件")
            for i, f in enumerate(output_files[:3], 1):  # 最多显示3个
                print(f"      {i}. {Path(f).name}")


        print("=" * 60)
        
        # 根据状态返回退出码
        if status == "success":
            sys.exit(0)
        elif status == "partial_success":
            sys.exit(1)  # 部分成功返回1
        else:
            sys.exit(2)  # 失败返回2
        
    except Exception as e:
        # 记录结束时间
        end_time = datetime.now()
        
        print(f"\n❌ Phylogeny分析失败: {e}")
        import traceback
        traceback.print_exc()
        
        try:
            # 创建基本的artifacts（可能部分成功）
            artifacts = {}
            failure_scope_info = locals().get("scope_info")
            if not failure_scope_info and args.scope_file:
                failure_scope_info = {
                    "enabled": True,
                    "scope_file": display_path(args.scope_file),
                }
                try:
                    scope_data = json.loads(Path(args.scope_file).read_text(encoding="utf-8"))
                    failure_scope_info.update({
                        "scope_id": scope_data.get("scope_id", ""),
                        "scope_mode": scope_data.get("scope_mode", ""),
                        "target_species": scope_data.get("target_species", []),
                        "reference_species": scope_data.get("reference_species", []),
                        "comparative_species": scope_data.get("comparative_species", []),
                        "selected_species": scope_data.get("selected_species", []),
                    })
                except Exception:
                    pass
            
            # 尝试写入失败manifest
            manifest_path = write_manifest(
                workspace=workspace,
                job_id=job_id,
                run_id=run_id,
                prefix=effective_prefix,
                start_time=start_time,
                end_time=end_time,
                artifacts=artifacts,
                metrics={},
                parameters={
                    "cpu": args.cpu, 
                    "prefix": effective_prefix,
                    "error": str(e)[:200],
                    "subfamily_file": "" if str(args.subfamily or "").strip().lower() in {"none", "skip", "disabled"} else (args.subfamily or "")
                },
                tool_versions=tool_versions if 'tool_versions' in locals() else {},
                cmd_history=[],
                identify_inputs=identify_inputs if 'identify_inputs' in locals() else [],
                status="failed",
                scope_info=failure_scope_info,
                output_prefix=locals().get("output_prefix_token", effective_prefix),
                analysis_route=args.analysis_route,
                member_set_id=args.member_set_id,
                subfamily_file="",
            )
            
            print(f"   失败Manifest已写入: {manifest_path}")
            
        except Exception as manifest_error:
            print(f"   写入失败Manifest时出错: {manifest_error}")
        
        sys.exit(2)


if __name__ == "__main__":
    main()