3155 lines
136 KiB
Plaintext
Executable File
3155 lines
136 KiB
Plaintext
Executable File
# 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
|