Dog brain representations of human facial expressions: Encoding happiness and differentiating specific negative expressions.
The 5 matches
- [1] § STAR★Methods › Method details › Experiment 2 › Quantification and statistical analysis ↔ preprocess_functions.py, lines 1524–1578 · score 0.70 · binary mask, FSL, FD, displacement, rotation, motion
- [2] § STAR★Methods › Method details › Experiment 2 › Quantification and statistical analysis ↔ tools/fsf_design.py, lines 1–50 · score 0.65 · convolving, swapped, FSL, Tool, split, filtered
- [3] § STAR★Methods › Method details › Experiment 1 › Quantification and statistical analysis ↔ utils.py, lines 179–313 · score 0.57 · anatomical image, FSL, kernel, smoothing, preprocessing, fMRI
- [4] § Results › Experiment 1: Happy faces elicit stronger right temporal-caudate responses than neutral faces in dogs ↔ export_static.py, lines 240–315 · score 0.52 · Sylvian gyrus, brain response, caudal, temporal, radius, dog
- [5] § STAR★Methods › Method details › Experiment 1 › Quantification and statistical analysis ↔ convert_coords.ipynb, lines 394–457 · score 0.51 · Dog Brain Toolkit, FSL, smoothing, atlas, fMRI
Paper
Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC
The paper is loaded when this pane is shown.
The authors' code
Python · 1,594 lines · 72 KB · no license · 1 match
- import os
- import utils
- import shutil
- from nilearn.plotting import plot_anat, show
- import ipywidgets as widgets
- import utils
- import os
- from importlib import reload
- from nilearn import plotting
- from nilearn import image as nli
- import nibabel as nib
- import numpy as np
- from ipywidgets import HBox, VBox
- import numpy as np
- from nilearn.plotting import plot_anat, show
- import ipywidgets as widgets
- import matplotlib.pyplot as plt
- from matplotlib.colors import ListedColormap
- from IPython.display import display, Markdown, clear_output
- # Description: Functions to preprocess fMRI data using FSL
- # Author: Raul Hernandez
- def bet_app(project_dict, sub_N, initial_params=None):
- global data_mask1, data_mask2, data_mask3
- # Create a custom colormaps
- # black
- black = np.array([0, 0, 0, 0]) # RGBA for black (last value is alpha)
- # red
- red = np.array([1, 0, 0, 1]) # RGBA for red
- # blue
- blue = np.array([0, 0, 1, 1]) # RGBA for blue
- # yellow
- yellow = np.array([1, 1, 0, 1]) # RGBA for yellow (last value is alpha)
- # create a red colormap
- red_cmap = ListedColormap([black, red])
- # create a blue colormap
- blue_cmap = ListedColormap([black, blue])
- # create a yellow colormap
- yellow_cmap = ListedColormap([black, yellow])
- dataset = project_dict['Dataset']
- session = project_dict['Session']
- task = project_dict['Task']
- specie = project_dict['Specie']
- datafolder = project_dict['Datafolder']
- # working directory
- workingdir = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2)
- mean_fct_file = workingdir + os.sep + specie + '-sub-' + str(sub_N).zfill(2)
- # adding session if there is one
- if session != '':
- mean_fct_file += '_ses-' + session
- cut_mean_fct_file = mean_fct_file + '_task-' + task + '_mean_fct.nii.gz'
- mask_file = mean_fct_file + '_task-' + task + '_mean_fct_mask.nii.gz'
- mask_file1 = mean_fct_file + '_task-' + task + '_mean_fct_mask1.nii.gz'
- mask_file2 = mean_fct_file + '_task-' + task + '_mean_fct_mask2.nii.gz'
- mask_file3 = mean_fct_file + '_task-' + task + '_mean_fct_mask3.nii.gz'
- mean_fct_file += '_task-' + task + '_mean_fct.nii.gz'
- masked_mean_fct_file = mean_fct_file[:-7] + '_brain.nii.gz'
- # load cut_mean_fct_file
- img = nib.load(cut_mean_fct_file)
- # Get the data from the image
- data = img.get_fdata()
- # if the mask files exist, load them
- if os.path.exists(mask_file1):
- mask_img1 = nib.load(mask_file1)
- data_mask1 = mask_img1.get_fdata()
- mask_img2 = nib.load(mask_file2)
- data_mask2 = mask_img2.get_fdata()
- mask_img3 = nib.load(mask_file3)
- data_mask3 = mask_img3.get_fdata()
- #params_file = mean_fct_file[:-7] + '_bet_params.txt'
- #param_dict = utils.read_params_file(params_file)
- # determine min and max values for each axis
- maxX,maxY,maxZ = img.shape
- def plot_bet(x,y,z, betx_1, bety_1, betz_1, betx_2, bety_2, betz_2, betx_3, bety_3, betz_3, plot_mask):
- global data_mask1, data_mask2, data_mask3
- """
- Plot slices from the sagittal, coronal, and axial views side by side.
- """
- fig, axes = plt.subplots(1, 3, figsize=(8, 5))
- # Sagittal
- sagittal_slice = data[x, :, :]
- axes[0].imshow(sagittal_slice.T, cmap='gray', origin='lower')
- axes[0].scatter(bety_1,betz_1,s=200, c='red')
- # plot blue only if betx_2 matches the current slice, else plot it with alpha=0.4
- if betx_2 == x:
- axes[0].scatter(bety_2,betz_2,s=200, c='blue')
- else:
- axes[0].scatter(bety_2,betz_2,s=200, c='blue', alpha=0.4)
- axes[0].scatter(bety_3,betz_3,s=200, c='yellow')
- axes[0].axis('off')
- # display mask if plot_mask is True
- if plot_mask:
- # get mask1 for sagital slice
- sagital_mask1 = data_mask1[x, :, :]
- # plot mask1 in red
- axes[0].imshow(sagital_mask1.T, cmap=red_cmap, origin='lower', alpha=0.5)
- # get mask2 for sagital slice
- sagital_mask2 = data_mask2[x, :, :]
- # plot mask2 in blue
- axes[0].imshow(sagital_mask2.T, cmap=blue_cmap, origin='lower', alpha=0.5)
- # get mask3 for sagital slice
- sagital_mask3 = data_mask3[x, :, :]
- # plot mask3 in yellow
- axes[0].imshow(sagital_mask3.T, cmap=yellow_cmap, origin='lower', alpha=0.5)
- # Coronal
- coronal_slice = data[:, y, :]
- axes[1].imshow(coronal_slice.T, cmap='gray', origin='lower')
- # plot red only if bety_1 matches the current slice, else plot it with alpha=0.4
- if bety_1 == y:
- axes[1].scatter(betx_1,betz_1,s=200, c='red')
- else:
- axes[1].scatter(betx_1,betz_1,s=200, c='red', alpha=0.4)
- axes[1].scatter(betx_2,betz_2,s=200, c='blue')
- axes[1].scatter(betx_3,betz_3,s=200, c='yellow')
- axes[1].axis('off')
- # display mask if plot_mask is True
- if plot_mask:
- # get mask1 for coronal slice
- coronal_mask1 = data_mask1[:, y, :]
- # plot mask1 in red
- axes[1].imshow(coronal_mask1.T, cmap=red_cmap, origin='lower', alpha=0.5)
- # get mask2 for coronal slice
- coronal_mask2 = data_mask2[:, y, :]
- # plot mask2 in blue
- axes[1].imshow(coronal_mask2.T, cmap=blue_cmap, origin='lower', alpha=0.5)
- # get mask3 for coronal slice
- coronal_mask3 = data_mask3[:, y, :]
- # plot mask3 in yellow
- axes[1].imshow(coronal_mask3.T, cmap=yellow_cmap, origin='lower', alpha=0.5)
- # Axial
- axial_slice = data[:, :, z]
- axes[2].imshow(axial_slice.T, cmap='gray', origin='lower')
- axes[2].scatter(betx_1,bety_1,s=200, c='red')
- axes[2].scatter(betx_2,bety_2,s=200, c='blue')
- # plot yellow only if betz_3 matches the current slice, else plot it with alpha=0.4
- if betz_3 == z:
- axes[2].scatter(betx_3,bety_3,s=200, c='yellow')
- else:
- axes[2].scatter(betx_3,bety_3,s=200, c='yellow', alpha=0.4)
- axes[2].axis('off')
- # display mask if plot_mask is True
- if plot_mask:
- # get mask1 for axial slice
- axial_mask1 = data_mask1[:, :, z]
- # plot mask1 in red
- axes[2].imshow(axial_mask1.T, cmap=red_cmap, origin='lower', alpha=0.5)
- # get mask2 for axial slice
- axial_mask2 = data_mask2[:, :, z]
- # plot mask2 in blue
- axes[2].imshow(axial_mask2.T, cmap=blue_cmap, origin='lower', alpha=0.5)
- # get mask3 for axial slice
- axial_mask3 = data_mask3[:, :, z]
- # plot mask3 in yellow
- axes[2].imshow(axial_mask3.T, cmap=yellow_cmap, origin='lower', alpha=0.5)
- # make a tight layout
- plt.tight_layout()
- # make the background black
- fig.patch.set_facecolor('black')
- # display the plot
- plt.show()
- # define function to apply BET
- def apply_bet_button(betx_1, bety_1, betz_1, thr_1, betx_2, bety_2, betz_2, thr_2, betx_3, bety_3, betz_3, thr_3):
- global data_mask1, data_mask2, data_mask3
- # determine name of params_file based on mean_fct_file
- params_file = mean_fct_file[:-7] + '_bet_params.txt'
- # create a dict with the parameters
- param_dict = {'betx_1':betx_1, 'bety_1':bety_1, 'betz_1':betz_1, 'thr_1':thr_1,
- 'betx_2':betx_2, 'bety_2':bety_2, 'betz_2':betz_2, 'thr_2':thr_2,
- 'betx_3':betx_3, 'bety_3':bety_3, 'betz_3':betz_3, 'thr_3':thr_3,
- 'output_file':mask_file, 'masked_file':masked_mean_fct_file,
- 'mask_file1':mask_file1, 'mask_file2':mask_file2, 'mask_file3':mask_file3,
- }
- # save the cutting parameters
- utils.write_params_file(params_file, param_dict)
- print('button pressed')
- command = f"./run_bet.sh {cut_mean_fct_file} {params_file}"
- os.system(command)
- # update button description and status
- col1.children[3].description = 'Display mask'
- col1.children[3].disabled = False
- # Load each of the three masks
- mask_img1 = nib.load(mask_file1)
- data_mask1 = mask_img1.get_fdata()
- mask_img2 = nib.load(mask_file2)
- data_mask2 = mask_img2.get_fdata()
- mask_img3 = nib.load(mask_file3)
- data_mask3 = mask_img3.get_fdata()
- print('Done')
- # create 4 sets of sliders, one for the slice, and 3 for the spheres marking the initial place of the BET sphere
- # check if mask file 1 already exist
- if os.path.exists(mask_file1):
- button_description = 'Display mask'
- button_disabled = False
- else:
- button_description = 'Mask not available'
- button_disabled = True
- # sliders for the slice
- col1 = widgets.VBox([widgets.IntSlider(min=0, max=maxX, step=1, value=np.round(maxX/2), description='X'),
- widgets.IntSlider(min=0, max=maxY, step=1, value=np.round(maxY/2), description='Y'),
- widgets.IntSlider(min=0, max=maxZ, step=1, value=np.round(maxZ/2), description='Z'),
- widgets.ToggleButton(value=False, description=button_description, button_style='info', disabled=button_disabled),
- ])
- # check if param_dict exists if yes, load it
- if os.path.exists(mean_fct_file[:-7] + '_bet_params.txt'):
- param_dict = utils.read_params_file(mean_fct_file[:-7] + '_bet_params.txt')
- else:
- if initial_params is not None:
- param_dict = {'betx_1':np.round(maxX/2),
- 'bety_1':initial_params['bety_1'],
- 'betz_1':initial_params['betz_1'],
- 'thr_1':initial_params['thr_1'],
- 'betx_2':np.round(maxX/2),
- 'bety_2':initial_params['bety_2'],
- 'betz_2':initial_params['betz_2'],
- 'thr_2':initial_params['thr_2'],
- 'betx_3':np.round(maxX/2),
- 'bety_3':initial_params['bety_3'],
- 'betz_3':initial_params['betz_3'],
- 'thr_3':initial_params['thr_3'],
- }
- else:
- param_dict = {'betx_1':np.round(maxX/2),
- 'bety_1':48,
- 'betz_1':11,
- 'thr_1':0.65,
- 'betx_2':np.round(maxX/2),
- 'bety_2':29,
- 'betz_2':11,
- 'thr_2':0.6,
- 'betx_3':np.round(maxX/2),
- 'bety_3':56,
- 'betz_3':2,
- 'thr_3':0.75,
- }
- # sliders for the spheres
- # first sphere, has 4 values, x,y,z and threshold
- col2 = widgets.VBox([widgets.IntSlider(min=0, max=maxX, step=1, value=param_dict['betx_1'], description='betx_1'),
- widgets.IntSlider(min=0, max=maxY, step=1, value=param_dict['bety_1'], description='bety_1'),
- widgets.IntSlider(min=0, max=maxZ, step=1, value=param_dict['betz_1'], description='betz_1'),
- widgets.FloatSlider(min=0, max=1, step=0.05, value=param_dict['thr_1'], description='thr_1'),
- ])
- # second sphere, has 4 values, x,y,z and threshold
- col3 = widgets.VBox([widgets.IntSlider(min=0, max=maxX, step=1, value=param_dict['betx_2'], description='betx_2'),
- widgets.IntSlider(min=0, max=maxY, step=1, value=param_dict['bety_2'], description='bety_2'),
- widgets.IntSlider(min=0, max=maxZ, step=1, value=param_dict['betz_2'], description='betz_2'),
- widgets.FloatSlider(min=0, max=1, step=0.05, value=param_dict['thr_2'], description='thr_2'),
- ])
- # third sphere, has 4 values, x,y,z and threshold
- col4 = widgets.VBox([widgets.IntSlider(min=0, max=maxX, step=1, value=param_dict['betx_3'], description='betx_3'),
- widgets.IntSlider(min=0, max=maxY, step=1, value=param_dict['bety_3'], description='bety_3'),
- widgets.IntSlider(min=0, max=maxZ, step=1, value=param_dict['betz_3'], description='betz_3'),
- widgets.FloatSlider(min=0, max=1, step=0.05, value=param_dict['thr_3'], description='thr_3'),
- ])
- # linkin the sliders to the function
- bet_app_out = widgets.interactive_output(plot_bet, {'x':col1.children[0], 'y':col1.children[1], 'z':col1.children[2],
- 'betx_1':col2.children[0], 'bety_1':col2.children[1], 'betz_1':col2.children[2],
- 'betx_2':col3.children[0], 'bety_2':col3.children[1], 'betz_2':col3.children[2],
- 'betx_3':col4.children[0], 'bety_3':col4.children[1], 'betz_3':col4.children[2],
- 'plot_mask':col1.children[3],
- })
- # button to apply the BET
- col5 = widgets.Button(description='Apply BET')
- # setting the button to call the function
- col5.on_click(lambda b: apply_bet_button(
- col2.children[0].value, col2.children[1].value, col2.children[2].value, col2.children[3].value,
- col3.children[0].value, col3.children[1].value, col3.children[2].value, col3.children[3].value,
- col4.children[0].value, col4.children[1].value, col4.children[2].value, col4.children[3].value,
- ))
- row1 = HBox([col1,col2,col3,col4])
- row2 = HBox([col5])
- bet_app_tab = VBox([row1,row2])
- return bet_app_tab, bet_app_out
- def crop_app(project_dict, sub_N):
- def plot_slices(x, y, z, x_lim1, y_lim1, z_lim1, x_lim2, y_lim2, z_lim2):
- """
- Plot slices from the sagittal, coronal, and axial views side by side.
- """
- # get the max in x
- maxX = data.shape[0]
- # flip x_lim1 and x_lim2
- x_lim1 = maxX - x_lim1
- x_lim2 = maxX - x_lim2
- fig, axes = plt.subplots(1, 3, figsize=(8, 5))
- # Sagittal
- sagittal_slice = data[x, :, :]
- axes[0].imshow(sagittal_slice.T, cmap='gray', origin='lower')
- axes[0].axis('off')
- # Coronal
- coronal_slice = data[:, y, :]
- axes[1].imshow(coronal_slice.T, cmap='gray', origin='lower')
- axes[1].axis('off')
- # Axial
- axial_slice = data[:, :, z]
- axes[2].imshow(axial_slice.T, cmap='gray', origin='lower')
- axes[2].axis('off')
- # make a tight layout
- plt.tight_layout()
- # make the background black
- fig.patch.set_facecolor('black')
- # Plotting red lines
- # plotting lines in sagital slice
- axes[0].axhline(y=z_lim1, color='red', lw=2)
- axes[0].axvline(x=y_lim1, color='red', lw=2)
- # plotting lines in coronal slice
- axes[1].axvline(x=x_lim1, color='red', lw=2)
- axes[1].axhline(y=z_lim1, color='red', lw=2)
- # plotting lines in axial slice
- axes[2].axvline(x=x_lim1, color='red', lw=2)
- axes[2].axhline(y=y_lim1, color='red', lw=2)
- # Plotting blue lines
- # plotting lines in sagital slice
- axes[0].axhline(y=z_lim2, color='blue', lw=2)
- axes[0].axvline(x=y_lim2, color='blue', lw=2)
- # plotting lines in coronal slice
- axes[1].axvline(x=x_lim2, color='blue', lw=2)
- axes[1].axhline(y=z_lim2, color='blue', lw=2)
- # plotting lines in axial slice
- axes[2].axvline(x=x_lim2, color='blue', lw=2)
- axes[2].axhline(y=y_lim2, color='blue', lw=2)
- plt.show()
- def apply_cut_button(x_lim1, y_lim1, z_lim1, x_lim2, y_lim2, z_lim2):
- # determine name of params_file based on mean_fct_file
- params_file = mean_fct_file[:-7] + '_cut_params.txt'
- # make sure that x_lim1 is smaller than x_lim2, if not invert them
- if x_lim1 > x_lim2:
- x_lim1, x_lim2 = x_lim2, x_lim1
- # make sure that y_lim1 is smaller than y_lim2, if not invert them
- if y_lim1 > y_lim2:
- y_lim1, y_lim2 = y_lim2, y_lim1
- # make sure that z_lim1 is smaller than z_lim2, if not invert them
- if z_lim1 > z_lim2:
- z_lim1, z_lim2 = z_lim2, z_lim1
- # create a dict with the cutting parameters
- param_dict = {'x_lim1':x_lim1, 'y_lim1':y_lim1,
- 'z_lim1':z_lim1, 'x_lim2':x_lim2,
- 'y_lim2':y_lim2, 'z_lim2':z_lim2,
- 'output_file':cut_mean_fct_file,
- }
- # save the cutting parameters
- utils.write_params_file(params_file, param_dict)
- print('button pressed')
- command = f"./remove_slices.sh {mean_fct_file} {params_file}"
- os.system(command)
- dataset = project_dict['Dataset']
- session = project_dict['Session']
- task = project_dict['Task']
- specie = project_dict['Specie']
- datafolder = project_dict['Datafolder']
- # working directory
- workingdir = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2)
- mean_fct_file = workingdir + os.sep + specie + '-sub-' + str(sub_N).zfill(2)
- # adding session if there is one
- if session != '':
- mean_fct_file += '_ses-' + session
- cut_mean_fct_file = mean_fct_file + '_task-' + task + '_mean_fct.nii.gz'
- mean_fct_file += '_task-' + task + '_mean_fct_uncut.nii.gz'
- # determine min and max values for each axis
- img = nib.load(mean_fct_file)
- # Get the data from the image
- data = img.get_fdata()
- slider_style = {'description_width': 'initial', 'width': '2px'}
- maxX,maxY,maxZ = img.shape
- col1 = widgets.VBox([widgets.IntSlider(min=0, max=maxX, step=1, value=np.round(maxX/2), description='X'),
- widgets.IntSlider(min=0, max=maxY, step=1, value=np.round(maxY/2), description='Y'),
- widgets.IntSlider(min=0, max=maxZ, step=1, value=np.round(maxZ/2), description='Z')])
- # check if param_dict exists if yes, load it
- if os.path.exists(mean_fct_file[:-7] + '_cut_params.txt'):
- param_dict = utils.read_params_file(mean_fct_file[:-7] + '_cut_params.txt')
- else:
- param_dict = {'x_lim1':0,
- 'y_lim1':0,
- 'z_lim1':0,
- 'x_lim2':maxX,
- 'y_lim2':maxY,
- 'z_lim2':maxZ,
- }
- col2 = widgets.VBox([widgets.IntSlider(min=0, max=maxX-1, step=1, value=param_dict['x_lim1'], description='lim X'),
- widgets.IntSlider(min=0, max=maxY-1, step=1, value=param_dict['y_lim1'], description='lim Y'),
- widgets.IntSlider(min=0, max=maxZ-1, step=1, value=param_dict['z_lim1'], description='lim Z')])
- col3 = widgets.VBox([widgets.IntSlider(min=0, max=maxX-1, step=1, value=param_dict['x_lim2'], description='lim X'),
- widgets.IntSlider(min=0, max=maxY-1, step=1, value=param_dict['y_lim2'], description='lim Y'),
- widgets.IntSlider(min=0, max=maxZ-1, step=1, value=param_dict['z_lim2'], description='lim Z')])
- # col2 = widgets.VBox([widgets.IntSlider(min=0, max=maxX-1, step=1, value= maxX - int(param_dict['x_lim2']), description='lim X'),
- # widgets.IntSlider(min=0, max=maxY-1, step=1, value=param_dict['y_lim1'], description='lim Y'),
- # widgets.IntSlider(min=0, max=maxZ-1, step=1, value=param_dict['z_lim1'], description='lim Z')])
- # col3 = widgets.VBox([widgets.IntSlider(min=0, max=maxX-1, step=1, value=maxX - int(param_dict['x_lim1']), description='lim X'),
- # widgets.IntSlider(min=0, max=maxY-1, step=1, value=param_dict['y_lim2'], description='lim Y'),
- # widgets.IntSlider(min=0, max=maxZ-1, step=1, value=param_dict['z_lim2'], description='lim Z')])
- col4 = widgets.Button(description='Apply cut')
- out = widgets.interactive_output(plot_slices, {'x':col1.children[0], 'y':col1.children[1], 'z':col1.children[2],
- 'x_lim1':col2.children[0], 'y_lim1':col2.children[1], 'z_lim1':col2.children[2],
- 'x_lim2':col3.children[0], 'y_lim2':col3.children[1], 'z_lim2':col3.children[2],
- })
- # setting the button to call the function
- col4.on_click(lambda b: apply_cut_button(
- col2.children[0].value, col2.children[1].value, col2.children[2].value,
- col3.children[0].value, col3.children[1].value, col3.children[2].value,
- ))
- tab_crop_app = VBox([HBox([col1,col2,col3]),col4])
- return tab_crop_app,out
- def check_job_status(job):
- # This function will check if the job finished, failed or is still running
- # Right now it will return 'Finished'
- job_status = 'Finished'
- return job_status
- def run_process(job):
- '''
- Will select the variables from the schedule_table and the project_dict to run the process
- Can run:
- preprocess_run
- get_mean_fct
- crop_interface
- bet_interface
- mean_to_STD
- run_to_STD
- '''
- # get the process to run
- # get the other variables of the job
- # user = job['User']
- dataset = job['Dataset']
- session = job['Session']
- task = job['Task']
- sub_N = job['sub_N']
- specie = job['Specie']
- process = job['Process']
- datafolder = job['Datafolder']
- job_status = job['Status']
- combination = job['Combination']
- smooth = job['Smooth']
- atlas_type = job['Atlas_type']
- base_run = 1
- img_type = 'brain2mm'
- run_prepro = job['Full_prepro']
- variation = job['Variation']
- session_and_run = job['session_and_run']
- first_time = job['first_time']
- use_anatomic = job['use_anatomic']
- print(f"Running process: {process}")
- # run the adecuate process
- if process == 'Preprocess':
- print('Running preprocess')
- for run_N, session in zip(job['run_N'], job['Sessions']):
- preprocess_run(
- sub_N, run_N, dataset, task, specie,
- datafolder, session, smooth,
- combination, run_prepro)
- elif process == 'Mean fct':
- print('Running get_mean_fct')
- runs_to_use = job['run_N']
- sessions_to_use = job['Sessions']
- session_and_run = job['session_and_run']
- get_mean_fct(
- sub_N, session_and_run, base_run, dataset,
- task, specie, datafolder, first_time=first_time)
- elif process == 'Mean to atlas':
- print('Running mean_to_STD')
- mean_to_STD(
- sub_N, dataset, task, specie, datafolder,
- atlas_type, img_type, variation=variation, use_anatomic=use_anatomic)
- elif process == 'Runs to atlas':
- print('Running run_to_STD')
- print('Variation:', variation)
- for run_N, session in zip(job['run_N'], job['Sessions']):
- run_to_STD(
- sub_N, run_N, dataset, task,
- specie, datafolder, atlas_type,
- img_type, session=session)
- else:
- print('Process not found')
- def check_file_status(project_dict, sub_N, run_N, session, process, verbose=False,
- model=None, dis_method=None, rsa_method=None, radius=None, rsa_model=None,
- stim_N=None, stim_1_name=None, stim_2_name=None, stim_types=None,
- reps=None, reps_group=None):
- '''
- Check which files are available for the process
- returns True if the files are available, False if not
- '''
- # print the inputs
- # print('Sub:', sub_N, 'Run:', run_N, 'Session:', session, 'Process:', process)
- # get specie from the project_dict
- specie = project_dict['Specie']
- # get dataset from the project_dict
- dataset = project_dict['Dataset']
- # get task from the project_dict
- task = project_dict['Task']
- # get datafolder from the project_dict
- datafolder = project_dict['Datafolder']
- if process == 'preprocess_run': # check if the preprocess files exist
- # create the filename
- if session != '':
- session = int(session)
- filename = (datafolder + os.sep + dataset + os.sep + 'BIDS' + os.sep +
- specie +'-sub-' + str(sub_N).zfill(2) + os.sep +
- specie + '-sub-' + str(sub_N).zfill(2) +
- '_ses-' + str(session).zfill(2) +
- '_task-' + task +
- '_run-' + str(run_N).zfill(2) +
- '_bold.nii.gz')
- filename_json = filename[:-7] + '.json'
- # check if filename and filename_json exist
- if os.path.exists(filename):
- if verbose:
- print('BIDS file exists: ' + filename)
- if os.path.exists(filename_json):
- if verbose:
- print('json file exists: ' + filename_json)
- return True, filename
- else:
- if verbose:
- print('BIDS file does not exist: ' + filename)
- return False, filename
- elif process == 'get_mean_fct':
- outputdir = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2)
- filename = specie + '-sub-' + str(sub_N).zfill(2)
- # adding session if there is one
- if session != '':
- session = int(session)
- filename += '_ses-' + f"{session:02d}"
- filename += '_task-' + task + '_run-' + str(run_N).zfill(2) + '_reoriented.nii.gz'
- # check if the file exists
- if os.path.exists(outputdir + os.sep + filename):
- if verbose:
- print('Reoriented file exists: ' + filename)
- return True, filename
- else:
- if verbose:
- print('Reoriented file does not exist: ' + filename)
- return False, filename
- elif process == 'mean_to_STD':
- # check if the mean file exists
- mean_fct_file = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2) + os.sep + specie + '-sub-' + str(sub_N).zfill(2)
- # adding session if there is one
- mean_fct_file += '_task-' + task + '_mean_fct_brain.nii.gz'
- # check if the file exists
- if os.path.exists(mean_fct_file):
- print('Mean functional file exists: ' + mean_fct_file)
- return True, mean_fct_file
- else:
- print('Mean functional file does not exist: ' + mean_fct_file)
- return False, mean_fct_file
- elif process == 'run_to_STD':
- preprocess_dir = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2)
- filename = (specie + '-sub-' + str(sub_N).zfill(2) +
- '_ses-' + session +
- '_task-' + task +
- '_run-' + str(run_N).zfill(2)
- )
- preprocessed_file = filename + '_mc.nii.gz'
- # check if the file exists
- if os.path.exists(preprocess_dir + os.sep + preprocessed_file):
- if verbose:
- print('Motion corrected file exists: ' + preprocessed_file)
- return True, preprocess_dir + os.sep + preprocessed_file
- else:
- if verbose:
- print('Motion corrected file does not exist: ' + preprocess_dir + os.sep + preprocessed_file)
- return False, preprocess_dir + os.sep + preprocessed_file
- elif process == 'BOLD': # check if the normalized BOLD file exist
- # if session is int convert to str with 2 digits
- if isinstance(session, int):
- session = str(session).zfill(2)
- filename = (datafolder + os.sep + dataset + os.sep + 'normalized' + os.sep +
- specie + '-sub-' + str(sub_N).zfill(2) + os.sep +
- specie + '-sub-' + str(sub_N).zfill(2) +
- '_ses-' + session +
- '_task-' + task +
- '_run-' + str(run_N).zfill(2) +
- '.nii.gz')
- # check if filename exists
- if os.path.exists(filename):
- if verbose:
- print('Normalized BOLD file exists: ' + filename)
- return True, filename
- else:
- if verbose:
- print('BOLD file does not exist: ' + filename)
- return False, filename
- elif process == 'feat_folder': # check if the first level GLM file exist
- # model = 'basic'
- # "P:\userdata\raulh87\data\EmoB\results\GLM\basic\D-sub-01\ses-01_task-EmoB_run-01.feat\stats\tstat1.nii.gz"
- # if session is int convert to str with 2 digits
- # check that all variables are available, if not indicate which is missing
- if model is None: # model is required
- print('Model is missing')
- return False, 'error in model'
- filename = (datafolder + os.sep + dataset + os.sep + 'results' + os.sep + 'GLM' + os.sep +
- model + os.sep + f"{specie}-sub-{sub_N:02d}" + os.sep + f"ses-{session}_task-{task}_run-{run_N:02d}.feat")
- # check if filename exist
- if os.path.exists(filename):
- if verbose:
- print('feat folder exists: ' + filename)
- return True, filename
- else:
- if verbose:
- print('feat folder does not exist: ' + filename)
- return False, filename
- # NOTE: the 'beta_map'/'beta_maps' branches below probe the *legacy* FEAT
- # layout (.feat/stats/pe*). They predate step 0.5 and have no callers today
- # (both check_file_status callers pass process='GLM'). If you revive them,
- # route through rsa_utils.resolve_beta_map instead -- these paths report
- # "missing" for any run whose .feat has been deleted after step 0.5.
- # They are not fixed in place because rsa_utils imports this module, so
- # importing it back here would be circular.
- elif process == 'beta_map': # check if the beta map exist
- # build the filename for the beta map
- # make sure that all variables are available, if not indicate which is missing
- list_vars = [model, stim_N]
- list_vars_name = ['model', 'stim_N']
- for var, var_name in zip(list_vars, list_vars_name):
- if var is None:
- print(f'{var_name} is missing')
- return False, f'error in {var_name}'
- # build the filename for the beta map
- filename = (datafolder + os.sep + dataset + os.sep + 'results' + os.sep + 'GLM' + os.sep +
- model + os.sep + f"{specie}-sub-{sub_N:02d}" + os.sep +
- f"ses-{session}_task-{task}_run-{run_N:02d}.feat" + os.sep + "stats" + os.sep + f"pe{stim_N*2 - 1}.nii.gz")
- # check if filename exist
- if os.path.exists(filename):
- if verbose:
- print('beta map exists: ' + filename)
- return True, filename
- else:
- if verbose:
- print('beta map does not exist: ' + filename)
- return False, filename
- elif process == 'beta_maps': # check if the beta map exist
- # build the filename for the beta map
- # make sure that all variables are available, if not indicate which is missing
- list_vars = [model, stim_types]
- list_vars_name = ['model', 'stim_types']
- for var, var_name in zip(list_vars, list_vars_name):
- if var is None:
- print(f'{var_name} is missing')
- return False, f'error in {var_name}'
- # initialize a list to store missing files
- files_missing, files_found = [], []
- for stim_Nx,_ in enumerate(stim_types):
- # starts at 1 and jumps of 2
- stim_N = stim_Nx + 1
- filename = (datafolder + os.sep + dataset + os.sep + 'results' + os.sep + 'GLM' + os.sep +
- model + os.sep + f"{specie}-sub-{sub_N:02d}" + os.sep +
- f"ses-{session}_task-{task}_run-{run_N:02d}.feat" + os.sep + "stats" + os.sep + f"pe{stim_N*2 - 1}.nii.gz")
- # add to files_found if exist
- if os.path.exists(filename):
- files_found.append(filename)
- else:
- files_missing.append(filename)
- if len(files_missing) == 0:
- if verbose:
- print('All beta maps exist')
- # return true and a list of all filenames
- return True, files_found
- else:
- if verbose:
- print('Some beta maps are missing')
- print('Missing files:', files_missing)
- return False, files_missing
- elif process == 'pairwise_similarity_maps':
- # build the filename for the pairwise similarity map
- # make sure that all variables are available, if not indicate which is missing
- list_vars = [model, dis_method, radius, stim_types]
- list_vars_name = ['model', 'dis_method', 'radius', 'stim_types']
- for var, var_name in zip(list_vars, list_vars_name):
- if var is None:
- print(f'{var_name} is missing')
- return False, f'error in {var_name}'
- # initialize a list to store missing files
- files_missing, files_found = [], []
- # check if "P:\userdata\raulh87\data\EmoB\results\RSA\basic-block\D-sub-01\ses-01_task-EmoB_run-01\r-3_correlation_A-1_A-2.nii.gz"
- for stim_1_N, stim_1_name in enumerate(stim_types):
- for stim_2_N, stim_2_name in enumerate(stim_types):
- if stim_2_N > stim_1_N: # only check for upper triangle
- filename = (datafolder + os.sep + dataset + os.sep + 'results' + os.sep + 'RSA' + os.sep +
- model + os.sep + f"{specie}-sub-{sub_N:02d}" + os.sep +
- f"ses-{session}_task-{task}_run-{run_N:02d}" + os.sep +
- f"r-{radius}_{dis_method}_{stim_1_name}_{stim_2_name}.nii.gz")
- # add to files_found if exist
- if os.path.exists(filename):
- if verbose:
- print('Exists: ' + filename + 'adding...')
- files_found.append(filename)
- else:
- if verbose:
- print('Missing: ' + filename + 'adding to missing...')
- files_missing.append(filename)
- if len(files_missing) == 0:
- if verbose:
- print('All pairwise similarity maps exist')
- # return true and a list of all filenames
- return True, files_found
- else:
- if verbose:
- print('Some pairwise similarity maps are missing')
- print('Missing files:', files_missing)
- return False, files_missing
- elif process == 'pairwise_similarity_map':
- # build the filename for the pairwise similarity map
- # make sure that all variables are available, if not indicate which is missing
- list_vars = [model, dis_method, radius, stim_1_name, stim_2_name]
- list_vars_name = ['model', 'dis_method', 'radius', 'stim_1_name', 'stim_2_name']
- for var, var_name in zip(list_vars, list_vars_name):
- if var is None:
- print(f'{var_name} is missing')
- return False, f'error in {var_name}'
- # build the filename for the pairwise similarity map
- filename = (datafolder + os.sep + dataset + os.sep + 'results' + os.sep + 'RSA' + os.sep +
- model + os.sep + f"{specie}-sub-{sub_N:02d}" + os.sep +
- f"ses-{session}_task-{task}_run-{run_N:02d}" + os.sep +
- f"r-{radius}_{dis_method}_{stim_1_name}_{stim_2_name}.nii.gz")
- # check if filename exist
- if os.path.exists(filename):
- if verbose:
- print('pairwise similarity map exists: ' + filename)
- return True, filename
- else:
- if verbose:
- print('pairwise similarity map does not exist: ' + filename)
- return False, filename
- elif process == 'model_similarity_map':
- # build the filename for the model similarity map
- # make sure that all variables are available, if not indicate which is missing
- list_vars = [model, dis_method, rsa_method, radius, rsa_model]
- list_vars_name = ['model', 'dis_method', 'rsa_method', 'radius', 'rsa_model'
- ]
- for var, var_name in zip(list_vars, list_vars_name):
- if var is None:
- print(f'{var_name} is missing')
- return False, f'error in {var_name}'
- filename = (datafolder + os.sep + dataset + os.sep + 'results' + os.sep + 'RSA' + os.sep +
- model + os.sep + rsa_model + os.sep + f"{specie}-sub-{sub_N:02d}" + os.sep + f"ses-{session}_task-{task}_run-{run_N:02d}" + os.sep +
- f"r-{radius}_{dis_method}_{rsa_method}.nii.gz")
- # check if filename exist
- if os.path.exists(filename):
- if verbose:
- print('model similarity map exists: ' + filename)
- return True, filename
- else:
- if verbose:
- print('model similarity map does not exist: ' + filename)
- return False, filename
- # write check for beta map
- elif process == 'mean_model_similarity_map':
- # build the filename for the mean model similarity map
- # make sure that all variables are available, if not indicate which is missing
- list_vars = [model, dis_method, rsa_method, radius, rsa_model]
- list_vars_name = ['model', 'dis_method', 'rsa_method', 'radius', 'rsa_model'
- ]
- for var, var_name in zip(list_vars, list_vars_name):
- if var is None:
- print(f'{var_name} is missing')
- return False, f'error in {var_name}'
- filename = (datafolder + os.sep + dataset + os.sep + 'results' + os.sep + 'RSA' + os.sep +
- model + os.sep + rsa_model + os.sep + 'mean' + os.sep +
- f"r-{radius}_{dis_method}_{rsa_method}_mean.nii.gz")
- # check if filename exist
- if os.path.exists(filename):
- if verbose:
- print('mean model similarity map exists: ' + filename)
- return True, filename
- else:
- if verbose:
- print('mean model similarity map does not exist: ' + filename)
- return False, filename
- elif process == 'model_similarity_maps_rnd':
- # checks if the permutations have been done:
- # true, files_found. If all files exist, false if one or more files are missing
- # false, missing_files. If one or more files are missing, list of missing files
- # make sure that all variables are available, if not indicate which is missing
- # reps - number of by participant repetitions
- list_vars = [model, dis_method, rsa_method, radius, rsa_model, reps]
- list_vars_name = ['model', 'dis_method', 'rsa_method', 'radius', 'rsa_model'
- 'reps']
- # check that all variables are available, if not indicate which is missing
- for var, var_name in zip(list_vars, list_vars_name):
- if var is None:
- print(f'{var_name} is missing')
- return False, f'error in {var_name}'
- # initialize a list to store missing files
- files_missing, files_found, file_num_missing = [], [], []
- # "P:\userdata\raulh87\data\EmoB\results\RSA_rnd\basic\emotion_valence\D-sub-01\ses-01_task-EmoB_run-01\r-3_pearson_kendall_0000.nii.gz"
- for rep in range(reps):
- filename = (datafolder + os.sep + dataset + os.sep + 'results' + os.sep + 'RSA_rnd' + os.sep +
- model + os.sep + rsa_model + os.sep + f"{specie}-sub-{sub_N:02d}" + os.sep + f"ses-{session}_task-{task}_run-{run_N:02d}" + os.sep +
- f"r-{radius}_{dis_method}_{rsa_method}_{str(rep).zfill(4)}.nii.gz")
- # add to files_found if exist
- if os.path.exists(filename):
- files_found.append(filename)
- else:
- files_missing.append(filename)
- file_num_missing.append(rep)
- if len(files_missing) == 0:
- if verbose:
- print('All model similarity rnd maps exist')
- # return true and a list of all filenames
- return True, files_found
- else:
- if verbose:
- print('Some model similarity rnd maps are missing')
- print('Missing files:', files_missing)
- return False, file_num_missing
- elif process == 'mean_model_similarity_maps_rnd':
- # checks if the permutations have been done:
- # true, files_found. If all files exist, false if one or more files are missing
- # false, missing_files. If one or more files are missing, list of missing files
- # make sure that all variables are available, if not indicate which is missing
- # reps - number of by participant repetitions
- list_vars = [model, dis_method, rsa_method, radius, rsa_model, reps, reps_group]
- list_vars_name = ['model', 'dis_method', 'rsa_method', 'radius', 'rsa_model',
- 'reps_group']
- # check that all variables are available, if not indicate which is missing
- for var, var_name in zip(list_vars, list_vars_name):
- if var is None:
- print(f'{var_name} is missing')
- return False, f'error in {var_name}'
- # initialize a list to store missing files
- files_missing, files_found, file_num_missing = [], [], []
- # "P:\userdata\raulh87\data\EmoB\results\RSA_rnd\basic\emotion_valence\r-r-3_pearson_kendall_mean_07181.nii.gz"
- for g_rep in range(reps_group):
- filename = (datafolder + os.sep + dataset + os.sep +
- 'results' + os.sep + 'RSA_rnd' + os.sep +
- model + os.sep + rsa_model + os.sep +
- 'r-' + str(radius) + '_' + dis_method + '_' + rsa_method + '_mean_' + str(g_rep).zfill(len(str(reps_group))) + '.nii.gz')
- # add to files_found if exist
- if os.path.exists(filename):
- files_found.append(filename)
- else:
- files_missing.append(filename)
- file_num_missing.append(rep)
- if len(files_missing) == 0:
- if verbose:
- print('All model similarity rnd maps exist')
- # return true and a list of all filenames
- return True, files_found
- else:
- if verbose:
- print('Some model similarity rnd maps are missing')
- print('Missing files:', files_missing)
- return False, file_num_missing
- elif process == 'get_vectors':
- # "P:\userdata\raulh87\data\EmoB\BIDS\sub-01\sub-01_ses-01_task-EmoB_run-01_events.csv"
- # if session is int convert to str with 2 digits
- if isinstance(session, int):
- session = str(session).zfill(2)
- # check if the events file exists
- filename = (datafolder + os.sep + dataset + os.sep + 'BIDS' + os.sep +
- specie + '-sub-' + str(sub_N).zfill(2) + os.sep +
- specie + '-sub-' + str(sub_N).zfill(2) +
- '_ses-' + session +
- '_task-' + task +
- '_run-' + str(run_N).zfill(2) + '_events.csv')
- # check if the file exists
- if os.path.exists(filename):
- if verbose:
- print('File exists: ' + filename)
- return True, filename
- else:
- if verbose:
- print('File does not exist: ' + filename)
- return False, filename
- else:
- # print process (process) not found
- print('Process not found: ' + process)
- return False, 'error in process'
- def preprocess_run(sub_N, run_N, dataset, task, specie, datafolder, session, smooth=0, combination=['-x','z','-y'], run_prepro=True):
- """
- Preprocesses a single run of a single subject.
- Reorients file
- """
- ## determine input file and output directory ##
- # input directory in BIDS format
- input_folder = datafolder + os.sep + dataset + os.sep + 'BIDS' + os.sep + specie + '-sub-' + str(sub_N).zfill(2)
- filename = input_folder + os.sep + specie + '-sub-' + str(sub_N).zfill(2) + '_ses-' + session + '_task-' + task + '_run-' + str(run_N).zfill(2) + '_bold.nii.gz'
- filename_json = input_folder + os.sep + specie + '-sub-' + str(sub_N).zfill(2) + '_ses-' + session + '_task-' + task + '_run-' + str(run_N).zfill(2) + '_bold.json'
- # print filename
- print('Input file: ' + filename)
- # get TR and number of volumes
- TR,volumes = utils.extract_params(filename)
- # create output directory, where the fsl output will be saved (preprocessed data)
- outputdir = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2)
- fsl_outputdir = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2) + os.sep + specie + '-sub-' + str(sub_N).zfill(2) + '_ses-' + session + '_task-' + task + '_run-' + str(run_N).zfill(2)
- slice_timming_path = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2) + os.sep + 'slice_timming_' + specie + '-sub-' + str(sub_N).zfill(2) + '_ses-' + session + '_task-' + task + '_run-' + str(run_N).zfill(2) + '.txt'
- # check that slice_timming_path folder exist
- if not os.path.exists(os.path.dirname(slice_timming_path)):
- os.makedirs(os.path.dirname(slice_timming_path))
- # get slice timing parameters
- slice_timming = utils.get_slice_timing(filename_json, slice_timming_path)
- # make sure the file was created
- if not os.path.exists(slice_timming_path):
- print('Slice timing file was not created: ' + slice_timming_path)
- return
- # check if slice_timming is not empty
- if slice_timming is None:
- print('No slice timing parameters found')
- print('File: ' + slice_timming_path + ' empty')
- return
- ## Filling out the design.fsf file ##
- # create list of labels to fill in the design.fsf file
- label_list = ['Outputdir', 'TR', 'Volumes', 'BET', 'Smooth', 'Input', 'SliceTimming']
- # create dictionary to fill in the design.fsf file
- to_fill_dict = dict()
- for label in label_list:
- to_fill_dict[label] = dict()
- if label == 'Outputdir':
- to_fill_dict[label]['string_to_find'] = 'set fmri(outputdir)'
- to_fill_dict[label]['string_to_replace'] = ('set fmri(outputdir) "' + fsl_outputdir + '"')
- elif label == 'TR':
- to_fill_dict[label]['string_to_find'] = 'set fmri(tr)'
- to_fill_dict[label]['string_to_replace'] = ('set fmri(tr) ' + str(TR))
- elif label == 'Volumes':
- to_fill_dict[label]['string_to_find'] = 'set fmri(npts)'
- to_fill_dict[label]['string_to_replace'] = ('set fmri(npts) ' + str(volumes))
- elif label == 'BET':
- to_fill_dict[label]['string_to_find'] = 'set fmri(bet_yn)'
- if specie == 'H':
- to_fill_dict[label]['string_to_replace'] = ('set fmri(bet_yn) 1')
- elif specie == 'D':
- to_fill_dict[label]['string_to_replace'] = ('set fmri(bet_yn) 0')
- elif label == 'Smooth':
- to_fill_dict[label]['string_to_find'] = 'set fmri(smooth)'
- to_fill_dict[label]['string_to_replace'] = ('set fmri(smooth) ' + str(smooth))
- elif label == 'Input':
- to_fill_dict[label]['string_to_find'] = 'set feat_files(1)'
- to_fill_dict[label]['string_to_replace'] = ('set feat_files(1) "' + filename + '"')
- elif label == 'SliceTimming':
- to_fill_dict[label]['string_to_find'] = 'set fmri(st_file)'
- to_fill_dict[label]['string_to_replace'] = ('set fmri(st_file) "' + slice_timming_path + '"')
- # fill in the design.fsf file
- design_path = os.path.join(os.getcwd(), 'FSL_designs' + os.sep + 'preprocess_slice_timming_from_JSON.fsf')
- # design_path = os.path.join(os.getcwd(), 'FSL_designs' + os.sep + 'preprocess_no-slice-timing.fsf')
- print('Design path: ' + design_path)
- design_modified_path = os.path.join(os.getcwd(), 'FSL_designs' + os.sep + 'preprocess_modified.fsf')
- if run_prepro:
- utils.fill_fsf(to_fill_dict, design_path, design_modified_path)
- # check if previous feat preprocessing directory exists, if so, delete it
- if os.path.exists(fsl_outputdir + '.feat'):
- shutil.rmtree(fsl_outputdir + '.feat')
- # run feat
- command = 'feat ' + design_modified_path
- # check if system is windows, if so, do not execute command
- print(command)
- if os.name == 'nt':
- print("The system is windows, command not executed.")
- elif os.name == 'posix':
- os.system(command)
- else:
- os.system(command)
- else:
- print('Preprocessing not run, skipping this step')
- ## reorient run ##
- base_filename =(
- specie + '-sub-' + str(sub_N).zfill(2) +
- '_ses-' + session +
- '_task-' + task +
- '_run-' + str(run_N).zfill(2))
- # non-oriented file
- non_oriented_file = base_filename + '_not-oriented.nii.gz'
- # oriented file
- reoriented_file = base_filename + '_reoriented.nii.gz'
- preprocessed_file = fsl_outputdir + '.feat' + os.sep + 'filtered_func_data.nii.gz'
- # check if the system is windows
- if os.name == 'nt':
- print('A copy should have been created, but this is Windows')
- print(non_oriented_file + ' a copy of this file here:')
- print('FSL output directory: ' + fsl_outputdir)
- else:
- #copy preprocessed_file to non_oriented_file
- shutil.copyfile(preprocessed_file, outputdir + os.sep + non_oriented_file)
- print(non_oriented_file + ' created')
- print('FSL output directory: ' + fsl_outputdir)
- # check if system is windows, if so, do not execute command
- if os.name == 'nt': # Windows
- print("system is windows, command not executed.")
- print("orientation to use: ",combination)
- print("non-oriented file: " + outputdir + os.sep + non_oriented_file)
- print("oriented file: " + outputdir + os.sep + reoriented_file)
- utils.reorient_file(outputdir + os.sep + non_oriented_file, outputdir + os.sep + reoriented_file, combination)
- def get_mean_fct(sub_N, session_and_run, base_run, dataset, task, specie, datafolder, first_time=True):
- """
- Calculates the mean functional image for a subject and a task.
- The mean functional image is calculated by averaging the mean images of each run.
- The mean image of each run is calculated by averaging all volumes of the run.
- The mean image of each run is calculated by averaging all volumes of the run.
- The motion is corrected for each run using the first volume of the first run as reference.
- The motion parameters are saved in a .par file.
- The mean image of each run
- """
- # Check if the system is windows
- if os.name == 'nt':
- print('The system is Windows, this is a test, no actual system or FSL commands will be run')
- # output directory where the fsl output will be saved (preprocessed data)
- outputdir = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2)
- # movement directory
- movementdir = datafolder + os.sep + dataset + os.sep + 'movement'
- # create movement directory if it does not exist
- if not os.path.exists(movementdir):
- os.makedirs(movementdir)
- # get initial run from session_and_run[0] = ['ses-01_run-1', 'ses-02_run-2']
- initial_run = int(session_and_run[0].split('_')[1].split('-')[1])
- initial_session = session_and_run[0].split('_')[0].split('-')[1]
- ## obtain volume to be used as base to correct all others ##
- filename = (specie + '-sub-' + str(sub_N).zfill(2) +
- '_ses-' + initial_session +
- '_task-' + task +
- '_run-' + str(initial_run).zfill(2) +
- '_reoriented.nii.gz')
- if first_time: # get the volume
- # get the first volume of the first run to use as base volume
- command = f"fslroi {outputdir + os.sep + filename} {outputdir + os.sep + 'base_vol.nii.gz'} 0 1"
- # print commmand
- print(command)
- # if the system is not windows, run the command
- if os.name != 'nt':
- os.system(command)
- else:
- # will use average created before
- print('Using existing base volume')
- # check if base_vol exists
- if not os.path.exists(outputdir + os.sep + 'base_vol.nii.gz'):
- print('base_vol.nii.gz does not exist, run the code with first_time = True')
- print('path: ' + outputdir + os.sep + 'base_vol.nii.gz')
- raise ValueError('base_vol.nii.gz does not exist, run the code with first_time = True')
- ## ----- ##
- # This string will be used to generate the mean image
- mean_images = ''
- ## calculate motion for each run and generate par file ##
- for n, file_ending in enumerate(session_and_run):
- print('processing ' + file_ending)
- session = file_ending.split('_')[0].split('-')[1]
- run_N = int(file_ending.split('_')[1].split('-')[1])
- print('processing ' + str(n+1) + ' of ' + str(len(session_and_run)) + ' runs')
- filename = (specie + '-sub-' + str(sub_N).zfill(2) +
- '_ses-' + session +
- '_task-' + task +
- '_run-' + str(run_N).zfill(2))
- # '_reoriented.nii.gz')
- # average the reoriented file to create a mean image
- print('calculating mean image...')
- command = f"fslmaths {outputdir + os.sep + filename + '_reoriented.nii.gz'} -Tmean {outputdir + os.sep + filename + '_mean_unaligned.nii.gz'}"
- # print commmand
- print(command)
- # if the system is not windows, run the command
- if os.name != 'nt':
- os.system(command)
- # command to calculate transformation matrix from mean image to base volume
- command = f"flirt -in {outputdir + os.sep + filename + '_mean_unaligned.nii.gz'} -ref {outputdir + os.sep + 'base_vol.nii.gz'} -out {outputdir + os.sep + filename + '_mean_aligned_tmp.nii.gz'} -omat {outputdir + os.sep + filename + '2base_vol.mat'}"
- # print commmand
- print(command)
- if os.name != 'nt':
- os.system(command)
- # apply the transformation matrix to the 4D file
- command = f"applyxfm4D {outputdir + os.sep + filename + '_reoriented.nii.gz'} {outputdir + os.sep + 'base_vol.nii.gz'} {outputdir + os.sep + filename + '_mc.nii.gz'} {outputdir + os.sep + filename + '2base_vol.mat'} -singlematrix"
- # print commmand
- print(command)
- # if the system is not windows, run the command
- if os.name != 'nt':
- os.system(command)
- # apply the transformation matrix to the file
- print('calculating mean image...')
- command = f"fslmaths {outputdir + os.sep + filename + '_mc.nii.gz'} -Tmean {outputdir + os.sep + filename + '_mean.nii.gz'}"
- # print commmand
- print(command)
- # if the system is not windows, run the command
- if os.name != 'nt':
- os.system(command)
- # if the system is not windows, run the command
- # print file saved
- print('aligned file saved as ' + filename + '_mc.nii.gz')
- # add filename to mean_images
- mean_images += outputdir + os.sep + filename + '_mean.nii.gz' + ' '
- if first_time: # if yes, calculate mean image
- mean_fct_file = (outputdir + os.sep +
- specie + '-sub-' + str(sub_N).zfill(2) +
- '_task-' + task + '_mean_fct_uncut.nii.gz'
- )
- # append mean images to a single 4D image
- command = f"fslmerge -t {mean_fct_file} {mean_images}"
- # print commmand
- print(command)
- # if the system is not windows, run the command
- if os.name != 'nt':
- os.system(command)
- print('mean fct file saved as ' + mean_fct_file)
- # calculate mean image
- command = f"fslmaths {mean_fct_file} -Tmean {mean_fct_file}"
- # print commmand
- print(command)
- # if the system is not windows, run the command
- if os.name != 'nt':
- os.system(command)
- print('done')
- def mean_to_STD(sub_N, dataset, task, specie, datafolder, atlas_type, img_type='brain2mm', variation='STD0', use_anatomic=False):
- """
- This function will take the mean functional image of a subject and transform it to the space of the atlas.
- sub_N: subject number
- dataset: dataset name
- task: task name
- specie: specie name
- datafolder: path to the data folder
- atlas_type: type of atlas to use
- img_type: type of image to use, default is 'brain2mm'
- """
- # working directory
- workingdir = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2)
- mean_fct_file = workingdir + os.sep + specie + '-sub-' + str(sub_N).zfill(2)
- cut_mean_fct_file = mean_fct_file + '_task-' + task + '_mean_fct.nii.gz'
- mean_fct_file += '_task-' + task + '_mean_fct.nii.gz'
- masked_mean_fct_file = mean_fct_file[:-7] + '_brain.nii.gz'
- mean_fct_file_STD = mean_fct_file[:-7] + '_STD.nii.gz'
- mean_fct2STD_mat = mean_fct_file[:-7] + '2STD.mat'
- #T1w_brain_file = datafolder + os.sep + dataset + os.sep + 'BIDS' + os.sep + 'sub-' + str(sub_N).zfill(2) + os.sep + 'sub-' + str(sub_N).zfill(2) + '_T1w_brain.nii.gz'
- T1w_brain_file = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2) + os.sep + specie + '-sub-' + str(sub_N).zfill(2) + '_T1w_brain.nii.gz'
- T1w_brain_file_unoriented = datafolder + os.sep + dataset + os.sep + 'BIDS' + os.sep + 'sub-' + str(sub_N).zfill(2) + os.sep + 'sub-' + str(sub_N).zfill(2) + '_T1w_brain.nii.gz'
- T1w_brain_STD_file = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2) + os.sep + specie + '-sub-' + str(sub_N).zfill(2) + '_T1w_brain_STD.nii.gz'
- T1w_brain2STD_mat = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2) + os.sep + specie + '-sub-' + str(sub_N).zfill(2)+ '_T1w_brain2STD.mat'
- mean_fct_T1w_file = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2) + os.sep + specie + '-sub-' + str(sub_N).zfill(2)+ '_mean_fct_T1w.nii.gz'
- mean_fct2T1w_mat = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2) + os.sep + specie + '-sub-' + str(sub_N).zfill(2)+ '_mean_fct2T1w.mat'
- if specie == 'D':
- specieS = 'Dog'
- elif specie == 'H':
- specieS = 'Hum'
- # generate path to atlas. The atlas is in the same folder as the script
- atlas_file = os.getcwd() + os.sep + "Atlas" + os.sep + specieS + os.sep + atlas_type + os.sep + img_type + ".nii.gz"
- # check which variation to use
- # A - mean
- # B - anatomic
- # C - STD
- if use_anatomic:
- # mean_to_anatomic(sub_N, dataset, task, specie, datafolder) # A -> B
- # anatomic_to_STD(sub_N, dataset, task, specie, datafolder, atlas_type, img_type=img_type) # B -> C
- # apply fslreorient2std to the T1w_brain file
- print('Reorienting T1w_brain file to standard orientation')
- command = f"fslreorient2std {T1w_brain_file_unoriented} {T1w_brain_file}"
- print(command)
- #if system is not windows, run the command
- if os.name != 'nt':
- os.system(command)
- print('Transforming masked mean functional image to T1w_brain A -> B')
- # masked_mean_fct to T1w_brain A -> B
- command = f"flirt -in {masked_mean_fct_file} -ref {T1w_brain_file} -out {mean_fct_T1w_file} -omat {mean_fct2T1w_mat} -bins 256 -cost corratio -searchrx -30 30 -searchry -30 30 -searchrz -30 30 -dof 7 -interp trilinear"
- # if the system is windows, don't run the command, just write it down
- print(command)
- if os.name == 'nt': # Windows
- print("System is Windows, command not executed")
- else:
- os.system(command)
- # T1w_brain to atlas B -> C T1w_brain_STD_file
- print('Transforming T1w_brain to atlas space B -> C')
- command = f"flirt -in {T1w_brain_file} -ref {atlas_file} -out {T1w_brain_STD_file} -omat {T1w_brain2STD_mat} -bins 256 -cost corratio -searchrx -45 45 -searchry -45 45 -searchrz -45 45 -dof 12 -interp trilinear"
- # flirt -in ${dataFolder}/data/${sub}/masks/${brain_nMeanfct}.nii.gz -ref ${atlasFile} -out ${dataFolder}/data/${sub}/masks/brainSTD.nii.gz -omat ${dataFolder}/data/${sub}/masks/brain2STD.mat -bins 256 -cost corratio -searchrx -45 45 -searchry -45 45 -searchrz -45 45 -dof 12 -interp trilinear
- print(command)
- if os.name != 'nt': # Windows
- os.system(command)
- print('Calculating transformation matrix from mean functional image to atlas A->B, B->C = A->C')
- # adding the matrices... #A->B + B->C = A->C
- command = f"convert_xfm -omat {mean_fct2STD_mat} -concat {T1w_brain2STD_mat} {mean_fct2T1w_mat}"
- print(command)
- if os.name == 'nt': # Windows
- print("System is Windows, command not executed")
- else:
- os.system(command)
- # apply the transformation matrix to the mean functional image
- print('Applying transformation matrix to masked mean functional image')
- command = f"flirt -in {masked_mean_fct_file} -ref {atlas_file} -out {mean_fct_file_STD} -applyxfm -init {mean_fct2STD_mat} -interp trilinear"
- else: # use mean functional image directly
- print('Transforming masked mean functional image directly to atlas space')
- print('variation: ' + variation)
- if variation == 'STD0':
- command = f"flirt -in {masked_mean_fct_file} -ref {atlas_file} -out {mean_fct_file_STD} -omat {mean_fct2STD_mat} -bins 256 -cost corratio -searchrx -90 90 -searchry -90 90 -searchrz -90 90 -dof 12 -interp trilinear"
- elif variation == 'STD1':
- command = f"flirt -in {masked_mean_fct_file} -ref {atlas_file} -out {mean_fct_file_STD} -omat {mean_fct2STD_mat} -bins 256 -cost corratio -searchrx -30 30 -searchry -30 30 -searchrz -30 30 -dof 12 -interp trilinear"
- else: #error, variation not found
- print('Variation not found')
- #generate error message
- return
- # if the system is windows, don't run the command, just write it down
- print(command)
- if os.name == 'nt': # Windows
- print("System is Windows, command not executed")
- else:
- os.system(command)
- def mean_to_anatomic(sub_N, dataset, task, specie, datafolder):
- '''
- input:
- \BIDS\sub-01_T1w_brain.nii.gz: anatomic file in BIDS format, skull removed, reoriented and labeled
- output:
- '''
- workingdir = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2)
- T1w_brain_file = datafolder + os.sep + dataset + os.sep + 'BIDS' + os.sep + 'sub-' + str(sub_N).zfill(2) + '_T1w_brain.nii.gz'
- mean_fct_T1w_file = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2) + '_mean_fct_T1w.nii.gz'
- mean_fct2T1w_mat = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2) + '_mean_fct2T1w.mat'
- mean_fct_file = workingdir + os.sep + specie + '-sub-' + str(sub_N).zfill(2) # like: D-sub-01_task-EmoB_run-01_mean_fct.nii.gz
- mean_fct_file += '_task-' + task + '_mean_fct.nii.gz'
- masked_mean_fct_file = mean_fct_file[:-7] + '_brain.nii.gz'
- # mean_fct_file_STD = mean_fct_file[:-7] + '_STD.nii.gz'
- # mean_fct2STD_mat = mean_fct_file[:-7] + '2STD.mat'
- if specie == 'D':
- specieS = 'Dog'
- elif specie == 'H':
- specieS = 'Hum'
- # masked_mean_fct to T1w_brain
- command = f"flirt -in {masked_mean_fct_file} -ref {T1w_brain_file} -out {mean_fct_T1w_file} -omat {mean_fct2T1w_mat} -bins 256 -cost corratio -searchrx -30 30 -searchry -30 30 -searchrz -30 30 -dof 7 -interp trilinear"
- # if the system is windows, don't run the command, just write it down
- print(command)
- if os.name == 'nt': # Windows
- print("System is Windows, command not executed")
- else:
- os.system(command)
- def anatomic_to_STD(sub_N, dataset, task, specie, datafolder, atlas_type, img_type='brain2mm'):
- '''
- input:
- \BIDS\sub-01_T1w_brain.nii.gz: anatomic file in BIDS format, skull removed, reoriented and labeled
- output:
- '''
- workingdir = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2)
- T1w_brain_file = datafolder + os.sep + dataset + os.sep + 'BIDS' + os.sep + 'sub-' + str(sub_N).zfill(2) + '_T1w_brain.nii.gz'
- T1w_brain_STD_file = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2) + '_T1w_brain_STD.nii.gz'
- T1w_brain2STD_mat = datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep + specie + '-sub-' + str(sub_N).zfill(2) + '_T1w_brain2STD.mat'
- if specie == 'D':
- specieS = 'Dog'
- elif specie == 'H':
- specieS = 'Hum'
- # generate path to atlas. The atlas is in the same folder as the script
- atlas_file = os.getcwd() + os.sep + "Atlas" + os.sep + specieS + os.sep + atlas_type + os.sep + img_type + ".nii.gz"
- # T1w_brain to atlas
- command = f"flirt -in {T1w_brain_file} -ref {atlas_file} -out {T1w_brain_STD_file} -omat {T1w_brain2STD_mat} -bins 256 -cost corratio -searchrx -30 30 -searchry -30 30 -searchrz -30 30 -dof 7 -interp trilinear"
- # if the system is windows, don't run the command, just write it down
- print(command)
- if os.name == 'nt': # Windows
- print("System is Windows, command not executed")
- else:
- os.system(command)
- def run_to_STD(sub_N, run_N, dataset, task, specie, datafolder, atlas_type, img_type='brain2mm', session=''):
- """
- This function will take a semi-processed run of a participant
- cut it, apply BET and transform it to the space of the atlas.
- sub_N: subject number
- run_N: run number
- dataset: dataset name
- task: task name
- specie: specie name
- datafolder: path to the data folder
- atlas_type: type of atlas to use
- img_type: type of image to use, default is 'brain2mm'
- session: session number (in case there is one)
- """
- # if the system is windows
- if os.name == 'nt':
- print('The system is Windows, FSL functions ans bash scripts will not be executed')
- # working directories
- std_dir = datafolder + os.sep + dataset + os.sep + 'normalized' + os.sep + specie + '-sub-' + str(sub_N).zfill(2)
- preprocess_dir = (datafolder + os.sep + dataset + os.sep + 'preprocessing' + os.sep +
- specie + '-sub-' + str(sub_N).zfill(2)
- )
- # cutting parameters file
- params_file = (preprocess_dir + os.sep +
- specie + '-sub-' + str(sub_N).zfill(2) +
- '_task-' + task + '_mean_fct_uncut_cut_params.txt'
- )
- filename = (specie + '-sub-' + str(sub_N).zfill(2) +
- '_ses-' + session +
- '_task-' + task +
- '_run-' + str(run_N).zfill(2))
- # name for reoriented and motion corrected file
- preprocessed_file = filename + '_mc.nii.gz'
- cut_file = filename + '_reoriented_mc_cut.nii.gz'
- # updating output params_file
- params_dict = utils.read_params_file(params_file)
- params_dict['output_file'] = preprocess_dir + os.sep + cut_file
- params_file_current = preprocess_dir + os.sep + filename + '_cut_params.txt'
- # params_file_current = params_file[:-30] + '_run-' + str(run_N).zfill(2) + '_cut_params.txt'
- # add mask file to parameters
- params_dict['mask_file'] = (preprocess_dir + os.sep +
- specie + '-sub-' + str(sub_N).zfill(2) +
- '_task-' + task + '_mean_fct_mask.nii.gz')
- # save the cutting parameters
- utils.write_params_file(params_file_current, params_dict)
- # cut preprocessed file and apply BET
- command = f"./remove_slices.sh {preprocess_dir + os.sep + preprocessed_file} {params_file_current}"
- print(command)
- if os.name != 'nt': # Windows
- os.system(command)
- # apply transformation to STD
- mean_fct2STD_mat = (
- preprocess_dir + os.sep +
- specie + '-sub-' + str(sub_N).zfill(2) +
- '_task-' + task + '_mean_fct2STD.mat'
- )
- # determine folder for the atlas
- if specie == 'D':
- specieS = 'Dog'
- elif specie == 'H':
- specieS = 'Hum'
- atlas_file = os.getcwd() + os.sep + "Atlas" + os.sep + specieS + os.sep + atlas_type + os.sep + img_type + ".nii.gz"
- # if normalized directory does not exist, create it
- if not os.path.exists(std_dir):
- os.makedirs(std_dir)
- print('Directory ' + std_dir + ' created')
- else:
- print('Directory ' + std_dir + ' already exists')
- command = f"flirt -in {preprocess_dir + os.sep + cut_file} -ref {atlas_file} -applyxfm -init {mean_fct2STD_mat} -out {std_dir + os.sep + filename + '.nii.gz'} -interp trilinear"
- print(command)
- if os.name != 'nt': # Windows
- os.system(command)
- print('done')
- from pathlib import Path
- import numpy as np
- from scipy.signal import detrend as _sp_detrend
- def fwd(par_file, radius=50.0, threshold=0.5, detrend_type="linear-demean", output_file=None, add_movement_params=True):
- """
- Compute framewise displacement (FD) and return a binary mask of frames above threshold.
- If written to disk, the file includes the 6 motion parameters followed by one
- additional censoring column per volume that exceeded the FD threshold.
- Parameters
- ----------
- par_file : str or Path
- Path to FSL .par file with 6 motion parameters (rotations in radians first).
- radius : float, optional
- Radius in mm used to convert rotations to displacement (default is 50).
- threshold : float, optional
- Threshold in mm for FD (default is 0.5).
- detrend_type : str, optional
- Detrending method: 'linear-demean', 'linear-nodemean', or 'none' (default is 'linear-demean').
- output_file : str or Path, optional
- If provided, saves a .txt file with the 6 motion parameters followed by one
- censoring column for each excluded volume.
- add_movement_params : bool, optional
- If True, includes motion parameters in the output file (default is True).
- Returns
- -------
- numpy.ndarray
- Array of 1s and 0s where 1 = frame above or equal to threshold, 0 = below.
- """
- par_file = Path(par_file)
- motion = np.loadtxt(par_file, ndmin=2, dtype=float)
- if motion.shape[1] != 6:
- raise ValueError(f"Expected 6 motion columns, got {motion.shape[1]}")
- motion_detrended = _detrend_columns(motion, detrend_type)
- motion_detrended[:, :3] *= radius # convert radians to mm
- d_motion = np.vstack([np.zeros((1, 6)), np.diff(motion_detrended, axis=0)])
- fd = np.sum(np.abs(d_motion), axis=1)
- mask = (fd >= threshold).astype(np.int8)
- if output_file:
- if add_movement_params:
- excluded_idx = np.flatnonzero(mask)
- censor_columns = np.zeros((motion.shape[0], excluded_idx.size), dtype=float)
- if excluded_idx.size > 0:
- censor_columns[excluded_idx, np.arange(excluded_idx.size)] = 1.0
- output_data = np.hstack([motion, censor_columns])
- np.savetxt(output_file, output_data, fmt="%.6f", delimiter="\t")
- print(f"Motion parameters and {excluded_idx.size} censoring columns saved to {output_file}")
- else:
- if mask.size == 0:
- np.savetxt(output_file, np.empty((0, 0)), fmt="%d", delimiter="\t")
- else:
- np.savetxt(output_file, mask, fmt="%d", delimiter="\t")
- print(f"FD mask saved to {output_file}")
- return mask
- def _detrend_columns(arr, kind="linear-demean"):
- if kind == "none":
- return arr.copy()
- out = arr.copy()
- demean = kind == "linear-demean"
- if kind in {"linear-demean", "linear-nodemean"}:
- for i in range(arr.shape[1]):
- col = arr[:, i]
- mu = col.mean() if demean else 0.0
- out[:, i] = _sp_detrend(col) + mu
- else:
- raise NotImplementedError("Only 'linear-demean', 'linear-nodemean', and 'none' are supported.")
- return out
preprocess_functions.py at commit 55ed729, no license · at the source
Overview
- Social, Cognitive and Affective Neuroscience Unit, Department of Cognition, Emotion, and Methods in Psychology, Faculty of Psychology, University of Vienna, Liebiggasse 5, 1010 Vienna, Austria
- Neuroethology of Communication Lab, Department of Ethology, Institute of Biology, Eötvös Loránd University, Pázmány Péter Sétány 1/C, 1117 Budapest, Hungary
- ELTE NAP Canine Brain Research Group, Pázmány Péter Sétány 1/C, 1117 Budapest, Hungary
- Department of Behavioral and Cognitive Neurobiology, Institute of Neurobiology, National Autonomous University of Mexico, Campus Juriquilla, Boulevard Juriquilla 3001, Querétaro 76230, México
- Faculty of Psychology, National Autonomous University of Mexico, Circuito Ciudad Universitaria, Mexico City 04510, México
Abstract
Dogs can distinguish human facial expressions, particularly happiness, yet brain processes remain unclear. Using fMRI, we conducted two experiments in awake pet dogs. In Experiment 1 (n = 8), happy faces elicited a stronger response than neutral faces in a right temporal cluster extending to the caudate nucleus, including the rostral Sylvian gyrus. In Experiment 2 (n = 12), dogs viewed faces expressing happiness, anger, fear, or sadness. Using the Experiment 1 cluster as a region of interest, a machine-learning classifier distinguished happiness from each negative facial expression, but not between negative pairs, showing differential BOLD responsiveness to happy faces. Whole-brain representational similarity analyses revealed activity patterns differentiating angry vs. fearful faces (right mid ectosylvian and left splenial gyri) and sad vs. fearful faces (right rostral suprasylvian gyrus) but not angry vs. sad faces. This provides direct evidence that dog brains can distinguish between two negative facial expressions.
Reproduced under the paper's license (CC BY), from the paper cited above.
Repository
Its files are read in the Code ↔ Paper reader above, with 5 matches between paragraphs and lines of code.
rhernandez00/dog_brain_toolkit
55ed729dab520cc0786c3c97e8a26d3f232abf5f, 22 September 2026Availability: 1 check, the latest on 27 September 2026: the link answers
- 27 September 2026: the link answers
137 files
- Atlas/
transform_mask.sh , Shell, 23 lines - Quality_check/
NIfTI_toolbox/ , MATLAB, 554 linesaffine.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 94 linesbipolar.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 189 linesbresenham_line3d.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 115 linesclip_nii.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 260 linescollapse_nii_scan.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 48 linesexpand_nii_scan.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 255 linesextra_nii_hdr.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 84 linesflip_lr.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 164 linesget_nii_frame.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 199 linesload_nii.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 207 linesload_nii_ext.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 280 linesload_nii_hdr.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 392 linesload_nii_img.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 200 linesload_untouch0_nii_hdr.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 187 linesload_untouch_header_only .m - Quality_check/
NIfTI_toolbox/ , MATLAB, 191 linesload_untouch_nii.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 217 linesload_untouch_nii_hdr.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 468 linesload_untouch_nii_img.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 210 linesmake_ana.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 256 linesmake_nii.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 83 linesmat_into_hdr.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 142 linespad_nii.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 321 linesreslice_nii.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 179 linesrri_file_menu.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 106 linesrri_orient.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 251 linesrri_orient_ui.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 636 linesrri_select_file.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 92 linesrri_xhair.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 33 linesrri_zoom_menu.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 286 linessave_nii.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 38 linessave_nii_ext.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 227 linessave_nii_hdr.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 219 linessave_untouch0_nii_hdr.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 71 linessave_untouch_header_only .m - Quality_check/
NIfTI_toolbox/ , MATLAB, 232 linessave_untouch_nii.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 207 linessave_untouch_nii_hdr.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 580 linessave_untouch_slice.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 40 linesunxform_nii.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 45 linesverify_nii_ext.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 4,873 linesview_nii.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 480 linesview_nii_menu.m - Quality_check/
NIfTI_toolbox/ , MATLAB, 521 linesxform_nii.m - Quality_check/
functions/ , MATLAB, 25 linesbramila_bold2perc.m - Quality_check/
functions/ , MATLAB, 109 linesbramila_dvars.m - Quality_check/
functions/ , MATLAB, 280 linesbramila_framewiseDisplac ement.m - Quality_check/
functions/ , MATLAB, 22 linescreateTxtColumns.m - Quality_check/
functions/ , MATLAB, 82 linesfuctionToTxt.m - Quality_check/
functions/ , MATLAB, 27 linesload_niiFolder.m - Quality_check/
functions/ , MATLAB, 280 linesload_nii_hdr.m - Quality_check/
functions/ , MATLAB, 191 linesload_untouch_nii.m - Quality_check/
functions/ , MATLAB, 217 linesload_untouch_nii_hdr.m - Quality_check/
functions/ , MATLAB, 35 lineswriteTxt.m - Quality_check/
run_quality_check.m , MATLAB, 117 lines - affine_correction.ipynb, Jupyter, 55 lines
- convert_coords.ipynb, Jupyter, 943 lines, 1 match
- docs/
app.js , JavaScript, 178 lines - export_static.py, Python, 620 lines, 1 match
- interface.ipynb, Jupyter, 897 lines
- job_status.py, Python, 91 lines
- old/
crop_and_bet.ipynb , Jupyter, 394 lines - old/
preprocess.ipynb , Jupyter, 1,529 lines - old/
running.ipynb , Jupyter, 73 lines - preprocess_CAPS-Knee.ipy
nb , Jupyter, 550 lines - preprocess_functions.py, Python, 1,594 lines, 1 match
- remove_slices.sh, Shell, 103 lines
- resample_atlas.sh, Shell, 26 lines
- rsa_model_builder.py, Python, 2,318 lines
- rsa_utils.py, Python, 4,180 lines
- run_GLM.py, Python, 239 lines
- run_bet.sh, Shell, 69 lines
- run_jobs.py, Python, 206 lines
- schedule_rsa.py, Python, 103 lines
- scheduler/
__init__.py , Python, 1 line - scheduler/
dag.py , Python, 271 lines - scheduler/
jobs.py , Python, 207 lines - scheduler/
paths.py , Python, 31 lines - searchlight.py, Python, 1,094 lines
- searchlight_difference.p
y , Python, 257 lines - step11.py, Python, 51 lines
- sync_rsa_to_drive.py, Python, 225 lines
- tests/
test_colab_regression.py , Python, 216 lines - tests/
test_colab_regression_in , Python, 124 linesference.py - tests/
test_generate_fsf.py , Python, 68 lines - tests/
test_regression_pipeline , Python, 240 lines.py - tools/
build_all_categories_gro , Python, 177 linesupings.py - tools/
build_models_manifest.py , Python, 151 lines - tools/
build_rsa_models.py , Python, 274 lines - tools/
bulk_check.py , Python, 518 lines - tools/
check_space.py , Python, 208 lines - tools/
colab_gpu/ , Jupyter, 82 linescolab_rsa.ipynb - tools/
colab_gpu/ , Jupyter, 53 linescolab_rsa4.ipynb - tools/
colab_gpu/ , Jupyter, 246 linescolab_rsa_group.ipynb - tools/
colab_gpu/ , Jupyter, 119 linescolab_rsa_regression.ipy nb - tools/
colab_gpu/ , Jupyter, 103 linescolab_rsa_regression_inf erence.ipynb - tools/
colab_gpu/ , Python, 2,068 linesgpu_group.py - tools/
colab_gpu/ , Python, 279 linesgpu_regression.py - tools/
colab_gpu/ , Python, 1,095 linesgpu_rsa.py - tools/
colab_gpu/ , Python, 199 linesmemory_estimate.py - tools/
colab_gpu/ , Python, 107 linesrefs/ build_refs.py - tools/
colab_gpu/ , Python, 161 linesrefs/ check_refs_mask.py - tools/
colab_gpu/ , Python, 191 linesrun_colab.py - tools/
colab_gpu/ , Python, 321 linesrun_colab_group.py - tools/
colab_gpu/ , Python, 408 linesrun_colab_regression.py - tools/
colab_gpu/ , Python, 199 linesrun_colab_regression_inf erence.py - tools/
colab_gpu/ , Python, 84 linestest_model_subsets.py - tools/
colab_gpu/ , Python, 297 linesvalidate_gpu.py - tools/
colab_gpu/ , Python, 647 linesvalidate_group.py - tools/
create_group_package.py , Python, 221 lines - tools/
create_package.py , Python, 571 lines - tools/
create_regression_infere , Python, 57 linesnce_package.py - tools/
create_regression_packag , Python, 164 linese.py - tools/
dashboard.py , Python, 133 lines - tools/
fsf_design.py , Python, 515 lines, 1 match - tools/
glm_designer.py , Python, 746 lines - tools/
hypothesis_explorer.py , Python, 3,881 lines - tools/
make_mask.py , Python, 232 lines - tools/
make_qr.py , Python, 50 lines - tools/
models_manifest.py , Python, 398 lines - tools/
pipeline_console.py , Python, 1,357 lines - tools/
pipeline_dashboard.py , Python, 1,457 lines - tools/
purge_step1_from_package , Python, 124 lines.py - tools/
queue_scan.py , Python, 198 lines - tools/
schedule_steps.py , Python, 164 lines - tools/
set_live.py , Python, 36 lines - tools/
unpack_results.py , Python, 461 lines - tools/
zmap_summary.py , Python, 1,294 lines - utils.py, Python, 1,005 lines, 1 match
- viz/
__init__.py , Python, 26 lines - viz/
datasource.py , Python, 183 lines - viz/
hypothesis_tree.py , Python, 320 lines - viz/
niftiutil.py , Python, 457 lines - viz/
planner_app.py , Python, 110 lines - viz/
scheduler_app.py , Python, 344 lines - viz/
stimuli.py , Python, 66 lines - viz/
viewer_app.py , Python, 337 lines - README.md, Text, 56 lines
The paper's code and data availability statement is in the Data section.
Tracing map
Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.
What the map holds:
- 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
- 136 scripts, each with its path and the digest of its content;
- 5 matches between paragraphs of the paper and lines of the code (method lexical-v1);
- neither the text of the paper nor the code itself.
Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.
Data
No dataset and no data link were found in the paper.
Data and code availability
The data associated with this manuscript have been deposited at Zenodo, and the custom preprocessing and analysis code used in this study is available through the Dog Brain Toolkit on GitHub. Both resources are publicly available, and their DOIs are listed in the key resources table. Any additional information required to reanalyze the data reported in this work is available from the lead contact upon request.
Reproduced under the paper's license (CC BY), from the paper cited above.
Versions
The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.
Version 3, 28 September 2026
- Authors: added Raúl Hernández-Pérez (0000-0003-3971-6236); Laura V. Cuaya (0000-0002-1073-7300); removed Raúl Hernández-Pérez; Laura V. Cuaya
Version 1, 27 September 2026: the first record
Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 6 keywords, 13 funders, 90 references.
Cite
This paper
Hernández-Pérez, R., Concha, L., Andics, A., Bernal-Gamboa, R., & Cuaya, L. V. (2026). Dog brain representations of human facial expressions: Encoding happiness and differentiating specific negative expressions. iScience, 29(8), 116900. https://
BibTeX
@article{hernandezperez2
author = {Hernández-Pérez, Raúl and Concha, Luis and Andics, Attila and Bernal-Gamboa, Rodolfo and Cuaya, Laura V.},
title = {{Dog brain representations of human facial expressions: Encoding happiness and differentiating specific negative expressions}},
journal = {iScience},
year = {2026},
month = aug,
volume = {29},
number = {8},
pages = {116900},
publisher = {Elsevier},
issn = {2589-0042},
doi = {10.1016/
url = {https://
pmid = {42666977},
pmcid = {PMC13523845}
}
RIS
TY - JOUR
AU - Hernández-Pérez, Raúl
AU - Concha, Luis
AU - Andics, Attila
AU - Bernal-Gamboa, Rodolfo
AU - Cuaya, Laura V.
TI - Dog brain representations of human facial expressions: Encoding happiness and differentiating specific negative expressions
T2 - iScience
J2 - iScience
PY - 2026
DA - 2026/
VL - 29
IS - 8
SP - 116900
SN - 2589-0042
PB - Elsevier
DO - 10.1016/
UR - https://
LA - en
ER -
CSL-JSON
{
"id": "10.1016/
"type": "article-journal",
"title": "Dog brain representations of human facial expressions: Encoding happiness and differentiating specific negative expressions",
"container-title": "iScience",
"author": [
{
"family": "Hernández-Pérez",
"given": "Raúl"
},
{
"family": "Concha",
"given": "Luis"
},
{
"family": "Andics",
"given": "Attila"
},
{
"family": "Bernal-Gamboa",
"given": "Rodolfo"
},
{
"family": "Cuaya",
"given": "Laura V."
}
],
"container-title-short":
"volume": "29",
"issue": "8",
"page": "116900",
"DOI": "10.1016/
"PMID": "42666977",
"PMCID": "PMC13523845",
"ISSN": "2589-0042",
"publisher": "Elsevier",
"URL": "https://
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
10
]
]
}
}
The tracing map gets a citation of its own once an author has validated it and it has a DOI.
Similar papers
The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.
- [1] doi:10.1038/s41467-026-73668-y [code]
- Convergent and divergent brain-cognition development in early adolescence.Journal: Nature communicationsIn common: Tools for NIfTI and ANALYZE image (MATLAB), FSL, Nilearn, 10 other tools, fMRI, 1 reference
- [2] doi:10.21203/rs.3.rs-9326213/v1 [code]
- Multi-task fMRI outperforms resting-state fMRI for revealing task-invariant organization of the human brainJournal: Research Square (preprint)In common: FSL, Nilearn, Plotly, 8 other tools, fMRI, cognitive, 2 references
- [3] doi:10.64898/2026.03.09.710558 [code]
- Multi-task fMRI outperforms resting-state fMRI for revealing task-invariant organization of the human brainJournal: bioRxiv (preprint)In common: FSL, Nilearn, Plotly, 8 other tools, fMRI, cognitive, 2 references
- [4] doi:10.1016/j.xcrm.2026.102943 [code]
- Parent-of-origin effects in Alzheimer's liability dissociate neurocognitive and cardiovascular traits in at-risk individuals.Journal: Cell reports. MedicineIn common: Tools for NIfTI and ANALYZE image (MATLAB), FSL, Nilearn, 10 other tools
- [5] doi:10.1002/hbm.70483 [code]
- Untamed: Unconstrained Tensor Decomposition and Graph Node Embedding for Cortical Parcellation.Journal: Human brain mappingIn common: Tools for NIfTI and ANALYZE image (MATLAB), FSL, Nilearn, 9 other tools, fMRI
- [6] doi:10.1038/s41586-026-10631-3 [code]
- A prognostic human brain network for diffuse midline glioma.Journal: NatureIn common: Tools for NIfTI and ANALYZE image (MATLAB), FSL, Nilearn, 9 other tools, 1 reference
- [7] doi:10.1016/j.neuron.2026.04.011 [code]
- Precision fMRI reveals densely interdigitated network patches with conserved motifs in the lateral prefrontal cortex.Journal: NeuronIn common: Tools for NIfTI and ANALYZE image (MATLAB), FSL, Nilearn, 8 other tools, fMRI, 1 reference
- [8] doi:10.3389/fnagi.2026.1742371 [code]
- Cross-sectional and longitudinal functional network alterations associated with subthreshold depressive symptoms in healthy older adults.Journal: Frontiers in aging neuroscienceIn common: Tools for NIfTI and ANALYZE image (MATLAB), FSL, Nilearn, 8 other tools, fMRI, 1 reference
- [9] doi:10.1038/s41592-026-03154-2 [code]
- Simultaneous single-cell calcium imaging of neuronal population activity and brain-wide BOLD fMRI.Journal: Nature methodsIn common: FSL, Nilearn, Image Processing Toolbox, 8 other tools, fMRI, 2 references
- [10] doi:10.1038/s41467-026-76452-0 [code]
- Music evokes shared neural representations of imagined narratives across sensory modalities.Journal: Nature communicationsIn common: Tools for NIfTI and ANALYZE image (MATLAB), FSL, Nilearn, 9 other tools, cognitive
Contribute
The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.
Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.
Claim this paper
Correct its record
Say what each link of this record is, remove the ones that are not the paper's, add the ones that are missing. The correction becomes a new version of the record, in its Versions section.
Validate its tracing map
You validate the map as this page shows it: 1 repository of the authors' code, each at its verified commit and with its license, 136 scripts, and 5 matches between paragraphs and code (see the Code and Map sections). It then receives a DOI on Zenodo, with you (your ORCID iD) and OSCR as its creators; the code itself is not deposited.
The map's fingerprint: sha256:df7ca280d33e3796…
Add the badge to its README
The badge links the code to this page. Copy one of these into the README of the paper's code: only you decide where it goes, and nothing is changed for you.
Markdown
[, paste the snippet at the top, then “Commit changes…” and, to review it first, “Create a new branch and start a pull request”. You open the pull request; OSCR asks for no permission.
Request its removal
To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).
Discussion, reproductions, activity
Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.
Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.
Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.
