Extraction of parcellated fMRI timeseries from the hierarchical downsampled DSURQE atlas

The following function will extract parcellated timeseries from a cleaned EPI file after confound regression and labeling of the image using the downsampled DSURQE atlas (/data/chamal/projects/Gabriel_DG/atlases/DSURQE_atlas/hierarchical_downsample/DSURQE_downsample.nii.gz). It will output a python dictionary with the label anatomical names and their associated timeseries.

def extract_timeseries(cleaned_file, atlas_path, out_name, label_filename='/data/chamal/projects/Gabriel_DG/atlases/DSURQE_atlas/hierarchical_downsample/DSURQE_downsample_labels.pkl', ABI_info='/data/chamal/projects/Gabriel_DG/atlases/AllenBrain/structures.csv', roi_exclusion_list=[]):
    #cleaned_file: path to cleaned EPI file after confound regression
    #atlas_path: path to atlas files directory, or path to a commonspace file
    #out_name: output file name, must end with .pkl
    #label_filename: .pkl file with info about the labels
    #ABI_info: .csv with ABI annotations
    #roi_exclusion_list: list of roi labels to omit
    import os
    import numpy as np
    import nibabel as nb
    from nilearn.input_data import NiftiMasker

    #prepare DSURQE downsampled atlas ids and ROI names
    import pickle
    with open(label_filename,'rb') as handle:
        labels=pickle.load(handle)

    roi_names=labels['roi']
    rh_ids=labels['rh']
    lh_ids=labels['lh']

    import pandas as pd
    ABI_df=pd.read_csv(ABI_info)
    ABI_names=list(ABI_df['name'])
    ABI_acronyms=list(ABI_df['acronym'])

    #replace the acronyms by the full names
    full_names=[]
    for roi_name in roi_names:
        if roi_name=='ILA-PL':
            full_names.append('Infra/Prelimbic area')
        elif roi_name=='AI-ORB':
            full_names.append('Orbital and Agranular insular areas')
        else:
            index=ABI_acronyms.index(roi_name)
            full_names.append(ABI_names[index])

    #determine indices to exclude based on roi_exclusion_list
    rh_rm_list=[]
    lh_rm_list=[]
    name_rm_list=[]
    for roi_name, rh, lh in zip(full_names, rh_ids, lh_ids):
        if rh in roi_exclusion_list or lh in roi_exclusion_list:
            print('removing '+roi_name)
            rh_rm_list.append(rh)
            lh_rm_list.append(lh)
            name_rm_list.append(roi_name)
    for roi_name, rh, lh in zip(name_rm_list, rh_rm_list, lh_rm_list):
        rh_ids.remove(rh)
        lh_ids.remove(lh)
        full_names.remove(roi_name)

    data_dict={}
    for roi_name, rh, lh in zip(full_names, rh_ids, lh_ids):
        exclude=False
        bilateral=False

        atlas=atlas_path
        atlas_img=nb.load(atlas)
        atlas_data=np.asarray(atlas_img.dataobj)
        if np.max(rh==atlas_data): #taking a ROI only if it has labeled voxels
            roi_mask=np.asarray(atlas_data==rh, dtype=int)
            roi_mask_nb=nb.Nifti1Image(roi_mask, nb.load(atlas).affine, nb.load(atlas).header)
            masker = NiftiMasker(mask_img=roi_mask_nb, standardize=False, verbose=0)
            voxel_timeseries = masker.fit_transform(cleaned_file) #extract the voxel timeseries within the mask
            roi_timeseries=np.mean(voxel_timeseries, axis=1) #take the mean ROI timeseries
            if lh==rh:
                label=roi_name
            else:
                bilateral=True
                label=roi_name+'_rh'
                if np.max(lh==atlas_data): #taking a ROI only if it has labeled voxels
                    roi_mask=np.asarray(atlas_data==lh, dtype=int)
                    roi_mask_nb=nb.Nifti1Image(roi_mask, nb.load(atlas).affine, nb.load(atlas).header)
                    masker = NiftiMasker(mask_img=roi_mask_nb, standardize=False, verbose=0)
                    voxel_timeseries = masker.fit_transform(cleaned_file) #extract the voxel timeseries within the mask
                    roi_timeseries_lh=np.mean(voxel_timeseries, axis=1) #take the mean ROI timeseries
                    lh_label=roi_name+'_lh'
                else:
                    exclude=True
                    print(roi_name+" excluded.")
        else:
            exclude=True
            print(roi_name+" excluded.")
        if exclude:
            continue #pass to the next label if it has to be excluded

        data_dict[str(label)]=roi_timeseries
        print(str(label))
        if bilateral:
            data_dict[str(lh_label)]=roi_timeseries_lh
            print(str(lh_label))

    import pickle

    with open(out_name, 'wb') as handle:
        pickle.dump(data_dict, handle, protocol=pickle.HIGHEST_PROTOCOL)
    return data_dict

Example code to derive a functional connectivity matrix from the parcellated timeseries:

import numpy as np

sub_timeseries = extract_timeseries('sub_cleaned_bold.nii.gz', 'sub_labels.nii.gz', 'sub_timeseries.pkl', label_filename='/data/chamal/projects/Gabriel_DG/atlases/DSURQE_atlas/hierarchical_downsample/DSURQE_downsample_labels.pkl', ABI_info='/data/chamal/projects/Gabriel_DG/atlases/AllenBrain/structures.csv')

roi_labels=sub_timeseries.keys()

#creating a 2D numpy array of ROI by timeseries
all_timeseries=[]
for roi in roi_labels:
    all_timeseries.append(sub_timeseries[roi])
all_timeseries=np.asarray(all_timeseries)

#derive the correlation matrix by calculating the pairwise correlation between every ROI timeseries
FC_matrix=np.corrcoef(all_timeseries)

#visualize the matrix
import matplotlib.pyplot as plt
plt.imshow(FC_matrix)

Using the python_FC_pkg (https://github.com/Gab-D-G/python_FC_pkg):

from python_FC_pkg import fc_utils
import numpy as np

sub_timeseries = extract_timeseries('sub_cleaned_bold.nii.gz', 'sub_labels.nii.gz', 'sub_timeseries.pkl', label_filename='/data/chamal/projects/Gabriel_DG/atlases/DSURQE_atlas/hierarchical_downsample/DSURQE_downsample_labels.pkl', ABI_info='/data/chamal/projects/Gabriel_DG/atlases/AllenBrain/structures.csv')

roi_labels=sub_timeseries.keys()

#creating a 2D numpy array of ROI by timeseries
all_timeseries=[]
for roi in roi_labels:
    all_timeseries.append(sub_timeseries[roi])
all_timeseries=np.asarray(all_timeseries)

#derive the correlation matrix by calculating the pairwise correlation between every ROI timeseries
FC_vector=fc_utils.get_FC_vector(all_timeseries.T, metric='r')
FC_matrix=fc_utils.recover_matrix(FC_vector, len(roi_labels))

#visualize the matrix
import matplotlib.pyplot as plt
plt.imshow(FC_matrix)