Source code for tvb.adapters.creators.siibra_base

# -*- coding: utf-8 -*-
#
#
# TheVirtualBrain-Framework Package. This package holds all Data Management, and
# Web-UI helpful to run brain-simulations. To use it, you also need to download
# TheVirtualBrain-Scientific Package (for simulators). See content of the
# documentation-folder for more details. See also http://www.thevirtualbrain.org
#
# (c) 2012-2023, Baycrest Centre for Geriatric Care ("Baycrest") and others
#
# This program is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software Foundation,
# either version 3 of the License, or (at your option) any later version.
# This program is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
# PARTICULAR PURPOSE.  See the GNU General Public License for more details.
# You should have received a copy of the GNU General Public License along with this
# program.  If not, see <http://www.gnu.org/licenses/>.
#
#
#   CITATION:
# When using The Virtual Brain for scientific publications, please cite it as explained here:
# https://www.thevirtualbrain.org/tvb/zwei/neuroscience-publications
#
#

"""
Utility functions for using siibra to extract Structural and Functional connectivities

.. moduleauthor:: Romina Baila <romina.baila@codemart.ro>
"""
import numpy as np
import siibra
from enum import Enum
from tvb.basic.logger.builder import get_logger
from tvb.datatypes import connectivity
from tvb.datatypes.graph import ConnectivityMeasure

LOGGER = get_logger(__name__)

# Concepts names
# atlases
HUMAN_ATLAS = 'Multilevel Human Atlas'  # DEFAULT, only this atlas has Struct. Conn.

# parcellations
JULICH_3_0 = 'Julich-Brain Cytoarchitectonic Atlas (v3.0.3)'  # DEFAULT
JULICH_2_9 = 'Julich-Brain Cytoarchitectonic Atlas (v2.9)'
parcellations = [JULICH_3_0, JULICH_2_9]

# cohorts
HCP_COHORT = 'HCP'  # DEFAULT
THOUSAND_BRAINS_COHORT = '1000BRAINS'


