Merge pull request #110 from FZJ-INM1-BDA/bugfix_ebrainsDatasetSpecies

bugfix: filter dataset by species

Author: dickscheid

pytest.ini

@@ -1,7 +1,4 @@
 [pytest]
-addopts = --maxfail=1 -rf
+addopts = -rf
 testpaths = 
-    test/core
-    test/retrieval
-    test/volumes
-    test/features
+    test/

siibra/core/atlas.py

@@ -32,20 +32,50 @@ class Atlas(
     spaces, as well as common functionalities of those.
     """
 
-    def __init__(self, identifier, name):
+    @staticmethod
+    def get_species_data(species_str: str):
+        if species_str == 'human':
+            return {
+                '@id': 'https://nexus.humanbrainproject.org/v0/data/minds/core/species/v1.0.0/0ea4e6ba-2681-4f7d-9fa9-49b915caaac9',
+                'name': 'Homo sapiens'
+            }
+        if species_str == 'rat':
+            return {
+                '@id': 'https://nexus.humanbrainproject.org/v0/data/minds/core/species/v1.0.0/f3490d7f-8f7f-4b40-b238-963dcac84412',
+                'name': 'Rattus norvegicus'
+            }
+        if species_str == 'mouse':
+            return {
+                '@id': 'https://nexus.humanbrainproject.org/v0/data/minds/core/species/v1.0.0/cfc1656c-67d1-4d2c-a17e-efd7ce0df88c',
+                'name': 'Mus musculus'
+            }
+        # TODO this may not be correct. Wait for feedback and get more accurate
+        if species_str == 'monkey':
+            return {
+                '@id': 'https://nexus.humanbrainproject.org/v0/data/minds/core/species/v1.0.0/3f75b0ad-dbcd-464e-b614-499a1b9ae86b',
+                'name': 'Primates'
+            }
+
+        raise ValueError(f'species with spec {species_str} cannot be decoded')
+
+    def __init__(self, identifier, name, species = None):
         """Construct an empty atlas object with a name and identifier."""
 
         AtlasConcept.__init__(self, identifier, name, dataset_specs=[])
 
         self._parcellations = []  # add with _add_parcellation
         self._spaces = []  # add with _add_space
+        if species is not None:
+            self.species = self.get_species_data(species)
 
     def _register_space(self, space):
         """Registers another reference space to the atlas."""
+        space.atlases.add(self)
         self._spaces.append(space)
 
     def _register_parcellation(self, parcellation):
         """Registers another parcellation to the atlas."""
+        parcellation.atlases.add(self)
         self._parcellations.append(parcellation)
 
     @property
@@ -75,7 +105,7 @@ def _from_json(cls, obj):
                 f"{cls.__name__} construction attempt from invalid json format (@type={obj.get('@type')}"
             )
         if all(["@id" in obj, "spaces" in obj, "parcellations" in obj]):
-            atlas = cls(obj["@id"], obj["name"])
+            atlas = cls(obj["@id"], obj["name"], species=obj["species"])
             for space_id in obj["spaces"]:
                 if not Space.REGISTRY.provides(space_id):
                     raise ValueError(

siibra/core/parcellation.py

@@ -156,6 +156,7 @@ def __init__(
             if self._datasets_cached is None:
                 self._datasets_cached = []
             self._datasets_cached.extend(maps)
+        self.atlases = set()
 
     @property
     def regiontree(self):

siibra/core/space.py

@@ -51,6 +51,7 @@ def __init__(
         AtlasConcept.__init__(self, identifier, name, dataset_specs)
         self.src_volume_type = src_volume_type
         self.type = template_type
+        self.atlases = set()
 
     def get_template(self):
         """

siibra/features/connectivity.py

@@ -128,9 +128,9 @@ class ConnectivityProfile(RegionalFeature):
 
     show_as_log = True
 
-    def __init__(self, regionspec: str, connectivitymatrix: ConnectivityMatrix, index):
+    def __init__(self, regionspec: str, connectivitymatrix: ConnectivityMatrix, index, **kwargs):
         assert regionspec is not None
-        RegionalFeature.__init__(self, regionspec)
+        RegionalFeature.__init__(self, regionspec, **kwargs)
         self._matrix_index = index
         self._matrix = connectivitymatrix
 
@@ -184,13 +184,14 @@ def __init__(self, **kwargs):
         for _, loader in self._QUERY.get_loaders("connectivity", ".json"):
             cm = ConnectivityMatrix._from_json(loader.data)
             for parcellation in cm.parcellations:
+                species = [atlas.species for atlas in parcellation.atlases]
                 for regionname in cm.regionnames:
                     region = parcellation.decode_region(regionname, build_group=False)
                     if region is None:
                         raise RuntimeError(
                             f"Could not decode region name {regionname} in {parcellation}"
                         )
-                    self.register(ConnectivityProfile(region, cm, regionname))
+                    self.register(ConnectivityProfile(region, cm, regionname, species=species))
 
 
 class ConnectivityMatrixQuery(FeatureQuery):

siibra/features/ebrains.py

@@ -22,8 +22,8 @@
 
 
 class EbrainsRegionalDataset(RegionalFeature, EbrainsDataset):
-    def __init__(self, regionspec, kg_id, name, embargo_status):
-        RegionalFeature.__init__(self, regionspec)
+    def __init__(self, regionspec, kg_id, name, embargo_status, species = []):
+        RegionalFeature.__init__(self, regionspec, species)
         EbrainsDataset.__init__(self, kg_id, name, embargo_status)
 
     @property
@@ -46,17 +46,26 @@ class EbrainsRegionalFeatureQuery(FeatureQuery):
     _FEATURETYPE = EbrainsRegionalDataset
 
     def __init__(self, **kwargs):
-
         FeatureQuery.__init__(self)
 
         loader = EbrainsRequest(
-            query_id="siibra-kg-feature-summary-0.0.1",
+            query_id="siibra-kg-feature-summary-0_0_4",
             schema="parcellationregion",
             params={"vocab": "https://schema.hbp.eu/myQuery/"},
         )
 
         for r in loader.data.get("results", []):
+
+            species_alt = []
+            # List, keys @id, name
             for dataset in r.get("datasets", []):
+                species_alt = [
+                    *species_alt,
+                    *dataset.get('ds_specimengroup_subject_species', []),
+                    *dataset.get('s_subject_species', []),
+                ]
+            for dataset in r.get("datasets", []):
+
                 ds_id = dataset.get("@id")
                 ds_name = dataset.get("name")
                 ds_embargo_status = dataset.get("embargo_status")
@@ -65,10 +74,28 @@ def __init__(self, **kwargs):
                         f"'{ds_name}' is not an interpretable dataset and will be skipped.\n(id:{ds_id})"
                     )
                     continue
