Upload New File

MRI_processing/fMRI_Preprocessing_Nback.m

+%% FUNCTION THAT PREPROCESSES FUNCTIONAL DATA FOR THE Nback TASK
+
+
+% VERSION_02 : This scrip applyies the following modules to process EPI images: 
+
+% 1- Calculate VDM
+% 2- Realign & Unwarp
+% 3- Brain extraction
+% 4- Smoothing. 
+
+% "subs_info": Contains subject ID and subjuct session name:
+% sub-001/ses-morning
+% "T": contains all info from the "HILIGHT_participants.xlsx" 
+% "paths": path to data and all toolboxes needed
+% "overwrite" : if set to "true" : the code will re-run the preprocessing of
+% functional data even if an output already exists ; if set to "false", 
+% the code will skip preprocessing if an output already exists.  
+
+% --> The preprocessed images are saved in the output_preprocessing folder,
+% in the subfolder corresponding to the task. For example, nback results
+% will be stored in output_preprocessing/nback.
+
+
+% --> BE CAREFUL : this code only reads "prepared" data. You need to have
+% the data downloaded on your computer in the prepared format:
+
+%..../prepared
+%         --> sub-008
+%               --> fMRI
+%               --> Phsio
+%               --> behav
+%               --> ses-morning
+%                    --> MRI
+%                        --> 005-b1map_sag_p2
+%                        --> 007-cmrr_mbep2d_p3_mb2_1.4iso_VolCheck 
+%                        --> ...
+%        --> sub-009
+
+
+% Written by R. Sharifpour- A-Berger
+% Cyclotron Research Centre, University of Liege, Belgium
+% 18th November 2022
+
+%%
+
+function fMRI_Preprocessing_Nback_v02(subs_info,T,paths,overwrite)
+for id = 1:length(subs_info)
+    
+    %this part of the code search for the functional files (nback task) of the
+    %current subject being processed
+    
+    subject_info = char(subs_info(id));
+    subject_ID = subject_info(5:7);
+    switch subject_info(13)
+        case 'm'
+            subject_session = 'M';
+        case 'e'
+            subject_session = 'E';
+    end
+    
+    
+    entry_table = ['HILIGHT_', subject_ID '_' subject_session] ;
+    index = strcmp(T.subject_Ids, entry_table) ;
+    patient_data = T(index,:) ;
+    
+    message = ['processing functional images of subject: sub-', subject_ID '/' subject_info(9:end)] ;
+    disp(message)
+    
+    fprintf('Looking for NIFTI functional files (TASK : NBACK)')
+    
+    %We verify that a sequence is encoded in the excel sheet, else next
+    %subject is preprocessed
+    nback_sequence = [str2num(patient_data.Nback{:})] ;
+    nback_sequence = isempty(nback_sequence) ;
+    if nback_sequence==1
+        message = ['(1) No data found for sub-', subject_ID '/' subject_info(9:end), 'or (2) No corresponding functional sequence found in the Excel table'] ;
+        disp(message)
+        continue
+    end
+    
+    prep_folder = dir(char([paths.data])) ;
+    
+    for iFile = 1:numel(prep_folder)
+        if(strncmpi(prep_folder(iFile).name,['sub-', subject_ID],7))
+            sub_folder = fullfile(char([paths.data]), prep_folder(iFile).name); % .../Prepared/sub-...
+        end
+    end
+    
+    mri_folder = dir(char([sub_folder subject_info(8:end) '/MRI'])) ; % .../Prepared/sub-.../ses-.../MRI
+    for iFile = 1:numel(mri_folder)
+        if (endsWith(mri_folder(iFile).name,'nBack'))
+            nback_folder = fullfile(char([sub_folder subject_info(8:end) '/MRI']), mri_folder(iFile).name); % .../Prepared/sub-.../ses-.../MRI/...-cmrr_mbep2d_p3_mb2_1.4iso_nBack
+        end
+    end
+    
+    %We extracted the name of the .nii files but since the input in the
+    %batch needs the entire path to the nifti files, we will change the
+    %content of the cell containing the names of the files to add the
+    %entire path before the name
+    
+    nback_folder_data = dir(nback_folder) ; % all files in the nback folder
+    nback_nii_data = dir(fullfile(nback_folder, '*.nii')) ; % all images in the nback folder (struct)
+    nback_nii_data = {nback_nii_data.name} ; % all images in the nback folder (Cell)
+    [r c] = size(nback_nii_data) ;
+    
+    %We save the entire path to the nifti files (with fullpath)
+    for i = 1:c
+        nback_nii_data{i} = char([nback_folder '/' nback_nii_data{i}]) ;
+    end
+    
+    %We will remove from this list the dummy scan at the beginning of the
+    %EPI acquisition. In order to do so, we have to load the excel table
+    %with the exact number of volumes that we won't preprocess.
+    
+    dummyscans_nback = patient_data.Var28;
+    dummyscans_nback = str2num(dummyscans_nback{:}) ;
+    
+    nback_nii_data(1:dummyscans_nback) = [] ;
+    firstEPI=nback_nii_data(1);
+    
+    %Here, we reshape the cell into a format readable by spm. The input is
+    %put in a single 1x1 cell containing #volume x 1 cells.
+    nback_nii_data = nback_nii_data' ;
+    nback_nii_data = {nback_nii_data} ;
+    fprintf('\n Functional data for the Nback task found \n ')
+    
+    %If not created yet, a new folder "ouput_preprocessing/nback" will be created
+    %in the ".../prepared/sub-ID/ses-sesNAME/MRI/"  path in which the preprocessed EPIs
+    %will be saved.
+    output_dir = char([sub_folder subject_info(8:end) '/MRI' '/output_preprocessing']) ;%we save that path to create an output folder
+    output_preprocessed_nback = char([output_dir '/nback']) ;
+    
+    if ~isempty(dir(output_preprocessed_nback)) && overwrite == false
+        message =  ['Output already exists and "overwrite" set to false. No preprocessing conducted for sub-', subject_ID] ;
+        disp(message)
+        continue
+    end
+    
+    if ~exist(output_dir)
+        mkdir(output_dir)
+        mkdir(output_preprocessed_nback)
+    elseif ~exist(output_preprocessed_nback)
+        mkdir(output_preprocessed_nback)
+    end
+    
+    %%
+ 
+    %STEP 1 : CALCULATE VDM
+    
+    %First, we access the sequence number corresponding to phase and magnitude in the
+    %Excel sheet and then extract the nifti files in the corresponding
+    %folders. Then the batch for VDM calculation will be created.
+    
+    phase_mag = patient_data.Gre_field_map ;
+    phase_mag = strsplit(phase_mag{1}, ',') ;
+    magnitude_index = str2num(phase_mag{1}) ;
+    phase_index = str2num(phase_mag{2}) ;
+    
+    magnitude_index = sprintf('%03d', magnitude_index) ;
+    phase_index = sprintf('%03d', phase_index) ;
+    
+    for iFile = 1:numel(mri_folder)
+        if(strncmpi(mri_folder(iFile).name,magnitude_index,numel(magnitude_index)) && endsWith(mri_folder(iFile).name,'eufind'))
+            File_mag = fullfile([mri_folder(iFile).folder,'/', mri_folder(iFile).name]) ;
+        end
+        if(strncmpi(mri_folder(iFile).name,phase_index,numel(phase_index)) && endsWith(mri_folder(iFile).name,'eufind'))
+            File_phase = fullfile([mri_folder(iFile).folder,'/', mri_folder(iFile).name]) ;
+        end
+    end
+    
+    magnitude = dir(File_mag) ;
+    phase = dir(File_phase) ;
+    string_0_2_nifti = '02.nii' ;
+    
+    for iFile = 1:numel(magnitude)
+        if(strncmp(magnitude(iFile).name(end:-1:1),string_0_2_nifti(end:-1:1),numel(string_0_2_nifti)) & strncmpi(magnitude(iFile).name,'s',numel('s')))
+            magnitude_nifti = fullfile(File_mag, magnitude(iFile).name) ;
+        end
+    end
+    
+    for iFile = 1:numel(phase)
+        if(strncmp(phase(iFile).name(end:-1:1),string_0_2_nifti(end:-1:1),numel(string_0_2_nifti)) & strncmpi(phase(iFile).name,'s',numel('s')))
+            phase_nifti = fullfile(File_phase, phase(iFile).name) ;
+        end
+    end
+    
+    first_EPI=nback_nii_data{1}{1};
+    clear matlabbatch
+    spm_jobman('initcfg') ;
+    
+    
+    matlabbatch{1}.spm.tools.fieldmap.calculatevdm.subj.data.presubphasemag.phase = cellstr([phase_nifti]);
+    matlabbatch{1}.spm.tools.fieldmap.calculatevdm.subj.data.presubphasemag.magnitude = cellstr([magnitude_nifti]);
+    matlabbatch{1}.spm.tools.fieldmap.calculatevdm.subj.defaults.defaultsfile = cellstr([paths.spm '/toolbox/FieldMap/pm_defaults_7T_cmrr_ss_v1.m']) ;
+    matlabbatch{1}.spm.tools.fieldmap.calculatevdm.subj.session.epi(1) = cellstr([first_EPI]);%{firstEPI};% cellstr([nback_nifti{1}{1}]);
+    matlabbatch{1}.spm.tools.fieldmap.calculatevdm.subj.matchvdm = 1;
+    matlabbatch{1}.spm.tools.fieldmap.calculatevdm.subj.sessname = 'session';
+    matlabbatch{1}.spm.tools.fieldmap.calculatevdm.subj.writeunwarped = 0;
+    matlabbatch{1}.spm.tools.fieldmap.calculatevdm.subj.anat = '';
+    matlabbatch{1}.spm.tools.fieldmap.calculatevdm.subj.matchanat = 0;
+    
+    spm_jobman('run', matlabbatch) % run batch
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    %%
+    %STEP 2 : REALIGN & UNWARP
+    
+    %Registration of EPIs to the "meanEPI" and correction for voxel displacement(VD)
+    % using the VDM map calculated on the previous step.
+
+    phase = dir(File_phase);
+    for iFile = 3:numel(phase)
+        if(strncmp(phase(iFile).name(1:4),'vdm5',4))
+            VDM = fullfile(File_phase, phase(iFile).name);
+        end
+    end
+    
+    clear matlabbatch
+    spm_jobman('initcfg') ;
+    
+    matlabbatch{1}.spm.spatial.realignunwarp.data.scans = nback_nii_data{1};
+    matlabbatch{1}.spm.spatial.realignunwarp.data.pmscan = cellstr([VDM]);
+    matlabbatch{1}.spm.spatial.realignunwarp.eoptions.quality = 0.9;
+    matlabbatch{1}.spm.spatial.realignunwarp.eoptions.sep = 4;
+    matlabbatch{1}.spm.spatial.realignunwarp.eoptions.fwhm = 5;
+    matlabbatch{1}.spm.spatial.realignunwarp.eoptions.rtm = 1; % register to mean
+    matlabbatch{1}.spm.spatial.realignunwarp.eoptions.einterp = 2;
+    matlabbatch{1}.spm.spatial.realignunwarp.eoptions.ewrap = [0 0 0];
+    matlabbatch{1}.spm.spatial.realignunwarp.eoptions.weight = '';
+    matlabbatch{1}.spm.spatial.realignunwarp.uweoptions.basfcn = [12 12];
+    matlabbatch{1}.spm.spatial.realignunwarp.uweoptions.regorder = 1;
+    matlabbatch{1}.spm.spatial.realignunwarp.uweoptions.lambda = 100000;
+    matlabbatch{1}.spm.spatial.realignunwarp.uweoptions.jm = 0;
+    matlabbatch{1}.spm.spatial.realignunwarp.uweoptions.fot = [4 5];
+    matlabbatch{1}.spm.spatial.realignunwarp.uweoptions.sot = [];
+    matlabbatch{1}.spm.spatial.realignunwarp.uweoptions.uwfwhm = 4;
+    matlabbatch{1}.spm.spatial.realignunwarp.uweoptions.rem = 1;
+    matlabbatch{1}.spm.spatial.realignunwarp.uweoptions.noi = 5;
+    matlabbatch{1}.spm.spatial.realignunwarp.uweoptions.expround = 'Average';
+    matlabbatch{1}.spm.spatial.realignunwarp.uwroptions.uwwhich = [2 1];
+    matlabbatch{1}.spm.spatial.realignunwarp.uwroptions.rinterp = 4;
+    matlabbatch{1}.spm.spatial.realignunwarp.uwroptions.wrap = [0 0 0];
+    matlabbatch{1}.spm.spatial.realignunwarp.uwroptions.mask = 1;
+    matlabbatch{1}.spm.spatial.realignunwarp.uwroptions.prefix = 'u';
+    
+    
+    
+    spm_jobman('run', matlabbatch) % run batch
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    %%
+    
+    %STEP 3 : BRAIN EXTRACTION
+    
+    % First, brain will be extracted from the "meanEPI" and then the brain mask
+    % will be applied on the rest of the EPIs.
+    
+    
+    % MeanEPI brain extraction
+    
+    input1= dir(fullfile(nback_folder, 'meanuf*.nii'));
+    input1= [input1.folder '/' input1.name];% input image
+    input2= [output_preprocessed_nback '/BET_FUNC.nii']; %  brain extracted image name
+    input3= [output_preprocessed_nback '/BET_FUNC_MASK.nii']; % brain mask name
+    cmdStr=[paths.synthstrip '/synthstrip-docker -i ' input1 ' -o ' input2 ' -m ' input3];
+    system(cmdStr)
+    
+    % Apply the brain mask on all EPIs.
+    
+    EPIs = dir(fullfile(nback_folder, 'uf*.nii')) ;
+    cd (output_preprocessed_nback);
+    for ii=1:numel(EPIs)
+        input= [EPIs(ii).folder '/' EPIs(ii).name];
+        output = ['BET_' EPIs(ii).name];
+        cmdStr=['fslmaths BET_FUNC_MASK.nii -mul ' input ' ' output];
+        system(cmdStr)
+        cmdStr=['gunzip ' output];
+        system(cmdStr)
+        
+    end
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    %%
+    
+    %STEP 4 : SMOOTHING
+    
+    % Find all brain extracted EPIs
+    
+    uf_files = dir(fullfile(output_preprocessed_nback, 'BET_uf*.nii')) ;
+    
+    uf_files = {uf_files.name} ;
+    [r c] = size(uf_files) ;
+    
+    for i = 1:c
+        uf_files{i} = char([output_preprocessed_nback '/' uf_files{i}]) ;
+    end
+    
+    uf_files = uf_files' ;
+    
+    % Smooth all EPIs with a 3mm gaussian kernel
+    
+    clear matlabbatch
+    spm_jobman('initcfg') ;
+    
+    matlabbatch{1}.spm.spatial.smooth.data = uf_files;
+    matlabbatch{1}.spm.spatial.smooth.fwhm = [3 3 3];
+    matlabbatch{1}.spm.spatial.smooth.dtype = 0;
+    matlabbatch{1}.spm.spatial.smooth.im = 0;
+    matlabbatch{1}.spm.spatial.smooth.prefix = 's';
+    
+    spm_jobman('run', matlabbatch) % run batch
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    %%
+    
+    % Here, the file cointains motion correction parameters and final preprocessed EPIs 
+    % created in the folder of the initial(unprocessed) images are moved 
+    % into the output_preprocessing/nback folder created at the beginning.
+    
+    old_path = fullfile(nback_folder, 'rp*');
+    movefile(old_path,output_preprocessed_nback)
+    old_path = fullfile(nback_folder, 'meanuf*.nii') ;
+    movefile(old_path,output_preprocessed_nback)
+    
+    % All other images will be deleted from the initial folder
+    
+    old_path = fullfile(nback_folder, '*uf*');
+    delete(old_path,output_preprocessed_nback)
+    
+    
+    fprintf('________________________________________________________')
+    
+    
+end
+fprintf('\n done : all subjects have been processed \n')
+end