Files
PlantAI-dev/apps/web/pipeline_backend.py.bak_species_select
2026-08-03 18:42:47 +08:00

3155 lines
136 KiB
Plaintext
Executable File
Raw Permalink Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
# web/pipeline_backend.py
import json
import os
import re
import shutil
import subprocess
import sys
import time
import traceback
from dataclasses import dataclass
from datetime import datetime, timedelta
from pathlib import Path
from typing import Dict, List, Optional, Tuple
# 可复现性支持
sys.path.insert(0, str(Path(__file__).parent.parent.parent))
from gfa_core.reproducibility_v2 import ReproducibilityManager
# API 客户端支持
from gfa_core.api_clients import (
MultiAPIQuery,
PlantCAREClient,
PfamClient,
NCBIClient,
EnsemblPlantsClient,
)
# 你要求:上传落到项目 input/;输出在 workspace/
INPUT_ROOT = Path("input")
WORKSPACE_ROOT = Path("workspace")
INPUT_ROOT.mkdir(exist_ok=True)
WORKSPACE_ROOT.mkdir(exist_ok=True)
RETENTION_DAYS = 7
IMAGE_SUFFIX = {".png", ".jpg", ".jpeg", ".svg", ".webp"}
PDF_SUFFIX = {".pdf"}
HTML_SUFFIX = {".html", ".htm"}
def _clean_prefix(value) -> str:
if value is None:
return ""
value = str(value).strip()
return value if value and value.lower() not in {"none", "null", "undefined"} else ""
def _read_prefix_from_manifest(path: Path) -> str:
try:
data = json.loads(path.read_text(encoding="utf-8"))
except Exception:
return ""
return _clean_prefix(data.get("prefix"))
def _gene_family_from_config_xlsx(config_path: Path) -> str:
if not config_path.exists():
return ""
try:
import pandas as pd
df = pd.read_excel(config_path)
except Exception as exc:
print(f"[WARN] Failed to read Gene_Family from {config_path}: {exc}")
return ""
if "Gene_Family" not in df.columns:
return ""
for value in df["Gene_Family"].tolist():
gene_family = _clean_prefix(value)
if gene_family:
return gene_family
return ""
def resolve_effective_prefix_for_job(job, prefix: str, stage: str) -> str:
prefix = _clean_prefix(prefix)
if prefix:
return prefix
in_dir = Path(job.input_dir)
ws = Path(job.workspace_dir)
species_config = in_dir / "genome_prepare" / "species_config.json"
try:
if species_config.exists():
data = json.loads(species_config.read_text(encoding="utf-8"))
target = data.get("target") if isinstance(data, dict) else None
if isinstance(target, dict):
resolved = _clean_prefix(target.get("prefix"))
if resolved and resolved.lower() != "multi_species":
return resolved
except Exception as exc:
print(f"[WARN] Failed to read prefix from {species_config}: {exc}")
for root in [in_dir / "final", ws / "final"]:
if not root.exists():
continue
for path in sorted(root.glob("*.genome_prepare.manifest.json"), key=lambda p: p.stat().st_mtime, reverse=True):
resolved = _read_prefix_from_manifest(path)
if resolved:
return resolved
for path in sorted(root.glob("Multi.*.identify.manifest.json"), key=lambda p: p.stat().st_mtime, reverse=True):
resolved = _read_prefix_from_manifest(path) or path.name[:-len(".identify.manifest.json")]
if resolved:
return resolved
for path in sorted(root.glob("Multi.*.final_members.pep.fa")):
resolved = path.name[:-len(".final_members.pep.fa")]
if resolved:
return resolved
gene_family = _gene_family_from_config_xlsx(in_dir / "config.xlsx")
if gene_family:
return f"Multi.{gene_family}"
for root in [in_dir, in_dir / "final", ws / "final"]:
if not root.exists():
continue
for path in sorted(root.glob("Ref_combined.*.pep.fasta")):
gene_family = path.name[len("Ref_combined."):-len(".pep.fasta")]
gene_family = _clean_prefix(gene_family)
if gene_family:
return f"Multi.{gene_family}"
return ""
def _has_non_comment_data(path: Path) -> bool:
"""Return True if a text table has at least one non-empty, non-comment row."""
try:
if not path or not path.exists() or path.stat().st_size == 0:
return False
with open(path, "r", encoding="utf-8", errors="ignore") as handle:
for raw in handle:
line = raw.strip()
if line and not line.startswith("#"):
return True
except Exception:
return False
return False
def _now_tag() -> str:
return datetime.now().strftime("%Y%m%d_%H%M%S")
def cleanup_old_jobs(retention_days: int = RETENTION_DAYS) -> str:
"""清理 input/<job> 和 workspace/<job> 中超过保留期的任务目录"""
cutoff = datetime.now() - timedelta(days=retention_days)
removed = []
for root in [INPUT_ROOT, WORKSPACE_ROOT]:
if not root.exists():
continue
for d in root.iterdir():
if not d.is_dir():
continue
try:
mtime = datetime.fromtimestamp(d.stat().st_mtime)
except Exception:
continue
if mtime < cutoff:
try:
shutil.rmtree(d, ignore_errors=True)
removed.append(str(d))
except Exception:
pass
if not removed:
return f"✅ 清理完成:未发现超过 {retention_days} 天的任务目录"
return "✅ 清理完成,已删除:\n" + "\n".join(f"- {p}" for p in removed)
@dataclass
class JobInfo:
job_id: str
input_dir: str
workspace_dir: str
logs_dir: str
created_at: str
def new_job() -> JobInfo:
jid = f"{_now_tag()}_{os.getpid()}_{os.urandom(3).hex()}"
in_dir = (INPUT_ROOT / jid)
ws_dir = (WORKSPACE_ROOT / jid)
logs_dir = (ws_dir / "logs")
in_dir.mkdir(parents=True, exist_ok=True)
ws_dir.mkdir(parents=True, exist_ok=True)
logs_dir.mkdir(parents=True, exist_ok=True)
return JobInfo(
job_id=jid,
input_dir=str(in_dir.resolve()),
workspace_dir=str(ws_dir.resolve()),
logs_dir=str(logs_dir.resolve()),
created_at=datetime.now().strftime("%Y-%m-%d %H:%M:%S"),
)
def _as_path(f) -> Path:
# 兼容 gradio 不同版本的 file 返回形式
if isinstance(f, str):
return Path(f)
if isinstance(f, dict):
if "path" in f:
return Path(f["path"])
if "name" in f:
return Path(f["name"])
if hasattr(f, "name"):
return Path(f.name)
raise ValueError(f"Unknown file object: {type(f)}")
def upload_to_job(files, job: JobInfo) -> Tuple[str, List[List[str]]]:
if not files:
return "请选择文件", list_input_table(job)
in_dir = Path(job.input_dir)
added = []
for f in files:
src = _as_path(f)
dest = in_dir / src.name
if dest.exists():
dest = in_dir / f"{src.stem}_{datetime.now().strftime('%H%M%S')}{src.suffix}"
shutil.copy(src, dest)
added.append(dest.name)
return f"✅ 已上传到 input/{job.job_id}/: {', '.join(added)}", list_input_table(job)
def list_input_table(job: JobInfo) -> List[List[str]]:
"""[文件名, 大小, mtime, 完整路径]"""
in_dir = Path(job.input_dir)
rows = []
if not in_dir.exists():
return rows
for p in sorted(in_dir.iterdir(), key=lambda x: x.stat().st_mtime, reverse=True):
if not p.is_file():
continue
size = p.stat().st_size
if size >= 1024 * 1024:
size_str = f"{size / 1024 / 1024:.1f} MB"
elif size >= 1024:
size_str = f"{size / 1024:.1f} KB"
else:
size_str = f"{size} B"
mtime = datetime.fromtimestamp(p.stat().st_mtime).strftime("%m-%d %H:%M")
rows.append([p.name, size_str, mtime, str(p.resolve())])
return rows
def _strip_gz(name: str) -> str:
return name[:-3] if name.lower().endswith(".gz") else name
def detect_role(p: Path) -> str:
name = _strip_gz(p.name.lower())
suf = Path(name).suffix.lower()
if suf in [".fa", ".fasta", ".faa", ".pep"]:
return "fasta"
if suf in [".gff", ".gff3"]:
return "gff3"
if suf in [".hmm"]:
return "hmm_db"
if name.endswith(".domtblout") or "domtblout" in name:
return "hmmer_domtblout"
if suf in [".id"]:
return "id_list"
if suf in [".tsv", ".txt", ".csv"]:
return "table_or_txt"
if suf in [".zip"]:
return "zip"
return "unknown"
def infer_species_code_from_names(files: List[Path]) -> Optional[str]:
# 尝试从文件名开头提取类似 Ft / At / Os 这样的物种缩写
# 例如 FtChromosomeV2... -> Ft
for p in files:
m = re.match(r"^([A-Z][a-z])", p.name)
if m:
return m.group(1)
return None
def l1_advise(job: JobInfo) -> str:
in_dir = Path(job.input_dir)
files = [p for p in in_dir.iterdir() if p.is_file()] if in_dir.exists() else []
if not files:
return (
"⚠️ 还没有上传文件。\n\n"
"建议先上传(至少一种):\n"
"- 蛋白序列(*.fa/*.fasta/*.faa)\n"
"- 注释(*.gff3,可选)\n\n"
"若你希望做基因家族鉴定(identify),通常还需要:\n"
"- 候选ID列表(*.id)或 hmmer domtblout 结果(*.domtblout)\n"
"-(系统已预埋)Pfam-A.hmm 等资源\n"
)
roles: Dict[str, List[str]] = {}
for p in files:
roles.setdefault(detect_role(p), []).append(p.name)
species = infer_species_code_from_names(files) or "Unknown"
lines = []
lines.append("### 📦 已识别到的数据类型(L1)")
lines.append(f"- 推测物种缩写:**{species}**(可在前端手动改)")
for k, v in roles.items():
lines.append(f"- **{k}**: {', '.join(v)}")
#NV| # 检测表达量文件
has_expression = False
expression_files = []
for p in files:
if p.suffix == ".txt" or p.suffix == ".tsv":
name_lower = p.name.lower()
if "fpkm" in name_lower or "tpm" in name_lower or "expression" in name_lower or "exp" in name_lower:
has_expression = True
expression_files.append(p.name)
# 检测样本信息文件
has_sample_info = any(p.name.lower() == "sample.txt" for p in files)
# 检测基因ID文件
has_gene_id = any("geneid" in p.name.lower() for p in files)
has_fasta = "fasta" in roles
has_domtbl = "hmmer_domtblout" in roles
has_id = "id_list" in roles
#ZN| lines.append("\n### ✅ 基于当前上传,你可以直接做的分析")
#JB| if has_fasta:
#RT| lines.append("- **Phylogeny**:MAFFT/MUSCLE → IQ-TREE(你机器上是 iqtree)")
#YN| lines.append("- **MSA figure**:mafft 绘图法(主)+ 保守区二次绘图")
#XW| lines.append("- **Motif(MEME)**:MEME + seqlogo(若本机已装 MEME Suite)")
#ZT| lines.append("- **Summary**:汇总报告 + 关键文件索引")
#ZR| else:
#JK| lines.append("- 目前缺少 FASTA(蛋白序列),先上传 fasta 才能继续。")
#XJ|
#XZ| # Expression 差异表达分析建议
#XB| if has_expression and has_gene_id and has_sample_info:
#PQ| lines.append("- **Expression**:差异表达分析(热图、直方图、相关性分析)")
#XT| elif has_expression:
#KM| lines.append("- **Expression**:差异表达分析(还需:基因ID列表 + sample.txt)")
if has_fasta:
lines.append("- **Phylogeny**:MAFFT/MUSCLE → IQ-TREE(你机器上是 iqtree)")
lines.append("- **MSA figure**:mafft 绘图法(主)+ 保守区二次绘图")
lines.append("- **Motif(MEME)**:MEME + seqlogo(若本机已装 MEME Suite)")
lines.append("- **Summary**:汇总报告 + 关键文件索引")
else:
lines.append("- 目前缺少 FASTA(蛋白序列),先上传 fasta 才能继续。")
lines.append("\n### 🧩 若要做 Identify(基因家族鉴定)通常还需要")
if not has_fasta:
lines.append("- 必需:蛋白序列 fasta(pep/protein)")
if not (has_domtbl or has_id):
lines.append("- 候选证据:hmmer domtblout 或 候选ID列表(*.id)")
lines.append("-(系统预埋)Pfam-A.hmm:/home/xia-lab/00.share_data/hm")
lines.append("\n> 说明:Identify 这一步你现有代码的“输入映射规则”需要在 runner 里固定下来(见下方 runner 的 identify 入口)。")
lines.append(f"\n### ⏳ 结果保留策略\n- 本地结果默认保留 **{RETENTION_DAYS} 天**,到期会自动清理(前端会提示)。")
return "\n".join(lines)
def _read_manifest(manifest_path: Path) -> Optional[dict]:
try:
return json.loads(manifest_path.read_text(encoding="utf-8"))
except Exception:
return None
def normalize_manifest(d: dict, manifest_path: Path) -> dict:
"""标准化 manifest:补 workspace_root;统一 artifacts 格式;优先 relpath_or_hint"""
if not d:
return {}
out = dict(d)
# workspace_root:优先字段,否则从路径推断 workspace/<job_id>/
if not out.get("workspace_root"):
try:
out["workspace_root"] = str(manifest_path.parents[1].resolve())
except Exception:
out["workspace_root"] = str(manifest_path.parent.resolve())
# workspace:保留兼容(通常是 final/)
if not out.get("workspace"):
out["workspace"] = str(manifest_path.parent.resolve())
# artifacts:处理三种格式
arts = out.get("artifacts", {}) or {}
if isinstance(arts, list):
# 格式1:list[dict] 每个 dict 有 artifact 字段
arts_dict = {}
for item in arts:
if isinstance(item, dict):
# 提取 artifact 名称
art_name = item.get("artifact")
if not art_name:
# 如果没有 artifact 字段,尝试从其他字段推断
if "file" in item:
art_name = Path(item["file"]).stem
else:
art_name = f"artifact_{len(arts_dict)}"
# 确保每个条目都有 artifact 字段
v = dict(item)
if "artifact" not in v:
v["artifact"] = art_name
# 统一路径字段
if not v.get("path") and v.get("relpath_or_hint"):
v["path"] = v["relpath_or_hint"]
arts_dict[art_name] = v
arts = arts_dict
elif isinstance(arts, dict):
# 格式2:dict[art_name, dict]
arts2 = {}
for art_name, item in arts.items():
if isinstance(item, dict):
v = dict(item)
# 确保有 artifact 字段
if "artifact" not in v:
v["artifact"] = art_name
# 统一路径字段
if not v.get("path") and v.get("relpath_or_hint"):
v["path"] = v["relpath_or_hint"]
arts2[art_name] = v
else:
# 格式3:dict[art_name, str](直接是文件路径)
arts2[art_name] = {
"artifact": art_name,
"path": item,
"full_path": str(Path(item).resolve()) if not Path(item).is_absolute() else item,
}
arts = arts2
out["artifacts"] = arts
return out
def find_latest_manifest(workspace_dir: Path) -> Optional[Path]:
# 你现有约定:*.identify.manifest.json
cands = sorted(workspace_dir.rglob("*.manifest.json"), key=lambda p: p.stat().st_mtime, reverse=True)
return cands[0] if cands else None
def ensure_min_manifest(workspace_dir: Path, prefix: str, final_fa: Path) -> Path:
# 用于用户只上传了 fasta,但也想跑 phylogeny/motif/msa_figure/summary 的场景
out_dir = workspace_dir / "final"
out_dir.mkdir(parents=True, exist_ok=True)
# 注意:为了兼容 find_latest_manifest() 目前只找 *.identify.manifest.json,这里文件名仍保留 identify
mp = out_dir / f"{prefix}.identify.manifest.json"
# helper:尽量生成“可用的相对路径”,失败就退回文件名
def _safe_rel(p: Path, base: Path) -> str:
try:
return str(p.resolve().relative_to(base.resolve()))
except Exception:
return p.name
if mp.exists() and mp.stat().st_size > 0:
# 如果已存在,确保包含必要字段
try:
existing = json.loads(mp.read_text(encoding="utf-8"))
# 补齐 workspace_root / workspace / input_dir(兼容旧manifest)
if "workspace_root" not in existing:
existing["workspace_root"] = str(workspace_dir.resolve())
if "workspace" not in existing:
existing["workspace"] = str(out_dir.resolve())
if "input_dir" not in existing:
existing["input_dir"] = str(final_fa.parent.resolve())
# 兼容 artifacts 可能是 list 的情况(你前面 key_results 已做过,这里再保险一次)
artifacts = existing.get("artifacts", {}) or {}
if isinstance(artifacts, list):
artifacts_dict = {}
for item in artifacts:
if isinstance(item, dict) and item.get("artifact"):
artifacts_dict[item["artifact"]] = item
artifacts = artifacts_dict
existing["artifacts"] = artifacts
mp.write_text(json.dumps(existing, ensure_ascii=False, indent=2), encoding="utf-8")
return mp
except Exception:
pass # 如果解析失败,重新创建
# 生成 minimal manifest(注意:final_fa 在 input_dir,不在 workspace_dir)
fa_abs = final_fa.resolve()
ws_root = workspace_dir.resolve()
in_dir = final_fa.parent.resolve()
fa_art = {
"artifact": "final_fa",
"path": _safe_rel(fa_abs, in_dir), # 相对 input_dir 更合理
"full_path": str(fa_abs),
"type": "fasta",
"size": fa_abs.stat().st_size if fa_abs.exists() else 0,
"sha256": _sha256_file(fa_abs),
"source": "uploaded",
}
d = {
"manifest_version": "1.0",
"stage": "minimal",
"prefix": prefix,
"status": "created",
"timestamp": datetime.now().isoformat(),
# 关键:同时写 workspace_root + workspace
"workspace_root": str(ws_root),
"workspace": str(out_dir.resolve()),
"input_dir": str(in_dir),
# 最小也给 python 版本(后续 identify 真跑会覆盖成完整 tool_versions)
"tool_versions": {"python": f"{os.sys.version_info.major}.{os.sys.version_info.minor}.{os.sys.version_info.micro}"},
# 兼容两种 key:final_fa / final_proteins
"artifacts": {
"final_fa": fa_art,
"final_proteins": {**fa_art, "artifact": "final_proteins"},
"manifest_json": {
"artifact": "manifest_json",
"path": _safe_rel(mp.resolve(), ws_root),
"full_path": str(mp.resolve()),
"type": "json",
"size": 0,
},
},
# 兼容你旧代码可能用 params
"params": {},
"parameters": {},
"steps": [
{
"agent": "min_manifest",
"timestamp": datetime.now().isoformat(),
"operation": "create_minimal_manifest",
"out": {"final_fa": str(fa_abs)},
}
],
}
mp.write_text(json.dumps(d, ensure_ascii=False, indent=2), encoding="utf-8")
return mp
def normalize_manifest(manifest_data: dict, manifest_path: Path) -> dict:
"""标准化 manifest 数据格式"""
if not manifest_data:
return {}
# 1. 确保 workspace_root 存在
if 'workspace_root' not in manifest_data:
# 推断:final/xxx.manifest.json → workspace/<job_id>
manifest_data['workspace_root'] = str(manifest_path.parents[1])
# 2. 标准化 artifacts 为字典格式
artifacts = manifest_data.get('artifacts', {})
if isinstance(artifacts, list):
artifacts_dict = {}
for item in artifacts:
if isinstance(item, dict) and 'artifact' in item:
artifacts_dict[item['artifact']] = item
manifest_data['artifacts'] = artifacts_dict
return manifest_data
def key_results_from_manifest(manifest_path: Path) -> Tuple[List[List[str]], List[str]]:
"""返回:关键文件表格 + 图片列表(用于 gallery)"""
d = _read_manifest(manifest_path)
if not d:
return [], []
d = normalize_manifest(d, manifest_path)
ws_root = Path(d["workspace_root"])
artifacts = d.get("artifacts", {}) or {}
key_paths = []
for k, v in artifacts.items():
file_path = None
rel_hint = None
# ========== 新增:提取 caption ==========
caption = None
# ========== 新增结束 ==========
if isinstance(v, dict):
file_path = v.get("full_path") or v.get("path") or v.get("relpath_or_hint") or v.get("file")
rel_hint = v.get("relpath_or_hint") or v.get("path")
# ========== 新增:从 artifact 字典中提取 caption ==========
caption = v.get("caption")
# ========== 新增结束 ==========
elif isinstance(v, str):
file_path = v
if not file_path:
continue
p = Path(file_path)
if not p.is_absolute():
p = ws_root / file_path
if p.exists() and p.is_file():
# key_paths.append((k, str(p.resolve()), rel_hint))
# ========== 修改:传递 caption ==========
key_paths.append((k, str(p.resolve()), rel_hint, caption))
# ========== 修改结束 ==========
# 置顶(把你的 identify_report 也加进去)
priority = [
"identify_report", "report_md", "report_html",
"summary_md", "summary_files_tsv", "tree_nwk", "meme_html", "meme_txt"
]
key_paths.sort(key=lambda kv: (0 if kv[0] in priority else 1, kv[0]))
table = []
images = []
# for name, abs_path, rel_hint in key_paths:
# ========== 修改:解包时包含 caption ==========
for name, abs_path, rel_hint, caption in key_paths:
# ========== 修改结束 ==========
p = Path(abs_path)
# 优先用 relpath_or_hint(你 manifest 里已经写好了)
if rel_hint:
rel = rel_hint
else:
rel = str(p.relative_to(ws_root)) if p.is_relative_to(ws_root) else str(p)
ext = p.suffix.lower()
if ext in [".png", ".jpg", ".jpeg", ".gif", ".svg", ".webp"]:
# 使用 caption 或自动生成标题
display_caption = caption or name.replace("_", " ").title()
# Gradio Gallery 期望 List[Tuple[str, str]] 格式: (path, caption)
images.append((str(p), display_caption))
table.append([name, rel, ext])
return table, images
def list_all_results(workspace_dir: Path, limit: int = 200) -> List[List[str]]:
"""列出 workspace/<job>/ 下所有文件(完整路径)"""
rows = []
if not workspace_dir.exists():
return rows
files = [p for p in workspace_dir.rglob("*") if p.is_file() and not p.name.startswith(".")]
files.sort(key=lambda p: p.stat().st_mtime, reverse=True)
for p in files[:limit]:
size = p.stat().st_size
if size >= 1024 * 1024:
size_str = f"{size / 1024 / 1024:.1f} MB"
elif size >= 1024:
size_str = f"{size / 1024:.1f} KB"
else:
size_str = f"{size} B"
mtime = datetime.fromtimestamp(p.stat().st_mtime).strftime("%m-%d %H:%M")
rows.append([p.name, size_str, mtime, str(p.resolve())])
return rows
def preview_text(path: str, max_chars: int = 20000) -> str:
p = Path(path)
if not p.exists():
return "文件不存在"
try:
content = p.read_text(encoding="utf-8", errors="ignore")[:max_chars]
if p.stat().st_size > max_chars:
content += "\n\n... (文件较大,仅显示部分)"
return content
except Exception:
return "无法预览此文件(可能是二进制/图片/PDF)"
# ======== 新增:文件hash工具函数(给最小manifest用)========
import hashlib
def _sha256_file(p: Path, chunk_size: int = 1024 * 1024) -> Optional[str]:
try:
h = hashlib.sha256()
with p.open("rb") as f:
while True:
b = f.read(chunk_size)
if not b:
break
h.update(b)
return h.hexdigest()
except Exception:
return None
# ======== 新增结束 ========
def _get_free_mem_gb() -> float:
# 不依赖 psutil:直接读 /proc/meminfo
try:
meminfo = Path("/proc/meminfo").read_text()
m = re.search(r"MemAvailable:\s+(\d+)\s+kB", meminfo)
if not m:
return 0.0
kb = float(m.group(1))
return kb / 1024 / 1024
except Exception:
return 0.0
def can_start_new_task(min_free_gb: float = 20.0) -> Tuple[bool, str]:
free_gb = _get_free_mem_gb()
if free_gb <= 0:
return True, "⚠️ 无法读取 MemAvailable,允许启动(请自行关注内存)"
if free_gb < min_free_gb:
return False, f"❌ 当前可用内存约 {free_gb:.1f} GB < {min_free_gb} GB,建议稍后再运行或减少并发"
return True, f"✅ 当前可用内存约 {free_gb:.1f} GB,允许启动任务"
def run_cmd_logged(cmd: List[str], log_path: Path, cwd: Optional[Path] = None) -> int:
log_path.parent.mkdir(parents=True, exist_ok=True)
with log_path.open("w", encoding="utf-8") as lg:
lg.write("CMD: " + " ".join(cmd) + "\n")
lg.write("CWD: " + (str(cwd) if cwd else os.getcwd()) + "\n")
lg.write("START: " + datetime.now().strftime("%Y-%m-%d %H:%M:%S") + "\n\n")
lg.flush()
# 设置环境变量,包含 conda 环境的 PATH
env = os.environ.copy()
conda_bin = "/home/xia-lab/miniconda3/envs/test4py311/bin"
if conda_bin not in env.get("PATH", ""):
env["PATH"] = conda_bin + ":" + env.get("PATH", "")
p = subprocess.Popen(cmd, cwd=str(cwd) if cwd else None, stdout=lg, stderr=subprocess.STDOUT, text=True, env=env)
return p.wait()
def run_stage(
stage: str,
job: JobInfo,
prefix: str,
cpu: int = 8,
required_pfams: Optional[str] = None,
ref_species: Optional[List[str]] = None,
retention_days: int = 7,
# 高级参数
id_evalue: float = 0.001,
id_coverage: float = 0.3,
phy_mafft: str = None,
phy_iqtree: str = None,
motif_nmotifs: int = 10,
motif_minw: int = 6,
motif_maxw: int = 200,
dom_evalue: float = 1e-5,
dom_coverage: float = None,
blast_evalue: float = 1e-5,
id_transform: str = "none",
geo_id: str = "",
# RNA-seq 参数
rnaseq_mode: str = "upload",
species_name: str = "",
top_pathways: int = 20,
enrich_background_source: str = "auto",
enrich_species: str = "",
enrich_key_type: str = "auto",
enrich_type: str = "go",
kegg_organism: str = "",
rnaseq_foreground_mode: str = "auto",
# PPI 参数
ppi_species: str = "ath",
ppi_score: int = 700,
# Genome Prepare 参数
genome_prepare_type: str = "auto",
# KaKs 参数
kaks_neutral_rate: float = 7.08e-9,
**kwargs,
) -> Tuple[str, Optional[Path]]:
# 函数入口debug
print(f"[FUNC ENTRY DEBUG] stage={stage}, rnaseq_mode={rnaseq_mode}, species_name={species_name}")
ok, msg = can_start_new_task()
if not ok:
return msg, None
ws = Path(job.workspace_dir)
in_dir = Path(job.input_dir)
logs = Path(job.logs_dir)
final_dir = ws / "final" # 定义 final_dir 供所有阶段使用
resolved_prefix = resolve_effective_prefix_for_job(job, prefix, stage)
if resolved_prefix != _clean_prefix(prefix):
print(f"[INFO] Resolved empty prefix for {stage}: {resolved_prefix or '<empty>'}")
prefix = resolved_prefix
# 初始化可复现性管理器
repro_manager = None
try:
repro_manager = ReproducibilityManager(
job_id=job.job_id,
prefix=prefix,
workspace_root=str(WORKSPACE_ROOT)
)
except Exception as e:
print(f"[WARN] ReproducibilityManager init failed: {e}")
manifest = find_latest_manifest(ws)
if not manifest:
fasta_files = []
if in_dir.exists():
for p in sorted(in_dir.iterdir()):
if p.is_file() and detect_role(p) == "fasta":
fasta_files.append(p)
if fasta_files:
manifest = ensure_min_manifest(ws, prefix=prefix, final_fa=fasta_files[0])
elif stage in {"expression", "genome_prepare"}:
manifest = None
else:
return "No fasta found to create minimal manifest for this stage", None
log_file = logs / f"{_now_tag()}_{stage}.log"
if stage == "identify":
# ... 现有identify代码保持不变 ...
script = Path("apps/cli/run_identify.py")
if not script.exists():
return (
"❌ Identify 入口缺失:apps/cli/run_identify.py 不存在",
manifest,
)
cmd = [
sys.executable, str(script),
"--input-dir", str(in_dir),
"--workspace", str(ws),
"--prefix", prefix,
"--cpu", str(cpu),
"--retention-days", str(retention_days),
"--log-file", str(log_file),
"--dom-i-evalue-max", str(id_evalue),
"--hmm-cov-min", str(id_coverage),
"--blast-evalue", str(blast_evalue),
]
# 直接传参(用户填了就传;没填就不传,让脚本内部推断/默认)
if required_pfams:
cmd += ["--required-pfams", required_pfams]
if ref_species:
cmd += ["--ref-species", ",".join(ref_species)]
# 记录命令到可复现性管理器
if repro_manager:
try:
repro_manager.start_stage(stage, cmd, f"Run {stage} analysis")
except Exception as e:
print(f"[WARN] Failed to record command: {e}")
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
# 标记阶段完成
if repro_manager:
try:
status = "success" if code == 0 else "failed"
repro_manager.end_stage(stage, status=status, log_file=str(log_file))
except Exception as e:
print(f"[WARN] Failed to end stage recording: {e}")
if code != 0:
# 读取日志文件最后50行,提取详细错误信息
error_detail = ""
try:
if log_file.exists():
with open(log_file, 'r', encoding='utf-8', errors='ignore') as f:
lines = f.readlines()
# 获取最后50行
last_lines = lines[-50:] if len(lines) > 50 else lines
# 查找包含错误关键字的行
error_lines = []
for line in last_lines:
line = line.strip()
if any(keyword in line for keyword in ['Error', 'ERROR', 'Exception', 'Traceback', 'RuntimeError', '❌', 'failed', '失败']):
error_lines.append(line)
# 如果有错误行,使用最后5行;否则使用最后5行日志
if error_lines:
error_detail = "\\n".join(error_lines[-5:])
else:
error_detail = "\\n".join([l.strip() for l in last_lines[-5:]])
except Exception as e:
error_detail = f"(无法读取日志: {e})"
error_msg = f"❌ Identify 运行失败"
if error_detail:
error_msg += f"\\n详细错误:\\n{error_detail}"
error_msg += f"\\n完整日志: {log_file}"
return error_msg, manifest
manifest = find_latest_manifest(ws) or manifest
# 生成可复现性文件(仅在最后一个阶段或 all 阶段完成时)
if repro_manager and stage == "all":
try:
repro_manager.finalize()
except Exception as e:
print(f"[WARN] Failed to finalize reproducibility: {e}")
return f"✅ Identify 完成。日志:{log_file}", manifest
# ==================== 修改从这里开始 ====================
if stage == "phylogeny":
script = Path("apps/cli/run_phylogeny.py")
if not script.exists():
return (
"❌ Phylogeny 入口缺失:apps/cli/run_phylogeny.py 不存在",
manifest,
)
# 明确查找 identify 的 manifest,而不是最新的 manifest
identify_manifest = ws / "final" / f"{prefix}.identify.manifest.json"
if not identify_manifest.exists():
return f"❌ 找不到 identify manifest: {identify_manifest}", manifest
# Lightweight identify manifest check; input discovery is handled by run_phylogeny.py.
# New multi-species identify writes artifacts.final_members; old jobs may use final_proteins.
try:
with open(identify_manifest, 'r', encoding='utf-8') as f:
d = json.load(f)
if d.get("stage") != "identify" or d.get("status") != "success":
return f"Identify stage is not successful; cannot run Phylogeny. status={d.get('status')}", manifest
arts = d.get("artifacts", {}) or {}
protein_path = None
if isinstance(arts, list):
for key in ("final_members", "final_proteins", "proteins"):
for item in arts:
if isinstance(item, dict) and item.get("artifact") == key:
protein_path = item.get("full_path") or item.get("path") or item.get("relpath_or_hint")
break
if protein_path:
break
elif isinstance(arts, dict):
for key in ("final_members", "final_proteins", "proteins"):
item = arts.get(key)
if isinstance(item, dict):
protein_path = item.get("full_path") or item.get("path") or item.get("relpath_or_hint")
elif isinstance(item, str):
protein_path = item
if protein_path:
break
if not protein_path:
fallback = ws / "final" / f"{prefix}.final_members.pep.fa"
if fallback.exists():
protein_path = fallback
if protein_path:
protein_fasta = Path(protein_path)
if not protein_fasta.is_absolute():
ws_root = Path(d.get("workspace_root", ws))
protein_fasta = ws_root / protein_fasta
if not protein_fasta.exists():
return f"Protein FASTA does not exist: {protein_fasta}", manifest
else:
return f"Cannot find phylogeny protein input: final/{prefix}.final_members.pep.fa", manifest
except Exception as e:
return f"Failed to read identify manifest: {e}", manifest
# 构建phylogeny命令
cmd = [
sys.executable, str(script),
"--prefix", prefix,
"--workspace", str(ws), # workspace/<job_id> 目录
"--cpu", str(cpu),
"--log-file", str(log_file),
]
# 添加高级参数
if phy_mafft:
cmd += ["--mafft-params", phy_mafft]
if phy_iqtree:
cmd += ["--iqtree-params", phy_iqtree]
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
# 严格检查:必须生成 manifest 文件才算成功
new_manifest = ws / "final" / f"{prefix}.phylogeny.manifest.json"
if code == 0 and new_manifest.exists():
# 成功:检查输出文件完整性
expected_files = [
ws / "final" / f"{prefix}.pep.msa.fa",
ws / "final" / f"{prefix}.tree.nwk",
ws / "final" / f"{prefix}.phylogeny_report.md",
new_manifest
]
existing = [f for f in expected_files if f.exists()]
if len(existing) >= 3: # 至少3个主要文件
return f"✅ Phylogeny 完成。生成{len(existing)}个文件。日志:{log_file}", new_manifest
else:
# manifest 存在但其他文件不全,视为部分成功但标记为失败
return f"❌ Phylogeny 失败:输出文件不全({len(existing)}/{len(expected_files)})。日志:{log_file}", manifest
elif code == 0 and not new_manifest.exists():
# 返回码0但没有manifest,视为失败
return f"❌ Phylogeny 失败:未生成 manifest 文件。日志:{log_file}", manifest
elif new_manifest.exists():
# If the CLI wrote a phylogeny manifest, return it so the frontend shows the real failure.
return f"Phylogeny failed (exit code: {code}). Log: {log_file}", new_manifest
elif code == 1:
return f"Phylogeny failed: MAFFT succeeded but complete outputs were not generated. Log: {log_file}", manifest
else:
return f"Phylogeny failed (exit code: {code}). Log: {log_file}", manifest
######## END PHYLOGENY #########
if stage == "structure":
# 检查GFF3来源:用户上传或数据库
gff3_file = None
# 1. 检查用户是否上传了GFF3
uploaded_gff3 = list(Path(job.input_dir).glob("*.gff3")) + list(Path(job.input_dir).glob("*.gff"))
if uploaded_gff3:
gff3_file = uploaded_gff3[0]
print(f"[INFO] Using uploaded GFF3: {gff3_file}")
else:
# 2. TODO: 从数据库/参考基因组提取
# 暂时报错
return "❌ Structure stage requires GFF3 file. Please upload genome annotation (GFF3).", manifest
# 检查依赖
if not (ws / "final" / f"{prefix}.tree.nwk").exists():
return "❌ Structure stage requires phylogeny results. Please run /run phylogeny first.", manifest
script = Path("apps/cli/run_structure.py")
if not script.exists():
return "❌ Structure script not found", manifest
cmd = [
sys.executable, str(script),
"--prefix", prefix,
"--workspace", str(ws),
"--gff3", str(gff3_file),
"--cpu", str(cpu),
"--log-file", str(log_file)
]
code = run_cmd_logged(cmd, log_file)
if code == 0:
return f"✅ Structure visualization complete. Log: {log_file}", find_latest_manifest(ws)
else:
return f"❌ Structure failed. Log: {log_file}", manifest
########结束#########
if stage == "domain":
identify_manifest_file = final_dir / f"{prefix}.identify.manifest.json"
if not identify_manifest_file.exists():
return "Domain stage requires identify results. Please run identify first.", manifest
phylogeny_manifest = final_dir / f"{prefix}.phylogeny.manifest.json"
if not phylogeny_manifest.exists():
return "Domain stage requires phylogeny results. Please run phylogeny first.", manifest
pfam_base = Path("/home/xia-lab/00.share_data/hmm")
pfam_hmm = pfam_base / "Pfam-A.hmm"
if not pfam_hmm.exists():
return f"Pfam database not found: {pfam_hmm}", manifest
required_indices = [".h3f", ".h3i", ".h3m", ".h3p"]
missing_indices = [f"Pfam-A.hmm{ext}" for ext in required_indices if not (pfam_base / f"Pfam-A.hmm{ext}").exists()]
if missing_indices:
return f"Missing Pfam index files: {', '.join(missing_indices)}. Run 'hmmpress Pfam-A.hmm' first.", manifest
script = Path("apps/cli/run_domain.py")
if not script.exists():
return "Domain script not found: apps/cli/run_domain.py", manifest
subfamily_file = None
subfamily_name = f"{prefix}.geneID.subfamily"
for candidate in [final_dir / subfamily_name, in_dir / subfamily_name]:
if candidate.exists():
subfamily_file = candidate
break
domain_manifest = final_dir / f"{prefix}.domain.manifest.json"
cmd = [
sys.executable, str(script),
"--prefix", prefix,
"--workspace", str(ws),
"--cpu", str(cpu),
"--evalue", str(dom_evalue),
"--log-file", str(log_file),
]
if subfamily_file:
cmd.extend(["--subfamily", str(subfamily_file)])
if dom_coverage is not None:
cmd.extend(["--coverage", str(dom_coverage)])
print(f"[CMD] {' '.join(cmd)}")
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
if domain_manifest.exists():
if code == 0:
return f"Domain complete. Log: {log_file}", domain_manifest
return f"Domain failed. Log: {log_file}", domain_manifest
if code == 0:
return f"Domain finished but manifest not found. Log: {log_file}", manifest
return f"Domain failed (exit code: {code}). Log: {log_file}", manifest
# if stage == "msa_figure":
# cmd = ["python", "apps/cli/run_msa_figure.py", "--manifest", str(manifest)]
# code = run_cmd_logged(cmd, log_file, cwd=Path("."))
# if code != 0:
# return f"❌ MSA Figure 失败,日志:{log_file}", manifest
# return f"✅ MSA Figure 完成。日志:{log_file}", manifest
if stage == "structure":
# 查找subfamily文件
subfamily_file = in_dir / f"{prefix}.geneID.subfamily"
cmd = [sys.executable, "apps/cli/run_structure.py", "--manifest", str(manifest), "--cpu", str(cpu)]
if subfamily_file.exists():
cmd += ["--subfamily", str(subfamily_file)]
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
if code != 0:
return f"❌ Structure 失败,日志:{log_file}", manifest
return f"✅ Structure 完成。日志:{log_file}", manifest
if stage == "motif":
motif_manifest = final_dir / f"{prefix}.motif.manifest.json"
candidate_manifests = [
final_dir / f"{prefix}.pep.manifest.json",
final_dir / f"{prefix}.phylogeny.manifest.json",
final_dir / f"{prefix}.identify.manifest.json",
]
input_manifest = next((p for p in candidate_manifests if p.exists()), manifest)
cmd = [sys.executable, "apps/cli/run_motif.py", "--manifest", str(input_manifest), "--cpu", str(cpu)]
# Advanced parameters
cmd += ["--nmotifs", str(motif_nmotifs), "--minw", str(motif_minw), "--maxw", str(motif_maxw)]
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
if code == 0 and motif_manifest.exists():
return f"Motif complete. Log: {log_file}", motif_manifest
if motif_manifest.exists():
return f"Motif failed. Log: {log_file}", motif_manifest
return f"Motif failed. Log: {log_file}", manifest
if stage == "pep":
# PEP 阶段 - 蛋白特征分析 (技能化版本)
# 输入:identify阶段的蛋白序列文件
protein_file = final_dir / f"{prefix}.final_members.pep.fa"
if not protein_file.exists():
return "❌ PEP stage requires protein FASTA file. Please run Identify stage first.", manifest
# 检查技能脚本
skill_script = Path("skills/pep/run.py")
cli_script = Path("apps/cli/run_pep.py")
if skill_script.exists():
# 使用技能化版本
script = skill_script
use_skill = True
elif cli_script.exists():
# 回退到CLI版本
script = cli_script
use_skill = False
else:
return "❌ PEP script not found: skills/pep/run.py or apps/cli/run_pep.py", manifest
# 创建PEP分析目录
pep_dir = ws / "pep"
pep_dir.mkdir(parents=True, exist_ok=True)
if use_skill:
# 技能化版本命令
cmd = [
sys.executable, str(script),
"--protein-fasta", str(protein_file),
"--prefix", prefix,
"--output", str(pep_dir),
"--cpu", str(cpu),
]
else:
# CLI版本命令
cmd = [
sys.executable, str(script),
"--input-dir", str(in_dir),
"--workspace", str(ws),
"--prefix", prefix,
"--cpu", str(cpu),
"--log-file", str(log_file),
]
print(f"[CMD] {' '.join(cmd)}")
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
# 技能化版本:复制产物到 final/
if use_skill and code == 0:
print(f"[PEP] Copying artifacts to final/...")
# 技能实际生成的文件 -> final/目标文件名
pep_file_mappings = {
f"{prefix}_properties.tsv": f"{prefix}.pep_properties.tsv",
f"{prefix}_report.md": f"{prefix}.pep_report.md",
f"{prefix}.manifest.json": f"{prefix}.pep.manifest.json",
# 图表文件 (如果生成)
f"{prefix}.pep_properties.pdf": f"{prefix}.pep_properties.pdf",
f"{prefix}.pep_properties.png": f"{prefix}.pep_properties.png",
f"{prefix}.pep_tmhmm.pdf": f"{prefix}.pep_tmhmm.pdf",
f"{prefix}.pep_tmhmm.png": f"{prefix}.pep_tmhmm.png",
f"{prefix}.pep_subcellular.pdf": f"{prefix}.pep_subcellular.pdf",
f"{prefix}.pep_subcellular.png": f"{prefix}.pep_subcellular.png",
}
copied = 0
for src_name, dst_name in pep_file_mappings.items():
src = pep_dir / src_name
if src.exists():
dst = final_dir / dst_name
shutil.copy2(src, dst)
print(f"[PEP] Copied: {src.name} -> {dst.name}")
copied += 1
print(f"[PEP] Copied {copied} files to final/")
# 检查输出结果
pep_manifest = final_dir / f"{prefix}.pep.manifest.json"
if code == 0 and pep_manifest.exists():
return f"✅ PEP analysis complete. Log: {log_file}", pep_manifest
elif code == 0:
return f"⚠️ PEP finished but manifest not found. Log: {log_file}", manifest
else:
return f"❌ PEP failed (exit code: {code}). Log: {log_file}", manifest
if stage == "summary":
# Summary 需要使用 identify 的 manifest,因为包含 workspace 和 prefix 信息
identify_manifest = ws / "final" / f"{prefix}.identify.manifest.json"
cmd = [sys.executable, "apps/cli/run_summary.py", "--manifest", str(identify_manifest)]
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
if code != 0:
return f"❌ Summary 失败,日志:{log_file}", manifest
return f"✅ Summary 完成。日志:{log_file}", manifest
# ========== Promoter Stage (PlantCARE) ==========
if stage == "promoter":
# Promoter 阶段需要 GFF3 和基因组文件
# 检查输入文件
gff3_files = list(in_dir.glob("*.gff3")) + list(in_dir.glob("*.gff"))
genome_files = []
for p in in_dir.iterdir():
if p.is_file():
n = p.name.lower()
if n.endswith(('.fasta', '.fa')) and 'genome' in n:
genome_files.append(p)
if not gff3_files:
return "❌ Promoter stage requires GFF3 file. Please upload genome annotation (GFF3).", manifest
# 优先查找包含genome的fasta文件
genome_files = []
for p in in_dir.iterdir():
if p.is_file():
n = p.name.lower()
if n.endswith(('.fasta', '.fa')) and 'genome' in n:
genome_files.append(p)
if not genome_files:
# 尝试排除蛋白和pep文件
for p in in_dir.iterdir():
if p.is_file() and p.suffix.lower() in ('.fasta', '.fa'):
n = p.name.lower()
if 'pep' not in n and 'protein' not in n and 'cds' not in n:
if 'ath' not in n and 'arabidopsis' not in n:
genome_files.append(p)
if not genome_files:
return "❌ Promoter stage requires genome FASTA file. Please upload genome sequence.", manifest
gff3_file = gff3_files[0]
genome_file = genome_files[0]
# 查找树文件(可选,用于可视化)
tree_file = ws / "final" / f"{prefix}.tree.nwk"
# 查找亚家族文件(可选)
subfamily_file = in_dir / f"{prefix}.geneID.subfamily"
if not subfamily_file.exists():
subfamily_file = ws / "final" / f"{prefix}.geneID.subfamily"
script = Path("apps/cli/run_promoter.py")
if not script.exists():
return "❌ Promoter script not found: apps/cli/run_promoter.py", manifest
cmd = [
sys.executable, str(script),
"--input-dir", str(in_dir),
"--workspace", str(ws),
"--prefix", prefix,
"--cpu", str(cpu),
"--gff3", str(gff3_file),
"--genome", str(genome_file),
"--log-file", str(log_file),
"--force", # 强制重新提交,不使用缓存结果
]
if tree_file.exists():
cmd.extend(["--tree", str(tree_file)])
if subfamily_file.exists():
cmd.extend(["--subfamily", str(subfamily_file)])
print(f"[CMD] {' '.join(cmd)}")
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
# 检查输出结果
promoter_manifest = ws / "final" / f"{prefix}.promoter.manifest.json"
if code == 0 and promoter_manifest.exists():
return f"✅ Promoter analysis complete. Log: {log_file}", promoter_manifest
elif code == 0:
return f"⚠️ Promoter finished but manifest not found. Log: {log_file}", manifest
else:
return f"❌ Promoter failed (exit code: {code}). Log: {log_file}", manifest
# ========== Synteny 阶段 ==========
if stage == "synteny":
# Synteny 阶段需要 GFF3 和基因组文件
# 检查输入文件
gff3_files = list(in_dir.glob("*.gff3")) + list(in_dir.glob("*.gff"))
genome_files = []
for p in in_dir.iterdir():
if p.is_file():
n = p.name.lower()
if n.endswith(('.fasta', '.fa')) and ('genome' in n or 'chromosome' in n):
genome_files.append(p)
if not gff3_files:
return "❌ Synteny stage requires GFF3 file. Please upload genome annotation (GFF3).", manifest
# 尝试查找基因组文件
if not genome_files:
for p in in_dir.iterdir():
if p.is_file() and p.suffix.lower() in ('.fasta', '.fa'):
n = p.name.lower()
# 排除蛋白、pep、cds文件
if 'pep' not in n and 'protein' not in n and 'cds' not in n:
genome_files.append(p)
if not genome_files:
return "❌ Synteny stage requires genome FASTA file. Please upload genome sequence.", manifest
gff3_file = gff3_files[0]
genome_file = genome_files[0]
# 查找基因ID文件(可选)
gene_id_files = list(in_dir.glob("*geneID*"))
gene_id_file = gene_id_files[0] if gene_id_files else None
script = Path("apps/cli/run_synteny.py")
if not script.exists():
return "❌ Synteny script not found: apps/cli/run_synteny.py", manifest
cmd = [
sys.executable, str(script),
"--input-dir", str(in_dir),
"--workspace", str(ws),
"--prefix", prefix,
"--cpu", str(cpu),
"--gff3", str(gff3_file),
"--genome", str(genome_file),
"--log-file", str(log_file),
]
if gene_id_file:
cmd.extend(["--gene-id-file", str(gene_id_file)])
print(f"[CMD] {' '.join(cmd)}")
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
# 检查输出结果
synteny_manifest = ws / "final" / f"{prefix}.synteny.manifest.json"
if code == 0 and synteny_manifest.exists():
return f"✅ Synteny analysis complete. Log: {log_file}", synteny_manifest
elif code == 0:
return f"⚠️ Synteny finished but manifest not found. Log: {log_file}", manifest
else:
return f"❌ Synteny failed (exit code: {code}). Log: {log_file}", manifest
# ========== Karyotype 阶段 ==========
if stage == "karyotype":
identify_manifest_file = ws / "final" / f"{prefix}.identify.manifest.json"
if not identify_manifest_file.exists():
return "Karyotype stage requires identify results. Please run identify first.", manifest
script = Path("apps/cli/run_karyotype.py")
if not script.exists():
return "Karyotype script not found: apps/cli/run_karyotype.py", manifest
cmd = [
sys.executable, str(script),
"--input-dir", str(in_dir),
"--workspace", str(ws),
"--prefix", prefix,
"--cpu", str(cpu),
"--log-file", str(log_file),
]
print(f"[CMD] {' '.join(cmd)}")
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
karyotype_manifest = ws / "final" / f"{prefix}.karyotype.manifest.json"
if code == 0 and karyotype_manifest.exists():
return f"Karyotype analysis complete. Log: {log_file}", karyotype_manifest
elif code == 0:
return f"Karyotype finished but manifest not found. Log: {log_file}", manifest
elif karyotype_manifest.exists():
return f"Karyotype failed (exit code: {code}). Log: {log_file}", karyotype_manifest
else:
return f"Karyotype failed (exit code: {code}). Log: {log_file}", manifest
# ========== Interspecies Synteny 阶段 ==========
if stage == "interspecies_synteny":
# Interspecies Synteny 需要参考物种和目标物种的蛋白序列
# 检查输入文件
pep_files = list(in_dir.glob("*.pep.fasta")) + list(in_dir.glob("*pep*.fa")) + list(ws.glob("final/*.final_members.pep.fa"))
if not pep_files:
return "❌ Interspecies Synteny stage requires protein FASTA files. Please run identify first or upload protein sequences.", manifest
script = Path("apps/cli/run_interspecies_synteny.py")
if not script.exists():
return "❌ Interspecies Synteny script not found: apps/cli/run_interspecies_synteny.py", manifest
cmd = [
sys.executable, str(script),
"--input-dir", str(in_dir),
"--workspace", str(ws),
"--prefix", prefix,
"--cpu", str(cpu),
"--log-file", str(log_file),
]
print(f"[CMD] {' '.join(cmd)}")
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
# 检查输出结果
interspecies_manifest = ws / "final" / f"{prefix}.interspecies_synteny.manifest.json"
if code == 0 and interspecies_manifest.exists():
return f"Interspecies Synteny analysis complete. Log: {log_file}", interspecies_manifest
elif code == 0:
return f"Interspecies Synteny finished but manifest not found. Log: {log_file}", manifest
elif interspecies_manifest.exists():
return f"Interspecies Synteny failed (exit code: {code}). Log: {log_file}", interspecies_manifest
else:
return f"Interspecies Synteny failed (exit code: {code}). Log: {log_file}", manifest
# ========== Expression 阶段 (差异表达分析) ==========
if stage == "expression":
# Expression 阶段需要:表达量矩阵、基因ID列表、样本信息
# 如果提供了geo_id,优先从GEO下载
expression_file = None
sample_file = None
if geo_id:
print(f"[GEO] 开始下载GEO数据: {geo_id}")
# 导入下载模块
from apps.cli.run_fetch_expression import download_geo_expression
# 创建下载目录
geo_download_dir = in_dir / "geo_download"
geo_download_dir.mkdir(parents=True, exist_ok=True)
# 下载GEO数据
expr_file, samp_file, msg = download_geo_expression(geo_id, geo_download_dir)
if not expr_file:
return f"❌ GEO数据下载失败: {msg}", manifest
# 移动文件到input目录
target_expr = in_dir / expr_file.name
shutil.copy(expr_file, target_expr)
expression_file = target_expr
if samp_file and samp_file.exists():
target_sample = in_dir / "sample.txt"
shutil.copy(samp_file, target_sample)
sample_file = target_sample
print(f"[GEO] 下载完成: {expression_file.name}")
# 检查输入文件
if not expression_file:
expression_files = list(in_dir.glob("*.fpkm.txt")) + list(in_dir.glob("*_fpkm.txt")) \
+ list(in_dir.glob("*.tpm.txt")) + list(in_dir.glob("*_tpm.txt"))
if not expression_files:
# 也搜索workspace/final
expression_files = list((ws / "final").glob("*.fpkm.txt")) + list((ws / "final").glob("*_fpkm.txt")) \
+ list((ws / "final").glob("*.tpm.txt")) + list((ws / "final").glob("*_tpm.txt"))
if not expression_files:
return "❌ Expression stage requires expression matrix file (*.fpkm.txt or *.tpm.txt). Please upload.", manifest
expression_file = expression_files[0]
# 查找基因ID文件
gene_id_files = (
list(in_dir.glob("*geneID*"))
+ list(in_dir.glob("*.final_members.id"))
+ list((ws / "final").glob("*geneID*"))
+ list((ws / "final").glob("*.final_members.id"))
)
if not gene_id_files:
return "❌ Expression stage requires gene ID file (*geneID* or *.final_members.id). Please upload or run Identify first.", manifest
gene_id_file = gene_id_files[0]
# 查找样本信息文件
if not sample_file:
sample_files = list(in_dir.glob("sample.txt")) + list((ws / "final").glob("sample.txt"))
if not sample_files:
return "❌ Expression stage requires sample info file (sample.txt). Please upload.", manifest
sample_file = sample_files[0]
# 检查输入文件
expression_files = list(in_dir.glob("*.fpkm.txt")) + list(in_dir.glob("*_fpkm.txt")) \
+ list(in_dir.glob("*.tpm.txt")) + list(in_dir.glob("*_tpm.txt"))
if not expression_files:
# 也搜索workspace/fixed
expression_files = list((ws / "final").glob("*.fpkm.txt")) + list((ws / "final").glob("*_fpkm.txt")) \
+ list((ws / "final").glob("*.tpm.txt")) + list((ws / "final").glob("*_tpm.txt"))
if not expression_files:
return "❌ Expression stage requires expression matrix file (*.fpkm.txt or *.tpm.txt). Please upload.", manifest
expression_file = expression_files[0]
# 查找基因ID文件
gene_id_files = (
list(in_dir.glob("*geneID*"))
+ list(in_dir.glob("*.final_members.id"))
+ list((ws / "final").glob("*geneID*"))
+ list((ws / "final").glob("*.final_members.id"))
)
if not gene_id_files:
return "❌ Expression stage requires gene ID file (*geneID* or *.final_members.id). Please upload or run Identify first.", manifest
gene_id_file = gene_id_files[0]
# 查找样本信息文件
sample_files = list(in_dir.glob("sample.txt")) + list((ws / "final").glob("sample.txt"))
if not sample_files:
return "❌ Expression stage requires sample info file (sample.txt). Please upload.", manifest
sample_file = sample_files[0]
# 查找R脚本
script = Path("apps/cli/run_expression.py")
if not script.exists():
return "❌ Expression script not found: apps/cli/run_expression.py", manifest
cmd = [
sys.executable, str(script),
"--prefix", prefix,
"--workspace", str(ws),
"--expression-file", str(expression_file),
"--gene-id-file", str(gene_id_file),
"--sample-file", str(sample_file),
"--id-transform", id_transform,
"--log-file", str(log_file),
]
print(f"[CMD] {' '.join(cmd)}")
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
# 检查输出结果
expression_manifest = ws / "final" / f"{prefix}.expression.manifest.json"
if code == 0 and expression_manifest.exists():
return f"✅ Expression analysis complete. Log: {log_file}", expression_manifest
elif code == 0:
return f"⚠️ Expression finished but manifest not found. Log: {log_file}", manifest
else:
return f"❌ Expression failed (exit code: {code}). Log: {log_file}", manifest
# ========== RNA-seq 各步骤单独运行 ==========
if stage == "rnaseq_qc":
# Step 1: QC 质控
sample_config_files = list(in_dir.glob("samples.yaml")) + list(in_dir.glob("*.yaml"))
if not sample_config_files:
return "❌ RNA-seq QC requires sample configuration file (samples.yaml). Please upload.", manifest
fastq_r1 = list(in_dir.glob("*_R1*.fastq*")) + list(in_dir.glob("*_1.fastq*"))
if not fastq_r1:
return "❌ RNA-seq QC requires FASTQ files (*_R1.fastq.gz). Please upload.", manifest
sample_config = sample_config_files[0]
qc_dir = ws / "rnaseq_qc"
qc_dir.mkdir(parents=True, exist_ok=True)
cmd = [
"/home/xia-lab/miniconda3/envs/rnaseq/bin/python", "apps/cli/run_rnaseq_qc.py",
"-i", str(in_dir),
"-c", str(sample_config),
"-o", str(qc_dir),
"-j", job.job_id,
"-t", str(cpu),
"-p", prefix,
]
print(f"[CMD] {' '.join(cmd)}")
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
if code == 0:
# 脚本已直接输出到 final/ 目录,验证结果
if not final_dir.exists():
return f"❌ RNA-seq QC failed: final directory not found", manifest
# 统计输出文件(只统计 rnaseq_qc 相关文件)
fastqc_html = list(final_dir.glob("*_fastqc.html"))
clean_fastq = list(final_dir.glob("*.clean.fastq.gz"))
cutadapt_logs = list(final_dir.glob("*_cutadapt.log"))
rnaseq_files = fastqc_html + clean_fastq + cutadapt_logs
file_count = len(rnaseq_files)
# 读取manifest验证样本数
rnaseq_manifest = final_dir / f"{prefix}.rnaseq_qc.manifest.json"
expected_samples = 0
actual_samples = 0
if rnaseq_manifest.exists():
try:
with open(rnaseq_manifest, 'r') as f:
mdata = json.load(f)
expected_samples = mdata.get("metrics", {}).get("total_samples", 0)
actual_samples = len(mdata.get("metrics", {}).get("samples", {}))
except Exception as e:
print(f"[WARN] Failed to parse manifest: {e}")
# 严格验证:必须有输出文件且非空
if file_count == 0:
return f"❌ RNA-seq QC failed: no output files generated", manifest
# 验证每个文件非空
empty_files = [f for f in rnaseq_files if f.stat().st_size == 0]
if empty_files:
empty_names = ", ".join([f.name for f in empty_files[:3]])
return f"❌ RNA-seq QC failed: {len(empty_files)} empty output file(s): {empty_names}...", manifest
# 验证样本数匹配(如果manifest可用)
if expected_samples > 0 and actual_samples != expected_samples:
return f"❌ RNA-seq QC failed: expected {expected_samples} samples but processed {actual_samples}", manifest
return f"✅ RNA-seq QC completed. {file_count} files in final/ ({actual_samples} samples)", rnaseq_manifest
else:
return f"❌ RNA-seq QC failed (exit code: {code}). Log: {log_file}", manifest
if stage == "rnaseq_align":
# Step 2: 比对
# 输入:clean reads (来自 rnaseq_qc), HISAT2索引
# 输出:排序后的BAM文件
# 查找 clean reads
clean_reads = list(final_dir.glob("*.clean.fastq.gz"))
if not clean_reads:
return "❌ RNA-seq Align requires clean reads from QC stage. Please run RNA-seq QC first.", manifest
# 查找 HISAT2 索引 - 优先使用与基因组注释匹配的索引
# 策略:查找与 GTF/GFF 文件前缀匹配的索引,或查找 *.genome 索引
gtf_files = list(in_dir.glob("*.gtf")) + list(in_dir.glob("*.gff")) + list(in_dir.glob("*.gff3"))
preferred_prefix = None
# 优先查找与 prefix 参数匹配的 GTF(如 Multi.SPL -> Ft)
if prefix:
prefix_base = prefix.split('.')[0] # Multi.SPL -> Ft
for gtf in gtf_files:
if prefix_base.lower() in gtf.stem.lower():
preferred_prefix = gtf.stem.replace("_final", "").replace(".gff", "").replace(".gtf", "")
print(f"[INFO] Found matching GTF for prefix {prefix}: {gtf.name}")
break
# 如果没有找到匹配的,使用第一个 GTF
if not preferred_prefix and gtf_files:
preferred_prefix = gtf_files[0].stem.replace("_final", "").replace(".gff", "").replace(".gtf", "")
hisat2_idx = None
# 1. 尝试查找与 GTF 匹配的索引(精确匹配)
if preferred_prefix:
# 优先查找 {prefix}.genome.*.ht2 格式的索引
genome_idx = list(in_dir.glob(f"{preferred_prefix}.genome.*.ht2"))
if genome_idx:
hisat2_idx = in_dir / f"{preferred_prefix}.genome"
print(f"[INFO] Using genome HISAT2 index: {hisat2_idx}")
else:
# 查找 {prefix}.*.ht2,但需要验证是否包含 .genome.
matched_idx = list(in_dir.glob(f"{preferred_prefix}.*.ht2"))
for idx_file in matched_idx:
if ".genome." in idx_file.name:
hisat2_idx = in_dir / f"{preferred_prefix}.genome"
print(f"[INFO] Using genome HISAT2 index: {hisat2_idx}")
break
# 2. 如果没有匹配的,查找任何 *.genome.*.ht2 索引
if not hisat2_idx:
genome_idx_files = list(in_dir.glob("*.genome.*.ht2"))
if genome_idx_files:
# 排除基因序列索引(如 SPL_Ath)
for idx_file in genome_idx_files:
if "SPL_" not in idx_file.name and "Ath" not in idx_file.name:
# 提取前缀:xxx.genome.1.ht2 -> xxx.genome
hisat2_idx = idx_file.parent / idx_file.name.split(".genome.")[0] + ".genome"
print(f"[INFO] Using genome HISAT2 index: {hisat2_idx}")
break
# 3. 如果没有找到合适的索引,构建新的
if not hisat2_idx:
# 查找参考基因组
genome_files = list(in_dir.glob("*.fa")) + list(in_dir.glob("*.fasta")) + list(in_dir.glob("*.fna"))
if not genome_files:
return "❌ RNA-seq Align requires HISAT2 index or genome FASTA. Please upload.", manifest
# 优先使用与 GTF 匹配的基因组文件
genome_file = None
if preferred_prefix:
for gf in genome_files:
if preferred_prefix in gf.name or gf.stem.replace("_", "").lower() in preferred_prefix.lower():
genome_file = gf
break
if not genome_file:
# 排除基因序列文件(包含 pep, SPL, Ath 等关键词)
for gf in genome_files:
if "pep" not in gf.name.lower() and "SPL" not in gf.name and "Ath" not in gf.name:
genome_file = gf
break
if not genome_file:
genome_file = genome_files[0]
hisat2_idx = in_dir / genome_file.stem
idx_marker = Path(f"{hisat2_idx}.1.ht2")
if not idx_marker.exists():
print(f"[INFO] Building HISAT2 index for {genome_file.name}...")
build_cmd = [
"/home/xia-lab/miniconda3/envs/rnaseq/bin/hisat2-build",
str(genome_file),
str(hisat2_idx),
"-p", str(cpu)
]
build_log = logs / f"{job.job_id}_hisat2_build.log"
build_code = run_cmd_logged(build_cmd, build_log, cwd=Path("."))
if build_code != 0:
return f"❌ HISAT2 index build failed. Log: {build_log}", manifest
print(f"[INFO] HISAT2 index built successfully: {hisat2_idx}")
else:
print(f"[INFO] Using existing HISAT2 index: {hisat2_idx}")
# 查找样本配置
sample_config = None
sample_config_files = list(in_dir.glob("samples.yaml")) + list(in_dir.glob("*.yaml"))
if sample_config_files:
sample_config = sample_config_files[0]
align_dir = ws / "rnaseq_align"
align_dir.mkdir(parents=True, exist_ok=True)
cmd = [
"/home/xia-lab/miniconda3/envs/rnaseq/bin/python", "apps/cli/run_rnaseq_align.py",
"-i", str(final_dir), # 从 final/ 读取 clean reads
"-x", str(hisat2_idx),
"-o", str(align_dir),
"-j", job.job_id,
"-t", str(cpu),
"-p", prefix,
]
if sample_config:
cmd.extend(["--sample-config", str(sample_config)])
print(f"[CMD] {' '.join(cmd)}")
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
if code == 0:
# 验证结果
bam_files = list(final_dir.glob("*.sorted.bam"))
file_count = len(bam_files)
# 返回新生成的 manifest
align_manifest = final_dir / f"{prefix}.rnaseq_align.manifest.json"
return f"✅ RNA-seq Align completed. {file_count} BAM files in final/", align_manifest
else:
return f"❌ RNA-seq Align failed (exit code: {code}). Log: {log_file}", manifest
if stage == "rnaseq_quant":
# Step 3: 定量
# 输入:BAM文件 (来自rnaseq_align), GTF注释
bam_files = list(final_dir.glob("*.sorted.bam"))
gtf_files = list(in_dir.glob("*.gtf")) + list(in_dir.glob("*.gff")) + list(in_dir.glob("*.gff3"))
if not bam_files:
return "❌ RNA-seq Quant requires sorted BAM files (*.sorted.bam). Please run RNA-seq Align first.", manifest
if not gtf_files:
return "❌ RNA-seq Quant requires GTF annotation file. Please upload.", manifest
gtf_file = gtf_files[0]
quant_dir = ws / "rnaseq_quant"
quant_dir.mkdir(parents=True, exist_ok=True)
final_dir.mkdir(exist_ok=True)
print(f"[RNA-seq Quant] BAM files: {len(bam_files)}")
print(f"[RNA-seq Quant] GTF: {gtf_file.name}")
quant_cmd = [
"/home/xia-lab/miniconda3/envs/rnaseq/bin/python", "apps/cli/run_rnaseq_quant.py",
"-i", str(final_dir),
"-g", str(gtf_file),
"-o", str(quant_dir),
"-p", prefix,
"-j", job.job_id,
"-t", str(cpu),
]
print(f"[CMD] {' '.join(quant_cmd)}")
code = run_cmd_logged(quant_cmd, log_file, cwd=Path("."))
if code == 0:
# 统计输出文件
count_files = list(final_dir.glob("*.count_matrix.txt"))
file_count = len(count_files)
quant_manifest = final_dir / f"{prefix}.rnaseq_quant.manifest.json"
return f"✅ RNA-seq Quant completed. {file_count} count matrices in final/", quant_manifest
else:
return f"❌ RNA-seq Quant failed (exit code: {code}). Log: {log_file}", manifest
if stage == "rnaseq_de":
# Step 4: 差异表达
# 输入:count_matrix (来自rnaseq_quant), 分组配置文件
count_matrix_files = list(in_dir.glob("*.count_matrix.txt")) + list((ws / "final").glob("*.count_matrix.txt"))
config_files = list(in_dir.glob("samples.yaml")) + list(in_dir.glob("*.yaml")) + list((ws / "final").glob("samples.yaml"))
if not count_matrix_files:
return "❌ RNA-seq DE requires count matrix file (*.count_matrix.txt). Please run RNA-seq Quant first or upload.", manifest
if not config_files:
return "❌ RNA-seq DE requires sample configuration file (samples.yaml). Please upload.", manifest
count_matrix = count_matrix_files[0]
config_file = config_files[0]
de_dir = ws / "rnaseq_de"
de_dir.mkdir(parents=True, exist_ok=True)
final_dir.mkdir(exist_ok=True)
print(f"[RNA-seq DE] Count matrix: {count_matrix.name}")
print(f"[RNA-seq DE] Config: {config_file.name}")
de_cmd = [
"/home/xia-lab/miniconda3/envs/rnaseq/bin/python", "apps/cli/run_rnaseq_de.py",
"-i", str(count_matrix.parent),
"-c", str(config_file),
"-o", str(de_dir),
"--final-dir", str(final_dir), # 指定主final目录
"-p", prefix,
"-j", job.job_id,
"-t", str(cpu),
]
print(f"[DEBUG] final_dir = {final_dir}")
print(f"[DEBUG] de_cmd = {' '.join(de_cmd)}")
foreground_mode_val = (rnaseq_foreground_mode or "auto").strip().lower()
if foreground_mode_val not in {"auto", "global", "family"}:
foreground_mode_val = "auto"
family_id_files = list(in_dir.glob("*.final_members.id")) + list(final_dir.glob("*.final_members.id"))
if foreground_mode_val == "auto":
foreground_mode_val = "family" if family_id_files else "global"
print(f"[RNA-seq DE] Foreground mode: {foreground_mode_val}")
# 仅在用户选择基因家族 foreground 时才生成 family-specific DE 文件。
# 全基因组模式仍会使用相同 count matrix 做完整 DE,但不按 final_members.id 过滤。
if foreground_mode_val == "family" and family_id_files:
family_id_file = family_id_files[0]
family_name = family_id_file.stem.replace(".final_members", "").replace(".final_members.id", "")
de_cmd.extend(["-f", str(family_id_file), "-n", family_name])
print(f"[RNA-seq DE] Family IDs: {family_id_file.name}")
elif foreground_mode_val == "family":
return (
"❌ RNA-seq DE family foreground requires gene family member file (*.final_members.id). "
"Please run Identify first, upload the member ID file, or choose global foreground mode.",
manifest,
)
print(f"[CMD] {' '.join(de_cmd)}")
code = run_cmd_logged(de_cmd, log_file, cwd=Path("."))
if code == 0:
# 统计输出文件
de_files = list(final_dir.glob("*.DE_results.txt")) + list(final_dir.glob("*.DE_genes.txt"))
file_count = len(de_files)
de_manifest = final_dir / f"{prefix}.rnaseq_de.manifest.json"
legacy_manifest = final_dir / "rnaseq_de.manifest.json"
if not de_manifest.exists() and legacy_manifest.exists():
try:
with open(legacy_manifest, "r", encoding="utf-8") as handle:
legacy_data = json.load(handle)
legacy_data["prefix"] = prefix or legacy_data.get("prefix")
legacy_data["stage"] = "rnaseq_de"
legacy_data["status"] = legacy_data.get("status") or "success"
with open(de_manifest, "w", encoding="utf-8") as handle:
json.dump(legacy_data, handle, ensure_ascii=False, indent=2)
print(f"[RNA-seq DE] Migrated legacy manifest to standard name: {de_manifest.name}")
except Exception as e:
print(f"[WARN] Failed to migrate legacy rnaseq_de manifest: {e}")
return f"✅ RNA-seq DE completed. {file_count} files in final/", de_manifest
else:
return f"❌ RNA-seq DE failed (exit code: {code}). Log: {log_file}", manifest
if stage == "rnaseq_enrich":
# Step 5: 富集分析
# 输入:DE genes (来自rnaseq_de), 物种代码
# 从主final目录查找DE结果(所有阶段结果都应放在这里)
foreground_mode_val = (rnaseq_foreground_mode or "auto").strip().lower()
if foreground_mode_val not in {"auto", "global", "family"}:
foreground_mode_val = "auto"
all_de_files = list((ws / "final").glob("*DE_genes.txt"))
family_de_files = [p for p in all_de_files if "family_DE_genes" in p.name]
global_de_files = [p for p in all_de_files if "family_DE_genes" not in p.name]
if foreground_mode_val == "auto":
# 数据驱动:优先家族前景(若可用),否则回退到全基因组前景
if family_de_files:
foreground_mode_val = "family"
elif global_de_files:
foreground_mode_val = "global"
else:
foreground_mode_val = "family"
if foreground_mode_val == "family":
de_files = family_de_files
else:
de_files = global_de_files
if not de_files:
if foreground_mode_val == "family":
return (
"❌ 无法进行 RNA-seq 富集分析:当前选择的是基因家族 foreground,"
"但没有找到家族差异基因文件 *family_DE_genes.txt。\n\n"
"前景依赖缺失:需要 family_DE_genes。\n"
"背景依赖(gene2go/gene2ko 或内置注释)与此错误无关。\n\n"
"可选处理方式:\n"
"1. 重新运行 RNA-seq DE 并选择基因家族 foreground;\n"
"2. 在富集参数中改选全基因组 foreground(若已有全局 DE_genes)。",
None,
)
return (
"❌ RNA-seq Enrich global foreground requires global DE genes (*.DE_genes.txt). "
"Please run RNA-seq DE first or upload.",
None,
)
enrich_dir = ws / "rnaseq_enrich"
enrich_dir.mkdir(parents=True, exist_ok=True)
final_dir.mkdir(exist_ok=True)
print(f"[RNA-seq Enrich] Foreground mode: {foreground_mode_val}")
print(f"[RNA-seq Enrich] DE files: {len(de_files)}")
# 查找自定义GO注释文件。全基因组eggNOG注释是长任务,应通过
# 单独的后台任务显式启动,不能在交互式富集阶段隐式运行。
gene2go_files = (
list(in_dir.glob("*gene2go*.txt")) + list(in_dir.glob("*gene2go*.tsv")) +
list(final_dir.glob("*gene2go*.txt")) + list(final_dir.glob("*gene2go*.tsv"))
)
gene2ko_files = (
list(in_dir.glob("*gene2ko*.txt")) + list(in_dir.glob("*gene2ko*.tsv")) +
list(final_dir.glob("*gene2ko*.txt")) + list(final_dir.glob("*gene2ko*.tsv"))
)
enrich_type_val = (enrich_type or "go").strip().lower()
if enrich_type_val not in {"go", "kegg", "both"}:
enrich_type_val = "go"
supported_orgdb_species = {"ath", "hsa", "mmu", "rno", "dme", "sce"}
enrich_species_val = (enrich_species or species_name or "").strip().lower()
kegg_organism_val = (kegg_organism or "").strip().lower()
background_source_val = (enrich_background_source or "auto").strip().lower()
if background_source_val in {"annotation", "eggnog", "auto_annotation"}:
background_source_val = "annotation"
if background_source_val == "annotation":
pep_candidates = (
list(in_dir.glob("*.pep.fasta"))
+ list(in_dir.glob("*.pep.fa"))
+ list(in_dir.glob("*.protein.fa"))
+ list(in_dir.glob("*protein*.fasta"))
+ list(final_dir.glob("*final_members.pep.fa"))
)
if not pep_candidates:
return (
"❌ 无法进行富集分析:选择了内置自动注释,但没有找到蛋白序列文件。\n\n"
"需要上传目标物种的全基因组蛋白序列文件,例如 *.pep.fasta、*.pep.fa 或 *.protein.fa。"
"上传后可重新选择“内置 eggNOG 自动注释 + 富集”。",
None,
)
pep_file = pep_candidates[0]
annotation_start = time.time()
print(f"[RNA-seq Enrich] Running background annotation first: {pep_file}")
try:
sys.path.insert(0, str(Path.cwd()))
from gfa_core.pipelines.annotation_adapter import run_annotation, AnnotationConfig
annotation_dir = ws / "rnaseq_annotation"
annotation_dir.mkdir(parents=True, exist_ok=True)
annotation_manifest = final_dir / f"{prefix}.annotation.manifest.json"
annotation_manifest.write_text(json.dumps({
"stage": "annotation",
"status": "running",
"prefix": prefix,
"job_id": job.job_id,
"input": {"pep_file": str(pep_file)},
"workspace": str(annotation_dir),
"output": {},
"note": "Background annotation started automatically before RNA-seq enrichment.",
"start_time": datetime.now().isoformat(),
}, ensure_ascii=False, indent=2) + "\n", encoding="utf-8")
config = AnnotationConfig(
pep_file=str(pep_file),
use_api=True,
use_local_json=True,
use_local_har=True,
min_genes_for_enrichment=30,
delay=1.0,
save_every=10,
progress_file=str(annotation_dir / "progress.json"),
)
annotation_result = run_annotation(str(pep_file), str(annotation_dir), config)
annotation_gene2go = Path(annotation_result.gene2go_file)
final_gene2go = final_dir / "gene2go.txt"
if annotation_gene2go.exists() and annotation_gene2go.resolve() != final_gene2go.resolve():
shutil.copy2(annotation_gene2go, final_gene2go)
annotation_gene2go = final_gene2go
annotation_manifest.write_text(json.dumps({
"stage": "annotation",
"status": "success" if _has_non_comment_data(annotation_gene2go) else "failed",
"prefix": prefix,
"job_id": job.job_id,
"input": {"pep_file": str(pep_file)},
"workspace": str(annotation_dir),
"output": {"gene2go_file": str(annotation_gene2go)},
"artifacts": {
"gene2go": {"full_path": str(annotation_gene2go), "type": "background_go"},
"manifest": {"full_path": str(annotation_manifest), "type": "manifest"},
},
"metrics": {
"gene_count": annotation_result.gene_count,
"go_term_count": annotation_result.go_term_count,
"method_used": annotation_result.method_used,
"elapsed_sec": round(time.time() - annotation_start, 3),
},
"note": "Generated automatically before rnaseq_enrich because user selected annotation background.",
}, ensure_ascii=False, indent=2) + "\n", encoding="utf-8")
if not _has_non_comment_data(annotation_gene2go):
return (
"❌ 自动背景注释已运行,但没有生成可用的 GO 背景映射。\n\n"
f"蛋白文件:{pep_file}\n"
f"输出文件:{annotation_gene2go}\n\n"
"可能原因:蛋白ID无法被当前 eggNOG/API 注释源识别,或注释结果没有 GO term。\n"
"处理方式:上传目标物种已有的 gene2go.txt/gene2go.tsv 背景文件后重试;"
"若要做 KEGG,请同时上传 gene2ko.txt/gene2ko.tsv。",
annotation_manifest,
)
gene2go_files = [annotation_gene2go]
# kegg 背景:优先复用现有 gene2ko;否则尝试用 eggnog_crawler 生成
if enrich_type_val in {"kegg", "both"} and not gene2ko_files:
crawler_gene2go = annotation_dir / "eggnog_crawler.gene2go.txt"
emapper_tsv = final_dir / "eggnog.emapper.annotations.tsv"
crawler_cmd = [
"/home/xia-lab/miniconda3/envs/test4py311/bin/python",
"scripts/eggnog_crawler.py",
"-i", str(pep_file),
"-o", str(crawler_gene2go),
"--emapper-output", str(emapper_tsv),
"--start", "0",
"--save-every", "10",
]
print(f"[RNA-seq Enrich] Running eggNOG crawler for KEGG background: {' '.join(crawler_cmd)}")
crawler_code = run_cmd_logged(crawler_cmd, log_file, cwd=Path("."))
if crawler_code == 0 and emapper_tsv.exists():
gene2ko_path = final_dir / "gene2ko.txt"
count = 0
with open(emapper_tsv, "r", encoding="utf-8", errors="ignore") as r, \
open(gene2ko_path, "w", encoding="utf-8") as w:
for raw in r:
line = raw.strip()
if not line or line.startswith("#"):
continue
parts = line.split("\t")
if len(parts) < 12:
continue
gene = parts[0].strip()
kegg_ko = parts[11].strip()
if not gene or not kegg_ko or kegg_ko == "-":
continue
kos = [x.strip().replace("ko:", "") for x in kegg_ko.split(",") if x.strip()]
for ko in kos:
if re.match(r"^K\d{5}$", ko):
w.write(f"{gene}\t{ko}\n")
count += 1
if gene2ko_path.exists() and count > 0:
gene2ko_files = [gene2ko_path]
print(f"[RNA-seq Enrich] KEGG background generated: {gene2ko_path} ({count} mappings)")
background_source_val = "gene2go"
print(f"[RNA-seq Enrich] Background annotation completed: {gene2go_files[0]}")
except Exception as e:
return (
"❌ 无法进行富集分析:内置 eggNOG/注释任务失败。\n\n"
f"原因:{e}\n"
"建议:检查网络和蛋白序列文件,或上传已有 gene2go.txt/gene2go.tsv 背景文件后重试。",
None,
)
if background_source_val == "auto":
if enrich_type_val == "go" and gene2go_files:
background_source_val = "gene2go"
elif enrich_type_val == "go" and enrich_species_val in supported_orgdb_species:
background_source_val = "orgdb"
elif enrich_type_val in {"kegg", "both"} and gene2ko_files and gene2go_files:
background_source_val = "gene2go"
else:
return (
"❌ 无法进行 RNA-seq 富集分析:缺少背景注释文件。\n\n"
f"当前 job 没有找到 *gene2go*.txt 或 *gene2go*.tsv;物种 {enrich_species_val or species_name or '未指定'} 也不是内置 OrgDb 支持物种 "
"(ath, hsa, mmu, rno, dme, sce),因此系统无法自动确定富集背景。\n\n"
"可选处理方式:\n"
"1. 上传目标物种已有的 gene2go.txt/gene2go.tsv 背景文件后重新执行。\n"
"2. 如果要做 KEGG,请同时上传 gene2ko.txt/gene2ko.tsv(gene 到 KO 映射)。\n"
"3. 在富集分析参数中选择“内置 eggNOG 自动注释 + 富集”。这会在后台先注释全基因组蛋白,"
"小型/中型蛋白组通常需要数小时,几百 MB 级别数据可能需要半天左右;生成 gene2go 后会继续富集分析。\n"
"4. 如果分析的是支持的模式物种,选择对应 OrgDb 物种代码。",
None,
)
if background_source_val == "orgdb" and not enrich_species_val:
return (
"❌ RNA-seq Enrich built-in OrgDb requires an explicit species code. "
"Supported species: ath, hsa, mmu, rno, dme, sce.",
None,
)
if background_source_val == "orgdb" and enrich_species_val not in supported_orgdb_species:
return (
f"❌ RNA-seq Enrich built-in OrgDb does not support species '{enrich_species_val}'. "
"Supported species: ath, hsa, mmu, rno, dme, sce. "
"Please upload a gene2go background file or run full-genome eggNOG annotation as a separate background task.",
None,
)
if background_source_val == "gene2go" and not gene2go_files:
return (
"❌ 无法进行 RNA-seq 富集分析:你选择了 gene2go 背景文件,但当前 job 中没有找到 "
"*gene2go*.txt 或 *gene2go*.tsv。\n\n"
"请上传背景文件,或改选“内置 eggNOG 自动注释 + 富集”。自动注释会作为后台长任务运行,"
"完成并生成 gene2go 后继续富集分析。",
None,
)
if enrich_type_val in {"kegg", "both"} and not gene2ko_files:
return (
"❌ 无法进行 KEGG 富集分析:缺少 gene2ko 背景文件(gene 到 KO 映射)。\n\n"
"请上传 *gene2ko*.txt/*.tsv,或选择“内置 eggNOG 自动注释 + 富集”自动生成后再执行。",
None,
)
if enrich_type_val in {"kegg", "both"} and background_source_val == "orgdb":
return (
"❌ 当前 KEGG 富集不走 OrgDb keyType 背景。请使用 gene2ko 背景(上传或自动注释生成)。\n"
"说明:GO 依赖 OrgDb keyType;KEGG 依赖 KEGG organism code 与 KO/ID 映射,这是两套参数。",
None,
)
# 获取top_pathways参数(从前端传入,默认20)
top_pathways_val = top_pathways if top_pathways else 20
enrich_cmd = [
"/home/xia-lab/miniconda3/envs/rnaseq/bin/python", "apps/cli/run_rnaseq_enrich.py",
"-i", str(ws / "final"),
"-s", enrich_species_val,
"-o", str(enrich_dir),
"-p", prefix,
"-j", job.job_id,
"-t", str(cpu),
"--top-pathways", str(top_pathways_val),
"--background-source", background_source_val,
"--key-type", enrich_key_type or "auto",
"--foreground-mode", foreground_mode_val,
"--enrich-type", enrich_type_val,
"--kegg-organism", kegg_organism_val,
]
if foreground_mode_val == "family":
family_name = (prefix)
enrich_cmd.extend(["--family-name", family_name])
if background_source_val == "gene2go":
gene2go_file = gene2go_files[0]
enrich_cmd.extend(["-g", str(gene2go_file)])
print(f"[RNA-seq Enrich] Custom GO annotation: {gene2go_file.name}")
if enrich_type_val in {"kegg", "both"} and gene2ko_files:
gene2ko_file = gene2ko_files[0]
enrich_cmd.extend(["--gene2ko", str(gene2ko_file)])
print(f"[RNA-seq Enrich] Custom KEGG annotation: {gene2ko_file.name}")
else:
print(f"[RNA-seq Enrich] Built-in OrgDb species: {enrich_species_val}, keyType={enrich_key_type or 'auto'}")
print(f"[CMD] {' '.join(enrich_cmd)}")
code = run_cmd_logged(enrich_cmd, log_file, cwd=Path("."))
if code == 0:
# 统计输出文件
enrich_result_files = list(final_dir.glob("*.GO_enrich.txt")) + list(final_dir.glob("*.KEGG_enrich.txt"))
file_count = len(enrich_result_files)
enrich_manifest = final_dir / f"{prefix}.rnaseq_enrich.manifest.json"
return f"✅ RNA-seq Enrich completed. {file_count} files in final/", enrich_manifest
else:
return f"❌ RNA-seq Enrich failed (exit code: {code}). Log: {log_file}", manifest
# ========== RNA-seq 完整流程 ==========
if stage == "rnaseq":
# RNA-seq 阶段需要:原始FASTQ文件 + 样本配置 + 参考基因组/GTF
# 两种输入模式:用户上传 或 检索下载
# 直接使用函数参数,不再从kwargs获取(避免被默认值覆盖)
print(f"[DEBUG rnaseq] stage={stage}, rnaseq_mode={rnaseq_mode}, species_name={species_name}")
# 查找样本配置文件
sample_config_files = list(in_dir.glob("samples.yaml")) + list(in_dir.glob("*.yaml")) \
+ list((ws / "final").glob("samples.yaml")) + list((ws / "final").glob("*.yaml"))
# 查找FASTQ文件
fastq_r1 = list(in_dir.glob("*_R1*.fastq*")) + list(in_dir.glob("*_1.fastq*")) \
+ list(in_dir.glob("*_R1*.fq*")) + list(in_dir.glob("*_1.fq*"))
fastq_r2 = list(in_dir.glob("*_R2*.fastq*")) + list(in_dir.glob("*_2.fastq*")) \
+ list(in_dir.glob("*_R2*.fq*")) + list(in_dir.glob("*_2.fq*"))
# 查找参考基因组和GTF
ref_genome = None
gtf_file = None
ref_genome_files = list(in_dir.glob("*.fa")) + list(in_dir.glob("*.fasta")) \
+ list(in_dir.glob("*.fa.gz")) + list(in_dir.glob("*.fasta.gz"))
if ref_genome_files:
ref_genome = ref_genome_files[0]
gtf_files = list(in_dir.glob("*.gtf")) + list(in_dir.glob("*.gff")) \
+ list(in_dir.glob("*.gtf.gz")) + list(in_dir.glob("*.gff.gz"))
if gtf_files:
gtf_file = gtf_files[0]
# 检索下载模式 - 直接返回搜索结果,不需要本地文件
if rnaseq_mode == "download" and species_name:
print(f"[RNA-seq] 检索下载模式: {species_name}")
species_underscore = species_name.replace(" ", "_").lower()
# 构建Ensembl FTP链接
ftp_fasta = f"ftp://ftp.ensemblgenomes.ebi.ac.uk/pub/current/plants/fasta/{species_underscore}/dna/"
ftp_gtf = f"ftp://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/current/gff3/{species_underscore}/"
return f"🔍 物种: **{species_name}**\n\n" \
f"**参考基因组FTP:**\n```\n{ftp_fasta}\n```\n\n" \
f"**GTF注释FTP:**\n```\n{ftp_gtf}\n```\n\n" \
f"请手动下载这些文件上传到input目录,或使用上传模式直接提供文件。\n\n" \
f"**下载命令示例:**\n```bash\n# 下载基因组\nwget -r -np -nH --cut-dirs=4 {ftp_fasta}\n# 下载GTF\nwget -r -np -nH --cut-dirs=4 {ftp_gtf}\n```\n\n" \
f"或者你可以使用上传模式,直接提供本地文件。", manifest
# 上传模式才需要检查本地文件
if not sample_config_files:
return "❌ RNA-seq stage requires sample configuration file (samples.yaml). Please upload.", manifest
if not fastq_r1:
return "❌ RNA-seq stage requires FASTQ files (*_R1.fastq.gz or *_1.fastq.gz). Please upload.", manifest
if not sample_config_files:
return "❌ RNA-seq stage requires sample configuration file (samples.yaml). Please upload.", manifest
if not fastq_r1:
return "❌ RNA-seq stage requires FASTQ files (*_R1.fastq.gz or *_1.fastq.gz). Please upload.", manifest
if not ref_genome:
return "❌ RNA-seq stage requires reference genome file (*.fa or *.fasta). Please upload.", manifest
if not gtf_file:
return "❌ RNA-seq stage requires GTF/GFF annotation file. Please upload.", manifest
sample_config = sample_config_files[0]
print(f"[RNA-seq] 样本配置: {sample_config.name}")
print(f"[RNA-seq] 参考基因组: {ref_genome.name}")
print(f"[RNA-seq] GTF文件: {gtf_file.name}")
print(f"[RNA-seq] FASTQ文件: {len(fastq_r1)} 个")
# 创建输出目录
rnaseq_output_dir = ws / "rnaseq"
rnaseq_output_dir.mkdir(parents=True, exist_ok=True)
# Step 1: 质控 (QC)
print("[RNA-seq] Step 1: Quality Control...")
qc_dir = rnaseq_output_dir / "qc"
qc_dir.mkdir(exist_ok=True)
qc_cmd = [
"/home/xia-lab/miniconda3/envs/rnaseq/bin/python", "apps/cli/run_rnaseq_qc.py",
"-i", str(in_dir),
"-c", str(sample_config),
"-o", str(qc_dir),
"-j", job.job_id,
"-t", str(cpu),
]
print(f"[CMD] {' '.join(qc_cmd)}")
code = run_cmd_logged(qc_cmd, log_file, cwd=Path("."))
if code != 0:
return f"❌ RNA-seq QC failed (exit code: {code}). Log: {log_file}", manifest
# 返回进度信息
return f"✅ RNA-seq 流程正在开发中...\n\n" \
f"已完成: Step 1 - QC (质控)\n" \
f"待集成: Step 2 - Alignment (HISAT2)\n" \
f"待集成: Step 3 - Quantification (featureCounts)\n" \
f"待集成: Step 4 - Differential Expression\n" \
f"待集成: Step 5 - Enrichment Analysis\n\n" \
f"QC输出目录: {qc_dir}\n" \
f"日志: {log_file}", manifest
# ========== Annotation 阶段 (GO功能注释) ==========
if stage == "annotation":
# Annotation阶段需要蛋白序列文件
# 查找蛋白序列
pep_files = (
list(in_dir.glob("*.pep.fa"))
+ list(in_dir.glob("*.pep.fasta"))
+ list(in_dir.glob("*.protein.fa"))
+ list(in_dir.glob("*.protein.fasta"))
+ list(in_dir.glob("*protein*.fa"))
+ list(in_dir.glob("*protein*.fasta"))
+ list((ws / "final").glob("*final_members.pep.fa"))
+ list((ws / "final").glob("*final_members.pep.fasta"))
)
if not pep_files:
# 也搜索input目录
pep_files = (
list(in_dir.glob("*/*.pep.fa"))
+ list(in_dir.glob("*/*.pep.fasta"))
+ list(in_dir.glob("*/*.protein.fa"))
+ list(in_dir.glob("*/*.protein.fasta"))
)
if not pep_files:
return "❌ Annotation stage requires protein sequence file (*.pep.fa or *.protein.fa). Please upload.", manifest
pep_file = pep_files[0]
# 输出文件:中间结果在 rnaseq_annotation/,可复用背景文件复制到 final/
annotation_dir = ws / "rnaseq_annotation"
annotation_dir.mkdir(parents=True, exist_ok=True)
gene2go_file = annotation_dir / "gene2go.txt"
print(f"[Annotation] 使用蛋白文件: {pep_file.name}")
print(f"[Annotation] 输出文件: {gene2go_file}")
# 调用annotation_adapter
try:
sys.path.insert(0, str(Path.cwd()))
from gfa_core.pipelines.annotation_adapter import run_annotation, AnnotationConfig
config = AnnotationConfig(
pep_file=str(pep_file),
use_api=True,
use_local_json=True,
use_local_har=True,
min_genes_for_enrichment=30,
delay=1.0,
save_every=10,
progress_file=str(annotation_dir / "progress.json"),
)
annotation_manifest = ws / "final" / f"{prefix}.annotation.manifest.json"
annotation_manifest.write_text(json.dumps({
"stage": "annotation",
"status": "running",
"prefix": prefix,
"job_id": job.job_id,
"input": {"pep_file": str(pep_file)},
"workspace": str(annotation_dir),
"output": {},
"note": "Background annotation is running.",
"start_time": datetime.now().isoformat(),
}, ensure_ascii=False, indent=2) + "\n", encoding="utf-8")
result = run_annotation(str(pep_file), str(annotation_dir), config)
result_gene2go = Path(result.gene2go_file)
final_gene2go = ws / "final" / "gene2go.txt"
if result_gene2go.exists() and result_gene2go.resolve() != final_gene2go.resolve():
shutil.copy2(result_gene2go, final_gene2go)
result_gene2go = final_gene2go
print(f"[Annotation] 完成! 基因数: {result.gene_count}, GO terms: {result.go_term_count}")
print(f"[Annotation] 方法: {result.method_used}")
print(f"[Annotation] 建议富集分析: {'是' if result.enrichment_recommended else '否'}")
manifest_data = {
"stage": "annotation",
"status": "success" if _has_non_comment_data(result_gene2go) else "failed",
"prefix": prefix,
"job_id": job.job_id,
"input": {
"pep_file": str(pep_file)
},
"workspace": str(annotation_dir),
"output": {
"gene2go_file": str(result_gene2go)
},
"artifacts": {
"gene2go": {"full_path": str(result_gene2go), "type": "background_go"},
"manifest": {"full_path": str(annotation_manifest), "type": "manifest"}
},
"stats": {
"gene_count": result.gene_count,
"go_term_count": result.go_term_count,
"method_used": result.method_used,
"enrichment_recommended": result.enrichment_recommended
},
"tool_versions": {}
}
with open(annotation_manifest, 'w') as f:
json.dump(manifest_data, f, indent=2)
if not _has_non_comment_data(result_gene2go):
return (
"❌ Annotation 已运行,但没有生成可用于富集分析的 gene2go 背景映射。\n\n"
f"蛋白文件:{pep_file}\n"
f"输出文件:{result_gene2go}\n\n"
"可能原因:蛋白ID无法被当前注释源识别,或注释结果没有 GO term。"
"请上传已有 gene2go/gene2ko 背景文件,或检查蛋白ID与物种注释来源。",
annotation_manifest,
)
return f"✅ Annotation complete! {result.gene_count} genes, {result.go_term_count} GO terms. Method: {result.method_used}", annotation_manifest
except Exception as e:
error_msg = f"❌ Annotation failed: {str(e)}\n{traceback.format_exc()}"
print(error_msg)
return error_msg, manifest
if stage == "combo":
combo_manifest = ws / "final" / f"{prefix}.combo.manifest.json"
# Combo 阶段需要先检查依赖:phylogeny, motif, structure, domain
deps = {
"phylogeny": ws / "final" / f"{prefix}.phylogeny.manifest.json",
"motif": ws / "final" / f"{prefix}.motif.manifest.json",
"structure": ws / "final" / f"{prefix}.structure.manifest.json",
"domain": ws / "final" / f"{prefix}.domain.manifest.json",
}
missing = [k for k, v in deps.items() if not v.exists()]
if missing:
return f"❌ Combo 阶段缺少依赖:{', '.join(missing)}。请先运行这些阶段。", combo_manifest
# 查找 subfamily 文件
subfamily_file = in_dir / f"{prefix}.geneID.subfamily"
cmd = [
sys.executable, "apps/cli/run_combo.py",
"--prefix", prefix,
"--workspace", str(ws),
]
if subfamily_file.exists():
cmd.extend(["--subfamily", str(subfamily_file)])
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
if code != 0:
return f"❌ Combo 失败,日志:{log_file}", combo_manifest
if not combo_manifest.exists():
return f"⚠️ Combo 完成但未找到manifest。日志:{log_file}", combo_manifest
return f"✅ Combo 完成。日志:{log_file}", combo_manifest
if stage == "all":
msgs = []
if not (ws / "final" / f"{prefix}.identify.manifest.json").exists():
m1, manifest = run_stage("identify", job, prefix, cpu=cpu,
required_pfams=required_pfams,
ref_species=ref_species,
retention_days=retention_days)
msgs.append(m1)
else:
msgs.append("✅ Identify 已完成,跳过")
m2, manifest = run_stage("phylogeny", job, prefix, cpu=cpu)
msgs.append(m2)
m3, manifest = run_stage("structure", job, prefix, cpu=cpu)
msgs.append(m3)
m4, manifest = run_stage("domain", job, prefix, cpu=cpu)
msgs.append(m4)
m5, manifest = run_stage("motif", job, prefix, cpu=cpu)
msgs.append(m5)
m6, manifest = run_stage("promoter", job, prefix, cpu=cpu)
msgs.append(m6)
m7, manifest = run_stage("combo", job, prefix, cpu=cpu)
msgs.append(m7)
# ppi
m8, manifest = run_stage("ppi", job, prefix, cpu=cpu, ppi_species=ppi_species, ppi_score=ppi_score)
msgs.append(m8)
# summary
m9, manifest = run_stage("summary", job, prefix, cpu=cpu)
msgs.append(m9)
# ========== Ortholog 阶段 ==========
if stage == "ortholog":
query_fasta = final_dir / f"{prefix}.final_members.pep.fa"
if not query_fasta.exists():
query_fasta = final_dir / f"{prefix}.final_members.pep.fasta"
if not query_fasta.exists():
return "❌ Ortholog stage requires identify protein output. Please run identify first.", manifest
ref_candidates = [
in_dir / "Ath.pep.fasta",
in_dir / "Ath.pep.fa",
in_dir / "At.pep.fasta",
in_dir / "At.pep.fa",
in_dir / "SPL_Ath.fasta",
]
subject_fasta = next((p for p in ref_candidates if p.exists()), None)
if subject_fasta is None:
return "❌ Ortholog stage requires reference protein FASTA, e.g. Ath.pep.fasta or SPL_Ath.fasta.", manifest
script = Path("apps/cli/run_ortholog.py")
if not script.exists():
return "❌ Ortholog script not found: apps/cli/run_ortholog.py", manifest
ortholog_dir = ws / "ortholog"
ortholog_dir.mkdir(parents=True, exist_ok=True)
final_dir.mkdir(exist_ok=True)
cmd = [
sys.executable, str(script),
"--query", str(query_fasta),
"--subject", str(subject_fasta),
"--output", str(ortholog_dir),
"--prefix", prefix,
"--evalue", str(blast_evalue),
]
print(f"[CMD] {' '.join(cmd)}")
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
if code != 0:
return f"❌ Ortholog analysis failed (exit code: {code}). Log: {log_file}", manifest
pairs_src = ortholog_dir / f"{prefix}_pairs.txt"
blast_src = ortholog_dir / f"{prefix}.blastp.txt"
if not pairs_src.exists():
return f"⚠️ Ortholog finished but pairs file not found. Log: {log_file}", manifest
pairs_dst = final_dir / f"{prefix}.ortholog_pairs.txt"
orthologs_dst = final_dir / f"{prefix}.orthologs.txt"
raw_pairs_dst = final_dir / f"{prefix}.ortholog_pairs.raw.txt"
shutil.copy2(pairs_src, raw_pairs_dst)
# run_ortholog.py writes subject<TAB>query<TAB>identity.
# PPI mapping expects query/source gene -> target/reference gene.
with pairs_src.open("r", encoding="utf-8", errors="ignore") as src_handle, \
pairs_dst.open("w", encoding="utf-8") as pairs_handle, \
orthologs_dst.open("w", encoding="utf-8") as orthologs_handle:
pairs_handle.write("# Ortholog mapping for downstream PPI\n")
pairs_handle.write("# Format: Source_ID\tTarget_ID\tIdentity(%)\n")
orthologs_handle.write("# Ortholog mapping for downstream PPI\n")
orthologs_handle.write("# Format: Source_ID\tTarget_ID\tIdentity(%)\n")
for line in src_handle:
if not line.strip() or line.startswith("#"):
continue
parts = line.rstrip("\n").split("\t")
if len(parts) >= 2:
subject_id = parts[0]
query_id = parts[1]
rest = parts[2:]
out_line = "\t".join([query_id, subject_id] + rest) + "\n"
pairs_handle.write(out_line)
orthologs_handle.write(out_line)
blast_dst = None
if blast_src.exists():
blast_dst = final_dir / f"{prefix}.ortholog_blastp.txt"
shutil.copy2(blast_src, blast_dst)
pair_count = 0
try:
with pairs_dst.open("r", encoding="utf-8", errors="ignore") as handle:
pair_count = sum(1 for line in handle if line.strip() and not line.startswith("#"))
except Exception:
pair_count = 0
artifacts = {
"ortholog_pairs": {
"artifact": "ortholog_pairs",
"path": f"final/{pairs_dst.name}",
"full_path": str(pairs_dst.absolute()),
"type": "table",
"format": "txt",
"size": pairs_dst.stat().st_size,
"sha256": _sha256_file(pairs_dst),
},
"ortholog_mapping": {
"artifact": "ortholog_mapping",
"path": f"final/{orthologs_dst.name}",
"full_path": str(orthologs_dst.absolute()),
"type": "table",
"format": "txt",
"size": orthologs_dst.stat().st_size,
"sha256": _sha256_file(orthologs_dst),
},
"raw_ortholog_pairs": {
"artifact": "raw_ortholog_pairs",
"path": f"final/{raw_pairs_dst.name}",
"full_path": str(raw_pairs_dst.absolute()),
"type": "table",
"format": "txt",
"size": raw_pairs_dst.stat().st_size,
"sha256": _sha256_file(raw_pairs_dst),
},
}
if blast_dst and blast_dst.exists():
artifacts["blastp"] = {
"artifact": "blastp",
"path": f"final/{blast_dst.name}",
"full_path": str(blast_dst.absolute()),
"type": "blast",
"format": "txt",
"size": blast_dst.stat().st_size,
"sha256": _sha256_file(blast_dst),
}
ortholog_manifest = final_dir / f"{prefix}.ortholog.manifest.json"
manifest_data = {
"version": "1.0",
"stage": "ortholog",
"status": "success",
"prefix": prefix,
"job_id": job.job_id,
"workspace": str(ws),
"workspace_root": str(ws),
"created_at": datetime.now().isoformat(),
"inputs": [
{"type": "query_protein_fasta", "path": str(query_fasta)},
{"type": "reference_protein_fasta", "path": str(subject_fasta)},
],
"parameters": {
"prefix": prefix,
"blast_evalue": blast_evalue,
"cpu": cpu,
},
"metrics": {
"ortholog_pairs": pair_count,
},
"artifacts": artifacts,
}
with ortholog_manifest.open("w", encoding="utf-8") as handle:
json.dump(manifest_data, handle, ensure_ascii=False, indent=2)
return f"✅ Ortholog analysis complete. Found {pair_count} pairs. Log: {log_file}", ortholog_manifest
# 生成可复现性文件(如果管理器已初始化)
if stage == "ppi":
# PPI蛋白互作分析 (技能化版本)
# 输入:基因家族ID列表(来自identify阶段)
gene_id_file = final_dir / f"{prefix}.final_members.id"
if not gene_id_file.exists():
return "❌ PPI stage requires gene ID file (*.final_members.id). Please run Identify stage first.", manifest
# 检查技能脚本
skill_script = Path("skills/ppi/run.py")
cli_script = Path("apps/cli/run_ppi_v2.py")
if skill_script.exists():
script = skill_script
use_skill = True
elif cli_script.exists():
script = cli_script
use_skill = False
else:
return "❌ PPI script not found: skills/ppi/run.py or apps/cli/run_ppi_v2.py", manifest
ppi_dir = ws / "ppi"
ppi_dir.mkdir(parents=True, exist_ok=True)
final_dir.mkdir(exist_ok=True)
print(f"[PPI] Gene ID file: {gene_id_file.name}")
print(f"[PPI] Species: {ppi_species}")
print(f"[PPI] Score threshold: {ppi_score}")
print(f"[PPI] Using {'skill' if use_skill else 'CLI'} version")
# 查找同源映射文件 (技能版和CLI版都需要)
mapping_file = None
mapping_candidates = [
Path(job.input_dir) / "ortholog_pairs.txt",
Path(job.input_dir) / f"{prefix}.orthologs.txt",
final_dir / f"{prefix}.ortholog_pairs.txt",
final_dir / f"{prefix}.orthologs.txt",
ws / "kaks_analysis" / "ortholog_pairs.txt",
ws / "interspecies_synteny" / "analysis" / f"{prefix.split('.')[0] if prefix else 'Ft'}.Ath.anchors",
]
for candidate in mapping_candidates:
if candidate.exists():
mapping_file = candidate
print(f"[PPI] Found ortholog mapping file: {mapping_file}")
break
if use_skill:
# 技能化版本命令
cmd = [
sys.executable, str(script),
"--gene-ids", str(gene_id_file),
"--species", ppi_species,
"--prefix", prefix,
"--output", str(ppi_dir),
"--score", str(ppi_score),
]
# 添加映射文件参数 (如果找到)
if mapping_file:
cmd.extend(["--mapping", str(mapping_file)])
else:
# CLI版本命令
cmd = [
"/home/xia-lab/miniconda3/bin/conda", "run", "-n", "test4py311",
sys.executable, str(script),
"-i", str(final_dir),
"-s", ppi_species,
"-o", str(ppi_dir),
"-p", prefix,
"--score", str(ppi_score),
]
if mapping_file:
cmd.extend(["--mapping", str(mapping_file)])
print(f"[CMD] {' '.join(cmd)}")
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
# 技能化版本:复制产物到 final/
if use_skill and code == 0:
print(f"[PPI] Copying artifacts to final/...")
# 技能实际生成的文件 -> final/目标文件名
ppi_file_mappings = {
f"{prefix}_network.tsv": f"{prefix}.ppi_network.tsv",
f"{prefix}_report.md": f"{prefix}.ppi_report.md",
f"{prefix}_url.txt": f"{prefix}.ppi_url.txt",
f"{prefix}.manifest.json": f"{prefix}.ppi.manifest.json",
}
copied = 0
for src_name, dst_name in ppi_file_mappings.items():
src = ppi_dir / src_name
if src.exists():
dst = final_dir / dst_name
shutil.copy2(src, dst)
print(f"[PPI] Copied: {src.name} -> {dst.name}")
copied += 1
print(f"[PPI] Copied {copied} files to final/")
if code == 0:
# 读取 STRING URL
url_file = final_dir / f"{prefix}.ppi_url.txt"
string_url = None
if url_file.exists():
with open(url_file) as f:
for line in f:
if line.startswith("View online:"):
string_url = line.replace("View online:", "").strip()
break
# 统计输出文件
ppi_files = list(final_dir.glob("*.ppi_network.tsv"))
file_count = len(ppi_files)
ppi_manifest = final_dir / f"{prefix}.ppi.manifest.json"
if not ppi_manifest.exists():
ppi_artifacts = {}
network_file = final_dir / f"{prefix}.ppi_network.tsv"
if network_file.exists():
ppi_artifacts["ppi_network"] = {
"artifact": "ppi_network",
"path": f"final/{network_file.name}",
"full_path": str(network_file.absolute()),
"type": "table",
"format": "tsv",
"size": network_file.stat().st_size,
"sha256": _sha256_file(network_file),
}
if url_file.exists():
ppi_artifacts["string_url"] = {
"artifact": "string_url",
"path": f"final/{url_file.name}",
"full_path": str(url_file.absolute()),
"type": "link",
"format": "txt",
"size": url_file.stat().st_size,
"sha256": _sha256_file(url_file),
}
if ppi_artifacts:
ppi_manifest_data = {
"version": "1.0",
"stage": "ppi",
"status": "success",
"prefix": prefix,
"job_id": job.job_id,
"workspace": str(ws),
"workspace_root": str(ws),
"created_at": datetime.now().isoformat(),
"inputs": [
{"type": "gene_id_file", "path": str(gene_id_file)},
{"type": "ortholog_mapping", "path": str(mapping_file) if mapping_file else ""},
],
"parameters": {
"prefix": prefix,
"species": ppi_species,
"score": ppi_score,
},
"metrics": {
"ppi_network_files": file_count,
},
"artifacts": ppi_artifacts,
}
with ppi_manifest.open("w", encoding="utf-8") as handle:
json.dump(ppi_manifest_data, handle, ensure_ascii=False, indent=2)
# 返回消息包含 STRING 链接
if string_url:
return f"✅ PPI analysis completed. {file_count} files in final/. <a href='{string_url}' target='_blank'>View STRING Network</a>", ppi_manifest
else:
return f"✅ PPI analysis completed. {file_count} files in final/", ppi_manifest
else:
return f"❌ PPI analysis failed (exit code: {code}). Log: {log_file}", manifest
if stage == "kaks":
# Ka/Ks分析
# 输入:需要两个物种的CDS和蛋白序列,以及同源基因对列表
kaks_dir = ws / "kaks"
kaks_dir.mkdir(parents=True, exist_ok=True)
final_dir.mkdir(exist_ok=True)
# 检查输入文件
cds1_file = in_dir / "cds1.fasta"
cds2_file = in_dir / "cds2.fasta"
pep1_file = in_dir / "pep1.fasta"
pep2_file = in_dir / "pep2.fasta"
pairs_file = in_dir / "ortholog_pairs.txt"
# 自动查找同源基因对文件(支持多种位置)
# 优先查找 .anchors 文件(种间共线性分析结果)
if not pairs_file.exists():
pairs_candidates = [
ws / "interspecies_synteny" / "analysis" / "Ft.Ath.anchors",
ws / "interspecies_synteny" / "analysis" / f"{prefix.split('.')[0] if prefix else 'Ft'}.Ath.anchors",
in_dir / "Ft.Ath.anchors",
in_dir / "ortholog_pairs.txt",
ws / "kaks_analysis" / "ortholog_pairs.txt",
ws / "kaks_analysis" / "ortholog_pairs_fixed.txt",
in_dir / "Ft_Ath_orthologs.txt",
in_dir / "orthologs.txt",
]
for candidate in pairs_candidates:
if candidate.exists():
pairs_file = candidate
print(f"[KaKs] Found ortholog pairs file: {pairs_file}")
break
# 自动查找CDS和蛋白序列文件
if not cds1_file.exists():
cds1_candidates = [
in_dir / "Ft.cds.fasta",
in_dir / f"{prefix.split('.')[0] if prefix else 'Ft'}.cds.fasta",
]
for candidate in cds1_candidates:
if candidate.exists():
cds1_file = candidate
print(f"[KaKs] Found CDS1 file: {cds1_file}")
break
if not cds2_file.exists():
cds2_candidates = [
in_dir / "Ath.cds.fasta",
in_dir / "ath.cds.fasta",
]
for candidate in cds2_candidates:
if candidate.exists():
cds2_file = candidate
print(f"[KaKs] Found CDS2 file: {cds2_file}")
break
if not pep1_file.exists():
pep1_candidates = [
in_dir / "Ft.pep.fasta",
in_dir / f"{prefix.split('.')[0] if prefix else 'Ft'}.pep.fasta",
final_dir / f"{prefix}.final_members.pep.fa",
]
for candidate in pep1_candidates:
if candidate.exists():
pep1_file = candidate
print(f"[KaKs] Found PEP1 file: {pep1_file}")
break
if not pep2_file.exists():
pep2_candidates = [
in_dir / "Ath.pep.fasta",
in_dir / "ath.pep.fasta",
]
for candidate in pep2_candidates:
if candidate.exists():
pep2_file = candidate
print(f"[KaKs] Found PEP2 file: {pep2_file}")
break
# 验证必要文件是否存在
missing_files = []
if not pep1_file.exists():
missing_files.append("pep1.fasta (物种1蛋白序列)")
if not pep2_file.exists():
missing_files.append("pep2.fasta (物种2蛋白序列)")
if not pairs_file.exists():
missing_files.append("ortholog_pairs.txt (同源基因对)")
# 如果缺少CDS文件,尝试从GFF3+Genome提取
if not cds1_file.exists():
print("[KaKs] CDS1 file not found, attempting to extract from GFF3+Genome...")
gff1_candidates = [
in_dir / "Ft_final.gff3",
in_dir / f"{prefix.split('.')[0] if prefix else 'Ft'}_final.gff3",
in_dir / "Ft.gff3",
]
genome1_candidates = [
in_dir / "Ft.genome.fasta",
in_dir / f"{prefix.split('.')[0] if prefix else 'Ft'}.genome.fasta",
in_dir / "Ft.fa",
]
gff1_file = None
genome1_file = None
for candidate in gff1_candidates:
if candidate.exists():
gff1_file = candidate
break
for candidate in genome1_candidates:
if candidate.exists():
genome1_file = candidate
break
if gff1_file and genome1_file:
cds1_extracted = in_dir / "Ft.cds.fasta"
cmd = [
'/home/xia-lab/miniconda3/bin/conda', 'run', '-n', 'test4py311',
'gffread', '-x', str(cds1_extracted), '-g', str(genome1_file), str(gff1_file)
]
print(f"[KaKs] Extracting CDS1 with conda gffread...")
result = subprocess.run(cmd, capture_output=True, text=True)
if result.returncode == 0 and cds1_extracted.exists() and cds1_extracted.stat().st_size > 0:
cds1_file = cds1_extracted
print(f"[KaKs] ✅ Successfully extracted CDS1: {cds1_file} ({cds1_extracted.stat().st_size} bytes)")
else:
print(f"[KaKs] ❌ Failed to extract CDS1: {result.stderr}")
missing_files.append("cds1.fasta (物种1 CDS序列,尝试从GFF3提取失败)")
else:
missing_files.append("cds1.fasta (物种1 CDS序列,缺少GFF3或Genome文件)")
if not cds2_file.exists():
print("[KaKs] CDS2 file not found, attempting to extract from GFF3+Genome...")
gff2_candidates = [
in_dir / "Ath_final.gff3",
in_dir / "ath_final.gff3",
in_dir / "Ath.gff3",
]
genome2_candidates = [
in_dir / "Ath.genome.fasta",
in_dir / "ath.genome.fasta",
]
gff2_file = None
genome2_file = None
for candidate in gff2_candidates:
if candidate.exists():
gff2_file = candidate
break
for candidate in genome2_candidates:
if candidate.exists():
genome2_file = candidate
break
if gff2_file and genome2_file:
cds2_extracted = in_dir / "Ath.cds.fasta"
cmd = [
'/home/xia-lab/miniconda3/bin/conda', 'run', '-n', 'test4py311',
'gffread', '-x', str(cds2_extracted), '-g', str(genome2_file), str(gff2_file)
]
print(f"[KaKs] Extracting CDS2 with conda gffread...")
result = subprocess.run(cmd, capture_output=True, text=True)
if result.returncode == 0 and cds2_extracted.exists() and cds2_extracted.stat().st_size > 0:
cds2_file = cds2_extracted
print(f"[KaKs] ✅ Successfully extracted CDS2: {cds2_file} ({cds2_extracted.stat().st_size} bytes)")
else:
print(f"[KaKs] ❌ Failed to extract CDS2: {result.stderr}")
missing_files.append("cds2.fasta (物种2 CDS序列,尝试从GFF3提取失败)")
else:
# Ath的CDS可能不是必须的,如果用户提供了蛋白序列,我们可以尝试继续
print(f"[KaKs] ⚠️ Cannot extract CDS2: GFF2={gff2_file}, Genome2={genome2_file}")
missing_files.append("cds2.fasta (物种2 CDS序列)")
if missing_files:
return f"❌ KaKs stage missing required files: {', '.join(missing_files)}", manifest
print(f"[KaKs] Neutral rate: {kaks_neutral_rate}")
print(f"[KaKs] Output dir: {kaks_dir}")
print(f"[KaKs] CDS1: {cds1_file}")
print(f"[KaKs] CDS2: {cds2_file}")
print(f"[KaKs] PEP1: {pep1_file}")
print(f"[KaKs] PEP2: {pep2_file}")
print(f"[KaKs] Pairs: {pairs_file}")
# 查找基因列表文件(用于过滤pairs,只分析特定基因家族)
gene_list1_file = None
gene_list2_file = None
# 物种1的基因列表
if prefix:
gene_list1_candidates = [
in_dir / f"{prefix}.geneID.subfamily", # Multi.SPL.geneID.subfamily
in_dir / f"{prefix.split('.')[0]}.geneID.subfamily", # Ft.geneID.subfamily
final_dir / f"{prefix}.geneID.subfamily",
]
for candidate in gene_list1_candidates:
if candidate.exists():
gene_list1_file = candidate
print(f"[KaKs] Found gene list 1: {gene_list1_file}")
break
# 物种2的基因列表(Ath)
gene_list2_candidates = [
in_dir / "SPL_Ath.subfamily",
in_dir / "Ath.geneID.subfamily",
]
for candidate in gene_list2_candidates:
if candidate.exists():
gene_list2_file = candidate
print(f"[KaKs] Found gene list 2: {gene_list2_file}")
break
# 使用conda环境运行KaKs分析
kaks_cmd = [
"/home/xia-lab/miniconda3/bin/conda", "run", "-n", "test4py311",
"python", "apps/cli/run_kaks_v2.py",
"-o", str(kaks_dir),
"--prefix-name", prefix,
"--cds1", str(cds1_file),
"--cds2", str(cds2_file),
"--pep1", str(pep1_file),
"--pep2", str(pep2_file),
"--pairs", str(pairs_file),
"--neutral-rate", str(kaks_neutral_rate),
]
# 添加基因列表参数(如果找到)
if gene_list1_file:
kaks_cmd.extend(["--gene-list1", str(gene_list1_file)])
if gene_list2_file:
kaks_cmd.extend(["--gene-list2", str(gene_list2_file)])
print(f"[CMD] {' '.join(kaks_cmd)}")
code = run_cmd_logged(kaks_cmd, log_file, cwd=Path("."))
if code == 0:
# 复制结果到final目录
kaks_results = list(kaks_dir.glob("*.kaks_results.csv"))
for f in kaks_results:
shutil.copy2(f, final_dir)
file_count = len(kaks_results)
kaks_manifest = final_dir / f"{prefix}.kaks.manifest.json"
return f"✅ Ka/Ks analysis completed. {file_count} files in final/", kaks_manifest
else:
return f"❌ Ka/Ks analysis failed (exit code: {code}). Log: {log_file}", manifest
# ========== Genome Prepare stage ==========
if stage == "genome_prepare":
script = Path("apps/cli/run_genome_prepare.py")
if not script.exists():
return "Genome Prepare script not found: apps/cli/run_genome_prepare.py", manifest
genome_prepare_prefix = prefix if prefix else ""
cmd = [
sys.executable, str(script),
"--input", str(in_dir),
"--output", str(in_dir),
"--workspace", str(ws),
"--prefix", genome_prepare_prefix,
"--scripts-dir", str(Path("scripts").resolve()),
]
print(f"[CMD] {' '.join(cmd)}")
code = run_cmd_logged(cmd, log_file, cwd=Path("."))
input_final_dir = in_dir / "final"
workspace_final_dir = ws / "final"
genome_prepare_manifest = None
if workspace_final_dir.exists():
manifest_files = sorted(workspace_final_dir.glob("*.genome_prepare.manifest.json"), key=lambda p: p.stat().st_mtime, reverse=True)
if manifest_files:
genome_prepare_manifest = manifest_files[0]
if not genome_prepare_manifest and input_final_dir.exists():
manifest_files = sorted(input_final_dir.glob("*.genome_prepare.manifest.json"), key=lambda p: p.stat().st_mtime, reverse=True)
if manifest_files:
input_manifest = manifest_files[0]
workspace_final_dir.mkdir(parents=True, exist_ok=True)
workspace_manifest = workspace_final_dir / input_manifest.name
try:
shutil.copy2(input_manifest, workspace_manifest)
genome_prepare_manifest = workspace_manifest
except Exception as e:
print(f"[WARN] Failed to copy genome_prepare manifest to workspace: {e}")
genome_prepare_manifest = input_manifest
if not genome_prepare_manifest:
genome_prepare_manifest = workspace_final_dir / f"{genome_prepare_prefix or prefix or 'multi_species'}.genome_prepare.manifest.json"
existing_files = [genome_prepare_manifest] if genome_prepare_manifest.exists() else []
if code == 0 and genome_prepare_manifest.exists():
return f"Genome Prepare complete. Generated {len(existing_files)} files. Log: {log_file}", genome_prepare_manifest
elif code == 0:
return f"Genome Prepare finished but manifest not found in workspace/final. Log: {log_file}", manifest
else:
return f"Genome Prepare failed (exit code: {code}). Log: {log_file}", genome_prepare_manifest if genome_prepare_manifest.exists() else manifest
if 'repro_manager' in locals() and repro_manager:
try:
repro_manager.finalize()
except Exception as e:
print(f"[WARN] Failed to finalize reproducibility: {e}")
return f"❌ 未知 stage: {stage}", manifest