Skip to content

D
DynaMate
Commit 63d7e209 by cassandra
Options

[fix] bug related to multiligand systems

parent 52cf521b
Showing with 162 additions and 24 deletions
src/agents/md_agent.py
......@@ -52,7 +52,71 @@ class MDAgent(BaseAgent):
self.completed_steps = []
self.completed_summary = ""
self.EXPECTED_FILES = ["md.tpr", "md.xtc", "md.edr", "md.log", "md.gro"]
self.EXPECTED_FILES_PROTEIN_TEMPLATE = [
"{pdb_id}.pdb",
"{pdb_id}_prepared.pdb",
"{pdb_id}_prepared_capped.pdb",
"{pdb_id}_prepared_capped_his.pdb",
"{pdb_id}.prmtop",
"{pdb_id}.inpcrd",
"{pdb_id}_tleap.pdb",
"topol.top",
"{pdb_id}.gro",
"topol_without_posre.top",
"em.gro",
"nvt.gro",
"nvt.xtc",
"npt.gro",
"npt.xtc",
"temperature.xvg",
"pressure.xvg",
"density.xvg",
"potential.xvg",
"md.gro",
"md.xtc",
"rmsd.xvg",
"rmsd_xtal.xvg",
"rmsf.xvg",
"gyrate.xvg",
"hbnum_prot_wat.xvg",
"hbnum_sidechain.xvg",
]
self.EXPECTED_FILES_TEMPLATE = [
"{pdb_id}.pdb",
"{pdb_id}_prepared.pdb",
"{ligand_name}.pdb",
["{ligand_name}_h.pdb", "{ligand_name}_1_h.pdb"],
"{pdb_id}_prepared_capped.pdb",
"{pdb_id}_prepared_capped_his.pdb",
["{ligand_name}_h.prepi", "{ligand_name}_1_h.prepi"],
"{ligand_name}.frcmod",
"complex.pdb",
"complex.prmtop",
"complex.inpcrd",
"complex_tleap.pdb",
"topol.top",
"complex.gro",
"topol_without_posre.top",
"{ligand_name}.gro",
"em.gro",
"nvt.gro",
"nvt.xtc",
"npt.gro",
"npt.xtc",
"temperature.xvg",
"pressure.xvg",
"density.xvg",
"potential.xvg",
"md.gro",
"md.xtc",
"rmsd.xvg",
"rmsd_xtal.xvg",
"rmsf.xvg",
"gyrate.xvg",
"hbnum_prot_lig.xvg",
"hbnum_prot_wat.xvg",
"hbnum_sidechain.xvg",
]
self.logger.info(f"MDAgent initialized.")
......@@ -60,17 +124,21 @@ class MDAgent(BaseAgent):
if (tool_name in ("gromacs_production", "gromacs_equil", "gromacs_analysis")) and (
" failed with return code " in tool_call
):
self.logger.error(f"{tool_call}")
return False
raise ToolOutputError(f"Gromacs tool execution failed: {tool_call}")
if tool_name in ("run_tleap", "run_tleap_ligand"):
if tool_name in ("run_tleap", "run_tleap_ligand", "param_ligand"):
if "tleap run failed with error:" in tool_call:
self.logger.error(f"{tool_call}")
return False
if "Ligand parameterization failed with error:" in tool_call:
self.logger.error(f"{tool_call}")
return False
raise ToolOutputError(f"TLEaP run failed {tool_call}")
if "ParmEd failed:" in tool_call:
self.logger.error(f"{tool_call}")
return False
raise ToolOutputError(f"ParmEd failed {tool_call}")
return True
......@@ -200,12 +268,56 @@ class MDAgent(BaseAgent):
return list(self.ESSENTIAL_STEPS)
def _format_expected_files(self, templates):
values = {
"pdb_id": self.pdb_id.upper(),
"ligand_name": self.ligand_name.upper(),
}
formatted = []
for item in templates:
if isinstance(item, (list, tuple)):
formatted.append(
[s.format(**values) for s in item]
)
else:
formatted.append(
item.format(**values)
)
return formatted
def _resolve_file(self, f):
if isinstance(f, (list, tuple)):
for candidate in f:
path = self.sandbox_dir / candidate
if path.exists():
return candidate
return f[0]
return f
def _pipeline_successful(self) -> bool:
"""Check whether the full MD pipeline completed successfully."""
missing = [f for f in self.EXPECTED_FILES if not (self.sandbox_dir / f).exists()]
expected_files = (
self._format_expected_files(self.EXPECTED_FILES_TEMPLATE)
if self.ligand_name
else self._format_expected_files(self.EXPECTED_FILES_PROTEIN_TEMPLATE)
)
if missing:
self.logger.error(f"Pipeline incomplete: missing final outputs {missing}")
missing_or_empty = []
for f in expected_files:
resolved = self._resolve_file(f)
path = self.sandbox_dir / resolved
if not path.exists() or path.stat().st_size == 0:
missing_or_empty.append(resolved)
if missing_or_empty:
self.logger.error(
f"Pipeline incomplete: missing final outputs {missing_or_empty}"
)
return False
return True
......@@ -260,7 +372,8 @@ class MDAgent(BaseAgent):
if self.ligand_name and success:
user_prompt = "\n==========\nWould you like me to calculate the free energy of binding for your protein-ligand system using the MMPBSA tool? (yes/no) \n\n"
user_answer = input(user_prompt).strip().lower()
# user_answer = input(user_prompt).strip().lower()
user_answer = "yes"
if user_answer in ("yes", "y"):
self.logger.info("Running MMPBSA calculation...")
......
src/constants.py
......@@ -16,4 +16,4 @@ DATA_DIR = Path(__file__).resolve().parent.parent / "sandbox"
AGENT_LOGS = Path(__file__).resolve().parent.parent / "agent_logs"
JSON_LOG_FILE = AGENT_LOGS / "agent_runs.jsonl"
MMPBSA_ENV_DIR = Path("/path/to/your/envs/mmpbsa")
\ No newline at end of file
MMPBSA_ENV_DIR = Path("/home/hackathon/miniforge3/envs/gmxMMPBSA/bin/gmx_MMPBSA")
\ No newline at end of file
src/scripts/equil_Gromacs.sh
......@@ -90,12 +90,33 @@ EOF
if grep -E "system1 +2" topol.top; then #special case for two identical chains named system1
echo "You have two identical chains named system1, therefore only one position restraint file for the first chain will be created." >> $LOG_FILE 2>&1
echo "Creating group for residues ${ranges[0]}" >> $LOG_FILE 2>&1
echo -e "ri ${ranges[0]}\n2 & \"r_${ranges[0]}\"\nq" | $GMX make_ndx -f em.gro -n index.ndx -o index.ndx >> $LOG_FILE 2>&1
if grep -Fq "[ r_${ranges[0]} ]" index.ndx; then # -F for fixed string (so can use [] without putting "^\[ r_${ranges[0]} \]", -q for quiet)
echo "Group r_${ranges[0]} already exists in index.ndx" >> $LOG_FILE 2>&1
else
echo "Adding group r_${ranges[0]} to index.ndx" >> $LOG_FILE 2>&1
echo -e "ri ${ranges[0]}\nq" | $GMX make_ndx -f em.gro -n index.ndx -o index.ndx >> $LOG_FILE 2>&1
fi
if grep -Fq "[ Protein-H_&_r_${ranges[0]} ]" index.ndx; then # -F for fixed string (so can use [] without putting "^\[ r_${ranges[0]} \]", -q for quiet)
echo "Group Protein-H_&_r_${ranges[0]} already exists in index.ndx" >> $LOG_FILE 2>&1
else
echo "Adding group Protein-H_&_r_${ranges[0]} to index.ndx" >> $LOG_FILE 2>&1
echo -e "2 & \"r_${ranges[0]}\"\nq" | $GMX make_ndx -f em.gro -n index.ndx -o index.ndx >> $LOG_FILE 2>&1
fi
else
i=1
for range in "${ranges[@]}"; do
echo "Creating group for residues $range..." >> $LOG_FILE 2>&1
echo -e "ri $range\n2 & \"r_${range}\"\nq" | $GMX make_ndx -f em.gro -n index.ndx -o index.ndx >> $LOG_FILE 2>&1
if grep -Fq "[ r_$range ]" index.ndx; then # -F for fixed string (so can use [] without putting "^\[ r_${ranges[0]} \]", -q for quiet)
echo "Group r_$range already exists in index.ndx" >> $LOG_FILE 2>&1
else
echo "Adding group r_$range to index.ndx" >> $LOG_FILE 2>&1
echo -e "ri $range\nq" | $GMX make_ndx -f em.gro -n index.ndx -o index.ndx >> $LOG_FILE 2>&1
fi
if grep -Fq "[ Protein-H_&_r_$range ]" index.ndx; then # -F for fixed string (so can use [] without putting "^\[ r_${ranges[0]} \]", -q for quiet)
echo "Group Protein-H_&_r_$range already exists in index.ndx" >> $LOG_FILE 2>&1
else
echo "Adding group Protein-H_&_r_$range to index.ndx" >> $LOG_FILE 2>&1
echo -e "2 & \"r_$range\"\nq" | $GMX make_ndx -f em.gro -n index.ndx -o index.ndx >> $LOG_FILE 2>&1
fi
((i++))
done
fi
......@@ -103,7 +124,10 @@ fi
# Step 6: generate posre.itp for each chain
if grep -E "system1 +2" topol.top; then #special case for two identical chains named system1
group_name="Protein-H_&_r_${ranges[0]}"
if ! ls posre.itp 1> /dev/null 2>&1;then
echo "Generating position restraints for chain1" >> $LOG_FILE 2>&1
echo "$group_name" | $GMX genrestr -f em.gro -n index.ndx -o "posre.itp" -fc 1000 1000 1000 >> $LOG_FILE 2>&1
fi
else
i=1
for range in "${ranges[@]}"; do
......
src/tools/ligand_tools.py
......@@ -143,7 +143,7 @@ def param_ligand(sandbox_dir: str, ligand_files: str | list[str], ligand_name: s
# Create frcmod file using parmchk2
cmd = shlex.split(f"parmchk2 -i {sandbox_dir}/{ligand_stem}.mol2 -f mol2 -o {sandbox_dir}/{ligand_stem}.frcmod")
run_5 = subprocess.run(cmd, cwd=sandbox_dir, capture_output=True, text=True, check=True)
run_5 = subprocess.run(cmd, cwd=sandbox_dir, capture_output=True, text=True)
if run_5.returncode != 0:
error_text = "\n".join(filter(None, [run_5.stderr, run_5.stdout]))
return f"Ligand parameterization failed with error: {error_text}"
......
src/tools/pdb_tools.py
......@@ -151,14 +151,14 @@ def prepare_pdb_file_ligand(sandbox_dir: str, pdb_id: str, ligand_name: str = No
# Extract ligand
if (ligand_name is not None) and (ligand_name != "XXX") and (ligand_name != "None") and (ligand_name != "None_h"):
# Count number of ligands
resnums = set()
ligand_keys = set()
with open(f"{sandbox_dir}/{pdb_id}.pdb") as f:
for line in f:
if line.startswith("HETATM") and ligand_name in line:
if line.startswith("HETATM") and line[17:20].strip() == ligand_name:
chain = line[21]
resnum = int(line[22:26])
resnums.add(resnum)
print("resnums is: ", resnums)
num_ligands = len(resnums)
ligand_keys.add((chain, resnum))
num_ligands = len(ligand_keys)
logger.info(f"IMPORTANT: Number of ligands {ligand_name} found: {num_ligands}")
if num_ligands == 0:
......@@ -183,12 +183,13 @@ def prepare_pdb_file_ligand(sandbox_dir: str, pdb_id: str, ligand_name: str = No
# Collect HETATM lines by residue number
with open(f"{sandbox_dir}/{pdb_id}.pdb", "r") as infile:
for line in infile:
if line.startswith("HETATM") and ligand_name in line:
resnum = int(line[22:26]) # residue number column
ligands[resnum].append(line)
if line.startswith("HETATM") and line[17:20].strip() == ligand_name:
chain = line[21].strip() or "_"
resnum = int(line[22:26])
ligands[(chain, resnum)].append(line)
# Write one file per ligand
for i, (resnum, atom_lines) in enumerate(ligands.items(), start=1):
for i, ((chain, resnum), atom_lines) in enumerate(ligands.items(), start=1):
ligand_pdb_file = f"{sandbox_dir}/{ligand_name}_{i}.pdb"
ligand_pdb_files_list.append(ligand_pdb_file)
with open(ligand_pdb_file, "w") as outfile:
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment