#!/usr/bin/env python3 """ Domain Analysis Skill - 蛋白结构域分析技能 使用HMMER对Pfam数据库进行蛋白结构域分析 Usage: python run.py --input protein.fa --pfam PF03110.hmm --output ./domain python run.py --config config.yaml --output ./domain python run.py --demo --output ./domain """ import argparse import hashlib import json import logging import os import shutil import subprocess import sys from datetime import datetime from pathlib import Path from typing import Dict, List, Optional import yaml # 设置日志 logging.basicConfig( level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s' ) logger = logging.getLogger(__name__) class DomainSkill: """蛋白结构域分析技能主类""" def __init__(self, config_path: Optional[str] = None): self.config = self._load_config(config_path) self.commands = [] # 记录所有执行的命令 self.checksums = {} # 记录文件校验 def _load_config(self, config_path: Optional[str]) -> Dict: """加载配置文件""" if config_path and os.path.exists(config_path): with open(config_path) as f: return yaml.safe_load(f) default_config = Path(__file__).parent / "config.yaml" with open(default_config) as f: return yaml.safe_load(f) def _which(self, cmd: str) -> Optional[str]: """检查命令是否存在""" return shutil.which(cmd) def _run_command(self, cmd: List[str], description: str): """运行命令并记录""" logger.info(f"Running: {description}") logger.info(f"Command: {' '.join(cmd)}") self.commands.append({ "description": description, "command": " ".join(cmd), "timestamp": datetime.now().isoformat() }) subprocess.run(cmd, check=True) def _compute_checksum(self, filepath: Path) -> str: """计算文件SHA256校验""" sha256 = hashlib.sha256() with open(filepath, "rb") as f: for chunk in iter(lambda: f.read(8192), b""): sha256.update(chunk) return sha256.hexdigest() def _record_file(self, filepath: Path, file_type: str = "output"): """记录文件信息""" if filepath.exists(): self.checksums[str(filepath)] = { "sha256": self._compute_checksum(filepath), "size": filepath.stat().st_size, "type": file_type } def run_hmmscan(self, protein_fasta: Path, pfam_hmm: Path, output_dir: Path) -> Dict: """运行hmmscan Pfam搜索""" params = self.config["parameters"]["hmmscan"] cpu = self.config["parameters"]["cpu"] hmmscan = self._which("hmmscan") if not hmmscan: raise RuntimeError("hmmscan not found. Please install HMMER.") output_dir.mkdir(parents=True, exist_ok=True) prefix = "domain_scan" domtblout = output_dir / f"{prefix}.pfam.domtblout" tblout = output_dir / f"{prefix}.pfam.tblout" hmmscan_out = output_dir / f"{prefix}.pfam.out" cmd = [ "hmmscan", "--cpu", str(cpu), "--domtblout", str(domtblout), "--tblout", str(tblout), "-o", str(hmmscan_out), ] if params.get("use_cut_ga"): cmd.append("--cut_ga") cmd.extend([str(pfam_hmm), str(protein_fasta)]) self._run_command(cmd, "HMMER hmmscan Pfam搜索") self._record_file(domtblout) self._record_file(tblout) self._record_file(hmmscan_out) return { "domtblout": str(domtblout), "tblout": str(tblout) } def generate_report(self, output_dir: Path, evalue: float, coverage: float) -> Path: """生成分析报告""" report_file = output_dir / "report.md" report_content = f"""# 蛋白结构域分析报告 ## 基本信息 - **E-value阈值**: {evalue} - **HMM覆盖度阈值**: {coverage} - **分析时间**: {datetime.now().strftime("%Y-%m-%d %H:%M:%S")} - **技能版本**: {self.config["skill"]["version"]} ## 分析结果 | 文件 | 描述 | |------|------| | domain_scan.pfam.domtblout | Domain表格输出 | | domain_scan.pfam.tblout | 表格输出 | ## 可复现性 所有分析命令已记录在 `commands.sh` 中， 环境配置见 `environment.yml`。 ## 校验信息 文件SHA256校验值见 `checksums.sha256`。 """ with open(report_file, "w") as f: f.write(report_content) self._record_file(report_file) return report_file def write_reproducibility_bundle(self, output_dir: Path): """写入可复现性文件包""" # commands.sh commands_file = output_dir / "commands.sh" with open(commands_file, "w") as f: f.write("#!/bin/bash\n# 蛋白结构域分析复现命令\n\n") for cmd in self.commands: f.write(f"# {cmd['description']}\n") f.write(f"{cmd['command']}\n\n") os.chmod(commands_file, 0o755) # environment.yml env_file = output_dir / "environment.yml" env_config = { "name": "domain-analysis", "channels": self.config["dependencies"]["conda"]["channels"], "dependencies": self.config["dependencies"]["conda"]["packages"] } with open(env_file, "w") as f: yaml.dump(env_config, f, default_flow_style=False) # checksums.sha256 checksums_file = output_dir / "checksums.sha256" with open(checksums_file, "w") as f: for filepath, info in self.checksums.items(): f.write(f"{info['sha256']} {filepath}\n") def run(self, protein_fasta: Path, pfam_hmm: Path, output_dir: Path) -> Dict: """运行完整的蛋白结构域分析流程""" logger.info(f"Starting domain analysis") # 记录输入文件 self._record_file(protein_fasta, "input") self._record_file(pfam_hmm, "input") # 运行hmmscan results = self.run_hmmscan(protein_fasta, pfam_hmm, output_dir) # 生成报告 params = self.config["parameters"]["hmmscan"] report = self.generate_report( output_dir, params["evalue"], params["coverage"] ) # 写入可复现性文件 self.write_reproducibility_bundle(output_dir) logger.info(f"Domain analysis completed. Report: {report}") return { "status": "success", "output_dir": str(output_dir), "report": str(report), "domtblout": results["domtblout"] } def run_demo(self, output_dir: Path) -> Dict: """运行Demo模式""" logger.info("Running in demo mode") # 创建模拟数据 demo_dir = output_dir / "demo_data" demo_dir.mkdir(parents=True, exist_ok=True) # 模拟蛋白质FASTA protein_fasta = demo_dir / "demo_proteins.fa" with open(protein_fasta, "w") as f: f.write(">FtGene001\nMKTLLLPLLLLVLLLLLVQGSISAQ\n" ">FtGene002\nMKTLLLPLLLLVLLLLLVQGSISAQ\n" ">AtGene001\nMKTLLLPLLLLVLLLLLVQGSISAQ\n") # 模拟PFAM HMM文件 pfam_hmm = demo_dir / "demo_pfam.hmm" pfam_hmm.write_text("HMMER3/f [3.3.2 | ]\n") # 创建模拟输出文件 output_dir.mkdir(parents=True, exist_ok=True) # domtblout (模拟) domtblout = output_dir / "domain_scan.pfam.domtblout" domtblout.write_text("""# --- full sequence --- --- best 1 domain --- --- domain number estimation --- # target name accession tlen query name accession qlen E-value score bias c-Evalue i-Evalue score bias from to from to from to acc PF03110.16 SBP 383 FtGene001 - 284 1.2e-45 156.8 0.0 1.1e-47 1.2e-45 156.8 0.0 52 376 52 376 1 1 0.89 PF03110.16 SBP 383 FtGene002 - 284 1.5e-42 148.2 0.0 1.4e-44 1.5e-42 148.2 0.0 55 380 55 380 1 1 0.87 """) # tblout (模拟) tblout = output_dir / "domain_scan.pfam.tblout" tblout.write_text("""# target name accession tlen query name accession qlen E-value score bias PF03110.16 SBP 383 FtGene001 - 284 1.2e-45 156.8 0.0 PF03110.16 SBP 383 FtGene002 - 284 1.5e-42 148.2 0.0 """) # report.md report_file = output_dir / "report.md" report_file.write_text(f"""# 蛋白结构域分析报告 ## 基本信息 - **E-value阈值**: 1e-5 (demo) - **HMM覆盖度阈值**: 0.5 (demo) - **分析时间**: {datetime.now().strftime("%Y-%m-%d %H:%M:%S")} - **技能版本**: {self.config["skill"]["version"]} ## 分析结果 | 文件 | 描述 | |------|------| | domain_scan.pfam.domtblout | Domain表格输出 | | domain_scan.pfam.tblout | 表格输出 | ## 可复现性 所有分析命令已记录在 `commands.sh` 中，环境配置见 `environment.yml`。 ## 校验信息 文件SHA256校验值见 `checksums.sha256`。 """) # commands.sh commands_file = output_dir / "commands.sh" commands_file.write_text("""#!/bin/bash # 蛋白结构域分析复现命令 (Demo模式) hmmscan --cpu 8 --domtblout domain.pfam.domtblout --tblout domain.pfam.tblout -o domain.pfam.out PF03110.hmm proteins.fa """) # environment.yml env_file = output_dir / "environment.yml" env_file.write_text("""name: domain-analysis channels: - bioconda - conda-forge dependencies: - hmmer=3.3.2 - biopython=1.81 """) # checksums.sha256 checksums_file = output_dir / "checksums.sha256" for f in [domtblout, tblout, report_file, commands_file, env_file]: if f.exists(): sha256 = hashlib.sha256() with open(f, "rb") as fp: for chunk in iter(lambda: fp.read(8192), b""): sha256.update(chunk) checksums_file.write_text(f"{sha256.hexdigest()} {f.name}\n") logger.info(f"Domain demo completed. Report: {report_file}") return { "status": "success", "output_dir": str(output_dir), "report": str(report_file), "domtblout": str(domtblout) } def main(): parser = argparse.ArgumentParser( description="Domain Analysis Skill - 蛋白结构域分析" ) parser.add_argument("--input", "-i", type=Path, help="蛋白序列FASTA文件") parser.add_argument("--pfam", type=Path, help="Pfam HMM模型文件") parser.add_argument("--evalue", type=float, help="E-value阈值") parser.add_argument("--coverage", type=float, help="HMM覆盖度阈值") parser.add_argument("--cpu", type=int, help="CPU线程数") parser.add_argument("--output", "-o", type=Path, required=True, help="输出目录") parser.add_argument("--config", "-c", type=Path, help="配置文件路径") parser.add_argument("--demo", action="store_true", help="运行Demo模式") args = parser.parse_args() # 初始化技能 skill = DomainSkill(args.config) # 更新配置参数 if args.evalue: skill.config["parameters"]["hmmscan"]["evalue"] = args.evalue if args.coverage: skill.config["parameters"]["hmmscan"]["coverage"] = args.coverage if args.cpu: skill.config["parameters"]["cpu"] = args.cpu # 运行 if args.demo: result = skill.run_demo(args.output) else: if not args.input or not args.pfam: parser.error("--input and --pfam are required (unless using --demo)") result = skill.run(args.input, args.pfam, args.output) # 输出结果 print(json.dumps(result, indent=2, ensure_ascii=False)) return 0 if result["status"] == "success" else 1 if __name__ == "__main__": sys.exit(main())