Merge pull request #140 from talagayev/add_last_frame

Add last frame

Author: talagayev

CHANGELOG.md

@@ -22,7 +22,7 @@ The rules for this file:
 talagayev, NDoering99
 
 ### Added
-- Changed to the use of `RDKitConverter` for Ligand recognition without `-l` flag in `OpenMMDL Analysis` (2025-04-10)
+- Addition of option to select final frame in OpenMMDL Analysis (2025-04-23, Issue #136, PR #140)
 - Addition of `MDontallo` Visualization (2025-04-11)
 - Addition of `SMIRNOFF` small molecule force field (2025-04-10, Issue #76)
 - Addition of `PyMOL` support for visualization (2025-04-10)
@@ -37,6 +37,8 @@ talagayev, NDoering99
 
 ### Changed
 - The class `TrajectorySaver` was moved to `trajectory_saving.py` (2025-04-24, Issue #141, PR #142)
+- Changed to the use of `RDKitConverter` for Ligand recognition
+  without `-l` flag in `OpenMMDL Analysis` (2025-04-10)
 
 ### Deprecated
 <!-- Soon-to-be removed features -->

docs/openmmdl_analysis.rst

@@ -32,6 +32,7 @@ Optional:
     -l = Ligand in SDF format
     -b = binding mode threshold. Is used to remove interactions under the defined procentual occurence from the binding mode generation. The default is 40% (accepted values: 0-100)
     -df = Dataframe (use if the interactions were already calculated, default name would be "interactions_gathered.csv")
+    -f = final frame of the analysis (if you want to analyze only a certain part of the trajectory). The default will be the full simulation trajectory analysis.
     -m = minimal transition threshold. Is used for the display of the binding mode transitions in the Markov state chains network figure. The default value is 1
     -c = CPU count, specify how many CPUs should be used, default is half of the CPU count.
     -p = Generate .pml files for pharmacophore visualization. The default is False (accepted values: True/False)

openmmdl/openmmdl_analysis/openmmdlanalysis.py

@@ -27,6 +27,7 @@
 from rdkit.Chem.Draw import rdMolDraw2D
 from plip.basic import config
 from MDAnalysis.analysis import rms
+from MDAnalysis.coordinates.DCD import DCDWriter
 from tqdm import tqdm
 
 
@@ -102,6 +103,12 @@ def main():
         help='Dataframe (use if the interactions were already calculated, default name would be "df_all.csv")',
         default=None,
     )
+    parser.add_argument(
+        "-f",
+        dest="frames",
+        help="final frame of the analysis (if you want to analyze only a certain part of the trajectory)",
+        default=None,
+    )
     parser.add_argument(
         "-m",
         dest="min_transition",
@@ -162,14 +169,12 @@ def main():
         help="Calculate the representative frame for each binding mode. Defaults to False",
         default=False,
     )
-
     parser.add_argument(
         "--watereps",
         dest="water_eps",
         help="Set the Eps for clustering, this defines how big clusters can be spatially in Angstrom",
         default=1.0,
     )
-
     parser.add_argument(
         "--figure",
         dest="figure_type",
@@ -268,6 +273,14 @@ def main():
     if not pdb_md:
         pdb_md = mda.Universe(topology, trajectory)
 
+    if frames != None:
+        print(f"\033[1mWriting out the first {frames} out for the analysis\033[0m")
+        ag = pdb_md.select_atoms("all")
+        with DCDWriter(f"{frames}.dcd", ag.n_atoms) as w:
+            for ts in pdb_md.trajectory[: int(frames) + 1]:
+                w.write(ag)
+        pdb_md = mda.Universe(topology, f"{frames}.dcd")
+
     trajsaver = TrajectorySaver(pdb_md, ligand, special_ligand, receptor_nucleic)
 
     # TODO maybe put this part into a function possibly in visualization_functions.py TrajectorySaver
@@ -815,6 +828,8 @@ def main():
             topology, "representative_waters.pdb", water_eps
         )
 
+    if frames != None:
+        os.remove(f"{frames}.dcd")
 
 if __name__ == "__main__":
     main()