%%% Pipeline MEP Mobi clear all close all %main_dir = fileparts(mfilename('fullpath')); % path of the script folder cd(main_dir) %% paths for PLUG-INs addpath([main_dir '\PLUG-INs\SOUND']);%SOUND addpath([main_dir '\PLUG-INs\TESA1.1.1']); %TESA addpath([main_dir '\PLUG-INs\M_functions']); % specific functions for the script %path of the raw data %addpath('/raw_data'); % initialize eeglab; addpath ('eeglab2022.1/') %% add path eeglab eeglab; close; %initialize eeg lab %initialize fieldtrip ft_defaults %% set folders for input and output: Y %%%% folder with data to be processed % datafolder_parent = %[main_dir '\mep_data']; % parent data folder includes all the subdir of MEP dataset to be processed % mkdir(main_dir, '\analysis_MEP'); % build output folder % output = [main_dir, '\analysis_MEP']; %% load tables to transcribe them Load_tables_MEP cd(main_dir) %% Variables ref_emg=[2,4]; epoching_long=[-0.2 0.5];% interval to epoch data, latency in sec relative to time-locking event epoching_long_VI=[-0.01 0.06];% interval to epoch data for visual inspection rmBase = [-100 -2]; % interval for baseline correction SR = 9600; %%step1.2 interp = [0 3]; %interval to cut and interpolate TMS pulse for TESA method, time in ms related to the event interp_Win = [0.35,0.35];%times before and after artifact window for fitting cubic function(TESA) HP_f = 1; % % HP filter frequency fftflag=1; % filterinf in the frequency domain, faster for high-order filters DS= 4800; %down sampling frequency tot_cond=7; %% 0. get data [subjects,n_subjects]= get_data_MOB_1(datafolder_parent,output,ref_emg,tot_cond); save('subjects.mat', 'subjects') save('n_subjects.mat', 'n_subjects') %% 0.1 HD OFF %uncomment these two lines if n_subjects already exist - and comment part %0. get data %load n_subjects.mat %load subjects.mat %% 1 for iSub = 1:n_subjects load n_subjects.mat load subjects.mat subj = subjects{iSub}; load([output '\' subj '\0_Data\Block_cond.mat']) for n_cond= 1:tot_cond name_condition = [subj '_' Block_cond{n_cond,2}]; name_folder= [subj '\' name_condition]; mkdir([output, '\' name_folder '\1_Preprocessing\step1_1']); current_output_folder = [output, '\' name_folder '\1_Preprocessing\step1_1']; if isfile([ current_output_folder '\' Block_cond{n_cond,1} '_EMG_step1_1.set']) ... %&& global_overwrite==0 && prep1_1_overwrite==0 warning(['skipping step1_1 for ' name_condition ' , file already present']) else load([output '\' subj '\0_Data\' subj '_' Block_cond{n_cond,2}]); data = eval([Block_cond{n_cond,2}]); EMG = pop_importdata... ('dataformat','matlab',... 'nbchan',height(data),... 'data', 'data',... 'setname', [subj '_rawmerged'],... 'srate',SR,... 'pnts',0,... 'xmin',0); EMG = eeg_checkset( EMG ); % using the trigger identified from TMS artifact % Read triggers from channel and then remove the channel and events % different from trigger EMG = pop_chanevent(EMG, length(EMG.data(:,1)),'edge',Block_cond{n_cond,4},'edgelen',0,'delchan','on','typename', num2str(Block_cond{n_cond,5})); EMG = pop_selectevent( EMG, 'type',str2num(Block_cond{n_cond,5}) ,'deleteevents','on'); %per S06 blocco 4: EMG = pop_selectevent( EMG, 'type',3 ,'deleteevents','on'); %S09 blocco 6: EMG = pop_selectevent( EMG, 'type',5 ,'deleteevents','on') EMG = pop_saveset(EMG, 'filename',[Block_cond{n_cond,1} '_EMG_step1_1.set'], 'filepath',current_output_folder ); end disp('Step 1.1 ended') % load EMG data x condition % interpolation %filters % downsampling % epoching %% STEP 1.2 Remove TMS-pulse artifact and interpolate, HP filter, dwnsampling,epoching previous_folder= [output, '\' name_folder '\1_Preprocessing\step1_1']; mkdir([output, '\' name_folder '\1_Preprocessing\step1_2']); current_output_folder = [output, '\' name_folder '\1_Preprocessing\step1_2']; if isfile([ current_output_folder '\' Block_cond{n_cond,1} '_step1_2.set']) warning(['skipping step1_2 for ' name_condition ' , file already present']) else EMG = pop_loadset('filename',[Block_cond{n_cond,1} '_EMG_step1_1.set'],'filepath', previous_folder); %output_folder = 'F:\Jack\MoBi\pipeline\analysis_MOBI_MEP'; % function to remove TMS pulse with TESA Method % inputs: EEG, interpolation interval,times before and after artifact window for fitting cubic function EMG= remove_TMS_TESA_G(EMG,interp,interp_Win); %plot with epoching and baseline correction only for visual purpose myplot_F(EMG,[-100 300],[-300 300],['After interpolation TESA'],epoching_long,rmBase); saveas(gcf, [current_output_folder '\after_interpolation_TESA']) close %%%%% High pass filter data % tic % EMG = high_pass_filter(EMG,HP_f,fftflag); % toc EMG = pop_eegfiltnew(EMG, 'locutoff',10, 'usefftfilt',1); EMG = pop_eegfiltnew(EMG, 'hicutoff',2500, 'usefftfilt',1); EMG = pop_eegfiltnew(EMG, 49.9,50.1,[],1,[],[],[],1); %plot after filtering without baseline correction myplot_F(EMG,[-100 300],[-500 500],['After filtering'],epoching_long); saveas(gcf, [current_output_folder '\after_filtering']) close % DownSample EMG = pop_resample(EMG, DS); myplot_F(EMG,[-100 300],[-500 500],['After resampling'],epoching_long); saveas(gcf, [current_output_folder '\after_resampling']) close % Epoching data EMG = pop_epoch(EMG, {}, epoching_long); %save figure myplot_F(EMG,[-100 300],[-500 500],['After Filtering']); saveas(gcf,[current_output_folder '\after_filtering']) close savename = [Block_cond{n_cond,1} '_step1_2.set']; EMG = pop_saveset( EMG, 'filename',savename,'filepath',current_output_folder); end % save set file which is closer to MEPs time range for latency visual inspection %if isfile([ current_output_folder '\' Block_cond{n_cond,1} '_step1_2_VisualInspection.set']) %warning(['skipping step1_2 for ' name_condition ' , file already present']) %else EMG = pop_loadset('filename',[Block_cond{n_cond,1} '_step1_2.set'],'filepath',current_output_folder); EMG2 = pop_epoch(EMG, {}, epoching_long_VI); savename = [Block_cond{n_cond,1} '_step1_2_VisualInspection.set']; EMG2 = pop_saveset( EMG2, 'filename',savename,'filepath',current_output_folder); % end disp('Step 1.2 ended') end end %% STEP 2. ESTIMATE MEP AMPLITUDE x CONDITIONS %TW = cell2table(cell(n_subjects,7), 'VariableNames', {'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'}); for iSub = 1:n_subjects load n_subjects.mat load subjects.mat subj = subjects{iSub} load([output '\' subj '\0_Data\Block_cond.mat']) for n_cond= 1:tot_cond name_condition = [subj '_' Block_cond{n_cond,2}]; name_folder= [subj '\' name_condition]; previous_folder= [output, '\' name_folder '\1_Preprocessing\step1_2']; mkdir([output, '\' name_folder '\2_MEP_AmplitudeP2P']); current_output_folder = [output, '\' name_folder '\2_MEP_AmplitudeP2P']; % if isfile([ previous_folder '\' Block_cond{n_cond,1} '_EMG_step1.set']) ... % %&& global_overwrite==0 && prep1_1_overwrite==0 % warning(['skipping step2 for ' name_condition ' , file already present']) % else EMG = pop_loadset('filename',[Block_cond{n_cond,1} '_step1_2.set'],'filepath',previous_folder); condition = Block_cond{n_cond,1}; numtrials = length(EMG.data(1,1,:)); nomep_trials = 0; handchan = 2; %RHAND p2p_mep = []; for itrial = 1:numtrials %% 3.4) amplitude estimation a = EMG.data(handchan,:,itrial); %channel x time point x current trial start_mep = 19; end_mep = 77; inizio_mep = find(EMG.times>=start_mep,1,"first"); fine_mep = find(EMG.times>=end_mep,1,"first"); z = a(inizio_mep:fine_mep); [minV,tpmin]= min(z); [maxV,tpmax] = max(z); %time window time_window = EMG.times(inizio_mep:fine_mep); %points in ms tmin = time_window(tpmin); tmax = time_window(tpmax); % peak2peak amplitude p2p = sprintf('%10.2f',(maxV - minV)); p2p(p2p==' ') =[]; p2p = str2num(p2p); %check trials without MEPs (exclude trial <50µV) p2p50 = p2p; if p2p50 <= 50 p2p50 = NaN; nomep_trials = nomep_trials + 1; end %amplitude x trial p2p_mep(itrial) = p2p50; p2p_all(itrial) = p2p; end p2p_meptrials = p2p_mep; p2p_allamp = p2p_all; %save([current_output_folder '\' Block_cond{iSub,1} '_' Block_cond{iSub,2} '_p2p_meptrials.mat'],'p2p_meptrials'); %save([current_output_folder '\' Block_cond{iSub,1} '_' Block_cond{iSub,2} '_p2p_allamp.mat'],'p2p_allamp'); save([current_output_folder '\p2p_meptrials.mat'],'p2p_meptrials'); save([current_output_folder '\p2p_allamp.mat'],'p2p_allamp'); %waitfor(gcf) % output tables if contains(condition,"Mo_PA") MEP_p2pAmp{iSub,1} = p2p_mep; nomep_sub(iSub, 1) = nomep_trials; MEP_mean(iSub,1) = nanmean(p2p_mep); elseif contains(condition,"Mo_AP") MEP_p2pAmp{iSub,2} = p2p_mep; nomep_sub(iSub, 2) = nomep_trials; MEP_mean(iSub,2) = nanmean(p2p_mep); elseif contains(condition,"Mo_LM") MEP_p2pAmp{iSub,3} = p2p_mep; nomep_sub(iSub, 3) = nomep_trials; MEP_mean(iSub,3) = nanmean(p2p_mep); elseif contains(condition,"Bi_PA_Contr") MEP_p2pAmp{iSub,4} = p2p_mep; nomep_sub(iSub, 4) = nomep_trials; MEP_mean(iSub,4) = nanmean(p2p_mep); elseif contains(condition,"Bi_AP") MEP_p2pAmp{iSub,5} = p2p_mep; nomep_sub(iSub, 5) = nomep_trials; MEP_mean(iSub,5) = nanmean(p2p_mep); elseif contains(condition,"Bi_LM") MEP_p2pAmp{iSub,6} = p2p_mep; nomep_sub(iSub, 6) = nomep_trials; MEP_mean(iSub,6) = nanmean(p2p_mep); elseif contains(condition,"Bi_PA") MEP_p2pAmp{iSub,7} = p2p_mep; nomep_sub(iSub, 7) = nomep_trials; MEP_mean(iSub,7) = nanmean(p2p_mep); %aggiungo questo per S38 % elseif contains(condition,"BI_PA") % MEP_p2pAmp{iSub,7} = p2p_mep; % nomep_sub(iSub, 7) = nomep_trials; % MEP_mean(iSub,7) = nanmean(p2p_mep); end %conditions in randomization order: MEP_p2pAmp_rc{iSub,n_cond} = p2p_mep; MEP_p2pAll_rc{iSub,n_cond} = p2p_all;%all amplitudes (also lower than 50uV) % %end end cd(output) TWMean(iSub,:) = array2table(MEP_mean(iSub,:),'VariableNames',{'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'}) TWP2P(iSub,:) = cell2table(MEP_p2pAmp(iSub,:),'VariableNames',{'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'}) TWNoM(iSub,:) = array2table(nomep_sub(iSub,:),'VariableNames',{'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'}) TQAMP(iSub,:) = cell2table(MEP_p2pAmp_rc(iSub,:),'VariableNames',{'B_1','B_2','B_3','B_4','B_5','B_6','B_7'}) TQALL(iSub,:) = cell2table(MEP_p2pAll_rc(iSub,:),'VariableNames',{'B_1','B_2','B_3','B_4','B_5','B_6','B_7'}) end %aggiungo parentesi e aggiungo nel loop % TWMean = array2table(MEP_mean,'VariableNames',{'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'}) % TWP2P = cell2table(MEP_p2pAmp,'VariableNames',{'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'}) % TWNoM = array2table(nomep_sub,'VariableNames',{'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'}) % TQAMP = cell2table(MEP_p2pAmp_rc,'VariableNames',{'B_1','B_2','B_3','B_4','B_5','B_6','B_7'}) % TQALL = cell2table(MEP_p2pAll_rc,'VariableNames',{'B_1','B_2','B_3','B_4','B_5','B_6','B_7'}) save([output '\nomep_sub.mat'],'TWNoM'); %numero trial senza MEP save([output '\MEP_p2pAmp.mat'],'TWP2P');% ampiezza mep picco-picco %ordine: {'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'} save([output '\MEP_mean.mat'],'TWMean'); %media dei trial x soggetto x condizione %ordine: {'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'} save([output '\MEP_p2pAmp_randomorder.mat'],'TQAMP'); % {'B_1','B_2''B_3''B_4''B_5''B_6''B_7'} save([output '\MEP_p2pAll_randomorder.mat'],'TQALL'); %% STEP 3. ESTIMATE MEP LATENCY x CONDITIONS % main_dir = 'C:\Users\delia\OneDrive\Documenti\MSc-cognitive neuroscience\Thesis\MEP'; %datafolder_parent = 'C:\Users\delia\OneDrive\Documenti\MSc-cognitive neuroscience\Thesis\MEP'; %'F:\Jack\MoBi\Data_MEP'; %[main_dir '\mep_data']; % parent data folder includes all the subdir of MEP dataset to be processed for iSub = 1:n_subjects load n_subjects.mat load subjects.mat subj = subjects{iSub}; load([output '\' subj '\0_Data\Block_cond.mat']) for n_cond = 1:tot_cond name_condition = [subj '_' Block_cond{n_cond,2}]; name_folder= [subj '\' name_condition]; previous_folder= [output, '\' name_folder '\1_Preprocessing\step1_2']; %load amplitude mat file to get NaN amplitudes amplitude_folder= [output, '\' name_folder '\2_MEP_AmplitudeP2P']; load([amplitude_folder,'\p2p_allamp.mat']) %create new folder for latency mkdir([output, '\' name_folder '\3_MEP_Latency']); current_output_folder = [output, '\' name_folder '\3_MEP_Latency']; %create new folder for plots mkdir([output, '\' name_folder '\3_MEP_Latency\plots']); plots_output_folder = [output, '\' name_folder '\3_MEP_Latency\plots']; % if isfile([ previous_folder '\' Block_cond{n_cond,1} '_EMG_step1.set']) ... % %&& global_overwrite==0 && prep1_1_overwrite==0 % warning(['skipping step2 for ' name_condition ' , file already present']) % else %load set file EMG = pop_loadset('filename',[Block_cond{n_cond,1} '_step1_2.set'],'filepath',previous_folder); %cut epochs from 10ms to 50ms epoching_long=[0.01 0.05];% interval to epoch data, latency in sec relative to time-locking event EMG = pop_epoch(EMG, {}, epoching_long); %variables condition = Block_cond{n_cond,1}; numtrials = length(EMG.data(1,1,:)); handchan = 2; lat_mep = []; lat_all = []; for itrial = 1:numtrials data = EMG.data(handchan,:,itrial); %channel x time point x current trial data=data'; time=EMG.times'; % reverse data if the first peak is negative [MA, IA] = max(data); [MI, II] = min(data); if II0; % to find points of positive derivative (0 non positive; 1 positive) % Label each separate region/cluster of ones. [labeledX, numRegions] = bwlabel(cons); % Get lengths of each region props = regionprops(labeledX, 'Area', 'PixelList'); regionLengths = [props.Area]; k = find(regionLengths==max(regionLengths)); % to find which is the biggest cluster %if one max cluster is found, calculate MEP onset and plot it if length(k)==1 fprintf('Region #%d is %d elements long and starts at element #%d\n',... k, regionLengths(k), props(k).PixelList(1,2)); start= props(k).PixelList(1,2); %index of start for longest cluster time(start) % MEP onset % to visualize MEPonset =zeros(size(z)); %to create a second channel MEPonset(start)=100; % to plot figure plot(cat(2, z, der1, MEPonset)); %to save the plot ax = gca; numtrial=num2str(itrial); %if the trial has no mep, name the image with "NO MEP" if p2p_allamp(itrial) <= 50 exportgraphics(ax,[plots_output_folder,'\latplot',numtrial,'_NO_MEP.jpg']) elseif p2p_allamp(itrial) > 50 exportgraphics(ax,[plots_output_folder,'\latplot',numtrial,'.jpg']) end v=EMG.times(1:IA); latency=v(start); disp([ ' MEP latency in condition ', condition, ', trial ', num2str(itrial), ', is ' num2str(latency), ' ms' ]) latency50=latency; %if it is not a MEP, save it as NaN if p2p_allamp(itrial) <= 50 latency50 = NaN end %if more than 1 max clusters are found, take into consideration the last one (in terms of time) elseif length(k)>1 k=max(k); fprintf('Region #%d is %d elements long and starts at element #%d\n',... k, regionLengths(k), props(k).PixelList(1,2)); start= props(k).PixelList(1,2); %index of start for longest cluster time(start) % MEP onset % to visualize MEPonset =zeros(size(z)); %to create a second channel MEPonset(start)=100; % to plot figure plot(cat(2, z, der1, MEPonset)); %to save the plot ax = gca; numtrial=num2str(itrial); %all images here are saved with "2P" (two peaks) to check them later on if p2p_allamp(itrial) <= 50 exportgraphics(ax,[plots_output_folder,'\latplot',numtrial,'_NO_MEP_2P.jpg']) elseif p2p_allamp(itrial) > 50 exportgraphics(ax,[plots_output_folder,'\latplot',numtrial,'_2P.jpg']) end v=EMG.times(1:IA); latency=v(start); disp([ ' MEP latency in condition ', condition, ', trial ', num2str(itrial), ', is ' num2str(latency), ' ms' ]) latency50=latency; if p2p_allamp(itrial) <= 50 latency50 = NaN end %if k is empty, the peak cannot be computed elseif isempty(k)==1 warning(['skipping latency estimation for ' name_condition ' , missing minimum input arguments']) latency=NaN; latency50=NaN; end lat_mep(itrial)=latency50; lat_all(itrial)=latency; end lat_meptrials=lat_mep; lat_alltrials=lat_all; save([current_output_folder '\lat_meptrials.mat'],'lat_meptrials'); save([current_output_folder '\lat_alltrials.mat'],'lat_alltrials'); if contains(condition,"Mo_PA") MEP_lat{iSub,1} = lat_mep; %nomep_sub(iSub, 1) = nomep_trials; MEAN_lat(iSub,1) = nanmean(lat_mep); elseif contains(condition,"Mo_AP") MEP_lat{iSub,2} = lat_mep; %nomep_sub(iSub, 2) = nomep_trials; MEAN_lat(iSub,2) = nanmean(lat_mep); elseif contains(condition,"Mo_LM") MEP_lat{iSub,3} = lat_mep; %nomep_sub(iSub, 3) = nomep_trials; MEAN_lat(iSub,3) = nanmean(lat_mep); elseif contains(condition,"Bi_PA_Contr") MEP_lat{iSub,4} = lat_mep; %nomep_sub(iSub, 4) = nomep_trials; MEAN_lat(iSub,4) = nanmean(lat_mep); elseif contains(condition,"Bi_AP") MEP_lat{iSub,5} = lat_mep; %nomep_sub(iSub, 5) = nomep_trials; MEAN_lat(iSub,5) = nanmean(lat_mep); elseif contains(condition,"Bi_LM") MEP_lat{iSub,6} = lat_mep; %nomep_sub(iSub, 6) = nomep_trials; MEAN_lat(iSub,6) = nanmean(lat_mep); elseif contains(condition,"Bi_PA") MEP_lat{iSub,7} = lat_mep; %nomep_sub(iSub, 7) = nomep_trials; MEAN_lat(iSub,7) = nanmean(lat_mep); %aggiungo questo per S38 elseif contains(condition,"BI_PA") MEP_lat{iSub,7} = lat_mep; %nomep_sub(iSub, 7) = nomep_trials; MEAN_lat(iSub,7) = nanmean(lat_mep); end %conditions in randomization order: MEP_lat_rc{iSub,n_cond} = lat_mep; MEP_latAll_rc{iSub,n_cond} = lat_all;%all amplitudes (also lower than 50uV) end cd(output) TOMeanLat(iSub,:) = array2table(MEAN_lat(iSub,:),'VariableNames',{'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'}) TOLat(iSub,:) = cell2table(MEP_lat(iSub,:),'VariableNames',{'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'}) TRLat(iSub,:) = cell2table(MEP_lat_rc(iSub,:),'VariableNames',{'B_1','B_2','B_3','B_4','B_5','B_6','B_7'}) TRALLat(iSub,:) = cell2table(MEP_latAll_rc(iSub,:),'VariableNames',{'B_1','B_2','B_3','B_4','B_5','B_6','B_7'}) end save([output '\MEP_Lat.mat'],'TOLat');% latenza mep %ordine: {'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'} save([output '\MEP_meanLat.mat'],'TOMeanLat'); %media dei trial x soggetto x condizione %ordine: {'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'} save([output '\MEP_Lat_randomorder.mat'],'TRLat'); % {'B_1','B_2''B_3''B_4''B_5''B_6''B_7'} save([output '\MEP_AlLat_randomorder.mat'],'TRALLat');