1716 lines
76 KiB
Python
Executable File
1716 lines
76 KiB
Python
Executable File
#!/usr/bin/env python3
|
||
"""
|
||
Interspecies Synteny阶段：物种间共线性分析（完整版）
|
||
|
||
输入：
|
||
- 参考物种蛋白序列 (final/{prefix}.final_members.pep.fa)
|
||
- 参考物种GFF3
|
||
- 目标物种蛋白序列
|
||
- 目标物种GFF3
|
||
|
||
输出：
|
||
- final/{prefix}.interspecies_tree.nwk (物种间进化树)
|
||
- final/{prefix}.interspecies_tree.pdf/png (进化树图)
|
||
- final/{prefix}.interspecies_synteny.pdf/png (共线性图)
|
||
- final/{prefix}.interspecies_synteny_report.md (分析报告)
|
||
- final/{prefix}.interspecies_synteny.manifest.json (机器清单)
|
||
"""
|
||
|
||
import os
|
||
import argparse
|
||
import json
|
||
import hashlib
|
||
import subprocess
|
||
import sys
|
||
import time
|
||
import shutil
|
||
import re
|
||
from datetime import datetime
|
||
from pathlib import Path
|
||
from typing import Dict, List, Any, Optional, Tuple
|
||
|
||
|
||
# ??conda??bin???PATH
|
||
CONDA_BIN = "/home/xia-lab/miniconda3/envs/test4py311/bin"
|
||
if CONDA_BIN not in os.environ.get("PATH", ""):
|
||
os.environ["PATH"] = CONDA_BIN + ":" + os.environ.get("PATH", "")
|
||
def pdf_to_png(pdf_path: Path, png_path: Path, dpi: int = 300) -> bool:
|
||
"""将PDF转换为PNG"""
|
||
try:
|
||
# 方法1: 使用pdftoppm (poppler-utils)
|
||
cmd = ["pdftoppm", "-png", "-r", str(dpi), "-singlefile",
|
||
str(pdf_path), str(png_path.with_suffix(''))]
|
||
result = subprocess.run(cmd, capture_output=True, text=True, timeout=60)
|
||
if result.returncode == 0:
|
||
return True
|
||
except Exception:
|
||
pass
|
||
|
||
try:
|
||
# 方法2: 使用convert (ImageMagick)
|
||
cmd = ["convert", "-density", str(dpi), str(pdf_path), "-quality", "100", str(png_path)]
|
||
result = subprocess.run(cmd, capture_output=True, text=True, timeout=60)
|
||
if result.returncode == 0:
|
||
return True
|
||
except Exception:
|
||
pass
|
||
|
||
return False
|
||
|
||
|
||
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 Exception:
|
||
return ""
|
||
|
||
|
||
def get_tool_versions() -> Dict[str, str]:
|
||
"""获取所需工具版本信息"""
|
||
versions = {}
|
||
|
||
for tool, cmd in [
|
||
("mafft", ["mafft", "--version"]),
|
||
("iqtree", ["iqtree", "--version"]),
|
||
("blastp", ["blastp", "-version"]),
|
||
("MCScanX", ["MCScanX", "-h"]),
|
||
]:
|
||
try:
|
||
result = subprocess.run(cmd, capture_output=True, text=True, check=False, timeout=5)
|
||
versions[tool] = "available" if result.returncode in [0, 1] else "not_found"
|
||
except Exception:
|
||
versions[tool] = "not_found"
|
||
|
||
return versions
|
||
|
||
|
||
def prepare_directories(workspace: Path) -> Dict[str, Path]:
|
||
"""创建目录结构"""
|
||
paths = {
|
||
"analysis_dir": workspace / "interspecies_synteny" / "analysis",
|
||
"final_dir": workspace / "final",
|
||
"logs_dir": workspace / "logs",
|
||
}
|
||
for name, path in paths.items():
|
||
path.mkdir(parents=True, exist_ok=True)
|
||
return paths
|
||
|
||
|
||
def merge_protein_sequences(
|
||
target_protein: Path,
|
||
ref_protein: Path,
|
||
output_fasta: Path
|
||
) -> Dict[str, Any]:
|
||
"""Merge target and reference family proteins for the interspecies tree."""
|
||
print("Merging family protein sequences...")
|
||
|
||
try:
|
||
with open(output_fasta, 'w') as outfile:
|
||
with open(target_protein, 'r') as f:
|
||
outfile.write(f.read().rstrip() + "\n")
|
||
with open(ref_protein, 'r') as f:
|
||
outfile.write(f.read().rstrip() + "\n")
|
||
|
||
target_count = sum(1 for line in open(target_protein) if line.startswith('>'))
|
||
ref_count = sum(1 for line in open(ref_protein) if line.startswith('>'))
|
||
|
||
print(f" merged target/ref family proteins: {target_count} + {ref_count} = {target_count + ref_count}")
|
||
|
||
return {
|
||
"success": True,
|
||
"output_file": output_fasta,
|
||
"target_count": target_count,
|
||
"ref_count": ref_count
|
||
}
|
||
except Exception as e:
|
||
return {"success": False, "stderr": str(e)}
|
||
|
||
|
||
def build_phylogenetic_tree(
|
||
merged_fasta: Path,
|
||
output_dir: Path,
|
||
prefix: str,
|
||
cpu: int = 8
|
||
) -> Dict[str, Any]:
|
||
"""构建系统发育树 (MAFFT + IQ-TREE)"""
|
||
print(f"构建进化树...")
|
||
|
||
# MAFFT比对
|
||
msa_file = output_dir / f"{prefix}_interspecies.mafft.msa.fa"
|
||
mafft_cmd = ["mafft", "--auto", "--thread", str(cpu), str(merged_fasta)]
|
||
|
||
try:
|
||
with open(msa_file, "w") as f_out:
|
||
result = subprocess.run(mafft_cmd, stdout=f_out, stderr=subprocess.PIPE,
|
||
text=True, check=False)
|
||
|
||
if result.returncode != 0:
|
||
return {"success": False, "error": f"MAFFT失败: {result.stderr}"}
|
||
|
||
# IQ-TREE建树
|
||
base_name = output_dir / f"{prefix}_interspecies.iqtree"
|
||
iqtree_cmd = [
|
||
"iqtree", "-s", str(msa_file), "-m", "MFP", "-bb", "1000",
|
||
"-nt", str(cpu), "-pre", str(base_name)
|
||
]
|
||
|
||
result = subprocess.run(iqtree_cmd, capture_output=True, text=True, check=False)
|
||
|
||
# IQ-TREE 生成的 treefile 是 {prefix}.iqtree.treefile，不是 {prefix}.treefile
|
||
tree_file = Path(f"{base_name}.treefile")
|
||
|
||
return {
|
||
"success": tree_file.exists(),
|
||
"tree_file": tree_file if tree_file.exists() else None,
|
||
"msa_file": msa_file
|
||
}
|
||
except Exception as e:
|
||
return {"success": False, "error": str(e)}
|
||
|
||
|
||
def filter_gff_files(
|
||
ref_gff3: Path,
|
||
target_gff3: Path,
|
||
output_dir: Path,
|
||
ref_species: str,
|
||
target_species: str
|
||
) -> Dict[str, Any]:
|
||
"""过滤GFF3文件 - 提取基因ID（去掉转录本版本后缀如 .1, .2）"""
|
||
print(f"过滤GFF3文件...")
|
||
|
||
import re
|
||
|
||
results = {}
|
||
for gff3_file, species, prefix in [
|
||
(ref_gff3, ref_species, "ref"),
|
||
(target_gff3, target_species, "target")
|
||
]:
|
||
output_gff = output_dir / f"{species}.gff"
|
||
|
||
try:
|
||
with open(gff3_file, 'r') as f:
|
||
lines = f.readlines()
|
||
|
||
count = 0
|
||
seen_genes = set() # 去重
|
||
|
||
with open(output_gff, 'w') as f:
|
||
for line in lines:
|
||
parts = line.strip().split('\t')
|
||
if len(parts) >= 9 and parts[2] == "mRNA":
|
||
attr = parts[8]
|
||
# 提取ID（如 ID=AT2G47070.1），保留版本号与BLAST结果匹配
|
||
gene_id = attr.split(';')[0].replace('ID=', '')
|
||
|
||
# 处理染色体名称：去掉 "Chr" 前缀（Ath用1,2,3而非Chr1,Chr2）
|
||
chr_name = parts[0]
|
||
if chr_name.startswith('Chr'):
|
||
chr_name = chr_name[3:] # Chr1 -> 1
|
||
|
||
# 去重：一个基因可能有多个转录本
|
||
if gene_id not in seen_genes:
|
||
seen_genes.add(gene_id)
|
||
f.write(f"{chr_name}\t{gene_id}\t{parts[3]}\t{parts[4]}\n")
|
||
count += 1
|
||
|
||
print(f" {species}: {count} 个基因（去重后）")
|
||
results[prefix] = {"success": True, "output_file": output_gff, "gene_count": count}
|
||
except Exception as e:
|
||
results[prefix] = {"success": False, "stderr": str(e)}
|
||
|
||
if results.get("ref", {}).get("success") and results.get("target", {}).get("success"):
|
||
return {
|
||
"success": True,
|
||
"ref_gff": results["ref"]["output_file"],
|
||
"target_gff": results["target"]["output_file"]
|
||
}
|
||
return {"success": False}
|
||
|
||
|
||
def run_reciprocal_blast(
|
||
ref_protein: Path,
|
||
target_protein: Path,
|
||
output_dir: Path,
|
||
ref_species: str,
|
||
target_species: str,
|
||
cpu: int = 8
|
||
) -> Dict[str, Any]:
|
||
"""双向BLAST比对 - 结果中基因名去掉版本后缀以匹配GFF"""
|
||
import re
|
||
print(f"双向BLAST比对...")
|
||
|
||
# 构建数据库
|
||
for protein_file in [ref_protein, target_protein]:
|
||
subprocess.run(["makeblastdb", "-in", str(protein_file), "-dbtype", "prot"],
|
||
capture_output=True, check=False)
|
||
|
||
# 双向blastp
|
||
blast1_out = output_dir / f"tmp.{ref_species}_{target_species}.blast"
|
||
subprocess.run([
|
||
"blastp", "-query", str(ref_protein), "-db", str(target_protein),
|
||
"-outfmt", "6", "-out", str(blast1_out), "-evalue", "1e-10",
|
||
"-num_alignments", "5", "-num_threads", str(cpu)
|
||
], capture_output=True, check=False)
|
||
|
||
blast2_out = output_dir / f"tmp.{target_species}_{ref_species}.blast"
|
||
subprocess.run([
|
||
"blastp", "-query", str(target_protein), "-db", str(ref_protein),
|
||
"-outfmt", "6", "-out", str(blast2_out), "-evalue", "1e-10",
|
||
"-num_alignments", "5", "-num_threads", str(cpu)
|
||
], capture_output=True, check=False)
|
||
|
||
# 合并BLAST结果（保留版本号与GFF匹配）
|
||
merged_blast = output_dir / f"{ref_species}_{target_species}.blast"
|
||
with open(merged_blast, 'w') as outfile:
|
||
for blast_file in [blast1_out, blast2_out]:
|
||
if blast_file.exists():
|
||
with open(blast_file, 'r') as infile:
|
||
for line in infile:
|
||
outfile.write(line) # 直接写入，保留版本号
|
||
|
||
return {"success": True, "output_file": merged_blast}
|
||
|
||
|
||
def run_mcscanx_interspecies(
|
||
work_dir: Path,
|
||
ref_species: str,
|
||
target_species: str
|
||
) -> Dict[str, Any]:
|
||
"""MCScanX种间共线性 - 需要合并的GFF文件"""
|
||
print(f"MCScanX种间共线性...")
|
||
|
||
prefix = f"{ref_species}_{target_species}"
|
||
|
||
# MCScanX 需要合并的 GFF 文件 (prefix.gff)
|
||
merged_gff = work_dir / f"{prefix}.gff"
|
||
ref_gff = work_dir / f"{ref_species}.gff"
|
||
target_gff = work_dir / f"{target_species}.gff"
|
||
|
||
if not merged_gff.exists() and ref_gff.exists() and target_gff.exists():
|
||
print(f" 合并 GFF 文件: {ref_gff.name} + {target_gff.name} -> {merged_gff.name}")
|
||
with open(merged_gff, 'w') as outfile:
|
||
for gff_file in [target_gff, ref_gff]: # 目标物种在前，参考物种在后
|
||
with open(gff_file, 'r') as infile:
|
||
outfile.write(infile.read())
|
||
|
||
cmd = ["MCScanX", f"./{prefix}"]
|
||
|
||
result = subprocess.run(cmd, capture_output=True, text=True,
|
||
cwd=str(work_dir), check=False)
|
||
|
||
collinearity_file = work_dir / f"{prefix}.collinearity"
|
||
|
||
return {
|
||
"success": collinearity_file.exists(),
|
||
"collinearity_file": collinearity_file if collinearity_file.exists() else None
|
||
}
|
||
|
||
|
||
|
||
|
||
def fasta_ids(fasta_file: Path) -> List[str]:
|
||
ids = []
|
||
with open(fasta_file, "r") as f:
|
||
for line in f:
|
||
if line.startswith(">"):
|
||
ids.append(line[1:].strip().split()[0])
|
||
return ids
|
||
|
||
|
||
def display_species_code(species: str, gene_id: str = "") -> str:
|
||
if re.match(r"^AT[1-5CM]G", gene_id, re.IGNORECASE):
|
||
return "AT"
|
||
return species
|
||
|
||
|
||
def species_for_gene(gene_id: str, species_list: List[str], default_species: str) -> str:
|
||
for species in sorted(species_list, key=len, reverse=True):
|
||
if gene_id.lower().startswith(species.lower()):
|
||
return display_species_code(species, gene_id)
|
||
if re.match(r"^AT[1-5CM]G", gene_id, re.IGNORECASE):
|
||
return "AT"
|
||
return display_species_code(default_species, gene_id)
|
||
|
||
|
||
def write_gene_species_map(target_fasta: Path, ref_fasta: Path, target_species: List[str], ref_species: List[str], output_file: Path) -> Path:
|
||
output_file.parent.mkdir(parents=True, exist_ok=True)
|
||
default_target = target_species[0] if target_species else "Target"
|
||
default_ref = ref_species[0] if ref_species else "Ref"
|
||
with open(output_file, "w") as out:
|
||
out.write("gene_id\tspecies\tgroup\n")
|
||
for gene_id in fasta_ids(target_fasta):
|
||
out.write(f"{gene_id}\t{species_for_gene(gene_id, target_species, default_target)}\ttarget\n")
|
||
for gene_id in fasta_ids(ref_fasta):
|
||
out.write(f"{gene_id}\t{species_for_gene(gene_id, ref_species, default_ref)}\tref\n")
|
||
return output_file
|
||
def normalize_seq_id(text: str) -> str:
|
||
"""Normalize FASTA/GFF IDs while preserving biologically meaningful versions."""
|
||
value = (text or "").strip()
|
||
if value.startswith(">"):
|
||
value = value[1:]
|
||
value = value.split()[0] if value.split() else value
|
||
if "|" in value:
|
||
value = value.split("|", 1)[0]
|
||
for prefix in ("gene:", "transcript:", "CDS:", "rna-", "gene-", "cds-"):
|
||
if value.startswith(prefix):
|
||
value = value[len(prefix):]
|
||
return value.strip().strip('"')
|
||
|
||
|
||
def id_variants(gene_id: str) -> set:
|
||
norm = normalize_seq_id(gene_id)
|
||
variants = {norm}
|
||
for prefix in ("gene:", "transcript:", "CDS:", "rna-", "gene-", "cds-"):
|
||
variants.add(normalize_seq_id(prefix + norm))
|
||
if "|" in gene_id:
|
||
variants.add(normalize_seq_id(gene_id.split("|", 1)[0]))
|
||
return {v for v in variants if v}
|
||
|
||
|
||
def fasta_ids(fasta_file: Path) -> List[str]:
|
||
ids = []
|
||
with open(fasta_file, "r") as f:
|
||
for line in f:
|
||
if line.startswith(">"):
|
||
seq_id = normalize_seq_id(line[1:])
|
||
if seq_id:
|
||
ids.append(seq_id)
|
||
return ids
|
||
|
||
|
||
def fasta_id_set(fasta_file: Path) -> set:
|
||
ids = set()
|
||
for seq_id in fasta_ids(fasta_file):
|
||
ids.update(id_variants(seq_id))
|
||
return ids
|
||
|
||
|
||
def display_species_code(species: str, gene_id: str = "") -> str:
|
||
if re.match(r"^AT[1-5CM]G", gene_id, re.IGNORECASE):
|
||
return "AT"
|
||
return species
|
||
|
||
|
||
def species_for_gene(gene_id: str, species_list: List[str], default_species: str, species_id_sets: Optional[Dict[str, set]] = None) -> str:
|
||
variants = id_variants(gene_id)
|
||
if species_id_sets:
|
||
for species, ids in species_id_sets.items():
|
||
if variants & ids:
|
||
return display_species_code(species, gene_id)
|
||
for species in sorted(species_list, key=len, reverse=True):
|
||
if normalize_seq_id(gene_id).lower().startswith(species.lower()):
|
||
return display_species_code(species, gene_id)
|
||
if re.match(r"^AT[1-5CM]G", gene_id, re.IGNORECASE):
|
||
return "AT"
|
||
return display_species_code(default_species, gene_id)
|
||
|
||
|
||
def collect_species_id_sets(groups: Dict[str, Any]) -> Dict[str, set]:
|
||
species_id_sets: Dict[str, set] = {}
|
||
for group_name in ("target", "ref"):
|
||
for species, item in groups.get(group_name, {}).items():
|
||
pep = item.get("resolved", {}).get("pep")
|
||
if pep and Path(pep).exists():
|
||
species_id_sets[species] = fasta_id_set(Path(pep))
|
||
return species_id_sets
|
||
|
||
|
||
def write_gene_species_map(
|
||
target_fasta: Path,
|
||
ref_fasta: Path,
|
||
target_species: List[str],
|
||
ref_species: List[str],
|
||
output_file: Path,
|
||
species_id_sets: Optional[Dict[str, set]] = None,
|
||
) -> Path:
|
||
output_file.parent.mkdir(parents=True, exist_ok=True)
|
||
default_target = target_species[0] if target_species else "Target"
|
||
default_ref = ref_species[0] if ref_species else "Ref"
|
||
with open(output_file, "w") as out:
|
||
out.write("gene_id\tspecies\tgroup\n")
|
||
for gene_id in fasta_ids(target_fasta):
|
||
species = species_for_gene(gene_id, target_species, default_target, species_id_sets)
|
||
out.write(f"{gene_id}\t{species}\ttarget\n")
|
||
for gene_id in fasta_ids(ref_fasta):
|
||
species = species_for_gene(gene_id, ref_species, default_ref, species_id_sets)
|
||
out.write(f"{gene_id}\t{species}\tref\n")
|
||
return output_file
|
||
|
||
|
||
def run_reciprocal_blast(
|
||
target_protein: Path,
|
||
ref_protein: Path,
|
||
output_dir: Path,
|
||
target_species: str,
|
||
ref_species: str,
|
||
cpu: int = 8
|
||
) -> Dict[str, Any]:
|
||
"""Run pairwise reciprocal BLASTP for one target-ref species pair."""
|
||
print(f"Reciprocal BLASTP: {target_species} vs {ref_species}")
|
||
|
||
pair_prefix = f"{target_species}.{ref_species}"
|
||
target_db = output_dir / f"{target_species}.blastdb"
|
||
ref_db = output_dir / f"{ref_species}.blastdb"
|
||
commands = []
|
||
|
||
for protein_file, db_prefix in [(target_protein, target_db), (ref_protein, ref_db)]:
|
||
cmd = ["makeblastdb", "-in", str(protein_file), "-dbtype", "prot", "-out", str(db_prefix)]
|
||
commands.append(" ".join(cmd))
|
||
result = subprocess.run(cmd, capture_output=True, text=True, check=False)
|
||
if result.returncode != 0:
|
||
return {"success": False, "stderr": result.stderr, "commands": commands}
|
||
|
||
target_to_ref = output_dir / f"tmp.{target_species}_{ref_species}.blast"
|
||
ref_to_target = output_dir / f"tmp.{ref_species}_{target_species}.blast"
|
||
for query, db_prefix, out_file in [(target_protein, ref_db, target_to_ref), (ref_protein, target_db, ref_to_target)]:
|
||
cmd = [
|
||
"blastp", "-query", str(query), "-db", str(db_prefix),
|
||
"-outfmt", "6", "-out", str(out_file), "-evalue", "1e-10",
|
||
"-num_alignments", "5", "-num_threads", str(cpu)
|
||
]
|
||
commands.append(" ".join(cmd))
|
||
result = subprocess.run(cmd, capture_output=True, text=True, check=False)
|
||
if result.returncode != 0:
|
||
return {"success": False, "stderr": result.stderr, "commands": commands}
|
||
|
||
merged_blast = output_dir / f"{pair_prefix}.blast"
|
||
with open(merged_blast, 'w') as outfile:
|
||
for blast_file in [target_to_ref, ref_to_target]:
|
||
if blast_file.exists():
|
||
with open(blast_file, 'r') as infile:
|
||
for line in infile:
|
||
if line.strip():
|
||
outfile.write(line)
|
||
|
||
return {"success": merged_blast.exists() and merged_blast.stat().st_size > 0, "output_file": merged_blast, "commands": commands}
|
||
|
||
|
||
def run_mcscanx_interspecies(work_dir: Path, target_species: str, ref_species: str) -> Dict[str, Any]:
|
||
"""Run MCScanX for one target-ref pair."""
|
||
pair_prefix = f"{target_species}.{ref_species}"
|
||
print(f"MCScanX pair: {pair_prefix}")
|
||
|
||
merged_gff = work_dir / f"{pair_prefix}.gff"
|
||
target_gff = work_dir / f"{target_species}.gff"
|
||
ref_gff = work_dir / f"{ref_species}.gff"
|
||
|
||
if target_gff.exists() and ref_gff.exists():
|
||
with open(merged_gff, 'w') as outfile:
|
||
for gff_file in [target_gff, ref_gff]:
|
||
with open(gff_file, 'r') as infile:
|
||
outfile.write(infile.read().rstrip() + "\n")
|
||
else:
|
||
return {"success": False, "error": f"Missing simplified GFF for {pair_prefix}"}
|
||
|
||
result = subprocess.run(["MCScanX", f"./{pair_prefix}"], capture_output=True, text=True, cwd=str(work_dir), check=False)
|
||
collinearity_file = work_dir / f"{pair_prefix}.collinearity"
|
||
return {
|
||
"success": collinearity_file.exists() and collinearity_file.stat().st_size > 0,
|
||
"collinearity_file": collinearity_file if collinearity_file.exists() else None,
|
||
"stdout": result.stdout,
|
||
"stderr": result.stderr,
|
||
"returncode": result.returncode,
|
||
}
|
||
|
||
|
||
def plot_interspecies_tree(
|
||
tree_file: Path,
|
||
output_dir: Path,
|
||
prefix: str,
|
||
target_species: str,
|
||
ref_species: str,
|
||
target_subfamily_file: Optional[Path] = None,
|
||
ref_subfamily_file: Optional[Path] = None,
|
||
gene_species_map: Optional[Path] = None
|
||
) -> Dict[str, Any]:
|
||
"""绘制物种间进化树图 - 使用专用的 interspecies_tree.R 脚本"""
|
||
print(f"绘制进化树图...")
|
||
|
||
# 查找专用的种间进化树R脚本
|
||
r_script = Path("scripts/interspecies_tree.R")
|
||
if not r_script.exists():
|
||
alt_script = Path("/home/xia-lab/wutx/00test/bbav2/scripts/interspecies_tree.R")
|
||
if alt_script.exists():
|
||
r_script = alt_script
|
||
else:
|
||
return {"success": False, "error": f"R脚本不存在: {r_script}"}
|
||
|
||
try:
|
||
# 使用专用的 interspecies_tree.R 脚本，传递物种参数和亚家族文件
|
||
# 使用conda环境的Rscript
|
||
conda_rscript = "/home/xia-lab/miniconda3/envs/test4py311/bin/Rscript"
|
||
cmd = [
|
||
conda_rscript, str(r_script),
|
||
str(tree_file),
|
||
str(output_dir),
|
||
target_species,
|
||
ref_species
|
||
]
|
||
|
||
# 添加亚家族文件参数（如果存在）
|
||
if target_subfamily_file and target_subfamily_file.exists():
|
||
cmd.append(str(target_subfamily_file))
|
||
print(f" 目标物种亚家族文件: {target_subfamily_file}")
|
||
else:
|
||
cmd.append("") # 空字符串表示无文件
|
||
|
||
if ref_subfamily_file and ref_subfamily_file.exists():
|
||
cmd.append(str(ref_subfamily_file))
|
||
print(f" 参考物种亚家族文件: {ref_subfamily_file}")
|
||
else:
|
||
cmd.append("") # 空字符串表示无文件
|
||
if gene_species_map and gene_species_map.exists():
|
||
cmd.append(str(gene_species_map))
|
||
print(f" Gene-to-species map: {gene_species_map}")
|
||
else:
|
||
cmd.append("")
|
||
result = subprocess.run(cmd, capture_output=True, text=True, check=False, timeout=120)
|
||
|
||
# 调试输出
|
||
print(f" R脚本返回码: {result.returncode}")
|
||
if result.stdout:
|
||
print(f" R脚本输出: {result.stdout[:500]}")
|
||
if result.stderr:
|
||
print(f" R脚本错误: {result.stderr[:500]}")
|
||
|
||
# 专用脚本输出固定为 interspecies_tree.pdf/png
|
||
output_pdf = output_dir / "interspecies_tree.pdf"
|
||
output_png = output_dir / "interspecies_tree.png"
|
||
|
||
# 重命名为规范名称
|
||
final_pdf = output_dir / f"{prefix}.interspecies_tree.pdf"
|
||
final_png = output_dir / f"{prefix}.interspecies_tree.png"
|
||
|
||
if output_pdf.exists():
|
||
shutil.move(output_pdf, final_pdf)
|
||
print(f" 进化树PDF: {final_pdf}")
|
||
if output_png.exists():
|
||
shutil.move(output_png, final_png)
|
||
print(f" 进化树PNG: {final_png}")
|
||
|
||
# PDF转PNG（如果没有生成PNG）
|
||
if final_pdf.exists() and not final_png.exists():
|
||
png_from_pdf = output_dir / "interspecies_tree.png"
|
||
if pdf_to_png(final_pdf, png_from_pdf):
|
||
final_png = output_dir / f"{prefix}.interspecies_tree.png"
|
||
shutil.move(png_from_pdf, final_png)
|
||
print(f" PDF转PNG: {final_png}")
|
||
|
||
return {
|
||
"success": final_pdf.exists() or final_png.exists(),
|
||
"pdf_file": final_pdf if final_pdf.exists() else None,
|
||
"png_file": final_png if final_png.exists() else None,
|
||
"stdout": result.stdout,
|
||
"stderr": result.stderr
|
||
}
|
||
except Exception as e:
|
||
return {"success": False, "error": str(e)}
|
||
|
||
|
||
def convert_blast_to_family_simple(
|
||
work_dir: Path,
|
||
blast_file: Path,
|
||
target_species: str,
|
||
ref_species: str,
|
||
) -> Dict[str, Any]:
|
||
"""Fallback: convert pairwise BLAST hits directly to JCVI anchors/simple.
|
||
|
||
This is not a strict MCScanX collinearity block. It is used when the ref side
|
||
only has family proteins, so MCScanX cannot satisfy whole-genome block rules.
|
||
"""
|
||
pair_prefix = f"{target_species}.{ref_species}"
|
||
target_bed = work_dir / f"{target_species}.bed"
|
||
ref_bed = work_dir / f"{ref_species}.bed"
|
||
anchors_file = work_dir / f"{pair_prefix}.anchors"
|
||
simple_file = work_dir / f"{pair_prefix}.anchors.simple"
|
||
|
||
valid_target = {line.rstrip("\n").split("\t")[3] for line in open(target_bed) if len(line.rstrip("\n").split("\t")) >= 4}
|
||
valid_ref = {line.rstrip("\n").split("\t")[3] for line in open(ref_bed) if len(line.rstrip("\n").split("\t")) >= 4}
|
||
seen = set()
|
||
links = []
|
||
if not blast_file.exists():
|
||
return {"success": False, "error": "BLAST file missing for family-only fallback"}
|
||
with open(blast_file, "r") as f:
|
||
for line in f:
|
||
if not line.strip() or line.startswith("#"):
|
||
continue
|
||
parts = line.rstrip("\n").split("\t")
|
||
if len(parts) < 2:
|
||
continue
|
||
a = normalize_seq_id(parts[0])
|
||
b = normalize_seq_id(parts[1])
|
||
target_gene = ref_gene = None
|
||
if a in valid_target and b in valid_ref:
|
||
target_gene, ref_gene = a, b
|
||
elif b in valid_target and a in valid_ref:
|
||
target_gene, ref_gene = b, a
|
||
if not target_gene or not ref_gene:
|
||
continue
|
||
key = (target_gene, ref_gene)
|
||
if key in seen:
|
||
continue
|
||
seen.add(key)
|
||
score_text = parts[11] if len(parts) > 11 else "1000"
|
||
try:
|
||
score = str(int(float(score_text)))
|
||
except Exception:
|
||
score = "1000"
|
||
links.append((target_gene, ref_gene, score))
|
||
if not links:
|
||
return {"success": False, "error": "No BLAST hits matched pair BED IDs for family-only fallback"}
|
||
with open(anchors_file, "w") as out:
|
||
out.write("##\n")
|
||
for target_gene, ref_gene, score in links:
|
||
out.write(f"{target_gene}\t{ref_gene}\t{score}\n")
|
||
with open(simple_file, "w") as out:
|
||
for target_gene, ref_gene, score in links:
|
||
out.write(f"{target_gene}\t{target_gene}\t{ref_gene}\t{ref_gene}\t{score}\t+\n")
|
||
return {
|
||
"success": True,
|
||
"anchors_file": anchors_file,
|
||
"simple_file": simple_file,
|
||
"anchor_count": len(links),
|
||
"fallback": "family_only_blast_links",
|
||
}
|
||
|
||
|
||
def convert_collinearity_to_simple(
|
||
work_dir: Path,
|
||
collinearity_file: Path,
|
||
target_species: str,
|
||
ref_species: str,
|
||
) -> Dict[str, Any]:
|
||
"""Convert one MCScanX .collinearity file into JCVI anchors/simple files."""
|
||
conda_python = "/home/xia-lab/miniconda3/envs/test4py311/bin/python"
|
||
pair_prefix = f"{target_species}.{ref_species}"
|
||
target_bed = work_dir / f"{target_species}.bed"
|
||
ref_bed = work_dir / f"{ref_species}.bed"
|
||
anchors_file = work_dir / f"{pair_prefix}.anchors"
|
||
anchors_new_file = work_dir / f"{pair_prefix}.anchors.new"
|
||
simple_file = work_dir / f"{pair_prefix}.anchors.simple"
|
||
|
||
valid_target = {line.rstrip("\n").split("\t")[3] for line in open(target_bed) if len(line.rstrip("\n").split("\t")) >= 4}
|
||
valid_ref = {line.rstrip("\n").split("\t")[3] for line in open(ref_bed) if len(line.rstrip("\n").split("\t")) >= 4}
|
||
|
||
current_block = []
|
||
anchors = 0
|
||
filtered = 0
|
||
with open(collinearity_file, "r") as f_in, open(anchors_file, "w") as f_out:
|
||
for line in f_in:
|
||
if not line.strip():
|
||
continue
|
||
if line.startswith("## Alignment"):
|
||
if current_block:
|
||
f_out.write("##\n")
|
||
f_out.writelines(current_block)
|
||
current_block = []
|
||
continue
|
||
if line.startswith("#"):
|
||
continue
|
||
parts = line.strip().split("\t")
|
||
if len(parts) < 3:
|
||
continue
|
||
gene_a = normalize_seq_id(parts[1])
|
||
gene_b = normalize_seq_id(parts[2])
|
||
target_gene = ref_gene = None
|
||
if gene_a in valid_target and gene_b in valid_ref:
|
||
target_gene, ref_gene = gene_a, gene_b
|
||
elif gene_b in valid_target and gene_a in valid_ref:
|
||
target_gene, ref_gene = gene_b, gene_a
|
||
if target_gene and ref_gene:
|
||
current_block.append(f"{target_gene} {ref_gene} 1000\n")
|
||
anchors += 1
|
||
else:
|
||
filtered += 1
|
||
if current_block:
|
||
f_out.write("##\n")
|
||
f_out.writelines(current_block)
|
||
|
||
if anchors == 0:
|
||
return {"success": False, "error": "No anchors matched pair BED IDs", "anchors_file": anchors_file, "filtered": filtered}
|
||
|
||
cmd_screen = [
|
||
conda_python, "-m", "jcvi.compara.synteny", "screen",
|
||
"--minspan=10", "--simple", str(anchors_file), str(anchors_new_file),
|
||
f"--qbed={target_bed.name}", f"--sbed={ref_bed.name}",
|
||
]
|
||
result = subprocess.run(cmd_screen, capture_output=True, text=True, cwd=str(work_dir), check=False, timeout=180)
|
||
for candidate in (work_dir / f"{pair_prefix}.anchors.new.simple", work_dir / f"{pair_prefix}.new.simple"):
|
||
if candidate.exists():
|
||
shutil.copy(candidate, simple_file)
|
||
break
|
||
if not simple_file.exists():
|
||
shutil.copy(anchors_file, simple_file)
|
||
|
||
return {
|
||
"success": simple_file.exists() and simple_file.stat().st_size > 0,
|
||
"anchors_file": anchors_file,
|
||
"anchors_new_file": anchors_new_file,
|
||
"simple_file": simple_file,
|
||
"anchor_count": anchors,
|
||
"filtered_count": filtered,
|
||
"screen_returncode": result.returncode,
|
||
"screen_stderr": result.stderr,
|
||
}
|
||
|
||
|
||
def _chr_sort_key(chrom: str):
|
||
cleaned = re.sub(r"^(chr|Chr|CHR)", "", chrom)
|
||
if cleaned.isdigit():
|
||
return (0, int(cleaned))
|
||
return (1, cleaned)
|
||
|
||
|
||
def _read_bed_chroms(bed_file: Path) -> List[str]:
|
||
skip = {"C", "M", "ChrC", "ChrM", "chrc", "chrm"}
|
||
chroms = set()
|
||
with open(bed_file, "r") as f:
|
||
for line in f:
|
||
parts = line.strip().split("\t")
|
||
if len(parts) >= 4 and parts[0] not in skip:
|
||
chroms.add(parts[0])
|
||
return sorted(chroms, key=_chr_sort_key)
|
||
|
||
|
||
def _bed_gene_ids(bed_file: Path) -> set:
|
||
ids = set()
|
||
if not bed_file.exists():
|
||
return ids
|
||
with open(bed_file, "r") as f:
|
||
for line in f:
|
||
parts = line.rstrip("\n").split("\t")
|
||
if len(parts) >= 4:
|
||
ids.add(normalize_seq_id(parts[3]))
|
||
return ids
|
||
|
||
|
||
def _is_family_gene(gene_id: str, family_ids: set) -> bool:
|
||
return bool(set(id_variants(gene_id)) & family_ids)
|
||
|
||
|
||
def _family_links_from_blast(
|
||
blast_file: Path,
|
||
family_ids: set,
|
||
target_bed_ids: set,
|
||
ref_bed_ids: set,
|
||
max_links: int = 300,
|
||
) -> List[Tuple[str, str, int]]:
|
||
if not blast_file or not blast_file.exists() or not family_ids:
|
||
return []
|
||
best: Dict[Tuple[str, str], float] = {}
|
||
with open(blast_file, "r") as f:
|
||
for line in f:
|
||
parts = line.rstrip("\n").split("\t")
|
||
if len(parts) < 12:
|
||
continue
|
||
q, s = normalize_seq_id(parts[0]), normalize_seq_id(parts[1])
|
||
if q not in target_bed_ids or s not in ref_bed_ids:
|
||
continue
|
||
if not (_is_family_gene(q, family_ids) and _is_family_gene(s, family_ids)):
|
||
continue
|
||
try:
|
||
score = float(parts[11])
|
||
except Exception:
|
||
score = 1.0
|
||
key = (q, s)
|
||
if score > best.get(key, -1):
|
||
best[key] = score
|
||
ranked = sorted(best.items(), key=lambda item: item[1], reverse=True)[:max_links]
|
||
return [(q, s, max(1, int(round(score)))) for (q, s), score in ranked]
|
||
|
||
|
||
def _write_marked_simple_file(
|
||
simple_file: Path,
|
||
marked_file: Path,
|
||
family_ids: set,
|
||
family_links: List[Tuple[str, str, int]],
|
||
) -> int:
|
||
if not simple_file.exists() or not family_ids:
|
||
return 0
|
||
marked_count = 0
|
||
seen_rows = set()
|
||
with open(simple_file, "r") as f_in, open(marked_file, "w") as f_out:
|
||
for line in f_in:
|
||
stripped = line.strip()
|
||
if not stripped or stripped.startswith("#"):
|
||
f_out.write(line)
|
||
continue
|
||
parts = stripped.split()
|
||
genes = [normalize_seq_id(p.split("*", 1)[-1]) for p in parts[:4]]
|
||
row_key = tuple(parts[:6]) if len(parts) >= 6 else tuple(parts)
|
||
seen_rows.add(row_key)
|
||
if any(_is_family_gene(g, family_ids) for g in genes) and "*" not in parts[0]:
|
||
f_out.write("#e31a1c*" + "\t".join(parts) + "\n")
|
||
marked_count += 1
|
||
else:
|
||
f_out.write(line)
|
||
for q, s, score in family_links:
|
||
row = (q, q, s, s, str(score), "+")
|
||
if row in seen_rows:
|
||
continue
|
||
f_out.write("#e31a1c*" + "\t".join(row) + "\n")
|
||
seen_rows.add(row)
|
||
marked_count += 1
|
||
return marked_count
|
||
|
||
|
||
def plot_multi_species_synteny(
|
||
work_dir: Path,
|
||
output_dir: Path,
|
||
prefix: str,
|
||
species_order: List[str],
|
||
pair_results: List[Dict[str, Any]],
|
||
family_ids: set,
|
||
) -> Dict[str, Any]:
|
||
"""Draw one multi-track JCVI karyotype figure from all successful target-ref pairs."""
|
||
print("Plotting multi-species JCVI synteny...")
|
||
conda_python = "/home/xia-lab/miniconda3/envs/test4py311/bin/python"
|
||
usable_pairs = [p for p in pair_results if p.get("simple_file") and Path(p["simple_file"]).exists()]
|
||
if not usable_pairs:
|
||
return {"success": False, "error": "No pairwise anchors.simple files available for JCVI"}
|
||
|
||
seqids_file = work_dir / f"{prefix}.seqids"
|
||
layout_file = work_dir / f"{prefix}.layout"
|
||
layout_mark_file = work_dir / f"{prefix}.layout_mark"
|
||
|
||
species_with_bed = []
|
||
with open(seqids_file, "w") as f:
|
||
for species in species_order:
|
||
bed = work_dir / f"{species}.bed"
|
||
if not bed.exists():
|
||
continue
|
||
chroms = _read_bed_chroms(bed)
|
||
if not chroms:
|
||
continue
|
||
species_with_bed.append(species)
|
||
f.write(",".join(chroms) + "\n")
|
||
|
||
if len(species_with_bed) < 2:
|
||
return {"success": False, "error": "Fewer than two species BED tracks are available"}
|
||
|
||
index = {species: i for i, species in enumerate(species_with_bed)}
|
||
n = len(species_with_bed)
|
||
y_values = [round(0.85 - (0.7 * i / max(n - 1, 1)), 3) for i in range(n)]
|
||
palette = ["#d73027", "#1a9850", "#4575b4", "#984ea3", "#ff7f00", "#4daf4a", "#377eb8", "#a65628"]
|
||
|
||
def write_layout(layout: Path, use_marked: bool) -> None:
|
||
with open(layout, "w") as f:
|
||
f.write("# y, xstart, xend, rotation, color, label, va, bed\n")
|
||
for i, species in enumerate(species_with_bed):
|
||
va = "top" if i == 0 else "bottom"
|
||
f.write(f"{y_values[i]}, .08, .92, 0, {palette[i % len(palette)]}, {species}, {va}, {species}.bed\n")
|
||
f.write("# edges\n")
|
||
for pair in usable_pairs:
|
||
target = pair["target_species"]
|
||
ref = pair["ref_species"]
|
||
if target not in index or ref not in index:
|
||
continue
|
||
edge = Path(pair.get("marked_simple_file") if use_marked and pair.get("marked_simple_file") else pair["simple_file"]).name
|
||
f.write(f"e, {index[target]}, {index[ref]}, {edge}\n")
|
||
|
||
for pair in usable_pairs:
|
||
simple = Path(pair["simple_file"])
|
||
marked = work_dir / f"{pair['target_species']}.{pair['ref_species']}.anchors.simple.mark"
|
||
target_bed_ids = _bed_gene_ids(work_dir / f"{pair['target_species']}.bed")
|
||
ref_bed_ids = _bed_gene_ids(work_dir / f"{pair['ref_species']}.bed")
|
||
family_links = _family_links_from_blast(Path(pair.get("blast_file", "")), family_ids, target_bed_ids, ref_bed_ids)
|
||
marked_count = _write_marked_simple_file(simple, marked, family_ids, family_links)
|
||
pair["family_highlight_links"] = len(family_links)
|
||
pair["marked_block_count"] = marked_count
|
||
if marked_count > 0:
|
||
pair["marked_simple_file"] = str(marked)
|
||
|
||
write_layout(layout_file, use_marked=False)
|
||
write_layout(layout_mark_file, use_marked=True)
|
||
|
||
def run_plot(layout: Path, final_pdf: Path, final_png: Path) -> Dict[str, Any]:
|
||
for stale in [work_dir / "karyotype.pdf", work_dir / "karyotype.png", final_pdf, final_png]:
|
||
if stale.exists():
|
||
stale.unlink()
|
||
cmd = [conda_python, "-m", "jcvi.graphics.karyotype", str(seqids_file.resolve()), str(layout.resolve()), "--notex"]
|
||
result = subprocess.run(cmd, capture_output=True, text=True, cwd=str(work_dir), check=False, timeout=180)
|
||
if (work_dir / "karyotype.pdf").exists():
|
||
shutil.copy(work_dir / "karyotype.pdf", final_pdf.resolve())
|
||
if (work_dir / "karyotype.png").exists():
|
||
shutil.copy(work_dir / "karyotype.png", final_png.resolve())
|
||
if final_pdf.exists() and (not final_png.exists() or final_png.stat().st_mtime < final_pdf.stat().st_mtime):
|
||
pdf_to_png(final_pdf, final_png)
|
||
return {"returncode": result.returncode, "stdout": result.stdout, "stderr": result.stderr}
|
||
|
||
final_pdf = output_dir / f"{prefix}.interspecies_synteny.pdf"
|
||
final_png = output_dir / f"{prefix}.interspecies_synteny.png"
|
||
final_mark_pdf = output_dir / f"{prefix}.interspecies_synteny.marked.pdf"
|
||
final_mark_png = output_dir / f"{prefix}.interspecies_synteny.marked.png"
|
||
|
||
base_result = run_plot(layout_file, final_pdf, final_png)
|
||
mark_result = run_plot(layout_mark_file, final_mark_pdf, final_mark_png)
|
||
|
||
return {
|
||
"success": final_pdf.exists() or final_png.exists(),
|
||
"pdf_file": final_pdf if final_pdf.exists() else None,
|
||
"png_file": final_png if final_png.exists() else None,
|
||
"marked_pdf_file": final_mark_pdf if final_mark_pdf.exists() else None,
|
||
"marked_png_file": final_mark_png if final_mark_png.exists() else None,
|
||
"seqids_file": seqids_file,
|
||
"layout_file": layout_file,
|
||
"layout_mark_file": layout_mark_file,
|
||
"base_plot": base_result,
|
||
"marked_plot": mark_result,
|
||
}
|
||
|
||
|
||
def write_report(workspace: Path, prefix: str, metrics: Dict, parameters: Dict,
|
||
tool_versions: Dict, start_time: datetime, end_time: datetime,
|
||
elapsed_sec: float, cmd: List[str]) -> Path:
|
||
"""写入报告"""
|
||
report_path = workspace / "final" / f"{prefix}.interspecies_synteny_report.md"
|
||
|
||
# 收集可复现性信息
|
||
reproducibility_info = []
|
||
reproducibility_info.append("## 可复现性信息\n\n")
|
||
reproducibility_info.append("### 执行环境\n")
|
||
reproducibility_info.append(f"- **分析时间**: {start_time.isoformat()}\n")
|
||
reproducibility_info.append(f"- **完成时间**: {end_time.isoformat()}\n")
|
||
reproducibility_info.append(f"- **运行时长**: {elapsed_sec:.2f} 秒\n")
|
||
reproducibility_info.append(f"- **工作目录**: {workspace}\n")
|
||
reproducibility_info.append(f"- **输入目录**: {parameters.get('input_dir', 'N/A')}\n\n")
|
||
|
||
reproducibility_info.append("### 工具版本\n")
|
||
for tool, version in tool_versions.items():
|
||
reproducibility_info.append(f"- **{tool}**: {version}\n")
|
||
reproducibility_info.append("\n")
|
||
|
||
reproducibility_info.append("### 执行命令\n")
|
||
reproducibility_info.append(f"```bash\n{' '.join(cmd)}\n```\n\n")
|
||
|
||
reproducibility_info.append("### 输入参数\n")
|
||
for key, value in parameters.items():
|
||
if key != 'input_dir':
|
||
reproducibility_info.append(f"- **{key}**: {value}\n")
|
||
reproducibility_info.append("\n")
|
||
|
||
with open(report_path, "w") as f:
|
||
f.write(f"# 物种间共线性分析 - {prefix}\n\n")
|
||
f.write(f"**时间**: {datetime.now().isoformat()}\n\n")
|
||
f.write(f"- 目标物种 (研究对象): {parameters.get('target_species')}\n")
|
||
f.write(f"- 参考物种 (对照): {parameters.get('ref_species')}\n")
|
||
f.write(f"- 目标蛋白文件: {parameters.get('target_protein_spl')}\n")
|
||
f.write(f"- 参考蛋白文件: {parameters.get('ref_protein_spl')}\n")
|
||
f.write(f"- 目标全蛋白组: {parameters.get('target_protein_all')}\n")
|
||
f.write(f"- 参考全蛋白组: {parameters.get('ref_protein_all')}\n")
|
||
f.write("\n")
|
||
f.writelines(reproducibility_info)
|
||
return report_path
|
||
|
||
|
||
def get_file_info(file_path: Path) -> Dict[str, Any]:
|
||
"""获取文件的详细信息，用于artifact规范"""
|
||
info = {
|
||
"size": 0,
|
||
"sha256": "",
|
||
"mtime": 0
|
||
}
|
||
if file_path.exists():
|
||
info["size"] = file_path.stat().st_size
|
||
info["mtime"] = file_path.stat().st_mtime
|
||
info["sha256"] = compute_file_hash(file_path)
|
||
return info
|
||
|
||
|
||
def write_manifest(workspace: Path, job_id: str, run_id: str, prefix: str,
|
||
start_time: datetime, end_time: datetime, artifacts: Dict,
|
||
metrics: Dict, parameters: Dict, tool_versions: Dict,
|
||
workspace_root: Path, input_dir: Path, cmd: List[str],
|
||
status: str = "success") -> Path:
|
||
"""写入manifest - 符合schema v1标准"""
|
||
elapsed_sec = (end_time - start_time).total_seconds()
|
||
|
||
# 构建符合规范的artifacts结构
|
||
formatted_artifacts = {}
|
||
for key, value in artifacts.items():
|
||
if isinstance(value, dict):
|
||
# 已经是字典格式，确保包含必需字段
|
||
formatted_artifacts[key] = value
|
||
elif isinstance(value, str):
|
||
# 字符串格式，转换为artifact对象
|
||
file_path = workspace / value
|
||
file_info = get_file_info(file_path)
|
||
formatted_artifacts[key] = {
|
||
"artifact": key,
|
||
"path": value,
|
||
"full_path": str(file_path.absolute()),
|
||
"type": "file",
|
||
"size": file_info["size"],
|
||
"sha256": file_info["sha256"],
|
||
"mtime": file_info["mtime"]
|
||
}
|
||
|
||
manifest_data = {
|
||
"manifest_version": "1.0",
|
||
"job_id": job_id,
|
||
"stage": "interspecies_synteny",
|
||
"run_id": run_id,
|
||
"prefix": prefix,
|
||
"timestamp": end_time.isoformat(),
|
||
"status": status,
|
||
"start_time": start_time.isoformat(),
|
||
"end_time": end_time.isoformat(),
|
||
"elapsed_sec": round(elapsed_sec, 2),
|
||
"cmd": cmd,
|
||
"workspace_root": str(workspace_root.absolute()),
|
||
"workspace": str(workspace.absolute()),
|
||
"input_dir": str(input_dir.absolute()),
|
||
"tool_versions": tool_versions,
|
||
"artifacts": formatted_artifacts,
|
||
"parameters": parameters,
|
||
"metrics": metrics
|
||
}
|
||
|
||
manifest_path = workspace / "final" / f"{prefix}.interspecies_synteny.manifest.json"
|
||
with open(manifest_path, "w", encoding='utf-8') as f:
|
||
json.dump(manifest_data, f, indent=2, ensure_ascii=False)
|
||
|
||
return manifest_path
|
||
|
||
|
||
|
||
def _clean_bool(value: Any) -> bool:
|
||
return str(value).strip().lower() in {"true", "1", "yes", "y"}
|
||
|
||
|
||
def _gene_family_from_prefix(prefix: str) -> str:
|
||
return prefix.split('.', 1)[1] if '.' in prefix else prefix
|
||
|
||
|
||
def _resolve_path(path_text: str, input_dir: Path) -> Optional[Path]:
|
||
if not path_text:
|
||
return None
|
||
path = Path(path_text)
|
||
if path.is_absolute():
|
||
return path if path.exists() else None
|
||
for candidate in (input_dir / path, input_dir / path.name):
|
||
if candidate.exists():
|
||
return candidate
|
||
return None
|
||
|
||
|
||
def _read_config_groups(input_dir: Path, default_family: str) -> Dict[str, Any]:
|
||
config_path = input_dir / "config.xlsx"
|
||
groups = {"target": {}, "ref": {}, "gene_family": default_family}
|
||
if not config_path.exists():
|
||
return groups
|
||
try:
|
||
import pandas as pd
|
||
df = pd.read_excel(config_path)
|
||
except Exception as exc:
|
||
print(f"[WARN] Failed to read config.xlsx: {exc}")
|
||
return groups
|
||
for _, row in df.iterrows():
|
||
species = str(row.get("Species_Code", "")).strip()
|
||
family = str(row.get("Gene_Family", default_family)).strip() or default_family
|
||
file_type = str(row.get("File_Type", "")).strip().lower()
|
||
original_name = str(row.get("Original_File_Name", "")).strip()
|
||
if not species:
|
||
continue
|
||
group_name = "ref" if _clean_bool(row.get("Ref", False)) else "target"
|
||
item = groups[group_name].setdefault(species, {"species": species, "gene_family": family, "files": {}})
|
||
item["gene_family"] = family
|
||
if file_type:
|
||
item["files"][file_type] = original_name
|
||
groups["gene_family"] = family or groups["gene_family"]
|
||
return groups
|
||
|
||
|
||
def _read_genome_prepare_artifacts(workspace: Path, input_dir: Path, prefix: str) -> Dict[str, Dict[str, Path]]:
|
||
result: Dict[str, Dict[str, Path]] = {}
|
||
for manifest_path in (workspace / "final" / f"{prefix}.genome_prepare.manifest.json", input_dir / "final" / f"{prefix}.genome_prepare.manifest.json"):
|
||
if not manifest_path.exists():
|
||
continue
|
||
try:
|
||
data = json.loads(manifest_path.read_text(encoding="utf-8"))
|
||
except Exception as exc:
|
||
print(f"[WARN] Failed to read genome_prepare manifest {manifest_path}: {exc}")
|
||
continue
|
||
for key, artifact in (data.get("artifacts") or {}).items():
|
||
kind = artifact.get("artifact")
|
||
species = key.rsplit("_", 1)[0]
|
||
found = _resolve_path(artifact.get("full_path") or artifact.get("path") or "", input_dir)
|
||
if found and kind:
|
||
result.setdefault(species, {})[kind] = found
|
||
return result
|
||
|
||
|
||
def _find_species_file(input_dir: Path, species: str, family: str, file_type: str, original_name: str = "") -> Optional[Path]:
|
||
if file_type == "all_pep":
|
||
patterns = [
|
||
f"{species}.{family}.all_pep.fasta", f"{species}.{family}.all_pep.fa",
|
||
f"{species}.{family}.protein.fasta", f"{species}.{family}.faa",
|
||
f"{species}.{family}.pep.fasta", f"{species}.{family}.pep.fa",
|
||
f"{species}.all_pep.fasta", f"{species}.pep.fasta", f"{species}.protein.fasta", f"{species}.faa",
|
||
]
|
||
elif file_type == "family_pep":
|
||
patterns = [
|
||
f"{species}.{family}.family_pep.fasta", f"{species}.{family}.family_pep.fa",
|
||
f"{species}.{family}.pep.fasta", f"{species}.{family}.pep.fa",
|
||
]
|
||
elif file_type == "pep":
|
||
patterns = [
|
||
f"{species}.{family}.pep.fasta", f"{species}.{family}.pep.fa",
|
||
f"{species}.{family}.protein.fasta", f"{species}.{family}.faa",
|
||
f"{species}.pep.fasta", f"{species}.pep.fa", f"{species}.protein.fasta", f"{species}.faa",
|
||
]
|
||
elif file_type == "gff":
|
||
patterns = [
|
||
f"{species}.{family}_final.gff3", f"{species}.{family}.gff3", f"{species}.{family}.gff",
|
||
f"{species}_final.gff3", f"{species}.gff3", f"{species}.gff",
|
||
]
|
||
else:
|
||
patterns = []
|
||
if original_name:
|
||
patterns.append(Path(original_name).name)
|
||
for pattern in patterns:
|
||
for candidate in (input_dir / pattern, input_dir / species / pattern):
|
||
if candidate.exists():
|
||
return candidate
|
||
return None
|
||
|
||
def _parse_species_selection(value: Optional[str]) -> List[str]:
|
||
if not value:
|
||
return []
|
||
return [item.strip() for item in str(value).split(',') if item.strip()]
|
||
|
||
|
||
def _filter_group_by_selection(group: Dict[str, Any], selected: List[str]) -> Dict[str, Any]:
|
||
if not selected:
|
||
return group
|
||
selected_set = set(selected)
|
||
return {species: item for species, item in group.items() if species in selected_set}
|
||
|
||
|
||
def _species_default_selection(available: List[str], preferred_ref: bool = False) -> List[str]:
|
||
if not available:
|
||
return []
|
||
if preferred_ref:
|
||
for preferred in ("at", "ath"):
|
||
for species in available:
|
||
if species.lower() == preferred:
|
||
return [species]
|
||
return [available[0]]
|
||
return list(available)
|
||
|
||
|
||
def _gene_matches_species(gene_id: str, species: str, species_id_sets: Dict[str, set]) -> bool:
|
||
variants = id_variants(gene_id)
|
||
if species in species_id_sets and variants & species_id_sets[species]:
|
||
return True
|
||
return normalize_seq_id(gene_id).lower().startswith(species.lower())
|
||
|
||
|
||
def _filter_fasta_by_species(source: Path, output_file: Path, selected_species: List[str], species_id_sets: Dict[str, set]) -> Path:
|
||
output_file.parent.mkdir(parents=True, exist_ok=True)
|
||
if not selected_species:
|
||
shutil.copy(source, output_file)
|
||
return output_file
|
||
keep = False
|
||
written = 0
|
||
with open(source, "r") as f_in, open(output_file, "w") as f_out:
|
||
for line in f_in:
|
||
if line.startswith(">"):
|
||
gene_id = normalize_seq_id(line[1:])
|
||
keep = any(_gene_matches_species(gene_id, species, species_id_sets) for species in selected_species)
|
||
if keep:
|
||
written += 1
|
||
f_out.write(line)
|
||
elif keep:
|
||
f_out.write(line)
|
||
if written == 0:
|
||
raise RuntimeError(f"Selected species {selected_species} produced an empty filtered FASTA from {source}")
|
||
return output_file
|
||
|
||
|
||
def discover_group_inputs(workspace: Path, input_dir: Path, prefix: str, args) -> Dict[str, Any]:
|
||
family = _gene_family_from_prefix(prefix)
|
||
groups = _read_config_groups(input_dir, family)
|
||
artifacts = _read_genome_prepare_artifacts(workspace, input_dir, prefix)
|
||
family = groups.get("gene_family") or family
|
||
|
||
available_target_species = list(groups["target"].keys())
|
||
available_ref_species = list(groups["ref"].keys())
|
||
requested_targets = _parse_species_selection(getattr(args, "target_species", None))
|
||
requested_refs = _parse_species_selection(getattr(args, "ref_species", None))
|
||
selected_targets = requested_targets or _species_default_selection(available_target_species, preferred_ref=False)
|
||
selected_refs = requested_refs or _species_default_selection(available_ref_species, preferred_ref=True)
|
||
|
||
unknown_targets = [sp for sp in selected_targets if sp not in groups["target"]]
|
||
unknown_refs = [sp for sp in selected_refs if sp not in groups["ref"]]
|
||
if unknown_targets:
|
||
raise RuntimeError(f"Selected target species not found in config Ref=False group: {', '.join(unknown_targets)}")
|
||
if unknown_refs:
|
||
raise RuntimeError(f"Selected reference species not found in config Ref=True group: {', '.join(unknown_refs)}")
|
||
if not selected_targets:
|
||
raise RuntimeError("No target species selected for interspecies analysis")
|
||
if not selected_refs:
|
||
raise RuntimeError("No reference species selected for interspecies analysis")
|
||
|
||
groups["target"] = _filter_group_by_selection(groups["target"], selected_targets)
|
||
groups["ref"] = _filter_group_by_selection(groups["ref"], selected_refs)
|
||
|
||
for group_name in ("target", "ref"):
|
||
for species, item in groups[group_name].items():
|
||
files = item.setdefault("resolved", {})
|
||
artifact_paths = artifacts.get(species, {})
|
||
original_files = item.get("files", {})
|
||
item_family = item.get("gene_family", family)
|
||
all_original = original_files.get("all_pep", "") or original_files.get("protein", "")
|
||
family_original = original_files.get("family_pep", "")
|
||
legacy_pep_original = original_files.get("pep", "")
|
||
if group_name == "ref" and legacy_pep_original and not family_original:
|
||
family_original = legacy_pep_original
|
||
if group_name == "target" and legacy_pep_original and not all_original:
|
||
all_original = legacy_pep_original
|
||
files["all_pep"] = (
|
||
_find_species_file(input_dir, species, item_family, "all_pep", all_original) or
|
||
artifact_paths.get("all_pep") or artifact_paths.get("ref_all_pep") or artifact_paths.get("pep")
|
||
)
|
||
files["family_pep"] = (
|
||
_find_species_file(input_dir, species, item_family, "family_pep", family_original) or
|
||
artifact_paths.get("family_pep") or artifact_paths.get("ref_family_pep") or artifact_paths.get("ref_pep")
|
||
)
|
||
files["pep"] = files.get("all_pep") or files.get("family_pep")
|
||
files["gff"] = (
|
||
artifact_paths.get("gff") or artifact_paths.get("ref_gff") or
|
||
_find_species_file(input_dir, species, item_family, "gff", original_files.get("gff", ""))
|
||
)
|
||
|
||
analysis_dir = workspace / "interspecies_synteny" / "analysis"
|
||
analysis_dir.mkdir(parents=True, exist_ok=True)
|
||
species_id_sets = collect_species_id_sets(groups)
|
||
|
||
raw_target_family_protein = Path(args.target_protein) if args.target_protein else workspace / "final" / f"{prefix}.final_members.pep.fa"
|
||
target_family_protein = raw_target_family_protein
|
||
if raw_target_family_protein and raw_target_family_protein.exists():
|
||
target_family_protein = _filter_fasta_by_species(
|
||
raw_target_family_protein,
|
||
analysis_dir / f"{prefix}.selected_target.family.pep.fa",
|
||
selected_targets,
|
||
species_id_sets,
|
||
)
|
||
|
||
ref_family_protein = Path(args.ref_protein) if args.ref_protein else None
|
||
if ref_family_protein is None:
|
||
ref_peps = [item.get("resolved", {}).get("family_pep") for item in groups["ref"].values() if item.get("resolved", {}).get("family_pep")]
|
||
if ref_peps:
|
||
ref_family_protein = analysis_dir / f"Ref_selected.{family}.pep.fasta"
|
||
merge_many_fastas(ref_peps, ref_family_protein)
|
||
else:
|
||
for candidate in [input_dir / f"Ref_combined.{family}.pep.fasta", input_dir / f"Ref_combined.{family}.pep.fa", input_dir / "Ref_combined.pep.fasta"]:
|
||
if candidate.exists():
|
||
ref_family_protein = candidate
|
||
break
|
||
if ref_family_protein and Path(ref_family_protein).exists():
|
||
ref_family_protein = _filter_fasta_by_species(
|
||
Path(ref_family_protein),
|
||
analysis_dir / f"{prefix}.selected_ref.family.pep.fa",
|
||
selected_refs,
|
||
species_id_sets,
|
||
)
|
||
|
||
target_all_peps = {sp: item.get("resolved", {}).get("all_pep") for sp, item in groups["target"].items() if item.get("resolved", {}).get("all_pep")}
|
||
ref_all_peps = {sp: item.get("resolved", {}).get("all_pep") for sp, item in groups["ref"].items() if item.get("resolved", {}).get("all_pep")}
|
||
target_gffs = {sp: item.get("resolved", {}).get("gff") for sp, item in groups["target"].items() if item.get("resolved", {}).get("gff")}
|
||
ref_gffs = {sp: item.get("resolved", {}).get("gff") for sp, item in groups["ref"].items() if item.get("resolved", {}).get("gff")}
|
||
|
||
return {
|
||
"gene_family": family,
|
||
"target_species": selected_targets,
|
||
"ref_species": selected_refs,
|
||
"available_target_species": available_target_species,
|
||
"available_ref_species": available_ref_species,
|
||
"selected_target_species": selected_targets,
|
||
"selected_ref_species": selected_refs,
|
||
"target_family_protein": target_family_protein,
|
||
"ref_family_protein": ref_family_protein,
|
||
"target_all_peps": target_all_peps,
|
||
"ref_all_peps": ref_all_peps,
|
||
"target_gffs": target_gffs,
|
||
"ref_gffs": ref_gffs,
|
||
"groups": groups,
|
||
}
|
||
|
||
|
||
def merge_many_fastas(files: List[Path], output_file: Path) -> Path:
|
||
output_file.parent.mkdir(parents=True, exist_ok=True)
|
||
with open(output_file, "w") as out:
|
||
for fasta in files:
|
||
with open(fasta, "r") as f:
|
||
out.write(f.read().rstrip() + "\n")
|
||
return output_file
|
||
|
||
|
||
def _parse_gff_attr(attr: str) -> Dict[str, str]:
|
||
result = {}
|
||
for item in attr.strip().strip(';').split(';'):
|
||
item = item.strip()
|
||
if not item:
|
||
continue
|
||
if '=' in item:
|
||
k, v = item.split('=', 1)
|
||
elif ' ' in item:
|
||
k, v = item.split(' ', 1)
|
||
else:
|
||
continue
|
||
result[k.strip()] = v.strip().strip('"')
|
||
return result
|
||
|
||
|
||
def _match_gff_id(attrs: Dict[str, str], protein_ids: set) -> Optional[str]:
|
||
candidates = [attrs.get("ID"), attrs.get("transcript_id"), attrs.get("Name"), attrs.get("Parent")]
|
||
for candidate in candidates:
|
||
if not candidate:
|
||
continue
|
||
for item in str(candidate).split(','):
|
||
norm = normalize_seq_id(item)
|
||
if norm in protein_ids:
|
||
return norm
|
||
for variant in id_variants(norm):
|
||
if variant in protein_ids:
|
||
return variant
|
||
return None
|
||
|
||
|
||
def build_species_position_files(species: str, gff3: Path, protein_fasta: Path, work_dir: Path) -> Dict[str, Any]:
|
||
"""Create MCScanX simplified GFF and JCVI BED for one species."""
|
||
simple_gff = work_dir / f"{species}.gff"
|
||
bed_file = work_dir / f"{species}.bed"
|
||
protein_ids = fasta_id_set(protein_fasta)
|
||
seen = set()
|
||
count = 0
|
||
features_seen = 0
|
||
work_dir.mkdir(parents=True, exist_ok=True)
|
||
with open(gff3, "r") as f_in, open(simple_gff, "w") as gff_out, open(bed_file, "w") as bed_out:
|
||
for line in f_in:
|
||
if not line.strip() or line.startswith('#'):
|
||
continue
|
||
parts = line.rstrip('\n').split(' ')
|
||
if len(parts) < 9 or parts[2] not in {"mRNA", "transcript", "gene"}:
|
||
continue
|
||
features_seen += 1
|
||
attrs = _parse_gff_attr(parts[8])
|
||
gene_id = _match_gff_id(attrs, protein_ids)
|
||
if not gene_id or gene_id in seen:
|
||
continue
|
||
seen.add(gene_id)
|
||
chrom = parts[0]
|
||
start, end = parts[3], parts[4]
|
||
gff_out.write(f"{chrom} {gene_id} {start} {end}\n")
|
||
bed_out.write(f"{chrom} {start} {end} {gene_id}\n")
|
||
count += 1
|
||
return {
|
||
"success": count > 0,
|
||
"species": species,
|
||
"gff": simple_gff,
|
||
"bed": bed_file,
|
||
"gene_count": count,
|
||
"features_seen": features_seen,
|
||
"protein_id_count": len(protein_ids),
|
||
"error": "No GFF IDs matched protein FASTA IDs" if count == 0 else "",
|
||
}
|
||
|
||
|
||
def collect_family_ids(target_family_fasta: Path, ref_family_fasta: Path, final_id_file: Optional[Path] = None) -> set:
|
||
ids = set()
|
||
for fasta in [target_family_fasta, ref_family_fasta]:
|
||
if fasta and Path(fasta).exists():
|
||
ids.update(fasta_id_set(Path(fasta)))
|
||
if final_id_file and final_id_file.exists():
|
||
with open(final_id_file, "r") as f:
|
||
for line in f:
|
||
if line.strip() and not line.startswith('#'):
|
||
ids.update(id_variants(line.strip().split()[0]))
|
||
return ids
|
||
|
||
|
||
def add_artifact(artifacts: Dict[str, Dict], key: str, file_path: Path, rel_path: str, artifact: str, file_type: str) -> None:
|
||
if not file_path or not file_path.exists():
|
||
return
|
||
info = get_file_info(file_path)
|
||
artifacts[key] = {
|
||
"artifact": artifact,
|
||
"path": rel_path,
|
||
"full_path": str(file_path.absolute()),
|
||
"type": file_type,
|
||
"size": info["size"],
|
||
"sha256": info["sha256"],
|
||
"mtime": info["mtime"],
|
||
}
|
||
|
||
|
||
def write_failure_manifest(workspace: Path, job_id: str, run_id: str, prefix: str, start_time: datetime,
|
||
workspace_root: Path, input_dir: Path, cmd: List[str], error: str) -> Path:
|
||
end_time = datetime.now()
|
||
return write_manifest(
|
||
workspace, job_id, run_id, prefix, start_time, end_time,
|
||
artifacts={}, metrics={"error": error}, parameters={"prefix": prefix, "error": error},
|
||
tool_versions={}, workspace_root=workspace_root, input_dir=input_dir, cmd=cmd, status="failed"
|
||
)
|
||
|
||
|
||
def main():
|
||
parser = argparse.ArgumentParser(description="????????")
|
||
parser.add_argument("--prefix", required=True, help="?????? Multi.SPL")
|
||
parser.add_argument("--workspace", required=True)
|
||
parser.add_argument("--input-dir")
|
||
parser.add_argument("--ref-species", default=None, help="???????? config.xlsx ??")
|
||
parser.add_argument("--target-species", default=None, help="???????? config.xlsx ??")
|
||
parser.add_argument("--ref-protein", help="?????????")
|
||
parser.add_argument("--target-protein", help="?????????")
|
||
parser.add_argument("--cpu", type=int, default=8)
|
||
parser.add_argument("--log-file", help="??????????????")
|
||
parser.add_argument("--workspace-root", help="??????????manifest?")
|
||
args = parser.parse_args()
|
||
|
||
start_time = datetime.now()
|
||
run_id = f"interspecies_synteny_{start_time.strftime('%Y%m%d_%H%M%S')}"
|
||
workspace = Path(args.workspace)
|
||
job_id = workspace.name
|
||
input_dir = Path(args.input_dir) if args.input_dir else workspace.parent.parent / "input" / job_id
|
||
workspace_root = Path(args.workspace_root) if args.workspace_root else workspace.parent
|
||
cmd = sys.argv
|
||
|
||
print("=" * 60)
|
||
print("Interspecies synteny analysis")
|
||
print("=" * 60)
|
||
|
||
try:
|
||
paths = prepare_directories(workspace)
|
||
tool_versions = get_tool_versions()
|
||
discovery = discover_group_inputs(workspace, input_dir, args.prefix, args)
|
||
target_group = "Target_combined"
|
||
ref_group = "Ref_combined"
|
||
print(f"Target species: {', '.join(discovery['target_species']) or 'NONE'}")
|
||
print(f"Reference species: {', '.join(discovery['ref_species']) or 'NONE'}")
|
||
print(f"Target family protein: {discovery['target_family_protein']}")
|
||
print(f"Ref family protein: {discovery['ref_family_protein']}")
|
||
|
||
if not discovery["target_family_protein"] or not Path(discovery["target_family_protein"]).exists():
|
||
raise RuntimeError(f"Target family protein missing: {discovery['target_family_protein']}")
|
||
if not discovery["ref_family_protein"] or not Path(discovery["ref_family_protein"]).exists():
|
||
raise RuntimeError("Reference family protein missing; expected Ref_combined.<Family>.pep.fasta or Ref species pep files")
|
||
|
||
artifacts = {}
|
||
metrics = {
|
||
"target_species": discovery["target_species"],
|
||
"ref_species": discovery["ref_species"],
|
||
"available_target_species": discovery.get("available_target_species", discovery["target_species"]),
|
||
"available_ref_species": discovery.get("available_ref_species", discovery["ref_species"]),
|
||
"selected_target_species": discovery.get("selected_target_species", discovery["target_species"]),
|
||
"selected_ref_species": discovery.get("selected_ref_species", discovery["ref_species"]),
|
||
"tree_success": False,
|
||
"synteny_success": False,
|
||
"synteny_skipped_reason": "",
|
||
}
|
||
|
||
print("\nMerging family proteins for tree...")
|
||
merged_fasta = paths["analysis_dir"] / f"{args.prefix}.target_ref.family.pep.fa"
|
||
merge_result = merge_protein_sequences(Path(discovery["target_family_protein"]), Path(discovery["ref_family_protein"]), merged_fasta)
|
||
if not merge_result.get("success"):
|
||
raise RuntimeError(f"Failed to merge family proteins: {merge_result.get('stderr')}")
|
||
metrics["target_family_count"] = merge_result.get("target_count")
|
||
metrics["ref_family_count"] = merge_result.get("ref_count")
|
||
|
||
print("\nBuilding interspecies tree...")
|
||
tree_result = build_phylogenetic_tree(merged_fasta, paths["analysis_dir"], args.prefix, args.cpu)
|
||
if not tree_result.get("tree_file"):
|
||
raise RuntimeError(f"Tree construction failed: {tree_result.get('error', 'unknown error')}")
|
||
final_tree = paths["final_dir"] / f"{args.prefix}.interspecies_tree.nwk"
|
||
shutil.copy(Path(tree_result["tree_file"]), final_tree)
|
||
metrics["tree_success"] = True
|
||
add_artifact(artifacts, "tree_nwk", final_tree, f"final/{args.prefix}.interspecies_tree.nwk", "tree_nwk", "tree")
|
||
|
||
gene_species_map = paths["analysis_dir"] / f"{args.prefix}.gene_to_species.tsv"
|
||
write_gene_species_map(Path(discovery["target_family_protein"]), Path(discovery["ref_family_protein"]), discovery["target_species"], discovery["ref_species"], gene_species_map, collect_species_id_sets(discovery["groups"]))
|
||
add_artifact(artifacts, "gene_to_species", gene_species_map, f"interspecies_synteny/analysis/{args.prefix}.gene_to_species.tsv", "gene_to_species", "tsv")
|
||
|
||
print("\nPlotting interspecies tree...")
|
||
target_subfamily_file = input_dir / f"{args.prefix}.geneID.subfamily"
|
||
if not target_subfamily_file.exists():
|
||
target_subfamily_file = workspace / "final" / f"{args.prefix}.geneID.subfamily"
|
||
ref_subfamily_file = input_dir / f"SPL_{discovery['ref_species'][0]}.subfamily" if discovery["ref_species"] else input_dir / "SPL_ref.subfamily"
|
||
if not ref_subfamily_file.exists():
|
||
candidates = sorted(input_dir.glob("SPL_*.subfamily")) + sorted(input_dir.glob("*.SPL.subfamily"))
|
||
ref_subfamily_file = candidates[0] if candidates else ref_subfamily_file
|
||
tree_plot_result = plot_interspecies_tree(
|
||
final_tree,
|
||
paths["final_dir"],
|
||
args.prefix,
|
||
target_group,
|
||
ref_group,
|
||
target_subfamily_file if target_subfamily_file.exists() else None,
|
||
ref_subfamily_file if ref_subfamily_file.exists() else None,
|
||
gene_species_map,
|
||
)
|
||
if tree_plot_result.get("pdf_file"):
|
||
add_artifact(artifacts, "tree_pdf", tree_plot_result["pdf_file"], f"final/{args.prefix}.interspecies_tree.pdf", "tree_pdf", "plot")
|
||
if tree_plot_result.get("png_file"):
|
||
add_artifact(artifacts, "tree_png", tree_plot_result["png_file"], f"final/{args.prefix}.interspecies_tree.png", "tree_png", "plot")
|
||
|
||
synteny_plot_result = {"success": False}
|
||
pair_results: List[Dict[str, Any]] = []
|
||
species_position_results: Dict[str, Dict[str, Any]] = {}
|
||
|
||
missing = []
|
||
if not discovery["target_all_peps"]:
|
||
missing.append("target all protein")
|
||
if not discovery["ref_all_peps"]:
|
||
missing.append("reference all protein")
|
||
if not discovery["target_gffs"]:
|
||
missing.append("target GFF")
|
||
if not discovery["ref_gffs"]:
|
||
missing.append("reference GFF")
|
||
|
||
if not missing:
|
||
print("\nPreparing per-species GFF/BED files for pairwise MCScanX...")
|
||
all_species = discovery["target_species"] + discovery["ref_species"]
|
||
for species in all_species:
|
||
pep = discovery["target_all_peps"].get(species) or discovery["ref_all_peps"].get(species)
|
||
gff = discovery["target_gffs"].get(species) or discovery["ref_gffs"].get(species)
|
||
if not pep or not gff:
|
||
species_position_results[species] = {"success": False, "error": "missing pep or gff"}
|
||
continue
|
||
result = build_species_position_files(species, Path(gff), Path(pep), paths["analysis_dir"])
|
||
species_position_results[species] = result
|
||
if result.get("success"):
|
||
add_artifact(artifacts, f"{species}_mcscanx_gff", result["gff"], f"interspecies_synteny/analysis/{species}.gff", "mcscanx_gff", "gff")
|
||
add_artifact(artifacts, f"{species}_bed", result["bed"], f"interspecies_synteny/analysis/{species}.bed", "jcvi_bed", "bed")
|
||
print(f" {species}: {result.get('gene_count', 0)} matched genes")
|
||
|
||
print("\nRunning pairwise target-ref MCScanX...")
|
||
for target_species in discovery["target_species"]:
|
||
for ref_species in discovery["ref_species"]:
|
||
pair_key = f"{target_species}.{ref_species}"
|
||
pair = {"pair": pair_key, "target_species": target_species, "ref_species": ref_species, "success": False}
|
||
target_ready = species_position_results.get(target_species, {}).get("success")
|
||
ref_ready = species_position_results.get(ref_species, {}).get("success")
|
||
target_pep = discovery["target_all_peps"].get(target_species)
|
||
ref_pep = discovery["ref_all_peps"].get(ref_species)
|
||
if not (target_ready and ref_ready and target_pep and ref_pep):
|
||
pair["error"] = "missing matched GFF/BED or protein FASTA"
|
||
pair_results.append(pair)
|
||
print(f" {pair_key}: skipped ({pair['error']})")
|
||
continue
|
||
blast_result = run_reciprocal_blast(Path(target_pep), Path(ref_pep), paths["analysis_dir"], target_species, ref_species, args.cpu)
|
||
pair["blast_success"] = bool(blast_result.get("success"))
|
||
if blast_result.get("output_file"):
|
||
pair["blast_file"] = str(blast_result["output_file"])
|
||
add_artifact(artifacts, f"{pair_key}_blast", blast_result["output_file"], f"interspecies_synteny/analysis/{pair_key}.blast", "blast", "data")
|
||
if not blast_result.get("success"):
|
||
pair["error"] = "BLAST failed or produced no hits"
|
||
pair_results.append(pair)
|
||
print(f" {pair_key}: {pair['error']}")
|
||
continue
|
||
mcscanx_result = run_mcscanx_interspecies(paths["analysis_dir"], target_species, ref_species)
|
||
pair_gff = paths["analysis_dir"] / f"{pair_key}.gff"
|
||
if pair_gff.exists():
|
||
pair["pair_gff"] = str(pair_gff)
|
||
add_artifact(artifacts, f"{pair_key}_gff", pair_gff, f"interspecies_synteny/analysis/{pair_key}.gff", "pair_mcscanx_gff", "gff")
|
||
pair["mcscanx_success"] = bool(mcscanx_result.get("success"))
|
||
if mcscanx_result.get("collinearity_file"):
|
||
pair["collinearity_file"] = str(mcscanx_result["collinearity_file"])
|
||
add_artifact(artifacts, f"{pair_key}_collinearity", mcscanx_result["collinearity_file"], f"interspecies_synteny/analysis/{pair_key}.collinearity", "collinearity", "data")
|
||
if not mcscanx_result.get("success"):
|
||
pair["error"] = "MCScanX did not produce collinearity"
|
||
pair_results.append(pair)
|
||
print(f" {pair_key}: {pair['error']}")
|
||
continue
|
||
simple_result = convert_collinearity_to_simple(paths["analysis_dir"], Path(mcscanx_result["collinearity_file"]), target_species, ref_species)
|
||
if not simple_result.get("success") and blast_result.get("output_file"):
|
||
fallback_result = convert_blast_to_family_simple(paths["analysis_dir"], Path(blast_result["output_file"]), target_species, ref_species)
|
||
if fallback_result.get("success"):
|
||
fallback_result["note"] = "MCScanX produced no whole-genome blocks; using family-only BLAST links for JCVI visualization"
|
||
simple_result = fallback_result
|
||
pair.update({k: str(v) if isinstance(v, Path) else v for k, v in simple_result.items()})
|
||
pair["success"] = bool(simple_result.get("success"))
|
||
if simple_result.get("anchors_file"):
|
||
add_artifact(artifacts, f"{pair_key}_anchors", simple_result["anchors_file"], f"interspecies_synteny/analysis/{pair_key}.anchors", "anchors", "data")
|
||
if simple_result.get("simple_file"):
|
||
add_artifact(artifacts, f"{pair_key}_anchors_simple", simple_result["simple_file"], f"interspecies_synteny/analysis/{pair_key}.anchors.simple", "anchors_simple", "data")
|
||
pair_results.append(pair)
|
||
print(f" {pair_key}: {'success' if pair['success'] else 'anchors conversion failed'}")
|
||
|
||
family_ids = collect_family_ids(Path(discovery["target_family_protein"]), Path(discovery["ref_family_protein"]), workspace / "final" / f"{args.prefix}.final_members.id")
|
||
synteny_plot_result = plot_multi_species_synteny(
|
||
paths["analysis_dir"], paths["final_dir"], args.prefix,
|
||
discovery["target_species"] + discovery["ref_species"], pair_results, family_ids,
|
||
)
|
||
successful_pairs = [p for p in pair_results if p.get("success")]
|
||
metrics["synteny_success"] = bool(synteny_plot_result.get("success")) and len(successful_pairs) == len(pair_results) and bool(pair_results)
|
||
if not successful_pairs:
|
||
metrics["synteny_skipped_reason"] = synteny_plot_result.get("error") or "No target-ref pair produced JCVI anchors"
|
||
else:
|
||
metrics["synteny_skipped_reason"] = "Missing full data for MCScanX: " + ", ".join(missing)
|
||
print(metrics["synteny_skipped_reason"])
|
||
|
||
metrics["species_position_results"] = {
|
||
sp: {k: str(v) if isinstance(v, Path) else v for k, v in result.items()}
|
||
for sp, result in species_position_results.items()
|
||
}
|
||
metrics["pairwise_synteny"] = pair_results
|
||
metrics["pair_success_count"] = sum(1 for p in pair_results if p.get("success"))
|
||
metrics["pair_total_count"] = len(pair_results)
|
||
|
||
if synteny_plot_result.get("pdf_file"):
|
||
add_artifact(artifacts, "synteny_pdf", synteny_plot_result["pdf_file"], f"final/{args.prefix}.interspecies_synteny.pdf", "synteny_pdf", "plot")
|
||
if synteny_plot_result.get("png_file"):
|
||
add_artifact(artifacts, "synteny_png", synteny_plot_result["png_file"], f"final/{args.prefix}.interspecies_synteny.png", "synteny_png", "plot")
|
||
if synteny_plot_result.get("marked_pdf_file"):
|
||
add_artifact(artifacts, "synteny_marked_pdf", synteny_plot_result["marked_pdf_file"], f"final/{args.prefix}.interspecies_synteny.marked.pdf", "synteny_marked_pdf", "plot")
|
||
if synteny_plot_result.get("marked_png_file"):
|
||
add_artifact(artifacts, "synteny_marked_png", synteny_plot_result["marked_png_file"], f"final/{args.prefix}.interspecies_synteny.marked.png", "synteny_marked_png", "plot")
|
||
if synteny_plot_result.get("seqids_file"):
|
||
add_artifact(artifacts, "jcvi_seqids", synteny_plot_result["seqids_file"], f"interspecies_synteny/analysis/{args.prefix}.seqids", "jcvi_seqids", "data")
|
||
if synteny_plot_result.get("layout_file"):
|
||
add_artifact(artifacts, "jcvi_layout", synteny_plot_result["layout_file"], f"interspecies_synteny/analysis/{args.prefix}.layout", "jcvi_layout", "data")
|
||
|
||
end_time = datetime.now()
|
||
status = "success" if metrics["tree_success"] and metrics["synteny_success"] else "partial_success"
|
||
parameters = {
|
||
"prefix": args.prefix,
|
||
"gene_family": discovery["gene_family"],
|
||
"target_group": target_group,
|
||
"ref_group": ref_group,
|
||
"target_species": discovery["target_species"],
|
||
"ref_species": discovery["ref_species"],
|
||
"available_target_species": discovery.get("available_target_species", discovery["target_species"]),
|
||
"available_ref_species": discovery.get("available_ref_species", discovery["ref_species"]),
|
||
"selected_target_species": discovery.get("selected_target_species", discovery["target_species"]),
|
||
"selected_ref_species": discovery.get("selected_ref_species", discovery["ref_species"]),
|
||
"target_protein_spl": str(discovery["target_family_protein"]),
|
||
"ref_protein_spl": str(discovery["ref_family_protein"]),
|
||
"input_dir": str(input_dir),
|
||
"cpu": args.cpu,
|
||
}
|
||
report_path = write_report(workspace, args.prefix, metrics, parameters, tool_versions, start_time, end_time, (end_time - start_time).total_seconds(), cmd)
|
||
add_artifact(artifacts, "report", report_path, f"final/{args.prefix}.interspecies_synteny_report.md", "report", "report")
|
||
manifest_path = write_manifest(workspace, job_id, run_id, args.prefix, start_time, end_time, artifacts, metrics, parameters, tool_versions, workspace_root, input_dir, cmd, status=status)
|
||
add_artifact(artifacts, "manifest", manifest_path, f"final/{args.prefix}.interspecies_synteny.manifest.json", "manifest", "json")
|
||
write_manifest(workspace, job_id, run_id, args.prefix, start_time, end_time, artifacts, metrics, parameters, tool_versions, workspace_root, input_dir, cmd, status=status)
|
||
print(f"Interspecies synteny finished with status: {status}")
|
||
print(f"Manifest: {manifest_path}")
|
||
return 0 if metrics["tree_success"] else 1
|
||
except Exception as exc:
|
||
print(f"Interspecies synteny failed: {exc}")
|
||
try:
|
||
write_failure_manifest(workspace, job_id, run_id, args.prefix, start_time, workspace_root, input_dir, cmd, str(exc))
|
||
except Exception as manifest_exc:
|
||
print(f"Failed to write failure manifest: {manifest_exc}")
|
||
return 1
|
||
|
||
|
||
if __name__ == "__main__":
|
||
sys.exit(main())
|