Relion5 input

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "pytom-match-pick"
-version = "0.7.3"
+version = "0.8.0"
 description = "PyTOM's GPU template matching module as an independent package"
 readme = "README.md"
 license = {file = "LICENSE"}

src/pytom_tm/entry_points.py

@@ -19,6 +19,7 @@
     ParseDefocus,
     BetweenZeroAndOne,
     ParseGPUIndices,
+    parse_relion5_star_data,
 )
 from pytom_tm.tmjob import load_json_to_tmjob
 from os import urandom
@@ -731,7 +732,7 @@ def match_template(argv=None):
         "--tilt-angles",
         nargs="+",
         type=str,
-        required=True,
+        required=False,
         action=ParseTiltAngles,
         help="Tilt angles of the tilt-series, either the minimum and maximum values of "
         "the tilts (e.g. --tilt-angles -59.1 60.1) or a .rawtlt/.tlt file with all the "
@@ -897,6 +898,17 @@ def match_template(argv=None):
         help="Specify a seed for the random number generator used for phase "
         "randomization for consistent results!",
     )
+    additional_group.add_argument(
+        "--relion5-tomograms-star",
+        type=pathlib.Path,
+        action=CheckFileExists,
+        required=False,
+        help="Here, you can provide a path to a RELION5 tomograms.star file (for "
+        "example "
+        "from a tomogram reconstruction job). pytom-match-pick will fetch all "
+        "the tilt-series metadata from this file and overwrite all other "
+        "metadata options.",
+    )
     device_group = parser.add_argument_group("Device control")
     device_group.add_argument(
         "-g",
@@ -922,6 +934,11 @@ def match_template(argv=None):
     args = parser.parse_args(argv)
     logging.basicConfig(level=args.log, force=True)
 
+    # parse CTF phase correction
+    phase_flip_correction = False
+    if args.tomogram_ctf_model is not None and args.tomogram_ctf_model == "phase-flip":
+        phase_flip_correction = True
+
     # combine ctf values to ctf_params list of dicts
     ctf_params = None
     if args.defocus is not None:
@@ -935,12 +952,6 @@ def match_template(argv=None):
                 "the required parameters (amplitude-contrast, "
                 "spherical-abberation or voltage) is/are missing."
             )