[docs] class Component2Modality(Enum): WEIGHTS = siibra.features.connectivity.StreamlineCounts TRACTS = siibra.features.connectivity.StreamlineLengths FUNCTIONAL_CONNECTIVITY = siibra.features.connectivity.FunctionalConnectivity
# ########################################## SIIBRA CREATOR INITIALIZATION #############################################
[docs] def get_cohorts_for_sc(parcellation_name): """ Given a parcellation name, return the name of all the cohorts related to it and containing data about Structural Connectivities. """ parcellation = siibra.parcellations[parcellation_name] features = siibra.features.get(parcellation, siibra.features.connectivity.StreamlineCounts) return [f.cohort.upper() for f in features]
# ######################################## SIIBRA PARAMETERS CONFIGURATION #############################################
[docs] def check_atlas_parcellation_compatible(atlas, parcellation): """ Given an atlas and a parcellation, verify that the atlas contains the parcellation, i.e. they are compatible """ return parcellation in list(atlas.parcellations)
[docs] def get_atlases_for_parcellation(parcelation): """ Given a parcelation, return all the atlases that contain it """ atlases = siibra.atlases atlases = [a for a in atlases if parcelation in list(a.parcellations)] return atlases
[docs] def get_parcellations_for_atlas(atlas): """ Given an atlas, return a list of all the parcellations inside it, which contain Structural conns. """ parcellations_available = [p for p in list(atlas.parcellations) if p.name in parcellations] return parcellations_available
[docs] def parse_subject_ids(subject_ids, cohort): """ Given a string representing subject ids or a range of subject ids; return the list containing all the included ids """ parsed_ids_as_str = [] individual_splits = subject_ids.split(';') zfill_value = 3 if cohort == HCP_COHORT else 4 # used in case of ranges for s in individual_splits: # if a range was specified if '-' in s: start, end = s.split('-') start_int = int(start) end_int = int(end) + 1 # so that the last element in range is also included ids_list_from_range = list(range(start_int, end_int)) parsed_ids_as_str.extend([str(id).zfill(zfill_value) for id in ids_list_from_range]) else: parsed_ids_as_str.append(s) parsed_ids_as_str = sorted(set(parsed_ids_as_str)) return parsed_ids_as_str
[docs] def init_siibra_params(atlas_name, parcellation_name, cohort_name, subject_ids): """ Initialize siibra parameters (if some were not given) and check the compatibility of the provided parameters. :param: atlas_name - name of atlas as str :param: parcellation_name - name of parcellation as str :param: cohort_name - name of cohort as str :param: subject_ids - list of unparsed subject ids given as str :return: (atlas, parcellation, cohort_name, subject_ids) - tuple containing a siibra atlas object, a siibra parcellation object and a cohort name, all compatible with each other, and a list of parsed ids """ # check that the atlas and the parcellation exist within siibra atlas = siibra.atlases[atlas_name] if atlas_name else None parcellation = siibra.parcellations[parcellation_name] if parcellation_name else None # check compatibility of atlas and parcellation if atlas and parcellation: compatible = check_atlas_parcellation_compatible(atlas, parcellation) if not compatible: raise ValueError(f'Atlas \'{atlas.name}\' does not contain parcellation \'{parcellation.name}\'. ' f'Please choose a different atlas and/or parcellation.') if atlas and not parcellation: LOGGER.warning(f'No parcellation was provided, so a default one will be selected.') parcellations = get_parcellations_for_atlas(atlas) no_parcellations = len(parcellations) if no_parcellations < 1: raise AttributeError(f'No default parcellation was found for atlas {atlas.name}!') elif no_parcellations > 1: LOGGER.info( f'Multiple parcellations were found for atlas {atlas.name}. An arbitrary one will be selected.') # select the newest parcellation version parcellation = [p for p in parcellations if p.is_newest_version][0] if not atlas and parcellation: LOGGER.warning('A parcellation was provided without an atlas, so a default atlas will be selected.') atlases = get_atlases_for_parcellation(parcellation) no_atlases = len(atlases) if no_atlases < 1: raise AttributeError(f'No default atlas containing parcellation {parcellation.name} was found!') elif no_atlases > 1: LOGGER.info( f'Multiple atlases containing parcellation {parcellation_name} were found. ' f'An arbitrary one will be selected') atlas = atlases[0] if not atlas and not parcellation: LOGGER.warning(f'No atlas and no parcellation were provided, so default ones will be selected.') atlas = siibra.atlases[HUMAN_ATLAS] parcellation = siibra.parcellations[JULICH_3_0] LOGGER.info(f'Using atlas {atlas.name} and parcellation {parcellation.name}') if cohort_name: # check the compatibility of cohort and parcellation cohort_options = get_cohorts_for_sc(parcellation.name) if cohort_name not in cohort_options: raise ValueError(f'The cohort \"{cohort_name}\" is not available for parcellation \"{parcellation.name}\"! ' f'Please choose one of the following cohorts: {cohort_options} or change ' f'the parcellation.') else: cohort_name = HCP_COHORT # compatible with all parcellations # check subject ids if not subject_ids: raise ValueError(f'Please provide at least one subject ID!') else: subject_ids = parse_subject_ids(subject_ids, cohort_name) return atlas, parcellation, cohort_name, subject_ids
# ############################################# CONNECTIVITY METHODS ###################################################
[docs] def get_hemispheres_for_regions(region_names): """ Given a list of region names, compute the hemispheres to which they belon to. 0 means the region belongs to the left hemisphere, 1 means it belongs to the right hemisphere (according to TVB convention). :param: region_names - list of str containing the names of the regions :return: hemi - list of ints indicating the hemisphere for each region in `region_names` """ LOGGER.info(f'Computing hemispheres for regions') hemi = [] for name in region_names: if 'right' in name: hemi.append(1) else: hemi.append(0) return hemi
[docs] def get_regions_positions(regions): """ Given a list of siibra regions, compute the positions of their centroids. :param: regions - list of siibra Regions :return: positions - list of tuples; each tuple represents the position of a region and contains 3 floating point coordinates """ LOGGER.info(f'Computing positions for regions') positions = [] for r in regions: space = r.supported_spaces[0] # choose first space that is available for that region centroid = r.spatial_props(space=space).components[0].centroid positions.append(centroid.coordinate) return positions
# ###################################### STRUCTURAL CONNECTIVITY METHODS ###############################################
[docs] def get_connectivity_matrix(parcellation, cohort, subjects, component): # type: (siibra.core.parcellation.Parcellation, str, list, Component2Modality) -> {} """ Retrieve the structural connectivity components (weights/tracts) for all the subjects provided, for the specified parcellation and cohort. The matrices are returned inside a dictionary, where the keys are the subject ids and the values represent the connectivity matrix. :param: parcellation - siibra Parcellation object for which we compute the connectivity matrices :param: cohort - name of cohort for which we compute the connectivity matrices :param: subjects - list containing the subject ids as strings :param: component - enum value specifying the connectivity component we want, weights or tracts return: conn_matrices - dict where key is the subject id and value is the conn. matrix """ modality = component.value features = siibra.features.get(parcellation, modality) conn_for_cohort = None # conn. obj. (StreamlineCounts/StreamlineLengths) for specified cohort conn_matrices = {} # dict containing connectivity matrices for the specified users # select the connectivities for the specific cohort for f in features: if f.cohort == cohort: conn_for_cohort = f break if conn_for_cohort is None: raise AttributeError(f"NO {modality} was found for cohort {cohort}") # for 1000BRAINS cohort, if the user did not specify a suffix (_1, _2), get all the possible ids for that subject if cohort == THOUSAND_BRAINS_COHORT and subjects is not None: all_subjects = sorted([conn.subject for conn in conn_for_cohort.elements]) subjects = [s for s in all_subjects if any(x in s for x in subjects)] # get the conn. matrices for each specified subject for s in subjects: conn = [c for c in conn_for_cohort if c.subject == s][0] matrix = conn.data conn_matrices[s] = matrix LOGGER.info(f'{component.name} for subject {s} retrieved successfully.') return conn_matrices
[docs] def create_tvb_structural_connectivity(weights_matrix, tracts_matrix, region_names, hemispheres, positions): # type: (pandas.DataFrame, pandas.DataFrame, list, list, list) -> tvb.datatypes.connectivity.Connectivity """ Create and configure a TVB Connectivity, based on its components obtained from siibra. :param: weights_matrix - pandas.DataFrame matrix for weights; obtained from siibra :param: tracts_matrix - pandas.DataFrame matrix for tracts; obtained from siibra :param: region_names - list of str containing the names of the regions for which the connectivity is computed :param: hemispheres - list of ints, corresponding to the hemisphere that each region from `region_names belongs to :param: positions - list of tuples, corresponding to region coordinates for each region from `region_names` :return: conn - a tvb.datatypes.connectivity.Connectivity object """ conn = connectivity.Connectivity() conn.weights = weights_matrix.to_numpy() conn.tract_lengths = tracts_matrix.to_numpy() conn.region_labels = np.array(region_names) conn.hemispheres = np.array(hemispheres, dtype=np.bool_) conn.centres = np.array(positions) conn.configure() return conn
[docs] def get_structural_connectivities_from_kg(atlas=None, parcellation=None, cohort=None, subject_ids=None): """ Return a dictionary of TVB Structural Connectivities using data from siibra and the KG, based on the specified atlas, parcelation and cohort, and for the specified subjects :param: atlas - str specifying the atlas name :param: parcellation - str specifying the parcellation name :param: cohort - str specifying the cohort name :param: subject_ids - unparsed str specifying the subject ids for which the connectivities will be retrived :return: connectivities - dict containing tvb structural Connectivities as values and the subject ids as keys """ atlas, parcellation, cohort, subject_ids = init_siibra_params(atlas, parcellation, cohort, subject_ids) connectivities = {} weights = get_connectivity_matrix(parcellation, cohort, subject_ids, Component2Modality.WEIGHTS) tracts = get_connectivity_matrix(parcellation, cohort, subject_ids, Component2Modality.TRACTS) # regions are the same for all weights and tract lengths matrices, so they can be computed only once regions = list(weights.values())[0].index.values region_names = [r.name for r in regions] hemi = get_hemispheres_for_regions(region_names) positions = get_regions_positions(regions) LOGGER.info(f'Computing TVB Connectivities') for subject, matrix in weights.items(): weights_matrix = matrix tracts_matrix = tracts[subject] tvb_conn = create_tvb_structural_connectivity(weights_matrix, tracts_matrix, region_names, hemi, positions) # structural connectivities stored as dict, where key is subject id, as we need it when computing connectivity # measures connectivities[subject] = tvb_conn return connectivities
# #################################### FUNCTIONAL CONNECTIVITY METHODS #################################################
[docs] def get_functional_connectivity_matrix(parcellation, cohort, subject): """ Get all the functional connectivities for the specified parcellation, cohort and just ONE specific subject; In v1.0a5 of siibra, for HCP cohort there are 5 groups of functional connectivities; each group contains 1 Functional connectivity for each subject addressed in the research :param: parcellation - siibra Parcellation object :param: cohort - str specifying the cohort name :param: subject - str specifying exactly one subject id :return: (fcs_list, fcs_names_list) - tuple containing 2 lists; `fcs_list` contains pandas.Dataframe matrices and `fcs_names_list` contains the name for each matrix from the previous list, obtained from the file they are stored in the KG """ modality = Component2Modality.FUNCTIONAL_CONNECTIVITY.value features = siibra.features.get(parcellation, modality) fcs_list = [] # the FC matrices fcs_names_list = [] # the name of the files containing the FC matrices in fcs_list for f in features: if f.cohort == cohort: f_conn = [c for c in f.elements if c.subject == subject][0] fc_matrix = f_conn.data fcs_names_list.append(f_conn.name) fcs_list.append(fc_matrix) return fcs_list, fcs_names_list
[docs] def create_tvb_connectivity_measure(siibra_fc, structural_connectivity, siibra_fc_name): """ Given a FunctionalConnectivity from siibra TVB Structural Connectivity (both for the same subject), return a TVB ConnectivityMeasure containing those 2 connectivities :param: siibra_fc - pandas.Dataframe matrix from siibra containing a functional connectivity :param: structural_connectivity - a TVB structural connectivity :param: siibra_fc_name - the name of the siibra functional connectivity object :return: conn_measure - tvb.datatypes.graph.ConnectivityMeasure representing a functional connectivity """ fc_matrix = siibra_fc.to_numpy() conn_measure = ConnectivityMeasure(array_data=fc_matrix, connectivity=structural_connectivity) conn_measure.title = siibra_fc_name return conn_measure
[docs] def get_connectivity_measures_from_kg(atlas=None, parcellation=None, cohort=None, subject_ids=None, structural_connectivities=None): """ Return a dictionary of TVB Connectivity Measures using data from siibra and the KG, based on the specified atlas, parcelation and cohort, and for the specified subjects :param: atlas - str specifying the atlas name :param: parcellation - str specifying the parcellation name :param: cohort - str specifying the cohort name :param: subject_ids - unparsed str specifying the subject ids for which the connectivities will be retrieved :param: structural_connectivities - dict of TVB Structural Connectivities computed for the subjects from `subject_ids`, where subject ids are keys and the structural connectivities are values :return: conn_measures - dict containing TVB Connectivity Measures as values and the subject ids as keys """ atlas, parcellation, cohort, subject_ids = init_siibra_params(atlas, parcellation, cohort, subject_ids) conn_measures = {} # for 1000BRAINS cohort, if the user did not specify a suffix (_1, _2), get all the possible ids for that subject if cohort == THOUSAND_BRAINS_COHORT and any('_' not in s for s in subject_ids): f = [f for f in siibra.features.get(parcellation, siibra.features.connectivity.FunctionalConnectivity) if f.cohort == THOUSAND_BRAINS_COHORT][0] # there is only one FC group for this 1000BRAINS subjects_for_cohort = [f_conn.subject for f_conn in f] subject_ids = [s for s in sorted(subjects_for_cohort) if any(x in s for x in subject_ids)] for s in subject_ids: conn_measures[s] = [] sc = structural_connectivities[s] fcs, fcs_names = get_functional_connectivity_matrix(parcellation, cohort, s) # create a tvb conn measure from a siibra fc and append it in the final dict for the specified subject for i in range(len(fcs)): conn_measure = create_tvb_connectivity_measure(fcs[i], sc, fcs_names[i]) conn_measures[s].append(conn_measure) return conn_measures
# ################################################# FINAL API ##########################################################
[docs] def get_connectivities_from_kg(atlas=None, parcellation=None, cohort=HCP_COHORT, subject_ids=None, compute_fc=False): """ Compute the TVB Structural Connectivities and optionally Functional Connectivities for the selected subjects :param: atlas - str specifying the atlas name :param: parcellation - str specifying the parcellation name :param: cohort - str specifying the cohort name :param: subject_ids - unparsed str specifying the subject ids for which the connectivities will be retrieved :param: compute_fc - boolean value indicating if for the specified subjects the functional connectivities should also be retrieved :return: (sc_dict, conn_measures_dict) - tuple containing 2 dictionaries: one with structural connectivities and one for functional connectivities; for each dictionary, the keys are the subject ids and the values are the connectivities """ conn_measures_dict = {} sc_dict = get_structural_connectivities_from_kg(atlas, parcellation, cohort, subject_ids) if compute_fc: conn_measures_dict = get_connectivity_measures_from_kg(atlas, parcellation, cohort, subject_ids, sc_dict) return sc_dict, conn_measures_dict