1041 lines
40 KiB
Python
Executable File
1041 lines
40 KiB
Python
Executable File
#!/usr/bin/env python3
|
||
"""
|
||
PEP阶段：蛋白特征分析（完整版）
|
||
|
||
输入：
|
||
- Identify阶段的蛋白序列文件 (final/{prefix}.final_members.pep.fa)
|
||
|
||
输出：
|
||
- final/{prefix}.pep_properties.pdf/png (理化性质图)
|
||
- final/{prefix}.pep_tmhmm.pdf/png (跨膜结构图)
|
||
- final/{prefix}.pep_subcellular.pdf/png (亚细胞定位图)
|
||
- final/{prefix}.pep_report.md (分析报告)
|
||
- final/{prefix}.pep.manifest.json (机器清单)
|
||
"""
|
||
|
||
import os
|
||
import argparse
|
||
import json
|
||
import hashlib
|
||
import subprocess
|
||
import sys
|
||
import time
|
||
import shutil
|
||
import re
|
||
from datetime import datetime
|
||
from pathlib import Path
|
||
from typing import Dict, List, Any, Optional
|
||
|
||
PROJECT_ROOT = Path(__file__).resolve().parents[2]
|
||
if str(PROJECT_ROOT) not in sys.path:
|
||
sys.path.insert(0, str(PROJECT_ROOT))
|
||
|
||
from gfa_core.gene_family_scope import read_scope, sanitize_token
|
||
|
||
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", "")
|
||
|
||
# Plant-mPLoc code is retained, but the default PEP subcellular localization
|
||
# method is WoLF PSORT because it preserves per-protein predictions reliably.
|
||
PLANT_MPLOC_ENABLED = False
|
||
|
||
|
||
def load_manifest(manifest_file: Path) -> Optional[Dict[str, Any]]:
|
||
try:
|
||
return json.loads(manifest_file.read_text(encoding="utf-8"))
|
||
except Exception:
|
||
return None
|
||
|
||
|
||
def count_fasta_records(fasta_file: Path) -> int:
|
||
count = 0
|
||
with open(fasta_file, encoding="utf-8", errors="ignore") as handle:
|
||
for line in handle:
|
||
if line.startswith(">"):
|
||
count += 1
|
||
return count
|
||
|
||
|
||
def ensure_protein_fasta_without_stop_markers(protein_fasta: Path) -> Path:
|
||
has_stop = False
|
||
with open(protein_fasta, encoding="utf-8", errors="ignore") as handle:
|
||
for line in handle:
|
||
if not line.startswith(">") and "*" in line:
|
||
has_stop = True
|
||
break
|
||
if not has_stop:
|
||
return protein_fasta
|
||
clean_name = re.sub(r"(\.fa(?:sta)?)$", r".clean\1", protein_fasta.name, flags=re.IGNORECASE)
|
||
if clean_name == protein_fasta.name:
|
||
clean_name = protein_fasta.name + ".clean.fa"
|
||
clean_fasta = protein_fasta.with_name(clean_name)
|
||
if clean_fasta.exists() and clean_fasta.stat().st_mtime >= protein_fasta.stat().st_mtime:
|
||
return clean_fasta
|
||
with open(protein_fasta, "r", encoding="utf-8", errors="ignore") as src, open(clean_fasta, "w", encoding="utf-8") as dst:
|
||
for line in src:
|
||
dst.write(line if line.startswith(">") else line.replace("*", ""))
|
||
print(f"[INFO] Created stop-marker-clean FASTA for PEP: {clean_fasta}")
|
||
return clean_fasta
|
||
|
||
|
||
def resolve_identify_protein_fasta_for_scope(workspace: Path, prefix: str, scope_data: Dict[str, Any]) -> Optional[Path]:
|
||
scope_mode = str((scope_data or {}).get("scope_mode") or "").strip()
|
||
if scope_mode not in {"", "target_only"}:
|
||
return None
|
||
final_dir = workspace / "final"
|
||
family = str((scope_data or {}).get("family") or (prefix.split(".")[-1] if "." in prefix else prefix)).strip()
|
||
targets = [str(x).strip() for x in (scope_data or {}).get("target_species") or [] if str(x).strip()]
|
||
candidates = [final_dir / f"{prefix}.final_members.pep.fa"]
|
||
if len(targets) == 1:
|
||
species = targets[0]
|
||
candidates = [
|
||
workspace / "identify" / species / "final" / f"{species}.final_gfmembers.pep.fa",
|
||
final_dir / f"Target.{species}.{family}.final_members.pep.fa",
|
||
final_dir / f"{species}.{family}.final_members.pep.fa",
|
||
*candidates,
|
||
]
|
||
for candidate in candidates:
|
||
if candidate.exists() and candidate.stat().st_size > 0:
|
||
return ensure_protein_fasta_without_stop_markers(candidate)
|
||
return None
|
||
|
||
|
||
def resolve_scope_protein_fasta(workspace: Path, prefix: str, scope_data: Dict[str, Any], cli_fasta: Optional[str] = None) -> Path:
|
||
if cli_fasta:
|
||
fasta_path = Path(cli_fasta)
|
||
if fasta_path.exists() and fasta_path.stat().st_size > 0:
|
||
return ensure_protein_fasta_without_stop_markers(fasta_path)
|
||
raise FileNotFoundError(f"Specified protein FASTA not found: {cli_fasta}")
|
||
|
||
final_dir = workspace / "final"
|
||
scope_id = str(scope_data.get("scope_id") or "")
|
||
manifest_candidates = []
|
||
if scope_id:
|
||
manifest_candidates.append(final_dir / f"{prefix}.{sanitize_token(scope_id)}.phylogeny.manifest.json")
|
||
manifest_candidates.append(final_dir / f"{prefix}.phylogeny.manifest.json")
|
||
|
||
for manifest_file in manifest_candidates:
|
||
if not manifest_file.exists():
|
||
continue
|
||
manifest_data = load_manifest(manifest_file)
|
||
if not manifest_data:
|
||
continue
|
||
fasta_text = manifest_data.get("parameters", {}).get("scope_input_fasta")
|
||
if fasta_text:
|
||
fasta_path = Path(fasta_text)
|
||
if fasta_path.exists() and fasta_path.stat().st_size > 0:
|
||
print(f"[INFO] Found scoped protein FASTA from phylogeny manifest: {fasta_path}")
|
||
return ensure_protein_fasta_without_stop_markers(fasta_path)
|
||
|
||
fallback = resolve_identify_protein_fasta_for_scope(workspace, prefix, scope_data)
|
||
if fallback:
|
||
print(f"[INFO] Using target-only Identify protein FASTA: {fallback}")
|
||
return fallback
|
||
|
||
raise FileNotFoundError(
|
||
"Cannot find scoped protein FASTA. Run phylogeny with the same scope first, "
|
||
"or provide --fasta explicitly."
|
||
)
|
||
|
||
|
||
def scope_parameters(scope_file: Optional[Path], scope_data: Optional[Dict[str, Any]]) -> Dict[str, Any]:
|
||
if not scope_file or not scope_data:
|
||
return {}
|
||
return {
|
||
"scope_file": str(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", ""),
|
||
}
|
||
|
||
|
||
def compute_file_hash(filepath: Path) -> str:
|
||
"""计算文件的SHA256哈希值"""
|
||
sha256_hash = hashlib.sha256()
|
||
try:
|
||
with open(filepath, "rb") as f:
|
||
for byte_block in iter(lambda: f.read(4096), b""):
|
||
sha256_hash.update(byte_block)
|
||
return sha256_hash.hexdigest()
|
||
except Exception:
|
||
return ""
|
||
|
||
|
||
def get_tool_versions() -> Dict[str, str]:
|
||
"""获取工具版本"""
|
||
versions = {}
|
||
|
||
try:
|
||
result = subprocess.run(["R", "--version"], capture_output=True, text=True, check=False, timeout=5)
|
||
versions["R"] = result.stdout.split('\n')[0] if result.stdout else "available"
|
||
except Exception:
|
||
versions["R"] = "not_found"
|
||
|
||
try:
|
||
biolib_bin = "/home/xzq/miniconda3/envs/test4py311/bin/biolib"
|
||
result = subprocess.run([biolib_bin, "--version"], capture_output=True, text=True, check=False, timeout=5)
|
||
versions["biolib"] = "available" if result.returncode == 0 else "not_found"
|
||
except Exception:
|
||
versions["biolib"] = "not_found"
|
||
|
||
return versions
|
||
|
||
|
||
def prepare_directories(workspace: Path) -> Dict[str, Path]:
|
||
"""创建目录"""
|
||
paths = {
|
||
"analysis_dir": workspace / "pep" / "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 parse_fasta(fasta_file: Path) -> Dict[str, str]:
|
||
"""解析FASTA文件"""
|
||
sequences = {}
|
||
current_id = None
|
||
current_seq = []
|
||
|
||
with open(fasta_file, 'r') as f:
|
||
for line in f:
|
||
line = line.strip()
|
||
if line.startswith('>'):
|
||
if current_id:
|
||
sequences[current_id] = ''.join(current_seq)
|
||
current_id = line[1:].split()[0]
|
||
current_seq = []
|
||
else:
|
||
current_seq.append(line)
|
||
|
||
if current_id:
|
||
sequences[current_id] = ''.join(current_seq)
|
||
|
||
return sequences
|
||
|
||
|
||
def external_tool_timeout(fasta_file: Path, seconds_per_sequence: int = 90, minimum: int = 300, maximum: int = 7200) -> int:
|
||
"""Estimate a conservative timeout for slow web/cloud predictors."""
|
||
try:
|
||
count = sum(1 for line in open(fasta_file, 'r') if line.startswith('>'))
|
||
except Exception:
|
||
count = 1
|
||
return max(minimum, min(maximum, count * seconds_per_sequence))
|
||
|
||
|
||
def write_clean_fasta(fasta_file: Path, output_file: Path) -> Path:
|
||
"""Write DeepTMHMM-friendly FASTA with short IDs and canonical amino acids."""
|
||
allowed = set("ACDEFGHIKLMNPQRSTVWY")
|
||
sequences = parse_fasta(fasta_file)
|
||
with open(output_file, "w") as f:
|
||
for seq_id, seq in sequences.items():
|
||
clean = "".join(aa for aa in seq.upper() if aa in allowed)
|
||
f.write(f">{seq_id}\n{clean}\n")
|
||
return output_file
|
||
|
||
|
||
def run_deeptmhmm_dtu_web(fasta_file: Path, output_dir: Path, timeout_sec: int) -> Dict[str, Any]:
|
||
"""Run DeepTMHMM through the DTU web form when BioLib file staging fails."""
|
||
from selenium import webdriver
|
||
from selenium.webdriver.common.by import By
|
||
from selenium.webdriver.firefox.options import Options
|
||
from selenium.webdriver.firefox.service import Service
|
||
from selenium.webdriver.support import expected_conditions as EC
|
||
from selenium.webdriver.support.ui import WebDriverWait
|
||
|
||
geckodriver = shutil.which("geckodriver")
|
||
if not geckodriver:
|
||
return {"success": False, "stderr": "geckodriver not found"}
|
||
|
||
clean_fasta = write_clean_fasta(fasta_file, output_dir / f"{fasta_file.stem}.deeptmhmm.clean.fa")
|
||
html_file = output_dir / "deeptmhmm_dtu.html"
|
||
topologies = output_dir / "predicted_topologies.3line"
|
||
gff_file = output_dir / "deeptmhmm_results.gff3"
|
||
|
||
options = Options()
|
||
options.add_argument("-headless")
|
||
driver = webdriver.Firefox(service=Service(geckodriver), options=options)
|
||
try:
|
||
driver.get("https://services.healthtech.dtu.dk/services/DeepTMHMM-1.0/")
|
||
WebDriverWait(driver, 30).until(EC.presence_of_element_located((By.NAME, "uploadfile")))
|
||
driver.find_element(By.NAME, "uploadfile").send_keys(str(clean_fasta.resolve()))
|
||
driver.find_element(By.CSS_SELECTOR, 'input[type="submit"]').click()
|
||
start = time.time()
|
||
while time.time() - start < timeout_sec:
|
||
body = driver.find_element(By.TAG_NAME, "body").text
|
||
if "Predicted Topologies" in body or "Predicted topologies" in body:
|
||
html_file.write_text(driver.page_source, encoding="utf-8")
|
||
lines = body.splitlines()
|
||
start_idx = None
|
||
for i, line in enumerate(lines):
|
||
if line.strip() == "Predicted Topologies":
|
||
start_idx = i + 1
|
||
break
|
||
if start_idx is None:
|
||
return {"success": False, "stderr": "Predicted Topologies section not found"}
|
||
topo_lines = []
|
||
gff_lines = []
|
||
in_gff = False
|
||
for line in lines[start_idx:]:
|
||
if line.startswith("##gff-version"):
|
||
in_gff = True
|
||
if line.startswith("DTU Health Tech") or line.startswith("ADDRESS & INFO"):
|
||
break
|
||
if in_gff:
|
||
if line.strip():
|
||
gff_lines.append(line)
|
||
else:
|
||
if line.strip():
|
||
topo_lines.append(line)
|
||
topologies.write_text("\n".join(topo_lines).strip() + "\n", encoding="utf-8")
|
||
gff_file.write_text("\n".join(gff_lines).strip() + "\n", encoding="utf-8")
|
||
return {
|
||
"success": topologies.exists() and topologies.stat().st_size > 0,
|
||
"results_md": html_file,
|
||
"topologies": topologies,
|
||
"gff": gff_file,
|
||
"method": "DTU DeepTMHMM web",
|
||
}
|
||
if "failed" in body.lower() or "error" in body.lower():
|
||
html_file.write_text(driver.page_source, encoding="utf-8")
|
||
return {"success": False, "stderr": body[:500]}
|
||
time.sleep(20)
|
||
html_file.write_text(driver.page_source, encoding="utf-8")
|
||
return {"success": False, "stderr": "DTU DeepTMHMM web job timed out"}
|
||
finally:
|
||
driver.quit()
|
||
|
||
|
||
def calculate_physicochemical_properties(sequences: Dict[str, str]) -> List[Dict]:
|
||
"""计算理化性质"""
|
||
print(f" 计算理化性质...")
|
||
|
||
# 氨基酸分子量
|
||
aa_weights = {
|
||
'A': 89.09, 'R': 174.20, 'N': 132.12, 'D': 133.10, 'C': 121.16,
|
||
'E': 147.13, 'Q': 146.15, 'G': 75.07, 'H': 155.16, 'I': 131.17,
|
||
'L': 131.17, 'K': 146.19, 'M': 149.21, 'F': 165.19, 'P': 115.13,
|
||
'S': 105.09, 'T': 119.12, 'W': 204.23, 'Y': 181.19, 'V': 117.15
|
||
}
|
||
|
||
# 氨基酸电荷
|
||
aa_charge = {
|
||
'D': -1, 'E': -1, 'K': 1, 'R': 1, 'H': 0.5
|
||
}
|
||
|
||
properties = []
|
||
|
||
for seq_id, sequence in sequences.items():
|
||
seq = sequence.upper()
|
||
length = len(seq)
|
||
|
||
# 分子量
|
||
mw = sum(aa_weights.get(aa, 110) for aa in seq) - (length - 1) * 18.02
|
||
|
||
# 等电点估算
|
||
acidic = seq.count('D') + seq.count('E')
|
||
basic = seq.count('K') + seq.count('R') + seq.count('H')
|
||
if basic > acidic:
|
||
pI = 7.0 + (basic - acidic) * 0.5
|
||
elif acidic > basic:
|
||
pI = 7.0 - (acidic - basic) * 0.3
|
||
else:
|
||
pI = 7.0
|
||
|
||
# 疏水性 (简化的GRAVY)
|
||
hydrophobic_aa = 'AILMFWV'
|
||
hydrophilic_aa = 'RNDQEK'
|
||
hydrophobic_count = sum(seq.count(aa) for aa in hydrophobic_aa)
|
||
hydrophilic_count = sum(seq.count(aa) for aa in hydrophilic_aa)
|
||
gravy = (hydrophobic_count - hydrophilic_count) / length if length > 0 else 0
|
||
|
||
# 脂肪族指数
|
||
aliphatic = (seq.count('A') + seq.count('V') + seq.count('I') + seq.count('L')) / length * 100
|
||
|
||
properties.append({
|
||
"id": seq_id,
|
||
"length": length,
|
||
"molecular_weight": mw,
|
||
"isoelectric_point": pI,
|
||
"gravy": gravy,
|
||
"aliphatic_index": aliphatic
|
||
})
|
||
|
||
return properties
|
||
|
||
|
||
def run_deeptmhmm(fasta_file: Path, output_dir: Path) -> Dict[str, Any]:
|
||
"""运行DeepTMHMM跨膜预测"""
|
||
print(f" 运行DeepTMHMM...")
|
||
|
||
results_md = output_dir / "deeptmhmm_results.md"
|
||
topologies = output_dir / "predicted_topologies.3line"
|
||
for old_file in (results_md, topologies):
|
||
if old_file.exists():
|
||
old_file.unlink()
|
||
timeout_sec = external_tool_timeout(fasta_file, seconds_per_sequence=120, minimum=900, maximum=10800)
|
||
|
||
web_result = run_deeptmhmm_dtu_web(fasta_file, output_dir, timeout_sec)
|
||
if web_result.get("success"):
|
||
return web_result
|
||
print(f" DTU网页DeepTMHMM失败，尝试BioLib: {web_result.get('stderr', '')[:200]}")
|
||
|
||
try:
|
||
biolib_bin = "/home/xzq/miniconda3/envs/test4py311/bin/biolib"
|
||
cmd = [biolib_bin, "run", "DTU/DeepTMHMM", "--fasta", str(fasta_file)]
|
||
result = subprocess.run(cmd, capture_output=True, text=True,
|
||
cwd=str(output_dir), timeout=timeout_sec)
|
||
|
||
if result.returncode != 0 or not topologies.exists():
|
||
stderr = result.stderr[:300] if result.stderr else "DeepTMHMM did not generate topology output"
|
||
print(f" DeepTMHMM未生成有效结果: {stderr}")
|
||
print(" 尝试通过DTU官方网页运行DeepTMHMM...")
|
||
return run_deeptmhmm_dtu_web(fasta_file, output_dir, timeout_sec)
|
||
|
||
return {
|
||
"success": True,
|
||
"results_md": results_md,
|
||
"topologies": topologies
|
||
}
|
||
|
||
except Exception as e:
|
||
print(f" DeepTMHMM错误: {e}")
|
||
return {"success": False, "stderr": str(e)}
|
||
|
||
|
||
def run_plant_mploc(fasta_file: Path, output_file: Path) -> Dict[str, Any]:
|
||
"""运行Plant-mPLoc"""
|
||
print(f" 运行Plant-mPLoc...")
|
||
|
||
try:
|
||
script = Path("scripts/get_Plant_mPLoc_out.py")
|
||
|
||
if script.exists():
|
||
cmd = [sys.executable, str(script), "--input", str(fasta_file), "--output", str(output_file)]
|
||
if output_file.exists():
|
||
output_file.unlink()
|
||
started_at = time.time()
|
||
timeout_sec = external_tool_timeout(fasta_file, seconds_per_sequence=90)
|
||
result = subprocess.run(cmd, capture_output=True, text=True, timeout=timeout_sec)
|
||
if result.returncode != 0:
|
||
return {"success": False, "stderr": result.stderr or result.stdout}
|
||
else:
|
||
return {"success": False, "stderr": f"Plant-mPLoc script not found: {script}"}
|
||
|
||
if not output_file.exists() or output_file.stat().st_size == 0 or output_file.stat().st_mtime < started_at:
|
||
return {"success": False, "stderr": "Plant-mPLoc did not generate output"}
|
||
|
||
return {"success": True, "output_file": output_file, "method": "Plant-mPLoc website automation"}
|
||
|
||
except Exception as e:
|
||
print(f" Plant-mPLoc错误: {e}")
|
||
return {"success": False, "stderr": str(e)}
|
||
|
||
|
||
def run_wolfpsort(fasta_file: Path, output_file: Path) -> Dict[str, Any]:
|
||
"""运行WoLF PSORT"""
|
||
print(f" 运行WoLF PSORT...")
|
||
|
||
try:
|
||
script = Path("scripts/wolfpsort_batch.py")
|
||
|
||
if script.exists():
|
||
cmd = [sys.executable, str(script), "--input", str(fasta_file), "--output", str(output_file)]
|
||
if output_file.exists():
|
||
output_file.unlink()
|
||
started_at = time.time()
|
||
timeout_sec = external_tool_timeout(fasta_file, seconds_per_sequence=90)
|
||
result = subprocess.run(cmd, capture_output=True, text=True, timeout=timeout_sec)
|
||
if result.returncode != 0:
|
||
return {"success": False, "stderr": result.stderr or result.stdout}
|
||
else:
|
||
return {"success": False, "stderr": f"WoLF PSORT script not found: {script}"}
|
||
|
||
if not output_file.exists() or output_file.stat().st_size == 0 or output_file.stat().st_mtime < started_at:
|
||
return {"success": False, "stderr": "WoLF PSORT did not generate output"}
|
||
|
||
return {"success": True, "output_file": output_file, "method": "WoLF PSORT website automation"}
|
||
|
||
except Exception as e:
|
||
print(f" WoLF PSORT错误: {e}")
|
||
return {"success": False, "stderr": str(e)}
|
||
|
||
|
||
def plot_properties_r(properties: List[Dict], output_pdf: Path, output_png: Path) -> bool:
|
||
"""使用R绘制理化性质图"""
|
||
print(f" 绘制理化性质图...")
|
||
|
||
# 准备数据文件
|
||
data_file = output_pdf.parent / "properties_data.tsv"
|
||
with open(data_file, 'w') as f:
|
||
f.write("id\tlength\tmw\tpI\tgravy\taliphatic\n")
|
||
for p in properties:
|
||
f.write(f"{p['id']}\t{p['length']}\t{p['molecular_weight']:.2f}\t"
|
||
f"{p['isoelectric_point']:.2f}\t{p['gravy']:.3f}\t{p['aliphatic_index']:.2f}\n")
|
||
|
||
r_script = f'''
|
||
library(ggplot2)
|
||
library(patchwork)
|
||
|
||
data <- read.delim("{data_file}", sep="\t", header=TRUE)
|
||
|
||
# 分子量分布
|
||
p1 <- ggplot(data, aes(x=mw/1000)) +
|
||
geom_histogram(bins=15, fill="steelblue", color="white") +
|
||
labs(title="Molecular Weight Distribution",
|
||
x="Molecular Weight (kDa)", y="Count") +
|
||
theme_minimal()
|
||
|
||
# 等电点分布
|
||
p2 <- ggplot(data, aes(x=pI)) +
|
||
geom_histogram(bins=15, fill="coral", color="white") +
|
||
labs(title="Isoelectric Point Distribution", x="pI", y="Count") +
|
||
theme_minimal()
|
||
|
||
# 序列长度分布
|
||
p3 <- ggplot(data, aes(x=length)) +
|
||
geom_histogram(bins=15, fill="seagreen", color="white") +
|
||
labs(title="Sequence Length Distribution", x="Length (aa)", y="Count") +
|
||
theme_minimal()
|
||
|
||
# 疏水性分布
|
||
p4 <- ggplot(data, aes(x=gravy)) +
|
||
geom_histogram(bins=15, fill="purple", color="white") +
|
||
labs(title="GRAVY (Hydrophobicity)", x="GRAVY Score", y="Count") +
|
||
theme_minimal()
|
||
|
||
combined <- (p1 + p2) / (p3 + p4)
|
||
|
||
ggsave("{output_pdf}", combined, width=12, height=10)
|
||
|
||
png("{output_png}", width=1200, height=1000, res=150)
|
||
print(combined)
|
||
dev.off()
|
||
'''
|
||
|
||
r_file = output_pdf.parent / "plot_properties.R"
|
||
with open(r_file, 'w') as f:
|
||
f.write(r_script)
|
||
|
||
try:
|
||
rscript = "/home/xzq/miniconda3/envs/test4py311/bin/Rscript"
|
||
result = subprocess.run([rscript, str(r_file)],
|
||
capture_output=True, text=True, timeout=120)
|
||
if result.returncode != 0:
|
||
print(f" R绘图错误: {result.stderr[:300]}")
|
||
return output_pdf.exists() and output_png.exists()
|
||
except Exception as e:
|
||
print(f" R绘图错误: {e}")
|
||
return False
|
||
|
||
|
||
def plot_tmhmm_results(topologies_file: Path, output_pdf: Path, output_png: Path) -> bool:
|
||
"""绘制跨膜结构图"""
|
||
print(f" 绘制跨膜结构图...")
|
||
|
||
# 解析拓扑结构
|
||
transmembrane_count = {}
|
||
|
||
if topologies_file.exists():
|
||
with open(topologies_file, 'r') as f:
|
||
current_id = None
|
||
for line in f:
|
||
line = line.strip()
|
||
if line.startswith('>'):
|
||
current_id = line[1:]
|
||
elif current_id and line and not line.startswith('>'):
|
||
# 统计M（跨膜）的数量
|
||
tm_count = line.count('M')
|
||
transmembrane_count[current_id] = tm_count
|
||
current_id = None
|
||
|
||
# 准备数据
|
||
data_file = output_pdf.parent / "tmhmm_data.tsv"
|
||
with open(data_file, 'w') as f:
|
||
f.write("id\ttm_count\n")
|
||
for seq_id, count in transmembrane_count.items():
|
||
f.write(f"{seq_id}\t{count}\n")
|
||
|
||
r_script = f'''
|
||
library(ggplot2)
|
||
|
||
data <- read.delim("{data_file}", sep="\t", header=TRUE)
|
||
|
||
# 跨膜螺旋数量分布
|
||
p <- ggplot(data, aes(x=tm_count)) +
|
||
geom_bar(fill="darkblue", color="white") +
|
||
labs(title="Transmembrane Helices Distribution",
|
||
x="Number of TM Helices", y="Count") +
|
||
theme_minimal()
|
||
|
||
ggsave("{output_pdf}", p, width=8, height=6)
|
||
|
||
png("{output_png}", width=800, height=600, res=150)
|
||
print(p)
|
||
dev.off()
|
||
'''
|
||
|
||
r_file = output_pdf.parent / "plot_tmhmm.R"
|
||
with open(r_file, 'w') as f:
|
||
f.write(r_script)
|
||
|
||
try:
|
||
rscript = "/home/xzq/miniconda3/envs/test4py311/bin/Rscript"
|
||
result = subprocess.run([rscript, str(r_file)],
|
||
capture_output=True, text=True, timeout=120)
|
||
if result.returncode != 0:
|
||
print(f" R绘图错误: {result.stderr[:300]}")
|
||
return output_pdf.exists() and output_png.exists()
|
||
except Exception as e:
|
||
print(f" R绘图错误: {e}")
|
||
return False
|
||
|
||
|
||
|
||
LOCATION_LABELS = {
|
||
"nucl": "Nucleus",
|
||
"nucleus": "Nucleus",
|
||
"cyto": "Cytoplasm",
|
||
"cytoplasm": "Cytoplasm",
|
||
"chlo": "Chloroplast",
|
||
"chloroplast": "Chloroplast",
|
||
"mito": "Mitochondrion",
|
||
"mitochondrion": "Mitochondrion",
|
||
"plas": "Plasma membrane",
|
||
"plasma membrane": "Plasma membrane",
|
||
"vacu": "Vacuole",
|
||
"vacuole": "Vacuole",
|
||
"golg": "Golgi apparatus",
|
||
"golgi": "Golgi apparatus",
|
||
"extr": "Extracellular",
|
||
"extracellular": "Extracellular",
|
||
"pero": "Peroxisome",
|
||
"peroxisome": "Peroxisome",
|
||
"e.r.": "Endoplasmic reticulum",
|
||
"er": "Endoplasmic reticulum",
|
||
"cysk": "Cytoskeleton",
|
||
"cytoskeleton": "Cytoskeleton",
|
||
}
|
||
|
||
|
||
def normalize_location_label(label: str) -> str:
|
||
clean = re.sub(r"\s+", " ", str(label or "").strip())
|
||
if not clean:
|
||
return ""
|
||
key = clean.lower().replace("_", " ")
|
||
if key in LOCATION_LABELS:
|
||
return LOCATION_LABELS[key]
|
||
if " " not in clean and "_" in clean:
|
||
parts = [normalize_location_label(part) for part in clean.split("_")]
|
||
parts = [part for part in parts if part]
|
||
if parts:
|
||
return " / ".join(parts)
|
||
return clean[:1].upper() + clean[1:]
|
||
|
||
|
||
def parse_plant_mploc_counts(plant_mploc_file: Path) -> tuple[Dict[str, int], int]:
|
||
counts: Dict[str, int] = {}
|
||
rows = 0
|
||
if not plant_mploc_file.exists() or plant_mploc_file.stat().st_size == 0:
|
||
return counts, rows
|
||
with open(plant_mploc_file, "r", encoding="utf-8", errors="ignore") as f:
|
||
for line in f:
|
||
line = line.strip()
|
||
if not line:
|
||
continue
|
||
parts = line.split("\t")
|
||
if len(parts) < 2:
|
||
continue
|
||
first = parts[0].strip().lower()
|
||
second = parts[1].strip()
|
||
if first in {"protein_id", "protein", "id", "gene_id", "sequence_id"}:
|
||
continue
|
||
loc = normalize_location_label(second)
|
||
if not loc:
|
||
continue
|
||
counts[loc] = counts.get(loc, 0) + 1
|
||
rows += 1
|
||
return counts, rows
|
||
|
||
|
||
def parse_wolfpsort_counts(wolfpsort_file: Path) -> tuple[Dict[str, int], int]:
|
||
counts: Dict[str, int] = {}
|
||
rows = 0
|
||
if not wolfpsort_file.exists() or wolfpsort_file.stat().st_size == 0:
|
||
return counts, rows
|
||
with open(wolfpsort_file, "r", encoding="utf-8", errors="ignore") as f:
|
||
for line in f:
|
||
line = line.strip()
|
||
if not line:
|
||
continue
|
||
parts = line.split("\t")
|
||
if len(parts) < 2 or parts[0].strip().lower() in {"protein_id", "protein", "id"}:
|
||
continue
|
||
prediction = parts[1].strip()
|
||
if not prediction:
|
||
continue
|
||
primary = prediction.split(",", 1)[0]
|
||
primary = primary.split(":", 1)[0].strip()
|
||
loc = normalize_location_label(primary)
|
||
if not loc:
|
||
continue
|
||
counts[loc] = counts.get(loc, 0) + 1
|
||
rows += 1
|
||
return counts, rows
|
||
|
||
|
||
def resolve_subcellular_counts(plant_mploc_file: Path, wolfpsort_file: Path,
|
||
allow_plant_fallback: bool = False) -> tuple[Dict[str, int], str, int]:
|
||
wolf_counts, wolf_rows = parse_wolfpsort_counts(wolfpsort_file)
|
||
if wolf_rows:
|
||
return wolf_counts, "WoLF PSORT", wolf_rows
|
||
|
||
if allow_plant_fallback:
|
||
plant_counts, plant_rows = parse_plant_mploc_counts(plant_mploc_file)
|
||
if plant_rows:
|
||
return plant_counts, "Plant-mPLoc", plant_rows
|
||
|
||
return {}, "WoLF PSORT", 0
|
||
|
||
def plot_subcellular_results(plant_mploc_file: Path, wolfpsort_file: Path,
|
||
output_pdf: Path, output_png: Path) -> bool:
|
||
"""绘制亚细胞定位图"""
|
||
print(f" 绘制亚细胞定位图...")
|
||
|
||
location_counts, source, parsed_rows = resolve_subcellular_counts(plant_mploc_file, wolfpsort_file)
|
||
if not location_counts:
|
||
print(" 亚细胞定位结果为空，跳过绘图")
|
||
return False
|
||
|
||
# 准备数据
|
||
data_file = output_pdf.parent / "subcellular_data.tsv"
|
||
with open(data_file, 'w') as f:
|
||
f.write("location\tcount\tsource\n")
|
||
for loc, count in sorted(location_counts.items(), key=lambda item: item[1], reverse=True):
|
||
f.write(f"{loc}\t{count}\t{source}\n")
|
||
|
||
r_script = f'''
|
||
library(ggplot2)
|
||
|
||
data <- read.delim("{data_file}", sep="\t", header=TRUE)
|
||
|
||
# 按count排序
|
||
data$location <- factor(data$location, levels=data$location[order(data$count, decreasing=TRUE)])
|
||
plot_source <- unique(data$source)[1]
|
||
|
||
p <- ggplot(data, aes(x=location, y=count, fill=location)) +
|
||
geom_bar(stat="identity", width=0.55, color="white") +
|
||
geom_text(aes(label=count), vjust=-0.3, size=4) +
|
||
labs(title=paste0("Subcellular Localization (", plot_source, ")"),
|
||
x="Location", y="Count") +
|
||
theme_minimal() +
|
||
theme(axis.text.x=element_text(angle=45, hjust=1), legend.position="none")
|
||
|
||
ggsave("{output_pdf}", p, width=8, height=6)
|
||
|
||
png("{output_png}", width=800, height=600, res=150)
|
||
print(p)
|
||
dev.off()
|
||
'''
|
||
|
||
r_file = output_pdf.parent / "plot_subcellular.R"
|
||
with open(r_file, 'w') as f:
|
||
f.write(r_script)
|
||
|
||
try:
|
||
rscript = "/home/xzq/miniconda3/envs/test4py311/bin/Rscript"
|
||
result = subprocess.run([rscript, str(r_file)],
|
||
capture_output=True, text=True, timeout=120)
|
||
if result.returncode != 0:
|
||
print(f" R绘图错误: {result.stderr[:300]}")
|
||
return output_pdf.exists() and output_png.exists()
|
||
except Exception as e:
|
||
print(f" R绘图错误: {e}")
|
||
return False
|
||
|
||
def write_report(workspace: Path, prefix: str, metrics: Dict, parameters: Dict,
|
||
tool_versions: Dict) -> Path:
|
||
"""写入报告"""
|
||
report_path = workspace / "final" / f"{prefix}.pep_report.md"
|
||
|
||
with open(report_path, "w") as f:
|
||
f.write(f"# 蛋白特征分析报告 - {prefix}\n\n")
|
||
f.write(f"**生成时间**: {datetime.now().isoformat()}\n\n")
|
||
|
||
f.write("## 1. 分析概述\n")
|
||
f.write("- **分析阶段**: PEP (Protein Features)\n")
|
||
f.write("- **分析内容**: \n")
|
||
f.write(" - 理化性质 (分子量、等电点、疏水性)\n")
|
||
f.write(" - 跨膜结构预测 (DeepTMHMM)\n")
|
||
f.write(" - 亚细胞定位 (WoLF PSORT；Plant-mPLoc代码保留但默认关闭)\n\n")
|
||
|
||
if "protein_count" in metrics:
|
||
f.write(f"- **蛋白数量**: {metrics['protein_count']}\n\n")
|
||
|
||
subcellular_data = workspace / "final" / "subcellular_data.tsv"
|
||
subcellular_rows = []
|
||
if subcellular_data.exists() and subcellular_data.stat().st_size > 0:
|
||
with open(subcellular_data, "r", encoding="utf-8", errors="ignore") as sf:
|
||
header = sf.readline().strip().split("\t")
|
||
for line in sf:
|
||
parts = line.rstrip("\n").split("\t")
|
||
if len(parts) < 2:
|
||
continue
|
||
row = dict(zip(header, parts))
|
||
subcellular_rows.append(row)
|
||
if subcellular_rows:
|
||
f.write("## 2. 亚细胞定位数量\n")
|
||
f.write("| Location | Count | Source |\n")
|
||
f.write("|---|---:|---|\n")
|
||
for row in subcellular_rows:
|
||
f.write(
|
||
f"| {row.get('location', '')} | {row.get('count', '')} | "
|
||
f"{row.get('source', '')} |\n"
|
||
)
|
||
f.write("\n")
|
||
|
||
f.write("## 3. 工具版本\n")
|
||
for tool, version in tool_versions.items():
|
||
f.write(f"- **{tool}**: {version}\n")
|
||
f.write("\n")
|
||
|
||
f.write("## 4. 输出文件\n")
|
||
f.write(f"- `{prefix}.pep_properties.pdf/png` - 理化性质分布图\n")
|
||
f.write(f"- `{prefix}.pep_tmhmm.pdf/png` - 跨膜结构预测图\n")
|
||
f.write(f"- `{prefix}.pep_subcellular.pdf/png` - 亚细胞定位图\n")
|
||
f.write("\n")
|
||
|
||
return report_path
|
||
|
||
|
||
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) -> Path:
|
||
"""写入manifest"""
|
||
elapsed_sec = (end_time - start_time).total_seconds()
|
||
|
||
manifest_data = {
|
||
"manifest_version": "1.0",
|
||
"job_id": job_id,
|
||
"stage": "pep",
|
||
"run_id": run_id,
|
||
"prefix": prefix,
|
||
"timestamp": end_time.isoformat(),
|
||
"status": "success",
|
||
"elapsed_sec": round(elapsed_sec, 2),
|
||
"tool_versions": tool_versions,
|
||
"artifacts": artifacts,
|
||
"parameters": parameters,
|
||
"metrics": metrics
|
||
}
|
||
|
||
final_dir = workspace / "final"
|
||
final_dir.mkdir(parents=True, exist_ok=True)
|
||
for key in (
|
||
"scope_id",
|
||
"scope_mode",
|
||
"scope_file",
|
||
"target_species",
|
||
"reference_species",
|
||
"comparative_species",
|
||
"selected_species",
|
||
"family",
|
||
):
|
||
if key in parameters:
|
||
manifest_data[key] = parameters[key]
|
||
|
||
manifest_path = final_dir / f"{prefix}.pep.manifest.json"
|
||
with open(manifest_path, "w", encoding='utf-8') as f:
|
||
json.dump(manifest_data, f, indent=2, ensure_ascii=False)
|
||
|
||
if parameters.get("scope_id"):
|
||
scoped_manifest = final_dir / f"{prefix}.{sanitize_token(str(parameters['scope_id']))}.pep.manifest.json"
|
||
manifest_data["scope_manifest"] = str(scoped_manifest.resolve())
|
||
with open(scoped_manifest, "w", encoding="utf-8") as f:
|
||
json.dump(manifest_data, f, indent=2, ensure_ascii=False)
|
||
|
||
return manifest_path
|
||
|
||
|
||
def main():
|
||
parser = argparse.ArgumentParser(description="蛋白特征分析")
|
||
parser.add_argument("--prefix", required=True)
|
||
parser.add_argument("--workspace", required=True)
|
||
parser.add_argument("--input-dir")
|
||
parser.add_argument("--fasta", help="蛋白FASTA文件；用于显式指定PEP输入")
|
||
parser.add_argument("--scope-file", help="target/reference/comparative scope JSON")
|
||
parser.add_argument("--cpu", type=int, default=8)
|
||
parser.add_argument("--log-file", help="日志文件路径（用于前端兼容）")
|
||
args = parser.parse_args()
|
||
|
||
start_time = datetime.now()
|
||
run_id = f"pep_{start_time.strftime('%Y%m%d_%H%M%S')}"
|
||
|
||
workspace = Path(args.workspace)
|
||
job_id = workspace.name
|
||
|
||
print("=" * 60)
|
||
print(f"🚀 蛋白特征分析")
|
||
print("=" * 60)
|
||
|
||
try:
|
||
scope_file = Path(args.scope_file) if args.scope_file else None
|
||
scope_data = None
|
||
if scope_file:
|
||
if not scope_file.exists():
|
||
raise FileNotFoundError(f"Scope file not found: {scope_file}")
|
||
scope_data = read_scope(scope_file)
|
||
|
||
# 1. 准备目录
|
||
print("\n📁 准备目录...")
|
||
paths = prepare_directories(workspace)
|
||
|
||
# 2. 获取工具版本
|
||
print("\n🔧 检查工具...")
|
||
tool_versions = get_tool_versions()
|
||
|
||
# 3. 查找输入文件
|
||
if scope_data:
|
||
protein_file = resolve_scope_protein_fasta(workspace, args.prefix, scope_data, args.fasta)
|
||
elif args.fasta:
|
||
protein_file = Path(args.fasta)
|
||
if not protein_file.exists() or protein_file.stat().st_size == 0:
|
||
print(f"❌ 找不到蛋白文件: {protein_file}")
|
||
return 1
|
||
protein_file = ensure_protein_fasta_without_stop_markers(protein_file)
|
||
else:
|
||
protein_file = workspace / "final" / f"{args.prefix}.final_members.pep.fa"
|
||
if not protein_file.exists() or protein_file.stat().st_size == 0:
|
||
print(f"❌ 找不到蛋白文件: {protein_file}")
|
||
return 1
|
||
protein_file = ensure_protein_fasta_without_stop_markers(protein_file)
|
||
|
||
print(f" 蛋白文件: {protein_file}")
|
||
|
||
# 解析序列
|
||
sequences = parse_fasta(protein_file)
|
||
print(f" 包含 {len(sequences)} 个序列")
|
||
|
||
artifacts = {}
|
||
metrics = {"protein_count": len(sequences)}
|
||
|
||
# 4. 理化性质分析
|
||
print("\n🧮 理化性质分析...")
|
||
properties = calculate_physicochemical_properties(sequences)
|
||
|
||
# 绘制理化性质图
|
||
props_pdf = paths["final_dir"] / f"{args.prefix}.pep_properties.pdf"
|
||
props_png = paths["final_dir"] / f"{args.prefix}.pep_properties.png"
|
||
|
||
if plot_properties_r(properties, props_pdf, props_png):
|
||
print(f" ✓ 理化性质图: {props_pdf}, {props_png}")
|
||
artifacts["properties_pdf"] = {"path": f"final/{args.prefix}.pep_properties.pdf"}
|
||
artifacts["properties_png"] = {"path": f"final/{args.prefix}.pep_properties.png"}
|
||
|
||
# 5. 跨膜结构预测
|
||
print("\n🌊 跨膜结构预测 (DeepTMHMM)...")
|
||
tmhmm_result = run_deeptmhmm(protein_file, paths["analysis_dir"])
|
||
|
||
metrics["tmhmm_success"] = bool(tmhmm_result.get("success"))
|
||
if not tmhmm_result.get("success"):
|
||
metrics["tmhmm_error"] = str(tmhmm_result.get("stderr") or tmhmm_result.get("error") or "DeepTMHMM unavailable")[:500]
|
||
if tmhmm_result["success"]:
|
||
# 绘制跨膜结构图
|
||
tmhmm_pdf = paths["final_dir"] / f"{args.prefix}.pep_tmhmm.pdf"
|
||
tmhmm_png = paths["final_dir"] / f"{args.prefix}.pep_tmhmm.png"
|
||
|
||
if plot_tmhmm_results(tmhmm_result["topologies"], tmhmm_pdf, tmhmm_png):
|
||
print(f" ✓ 跨膜结构图: {tmhmm_pdf}, {tmhmm_png}")
|
||
artifacts["tmhmm_pdf"] = {"path": f"final/{args.prefix}.pep_tmhmm.pdf"}
|
||
artifacts["tmhmm_png"] = {"path": f"final/{args.prefix}.pep_tmhmm.png"}
|
||
|
||
# 6. 亚细胞定位预测
|
||
print("\n🎯 亚细胞定位预测...")
|
||
|
||
# Plant-mPLoc code is retained but disabled by default.
|
||
plant_mploc_file = paths["analysis_dir"] / f"{args.prefix}_plant_mploc.txt"
|
||
if PLANT_MPLOC_ENABLED:
|
||
plant_result = run_plant_mploc(protein_file, plant_mploc_file)
|
||
else:
|
||
print(" Plant-mPLoc已关闭；亚细胞定位默认使用WoLF PSORT")
|
||
plant_result = {
|
||
"success": False,
|
||
"disabled": True,
|
||
"method": "Plant-mPLoc website automation",
|
||
"stderr": "Plant-mPLoc disabled; WoLF PSORT is the configured subcellular localization method",
|
||
}
|
||
|
||
# WoLF PSORT
|
||
wolfpsort_file = paths["analysis_dir"] / f"{args.prefix}_wolfpsort.txt"
|
||
wolf_result = run_wolfpsort(protein_file, wolfpsort_file)
|
||
metrics["subcellular_method"] = "WoLF PSORT"
|
||
metrics["plant_mploc_enabled"] = PLANT_MPLOC_ENABLED
|
||
metrics["plant_mploc_success"] = bool(plant_result.get("success"))
|
||
metrics["wolfpsort_success"] = bool(wolf_result.get("success"))
|
||
if plant_result.get("disabled"):
|
||
metrics["plant_mploc_status"] = "disabled"
|
||
elif not plant_result.get("success"):
|
||
metrics["plant_mploc_error"] = str(plant_result.get("stderr") or plant_result.get("error") or "Plant-mPLoc unavailable")[:500]
|
||
if not wolf_result.get("success"):
|
||
metrics["wolfpsort_error"] = str(wolf_result.get("stderr") or wolf_result.get("error") or "WoLF PSORT unavailable")[:500]
|
||
|
||
# 绘制亚细胞定位图；默认只接受WoLF PSORT逐蛋白预测作为数量来源。
|
||
if wolf_result["success"]:
|
||
subcell_pdf = paths["final_dir"] / f"{args.prefix}.pep_subcellular.pdf"
|
||
subcell_png = paths["final_dir"] / f"{args.prefix}.pep_subcellular.png"
|
||
|
||
if plot_subcellular_results(plant_mploc_file, wolfpsort_file, subcell_pdf, subcell_png):
|
||
sub_counts, sub_source, sub_rows = resolve_subcellular_counts(plant_mploc_file, wolfpsort_file)
|
||
metrics["subcellular_source"] = sub_source
|
||
metrics["subcellular_prediction_count"] = sub_rows
|
||
metrics["subcellular_counts"] = sub_counts
|
||
print(f" ✓ 亚细胞定位图: {subcell_pdf}, {subcell_png}")
|
||
artifacts["subcellular_pdf"] = {"path": f"final/{args.prefix}.pep_subcellular.pdf"}
|
||
artifacts["subcellular_png"] = {"path": f"final/{args.prefix}.pep_subcellular.png"}
|
||
artifacts["subcellular_data"] = {"path": "final/subcellular_data.tsv"}
|
||
|
||
# 7. 生成报告
|
||
end_time = datetime.now()
|
||
|
||
print("\n📝 生成报告...")
|
||
parameters = {
|
||
"prefix": args.prefix,
|
||
"cpu": args.cpu,
|
||
"protein_fasta": str(protein_file),
|
||
"protein_count": count_fasta_records(protein_file),
|
||
}
|
||
parameters.update(scope_parameters(scope_file, scope_data))
|
||
|
||
report_path = write_report(workspace, args.prefix, metrics, parameters, tool_versions)
|
||
artifacts["report"] = {"path": f"final/{args.prefix}.pep_report.md"}
|
||
|
||
# 8. 生成manifest
|
||
print("\n📋 生成Manifest...")
|
||
manifest_path = write_manifest(
|
||
workspace, job_id, run_id, args.prefix,
|
||
start_time, end_time, artifacts, metrics, parameters, tool_versions
|
||
)
|
||
|
||
print("\n" + "=" * 60)
|
||
print(f"✅ 蛋白特征分析完成!")
|
||
print("=" * 60)
|
||
|
||
return 0
|
||
|
||
except Exception as e:
|
||
print(f"\n❌ 失败: {e}")
|
||
import traceback
|
||
traceback.print_exc()
|
||
return 1
|
||
|
||
|
||
if __name__ == "__main__":
|
||
sys.exit(main())
|