2804 lines
129 KiB
Python
Executable File
2804 lines
129 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
|
||
import gzip
|
||
from datetime import datetime
|
||
from pathlib import Path
|
||
from typing import Dict, List, Any, Optional, Tuple
|
||
|
||
try:
|
||
from gfa_core.gene_family_scope import read_scope, sanitize_token
|
||
from gfa_core.gene_family_stage_context import apply_stage_context, build_stage_context, infer_analysis_route, member_set_id_for
|
||
except Exception:
|
||
def read_scope(scope_path):
|
||
return json.loads(Path(scope_path).read_text(encoding="utf-8"))
|
||
|
||
def sanitize_token(value: str) -> str:
|
||
token = re.sub(r"[^A-Za-z0-9_.-]+", "-", str(value or "").strip())
|
||
token = token.strip(".-_")
|
||
return token or "NA"
|
||
|
||
def infer_analysis_route(stage: str, scope_data: Dict[str, Any] | None = None) -> str:
|
||
scope_data = scope_data or {}
|
||
return "multi_species" if scope_data.get("reference_species") else "unscoped"
|
||
|
||
def member_set_id_for(prefix: str, scope_data: Dict[str, Any] | None, analysis_route: str) -> str:
|
||
scope_id = str((scope_data or {}).get("scope_id") or "").strip()
|
||
return f"{sanitize_token(scope_id or prefix)}:{analysis_route}:all_members"
|
||
|
||
def build_stage_context(*args, **kwargs):
|
||
return None
|
||
|
||
def apply_stage_context(manifest: dict, context: Any) -> dict:
|
||
return manifest
|
||
|
||
|
||
# ??conda??bin???PATH
|
||
CONDA_BIN = "/home/xzq/miniconda3/envs/test4py311/bin"
|
||
if CONDA_BIN not in os.environ.get("PATH", ""):
|
||
os.environ["PATH"] = CONDA_BIN + ":" + os.environ.get("PATH", "")
|
||
|
||
|
||
def open_text(path: Path, mode: str = "rt"):
|
||
path = Path(path)
|
||
if path.suffix == ".gz":
|
||
return gzip.open(path, mode, encoding="utf-8", errors="ignore")
|
||
return open(path, mode.replace("t", ""), encoding="utf-8", errors="ignore")
|
||
|
||
|
||
def materialize_plain_fasta(source: Path, output_dir: Path, label: str) -> Path:
|
||
"""Return a plain, search-safe protein FASTA for DIAMOND/BLASTP."""
|
||
source = Path(source)
|
||
output_dir.mkdir(parents=True, exist_ok=True)
|
||
stem = source.name[:-3] if source.suffix.lower() == ".gz" else source.name
|
||
stem = stem or f"{sanitize_token(label)}.fa"
|
||
dest = output_dir / f"{sanitize_token(label)}.{stem}"
|
||
if dest.exists() and dest.stat().st_size > 0 and dest.stat().st_mtime >= source.stat().st_mtime:
|
||
return dest
|
||
print(f" Materializing search-safe protein FASTA: {source} -> {dest}")
|
||
replacements = 0
|
||
sequence_lines = 0
|
||
valid_letters = set("ABCDEFGHIJKLMNOPQRSTUVWXYZ")
|
||
with open_text(source, "rt") as src, dest.open("w", encoding="utf-8") as dst:
|
||
for line in src:
|
||
if line.startswith(">"):
|
||
dst.write(line.rstrip("\n") + "\n")
|
||
continue
|
||
cleaned = []
|
||
for char in line.strip().upper():
|
||
if not char:
|
||
continue
|
||
if char in valid_letters:
|
||
cleaned.append(char)
|
||
else:
|
||
cleaned.append("X")
|
||
replacements += 1
|
||
if cleaned:
|
||
dst.write("".join(cleaned) + "\n")
|
||
sequence_lines += 1
|
||
if not dest.exists() or dest.stat().st_size == 0:
|
||
raise RuntimeError(f"Failed to materialize protein FASTA: {source}")
|
||
if replacements:
|
||
print(f" Replaced {replacements} non-letter protein character(s) with X in {dest.name}")
|
||
if sequence_lines == 0:
|
||
raise RuntimeError(f"No sequence lines found in protein FASTA: {source}")
|
||
return dest
|
||
|
||
|
||
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 current_input_dir(workspace: Path) -> Path:
|
||
return workspace.parent.parent / "input" / workspace.name
|
||
|
||
|
||
def looks_like_path_candidate(value: str) -> bool:
|
||
text = str(value or "").strip()
|
||
if not text or len(text) > 2048 or any(ch in text for ch in "\r\n\t"):
|
||
return False
|
||
if text.startswith(("/", "./", "../", "~")) or "/" in text or "\\" in text:
|
||
return True
|
||
if any(ch.isspace() for ch in text):
|
||
return False
|
||
return bool(Path(text).suffix)
|
||
|
||
|
||
def rebase_current_job_path(workspace: Path, raw_path: Any, input_dir: Optional[Path] = None) -> Optional[Path]:
|
||
path_text = str(raw_path or "").strip()
|
||
if not looks_like_path_candidate(path_text):
|
||
return None
|
||
path = Path(path_text)
|
||
input_dir = input_dir or current_input_dir(workspace)
|
||
parts = list(path.parts)
|
||
for marker, root in (("input", input_dir), ("workspace", workspace)):
|
||
if marker in parts:
|
||
idx = parts.index(marker)
|
||
if idx + 1 < len(parts) and parts[idx + 1] == workspace.name:
|
||
candidate = root.joinpath(*parts[idx + 2 :])
|
||
if candidate.exists():
|
||
return candidate
|
||
if not path.is_absolute():
|
||
for candidate in (Path.cwd() / path, input_dir / path, input_dir / path.name, workspace / path):
|
||
if candidate.exists():
|
||
return candidate
|
||
for candidate in (
|
||
input_dir / path.name,
|
||
workspace / "final" / path.name,
|
||
workspace / "interspecies_synteny" / "analysis" / path.name,
|
||
):
|
||
if candidate.exists():
|
||
return candidate
|
||
if path.exists():
|
||
return path
|
||
return None
|
||
|
||
|
||
def rebase_job_path_text(workspace: Path, value: str) -> str:
|
||
if not value:
|
||
return value
|
||
if looks_like_path_candidate(value):
|
||
direct = rebase_current_job_path(workspace, value)
|
||
if direct:
|
||
return str(direct.absolute())
|
||
job_id = re.escape(workspace.name)
|
||
pattern = re.compile(rf"/[^\s\"']*?(?:input|workspace)/{job_id}/[^\s\"']+")
|
||
|
||
def replace_match(match: re.Match[str]) -> str:
|
||
original = match.group(0).rstrip(",)]}")
|
||
suffix = match.group(0)[len(original):]
|
||
rebased = rebase_current_job_path(workspace, original)
|
||
if rebased:
|
||
return str(rebased.absolute()) + suffix
|
||
return match.group(0)
|
||
|
||
return pattern.sub(replace_match, value)
|
||
|
||
|
||
def rebase_manifest_payload_paths(workspace: Path, value: Any) -> Any:
|
||
if isinstance(value, dict):
|
||
return {key: rebase_manifest_payload_paths(workspace, item) for key, item in value.items()}
|
||
if isinstance(value, list):
|
||
return [rebase_manifest_payload_paths(workspace, item) for item in value]
|
||
if isinstance(value, str):
|
||
return rebase_job_path_text(workspace, value)
|
||
return value
|
||
|
||
|
||
def summarize_process_result(result: subprocess.CompletedProcess) -> Dict[str, Any]:
|
||
"""Keep process diagnostics compact and avoid persisting wrapped absolute paths."""
|
||
def clean_lines(text: str) -> List[str]:
|
||
cleaned: List[str] = []
|
||
for line in (text or "").splitlines():
|
||
if any(marker in line for marker in ("/data_old", ".bak", "bbav5_0502_Yang_Tihuan.7880")):
|
||
continue
|
||
stripped = line.strip()
|
||
if stripped:
|
||
cleaned.append(stripped)
|
||
return cleaned[-20:]
|
||
|
||
stdout = "\n".join(clean_lines(result.stdout))
|
||
stderr = "\n".join(clean_lines(result.stderr))
|
||
return {
|
||
"returncode": result.returncode,
|
||
"stdout_tail": stdout[:2000],
|
||
"stderr_tail": stderr[:2000],
|
||
}
|
||
|
||
|
||
def resolve_iqtree_executable() -> str:
|
||
"""Resolve IQ-TREE consistently with the phylogeny stage."""
|
||
candidates = [
|
||
Path(sys.executable).resolve().parent / "iqtree3",
|
||
Path(sys.executable).resolve().parent / "iqtree2",
|
||
shutil.which("iqtree3"),
|
||
shutil.which("iqtree2"),
|
||
shutil.which("iqtree"),
|
||
]
|
||
for candidate in candidates:
|
||
if not candidate:
|
||
continue
|
||
path = Path(str(candidate))
|
||
if path.exists() and os.access(path, os.X_OK):
|
||
return str(path)
|
||
if isinstance(candidate, str) and shutil.which(candidate):
|
||
return str(candidate)
|
||
return "iqtree3"
|
||
|
||
|
||
def get_tool_versions() -> Dict[str, str]:
|
||
"""获取所需工具版本信息"""
|
||
versions = {}
|
||
|
||
for tool, cmd in [
|
||
("mafft", ["mafft", "--version"]),
|
||
("iqtree", [resolve_iqtree_executable(), "--version"]),
|
||
("diamond", ["diamond", "--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_text(target_protein, 'rt') as f:
|
||
outfile.write(f.read().rstrip() + "\n")
|
||
with open_text(ref_protein, 'rt') as f:
|
||
outfile.write(f.read().rstrip() + "\n")
|
||
|
||
target_count = sum(1 for line in open_text(target_protein, 'rt') if line.startswith('>'))
|
||
ref_count = sum(1 for line in open_text(ref_protein, 'rt') 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"构建进化树...")
|
||
expected_tips = sum(1 for line in open_text(merged_fasta, "rt") if line.startswith(">"))
|
||
|
||
# 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"
|
||
tree_file = Path(f"{base_name}.treefile")
|
||
existing_tips = _read_tree_tip_labels(tree_file) if tree_file.exists() else set()
|
||
if existing_tips and len(existing_tips) != expected_tips:
|
||
print(
|
||
f" Removing stale IQ-TREE outputs ({len(existing_tips)} tips != {expected_tips} merged sequences)"
|
||
)
|
||
for stale in output_dir.glob(f"{prefix}_interspecies.iqtree*"):
|
||
try:
|
||
stale.unlink()
|
||
except OSError:
|
||
pass
|
||
iqtree_cmd = [
|
||
resolve_iqtree_executable(), "-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)
|
||
if result.returncode != 0:
|
||
return {"success": False, "error": f"IQ-TREE失败: {result.stderr or result.stdout}"}
|
||
|
||
if not tree_file.exists():
|
||
return {"success": False, "error": "IQ-TREE未生成 treefile"}
|
||
built_tips = _read_tree_tip_labels(tree_file)
|
||
if len(built_tips) != expected_tips:
|
||
return {
|
||
"success": False,
|
||
"error": f"IQ-TREE tree tip count ({len(built_tips)}) != merged sequence count ({expected_tips})",
|
||
}
|
||
|
||
return {
|
||
"success": True,
|
||
"tree_file": tree_file,
|
||
"msa_file": msa_file,
|
||
"tip_count": len(built_tips),
|
||
"iqtree_executable": iqtree_cmd[0],
|
||
"command": " ".join(iqtree_cmd),
|
||
}
|
||
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_text(gff3_file, 'rt') 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_text(fasta_file, "rt") 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)
|
||
seeds = {norm}
|
||
if "|" in (gene_id or ""):
|
||
seeds.add(normalize_seq_id(gene_id.split("|", 1)[0]))
|
||
expanded = set()
|
||
for seed in seeds:
|
||
if not seed:
|
||
continue
|
||
expanded.update({seed, seed.upper(), seed.lower()})
|
||
if re.search(r"\.\d+$", seed) and not seed.lower().startswith(("evm.model.", "evm.tu.")):
|
||
stripped = re.sub(r"\.\d+$", "", seed)
|
||
expanded.update({stripped, stripped.upper(), stripped.lower()})
|
||
suffix_stripped = re.sub(r"_[PT]\d+$", "", seed)
|
||
expanded.update({suffix_stripped, suffix_stripped.upper(), suffix_stripped.lower()})
|
||
prefixed = set(expanded)
|
||
for value in expanded:
|
||
for prefix in ("gene:", "transcript:", "CDS:", "rna-", "gene-", "cds-"):
|
||
prefixed.add(normalize_seq_id(prefix + value))
|
||
return {v for v in prefixed if v}
|
||
|
||
|
||
def fasta_ids(fasta_file: Path) -> List[str]:
|
||
ids = []
|
||
with open_text(fasta_file, "rt") 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 _merge_blast_outputs(blast_files: List[Path], merged_blast: Path) -> bool:
|
||
with open(merged_blast, "w") as outfile:
|
||
for blast_file in blast_files:
|
||
if blast_file.exists():
|
||
with open(blast_file, "r") as infile:
|
||
for line in infile:
|
||
if line.strip():
|
||
outfile.write(line)
|
||
return merged_blast.exists() and merged_blast.stat().st_size > 0
|
||
|
||
|
||
def _run_reciprocal_blastp_engine(
|
||
target_protein: Path,
|
||
ref_protein: Path,
|
||
output_dir: Path,
|
||
target_species: str,
|
||
ref_species: str,
|
||
cpu: int,
|
||
evalue: str,
|
||
max_target_seqs: int,
|
||
) -> Dict[str, Any]:
|
||
pair_prefix = f"{target_species}.{ref_species}"
|
||
target_db = output_dir / f"{target_species}.blastdb"
|
||
ref_db = output_dir / f"{ref_species}.blastdb"
|
||
commands = []
|
||
stderr_items = []
|
||
start_time = time.time()
|
||
|
||
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))
|
||
try:
|
||
result = subprocess.run(cmd, capture_output=True, text=True, check=False)
|
||
except FileNotFoundError:
|
||
return {"success": False, "stderr": "makeblastdb not found", "commands": commands, "engine_used": "blastp"}
|
||
if result.stderr:
|
||
stderr_items.append(result.stderr)
|
||
if result.returncode != 0:
|
||
return {"success": False, "stderr": "\n".join(stderr_items), "commands": commands, "engine_used": "blastp"}
|
||
|
||
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", str(evalue),
|
||
"-num_alignments", str(max_target_seqs), "-num_threads", str(cpu),
|
||
]
|
||
commands.append(" ".join(cmd))
|
||
try:
|
||
result = subprocess.run(cmd, capture_output=True, text=True, check=False)
|
||
except FileNotFoundError:
|
||
return {"success": False, "stderr": "blastp not found", "commands": commands, "engine_used": "blastp"}
|
||
if result.stderr:
|
||
stderr_items.append(result.stderr)
|
||
if result.returncode != 0:
|
||
return {"success": False, "stderr": "\n".join(stderr_items), "commands": commands, "engine_used": "blastp"}
|
||
|
||
merged_blast = output_dir / f"{pair_prefix}.blast"
|
||
success = _merge_blast_outputs([target_to_ref, ref_to_target], merged_blast)
|
||
return {
|
||
"success": success,
|
||
"output_file": merged_blast if success else None,
|
||
"commands": commands,
|
||
"stderr": "\n".join(stderr_items),
|
||
"output_stats": {
|
||
str(path.name): path.stat().st_size if path.exists() else -1
|
||
for path in [target_to_ref, ref_to_target, merged_blast]
|
||
},
|
||
"engine_used": "blastp",
|
||
"elapsed_sec": time.time() - start_time,
|
||
}
|
||
|
||
|
||
def _run_reciprocal_diamond_engine(
|
||
target_protein: Path,
|
||
ref_protein: Path,
|
||
output_dir: Path,
|
||
target_species: str,
|
||
ref_species: str,
|
||
cpu: int,
|
||
evalue: str,
|
||
max_target_seqs: int,
|
||
) -> Dict[str, Any]:
|
||
pair_prefix = f"{target_species}.{ref_species}"
|
||
target_db = output_dir / f"{target_species}.diamond_db"
|
||
ref_db = output_dir / f"{ref_species}.diamond_db"
|
||
commands = []
|
||
stderr_items = []
|
||
start_time = time.time()
|
||
|
||
for protein_file, db_prefix in [(target_protein, target_db), (ref_protein, ref_db)]:
|
||
cmd = ["diamond", "makedb", "--in", str(protein_file), "--db", str(db_prefix)]
|
||
commands.append(" ".join(cmd))
|
||
try:
|
||
result = subprocess.run(cmd, capture_output=True, text=True, check=False)
|
||
except FileNotFoundError:
|
||
return {"success": False, "stderr": "diamond not found", "commands": commands, "engine_used": "diamond"}
|
||
if result.stderr:
|
||
stderr_items.append(result.stderr)
|
||
if result.returncode != 0:
|
||
return {"success": False, "stderr": "\n".join(stderr_items), "commands": commands, "engine_used": "diamond"}
|
||
|
||
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 = [
|
||
"diamond", "blastp", "--query", str(query), "--db", str(db_prefix),
|
||
"--out", str(out_file), "--outfmt", "6", "--max-target-seqs", str(max_target_seqs),
|
||
"--evalue", str(evalue), "--threads", str(cpu),
|
||
]
|
||
commands.append(" ".join(cmd))
|
||
try:
|
||
result = subprocess.run(cmd, capture_output=True, text=True, check=False)
|
||
except FileNotFoundError:
|
||
return {"success": False, "stderr": "diamond not found", "commands": commands, "engine_used": "diamond"}
|
||
if result.stderr:
|
||
stderr_items.append(result.stderr)
|
||
if result.returncode != 0:
|
||
return {"success": False, "stderr": "\n".join(stderr_items), "commands": commands, "engine_used": "diamond"}
|
||
|
||
merged_blast = output_dir / f"{pair_prefix}.blast"
|
||
success = _merge_blast_outputs([target_to_ref, ref_to_target], merged_blast)
|
||
return {
|
||
"success": success,
|
||
"output_file": merged_blast if success else None,
|
||
"commands": commands,
|
||
"stderr": "\n".join(stderr_items),
|
||
"output_stats": {
|
||
str(path.name): path.stat().st_size if path.exists() else -1
|
||
for path in [target_to_ref, ref_to_target, merged_blast]
|
||
},
|
||
"engine_used": "diamond",
|
||
"elapsed_sec": time.time() - start_time,
|
||
}
|
||
|
||
|
||
def run_reciprocal_blast(
|
||
target_protein: Path,
|
||
ref_protein: Path,
|
||
output_dir: Path,
|
||
target_species: str,
|
||
ref_species: str,
|
||
cpu: int = 8,
|
||
engine: str = "diamond",
|
||
evalue: str = "1e-10",
|
||
max_target_seqs: int = 5,
|
||
) -> Dict[str, Any]:
|
||
"""Run pairwise reciprocal protein search for one target-ref species pair."""
|
||
target_protein = materialize_plain_fasta(target_protein, output_dir, f"{target_species}.search")
|
||
ref_protein = materialize_plain_fasta(ref_protein, output_dir, f"{ref_species}.search")
|
||
requested_engine = (engine or "diamond").strip().lower()
|
||
if requested_engine not in {"diamond", "blastp"}:
|
||
requested_engine = "diamond"
|
||
print(f"Reciprocal protein search ({requested_engine}): {target_species} vs {ref_species}")
|
||
|
||
if requested_engine == "blastp":
|
||
result = _run_reciprocal_blastp_engine(target_protein, ref_protein, output_dir, target_species, ref_species, cpu, evalue, max_target_seqs)
|
||
result["engine_requested"] = requested_engine
|
||
result["fallback_from"] = None
|
||
return result
|
||
|
||
result = _run_reciprocal_diamond_engine(target_protein, ref_protein, output_dir, target_species, ref_species, cpu, evalue, max_target_seqs)
|
||
if result.get("success"):
|
||
result["engine_requested"] = requested_engine
|
||
result["fallback_from"] = None
|
||
return result
|
||
|
||
diamond_stderr = result.get("stderr", "")
|
||
diamond_commands = list(result.get("commands", []))
|
||
print(f" Diamond failed for {target_species}.{ref_species}; falling back to NCBI blastp")
|
||
fallback = _run_reciprocal_blastp_engine(target_protein, ref_protein, output_dir, target_species, ref_species, cpu, evalue, max_target_seqs)
|
||
fallback["engine_requested"] = requested_engine
|
||
fallback["fallback_from"] = "diamond"
|
||
fallback["diamond_stderr"] = diamond_stderr
|
||
fallback["commands"] = diamond_commands + list(fallback.get("commands", []))
|
||
return fallback
|
||
|
||
|
||
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 _read_tree_tip_labels(tree_path: Path) -> set[str]:
|
||
try:
|
||
content = tree_path.read_text(encoding="utf-8", errors="ignore")
|
||
except Exception:
|
||
return set()
|
||
return {m.group(1).strip().strip("'").strip('"') for m in re.finditer(r"(?<=[(,])([^():,;]+):", content) if m.group(1).strip()}
|
||
|
||
|
||
def _merge_matching_subfamily_files(candidates: List[Path], tree_file: Path, output_dir: Path, prefix: str) -> Optional[Path]:
|
||
tips = _read_tree_tip_labels(tree_file)
|
||
if not tips:
|
||
return None
|
||
merged: Dict[str, str] = {}
|
||
used_files: List[str] = []
|
||
seen: set[Path] = set()
|
||
for candidate in candidates:
|
||
if not candidate or candidate in seen or not candidate.exists() or not candidate.is_file():
|
||
continue
|
||
seen.add(candidate)
|
||
matched_here = False
|
||
try:
|
||
with candidate.open("r", encoding="utf-8", errors="ignore") as handle:
|
||
for raw in handle:
|
||
parts = raw.strip().split()
|
||
if len(parts) < 2:
|
||
continue
|
||
gene_id, subfamily = parts[0], parts[1]
|
||
if gene_id in tips and gene_id not in merged:
|
||
merged[gene_id] = subfamily
|
||
matched_here = True
|
||
except Exception:
|
||
continue
|
||
if matched_here:
|
||
used_files.append(str(candidate))
|
||
if not merged:
|
||
return None
|
||
output_dir.mkdir(parents=True, exist_ok=True)
|
||
merged_path = output_dir / f"{prefix}.merged.geneID.subfamily"
|
||
with merged_path.open("w", encoding="utf-8") as handle:
|
||
for gene_id in sorted(merged):
|
||
handle.write(f"{gene_id}\t{merged[gene_id]}\n")
|
||
print(f" 合并亚家族文件: {merged_path} ({len(merged)}/{len(tips)} tips; sources={len(used_files)})")
|
||
return merged_path
|
||
|
||
|
||
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,
|
||
subfamily_arc_enabled: bool = True,
|
||
subfamily_arc_gap: float = 1.2,
|
||
subfamily_arc_label_gap: float = 0.6,
|
||
subfamily_arc_label_size: float = 0.0,
|
||
) -> Dict[str, Any]:
|
||
"""绘制物种间进化树图 - 使用专用的 interspecies_tree.R 脚本"""
|
||
print(f"绘制进化树图...")
|
||
|
||
# 查找专用的种间进化树R脚本
|
||
r_script = Path("scripts/interspecies_tree.R")
|
||
if not r_script.exists():
|
||
return {"success": False, "error": f"R脚本不存在: {r_script}"}
|
||
|
||
try:
|
||
# 使用专用的 interspecies_tree.R 脚本，传递物种参数和亚家族文件
|
||
# 使用conda环境的Rscript
|
||
conda_rscript = "/home/xzq/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("")
|
||
cmd.extend([
|
||
str(0.4),
|
||
str(2.8),
|
||
"black",
|
||
str(0.5),
|
||
"triangle",
|
||
"circle",
|
||
str(3.0),
|
||
str(2.0),
|
||
"true" if subfamily_arc_enabled else "false",
|
||
str(subfamily_arc_gap),
|
||
str(subfamily_arc_label_gap),
|
||
str(subfamily_arc_label_size if subfamily_arc_label_size and subfamily_arc_label_size > 0 else ""),
|
||
])
|
||
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"
|
||
output_circular_pdf = output_dir / "interspecies_tree.circular_tree.pdf"
|
||
output_circular_png = output_dir / "interspecies_tree.circular_tree.png"
|
||
|
||
# 重命名为规范名称
|
||
final_pdf = output_dir / f"{prefix}.interspecies_tree.pdf"
|
||
final_png = output_dir / f"{prefix}.interspecies_tree.png"
|
||
final_circular_pdf = output_dir / f"{prefix}.interspecies_tree.circular_tree.pdf"
|
||
final_circular_png = output_dir / f"{prefix}.interspecies_tree.circular_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}")
|
||
if output_circular_pdf.exists():
|
||
shutil.move(output_circular_pdf, final_circular_pdf)
|
||
print(f" 进化树圈图PDF: {final_circular_pdf}")
|
||
if output_circular_png.exists():
|
||
shutil.move(output_circular_png, final_circular_png)
|
||
print(f" 进化树圈图PNG: {final_circular_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,
|
||
"circular_pdf_file": final_circular_pdf if final_circular_pdf.exists() else None,
|
||
"circular_png_file": final_circular_png if final_circular_png.exists() else None,
|
||
"stdout": result.stdout,
|
||
"stderr": result.stderr,
|
||
"subfamily_arc_enabled": subfamily_arc_enabled,
|
||
"subfamily_arc_gap": subfamily_arc_gap,
|
||
"subfamily_arc_label_gap": subfamily_arc_label_gap,
|
||
"subfamily_arc_label_size": subfamily_arc_label_size,
|
||
}
|
||
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/xzq/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.absolute()), str(anchors_new_file.absolute()),
|
||
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)
|
||
screen_summary = summarize_process_result(result)
|
||
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": screen_summary.get("stderr_tail", ""),
|
||
}
|
||
|
||
|
||
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 _family_block_signatures_from_anchors(anchors_file: Path, family_ids: set) -> set:
|
||
"""Return JCVI simple-row signatures for blocks containing family genes.
|
||
|
||
JCVI's anchors.new keeps full block membership, while anchors.simple stores
|
||
only block endpoints. Published-style family highlighting should mark the
|
||
whole block when any family member occurs inside the block, not only when the
|
||
member happens to be one of the simple endpoints.
|
||
"""
|
||
signatures = set()
|
||
if not anchors_file or not anchors_file.exists() or not family_ids:
|
||
return signatures
|
||
|
||
def flush(block: List[Tuple[str, str]]) -> None:
|
||
if not block:
|
||
return
|
||
contains_family = any(_is_family_gene(q, family_ids) or _is_family_gene(s, family_ids) for q, s in block)
|
||
if not contains_family:
|
||
return
|
||
q1, s1 = block[0]
|
||
q2, s2 = block[-1]
|
||
signatures.add((frozenset((q1, q2)), frozenset((s1, s2))))
|
||
|
||
block: List[Tuple[str, str]] = []
|
||
with open(anchors_file, "r") as f:
|
||
for raw in f:
|
||
line = raw.strip()
|
||
if not line:
|
||
continue
|
||
if line.startswith("###") or line.startswith("#"):
|
||
flush(block)
|
||
block = []
|
||
continue
|
||
parts = line.split()
|
||
if len(parts) < 2:
|
||
continue
|
||
block.append((normalize_seq_id(parts[0]), normalize_seq_id(parts[1])))
|
||
flush(block)
|
||
return signatures
|
||
|
||
|
||
def _simple_block_signature(parts: List[str]) -> Optional[Tuple[frozenset, frozenset]]:
|
||
if len(parts) < 4:
|
||
return None
|
||
q1, q2, s1, s2 = [normalize_seq_id(p.split("*", 1)[-1]) for p in parts[:4]]
|
||
return (frozenset((q1, q2)), frozenset((s1, s2)))
|
||
|
||
|
||
def _write_marked_simple_file(
|
||
simple_file: Path,
|
||
marked_file: Path,
|
||
family_ids: set,
|
||
family_links: List[Tuple[str, str, int]],
|
||
anchors_new_file: Optional[Path] = None,
|
||
) -> int:
|
||
if not simple_file.exists() or not family_ids:
|
||
return 0
|
||
family_block_signatures = _family_block_signatures_from_anchors(anchors_new_file, family_ids) if anchors_new_file else set()
|
||
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)
|
||
row_sig = _simple_block_signature(parts)
|
||
seen_rows.add(row_key)
|
||
direct_endpoint_hit = any(_is_family_gene(g, family_ids) for g in genes)
|
||
block_hit = row_sig in family_block_signatures if row_sig is not None else False
|
||
if (direct_endpoint_hit or block_hit) 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/xzq/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, Path(pair.get("anchors_new_file", "")))
|
||
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.absolute()), str(layout.absolute()), "--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.absolute())
|
||
if (work_dir / "karyotype.png").exists():
|
||
shutil.copy(work_dir / "karyotype.png", final_png.absolute())
|
||
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 summarize_process_result(result)
|
||
|
||
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 plot_pairwise_synteny(
|
||
work_dir: Path,
|
||
output_dir: Path,
|
||
prefix: str,
|
||
pair_results: List[Dict[str, Any]],
|
||
family_ids: set,
|
||
) -> Dict[str, Any]:
|
||
"""Draw one JCVI karyotype figure for each successful target-ref pair."""
|
||
print("Plotting pairwise JCVI synteny...")
|
||
conda_python = "/home/xzq/miniconda3/envs/test4py311/bin/python"
|
||
usable_pairs = [p for p in pair_results if p.get("success") and p.get("simple_file") and Path(p["simple_file"]).exists()]
|
||
if not usable_pairs:
|
||
return {"success": False, "error": "No successful pairwise anchors.simple files available for JCVI", "pairs": []}
|
||
|
||
pair_outputs: List[Dict[str, Any]] = []
|
||
palette = ["#d73027", "#4575b4"]
|
||
|
||
def run_plot(seqids: Path, 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.absolute()), str(layout.absolute()), "--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.absolute())
|
||
if (work_dir / "karyotype.png").exists():
|
||
shutil.copy(work_dir / "karyotype.png", final_png.absolute())
|
||
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 summarize_process_result(result)
|
||
|
||
for pair in usable_pairs:
|
||
target = pair["target_species"]
|
||
ref = pair["ref_species"]
|
||
pair_key = f"{target}.{ref}"
|
||
target_bed = work_dir / f"{target}.bed"
|
||
ref_bed = work_dir / f"{ref}.bed"
|
||
simple = Path(pair["simple_file"])
|
||
if not target_bed.exists() or not ref_bed.exists():
|
||
pair_outputs.append({"pair": pair_key, "success": False, "error": "missing pair BED file"})
|
||
continue
|
||
|
||
marked = work_dir / f"{pair_key}.anchors.simple.mark"
|
||
target_bed_ids = _bed_gene_ids(target_bed)
|
||
ref_bed_ids = _bed_gene_ids(ref_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, Path(pair.get("anchors_new_file", "")))
|
||
pair["pairwise_family_highlight_links"] = len(family_links)
|
||
pair["pairwise_marked_block_count"] = marked_count
|
||
if marked_count > 0:
|
||
pair["marked_simple_file"] = str(marked)
|
||
|
||
seqids_file = work_dir / f"{pair_key}.seqids"
|
||
layout_file = work_dir / f"{pair_key}.layout"
|
||
layout_mark_file = work_dir / f"{pair_key}.layout_mark"
|
||
species_rows = [(target, target_bed), (ref, ref_bed)]
|
||
with open(seqids_file, "w") as f:
|
||
for _, bed in species_rows:
|
||
chroms = _read_bed_chroms(bed)
|
||
if not chroms:
|
||
pair_outputs.append({"pair": pair_key, "success": False, "error": f"no chromosomes found in {bed.name}"})
|
||
break
|
||
f.write(",".join(chroms) + "\n")
|
||
if pair_outputs and pair_outputs[-1].get("pair") == pair_key and not pair_outputs[-1].get("success") and pair_outputs[-1].get("error", "").startswith("no chromosomes"):
|
||
continue
|
||
|
||
def write_layout(layout: Path, use_marked: bool) -> None:
|
||
edge_file = Path(pair.get("marked_simple_file") if use_marked and pair.get("marked_simple_file") else simple).name
|
||
with open(layout, "w") as f:
|
||
f.write("# y, xstart, xend, rotation, color, label, va, bed\n")
|
||
f.write(f"0.75, .08, .92, 0, {palette[0]}, {target}, top, {target}.bed\n")
|
||
f.write(f"0.25, .08, .92, 0, {palette[1]}, {ref}, bottom, {ref}.bed\n")
|
||
f.write("# edges\n")
|
||
f.write(f"e, 0, 1, {edge_file}\n")
|
||
|
||
write_layout(layout_file, use_marked=False)
|
||
write_layout(layout_mark_file, use_marked=True)
|
||
|
||
final_pdf = output_dir / f"{prefix}.{pair_key}.interspecies_synteny.pdf"
|
||
final_png = output_dir / f"{prefix}.{pair_key}.interspecies_synteny.png"
|
||
final_mark_pdf = output_dir / f"{prefix}.{pair_key}.interspecies_synteny.marked.pdf"
|
||
final_mark_png = output_dir / f"{prefix}.{pair_key}.interspecies_synteny.marked.png"
|
||
base_result = run_plot(seqids_file, layout_file, final_pdf, final_png)
|
||
mark_result = run_plot(seqids_file, layout_mark_file, final_mark_pdf, final_mark_png)
|
||
success = final_pdf.exists() or final_png.exists()
|
||
pair_outputs.append({
|
||
"pair": pair_key,
|
||
"target_species": target,
|
||
"ref_species": ref,
|
||
"success": success,
|
||
"pdf_file": str(final_pdf) if final_pdf.exists() else "",
|
||
"png_file": str(final_png) if final_png.exists() else "",
|
||
"marked_pdf_file": str(final_mark_pdf) if final_mark_pdf.exists() else "",
|
||
"marked_png_file": str(final_mark_png) if final_mark_png.exists() else "",
|
||
"seqids_file": str(seqids_file),
|
||
"layout_file": str(layout_file),
|
||
"layout_mark_file": str(layout_mark_file),
|
||
"base_plot": base_result,
|
||
"marked_plot": mark_result,
|
||
})
|
||
|
||
return {
|
||
"success": any(item.get("success") for item in pair_outputs),
|
||
"pairs": pair_outputs,
|
||
"pairwise_plot_count": sum(1 for item in pair_outputs if item.get("success")),
|
||
}
|
||
|
||
|
||
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(f"- BLAST engine requested: {parameters.get('blast_engine_requested')}\n")
|
||
f.write(f"- BLAST engine used: {metrics.get('blast_engine_used')}\n")
|
||
f.write(f"- BLAST fallback count: {metrics.get('blast_fallback_count', 0)}\n")
|
||
species_positions = metrics.get("species_position_results") or {}
|
||
if species_positions:
|
||
f.write("- 染色体背景: 保留主核染色体/染色体尺度序列，过滤疑似细胞器和短 contig/scaffold；无法识别主序列时回退全量。\n")
|
||
for species, info in species_positions.items():
|
||
f.write(
|
||
f" - {species}: plotted {info.get('chromosomes_plotted', 'N/A')}/"
|
||
f"{info.get('chromosomes_total', 'N/A')} seqids; "
|
||
f"filtered genes {info.get('genes_filtered_by_chromosome', 'N/A')} "
|
||
f"({info.get('chromosomes_filter_mode', 'N/A')})\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"]
|
||
}
|
||
|
||
formatted_artifacts = rebase_manifest_payload_paths(workspace, formatted_artifacts)
|
||
rebased_metrics = rebase_manifest_payload_paths(workspace, metrics)
|
||
rebased_parameters = rebase_manifest_payload_paths(workspace, parameters)
|
||
|
||
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": rebase_manifest_payload_paths(workspace, 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": rebased_parameters,
|
||
"metrics": rebased_metrics
|
||
}
|
||
for key in (
|
||
"scope_id",
|
||
"scope_mode",
|
||
"scope_file",
|
||
"target_species",
|
||
"reference_species",
|
||
"comparative_species",
|
||
"selected_species",
|
||
"family",
|
||
"analysis_route",
|
||
"member_set_id",
|
||
):
|
||
if key in rebased_parameters:
|
||
manifest_data[key] = rebased_parameters[key]
|
||
|
||
if rebased_parameters.get("scope_id"):
|
||
scope_context_data = {
|
||
"scope_id": rebased_parameters.get("scope_id", ""),
|
||
"scope_mode": rebased_parameters.get("scope_mode", ""),
|
||
"family": rebased_parameters.get("family", ""),
|
||
"target_species": rebased_parameters.get("target_species", []),
|
||
"reference_species": rebased_parameters.get("reference_species", []),
|
||
"comparative_species": rebased_parameters.get("comparative_species", []),
|
||
"selected_species": rebased_parameters.get("selected_species", []),
|
||
}
|
||
context = build_stage_context(
|
||
stage="interspecies_synteny",
|
||
prefix=prefix,
|
||
scope_data=scope_context_data,
|
||
scope_file=rebased_parameters.get("scope_file", ""),
|
||
output_prefix=prefix,
|
||
analysis_route=rebased_parameters.get("analysis_route", ""),
|
||
member_set_id=rebased_parameters.get("member_set_id", ""),
|
||
)
|
||
if context:
|
||
manifest_data = apply_stage_context(manifest_data, context)
|
||
|
||
manifest_path = workspace / "final" / f"{prefix}.interspecies_synteny.manifest.json"
|
||
manifest_path.parent.mkdir(parents=True, exist_ok=True)
|
||
with open(manifest_path, "w", encoding='utf-8') as f:
|
||
json.dump(manifest_data, f, indent=2, ensure_ascii=False)
|
||
if rebased_parameters.get("scope_id"):
|
||
scoped_manifest_path = workspace / "final" / f"{prefix}.{sanitize_token(rebased_parameters.get('scope_id', ''))}.interspecies_synteny.manifest.json"
|
||
manifest_data["scope_manifest"] = str(scoped_manifest_path.absolute())
|
||
with open(scoped_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 _clean_text(value: Any) -> str:
|
||
if value is None:
|
||
return ""
|
||
try:
|
||
import pandas as pd
|
||
if pd.isna(value):
|
||
return ""
|
||
except Exception:
|
||
pass
|
||
text = str(value).strip()
|
||
return "" if text.lower() in {"nan", "none", "null", "undefined", "-"} else text
|
||
|
||
|
||
def _gene_family_from_prefix(prefix: str) -> str:
|
||
return prefix.split('.', 1)[1] if '.' in prefix else prefix
|
||
|
||
|
||
def _scope_file_record(workspace: Path, input_dir: Path, scope_data: Dict[str, Any], species: str, file_type: str) -> Optional[Path]:
|
||
for record in ((scope_data.get("files") or {}).get(species) or {}).get(file_type, []) or []:
|
||
path_text = record.get("path") or record.get("full_path") or ""
|
||
if path_text:
|
||
path = rebase_current_job_path(workspace, path_text, input_dir=input_dir) or Path(path_text)
|
||
if path.exists():
|
||
return path
|
||
return None
|
||
|
||
|
||
def _add_scope_species_to_groups(workspace: Path, input_dir: Path, groups: Dict[str, Any], scope_data: Dict[str, Any]) -> None:
|
||
family = scope_data.get("family") or groups.get("gene_family") or ""
|
||
target_set = set(scope_data.get("target_species") or [])
|
||
reference_set = set(scope_data.get("reference_species") or [])
|
||
comparative_set = set(scope_data.get("comparative_species") or [])
|
||
for species in target_set | reference_set | comparative_set:
|
||
group_name = "target" if species in target_set else "ref"
|
||
item = groups[group_name].setdefault(species, {"species": species, "gene_family": family, "files": {}})
|
||
item["gene_family"] = item.get("gene_family") or family
|
||
item["role"] = "target" if species in target_set else ("reference" if species in reference_set else "comparative")
|
||
files = item.setdefault("files", {})
|
||
for file_type in ("pep", "family_pep", "gff"):
|
||
path = _scope_file_record(workspace, input_dir, scope_data, species, file_type)
|
||
if path:
|
||
key = "all_pep" if file_type == "pep" else file_type
|
||
files[key] = str(path)
|
||
|
||
|
||
def _apply_scope_selection_to_args(args, scope_data: Dict[str, Any]) -> None:
|
||
targets = [str(x).strip() for x in (scope_data.get("target_species") or []) if str(x).strip()]
|
||
selected = [str(x).strip() for x in (scope_data.get("selected_species") or []) if str(x).strip()]
|
||
refs = [str(x).strip() for x in (scope_data.get("reference_species") or []) if str(x).strip()]
|
||
comps = [str(x).strip() for x in (scope_data.get("comparative_species") or []) if str(x).strip()]
|
||
selected_ref_side = [sp for sp in selected if sp not in set(targets)]
|
||
if not selected_ref_side:
|
||
selected_ref_side = refs + comps
|
||
if targets and not getattr(args, "target_species", None):
|
||
args.target_species = ",".join(targets)
|
||
if selected_ref_side and not getattr(args, "ref_species", None):
|
||
args.ref_species = ",".join(dict.fromkeys(selected_ref_side))
|
||
|
||
|
||
def _resolve_path(path_text: str, input_dir: Path, workspace: Optional[Path] = None) -> Optional[Path]:
|
||
if not path_text:
|
||
return None
|
||
if workspace:
|
||
rebased = rebase_current_job_path(workspace, path_text, input_dir=input_dir)
|
||
if rebased:
|
||
return rebased
|
||
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": {}, "seed_reference": {}, "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 = _clean_text(row.get("Species_Code", ""))
|
||
row_family = _clean_text(row.get("Gene_Family", ""))
|
||
family = row_family or default_family
|
||
raw_file_type = _clean_text(row.get("File_Type", "")).lower()
|
||
file_scope = _clean_text(row.get("File_Scope", "")).lower()
|
||
file_type = raw_file_type
|
||
original_name = _clean_text(row.get("Original_File_Name", ""))
|
||
role = _clean_text(row.get("Species_Role", "")).lower()
|
||
if not species:
|
||
continue
|
||
if role in {"seed", "seed_reference", "known_family_reference"}:
|
||
group_name = "seed_reference"
|
||
elif role in {"comparative", "comparison", "outgroup", "reference", "ref", "comparison_reference"}:
|
||
group_name = "ref"
|
||
else:
|
||
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
|
||
item["role"] = role or group_name
|
||
if raw_file_type in {"pep", "protein"}:
|
||
if file_scope in {"known_family", "family", "seed"} or group_name == "seed_reference":
|
||
file_type = "family_pep"
|
||
else:
|
||
file_type = "all_pep"
|
||
if file_type:
|
||
item["files"][file_type] = original_name
|
||
if row_family:
|
||
groups["gene_family"] = row_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")
|
||
if "_ref_" in key:
|
||
species = key.split("_ref_", 1)[0]
|
||
else:
|
||
species = key.split("_", 1)[0]
|
||
found = _resolve_path(artifact.get("full_path") or artifact.get("path") or "", input_dir, workspace=workspace)
|
||
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}.all_pep.fasta", f"{species}.pep.fasta", f"{species}.protein.fasta", f"{species}.faa",
|
||
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",
|
||
]
|
||
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 _find_species_family_pep(
|
||
workspace: Path,
|
||
input_dir: Path,
|
||
species: str,
|
||
family: str,
|
||
original_name: str = "",
|
||
) -> Optional[Path]:
|
||
"""Resolve per-species family protein FASTA for interspecies tree building."""
|
||
found = _find_species_file(input_dir, species, family, "family_pep", original_name)
|
||
if found:
|
||
return found
|
||
for candidate in (
|
||
workspace / "final" / f"{species}.{family}.final_members.pep.fa",
|
||
workspace / "final" / f"{species}.{family}.curated_members.pep.fa",
|
||
input_dir / "final" / f"{species}.{family}.final_members.pep.fa",
|
||
input_dir / f"{species}.{family}.pep.fa",
|
||
input_dir / f"{species}.{family}.pep.fasta",
|
||
):
|
||
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_text(source, "rt") 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)
|
||
scope_data = getattr(args, "scope_data", {}) or {}
|
||
if scope_data:
|
||
groups["gene_family"] = scope_data.get("family") or groups.get("gene_family") or family
|
||
_add_scope_species_to_groups(workspace, input_dir, groups, scope_data)
|
||
# Interspecies synteny compares target species against the reference side.
|
||
# Known-family seed references such as At.SPL are valid reference inputs even
|
||
# when config marks them as seed_reference instead of comparative/ref.
|
||
for species, item in (groups.get("seed_reference") or {}).items():
|
||
groups["ref"].setdefault(species, item)
|
||
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
|
||
scope_all_pep = _scope_file_record(workspace, input_dir, scope_data, species, "pep") if scope_data else None
|
||
scope_family_pep = _scope_file_record(workspace, input_dir, scope_data, species, "family_pep") if scope_data else None
|
||
scope_gff = _scope_file_record(workspace, input_dir, scope_data, species, "gff") if scope_data else None
|
||
prepared_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["all_pep"] = prepared_all_pep or scope_all_pep
|
||
files["family_pep"] = (
|
||
scope_family_pep or
|
||
_find_species_family_pep(workspace, input_dir, species, item_family, 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", "")) or
|
||
scope_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 None
|
||
if raw_target_family_protein is None:
|
||
target_peps = [item.get("resolved", {}).get("family_pep") for item in groups["target"].values() if item.get("resolved", {}).get("family_pep")]
|
||
if target_peps:
|
||
raw_target_family_protein = analysis_dir / f"Target_selected.{family}.pep.fasta"
|
||
merge_many_fastas(target_peps, raw_target_family_protein)
|
||
else:
|
||
raw_target_family_protein = 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():
|
||
selected_target_output = analysis_dir / f"{prefix}.selected_target.family.pep.fa"
|
||
prefix_species = prefix.split(".", 1)[0] if "." in prefix else ""
|
||
single_species_family_fasta = (
|
||
len(selected_targets) == 1
|
||
and prefix_species == selected_targets[0]
|
||
and raw_target_family_protein.name == f"{prefix}.final_members.pep.fa"
|
||
)
|
||
if single_species_family_fasta:
|
||
# Species-specific Identify output may use transcript IDs while the
|
||
# standardized whole-proteome ID set uses protein IDs. It is already
|
||
# the selected species, so do not split it again by species IDs.
|
||
shutil.copy2(raw_target_family_protein, selected_target_output)
|
||
target_family_protein = selected_target_output
|
||
else:
|
||
target_family_protein = _filter_fasta_by_species(
|
||
raw_target_family_protein,
|
||
selected_target_output,
|
||
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 is None and raw_target_family_protein and raw_target_family_protein.exists() and selected_refs:
|
||
# Comparative species often have no curated family FASTA. If identify was
|
||
# run on their whole proteomes, split their family members from the same
|
||
# final_members file used for targets.
|
||
try:
|
||
ref_family_protein = _filter_fasta_by_species(
|
||
raw_target_family_protein,
|
||
analysis_dir / f"{prefix}.selected_ref.family.pep.fa",
|
||
selected_refs,
|
||
species_id_sets,
|
||
)
|
||
except Exception as exc:
|
||
print(f"[WARN] Could not derive reference family proteins from final_members: {exc}")
|
||
ref_family_protein = None
|
||
selected_ref_output = analysis_dir / f"{prefix}.selected_ref.family.pep.fa"
|
||
if ref_family_protein and Path(ref_family_protein).exists() and Path(ref_family_protein).resolve() != selected_ref_output.resolve():
|
||
ref_family_protein = _filter_fasta_by_species(
|
||
Path(ref_family_protein),
|
||
selected_ref_output,
|
||
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,
|
||
"scope": scope_data,
|
||
}
|
||
|
||
|
||
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_text(fasta, "rt") 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 _is_primary_chromosome_name(chrom: str) -> bool:
|
||
name = str(chrom or "").strip()
|
||
if re.fullmatch(r"\d+", name):
|
||
return True
|
||
return bool(re.fullmatch(r"chr(?:omosome)?[\s_.-]*\d+[A-Za-z]?", name, flags=re.IGNORECASE))
|
||
|
||
|
||
def _is_obvious_non_nuclear_sequence(chrom: str) -> bool:
|
||
name = str(chrom or "").strip().lower()
|
||
if re.search(r"(chloroplast|plastid|mitochond|mitoch|scaffold|contig|unplaced|unlocalized|random)", name):
|
||
return True
|
||
return bool(re.fullmatch(r"(chrm|chrmt|mt|mito|chrc|chrpt|pt|cp)", name))
|
||
|
||
|
||
def _select_primary_chromosomes_from_lengths(chrom_lengths: Dict[str, int]) -> Tuple[set, Dict[str, Any]]:
|
||
"""Select nuclear chromosome-scale records while falling back on unusual assemblies."""
|
||
records = [(chrom, int(length)) for chrom, length in chrom_lengths.items() if chrom and int(length) > 0]
|
||
total = len(records)
|
||
if not records:
|
||
return set(), {
|
||
"chromosomes_total": 0,
|
||
"chromosomes_plotted": 0,
|
||
"chromosomes_filtered": 0,
|
||
"chromosomes_filter_mode": "empty",
|
||
"chromosomes_filtered_names": [],
|
||
}
|
||
|
||
selected: List[Tuple[str, int]]
|
||
mode: str
|
||
named_primary = [(chrom, length) for chrom, length in records if _is_primary_chromosome_name(chrom)]
|
||
if named_primary:
|
||
selected = named_primary
|
||
mode = "primary_name"
|
||
else:
|
||
nuclear_candidates = [
|
||
(chrom, length)
|
||
for chrom, length in records
|
||
if not _is_obvious_non_nuclear_sequence(chrom)
|
||
]
|
||
length_source = nuclear_candidates or records
|
||
sorted_lengths = sorted((length for _, length in length_source), reverse=True)
|
||
median_length = sorted_lengths[len(sorted_lengths) // 2]
|
||
max_length = sorted_lengths[0]
|
||
threshold = max(1_000_000, int(median_length * 0.10), int(max_length * 0.02))
|
||
selected = [(chrom, length) for chrom, length in length_source if length >= threshold]
|
||
mode = "length_distribution"
|
||
|
||
if not selected:
|
||
selected = records
|
||
mode = "fallback_all_records"
|
||
|
||
selected_names = {chrom for chrom, _ in selected}
|
||
filtered_names = sorted(chrom for chrom, _ in records if chrom not in selected_names)
|
||
return selected_names, {
|
||
"chromosomes_total": total,
|
||
"chromosomes_plotted": len(selected_names),
|
||
"chromosomes_filtered": max(0, total - len(selected_names)),
|
||
"chromosomes_filter_mode": mode,
|
||
"chromosomes_filtered_names": filtered_names,
|
||
}
|
||
|
||
|
||
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 filtered 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()
|
||
records: List[Tuple[str, str, str, str]] = []
|
||
chrom_lengths: Dict[str, int] = {}
|
||
features_seen = 0
|
||
work_dir.mkdir(parents=True, exist_ok=True)
|
||
with open_text(gff3, "rt") as f_in:
|
||
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]
|
||
try:
|
||
chrom_lengths[chrom] = max(chrom_lengths.get(chrom, 0), int(float(end)))
|
||
except (TypeError, ValueError):
|
||
pass
|
||
records.append((chrom, start, end, gene_id))
|
||
|
||
selected_chroms, filter_metrics = _select_primary_chromosomes_from_lengths(chrom_lengths)
|
||
if not selected_chroms:
|
||
selected_chroms = set(chrom_lengths)
|
||
filtered_records = [record for record in records if record[0] in selected_chroms]
|
||
if not filtered_records and records:
|
||
filtered_records = records
|
||
filter_metrics["chromosomes_filter_mode"] = "fallback_all_records"
|
||
filter_metrics["chromosomes_plotted"] = filter_metrics["chromosomes_total"]
|
||
filter_metrics["chromosomes_filtered"] = 0
|
||
filter_metrics["chromosomes_filtered_names"] = []
|
||
|
||
with open(simple_gff, "w") as gff_out, open(bed_file, "w") as bed_out:
|
||
for chrom, start, end, gene_id in filtered_records:
|
||
gff_out.write(f"{chrom} {gene_id} {start} {end}\n")
|
||
bed_out.write(f"{chrom} {start} {end} {gene_id}\n")
|
||
|
||
count = len(filtered_records)
|
||
filter_metrics["genes_before_chromosome_filter"] = len(records)
|
||
filter_metrics["genes_after_chromosome_filter"] = count
|
||
filter_metrics["genes_filtered_by_chromosome"] = max(0, len(records) - count)
|
||
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 "",
|
||
**filter_metrics,
|
||
}
|
||
|
||
|
||
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,
|
||
scope_data: Optional[Dict[str, Any]] = None, scope_file: str = "",
|
||
analysis_route: str = "", member_set_id: str = "") -> Path:
|
||
end_time = datetime.now()
|
||
scope_data = scope_data or {}
|
||
parameters = {
|
||
"prefix": prefix,
|
||
"error": error,
|
||
"scope_file": scope_file,
|
||
"scope_id": scope_data.get("scope_id", ""),
|
||
"scope_mode": scope_data.get("scope_mode", ""),
|
||
"target_species": scope_data.get("target_species", []),
|
||
"reference_species": scope_data.get("reference_species", []),
|
||
"comparative_species": scope_data.get("comparative_species", []),
|
||
"selected_species": scope_data.get("selected_species", []),
|
||
"family": scope_data.get("family", ""),
|
||
}
|
||
if scope_data:
|
||
route = analysis_route or infer_analysis_route("interspecies_synteny", scope_data)
|
||
parameters["analysis_route"] = route
|
||
parameters["member_set_id"] = member_set_id or member_set_id_for(prefix, scope_data, route)
|
||
return write_manifest(
|
||
workspace, job_id, run_id, prefix, start_time, end_time,
|
||
artifacts={}, metrics={"error": error}, parameters=parameters,
|
||
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("--blast-engine", choices=["diamond", "blastp"], default="diamond", help="protein search engine for MCScanX pairwise BLAST file")
|
||
parser.add_argument("--plot-mode", choices=["multi", "pairwise", "both"], default="both", help="synteny plot outputs to generate: multi species, pairwise, or both")
|
||
parser.add_argument("--blast-evalue", default="1e-10", help="evalue cutoff for diamond/blastp search")
|
||
parser.add_argument("--max-target-seqs", type=int, default=5, help="maximum target sequences per query")
|
||
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?")
|
||
parser.add_argument("--reuse-tree", action="store_true", help="Reuse existing final interspecies tree if present")
|
||
parser.add_argument("--scope-file", help="Gene-family scope JSON from resolve_scope.py")
|
||
parser.add_argument("--analysis-route", help="Scoped analysis route recorded in manifest")
|
||
parser.add_argument("--member-set-id", help="Scoped member set identity recorded in manifest")
|
||
parser.add_argument("--tree-subfamily-arc-enabled", choices=["true", "false"], default="true",
|
||
help="Enable dynamic subfamily outer arcs in circular interspecies tree plots")
|
||
parser.add_argument("--tree-subfamily-arc-gap", type=float, default=1.2,
|
||
help="Extra gap between gene labels and circular subfamily arcs")
|
||
parser.add_argument("--tree-subfamily-arc-label-gap", type=float, default=0.6,
|
||
help="Extra gap between circular subfamily arcs and labels")
|
||
parser.add_argument("--tree-subfamily-arc-label-size", type=float, default=0.0,
|
||
help="Subfamily arc label size; 0 means automatic")
|
||
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
|
||
scope_data: Dict[str, Any] = {}
|
||
scope_file = ""
|
||
analysis_route = str(args.analysis_route or "").strip()
|
||
member_set_id = str(args.member_set_id or "").strip()
|
||
|
||
print("=" * 60)
|
||
print("Interspecies synteny analysis")
|
||
print("=" * 60)
|
||
|
||
try:
|
||
if args.scope_file:
|
||
scope_path = Path(args.scope_file)
|
||
if not scope_path.exists():
|
||
raise FileNotFoundError(f"Scope file not found: {args.scope_file}")
|
||
scope_file = str(scope_path.absolute())
|
||
scope_data = read_scope(scope_path)
|
||
args.scope_data = scope_data
|
||
_apply_scope_selection_to_args(args, scope_data)
|
||
if not analysis_route:
|
||
analysis_route = infer_analysis_route("interspecies_synteny", scope_data)
|
||
if not member_set_id and scope_data.get("scope_id") and analysis_route:
|
||
member_set_id = member_set_id_for(args.prefix, scope_data, analysis_route)
|
||
else:
|
||
args.scope_data = {}
|
||
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": "",
|
||
"blast_engine_requested": args.blast_engine,
|
||
"plot_mode": args.plot_mode,
|
||
"multi_plot_generated": False,
|
||
"pairwise_plot_count": 0,
|
||
"pairwise_plot_success": False,
|
||
"blast_engine_used": [],
|
||
"blast_fallback_count": 0,
|
||
}
|
||
|
||
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")
|
||
metrics["tree_merged_sequence_count"] = (merge_result.get("target_count") or 0) + (merge_result.get("ref_count") or 0)
|
||
ref_family_sources = [
|
||
str(item.get("resolved", {}).get("family_pep"))
|
||
for item in discovery["groups"].get("ref", {}).values()
|
||
if item.get("resolved", {}).get("family_pep")
|
||
]
|
||
metrics["tree_ref_family_sources"] = ref_family_sources
|
||
|
||
print("\nBuilding interspecies tree...")
|
||
final_tree = paths["final_dir"] / f"{args.prefix}.interspecies_tree.nwk"
|
||
reuse_tree = bool(args.reuse_tree and final_tree.exists())
|
||
if reuse_tree:
|
||
merged_tip_count = sum(1 for line in open_text(merged_fasta, "rt") if line.startswith(">"))
|
||
existing_tip_count = len(_read_tree_tip_labels(final_tree))
|
||
if existing_tip_count != merged_tip_count:
|
||
print(
|
||
f"Existing tree tip count ({existing_tip_count}) != merged family count ({merged_tip_count}); rebuilding tree"
|
||
)
|
||
reuse_tree = False
|
||
if reuse_tree:
|
||
print(f"Reusing existing interspecies tree: {final_tree}")
|
||
tree_result = {"tree_file": final_tree}
|
||
else:
|
||
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')}")
|
||
if Path(tree_result["tree_file"]).resolve() != final_tree.resolve():
|
||
shutil.copy(Path(tree_result["tree_file"]), final_tree)
|
||
metrics["tree_success"] = True
|
||
if tree_result.get("iqtree_executable"):
|
||
metrics["iqtree_executable"] = tree_result.get("iqtree_executable")
|
||
if tree_result.get("command"):
|
||
metrics["iqtree_command"] = tree_result.get("command")
|
||
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
|
||
|
||
subfamily_candidates: List[Path] = [
|
||
target_subfamily_file,
|
||
ref_subfamily_file,
|
||
input_dir / f"{args.prefix}.target_plus_reference.geneID.subfamily",
|
||
workspace / "final" / f"{args.prefix}.target_plus_reference.geneID.subfamily",
|
||
]
|
||
for species_value in list(discovery.get("target_species") or []) + list(discovery.get("ref_species") or []):
|
||
species_stem = str(species_value).strip()
|
||
if not species_stem:
|
||
continue
|
||
subfamily_candidates.extend([
|
||
input_dir / f"{species_stem}.SPL.target_plus_reference.geneID.subfamily",
|
||
workspace / "final" / f"{species_stem}.SPL.target_plus_reference.geneID.subfamily",
|
||
input_dir / f"{species_stem}.SPL.geneID.subfamily",
|
||
workspace / "final" / f"{species_stem}.SPL.geneID.subfamily",
|
||
])
|
||
subfamily_candidates.extend(sorted((workspace / "final").glob(f"{species_stem}.*.geneID.subfamily"), key=lambda p: p.stat().st_mtime, reverse=True))
|
||
merged_subfamily_file = _merge_matching_subfamily_files(subfamily_candidates, final_tree, paths["analysis_dir"], args.prefix)
|
||
if merged_subfamily_file:
|
||
target_subfamily_file = merged_subfamily_file
|
||
ref_subfamily_file = None
|
||
|
||
tree_plot_result = plot_interspecies_tree(
|
||
final_tree,
|
||
paths["final_dir"],
|
||
args.prefix,
|
||
target_group,
|
||
ref_group,
|
||
target_subfamily_file if target_subfamily_file and target_subfamily_file.exists() else None,
|
||
ref_subfamily_file if ref_subfamily_file and ref_subfamily_file.exists() else None,
|
||
gene_species_map,
|
||
subfamily_arc_enabled=str(args.tree_subfamily_arc_enabled).lower() != "false",
|
||
subfamily_arc_gap=args.tree_subfamily_arc_gap,
|
||
subfamily_arc_label_gap=args.tree_subfamily_arc_label_gap,
|
||
subfamily_arc_label_size=args.tree_subfamily_arc_label_size,
|
||
)
|
||
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")
|
||
if tree_plot_result.get("circular_pdf_file"):
|
||
add_artifact(artifacts, "tree_circular_pdf", tree_plot_result["circular_pdf_file"], f"final/{args.prefix}.interspecies_tree.circular_tree.pdf", "tree_circular_pdf", "plot")
|
||
if tree_plot_result.get("circular_png_file"):
|
||
add_artifact(artifacts, "tree_circular_png", tree_plot_result["circular_png_file"], f"final/{args.prefix}.interspecies_tree.circular_tree.png", "tree_circular_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 after chromosome filter "
|
||
f"({result.get('chromosomes_plotted', 'N/A')}/{result.get('chromosomes_total', 'N/A')} seqids plotted; "
|
||
f"filtered genes={result.get('genes_filtered_by_chromosome', 0)})"
|
||
)
|
||
|
||
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, args.blast_engine, args.blast_evalue, args.max_target_seqs)
|
||
pair["blast_success"] = bool(blast_result.get("success"))
|
||
pair["blast_engine_requested"] = blast_result.get("engine_requested", args.blast_engine)
|
||
pair["blast_engine_used"] = blast_result.get("engine_used")
|
||
pair["fallback_from"] = blast_result.get("fallback_from")
|
||
pair["blast_commands"] = blast_result.get("commands", [])
|
||
pair["blast_stderr"] = blast_result.get("stderr", "")
|
||
pair["diamond_stderr"] = blast_result.get("diamond_stderr", "")
|
||
pair["blast_output_stats"] = blast_result.get("output_stats", {})
|
||
if blast_result.get("engine_used"):
|
||
metrics.setdefault("blast_engine_used", []).append(blast_result.get("engine_used"))
|
||
if blast_result.get("fallback_from"):
|
||
metrics["blast_fallback_count"] = metrics.get("blast_fallback_count", 0) + 1
|
||
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")
|
||
successful_pairs = [p for p in pair_results if p.get("success")]
|
||
pairwise_plot_result = {"success": False, "pairs": [], "pairwise_plot_count": 0}
|
||
if args.plot_mode in {"multi", "both"}:
|
||
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,
|
||
)
|
||
if args.plot_mode in {"pairwise", "both"}:
|
||
pairwise_plot_result = plot_pairwise_synteny(
|
||
paths["analysis_dir"], paths["final_dir"], args.prefix, pair_results, family_ids,
|
||
)
|
||
metrics["multi_plot_generated"] = bool(synteny_plot_result.get("success"))
|
||
metrics["pairwise_plot_success"] = bool(pairwise_plot_result.get("success"))
|
||
metrics["pairwise_plot_count"] = int(pairwise_plot_result.get("pairwise_plot_count") or 0)
|
||
metrics["pairwise_plots"] = pairwise_plot_result.get("pairs", [])
|
||
requested_plot_success = True
|
||
if args.plot_mode in {"multi", "both"}:
|
||
requested_plot_success = requested_plot_success and bool(synteny_plot_result.get("success"))
|
||
if args.plot_mode in {"pairwise", "both"}:
|
||
requested_plot_success = requested_plot_success and bool(pairwise_plot_result.get("success"))
|
||
metrics["synteny_success"] = requested_plot_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 pairwise_plot_result.get("error") or "No target-ref pair produced JCVI anchors"
|
||
elif not requested_plot_success:
|
||
metrics["synteny_skipped_reason"] = synteny_plot_result.get("error") or pairwise_plot_result.get("error") or "Requested synteny plot output was not generated"
|
||
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")
|
||
for pair_plot in (metrics.get("pairwise_plots") or []):
|
||
pair_key = sanitize_token(pair_plot.get("pair") or "pair")
|
||
for artifact_key, file_key, artifact_name, file_type in [
|
||
("pairwise_pdf", "pdf_file", "pairwise_synteny_pdf", "plot"),
|
||
("pairwise_png", "png_file", "pairwise_synteny_png", "plot"),
|
||
("pairwise_marked_pdf", "marked_pdf_file", "pairwise_synteny_marked_pdf", "plot"),
|
||
("pairwise_marked_png", "marked_png_file", "pairwise_synteny_marked_png", "plot"),
|
||
("pairwise_seqids", "seqids_file", "pairwise_jcvi_seqids", "data"),
|
||
("pairwise_layout", "layout_file", "pairwise_jcvi_layout", "data"),
|
||
("pairwise_layout_mark", "layout_mark_file", "pairwise_jcvi_layout_mark", "data"),
|
||
]:
|
||
file_text = pair_plot.get(file_key) or ""
|
||
if not file_text:
|
||
continue
|
||
file_path = Path(file_text)
|
||
if not file_path.exists():
|
||
continue
|
||
rel = f"final/{file_path.name}" if file_path.parent == paths["final_dir"] else f"interspecies_synteny/analysis/{file_path.name}"
|
||
add_artifact(artifacts, f"{pair_key}_{artifact_key}", file_path, rel, artifact_name, file_type)
|
||
|
||
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,
|
||
"blast_engine_requested": args.blast_engine,
|
||
"plot_mode": args.plot_mode,
|
||
"blast_evalue": args.blast_evalue,
|
||
"max_target_seqs": args.max_target_seqs,
|
||
"tree_subfamily_arc_enabled": str(args.tree_subfamily_arc_enabled).lower() != "false",
|
||
"tree_subfamily_arc_gap": args.tree_subfamily_arc_gap,
|
||
"tree_subfamily_arc_label_gap": args.tree_subfamily_arc_label_gap,
|
||
"tree_subfamily_arc_label_size": args.tree_subfamily_arc_label_size,
|
||
}
|
||
if scope_data:
|
||
parameters.update({
|
||
"scope_file": scope_file,
|
||
"scope_id": scope_data.get("scope_id", ""),
|
||
"scope_mode": scope_data.get("scope_mode", ""),
|
||
"target_species": scope_data.get("target_species", []),
|
||
"reference_species": scope_data.get("reference_species", []),
|
||
"comparative_species": scope_data.get("comparative_species", []),
|
||
"selected_species": scope_data.get("selected_species", []),
|
||
"family": scope_data.get("family", ""),
|
||
})
|
||
route = analysis_route or infer_analysis_route("interspecies_synteny", scope_data)
|
||
parameters["analysis_route"] = route
|
||
parameters["member_set_id"] = member_set_id or member_set_id_for(args.prefix, scope_data, route)
|
||
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),
|
||
scope_data=scope_data, scope_file=scope_file,
|
||
analysis_route=analysis_route, member_set_id=member_set_id,
|
||
)
|
||
except Exception as manifest_exc:
|
||
print(f"Failed to write failure manifest: {manifest_exc}")
|
||
return 1
|
||
|
||
|
||
if __name__ == "__main__":
|
||
sys.exit(main())
|