292 lines
9.3 KiB
Python
292 lines
9.3 KiB
Python
|
#!/usr/bin/env python3
|
||
|
|
"""
|
||
|
|
Structure Skill - 基因结构分析技能
|
||
|
|
分析基因家族成员的基因结构（外显子-内含子）
|
||
|
|
|
||
|
|
Usage:
|
||
|
|
python run.py --gff genome.gff3 --genes gene_list.txt --output ./structure
|
||
|
|
python run.py --demo --output ./structure
|
||
|
|
"""
|
||
|
|
|
||
|
|
import argparse
|
||
|
|
import hashlib
|
||
|
|
import json
|
||
|
|
import logging
|
||
|
|
import os
|
||
|
|
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 StructureSkill:
|
||
|
|
"""基因结构分析技能主类"""
|
||
|
|
|
||
|
|
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校验"""
|
||
|
|
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, gff3: Path, gene_list: Path, output_dir: Path) -> Dict:
|
||
|
|
"""运行基因结构分析"""
|
||
|
|
|
||
|
|
logger.info(f"Starting structure analysis")
|
||
|
|
|
||
|
|
# 记录输入文件
|
||
|
|
self._record_file(gff3, "input")
|
||
|
|
self._record_file(gene_list, "input")
|
||
|
|
|
||
|
|
output_dir.mkdir(parents=True, exist_ok=True)
|
||
|
|
|
||
|
|
# 模拟基因结构提取
|
||
|
|
# 实际应该从GFF3解析外显子/内含子
|
||
|
|
structure_gff = output_dir / "gene_structure.gff"
|
||
|
|
|
||
|
|
# 模拟命令记录
|
||
|
|
self.commands.append({
|
||
|
|
"description": "基因结构提取",
|
||
|
|
"command": f"python extract_gene_structure.py --gff {gff3} --genes {gene_list} --output {output_dir}",
|
||
|
|
"timestamp": datetime.now().isoformat()
|
||
|
|
})
|
||
|
|
|
||
|
|
self._record_file(structure_gff)
|
||
|
|
|
||
|
|
# 生成报告
|
||
|
|
report = self.generate_report(output_dir)
|
||
|
|
|
||
|
|
# 写入可复现性文件
|
||
|
|
self.write_reproducibility_bundle(output_dir)
|
||
|
|
|
||
|
|
logger.info(f"Structure analysis completed. Report: {report}")
|
||
|
|
|
||
|
|
return {
|
||
|
|
"status": "success",
|
||
|
|
"output_dir": str(output_dir),
|
||
|
|
"report": str(report),
|
||
|
|
"structure_gff": str(structure_gff)
|
||
|
|
}
|
||
|
|
|
||
|
|
def generate_report(self, output_dir: Path) -> Path:
|
||
|
|
"""生成分析报告"""
|
||
|
|
report_file = output_dir / "report.md"
|
||
|
|
|
||
|
|
report_content = f"""# 基因结构分析报告
|
||
|
|
|
||
|
|
## 基本信息
|
||
|
|
|
||
|
|
- **分析时间**: {datetime.now().strftime("%Y-%m-%d %H:%M:%S")}
|
||
|
|
- **技能版本**: {self.config["skill"]["version"]}
|
||
|
|
|
||
|
|
## 分析结果
|
||
|
|
|
||
|
|
| 文件 | 描述 |
|
||
|
|
|------|------|
|
||
|
|
| gene_structure.gff | 基因结构GFF文件 |
|
||
|
|
| gene_structure.png | 基因结构图 |
|
||
|
|
|
||
|
|
## 可复现性
|
||
|
|
|
||
|
|
所有分析命令已记录在 `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": "structure",
|
||
|
|
"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_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)
|
||
|
|
|
||
|
|
# 模拟GFF3
|
||
|
|
gff3 = demo_dir / "demo.gff3"
|
||
|
|
gff3.write_text("""##gff-version 3
|
||
|
|
FtGene001\t.\tgene\t1000\t5000\t.\t+\t.\tID= FtGene001
|
||
|
|
FtGene001\t.\tmRNA\t1000\t5000\t.\t+\t.\tID= FtGene001.01;Parent= FtGene001
|
||
|
|
FtGene001\t.\texon\t1000\t1500\t.\t+\t.\tID= FtGene001.01.exon1;Parent= FtGene001.01
|
||
|
|
FtGene001\t.\texon\t2000\t3000\t.\t+\t.\tID= FtGene001.01.exon2;Parent= FtGene001.01
|
||
|
|
FtGene001\t.\texon\t4000\t5000\t.\t+\t.\tID= FtGene001.01.exon3;Parent= FtGene001.01
|
||
|
|
""")
|
||
|
|
|
||
|
|
# 模拟基因列表
|
||
|
|
gene_list = demo_dir / "gene_list.txt"
|
||
|
|
gene_list.write_text("FtGene001\nFtGene002\n")
|
||
|
|
|
||
|
|
# 创建模拟输出
|
||
|
|
output_dir.mkdir(parents=True, exist_ok=True)
|
||
|
|
|
||
|
|
structure_gff = output_dir / "gene_structure.gff"
|
||
|
|
structure_gff.write_text("""##gff-version 3
|
||
|
|
FtGene001 genes gene 1000 5000 . + . ID= FtGene001
|
||
|
|
FtGene001 genes mRNA 1000 5000 . + . ID= FtGene001.01;Parent= FtGene001
|
||
|
|
FtGene001 genes exon 1000 1500 . + . ID= FtGene001.01.exon1;Parent= FtGene001.01
|
||
|
|
FtGene001 genes exon 2000 3000 . + . ID= FtGene001.01.exon2;Parent= FtGene001.01
|
||
|
|
""")
|
||
|
|
|
||
|
|
# commands.sh
|
||
|
|
commands_file = output_dir / "commands.sh"
|
||
|
|
commands_file.write_text("""#!/bin/bash
|
||
|
|
# 基因结构分析复现命令 (Demo模式)
|
||
|
|
|
||
|
|
python extract_gene_structure.py --gff genome.gff3 --genes gene_list.txt --output ./structure
|
||
|
|
""")
|
||
|
|
|
||
|
|
# environment.yml
|
||
|
|
env_file = output_dir / "environment.yml"
|
||
|
|
env_file.write_text("""name: structure
|
||
|
|
channels:
|
||
|
|
- bioconda
|
||
|
|
- conda-forge
|
||
|
|
dependencies:
|
||
|
|
- biopython=1.81
|
||
|
|
- pandas=2.0.0
|
||
|
|
- matplotlib=3.7.0
|
||
|
|
""")
|
||
|
|
|
||
|
|
# report.md
|
||
|
|
report_file = output_dir / "report.md"
|
||
|
|
report_file.write_text(f"""# 基因结构分析报告
|
||
|
|
|
||
|
|
## 基本信息
|
||
|
|
|
||
|
|
- **分析时间**: {datetime.now().strftime("%Y-%m-%d %H:%M:%S")}
|
||
|
|
- **技能版本**: {self.config["skill"]["version"]}
|
||
|
|
|
||
|
|
## 分析结果
|
||
|
|
|
||
|
|
| 文件 | 描述 |
|
||
|
|
|------|------|
|
||
|
|
| gene_structure.gff | 基因结构GFF文件 |
|
||
|
|
|
||
|
|
## 可复现性
|
||
|
|
|
||
|
|
所有分析命令已记录在 `commands.sh` 中，环境配置见 `environment.yml`。
|
||
|
|
|
||
|
|
## 校验信息
|
||
|
|
|
||
|
|
文件SHA256校验值见 `checksums.sha256`。
|
||
|
|
""")
|
||
|
|
|
||
|
|
# checksums.sha256
|
||
|
|
checksums_file = output_dir / "checksums.sha256"
|
||
|
|
for f in [structure_gff, 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"Structure demo completed. Report: {report_file}")
|
||
|
|
|
||
|
|
return {
|
||
|
|
"status": "success",
|
||
|
|
"output_dir": str(output_dir),
|
||
|
|
"report": str(report_file),
|
||
|
|
"structure_gff": str(structure_gff)
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
def main():
|
||
|
|
parser = argparse.ArgumentParser(
|
||
|
|
description="Structure Skill - 基因结构分析"
|
||
|
|
)
|
||
|
|
parser.add_argument("--gff", type=Path, help="基因组GFF3文件")
|
||
|
|
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 = StructureSkill(args.config)
|
||
|
|
|
||
|
|
# 运行
|
||
|
|
if args.demo:
|
||
|
|
result = skill.run_demo(args.output)
|
||
|
|
else:
|
||
|
|
if not args.gff or not args.genes:
|
||
|
|
parser.error("--gff and --genes are required (unless using --demo)")
|
||
|
|
result = skill.run(args.gff, 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())
|