-                regionname = r.get("name", None)
+                regionname: str = r.get("name", None)
+                alias: str = r.get("alias", None)
+
+                # species defined for the current dataset
+                dataset_species = [
+                    *dataset.get('ds_specimengroup_subject_species', []),
+                    *dataset.get('s_subject_species', []),
+                ]
+
+                # if the current dataset has species defined, use the current species, else use the general speices
+                species = [*r.get("species", []), *(dataset_species if dataset_species else species_alt)] # list with keys @id, identifier, name
+
+                # filter species by @id attribute
+                unique_species = []
+                for sp in species:
+                    if sp.get('@id') in [s.get('@id') for s in unique_species]:
+                        continue
+                    unique_species.append(sp)
+
                 self.register(
                     EbrainsRegionalDataset(
-                        regionname, ds_id, ds_name, ds_embargo_status
+                        alias or regionname, ds_id, ds_name, ds_embargo_status, unique_species
                     )
                 )
 

siibra/features/feature.py

@@ -181,7 +181,7 @@ class RegionalFeature(Feature):
     TODO store region as an object that has a link to the parcellation
     """
 
-    def __init__(self, regionspec: Tuple[str, Region]):
+    def __init__(self, regionspec: Tuple[str, Region], species = [], **kwargs):
         """
         Parameters
         ----------
