1322 lines
51 KiB
Python
Executable File
1322 lines
51 KiB
Python
Executable File
|
||
|
||
#!/usr/bin/env python3
|
||
"""
|
||
Domain stage: 蛋白结构域分析 (Pfam HMM 扫描 + 可视化)
|
||
符合 workspace_contract_v1.md 和 manifest_v1.schema.json 契约
|
||
|
||
执行路径:
|
||
1. 单独触发: 独立执行，自动从上游 manifest 查找输入
|
||
2. Run All: 顺序执行，依赖 Identify 和 Phylogeny 阶段产物
|
||
|
||
目录结构适配:
|
||
- 前端上传: input/{job_id}/ (如 input/20260212_143523_2480772_ec5c6d/)
|
||
- 工作目录: workspace/{job_id}/ (如 workspace/20260212_143523_2480772_ec5c6d/)
|
||
"""
|
||
|
||
import argparse
|
||
import json
|
||
import subprocess
|
||
import sys
|
||
import os
|
||
import re
|
||
import uuid
|
||
from datetime import datetime
|
||
from pathlib import Path
|
||
from typing import Dict, List, Optional, Tuple
|
||
import hashlib
|
||
|
||
PROJECT_ROOT = Path(__file__).absolute().parents[2]
|
||
if str(PROJECT_ROOT) not in sys.path:
|
||
sys.path.insert(0, str(PROJECT_ROOT))
|
||
|
||
from gfa_core.gene_family_scope import read_scope, sanitize_token
|
||
from gfa_core.gene_family_stage_context import apply_stage_context, build_stage_context, user_output_prefix_for, write_legacy_manifest_copy
|
||
from gfa_core.mas.effective_members import MemberResolutionError, resolve_effective_members
|
||
|
||
|
||
# ===================== 配置常量 =====================
|
||
|
||
# Pfam 数据库项目内路径（可配置）
|
||
DEFAULT_PFAM_BASE = str(PROJECT_ROOT / "data" / "hmm")
|
||
DEFAULT_PFAM_HMM = "Pfam-A.hmm"
|
||
|
||
# E-value 默认阈值
|
||
DEFAULT_EVALUE = 1e-5
|
||
|
||
# R 脚本路径（相对于项目根目录）
|
||
DEFAULT_RSCRIPT = "scripts/tree_domain.R"
|
||
|
||
CONDA_BIN = "/home/xzq/miniconda3/envs/test4py311/bin"
|
||
if CONDA_BIN not in os.environ.get("PATH", ""):
|
||
os.environ["PATH"] = CONDA_BIN + ":" + os.environ.get("PATH", "")
|
||
|
||
|
||
# ===================== 工具函数 =====================
|
||
|
||
def display_path(path: Path | str) -> str:
|
||
"""Return the user-facing symlink path instead of resolving data-root links."""
|
||
return str(Path(path).absolute())
|
||
|
||
|
||
def current_input_dir(workspace: Path) -> Path:
|
||
return workspace.parent.parent / "input" / workspace.name
|
||
|
||
|
||
def workspace_root_dir(workspace: Path) -> Path:
|
||
return workspace.parent
|
||
|
||
|
||
def inside_current_job(workspace: Path, path: Path) -> bool:
|
||
input_dir = current_input_dir(workspace)
|
||
try:
|
||
resolved = path.resolve()
|
||
return resolved.is_relative_to(workspace.resolve()) or resolved.is_relative_to(input_dir.resolve())
|
||
except OSError:
|
||
return False
|
||
|
||
|
||
def rebase_current_job_path(workspace: Path, path: Path) -> Optional[Path]:
|
||
"""Map copied legacy workspace/input paths back onto the current job roots."""
|
||
parts = path.parts
|
||
for root_name, current_root in (("workspace", workspace), ("input", current_input_dir(workspace))):
|
||
for idx, part in enumerate(parts):
|
||
if part != root_name:
|
||
continue
|
||
if idx + 1 >= len(parts) or parts[idx + 1] != workspace.name:
|
||
continue
|
||
suffix = Path(*parts[idx + 2:]) if idx + 2 < len(parts) else Path()
|
||
candidate = current_root / suffix
|
||
if candidate.exists():
|
||
return candidate
|
||
return None
|
||
|
||
|
||
def resolve_job_path(workspace: Path, raw_path: str | Path, manifest_path: Optional[Path] = None) -> Optional[Path]:
|
||
"""Resolve a manifest path without trusting stale absolute full_path values."""
|
||
raw_text = str(raw_path or "").strip()
|
||
if not raw_text:
|
||
return None
|
||
path = Path(raw_text)
|
||
if path.is_absolute():
|
||
rebased = rebase_current_job_path(workspace, path)
|
||
if rebased:
|
||
return rebased.absolute()
|
||
if inside_current_job(workspace, path):
|
||
return path.absolute()
|
||
return None
|
||
|
||
candidates: List[Path] = [workspace / path]
|
||
if manifest_path:
|
||
candidates.extend([manifest_path.parent / path, manifest_path.parent / path.name])
|
||
if str(path).startswith("final/") and manifest_path:
|
||
candidates.append(manifest_path.parent / path.name)
|
||
for candidate in candidates:
|
||
if candidate.exists():
|
||
return candidate
|
||
return None
|
||
|
||
|
||
def manifest_artifact_path(workspace: Path, manifest_path: Path, artifact: Dict) -> Optional[Path]:
|
||
if not isinstance(artifact, dict):
|
||
return None
|
||
for key in ("path", "relpath_or_hint", "file", "full_path"):
|
||
resolved = resolve_job_path(workspace, artifact.get(key) or "", manifest_path)
|
||
if resolved:
|
||
return resolved
|
||
return None
|
||
|
||
def get_tool_versions() -> Dict[str, str]:
|
||
"""获取工具版本信息"""
|
||
versions = {}
|
||
|
||
# HMMER 版本
|
||
try:
|
||
r = subprocess.run(["hmmscan", "-h"], capture_output=True, text=True, check=False)
|
||
if r.returncode == 0 and r.stdout:
|
||
first_line = r.stdout.strip().split('\\n')[0]
|
||
versions["hmmer"] = first_line
|
||
except:
|
||
versions["hmmer"] = "not_found"
|
||
|
||
# R 版本
|
||
try:
|
||
r = subprocess.run(["R", "--version"], capture_output=True, text=True, check=False)
|
||
if r.returncode == 0 and r.stdout:
|
||
versions["R"] = r.stdout.split("\\n")[0]
|
||
except:
|
||
versions["R"] = "not_found"
|
||
|
||
# Biopython 版本
|
||
try:
|
||
import Bio
|
||
versions["biopython"] = Bio.__version__
|
||
except:
|
||
versions["biopython"] = "not_found"
|
||
|
||
return versions
|
||
|
||
|
||
def prepare_directories(workspace: Path) -> Dict[str, Path]:
|
||
"""创建Domain阶段的标准目录结构"""
|
||
paths = {
|
||
"analysis_dir": workspace / "domain" / "analysis",
|
||
"final_dir": workspace / "final",
|
||
"logs_dir": workspace / "logs",
|
||
}
|
||
|
||
print("创建目录结构:")
|
||
for name, path in paths.items():
|
||
path.mkdir(parents=True, exist_ok=True)
|
||
print(f" {name}: {path}")
|
||
|
||
return paths
|
||
|
||
|
||
def compute_file_hash(filepath: Path) -> str:
|
||
"""计算文件 SHA256 哈希"""
|
||
sha256_hash = hashlib.sha256()
|
||
try:
|
||
with open(filepath, "rb") as f:
|
||
for byte_block in iter(lambda: f.read(4096), b""):
|
||
sha256_hash.update(byte_block)
|
||
return sha256_hash.hexdigest()
|
||
except:
|
||
return ""
|
||
|
||
|
||
def find_latest_manifest(final_dir: Path, stage: str) -> Optional[Path]:
|
||
"""查找指定阶段的最新 manifest 文件"""
|
||
pattern = f"*.{stage}.manifest.json"
|
||
manifests = list(final_dir.glob(pattern))
|
||
if not manifests:
|
||
return None
|
||
return max(manifests, key=lambda p: p.stat().st_mtime)
|
||
|
||
|
||
def read_manifest(manifest_file: Path) -> Dict:
|
||
"""读取 manifest 文件"""
|
||
with open(manifest_file, 'r') as f:
|
||
return json.load(f)
|
||
|
||
|
||
def resolve_scope_tree_file(ws: Path, prefix: str, scope_data: Optional[Dict], cli_tree: Optional[str] = None) -> Path:
|
||
final_dir = ws / "final"
|
||
if cli_tree:
|
||
tree_path = resolve_job_path(ws, cli_tree) or Path(cli_tree)
|
||
if tree_path.exists() and tree_path.stat().st_size > 0:
|
||
return tree_path.absolute()
|
||
raise FileNotFoundError(f"Specified tree file not found: {cli_tree}")
|
||
|
||
manifest_candidates: List[Path] = []
|
||
scope_id = str((scope_data or {}).get("scope_id") or "").strip()
|
||
if scope_id:
|
||
manifest_candidates.append(final_dir / f"{prefix}.{sanitize_token(scope_id)}.phylogeny.manifest.json")
|
||
else:
|
||
manifest_candidates.append(final_dir / f"{prefix}.phylogeny.manifest.json")
|
||
|
||
for manifest_file in manifest_candidates:
|
||
if not manifest_file.exists() or manifest_file.stat().st_size == 0:
|
||
continue
|
||
try:
|
||
manifest_data = read_manifest(manifest_file)
|
||
except Exception:
|
||
continue
|
||
artifact = (manifest_data.get("artifacts", {}) or {}).get("tree")
|
||
tree_path = None
|
||
if isinstance(artifact, dict):
|
||
tree_path = manifest_artifact_path(ws, manifest_file, artifact)
|
||
elif isinstance(artifact, str):
|
||
tree_path = resolve_job_path(ws, artifact, manifest_file)
|
||
if tree_path and tree_path.exists() and tree_path.stat().st_size > 0:
|
||
return tree_path.absolute()
|
||
|
||
fallback = final_dir / f"{prefix}.tree.nwk"
|
||
if not scope_id and fallback.exists() and fallback.stat().st_size > 0:
|
||
return fallback.absolute()
|
||
raise FileNotFoundError(f"Tree file not found for prefix {prefix}. Please run phylogeny first or specify --tree.")
|
||
|
||
|
||
def count_fasta_records(fasta_file: Path) -> int:
|
||
count = 0
|
||
with open(fasta_file, encoding="utf-8", errors="ignore") as handle:
|
||
for line in handle:
|
||
if line.startswith(">"):
|
||
count += 1
|
||
return count
|
||
|
||
|
||
def ensure_protein_fasta_without_stop_markers(protein_fasta: Path) -> Path:
|
||
"""Use a job-local derived FASTA without '*' stop markers when needed."""
|
||
has_stop = False
|
||
with open(protein_fasta, encoding="utf-8", errors="ignore") as handle:
|
||
for line in handle:
|
||
if not line.startswith(">") and "*" in line:
|
||
has_stop = True
|
||
break
|
||
if not has_stop:
|
||
return protein_fasta
|
||
clean_name = re.sub(r"(\.fa(?:sta)?)$", r".clean\1", protein_fasta.name, flags=re.IGNORECASE)
|
||
if clean_name == protein_fasta.name:
|
||
clean_name = protein_fasta.name + ".clean.fa"
|
||
clean_fasta = protein_fasta.with_name(clean_name)
|
||
if clean_fasta.exists() and clean_fasta.stat().st_mtime >= protein_fasta.stat().st_mtime:
|
||
return clean_fasta
|
||
with open(protein_fasta, "r", encoding="utf-8", errors="ignore") as src, open(clean_fasta, "w", encoding="utf-8") as dst:
|
||
for line in src:
|
||
dst.write(line if line.startswith(">") else line.replace("*", ""))
|
||
print(f"[INFO] Created stop-marker-clean FASTA for hmmscan: {clean_fasta}")
|
||
return clean_fasta
|
||
|
||
|
||
def resolve_scope_protein_fasta(ws: Path, prefix: str, scope_data: Dict, cli_proteins: Optional[str] = None) -> Path:
|
||
if cli_proteins:
|
||
protein_path = resolve_job_path(ws, cli_proteins) or Path(cli_proteins)
|
||
if protein_path.exists() and protein_path.stat().st_size > 0:
|
||
return ensure_protein_fasta_without_stop_markers(protein_path)
|
||
raise FileNotFoundError(f"Specified protein file not found: {protein_path}")
|
||
|
||
final_dir = ws / "final"
|
||
scope_id = str(scope_data.get("scope_id") or "")
|
||
manifest_candidates = []
|
||
if scope_id:
|
||
manifest_candidates.append(final_dir / f"{prefix}.{sanitize_token(scope_id)}.phylogeny.manifest.json")
|
||
else:
|
||
manifest_candidates.append(final_dir / f"{prefix}.phylogeny.manifest.json")
|
||
|
||
for manifest_file in manifest_candidates:
|
||
if not manifest_file.exists():
|
||
continue
|
||
try:
|
||
manifest_data = read_manifest(manifest_file)
|
||
except Exception:
|
||
continue
|
||
fasta_text = manifest_data.get("parameters", {}).get("scope_input_fasta")
|
||
if fasta_text:
|
||
fasta_path = resolve_job_path(ws, fasta_text, manifest_file)
|
||
if fasta_path and fasta_path.exists() and fasta_path.stat().st_size > 0:
|
||
print(f"[INFO] Found scoped protein FASTA from phylogeny manifest: {fasta_path}")
|
||
return ensure_protein_fasta_without_stop_markers(fasta_path)
|
||
|
||
raise FileNotFoundError(
|
||
"Cannot find scoped protein FASTA. Run phylogeny with the same scope first, "
|
||
"or provide --proteins explicitly."
|
||
)
|
||
|
||
|
||
def resolve_project_root(ws: Path) -> Path:
|
||
"""
|
||
解析项目根目录
|
||
从 workspace/{job_id}/ 向上回溯到项目根目录
|
||
假设结构: .../bba/workspace/{job_id}/ -> .../bba/
|
||
"""
|
||
# ws 是 workspace/{job_id}/，parent 是 workspace/，parent.parent 才是项目根目录
|
||
return ws.parent.parent
|
||
|
||
|
||
def resolve_subfamily_file(ws: Path, prefix: str, job_id: str, cli_subfamily: Optional[str] = None, output_prefix: Optional[str] = None) -> Optional[Path]:
|
||
"""
|
||
解析亚家族文件路径
|
||
命名规则: {prefix}.confirmed_subfamily.tsv (如 Ft.SPL.confirmed_subfamily.tsv)
|
||
优先级: final/ > input/{job_id}/ (项目根目录)
|
||
|
||
实际目录结构：
|
||
- workspace/{job_id}/final/ <- 阶段产物
|
||
- input/{job_id}/ <- 前端上传（项目根目录下）
|
||
"""
|
||
# 方式一：命令行指定
|
||
if cli_subfamily:
|
||
if str(cli_subfamily).strip().lower() in {"none", "skip", "disabled"}:
|
||
print("[INFO] Subfamily coloring explicitly disabled")
|
||
return None
|
||
subfamily_path = resolve_job_path(ws, cli_subfamily) or Path(cli_subfamily)
|
||
if subfamily_path.exists():
|
||
print(f"[INFO] Using specified subfamily file: {subfamily_path}")
|
||
return subfamily_path
|
||
else:
|
||
print(f"[WARNING] Specified subfamily file not found: {cli_subfamily}")
|
||
return None
|
||
|
||
# 自动查找只接受 confirmed_subfamily；prediction/geneID 必须由 --subfamily 显式传入。
|
||
names = []
|
||
for base in [output_prefix, prefix]:
|
||
if not base or base in names:
|
||
continue
|
||
names.extend([
|
||
f"{base}.confirmed_subfamily.tsv",
|
||
])
|
||
subfamily_name = f"{prefix}.confirmed_subfamily.tsv"
|
||
|
||
# 项目根目录（workspace 的父目录的父目录）
|
||
project_root = resolve_project_root(ws)
|
||
|
||
# 搜索路径优先级：
|
||
search_paths = []
|
||
|
||
# 1. final/ 目录（阶段产物，优先级最高）
|
||
final_dir = ws / "final"
|
||
if final_dir.exists():
|
||
search_paths.extend(final_dir / name for name in names)
|
||
|
||
# 2. 项目根目录的 input/{job_id}/ 目录（前端上传的实际位置）
|
||
# 修正：input/ 直接在项目根目录下，不是 workspace/ 下
|
||
search_paths.extend(
|
||
project_root / "input" / job_id / name
|
||
for name in names
|
||
|
||
)
|
||
|
||
# 3. workspace 下的 input/ 目录（兼容性检查，某些旧版本可能用这里）
|
||
workspace_input = ws / "input"
|
||
if workspace_input.exists():
|
||
search_paths.extend(workspace_input / name for name in names)
|
||
|
||
for candidate in search_paths:
|
||
if candidate.exists():
|
||
print(f"[INFO] Found subfamily file: {candidate}")
|
||
return candidate
|
||
|
||
print(f"[INFO] No subfamily file found ({subfamily_name})")
|
||
print(f"[DEBUG] Searched paths:")
|
||
for p in search_paths:
|
||
print(f" - {p} (exists: {p.exists()})")
|
||
return None
|
||
|
||
def prepare_pfam_database(domain_dir: Path, pfam_base: str = DEFAULT_PFAM_BASE) -> Tuple[Path, bool, str]:
|
||
"""
|
||
准备 Pfam 数据库（软连接方式）
|
||
|
||
Returns:
|
||
(pfam_hmm_path, success, error_message)
|
||
"""
|
||
pfam_source = Path(pfam_base) / DEFAULT_PFAM_HMM
|
||
|
||
# 检查源文件是否存在
|
||
if not pfam_source.exists():
|
||
return None, False, f"Pfam source not found: {pfam_source}"
|
||
|
||
# 检查索引文件
|
||
index_extensions = ['.h3f', '.h3i', '.h3m', '.h3p']
|
||
missing_indices = []
|
||
|
||
for ext in index_extensions:
|
||
idx_file = Path(pfam_base) / f"{DEFAULT_PFAM_HMM}{ext}"
|
||
if not idx_file.exists():
|
||
missing_indices.append(str(idx_file))
|
||
|
||
if missing_indices:
|
||
return None, False, f"Missing Pfam index files: {', '.join(missing_indices)}. Run 'hmmpress {pfam_source}' first."
|
||
|
||
# 创建软连接
|
||
pfam_link = domain_dir / DEFAULT_PFAM_HMM
|
||
|
||
try:
|
||
if pfam_link.exists() or pfam_link.is_symlink():
|
||
pfam_link.unlink()
|
||
os.symlink(pfam_source, pfam_link)
|
||
print(f"[INFO] Linked Pfam database: {pfam_source} -> {pfam_link}")
|
||
|
||
# 链接索引文件
|
||
for ext in index_extensions:
|
||
src_idx = Path(pfam_base) / f"{DEFAULT_PFAM_HMM}{ext}"
|
||
link_idx = domain_dir / f"{DEFAULT_PFAM_HMM}{ext}"
|
||
if link_idx.exists() or link_idx.is_symlink():
|
||
link_idx.unlink()
|
||
os.symlink(src_idx, link_idx)
|
||
|
||
print(f"[INFO] Linked {len(index_extensions)} index files")
|
||
return pfam_link, True, ""
|
||
|
||
except Exception as e:
|
||
return None, False, f"Failed to create Pfam symlinks: {str(e)}"
|
||
|
||
|
||
def run_hmmscan(pfam_hmm: Path, protein_fasta: Path, domain_dir: Path,
|
||
output_prefix: str = "", cpu: int = 8, evalue: float = DEFAULT_EVALUE) -> Tuple[bool, Dict]:
|
||
"""
|
||
运行 hmmscan 进行结构域搜索
|
||
|
||
Returns:
|
||
(success, result_dict)
|
||
"""
|
||
basename = sanitize_token(output_prefix) if output_prefix else "results"
|
||
domtblout = domain_dir / f"{basename}.results.domtblout"
|
||
tblout = domain_dir / f"{basename}.results.tblout"
|
||
output_log = domain_dir / f"{basename}.results.out"
|
||
|
||
cmd = [
|
||
"hmmscan",
|
||
"--cpu", str(cpu),
|
||
"--domtblout", str(domtblout),
|
||
"--tblout", str(tblout),
|
||
"-o", str(output_log),
|
||
"-E", str(evalue), # 报告阈值
|
||
"--domE", str(evalue), # 结构域报告阈值
|
||
str(pfam_hmm),
|
||
str(protein_fasta)
|
||
]
|
||
|
||
print(f"[CMD] {' '.join(cmd)}")
|
||
|
||
start_time = datetime.now()
|
||
|
||
try:
|
||
result = subprocess.run(cmd, capture_output=True, text=True, check=False)
|
||
end_time = datetime.now()
|
||
|
||
success = result.returncode == 0
|
||
elapsed = (end_time - start_time).total_seconds()
|
||
|
||
# 检查输出文件
|
||
domtblout_exists = domtblout.exists() and domtblout.stat().st_size > 0
|
||
|
||
if success and not domtblout_exists:
|
||
success = False
|
||
error_msg = "hmmscan returned 0 but domtblout is empty"
|
||
elif not success:
|
||
error_msg = result.stderr[-2000:] if result.stderr else "hmmscan failed"
|
||
else:
|
||
error_msg = ""
|
||
|
||
return success, {
|
||
"domtblout": domtblout if domtblout.exists() else None,
|
||
"tblout": tblout if tblout.exists() else None,
|
||
"output_log": output_log if output_log.exists() else None,
|
||
"returncode": result.returncode,
|
||
"stderr": result.stderr[-2000:] if result.stderr else "",
|
||
"elapsed_sec": elapsed,
|
||
"error": error_msg
|
||
}
|
||
|
||
except Exception as e:
|
||
return False, {
|
||
"domtblout": None,
|
||
"tblout": None,
|
||
"output_log": None,
|
||
"returncode": -1,
|
||
"stderr": str(e),
|
||
"elapsed_sec": 0,
|
||
"error": str(e)
|
||
}
|
||
|
||
|
||
def parse_domtblout(domtblout: Path, evalue_cutoff: float = DEFAULT_EVALUE) -> Tuple[List[Dict], Dict]:
|
||
"""
|
||
解析 hmmscan domtblout 文件
|
||
|
||
Returns:
|
||
(domains_list, stats_dict)
|
||
"""
|
||
domains = []
|
||
stats = {
|
||
"total_hits": 0,
|
||
"significant_hits": 0,
|
||
"unique_proteins": set(),
|
||
"unique_domains": set(),
|
||
"protein_domain_counts": {}
|
||
}
|
||
|
||
if not domtblout.exists():
|
||
return [], stats
|
||
|
||
with open(domtblout, 'r') as f:
|
||
for line in f:
|
||
line = line.strip()
|
||
if not line or line.startswith('#'):
|
||
continue
|
||
|
||
stats["total_hits"] += 1
|
||
|
||
# domtblout 格式（空格分隔，22列）
|
||
parts = line.split()
|
||
if len(parts) < 22:
|
||
continue
|
||
|
||
try:
|
||
# 解析关键字段
|
||
target_name = parts[0] # Domain name
|
||
t_accession = parts[1] # Domain accession
|
||
tlen = int(parts[2]) # Domain length
|
||
query_name = parts[3] # Protein ID
|
||
q_accession = parts[4] # Protein accession
|
||
qlen = int(parts[5]) # Protein length
|
||
|
||
# E-values and scores
|
||
E_full = float(parts[6]) # Full sequence E-value
|
||
score_full = float(parts[7]) # Full sequence score
|
||
|
||
# Domain coordinates
|
||
dom_index = int(parts[9]) # Domain index
|
||
dom_of = int(parts[10]) # Total domains
|
||
|
||
c_evalue = float(parts[11]) # Conditional E-value
|
||
i_evalue = float(parts[12]) # Independent E-value
|
||
dom_score = float(parts[13]) # Domain score
|
||
|
||
hmm_from = int(parts[15]) # HMM start
|
||
hmm_to = int(parts[16]) # HMM end
|
||
ali_from = int(parts[17]) # Alignment start (protein)
|
||
ali_to = int(parts[18]) # Alignment end (protein)
|
||
env_from = int(parts[19]) # Envelope start
|
||
env_to = int(parts[20]) # Envelope end
|
||
acc = float(parts[21]) # Average posterior probability
|
||
|
||
# 应用 E-value 过滤（使用 i-Evalue，更严格）
|
||
if i_evalue > evalue_cutoff:
|
||
continue
|
||
|
||
stats["significant_hits"] += 1
|
||
stats["unique_proteins"].add(query_name)
|
||
stats["unique_domains"].add(target_name)
|
||
|
||
if query_name not in stats["protein_domain_counts"]:
|
||
stats["protein_domain_counts"][query_name] = 0
|
||
stats["protein_domain_counts"][query_name] += 1
|
||
|
||
domains.append({
|
||
"target_name": target_name,
|
||
"target_accession": t_accession,
|
||
"query_name": query_name,
|
||
"query_length": qlen,
|
||
"i_evalue": i_evalue,
|
||
"domain_score": dom_score,
|
||
"hmm_from": hmm_from,
|
||
"hmm_to": hmm_to,
|
||
"ali_from": ali_from,
|
||
"ali_to": ali_to,
|
||
"env_from": env_from,
|
||
"env_to": env_to,
|
||
"accuracy": acc
|
||
})
|
||
|
||
except (ValueError, IndexError) as e:
|
||
continue
|
||
|
||
# 转换 set 为 count
|
||
stats["unique_proteins"] = len(stats["unique_proteins"])
|
||
stats["unique_domains"] = len(stats["unique_domains"])
|
||
|
||
return domains, stats
|
||
|
||
|
||
def run_r_plotting(tree_file: Path, subfamily_file: Optional[Path],
|
||
domtblout_file: Path, output_prefix: Path,
|
||
rscript_path: Optional[Path] = None,
|
||
evalue_cutoff: float = DEFAULT_EVALUE,
|
||
coverage_cutoff: Optional[float] = None) -> Dict:
|
||
"""
|
||
运行 R 绘图脚本 tree_domain.R
|
||
|
||
Args:
|
||
tree_file: 系统发育树文件 (Newick)
|
||
subfamily_file: 亚家族分类文件 (可选)
|
||
domtblout_file: hmmscan domtblout 文件
|
||
output_prefix: 输出文件前缀（不含扩展名）
|
||
rscript_path: R 脚本路径
|
||
evalue_cutoff: E-value 阈值
|
||
coverage_cutoff: 覆盖度阈值 (百分比，如 50 表示 50%)
|
||
|
||
Returns:
|
||
Dict with success, output_files, etc.
|
||
"""
|
||
if rscript_path is None:
|
||
# 从项目根目录推断
|
||
rscript_path = Path(__file__).parent.parent.parent / "scripts" / "tree_domain.R"
|
||
|
||
if not rscript_path.exists():
|
||
return {
|
||
"success": False,
|
||
"error": f"R script not found: {rscript_path}"
|
||
}
|
||
|
||
# 构建命令 - 适配 tree_domain.R 的参数格式
|
||
# 用法: Rscript tree_domain.R <treefile> <subfamily_file> <domtblout_file> [evalue] [coverage] [output_prefix]
|
||
cmd = [
|
||
"Rscript", str(rscript_path),
|
||
str(tree_file),
|
||
str(subfamily_file) if subfamily_file else "NONE", # R 脚本需要处理 NONE
|
||
str(domtblout_file),
|
||
str(evalue_cutoff),
|
||
str(coverage_cutoff) if coverage_cutoff else "", # 空字符串表示不启用
|
||
str(output_prefix)
|
||
]
|
||
|
||
print(f"[CMD] {' '.join(cmd)}")
|
||
|
||
start_time = datetime.now()
|
||
|
||
try:
|
||
result = subprocess.run(cmd, capture_output=True, text=True, check=False, timeout=600)
|
||
end_time = datetime.now()
|
||
|
||
success = result.returncode == 0
|
||
|
||
# 检查输出文件（R 脚本输出: {prefix}.pdf 和 {prefix}.png）
|
||
output_files = []
|
||
for ext in ['.pdf', '.png']:
|
||
f = Path(f"{output_prefix}{ext}")
|
||
if f.exists():
|
||
output_files.append(str(f))
|
||
|
||
return {
|
||
"success": success and len(output_files) > 0,
|
||
"returncode": result.returncode,
|
||
"stdout": result.stdout[-2000:],
|
||
"stderr": result.stderr[-2000:],
|
||
"output_files": output_files,
|
||
"elapsed_sec": (end_time - start_time).total_seconds()
|
||
}
|
||
|
||
except subprocess.TimeoutExpired:
|
||
return {
|
||
"success": False,
|
||
"error": "R plotting timeout (>600s)",
|
||
"elapsed_sec": 600
|
||
}
|
||
except Exception as e:
|
||
return {
|
||
"success": False,
|
||
"error": str(e),
|
||
"elapsed_sec": 0
|
||
}
|
||
|
||
|
||
def generate_report(prefix: str, metrics: Dict, parameters: Dict,
|
||
output_file: Path, artifacts: Dict, plot_success: bool) -> None:
|
||
"""生成 domain_report.md 报告"""
|
||
|
||
report_lines = [
|
||
f"# Domain Analysis Report: {prefix}",
|
||
"",
|
||
f"**Generated:** {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}",
|
||
f"**Status:** {'Success' if plot_success else 'Partial Success (Plotting Failed)'}",
|
||
"",
|
||
"## Summary",
|
||
"",
|
||
f"- **Total Proteins:** {metrics['proteins_total']}",
|
||
f"- **Proteins with Domains:** {metrics['proteins_with_domains']} ({metrics['domain_coverage_rate']*100:.1f}%)",
|
||
f"- **Total Domain Hits:** {metrics['total_domain_hits']}",
|
||
f"- **Significant Hits (E < {parameters.get('evalue_cutoff', DEFAULT_EVALUE)}):** {metrics['significant_domain_hits']}",
|
||
f"- **Unique Domain Types:** {metrics['unique_domain_types']}",
|
||
f"- **Average Domains per Protein:** {metrics['avg_domains_per_protein']:.2f}",
|
||
"",
|
||
"## Domain Statistics",
|
||
"",
|
||
]
|
||
|
||
# 添加 top domains 信息
|
||
if "top_domains" in metrics and metrics["top_domains"]:
|
||
report_lines.append("### Top Domains by Frequency")
|
||
report_lines.append("")
|
||
for domain, count in metrics["top_domains"]:
|
||
report_lines.append(f"- **{domain}:** {count} hits")
|
||
report_lines.append("")
|
||
|
||
report_lines.extend([
|
||
"## Parameters",
|
||
"",
|
||
f"- **E-value Cutoff:** {parameters.get('evalue_cutoff', DEFAULT_EVALUE)}",
|
||
f"- **Coverage Cutoff:** {parameters.get('coverage_cutoff', 'Not applied')}",
|
||
f"- **CPU Cores:** {parameters.get('cpu', 'N/A')}",
|
||
f"- **Pfam Database:** {parameters.get('pfam_hmm', 'N/A')}",
|
||
f"- **Protein Source:** {parameters.get('protein_source', 'N/A')}",
|
||
f"- **Tree File:** {parameters.get('tree_file', 'N/A')}",
|
||
f"- **Subfamily File:** {parameters.get('subfamily_file', 'Not provided')}",
|
||
"",
|
||
"## Output Files",
|
||
""
|
||
])
|
||
|
||
for key, artifact in artifacts.items():
|
||
if isinstance(artifact, dict):
|
||
report_lines.append(f"- **{key}:** `{artifact.get('path', 'N/A')}` ({artifact.get('type', 'unknown')})")
|
||
|
||
report_lines.extend([
|
||
"",
|
||
"## Visualization",
|
||
""
|
||
])
|
||
|
||
if plot_success:
|
||
report_lines.append("Protein domain architecture visualization was successfully generated with three panels:")
|
||
report_lines.append("1. **Phylogenetic Tree** - Gene relationships with subfamily coloring")
|
||
report_lines.append("2. **Subfamily Bars** - Visual connection between tree and domains")
|
||
report_lines.append("3. **Domain Architecture** - Pfam domain positions along protein sequences")
|
||
report_lines.append("")
|
||
report_lines.append("Output formats:")
|
||
report_lines.append("- **PDF**: Vector graphics for publication")
|
||
report_lines.append("- **PNG**: 300 DPI raster image for presentation")
|
||
else:
|
||
report_lines.append("**Warning:** Domain visualization failed. See logs for details.")
|
||
report_lines.append("")
|
||
report_lines.append("Raw domain data is available in the domtblout file for manual analysis.")
|
||
|
||
report_lines.extend([
|
||
"",
|
||
"---",
|
||
"*Generated by Gene Family Analysis Pipeline - Domain Stage*"
|
||
])
|
||
|
||
with open(output_file, 'w') as f:
|
||
f.write('\\n'.join(report_lines))
|
||
|
||
print(f"[INFO] Report generated: {output_file}")
|
||
|
||
|
||
def write_manifest(manifest: Dict, manifest_file: Path) -> None:
|
||
"""写入 manifest 文件"""
|
||
with open(manifest_file, 'w') as f:
|
||
json.dump(manifest, f, indent=2)
|
||
print(f"[INFO] Manifest: {manifest_file}")
|
||
|
||
|
||
def write_failure_manifest(workspace: Path, final_dir: Path, prefix: str,
|
||
job_id: str, run_id: str, start_time: datetime,
|
||
error_msg: str, partial_artifacts: Dict = None) -> None:
|
||
"""写入失败状态的 manifest"""
|
||
end_time = datetime.now()
|
||
|
||
manifest = {
|
||
"manifest_version": "1.0",
|
||
"job_id": job_id,
|
||
"stage": "domain",
|
||
"run_id": run_id,
|
||
"prefix": prefix,
|
||
"timestamp": end_time.isoformat(),
|
||
"status": "failed",
|
||
"start_time": start_time.isoformat(),
|
||
"end_time": end_time.isoformat(),
|
||
"elapsed_sec": (end_time - start_time).total_seconds(),
|
||
"cmd": [str(Path(__file__).name)],
|
||
"workspace_root": display_path(workspace_root_dir(workspace)),
|
||
"workspace": display_path(workspace),
|
||
"input_dir": display_path(current_input_dir(workspace)),
|
||
"artifacts": partial_artifacts or {},
|
||
"parameters": {
|
||
"stage": "domain",
|
||
"error": error_msg
|
||
},
|
||
"metrics": {
|
||
"error_type": "execution_failed"
|
||
}
|
||
}
|
||
|
||
manifest_file = final_dir / f"{prefix}.domain.manifest.json"
|
||
write_manifest(manifest, manifest_file)
|
||
print(f"[ERROR] Failure manifest written: {manifest_file}")
|
||
|
||
|
||
def quarantine_stale_unscoped_failure(final_dir: Path, prefix: str, stage: str, canonical_manifest: Path) -> None:
|
||
"""Move stale unscoped failure manifests aside after a scoped success is available."""
|
||
stale_manifest = final_dir / f"{prefix}.{stage}.manifest.json"
|
||
if stale_manifest == canonical_manifest or not stale_manifest.exists():
|
||
return
|
||
try:
|
||
with open(stale_manifest, "r", encoding="utf-8") as handle:
|
||
stale_data = json.load(handle)
|
||
except Exception:
|
||
return
|
||
if stale_data.get("status") != "failed":
|
||
return
|
||
if stale_data.get("scope_id") or stale_data.get("stage_context"):
|
||
return
|
||
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
|
||
archived = final_dir / f"{prefix}.{stage}.stale_failed.{timestamp}.json"
|
||
stale_manifest.rename(archived)
|
||
print(f"[INFO] Quarantined stale unscoped failure manifest: {archived}")
|
||
|
||
|
||
def calculate_domain_metrics(domains: List[Dict], total_proteins: int,
|
||
hmmscan_stats: Dict) -> Dict:
|
||
"""计算结构域相关的 metrics"""
|
||
|
||
# 统计每个 domain 的频率
|
||
domain_freq = {}
|
||
for d in domains:
|
||
name = d["target_name"]
|
||
domain_freq[name] = domain_freq.get(name, 0) + 1
|
||
|
||
# Top 10 domains
|
||
top_domains = sorted(domain_freq.items(), key=lambda x: x[1], reverse=True)[:10]
|
||
|
||
# 计算有 domain 的蛋白数量
|
||
proteins_with_domains = len(set(d["query_name"] for d in domains))
|
||
|
||
metrics = {
|
||
"proteins_total": total_proteins,
|
||
"proteins_with_domains": proteins_with_domains,
|
||
"domain_coverage_rate": round(proteins_with_domains / max(total_proteins, 1), 4),
|
||
"total_domain_hits": hmmscan_stats.get("total_hits", 0),
|
||
"significant_domain_hits": len(domains),
|
||
"unique_domain_types": len(domain_freq),
|
||
"avg_domains_per_protein": round(len(domains) / max(proteins_with_domains, 1), 2),
|
||
"top_domains": top_domains
|
||
}
|
||
|
||
return metrics
|
||
|
||
|
||
# ===================== 主函数 =====================
|
||
|
||
def main():
|
||
parser = argparse.ArgumentParser(
|
||
description='Domain stage - Protein domain analysis (Pfam HMM scan + visualization)'
|
||
)
|
||
parser.add_argument('--prefix', required=True,
|
||
help='Analysis prefix (e.g., Ft.SPL)')
|
||
parser.add_argument('--workspace', required=True,
|
||
help='Workspace directory (e.g., workspace/20260212_143523_2480772_ec5c6d)')
|
||
parser.add_argument('--proteins',
|
||
help='Protein FASTA file (auto-detect from Identify manifest if not specified)')
|
||
parser.add_argument('--scope-file',
|
||
help='Gene-family scope JSON from resolve_scope.py')
|
||
parser.add_argument('--analysis-route', default='', help='Explicit analysis route identity for manifest/status')
|
||
parser.add_argument('--member-set-id', default='', help='Explicit member-set identity for manifest/status')
|
||
parser.add_argument('--tree',
|
||
help='Tree file (auto-detect from Phylogeny manifest if not specified)')
|
||
parser.add_argument('--subfamily',
|
||
help='Subfamily classification file (explicit path; auto-detect only uses confirmed_subfamily.tsv)')
|
||
parser.add_argument('--pfam-base',
|
||
default=DEFAULT_PFAM_BASE,
|
||
help=f'Pfam database base directory (default: {DEFAULT_PFAM_BASE})')
|
||
parser.add_argument('--evalue', type=float,
|
||
default=DEFAULT_EVALUE,
|
||
help=f'E-value cutoff (default: {DEFAULT_EVALUE})')
|
||
parser.add_argument('--coverage', type=float,
|
||
help='Coverage cutoff percentage (optional, e.g., 50 for 50%%)')
|
||
parser.add_argument('--cpu', type=int, default=8,
|
||
help='CPU cores for hmmscan (default: 8)')
|
||
parser.add_argument('--rscript',
|
||
help='Path to tree_domain.R script (auto-detect if not specified)')
|
||
parser.add_argument('--log-file',
|
||
help='Log file path')
|
||
parser.add_argument('--run-id',
|
||
help='Run ID (auto-generated if not provided)')
|
||
|
||
args = parser.parse_args()
|
||
|
||
run_id = args.run_id or str(uuid.uuid4())[:8]
|
||
job_id = Path(args.workspace).name
|
||
|
||
start_time = datetime.now()
|
||
ws = Path(args.workspace)
|
||
|
||
# 验证 workspace 路径
|
||
if not ws.exists():
|
||
print(f"[ERROR] Workspace directory not found: {ws}")
|
||
sys.exit(1)
|
||
|
||
# 1. 创建标准目录结构（契约要求）
|
||
domain_dir = ws / "domain" / "analysis"
|
||
domain_dir.mkdir(parents=True, exist_ok=True)
|
||
final_dir = ws / "final"
|
||
final_dir.mkdir(parents=True, exist_ok=True)
|
||
|
||
scope_data: Optional[Dict] = None
|
||
scope_file: Optional[Path] = None
|
||
if args.scope_file:
|
||
scope_file = Path(args.scope_file)
|
||
if not scope_file.exists():
|
||
error_msg = f"Scope file not found: {args.scope_file}"
|
||
print(f"[ERROR] {error_msg}")
|
||
write_failure_manifest(ws, final_dir, args.prefix, job_id, run_id, start_time, error_msg)
|
||
sys.exit(1)
|
||
scope_data = read_scope(scope_file)
|
||
scope_id = str((scope_data or {}).get("scope_id") or "").strip()
|
||
output_prefix = user_output_prefix_for(
|
||
"domain",
|
||
args.prefix,
|
||
scope_data or {},
|
||
analysis_route=args.analysis_route,
|
||
member_set_id=args.member_set_id,
|
||
) if scope_id else args.prefix
|
||
|
||
print(f"[INFO] Domain stage started for {args.prefix}")
|
||
print(f"[INFO] Job ID: {job_id}, Run ID: {run_id}")
|
||
print(f"[INFO] Workspace: {display_path(ws)}")
|
||
|
||
# 2. 解析输入文件（支持命令行指定或自动查找）
|
||
|
||
# 2.1 蛋白序列文件
|
||
if scope_data:
|
||
try:
|
||
protein_file = resolve_scope_protein_fasta(ws, args.prefix, scope_data, args.proteins)
|
||
except FileNotFoundError as e:
|
||
error_msg = str(e)
|
||
print(f"[ERROR] {error_msg}")
|
||
write_failure_manifest(ws, final_dir, args.prefix, job_id, run_id, start_time, error_msg)
|
||
sys.exit(1)
|
||
elif args.proteins:
|
||
protein_file = resolve_job_path(ws, args.proteins) or Path(args.proteins)
|
||
if not protein_file.exists():
|
||
error_msg = f"Specified protein file not found: {protein_file}"
|
||
print(f"[ERROR] {error_msg}")
|
||
write_failure_manifest(ws, final_dir, args.prefix, job_id, run_id, start_time, error_msg)
|
||
sys.exit(1)
|
||
else:
|
||
phase = "subfamily_evidence_required" if args.analysis_route == "subfamily_evidence" else "post_confirm_required"
|
||
try:
|
||
resolved = resolve_effective_members(
|
||
workspace=ws,
|
||
final_dir=final_dir,
|
||
input_dir=current_input_dir(ws),
|
||
prefix=args.prefix,
|
||
phase=phase, # type: ignore[arg-type]
|
||
require="pep",
|
||
member_set_id=args.member_set_id,
|
||
)
|
||
if resolved.path_pep:
|
||
protein_file = resolved.path_pep
|
||
print(f"[INFO] Protein file from {resolved.source_kind}: {protein_file}")
|
||
else:
|
||
protein_file = Path("")
|
||
except MemberResolutionError as e:
|
||
if args.analysis_route == "subfamily_evidence":
|
||
error_msg = f"Complete identify member review before Domain subfamily evidence: {e}"
|
||
else:
|
||
error_msg = f"Complete subfamily confirmation before regular Domain analysis: {e}"
|
||
print(f"[ERROR] {error_msg}")
|
||
write_failure_manifest(ws, final_dir, args.prefix, job_id, run_id, start_time, error_msg)
|
||
sys.exit(1)
|
||
if protein_file.exists() and protein_file.stat().st_size > 0:
|
||
protein_file = ensure_protein_fasta_without_stop_markers(protein_file)
|
||
print(f"[INFO] Protein file: {protein_file}")
|
||
protein_record_count = count_fasta_records(protein_file)
|
||
else:
|
||
protein_file = Path("")
|
||
|
||
if scope_data or args.proteins:
|
||
protein_file = ensure_protein_fasta_without_stop_markers(protein_file)
|
||
print(f"[INFO] Protein file: {protein_file}")
|
||
protein_record_count = count_fasta_records(protein_file)
|
||
|
||
if not scope_data and not args.proteins and not protein_file:
|
||
candidates = [final_dir / f"{args.prefix}.final_members.pep.fa"]
|
||
if "." in args.prefix:
|
||
species, family = args.prefix.split(".", 1)
|
||
if species not in {"Target", "Multi", "Reference", "Comparative"} and family:
|
||
candidates.append(final_dir / f"Target.{species}.{family}.final_members.pep.fa")
|
||
# 从 Identify manifest 查找
|
||
identify_manifest = find_latest_manifest(final_dir, "identify")
|
||
if identify_manifest:
|
||
try:
|
||
manifest_data = read_manifest(identify_manifest)
|
||
artifacts = manifest_data.get("artifacts", {}) if isinstance(manifest_data.get("artifacts", {}), dict) else {}
|
||
for key in ("final_members", "final_proteins", "proteins"):
|
||
item = artifacts.get(key)
|
||
protein_path = manifest_artifact_path(ws, identify_manifest, item) if isinstance(item, dict) else None
|
||
if protein_path:
|
||
candidates.append(protein_path)
|
||
except Exception:
|
||
pass
|
||
protein_file = next((path for path in candidates if path.exists() and path.stat().st_size > 0), candidates[0])
|
||
|
||
if not protein_file.exists():
|
||
error_msg = "Protein file not found. Checked: " + ", ".join(str(path) for path in candidates) + ". Please run Identify stage first or specify --proteins."
|
||
print(f"[ERROR] {error_msg}")
|
||
write_failure_manifest(ws, final_dir, args.prefix, job_id, run_id, start_time, error_msg)
|
||
sys.exit(1)
|
||
protein_file = ensure_protein_fasta_without_stop_markers(protein_file)
|
||
print(f"[INFO] Protein file: {protein_file}")
|
||
protein_record_count = count_fasta_records(protein_file)
|
||
|
||
# 2.2 树文件
|
||
try:
|
||
tree_file = resolve_scope_tree_file(ws, args.prefix, scope_data, args.tree)
|
||
except FileNotFoundError as e:
|
||
error_msg = str(e)
|
||
print(f"[ERROR] {error_msg}")
|
||
write_failure_manifest(ws, final_dir, args.prefix, job_id, run_id, start_time, error_msg)
|
||
sys.exit(1)
|
||
|
||
print(f"[INFO] Tree file: {tree_file}")
|
||
|
||
# 2.3 亚家族文件（可选，但 R 脚本需要）
|
||
# 自动查找只接受 confirmed_subfamily；prediction/geneID 必须显式传入
|
||
subfamily_file = resolve_subfamily_file(ws, args.prefix, job_id, args.subfamily, output_prefix=output_prefix)
|
||
|
||
# 3. 准备 Pfam 数据库
|
||
print("[INFO] Preparing Pfam database...")
|
||
pfam_hmm, pfam_success, pfam_error = prepare_pfam_database(domain_dir, args.pfam_base)
|
||
|
||
if not pfam_success:
|
||
print(f"[ERROR] {pfam_error}")
|
||
write_failure_manifest(ws, final_dir, args.prefix, job_id, run_id, start_time, pfam_error)
|
||
sys.exit(1)
|
||
|
||
print(f"[INFO] Pfam ready: {pfam_hmm}")
|
||
|
||
# 4. 运行 hmmscan
|
||
print(f"[INFO] Running hmmscan (this may take a while)...")
|
||
hmm_success, hmm_result = run_hmmscan(
|
||
pfam_hmm, protein_file, domain_dir,
|
||
output_prefix=output_prefix,
|
||
cpu=args.cpu, evalue=args.evalue
|
||
)
|
||
|
||
if not hmm_success:
|
||
error_msg = f"hmmscan failed: {hmm_result.get('error', 'unknown error')}"
|
||
print(f"[ERROR] {error_msg}")
|
||
print(f"[ERROR] stderr: {hmm_result.get('stderr', '')}")
|
||
write_failure_manifest(ws, final_dir, args.prefix, job_id, run_id, start_time, error_msg)
|
||
sys.exit(1)
|
||
|
||
print(f"[INFO] hmmscan completed in {hmm_result['elapsed_sec']:.1f}s")
|
||
print(f"[INFO] domtblout: {hmm_result['domtblout']}")
|
||
|
||
# 5. 解析 domtblout
|
||
print(f"[INFO] Parsing domain hits (E < {args.evalue})...")
|
||
domains, hmm_stats = parse_domtblout(hmm_result['domtblout'], args.evalue)
|
||
|
||
print(f"[INFO] - Total hits: {hmm_stats['total_hits']}")
|
||
print(f"[INFO] - Significant hits: {hmm_stats['significant_hits']}")
|
||
print(f"[INFO] - Unique proteins: {hmm_stats['unique_proteins']}")
|
||
print(f"[INFO] - Unique domains: {hmm_stats['unique_domains']}")
|
||
|
||
if not domains:
|
||
print(f"[WARNING] No significant domain hits found with E-value < {args.evalue}")
|
||
# 不退出，继续生成空报告
|
||
|
||
# 6. 计算 metrics
|
||
# 获取总蛋白数量（从蛋白文件）
|
||
try:
|
||
from Bio import SeqIO
|
||
total_proteins = len(list(SeqIO.parse(protein_file, "fasta")))
|
||
except:
|
||
# 回退：从 domains 估算
|
||
total_proteins = hmm_stats.get("unique_proteins", 0)
|
||
|
||
domain_metrics = calculate_domain_metrics(domains, total_proteins, hmm_stats)
|
||
|
||
print(f"[INFO] Metrics calculated:")
|
||
for k, v in domain_metrics.items():
|
||
if k != "top_domains":
|
||
print(f" - {k}: {v}")
|
||
|
||
# 7. R 绘图（输出到 final/）
|
||
# 使用标准命名：{prefix}.domain（不含扩展名，R脚本会自动添加 .pdf 和 .png）
|
||
output_path_prefix = final_dir / f"{output_prefix}.domain"
|
||
|
||
# 确定 R 脚本路径
|
||
if args.rscript:
|
||
rscript_path = Path(args.rscript)
|
||
else:
|
||
rscript_path = Path(__file__).parent.parent.parent / "scripts" / "tree_domain.R"
|
||
|
||
plot_result = run_r_plotting(
|
||
tree_file,
|
||
subfamily_file,
|
||
hmm_result['domtblout'],
|
||
output_path_prefix,
|
||
rscript_path=rscript_path,
|
||
evalue_cutoff=args.evalue,
|
||
coverage_cutoff=args.coverage
|
||
)
|
||
|
||
plot_success = plot_result.get("success", False)
|
||
if not plot_success:
|
||
print(f"[WARNING] R plotting failed: {plot_result.get('error', plot_result.get('stderr', 'unknown'))}")
|
||
else:
|
||
print(f"[INFO] Plot files generated: {plot_result.get('output_files', [])}")
|
||
|
||
# 8. 准备 artifacts
|
||
end_time = datetime.now()
|
||
artifacts = {}
|
||
|
||
# 中间文件（domtblout, tblout）
|
||
if hmm_result['domtblout'] and hmm_result['domtblout'].exists():
|
||
artifacts["domtblout"] = {
|
||
"artifact": "domtblout",
|
||
"path": str(hmm_result['domtblout'].relative_to(ws)),
|
||
"full_path": display_path(hmm_result['domtblout']),
|
||
"type": "txt",
|
||
"size": hmm_result['domtblout'].stat().st_size,
|
||
"sha256": compute_file_hash(hmm_result['domtblout']),
|
||
"description": "hmmscan domain table output (raw)"
|
||
}
|
||
|
||
if hmm_result['tblout'] and hmm_result['tblout'].exists():
|
||
artifacts["tblout"] = {
|
||
"artifact": "tblout",
|
||
"path": str(hmm_result['tblout'].relative_to(ws)),
|
||
"full_path": display_path(hmm_result['tblout']),
|
||
"type": "txt",
|
||
"size": hmm_result['tblout'].stat().st_size,
|
||
"description": "hmmscan per-sequence hits (raw)"
|
||
}
|
||
|
||
# 最终产物（图片）- 同时支持 PNG 和 PDF
|
||
for img_file in plot_result.get("output_files", []):
|
||
img_path = Path(img_file)
|
||
if img_path.exists():
|
||
key = f"domain_plot_{img_path.suffix[1:]}" # domain_plot_png, domain_plot_pdf
|
||
artifacts[key] = {
|
||
"artifact": "domain_plot",
|
||
"path": str(img_path.relative_to(ws)),
|
||
"full_path": display_path(img_path),
|
||
"type": "image",
|
||
"format": img_path.suffix[1:],
|
||
"size": img_path.stat().st_size,
|
||
"sha256": compute_file_hash(img_path),
|
||
"description": f"Domain architecture visualization ({img_path.suffix[1:].upper()})"
|
||
}
|
||
|
||
# 9. 生成报告
|
||
report_file = final_dir / f"{output_prefix}.domain_report.md"
|
||
|
||
parameters = {
|
||
"cpu": args.cpu,
|
||
"prefix": args.prefix,
|
||
"output_prefix": output_prefix,
|
||
"stage": "domain",
|
||
"job_id": job_id,
|
||
"run_id": run_id,
|
||
"evalue_cutoff": args.evalue,
|
||
"coverage_cutoff": args.coverage,
|
||
"pfam_hmm": str(pfam_hmm),
|
||
"pfam_base": args.pfam_base,
|
||
"protein_source": str(protein_file),
|
||
"protein_sha256": compute_file_hash(protein_file),
|
||
"protein_record_count": protein_record_count,
|
||
"tree_file": str(tree_file),
|
||
"subfamily_file": str(subfamily_file) if subfamily_file else None,
|
||
"hmmscan_elapsed_sec": hmm_result.get("elapsed_sec", 0),
|
||
"Rscript_path": str(rscript_path),
|
||
"tool_versions": get_tool_versions(),
|
||
"timestamp_start": start_time.isoformat(),
|
||
"timestamp_end": end_time.isoformat(),
|
||
"workspace_root": display_path(workspace_root_dir(ws)),
|
||
"workspace": display_path(ws),
|
||
"input_dir": display_path(current_input_dir(ws)),
|
||
"domain_analysis_dir": display_path(domain_dir),
|
||
"final_dir": display_path(final_dir)
|
||
}
|
||
if scope_data:
|
||
parameters["scope_file"] = str(scope_file)
|
||
parameters["scope_id"] = scope_data.get("scope_id", "")
|
||
parameters["scope_mode"] = scope_data.get("scope_mode", "")
|
||
parameters["target_species"] = scope_data.get("target_species", [])
|
||
parameters["reference_species"] = scope_data.get("reference_species", [])
|
||
parameters["comparative_species"] = scope_data.get("comparative_species", [])
|
||
parameters["selected_species"] = scope_data.get("selected_species", [])
|
||
|
||
generate_report(output_prefix, domain_metrics, parameters, report_file, artifacts, plot_success)
|
||
|
||
if report_file.exists():
|
||
artifacts["domain_report"] = {
|
||
"artifact": "domain_report",
|
||
"path": str(report_file.relative_to(ws)),
|
||
"full_path": display_path(report_file),
|
||
"type": "report",
|
||
"format": "markdown",
|
||
"size": report_file.stat().st_size,
|
||
"sha256": compute_file_hash(report_file),
|
||
"description": "Human-readable analysis report"
|
||
}
|
||
|
||
# 10. 确定状态
|
||
# if not domains:
|
||
# status = "partial_success" # hmmscan 成功但无显著 hits
|
||
# elif plot_success:
|
||
# status = "success"
|
||
# else:
|
||
# status = "partial_success
|
||
# 检查输出文件是否存在
|
||
expected_png = final_dir / f"{output_prefix}.domain.png"
|
||
expected_pdf = final_dir / f"{output_prefix}.domain.pdf"
|
||
|
||
if plot_success and (expected_png.exists() or expected_pdf.exists()):
|
||
status = "success" # 只要绘图成功，就视为成功
|
||
else:
|
||
status = "partial_success"
|
||
|
||
# 11. 构建完整 manifest
|
||
manifest = {
|
||
"manifest_version": "1.0",
|
||
"job_id": job_id,
|
||
"stage": "domain",
|
||
"run_id": run_id,
|
||
"prefix": args.prefix,
|
||
"output_prefix": output_prefix,
|
||
"timestamp": end_time.isoformat(),
|
||
"status": status,
|
||
"start_time": start_time.isoformat(),
|
||
"end_time": end_time.isoformat(),
|
||
"elapsed_sec": (end_time - start_time).total_seconds(),
|
||
"cmd": [
|
||
str(Path(__file__).name),
|
||
"--prefix", args.prefix,
|
||
"--workspace", str(ws),
|
||
"--proteins", str(protein_file),
|
||
"--tree", str(tree_file),
|
||
"--pfam-base", args.pfam_base,
|
||
"--evalue", str(args.evalue),
|
||
"--cpu", str(args.cpu),
|
||
"--run-id", run_id
|
||
] + (["--scope-file", str(scope_file)] if scope_file else [])
|
||
+ (["--subfamily", str(subfamily_file)] if subfamily_file else [])
|
||
+ (["--coverage", str(args.coverage)] if args.coverage else [])
|
||
+ (["--rscript", str(rscript_path)] if args.rscript else []),
|
||
"workspace_root": display_path(workspace_root_dir(ws)),
|
||
"workspace": display_path(ws),
|
||
"input_dir": display_path(current_input_dir(ws)),
|
||
"artifacts": artifacts,
|
||
"parameters": parameters,
|
||
"metrics": domain_metrics
|
||
}
|
||
manifest["metrics"]["protein_record_count"] = protein_record_count
|
||
if scope_data:
|
||
manifest["scope_id"] = scope_data.get("scope_id", "")
|
||
manifest["scope_mode"] = scope_data.get("scope_mode", "")
|
||
manifest["scope_file"] = str(scope_file)
|
||
manifest["target_species"] = scope_data.get("target_species", [])
|
||
manifest["reference_species"] = scope_data.get("reference_species", [])
|
||
manifest["comparative_species"] = scope_data.get("comparative_species", [])
|
||
manifest["selected_species"] = scope_data.get("selected_species", [])
|
||
manifest["scope"] = {
|
||
"family": scope_data.get("family", ""),
|
||
"prefix": scope_data.get("prefix", ""),
|
||
"scope_id": scope_data.get("scope_id", ""),
|
||
"scope_mode": scope_data.get("scope_mode", ""),
|
||
"target_species": scope_data.get("target_species", []),
|
||
"reference_species": scope_data.get("reference_species", []),
|
||
"comparative_species": scope_data.get("comparative_species", []),
|
||
"selected_species": scope_data.get("selected_species", []),
|
||
}
|
||
context = build_stage_context(
|
||
stage="domain",
|
||
prefix=args.prefix,
|
||
scope_data=scope_data or {},
|
||
scope_file=scope_file or "",
|
||
output_prefix=output_prefix,
|
||
subfamily_file=subfamily_file,
|
||
analysis_route=args.analysis_route,
|
||
member_set_id=args.member_set_id,
|
||
)
|
||
manifest["tree_source"] = str(tree_file)
|
||
manifest["protein_fasta"] = str(protein_file)
|
||
apply_stage_context(manifest, context)
|
||
|
||
manifest_file = final_dir / f"{output_prefix}.domain.manifest.json"
|
||
write_manifest(manifest, manifest_file)
|
||
if scope_data and scope_data.get("scope_id"):
|
||
write_legacy_manifest_copy(manifest_file, args.prefix, str(scope_data.get("scope_id") or ""), "domain")
|
||
quarantine_stale_unscoped_failure(final_dir, args.prefix, "domain", manifest_file)
|
||
|
||
# 12. 完成
|
||
if status == "success":
|
||
print(f"[INFO] Domain stage completed successfully")
|
||
else:
|
||
print(f"[WARNING] Domain stage completed with status: {status}")
|
||
|
||
print("\\n=== OUTPUT_FILES ===")
|
||
for k, v in artifacts.items():
|
||
if isinstance(v, dict) and v.get("type") == "image":
|
||
print(f"IMAGE: {v['full_path']}")
|
||
|
||
sys.exit(0 if status == "success" else 1)
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|