437 lines
15 KiB
Python
Executable File
437 lines
15 KiB
Python
Executable File
#!/usr/bin/env python3
|
|
"""
|
|
Phylogeny Skill - 系统发育分析技能
|
|
支持MAFFT/Muscle多序列比对和IQ-TREE/FastTree建树
|
|
|
|
Usage:
|
|
python run.py --input protein.fa --output ./phylogeny
|
|
python run.py --config config.yaml --output ./phylogeny
|
|
python run.py --demo --output ./phylogeny
|
|
"""
|
|
|
|
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 PhylogenySkill:
|
|
"""系统发育分析技能主类"""
|
|
|
|
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,
|
|
dry_run: bool = False, capture: bool = False) -> str:
|
|
"""运行命令并记录"""
|
|
logger.info(f"Running: {description}")
|
|
logger.info(f"Command: {' '.join(cmd)}")
|
|
|
|
self.commands.append({
|
|
"description": description,
|
|
"command": " ".join(cmd),
|
|
"timestamp": datetime.now().isoformat()
|
|
})
|
|
|
|
if dry_run:
|
|
logger.info("Dry run mode - command not executed")
|
|
return ""
|
|
|
|
if capture:
|
|
result = subprocess.run(
|
|
cmd,
|
|
capture_output=True,
|
|
text=True,
|
|
check=True
|
|
)
|
|
return result.stdout
|
|
else:
|
|
subprocess.run(cmd, check=True)
|
|
return ""
|
|
|
|
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_alignment(self, protein_fasta: Path, output_dir: Path,
|
|
dry_run: bool = False) -> Path:
|
|
"""运行多序列比对"""
|
|
aligner = self.config["parameters"]["aligner"]["tool"]
|
|
aln_out = output_dir / "alignment.aln"
|
|
|
|
params = self.config["parameters"]
|
|
cpu = params["cpu"]
|
|
|
|
if aligner == "mafft":
|
|
mafft = self._which("mafft")
|
|
if not mafft:
|
|
raise RuntimeError("mafft not found. Please install mafft.")
|
|
|
|
cmd = [mafft, "--auto", "--thread", str(cpu), str(protein_fasta)]
|
|
self._run_command(cmd, "MAFFT多序列比对", dry_run)
|
|
|
|
if dry_run:
|
|
# 模拟比对输出
|
|
aln_out.write_text(f">gene1\nMKTLLLPLLL-\n>gene2\nMKTLLLPLLL-\n>gene3\nMKTLLLPLLL-\n")
|
|
else:
|
|
with open(aln_out, "w") as w:
|
|
subprocess.run(cmd, check=True, stdout=w)
|
|
|
|
elif aligner == "muscle":
|
|
muscle = self._which("muscle")
|
|
if not muscle:
|
|
raise RuntimeError("muscle not found. Please install muscle.")
|
|
|
|
cmd = [muscle, "-align", str(protein_fasta), "-output", str(aln_out)]
|
|
self._run_command(cmd, "Muscle多序列比对", dry_run)
|
|
|
|
if dry_run:
|
|
aln_out.write_text(f">gene1\nMKTLLLPLLL-\n>gene2\nMKTLLLPLLL-\n>gene3\nMKTLLLPLLL-\n")
|
|
else:
|
|
raise ValueError(f"Unknown aligner: {aligner}")
|
|
|
|
self._record_file(aln_out)
|
|
return aln_out
|
|
|
|
def run_tree(self, aln_file: Path, output_dir: Path,
|
|
dry_run: bool = False) -> Path:
|
|
"""运行系统发育树构建"""
|
|
tree_builder = self.config["parameters"]["tree_builder"]["tool"]
|
|
tree_out = output_dir / "phylogeny.tree"
|
|
|
|
params = self.config["parameters"]
|
|
cpu = params["cpu"]
|
|
|
|
if tree_builder == "iqtree":
|
|
iq = self._which("iqtree2") or self._which("iqtree")
|
|
if not iq:
|
|
raise RuntimeError("iqtree2/iqtree not found. Please install IQ-TREE.")
|
|
|
|
prefix = str(output_dir / "phylogeny.iqtree")
|
|
cmd = [iq, "-s", str(aln_file), "-nt", str(cpu), "-pre", prefix]
|
|
self._run_command(cmd, "IQ-TREE系统发育树构建", dry_run)
|
|
|
|
if dry_run:
|
|
# 模拟树输出
|
|
tree_out.write_text("(gene1:0.1,gene2:0.1,gene3:0.1);\n")
|
|
else:
|
|
subprocess.run(cmd, check=True)
|
|
treefile = Path(prefix + ".treefile")
|
|
if treefile.exists():
|
|
tree_out.write_text(treefile.read_text())
|
|
else:
|
|
raise RuntimeError(f"IQ-TREE treefile not found: {treefile}")
|
|
|
|
elif tree_builder == "fasttree":
|
|
ft = self._which("FastTree") or self._which("fasttree")
|
|
if not ft:
|
|
raise RuntimeError("FastTree not found. Please install FastTree.")
|
|
|
|
cmd = [ft, str(aln_file)]
|
|
self._run_command(cmd, "FastTree系统发育树构建", dry_run)
|
|
|
|
if dry_run:
|
|
tree_out.write_text("(gene1:0.1,gene2:0.1,gene3:0.1);\n")
|
|
else:
|
|
with open(tree_out, "w") as w:
|
|
subprocess.run(cmd, check=True, stdout=w)
|
|
else:
|
|
raise ValueError(f"Unknown tree_builder: {tree_builder}")
|
|
|
|
self._record_file(tree_out)
|
|
return tree_out
|
|
|
|
def generate_report(self, output_dir: Path, aligner: str,
|
|
tree_builder: str) -> Path:
|
|
"""生成分析报告"""
|
|
report_file = output_dir / "report.md"
|
|
|
|
report_content = f"""# 系统发育分析报告
|
|
|
|
## 基本信息
|
|
|
|
- **比对工具**: {aligner}
|
|
- **建树工具**: {tree_builder}
|
|
- **分析时间**: {datetime.now().strftime("%Y-%m-%d %H:%M:%S")}
|
|
- **技能版本**: {self.config["skill"]["version"]}
|
|
|
|
## 分析结果
|
|
|
|
| 文件 | 描述 |
|
|
|------|------|
|
|
| alignment.aln | 多序列比对结果 |
|
|
| phylogeny.tree | 系统发育树 (Newick格式) |
|
|
|
|
## 可复现性
|
|
|
|
所有分析命令已记录在 `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": "phylogeny",
|
|
"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, output_dir: Path) -> Dict:
|
|
"""运行完整的系统发育分析流程"""
|
|
|
|
logger.info(f"Starting phylogeny analysis")
|
|
|
|
# 创建输出目录
|
|
output_dir.mkdir(parents=True, exist_ok=True)
|
|
|
|
# 记录输入文件
|
|
self._record_file(protein_fasta, "input")
|
|
|
|
# 1. 多序列比对
|
|
aligner = self.config["parameters"]["aligner"]["tool"]
|
|
dry_run = aligner == "mafft" and not self._which("mafft")
|
|
aln_file = self.run_alignment(protein_fasta, output_dir, dry_run)
|
|
|
|
# 2. 建树
|
|
tree_builder = self.config["parameters"]["tree_builder"]["tool"]
|
|
dry_run = dry_run or (tree_builder == "iqtree" and not (self._which("iqtree2") or self._which("iqtree")))
|
|
tree_file = self.run_tree(aln_file, output_dir, dry_run)
|
|
|
|
# 3. 生成报告
|
|
report = self.generate_report(output_dir, aligner, tree_builder)
|
|
|
|
# 4. 写入可复现性文件
|
|
self.write_reproducibility_bundle(output_dir)
|
|
|
|
logger.info(f"Phylogeny analysis completed. Report: {report}")
|
|
|
|
return {
|
|
"status": "success",
|
|
"output_dir": str(output_dir),
|
|
"report": str(report),
|
|
"alignment": str(aln_file),
|
|
"tree": str(tree_file)
|
|
}
|
|
|
|
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"
|
|
">AtGene002\nMKTLLLPLLLLVLLLLLVQGSISAQ\n")
|
|
|
|
# 创建输出目录
|
|
output_dir.mkdir(parents=True, exist_ok=True)
|
|
|
|
# 直接创建模拟输出文件（不依赖外部工具）
|
|
# alignment.aln
|
|
aln_file = output_dir / "alignment.aln"
|
|
aln_file.write_text(">FtGene001\nMKTLLLPLLLLVLLLLLVQGSISAQ\n"
|
|
">FtGene002\nMKTLLLPLLLLVLLLLLVQGSISAQ\n"
|
|
">AtGene001\nMKTLLLPLLLLVLLLLLVQGSISAQ\n"
|
|
">AtGene002\nMKTLLLPLLLLVLLLLLVQGSISAQ\n")
|
|
|
|
# phylogeny.tree
|
|
tree_file = output_dir / "phylogeny.tree"
|
|
tree_file.write_text("((FtGene001:0.1,FtGene002:0.1):0.05,(AtGene001:0.1,AtGene002:0.1):0.05);\n")
|
|
|
|
# report.md
|
|
report_file = output_dir / "report.md"
|
|
report_file.write_text(f"""# 系统发育分析报告
|
|
|
|
## 基本信息
|
|
|
|
- **比对工具**: mafft (demo)
|
|
- **建树工具**: iqtree (demo)
|
|
- **分析时间**: {datetime.now().strftime("%Y-%m-%d %H:%M:%S")}
|
|
- **技能版本**: {self.config["skill"]["version"]}
|
|
|
|
## 分析结果
|
|
|
|
| 文件 | 描述 |
|
|
|------|------|
|
|
| alignment.aln | 多序列比对结果 |
|
|
| phylogeny.tree | 系统发育树 (Newick格式) |
|
|
|
|
## 可复现性
|
|
|
|
所有分析命令已记录在 `commands.sh` 中，环境配置见 `environment.yml`。
|
|
|
|
## 校验信息
|
|
|
|
文件SHA256校验值见 `checksums.sha256`。
|
|
""")
|
|
|
|
# commands.sh (模拟)
|
|
commands_file = output_dir / "commands.sh"
|
|
commands_file.write_text("""#!/bin/bash
|
|
# 系统发育分析复现命令 (Demo模式)
|
|
|
|
# MAFFT多序列比对
|
|
mafft --auto --thread 8 input_proteins.fa > alignment.aln
|
|
|
|
# IQ-TREE系统发育树构建
|
|
iqtree -s alignment.aln -nt 8 -pre phylogeny.iqtree
|
|
""")
|
|
|
|
# environment.yml (模拟)
|
|
env_file = output_dir / "environment.yml"
|
|
env_file.write_text("""name: phylogeny
|
|
channels:
|
|
- bioconda
|
|
- conda-forge
|
|
dependencies:
|
|
- mafft=7.505
|
|
- iqtree=2.2.0
|
|
- biopython=1.81
|
|
- pandas=2.0.0
|
|
""")
|
|
|
|
# checksums.sha256 (模拟)
|
|
checksums_file = output_dir / "checksums.sha256"
|
|
for f in [aln_file, tree_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"Phylogeny demo completed. Report: {report_file}")
|
|
|
|
return {
|
|
"status": "success",
|
|
"output_dir": str(output_dir),
|
|
"report": str(report_file),
|
|
"alignment": str(aln_file),
|
|
"tree": str(tree_file)
|
|
}
|
|
|
|
|
|
def main():
|
|
parser = argparse.ArgumentParser(
|
|
description="Phylogeny Skill - 系统发育分析"
|
|
)
|
|
parser.add_argument("--input", "-i", type=Path, help="蛋白序列FASTA文件")
|
|
parser.add_argument("--aligner", choices=["mafft", "muscle"],
|
|
help="比对工具")
|
|
parser.add_argument("--tree-builder", choices=["iqtree", "fasttree"],
|
|
help="建树工具")
|
|
parser.add_argument("--cpu", type=int, default=8, 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 = PhylogenySkill(args.config)
|
|
|
|
# 更新配置参数
|
|
if args.aligner:
|
|
skill.config["parameters"]["aligner"]["tool"] = args.aligner
|
|
if args.tree_builder:
|
|
skill.config["parameters"]["tree_builder"]["tool"] = args.tree_builder
|
|
if args.cpu:
|
|
skill.config["parameters"]["cpu"] = args.cpu
|
|
|
|
# 运行
|
|
if args.demo:
|
|
result = skill.run_demo(args.output)
|
|
else:
|
|
if not args.input:
|
|
parser.error("--input is required (unless using --demo)")
|
|
result = skill.run(args.input, 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()) |