@@ -194,6 +194,11 @@ def __init__(self, regionspec: Tuple[str, Region]):
             )
         Feature.__init__(self)
         self.regionspec = regionspec
+        self.species = species
+
+    @property
+    def species_ids(self):
+        return [s.get('@id') for s in self.species]
 
     def match(self, concept):
         """
@@ -208,6 +213,24 @@ def match(self, concept):
         -------
         True, if match was successful, otherwise False
         """
+
+        # first check if any of
+        try:
+            if isinstance(concept, Region):
+                atlases = concept.parcellation.atlases
+            if isinstance(concept, Parcellation):
+                atlases = concept.atlases
+            if isinstance(concept, Atlas):
+                atlases = {concept}
+            if atlases:
+                # if self.species_ids is defined, and the concept is explicitly not in
+                # return False
+                if all(atlas.species.get('@id') not in self.species_ids for atlas in atlases):
+                    return False
+        # for backwards compatibility. If any attr is not found, pass
+        except AttributeError:
+            pass
+
         self._match = None
 
         # regionspec might be a specific region, then we can
@@ -227,13 +250,15 @@ def match(self, concept):
             for w in concept.key.split('_'):
                 spec = spec.replace(w.lower(), '')
             for match in concept.regiontree.find(spec):
+                # TODO what's with the mutation here?
                 self._match = match
                 return True
 
         elif isinstance(concept, Region):
             for w in concept.parcellation.key.split('_'):
                 spec = spec.replace(w.lower(), '')
             for match in concept.find(spec):
+                # TODO what's with the mutation here?
                 self._match = match
                 return True
 
@@ -245,6 +270,7 @@ def match(self, concept):
                 spec = spec.replace(w.lower(), '')
             for p in concept.parcellations:
                 for match in p.regiontree.find(spec):
+                    # TODO what's with the mutation here?
                     self._match = match
                     return True
         else:

siibra/features/receptors.py

