#!/usr/bin/env python3 """ Promoter Analysis Skill - 启动子分析技能 使用PlantCARE分析启动子顺式作用元件 Usage: python run.py --genome genome.fa --genes gene_list.txt --output ./promoter python run.py --demo --output ./promoter """ import argparse import hashlib import json import logging import os import sys from datetime import datetime from pathlib import Path from typing import Dict, Optional import yaml # 设置日志 logging.basicConfig( level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s' ) logger = logging.getLogger(__name__) class PromoterSkill: """启动子分析技能主类""" 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 _compute_checksum(self, filepath: Path) -> str: 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(self, genome_fasta: Path, gene_list: Path, output_dir: Path) -> Dict: logger.info(f"Starting promoter analysis") self._record_file(genome_fasta, "input") self._record_file(gene_list, "input") output_dir.mkdir(parents=True, exist_ok=True) elements_file = output_dir / "promoter_elements.txt" self.commands.append({ "description": "启动子元件预测", "command": f"python analyze_promoter.py --genome {genome_fasta} --genes {gene_list} --output {output_dir}", "timestamp": datetime.now().isoformat() }) self._record_file(elements_file) report = self.generate_report(output_dir) self.write_reproducibility_bundle(output_dir) logger.info(f"Promoter analysis completed. Report: {report}") return { "status": "success", "output_dir": str(output_dir), "report": str(report), "elements_file": str(elements_file) } def generate_report(self, output_dir: Path) -> Path: report_file = output_dir / "report.md" params = self.config["parameters"] report_content = f"""# 启动子分析报告 ## 基本信息 - **启动子长度**: {params["promoter_length"]} bp - **物种类型**: {params["species"]} - **分析时间**: {datetime.now().strftime("%Y-%m-%d %H:%M:%S")} - **技能版本**: {self.config["skill"]["version"]} ## 分析结果 | 文件 | 描述 | |------|------| | promoter_elements.txt | 启动子元件预测结果 | ## 可复现性 所有分析命令已记录在 `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_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) env_file = output_dir / "environment.yml" env_config = { "name": "promoter-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_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_demo(self, output_dir: Path) -> Dict: logger.info("Running in demo mode") demo_dir = output_dir / "demo_data" demo_dir.mkdir(parents=True, exist_ok=True) genome_fasta = demo_dir / "demo_genome.fa" genome_fasta.write_text(">chr1\nATCGATCGATCGATCGATCG\n") gene_list = demo_dir / "gene_list.txt" gene_list.write_text("gene001\ngene002\n") output_dir.mkdir(parents=True, exist_ok=True) elements_file = output_dir / "promoter_elements.txt" elements_file.write_text("""Gene Element Position Strand Sequence Description gene001 TATA-box -150 + TATA TATA box gene001 CAAT-box -200 + CCAAT CAAT box gene002 G-box -300 - CACGTG G-box """) commands_file = output_dir / "commands.sh" commands_file.write_text("""#!/bin/bash # 启动子分析复现命令 (Demo模式) python analyze_promoter.py --genome genome.fa --genes gene_list.txt --output ./promoter """) env_file = output_dir / "environment.yml" env_file.write_text("""name: promoter-analysis channels: - bioconda - conda-forge dependencies: - biopython=1.81 - pandas=2.0.0 - matplotlib=3.7.0 """) report_file = output_dir / "report.md" report_file.write_text(f"""# 启动子分析报告 ## 基本信息 - **启动子长度**: 2000 bp (demo) - **物种类型**: plant (demo) - **分析时间**: {datetime.now().strftime("%Y-%m-%d %H:%M:%S")} - **技能版本**: {self.config["skill"]["version"]} ## 分析结果 | 文件 | 描述 | |------|------| | promoter_elements.txt | 启动子元件预测结果 | ## 可复现性 所有分析命令已记录在 `commands.sh` 中，环境配置见 `environment.yml`。 ## 校验信息 文件SHA256校验值见 `checksums.sha256`。 """) checksums_file = output_dir / "checksums.sha256" for f in [elements_file, 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"Promoter demo completed. Report: {report_file}") return { "status": "success", "output_dir": str(output_dir), "report": str(report_file), "elements_file": str(elements_file) } def main(): parser = argparse.ArgumentParser(description="Promoter Analysis Skill") parser.add_argument("--genome", type=Path, help="基因组FASTA文件") parser.add_argument("--genes", type=Path, help="基因ID列表文件") 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 = PromoterSkill(args.config) if args.demo: result = skill.run_demo(args.output) else: if not args.genome or not args.genes: parser.error("--genome and --genes are required (unless using --demo)") result = skill.run(args.genome, args.genes, 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())