-        phase_flip_correction = False
-        if (
-            args.tomogram_ctf_model is not None
-            and args.tomogram_ctf_model == "phase-flip"
-        ):
-            phase_flip_correction = True
         ctf_params = [
             {
                 "defocus": defocus * 1e-6,
@@ -953,6 +964,28 @@ def match_template(argv=None):
             for defocus in args.defocus
         ]
 
+    if args.relion5_tomograms_star is not None:
+        voxel_size, tilt_angles, dose_accumulation, ctf_params, defocus_handedness = (
+            parse_relion5_star_data(
+                args.relion5_tomograms_star,
+                args.tomogram,
+                phase_flip_correction=phase_flip_correction,
+                phase_shift=args.phase_shift,
+            )
+        )
+        per_tilt_weighting = True
+    else:
+        if args.tilt_angles is None:
+            raise ValueError(
+                "Without tilt angles the missing wedge cannot be calculated. A "
+                "minimal run requires tilt angles."
+            )
+        voxel_size = args.voxel_size_angstrom
+        defocus_handedness = args.defocus_handedness
+        tilt_angles = args.tilt_angles
+        dose_accumulation = args.dose_accumulation
+        per_tilt_weighting = args.per_tilt_weighting
+
     if args.angular_search is None and args.particle_diameter is None:
         raise ValueError(
             "Either the angular search should be specifically set or a particle "
@@ -970,23 +1003,23 @@ def match_template(argv=None):
         mask_is_spherical=True
         if args.non_spherical_mask is None
         else (not args.non_spherical_mask),
-        tilt_angles=args.tilt_angles,
-        tilt_weighting=args.per_tilt_weighting,
+        tilt_angles=tilt_angles,
+        tilt_weighting=per_tilt_weighting,
         search_x=args.search_x,
         search_y=args.search_y,
         search_z=args.search_z,
         tomogram_mask=args.tomogram_mask,
-        voxel_size=args.voxel_size_angstrom,
+        voxel_size=voxel_size,
         low_pass=args.low_pass,
         high_pass=args.high_pass,
-        dose_accumulation=args.dose_accumulation,
+        dose_accumulation=dose_accumulation,
         ctf_data=ctf_params,
         whiten_spectrum=args.spectral_whitening,
         rotational_symmetry=args.z_axis_rotational_symmetry,
         particle_diameter=args.particle_diameter,
         random_phase_correction=args.random_phase_correction,
         rng_seed=args.rng_seed,
-        defocus_handedness=args.defocus_handedness,
+        defocus_handedness=defocus_handedness,
         output_dtype=np.float16 if args.half_precision else np.float32,
     )
 

src/pytom_tm/extract.py

@@ -276,6 +276,11 @@ def extract_particles(
     pixel_size = job.voxel_size
     tomogram_id = job.tomo_id
 
+    # remove relion5 reconstructed tomogram name as it messes with linking the tilt
+    # series id when extracting subtomos
+    if relion5_compat and tomogram_id.startswith("rec_"):
+        tomogram_id = tomogram_id[4:]
+
     data = []
     scores = []
 

src/pytom_tm/io.py

@@ -4,6 +4,7 @@
 import logging
 import numpy.typing as npt
 import numpy as np
+import starfile
 from contextlib import contextmanager
 from operator import attrgetter
 
@@ -481,3 +482,103 @@ def read_defocus_file(file_name: pathlib.Path) -> list[float, ...]:
         return read_txt_file(file_name)
     else:
         raise ValueError("Defocus file needs to have format .defocus or .txt")
+
+
+def parse_relion5_star_data(
+    tomograms_star_path: pathlib.Path,
+    tomogram_path: pathlib.Path,
+    phase_flip_correction: bool = False,
+    phase_shift: float = 0.0,
+) -> tuple[float, list[float, ...], list[float, ...], list[dict, ...], int]:
+    """Read RELION5 metadata from a project directory.
+
+    Parameters
+    ----------
+    tomograms_star_path: pathlib.Path
+        the tomograms.star from a RELION5 reconstruct job contains invariable metadata
+        and points to a tilt series star file with fitted values
+    tomogram_path: pathlib.Path
+        path to the tomogram for template matching; we use the name to pattern match in
+        the RELION5 star file
+    phase_flip_correction: bool, default False
+    phase_shift: float, default 0.0
+
+    Returns
+    -------
+    tomogram_voxel_size, tilt_angles, dose_accumulation, ctf_params, defocus_handedness:
+        tuple[float, list[float, ...], list[float, ...], list[dict, ...], int]
+    """
+    tomogram_id = tomogram_path.stem
+    tomograms_star_data = starfile.read(tomograms_star_path)
+
+    # match the tomo_id and check if viable
+    matches = [
+        (i, x)
+        for i, x in enumerate(tomograms_star_data["rlnTomoName"])
+        if x in tomogram_id
+    ]
+    if len(matches) == 1:
+        tomogram_meta_data = tomograms_star_data.loc[matches[0][0]]
+    else:
+        raise ValueError(
+            "Multiple or zero matches of tomogram id in RELION5 STAR file. Aborting..."
+        )
+
+    # grab the path to tilt series star for tilt angles, defocus and dose
+    tilt_series_star_path = pathlib.Path(
+        tomogram_meta_data["rlnTomoTiltSeriesStarFile"]
+    )
+    # update the path to a location we can actually find from CD
+    tilt_series_star_path = tomograms_star_path.parent.joinpath("tilt_series").joinpath(
+        tilt_series_star_path.name
+    )
+    tilt_series_star_data = starfile.read(tilt_series_star_path)
+
+    # we extract tilt angles, dose accumulation and ctf params
+    # TODO We need to have an internal structure for tilt series meta data
+    tilt_angles = list(tilt_series_star_data["rlnTomoNominalStageTiltAngle"])
+    dose_accumulation = list(tilt_series_star_data["rlnMicrographPreExposure"])
+
+    # _rlnTomoName  # 1
+    # _rlnVoltage  # 2
+    # _rlnSphericalAberration  # 3
+    # _rlnAmplitudeContrast  # 4
+    # _rlnMicrographOriginalPixelSize  # 5
+    # _rlnTomoHand  # 6
+    # _rlnTomoTiltSeriesPixelSize  # 7
+    # _rlnTomoTiltSeriesStarFile  # 8
+    # _rlnTomoTomogramBinning  # 9
+    # _rlnTomoSizeX  # 10
+    # _rlnTomoSizeY  # 11
+    # _rlnTomoSizeZ  # 12
+    # _rlnTomoReconstructedTomogram  # 13
+    tomogram_voxel_size = float(
+        tomogram_meta_data["rlnTomoTiltSeriesPixelSize"]
+        * tomogram_meta_data["rlnTomoTomogramBinning"]
+    )
+    defocus_handedness = int(tomogram_meta_data["rlnTomoHand"])
+
+    ctf_params = [
+        {
+            "defocus": defocus * 1e-6,
+            "amplitude_contrast": tomogram_meta_data["rlnAmplitudeContrast"],
+            "voltage": tomogram_meta_data["rlnVoltage"] * 1e3,
+            "spherical_aberration": tomogram_meta_data["rlnSphericalAberration"] * 1e-3,
+            "flip_phase": phase_flip_correction,
+            "phase_shift_deg": phase_shift,  # Unclear to me where the phase
+            # shifts are annotated in RELION5 tilt series file
+        }
+        for defocus in (
+            tilt_series_star_data.rlnDefocusV + tilt_series_star_data.rlnDefocusU
+        )
+        / 2
+        * 1e-4
+    ]
+
+    return (
+        tomogram_voxel_size,
+        tilt_angles,
+        dose_accumulation,
+        ctf_params,
+        defocus_handedness,
+    )