@@ -357,9 +357,9 @@ class ReceptorDistribution(RegionalFeature, EbrainsDataset):
     TODO lazy loading could be more elegant.
     """
 
-    def __init__(self, region, kg_result):
+    def __init__(self, region, kg_result, **kwargs):
 
-        RegionalFeature.__init__(self, region)
+        RegionalFeature.__init__(self, region, **kwargs)
         EbrainsDataset.__init__(self, kg_result["identifier"], kg_result["name"])
 
         self.info = kg_result["description"]
@@ -488,13 +488,18 @@ class ReceptorQuery(FeatureQuery):
 
     def __init__(self,**kwargs):
         FeatureQuery.__init__(self)
-        kg_query = EbrainsRequest(query_id="siibra_receptor_densities").get()
-        # kg_query = ebrains.execute_query_by_id('minds', 'core', 'dataset', 'v1.0.0', )
+        kg_req = EbrainsRequest(
+            query_id="siibra_receptor_densities-0_0_2",
+            params={'vocab': 'https://schema.hbp.eu/myQuery/' }
+        )
+        kg_query = kg_req.get()
+        
         not_used = 0
         for kg_result in kg_query["results"]:
-            region_names = [e["name"] for e in kg_result["region"]]
+            region_names = [p_region["name"] for p_region in kg_result["parcellationRegion"]]
+            species = kg_result.get('species', [])
             for region_name in region_names:
-                f = ReceptorDistribution(region_name, kg_result)
+                f = ReceptorDistribution(region_name, kg_result, species=species)
                 if f.fingerprint is None:
                     not_used += 1
                 else:

test/features/test_ebrainsquery.py

@@ -1,6 +1,9 @@
+from typing import List
 import unittest
 import siibra
-
+import pytest
+from siibra.core import Parcellation, Atlas, Region
+from siibra.features.feature import Feature
 
 class TestEbrainsQuery(unittest.TestCase):
     @classmethod
@@ -17,5 +20,60 @@ def test_no_duplicates_returned(self):
         assert len(self.feat) == len(list(set(ids)))
 
 
+parameter = [
+    ('rat', 'v3', 'neocortex', {
+        'exclude': [
+            # buggy. hippocampus + difumo 512 is a big issue
+            'DiFuMo atlas (512 dimensions)'
+        ],
+        'include': [
+            # some of these still clearly doesn't look right
+            # for e.g. some v1/v2 are in here
+            # but one step at a time...
+            'Large scale multi-channel EEG in rats',
+            'Wistar rat hippocampus CA1 pyramidal cell morphologies – Extension with additional reconstructions',
+            'Auditory stimulation during the sleep-wake cycle in the freely moving rat',
+            '3D high resolution SRXTM image data of cortical vasculature of rat brain.',
+            'Density measurements of different receptors for CA1 (Hippocampus) [rat, v2.0]',
+            'Visualization of projections from insular cortex in rat with the use of anterograde tracers',
+            'Density measurements of different receptors for CA, stratum moleculare (Hippocampus) [rat, v2.0]',
+            'Density measurements of different receptors for CA2 (Hippocampus) [rat, v2.0]',
+            'PCI-like measure in rodents',
+            'Density measurements of different receptors for CA3 (Hippocampus) [rat, v2.0]',
+            'Dose-dependent effects of ketamine on spontaneous and evoked EEG activity in rats',
+            'Detailed dynamical laminar organisation in different cortical areas (in rats in vivo)',
+            'Density measurements of different receptors for CA, stratum cellulare (Hippocampus) [rat, v2.0]',
+            '3D reconstructions of pyramidal cells in rat hippocampal CA1 region',
+            'Electrophysiological data of cortical layer 6 neurons and synaptically coupled neuronal pairs',
+            'Density measurements of different receptors for DG (Hippocampus) [rat, v1.0]',
+            'Test of consciousness metrics in rodents',
+            'Morphological data of cortical layer 6 neurons and synaptically coupled neuronal pairs',
+            'Visualization of projections from posterior parietal cortex in rat with the use of anterograde tracers',
+            'Immunofluorescence data of cortical layer 6 neurons',
+            'Density measurements of different receptors for DG (Hippocampus) [rat, v2.0]',
+            'Graphical representation of rat cortical vasculature reconstructed from high resolution 3D SRXTM data.',
+            'Density measurements of different receptors for CA, stratum cellulare (Hippocampus) [rat, v1.0]',
+            'Wistar rat hippocampus CA1 pyramidal cell morphologies',
+            'Density measurements of different receptors for CA3 (Hippocampus) [rat, v1.0]',
+            'Density measurements of different receptors for CA1 (Hippocampus) [rat, v1.0]',
+            'Density measurements of different receptors for CA, stratum moleculare (Hippocampus) [rat, v1.0]',
+            'Density measurements of different receptors for CA2 (Hippocampus) [rat, v1.0]',
+            'Multi-area recordings from visual and somatosensory cortices, perirhinal cortex and hippocampal CA1']
+    })
+]
+
+@pytest.mark.parametrize('atlas_id,parc_id,region_id,inc_exc', parameter)
+def test_species(atlas_id,parc_id,region_id,inc_exc):
+    atlas:Atlas = siibra.atlases[atlas_id]
+    parc:Parcellation = atlas.parcellations[parc_id]
+    r:Region = parc.decode_region(region_id)
+    features: List[Feature] = siibra.get_features(r, 'ebrains')
+    feature_names = [f.name for f in features]
+
+    excludes: List[str] = inc_exc.get('exclude')
+    includes: List[str] = inc_exc.get('include')
+    assert all(exc not in feature_names for exc in excludes)
+    assert all(inc in feature_names for inc in includes)
+
 if __name__ == "__main__":
     unittest.main()