v0.10.3.dev0: convert multiple .star files into a single EMDB-SFF f…

…ile; new class `sfftk.formats.star.RelionMultiStarSegmentation`; however, only one mask/particle will be accepted for all, which is handy when displaying the particles with different colours on the same viz.

sfftk/__init__.py

@@ -3,4 +3,4 @@
 
 BASE_DIR = os.path.dirname(__file__)
 
-SFFTK_VERSION = 'v0.10.2.dev1'
+SFFTK_VERSION = 'v0.10.3.dev0'

sfftk/core/parser.py

@@ -29,7 +29,7 @@
 __updated__ = '2018-02-14'
 
 VERBOSITY_RANGE = _xrange(4)
-MULTI_FILE_FORMATS = ['stl', 'map', 'mrc', 'rec']
+MULTI_FILE_FORMATS = ['stl', 'map', 'mrc', 'rec', 'star']
 PREPABLE_FILE_FORMATS = ['mrc', 'map', 'rec']
 RESCALABLE_FILE_FORMATS = ['stl']
 # some file extensions are used by multiple file types
@@ -497,9 +497,13 @@
          "1) the mask labels (key: 'mask_to_label') and "
          "2) the hierarchical relationship between labels (key: 'label_tree')"
 )
-multi_or_label_mutex_parser.add_argument(
+# multi_or_label_mutex_parser.add_argument(
+#     '--subtomogram-average',
+#     help="the result of subtomogram averaging in CCP4 format (.mrc, .map, .rec)"
+# )
+convert_parser.add_argument(
     '--subtomogram-average',
-    help="the result of subtomogram averaging in CCP4 format (.mrc, .map, .rec)"
+    help="the result of subtomogram averaging or a particle mask for visualisation in CCP4 format (.mrc, .map, .rec)"
 )
 convert_parser.add_argument(
     '--image-name-field',

sfftk/formats/star.py

@@ -15,6 +15,7 @@
 
 """
 import numpy
+import pathlib
 import sfftkrw.schema.adapter_v0_8_0_dev1 as schema
 from sfftkrw.core.print_tools import print_date
 
@@ -30,30 +31,35 @@ class RelionStarHeader(mapformat.MaskHeader):
 class RelionStarSegment(Segment):
     """Class representing a Relion STAR file segment"""
 
-    def __init__(self, particles: starreader.StarTable, euler_angle_convention='ZYZ', degrees=True, verbose=False):
+    def __init__(self, particles: starreader.StarTable, euler_angle_convention='ZYZ', degrees=True, verbose=False,
+                 name="Particle refined using subtomogram averaging"):
         self._particles = particles
         self._euler_angle_convention = euler_angle_convention
         self._degrees = degrees
         self._verbose = verbose
+        self._name = name
 
     def convert(self, **kwargs):
         """Convert the segment to an EMDB-SFF segment"""
         segment = schema.SFFSegment()
         # metadata
         segment.biological_annotation = schema.SFFBiologicalAnnotation()
-        segment.biological_annotation.name = "Particle refined using subtomogram averaging"
+        segment.biological_annotation.name = self._name
         segment.colour = schema.SFFRGBA(random_colour=True)
         segment.shape_primitive_list = schema.SFFShapePrimitiveList()
-        transforms = schema.SFFTransformList()
+        if 'transforms' in kwargs:
+            transforms = kwargs['transforms']
+        else:
+            transforms = schema.SFFTransformList()
         if self._verbose:
             print_date(f"Using Euler angle convention: {self._euler_angle_convention}")
             print_date(f"Euler angles in degrees: {not self._degrees}")
-        for id, particle in enumerate(self._particles, start=1):
+        for particle in self._particles:
             transform = schema.SFFTransformationMatrix.from_array(
                 particle.to_affine_transform(
                     axes=self._euler_angle_convention,
                     degrees=self._degrees
-                ), id=id
+                )
             )
             shape = schema.SFFSubtomogramAverage(
                 lattice_id=kwargs.get('lattice_id'),
@@ -81,7 +87,8 @@ def __init__(self, fn, particle_fn, euler_angle_convention='ZYZ', degrees=True,
                 self._segmentation.tables['_rln'],
                 euler_angle_convention=self._euler_angle_convention,
                 degrees=self._degrees,
-                verbose=_kwargs.get('verbose', False)
+                verbose=_kwargs.get('verbose', False),
+                name=pathlib.Path(self._fn).name
             )]
 
     @property
@@ -143,8 +150,96 @@ def convert(self, name=None, software_version=None, processing_details=None, det
         segmentation.lattice_list.append(lattice)
         segmentation.segment_list = schema.SFFSegmentList()
         for segment in self.segments:
-            _segment, _transforms = segment.convert(lattice_id=lattice.id)
+            _segment, _transforms = segment.convert(lattice_id=lattice.id, transforms=segmentation.transform_list)
+            segmentation.segment_list.append(_segment)
+            segmentation.transform_list = _transforms
+        segmentation.details = details
+        return segmentation
+
+
+class RelionMultiStarSegmentation(Segmentation):
+    def __init__(self, fn_list, particle_fn, euler_angle_convention='ZYZ', degrees=True, *_args, **_kwargs):
+        """Initialise the segmentation"""
+        self._fns = fn_list
+        self._particle_fn = particle_fn
+        self._euler_angle_convention = euler_angle_convention
+        self._degrees = degrees
+        self._segmentation = list()
+        for fn in self._fns:
+            self._segmentation.append(starreader.get_data(fn, *_args, **_kwargs))
+        self._density = mapreader.get_data(self._particle_fn, *_args, **_kwargs)
+        self._segments = list()
+        for index, _segmentation in enumerate(self._segmentation):
+            self._segments.append(
+                RelionStarSegment(
+                    _segmentation.tables['_rln'],
+                    euler_angle_convention=self._euler_angle_convention,
+                    degrees=self._degrees,
+                    verbose=_kwargs.get('verbose', False),
+                    name=pathlib.Path(self._fns[index]).name
+                ))
+
+    @property
+    def header(self, ):
+        """Return the header"""
+        return RelionStarHeader(self._density)
+
+    @property
+    def segments(self):
+        """Return the segments"""
+        return self._segments
+
+    def convert(self, name=None, software_version=None, processing_details=None, details=None, verbose=False,
+                transform=None):
+        """Convert the segmentation to an EMDB-SFF segmentation"""
+        segmentation = schema.SFFSegmentation()
+        # metadata
+        segmentation.name = name if name is not None else "RELION Subtomogram Average"
+        segmentation.primary_descriptor = "shape_primitive_list"
+        segmentation.software_list = schema.SFFSoftwareList()
+        segmentation.software_list.append(
+            schema.SFFSoftware(
+                name='RELION',
+                version=software_version if software_version is not None else 'v4.0',
+                processing_details=processing_details
+            )
+        )
+        segmentation.details = details
+        # transforms
+        segmentation.transform_list = schema.SFFTransformList()
+        if transform is not None:
+            _transform = schema.SFFTransformationMatrix.from_array(transform)
+            segmentation.transform_list.append(
+                _transform
+            )
+        else:
+            _transform = schema.SFFTransformationMatrix.from_array(
+                numpy.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], ]))
+            segmentation.transform_list.append(
+                _transform
+            )
+        # lattice: we need to know the lattice because we reference it in the segment
+        segmentation.lattice_list = schema.SFFLatticeList()
+        lattice = schema.SFFLattice(
+            mode=self.header.mode,
+            endinaness=self.header.endianness,
+            size=schema.SFFVolumeStructure(
+                cols=self.header.cols,
+                rows=self.header.rows,
+                sections=self.header.sections
+            ),
+            start=schema.SFFVolumeIndex(
+                cols=self.header.start_cols,
+                rows=self.header.start_rows,
+                sections=self.header.start_sections
+            ),
+            data=self._density.voxels
+        )
+        segmentation.lattice_list.append(lattice)
+        segmentation.segment_list = schema.SFFSegmentList()
+        for segment in self.segments:
+            _segment, _transforms = segment.convert(lattice_id=lattice.id, transforms=segmentation.transform_list)
             segmentation.segment_list.append(_segment)
-            segmentation.transform_list.extend(_transforms)
+            segmentation.transform_list = _transforms
         segmentation.details = details
         return segmentation

sfftk/sff.py

@@ -125,6 +125,15 @@ def handle_convert(args, configs):  # @UnusedVariable
             elif re.match(r'.*\.stl$', args.from_file[0], re.IGNORECASE):
                 from .formats.stl import STLSegmentation
                 seg = STLSegmentation(args.from_file)
+            elif re.match(r'.*\.star$', args.from_file[0], re.IGNORECASE):
+                from .formats.star import RelionMultiStarSegmentation
+                seg = RelionMultiStarSegmentation(
+                    args.from_file, args.subtomogram_average,
+                    euler_angle_convention=args.euler_angle_convention,
+                    degrees=not args.radians,
+                    image_name_field=args.image_name_field,
+                    verbose=args.verbose
+                )
             else:
                 raise ValueError("Unknown file type '{}'".format(', '.join(args.from_file)))
         else:

sfftk/test_data/segmentations/test_data11.star

@@ -0,0 +1,25 @@
+
+data_refined_volume
+
+loop_
+_rlnMagnification #1
+_rlnDetectorPixelSize #2
+_rlnCoordinateX #3
+_rlnCoordinateY #4
+_rlnCoordinateZ #5
+_rlnAngleRot #6
+_rlnAngleTilt #7
+_rlnAnglePsi #8
+_rlnTomoName #9
+_rlnCtfImage #10
+_rlnRandomSubset #11
+_rlnPixelSize #12
+_rlnVoltage #13
+_rlnSphericalAberration #14
+_rlnMicrographName #15
+  10000.0  1.96000   487.84700  2451.94100   917.89670   -31.19659   48.16179   -60.32243  file0_1.96A.mrc  file0_ctf_1.96A.mrc  1  1.96000  300.0  1.470  file.tomostar
+  10000.0  1.96000  3404.12600  1618.88500  1425.73000  -139.90820  115.07800   -94.01349  file0_1.96A.mrc  file1_ctf_1.96A.mrc  2  1.96000  300.0  1.470  file.tomostar
+  10000.0  1.96000  2956.05500  1548.24100  1593.43900    60.32928   62.41777   -67.58646  file0_1.96A.mrc  file2_ctf_1.96A.mrc  1  1.96000  300.0  1.470  file.tomostar
+  10000.0  1.96000  1301.72100   810.92630  1071.27200   178.48200   85.17842  -106.03530  file0_1.96A.mrc  file3_ctf_1.96A.mrc  2  1.96000  300.0  1.470  file.tomostar
+  10000.0  1.96000  3215.36000   618.29010  1319.56000   164.71390  108.29200    86.71450  file0_1.96A.mrc  file4_ctf_1.96A.mrc  1  1.96000  300.0  1.470  file.tomostar
+  10000.0  1.96000  2902.59800  3907.25900   612.56200  -168.82450  131.84540   -40.54673  file0_1.96A.mrc  file5_ctf_1.96A.mrc  2  1.96000  300.0  1.470  file.tomostar

sfftk/test_data/sff/v0.8/output_emd_1181.json

@@ -45,15 +45,15 @@
       "id": 15559,
       "parent_id": 0,
       "biological_annotation": {
-        "name": "gleam center ray",
-        "description": "Eirmod rhoncusmaecenas himenaeos ipsumcurabitur, clita magnainteger proin. Aptent nibh antesuspendisse rebum, himenaeos quam nec lacinia nihil lacus turpis iaculis dui.",
-        "number_of_instances": 98,
+        "name": "bush creek soil",
+        "description": "Imperdiet netus volutpat nostra vel tortorvestibulum. Clita phasellus conguenulla nostrud facilisicurabitur sea duo, donec aliquip mus aliquet habitasse urnamorbi vivamus.",
+        "number_of_instances": 654,
         "external_references": [
           {
             "id": 0,
-            "resource": "web",
-            "url": "passbook",
-            "accession": "moves",
+            "resource": "dirt",
+            "url": "acre",
+            "accession": "adhesives",
             "label": null,
             "description": null
           }

sfftk/test_data/sff/v0.8/output_emd_1181.sff

@@ -28,12 +28,12 @@
     <segment_list>
         <segment id="15559" parent_id="0">
             <biological_annotation>
-                <name>construction tubing debt</name>
-                <description>Ametduis vestibulumnulla taciti sea, euismod quisque nec consetetur. Doming magnainteger dui iaculis justo, no muspellentesque tempus.</description>
+                <name>barometer variation fake</name>
+                <description>Sea fermentum bibendum faucibusvestibulum no, rhoncusmaecenas feugait in et blandit. Quisque imperdietaliquam velit sapien.</description>
                 <external_references>
-                    <ref id="0" resource="provisions" url="analyzers" accession="habits"/>
+                    <ref id="0" resource="hairpins" url="operand" accession="baud"/>
                 </external_references>
-                <number_of_instances>384</number_of_instances>
+                <number_of_instances>800</number_of_instances>
             </biological_annotation>
             <colour>
                 <red>0.921817600727081</red>

sfftk/unittests/test_formats.py

@@ -109,6 +109,17 @@ def read_star(self):
                 image_name_field='_rlnTomoName'
             )
 
+    def read_star_multi(self):
+        """Read RELION .star multi files"""
+        if not hasattr(self, 'star_multi_file'):
+            self.star_multi0_file = os.path.join(self.segmentations_path, 'test_data8.star')
+            self.star_multi1_file = os.path.join(self.segmentations_path, 'test_data11.star')
+            self.particle_file = os.path.join(self.segmentations_path, 'test_data.map')
+            self.image_file = os.path.join(self.segmentations_path, 'test_data.map')
+            self.star_multi_segmentation = star.RelionMultiStarSegmentation(
+                [self.star_multi0_file, self.star_multi1_file], self.particle_file,
+            )
+
     def read_stl(self):
         """Read .stl files"""
         if not hasattr(self, 'stl_file'):
@@ -210,6 +221,15 @@ def test_star_read(self):
         self.assertIsInstance(self.star_segmentation.segments, list)
         self.assertIsInstance(self.star_segmentation.segments[0], star.RelionStarSegment)
 
+    def test_multiple_star_read(self):
+        """Read multiple RELION (.star) segmentation"""
+        self.read_star_multi()
+        # assertions
+        self.assertIsInstance(self.star_multi_segmentation.header, star.RelionStarHeader)
+        self.assertIsInstance(self.star_multi_segmentation.segments, list)
+        self.assertEqual(2, len(self.star_multi_segmentation.segments))
+        self.assertIsInstance(self.star_multi_segmentation.segments[0], star.RelionStarSegment)
+
     def test_stl_read(self):
         """Read a Stereo Lithography (.stl) segmentation"""
         self.read_stl()
@@ -513,7 +533,7 @@ def test_star_convert(self):
         self.assertEqual("Something interesting", seg.details)
         segment = seg.segment_list[0]  # there is only one segment anyway
         self.assertIsNotNone(segment.biological_annotation)
-        self.assertIsNotNone(segment.biological_annotation.name)
+        self.assertEqual("test_data8.star", segment.biological_annotation.name)
         self.assertGreaterEqual(segment.biological_annotation.number_of_instances, 1)
         self.assertIsNotNone(segment.colour)
         self.assertEqual(0, len(segment.mesh_list))
@@ -533,6 +553,58 @@ def test_star_convert(self):
         self.assertEqual(1.0, sta.value)
         self.assertEqual(1, sta.transform_id)
 
+    def test_star_multi_convert(self):
+        """Convert a segmentation from multiple RELION .star files to an SFFSegmentation object"""
+        self.read_star_multi()
+        args, configs = parse_args(
+            f"convert {self.star_multi0_file} {self.star_multi1_file} --multi-file --image {self.image_file} "
+            f"--image-name-field _rlnTomoName --details 'Something interesting' "
+            f"--subtomogram-average {self.particle_file}",
+            use_shlex=True
+        )
+        transform = mapreader.compute_transform(args.image)
+        seg = self.star_multi_segmentation.convert(details=args.details, transform=transform)
+        # assertions
+        self.assertIsInstance(seg, schema.SFFSegmentation)
+        self.assertEqual("RELION Subtomogram Average", seg.name)
+        self.assertEqual(seg.version, self.schema_version)
+        self.assertEqual(seg.software_list[0].name, 'RELION')
+        self.assertEqual("shape_primitive_list", seg.primary_descriptor)
+        self.assertEqual(seg.transform_list[0].id, 0)
+        self.assertEqual(
+            '18.20999987022425 0.0 0.0 0.0 0.0 18.209998684928305 0.0 0.0 0.0 0.0 '
+            '18.209998762103872 0.0',
+            seg.transform_list[0].data
+        )
+        self.assertEqual(13, len(seg.transform_list))
+        self.assertEqual(12, seg.transform_list[-1].id)
+        self.assertEqual("Something interesting", seg.details)
+        self.assertEqual(2, len(seg.segment_list))
+        segment1, segment2 = seg.segment_list
+        self.assertIsNotNone(segment1.biological_annotation)
+        self.assertEqual("test_data8.star", segment1.biological_annotation.name)
+        self.assertEqual("test_data11.star", segment2.biological_annotation.name)
+        self.assertGreaterEqual(segment1.biological_annotation.number_of_instances, 1)
+        self.assertIsNotNone(segment1.colour)
+        self.assertEqual(0, len(segment1.mesh_list))
+        self.assertEqual(6, segment1.shape_primitive_list[-1].transform_id)
+        self.assertEqual(12, segment2.shape_primitive_list[-1].transform_id)
+        # lattice
+        lattice = seg.lattice_list[0]
+        self.assertIsInstance(lattice, schema.SFFLattice)
+        self.assertEqual(0, lattice.id)
+        self.assertEqual('float32', lattice.mode)
+        self.assertEqual('little', lattice.endianness)
+        self.assertIsInstance(lattice.size, schema.SFFVolumeStructure)
+        self.assertIsInstance(lattice.start, schema.SFFVolumeIndex)
+        self.assertIsNotNone(lattice.data)
+        sta = segment1.shape_primitive_list[0]
+        self.assertIsInstance(sta, schema.SFFSubtomogramAverage)
+        self.assertEqual(0, sta.id)
+        self.assertEqual(lattice.id, sta.lattice_id)
+        self.assertEqual(1.0, sta.value)
+        self.assertEqual(1, sta.transform_id)
+
     def test_stl_convert(self):
         """Convert a segmentation from an Stereo Lithography file to an SFFSegmentation object"""
         self.read_stl()