gin commit from psy497.lan

New files: 7

Pilot_analysis/GB_gtec_files_concatenation.m

@@ -0,0 +1,39 @@
+function  [alt_y, SR, MetaData] = GB_gtec_files_concatenation(file_tbc)
+
+    % get the last session (final number +1)
+    last_section = size(file_tbc,1);
+    
+    % 
+    alt_y = [];
+    
+    % for each chunk
+    for k = 1:last_section
+        % load the file
+        sect{k} = load(file_tbc(k).name);
+        
+        alt_y = horzcat(alt_y, sect{k}.y);
+        
+    end
+    
+    
+    % time channel is long as the sum of the timings
+    alt_y(1,:) = [1:size(alt_y,2)];
+    
+%     new_file = struct; 
+%     new_file.y = alt_y;
+%     new_file.SR = SR;
+%     
+%     if isfield(sect{end}, 'MetaData')
+%     new_file.MetaData = sect{end}.MetaData;
+%     end
+
+SR = sect{last_section}.SR; 
+
+if isfield(sect{last_section}, 'MetaData')
+    MetaData = sect{last_section}.MetaData;
+else
+    MetaData = '';
+end
+
+end
+

Pilot_analysis/JA_pilot_analysis.R

@@ -0,0 +1,127 @@
+
+rm(list = ls())
+
+
+library(openxlsx)
+library(dplyr) 
+library(ggplot2)
+library(tidyverse)
+library(ez)
+library(plotly)
+library(ggsignif)
+library(cowplot)
+library(multipanelfigure)
+library(ggh4x)
+library(svglite)
+library(BayesFactor)
+
+
+setwd("/Users/b1037210/Desktop/Progetti/JA_Stage1/")
+
+df = read.xlsx("final_table_for_R_10-2500_def_o1_rep.xlsx")
+
+
+df$TMS_timings_rep = as.factor(df$TMS_timings_rep)
+levels(df$TMS_timings_rep) = c("Late TMS","Early TMS")
+df$TMS_timings_rep <- factor(df$TMS_timings_rep, levels = c("Early TMS","Late TMS"))
+
+df$subject = as.factor(df$subject)
+df$CF_colors_order = as.factor(df$CF_colors_order)
+levels(df$CF_colors_order) = c("Lift", "Catch", "Press")
+
+df <- subset(df, CF_colors_order == "Press"| CF_colors_order == "Lift" ) 
+
+df$CF_colors_order = droplevels(df$CF_colors_order)
+
+groups <- group_by(df, CF_colors_order, subject, TMS_timings_rep) 
+
+# summarize by using the medians of each subject --> have a table with a "mean" column for each cell for each participant
+data_avgs <- summarise(groups,
+                       mean = median(ECD_MEP_amplitudes, na.rm=TRUE))
+
+# take a peek
+data_avgs = as.data.frame(data_avgs)
+
+df_final = data_avgs
+
+attach(df_final)
+
+
+## descriptive statistics
+s_alt = ezStats(data=df_final,
+                dv = .(mean),
+                wid = .(subject),
+                within = .(CF_colors_order, TMS_timings_rep))
+
+## plot now
+pd = position_dodge(width=0.0)
+
+# create basic plot with x axis as Cngruendy and y as mean (which is median actually)
+g = ggplot(df_final, aes(x = CF_colors_order, y = mean))
+
+# create Facets, the level of one variable above, the name of the levels of the other on the new line (\n)
+g = g+facet_grid(col = vars(interaction(TMS_timings_rep, sep = "\n",lex.order = TRUE)))
+
+# create bar plot 
+g = g+ geom_bar(aes(fill = CF_colors_order), stat = "summary", fun = "mean")
+
+
+g = g+ geom_line(aes(group = subject)
+                 , stat = "summary"
+                 , fun = "mean" 
+                 ,size= .3
+                 , position = pd)
+
+
+#new filling colors for the bars --> values = c("red","orange","blue","green"))
+g = g+scale_fill_manual(values = c("red","blue","tomato","dodgerblue"))
+
+g = g+theme_classic()
+
+g = g+theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) #top, right, bottom, left
+
+g = g+theme(strip.text.x = element_text(size=20, angle=0, vjust = 0),
+            strip.background = element_rect(colour="white", fill="white"))
+
+g = g + theme(
+  axis.title.x = element_text( size = 20, face = "bold", margin = margin(t = 10, r = 0, b = 0, l = 0)),
+  axis.title.y = element_text(size = 20, face = "bold", margin = margin(t = 0, r = 10, b = 0, l = 0)),
+  
+  axis.text.x = element_text(size = 15),# family = 'Helvetica',
+  axis.text.y = element_text(size = 15))  
+
+g = g+xlab("Cue instruction") + ylab("peak-to-peak amplitude (µV)") #ylab("TMS-evoked pressure (a.u.)") #
+
+g = g+theme(
+  legend.position = "none")
+
+
+file_name = "ECU" #"Pressure (z)"#"ECU(z) - FDS(z)"#"Flexor Digitorum Superficialis (Z-scores)\n"
+g = g+ggtitle(file_name)
+g = g+theme(plot.title = element_text(hjust = 0.5, size = 35, face = "bold"))
+
+
+g
+
+
+## stats
+late_lift = subset(data_avgs, TMS_timings_rep == "Late TMS" & CF_colors_order == "Lift")
+late_press = subset(data_avgs, TMS_timings_rep == "Late TMS" & CF_colors_order == "Press")
+
+early_lift = subset(data_avgs, TMS_timings_rep == "Early TMS" & CF_colors_order == "Lift")
+early_press = subset(data_avgs, TMS_timings_rep == "Early TMS" & CF_colors_order == "Press")
+
+shapiro.test(late_lift$mean)
+shapiro.test(late_press$mean)
+
+shapiro.test(early_lift$mean)
+shapiro.test(early_press$mean)
+
+t.test(early_lift$mean, early_press$mean, paired = TRUE, alternative = "greater")
+t.test(late_lift$mean, late_press$mean, paired = TRUE, alternative = "greater")
+
+
+
+## save
+file_title = paste0(file_name, ".jpg");
+ggsave(file_title,g,width = 10, height = 7, device = "jpg")

Pilot_analysis/JA_pilot_analysis_DEF_rep.m

@@ -0,0 +1,366 @@
+%%
+clear
+general_directory = '/Users/b1037210/Desktop/Progetti/JA_Stage1/'
+cd(general_directory)
+addpath(general_directory)
+%%
+addpath(genpath('/Users/b1037210/Documents/MATLAB/fieldtrip-20161117'))
+
+ft_defaults
+%%
+
+subjects = ...
+    {'S01_pilot',...
+    'S02_pilot',...
+    'S03_pilot',...
+    'S04_pilot',...
+    'S05_pilot',...
+    'S06_pilot',...
+    'S07_pilot',...
+    'S08_pilot',...
+    'S09_pilot',...
+    'S10_pilot',...
+    'S11_pilot',...
+    'S12_pilot'}
+
+
+threshold = {[0.018 0.032],...
+    [0.016 0.035],...
+    [0.018 0.040],...
+    [0.019 0.039],...
+    [0.017 0.037],...
+    [0.016 0.036],...
+    [0.019 0.042],...
+    [0.018 0.039],...
+    [0.016 0.035],...
+    [0.015 0.035],...
+    [0.02 0.03],...
+    [0.02 0.04]}
+
+
+
+TMS_timings = {'-0.05', '0.0'};
+final_table = table;
+%%
+for s = 1:length(subjects)
+    
+    %% cd directory
+    cd([general_directory,subjects{s}])
+    
+    %% ...for each tms timing 
+    for t = 1:length(TMS_timings)
+        
+        %% file loading
+        file_tmp = (dir(strcat(subjects{s},'_', TMS_timings{t},'_TMS*')))
+        
+        
+        if size(file_tmp,1)>1
+        
+        disp(strcat('concatenate ',file_tmp(1).name))
+        [y SR MetaData] = GB_gtec_files_concatenation(file_tmp);
+        
+        else
+        
+        disp(strcat('loading ', file_tmp.name)) 
+        load(file_tmp.name)
+
+        end
+        
+        clear file_tmp
+        %% name the trigger channel and time channel, then erase time from y
+        
+        % select trigger row and time row
+        trigger = y(end,:); % trigger channel for glove touch detection
+        time = y(1,:); % time-series in seconds
+        
+        % new data without time, so channels in the gtect settings and data
+        % analysis have the same nnumber
+        y = y(2:end,:);
+        
+        %% get triggers and sample points at which each event of interest happens
+        
+        % get were variations happen within trigger channel
+        diff_trigger = diff(trigger);
+        
+        % get trigger codes wihtin trgger line
+        trigger_numbers = unique(trigger);
+        
+        % find were the wash_out ends
+        end_wash_out = find(diff_trigger == 1)+1;
+        end_wash_out = end_wash_out(1);
+        
+        % for each code present, get the sample were the rising triggers are present
+        for x = 1:length(trigger_numbers)
+            
+            % find rising triggers and add 1 sampling point (due to diff which removes one)
+            tmp = find(diff_trigger == trigger_numbers(x))+1;
+            
+            % based on the code --> assign the tmp numbers to the correct variable
+            switch trigger_numbers(x)
+                
+%                 case 14 % color press CF
+%                     sp_slow_id = tmp(tmp>end_wash_out);
+%                     
+%                 case 15 % color press CF
+%                     sp_fast_id = tmp(tmp>end_wash_out);
+                
+                case 16 % color press CF
+                    sp_out_slow_id = tmp(tmp>end_wash_out);
+                
+                case 17 % color press CF
+                    sp_out_fast_id = tmp(tmp>end_wash_out);
+                    
+%                 case 18 % color press CF
+%                     sp_out_ok_catch_id = tmp(tmp>end_wash_out);
+%                      
+%                 case 19 % color press CF
+%                     sp_out_failed_catch_id = tmp(tmp>end_wash_out);
+%                 
+                    
+                    
+                    
+                    
+                case 20 % color press CF
+                    sp_color_CF_press = tmp(tmp>end_wash_out);
+                    
+                case 21 % color lift CF
+                    sp_color_CF_lift = tmp(tmp>end_wash_out);
+                    
+                case 22 % color catch CF
+                    sp_color_CF_catch = tmp(tmp>end_wash_out);
+                    
+                case 30 % color press SBJ
+                    sp_color_SBJ_press = tmp(tmp>end_wash_out);
+                    
+                case 31 % color_lift SBJ
+                    sp_color_SBJ_lift = tmp(tmp>end_wash_out);
+                    
+                case 40 % CF hitting
+                    sp_CF_hitting = tmp(tmp>end_wash_out);
+                    
+                case 41 % SBJ hitting
+                    sp_SBJ_hitting = tmp(tmp>end_wash_out);
+                    
+                case 128 % TMS pulse
+                    sp_TMS = tmp(tmp>end_wash_out);
+            end
+            
+        end
+        
+        %%
+        
+        % check dimensions
+        size(sp_color_CF_press)
+        size(sp_color_CF_lift)
+        size(sp_color_CF_catch)
+        size(sp_TMS)
+        
+        %
+        [CF_colors_idx] = sortrows(vertcat(sp_color_CF_press',sp_color_CF_lift',sp_color_CF_catch'));
+        CF_colors_order = CF_colors_idx;
+        
+        % assign codes to sample points corresponding to press lift and catch
+        CF_colors_order(ismember(CF_colors_idx, sp_color_CF_press)) = 1;
+        CF_colors_order(ismember(CF_colors_idx, sp_color_CF_lift)) = -1;
+        CF_colors_order(ismember(CF_colors_idx, sp_color_CF_catch)) = 0;
+        
+        %%
+        
+
+        size(sp_out_slow_id)
+        size(sp_out_fast_id)
+
+         
+  [trials_idx] = sortrows(vertcat(sp_out_slow_id',sp_out_fast_id'))
+
+ 
+         
+       %%
+       for ex = 1:length(trials_idx)
+          
+           [v_abs(ex) p_abs(ex)] = min(abs(CF_colors_idx - trials_idx(ex)))
+          ddd = CF_colors_idx - trials_idx(ex);
+           exclude(ex) = p_abs(ex);
+           is_exclude_positive(ex) = ddd(exclude(ex))
+       end
+       
+       clear v_abs
+       clear ddd
+       clear is_exclude_positive
+       
+       
+        
+        %% EMG preprocessing
+        
+        ECD_right = y(2,:);
+        FDS_right = y(4,:);
+        ECD_response = y(6,:);
+        
+        
+        
+        Fsample = SR;
+        hp_filter = 10
+        lp_filter = 2500
+        
+   
+        ECD_right_filt = ft_preproc_highpassfilter(ECD_right, Fsample, hp_filter,1);
+        FDS_right_filt = ft_preproc_highpassfilter(FDS_right, Fsample, hp_filter,1);
+        
+        ECD_right_filt = ft_preproc_lowpassfilter(ECD_right_filt, Fsample, lp_filter,1);
+        FDS_right_filt = ft_preproc_lowpassfilter(FDS_right_filt, Fsample, lp_filter,1);
+        
+        
+        ECD_right_filt = ft_preproc_dftfilter(ECD_right_filt, Fsample);
+        FDS_right_filt = ft_preproc_dftfilter(FDS_right_filt, Fsample);
+        
+        
+ %% take away the decay from MEPs
+        
+cfg = [];
+cfg.start_time = 0.005;
+cfg.end_time = 0.15;
+cfg.SR = SR;
+cfg.sp_trigger = sp_TMS;
+
+[ECD_sweep_cuts_tmp, ECD_sweep_time, ECD_sample_points] = gb_cutting_sweeps(cfg, ECD_right_filt);
+[FDS_sweep_cuts_tmp, FDS_sweep_time, FDS_sample_points] = gb_cutting_sweeps(cfg, FDS_right_filt);
+
+
+
+
+        %% get MEPs p2p amplitudes
+        cfg = [];
+        cfg.start_time = threshold{s}(1);
+        cfg.end_time = threshold{s}(2);
+        cfg.SR = SR;
+        cfg.sp_trigger = sp_TMS;
+        
+        [ECD_MEP_amplitudes] = gb_MEP_p2p(cfg, ECD_right_filt);
+        [FDS_MEP_amplitudes] = gb_MEP_p2p(cfg, FDS_right_filt);
+        
+        
+%% fit the decay and remove it
+
+cfg = [];
+cfg.start_time = -0.15;
+cfg.end_time = -0.002;
+cfg.SR = SR;
+cfg.sp_trigger = sp_TMS;
+
+[ECD_sweep_cuts_bkgrnd_tmp, ECD_sweep_time_bkgrnd, ECD_sample_points_bkgrnd] = gb_cutting_sweeps(cfg, ECD_right_filt)
+
+[FDS_sweep_cuts_bkgrnd_tmp, FDS_sweep_time_bkgrnd, FDS_sample_points_bkgrnd] = gb_cutting_sweeps(cfg, FDS_right_filt)
+
+%%
+for x = 1:length(ECD_sweep_cuts_bkgrnd_tmp)
+[~, ~, resECD{x}] = fit(ECD_sweep_time_bkgrnd',ECD_sweep_cuts_bkgrnd_tmp{x}','poly3'); % time in ms
+ECD_sweep_cuts_bkgrnd{x} =resECD{x}.residuals';
+
+[~, ~, resFDS{x}] = fit(FDS_sweep_time_bkgrnd',FDS_sweep_cuts_bkgrnd_tmp{x}','poly3'); % time in ms
+FDS_sweep_cuts_bkgrnd{x} =resFDS{x}.residuals';
+
+end
+clear resECD resFDS
+
+%%
+for x = 1:length(ECD_sample_points_bkgrnd)
+    
+ECD_right_filt(ECD_sample_points_bkgrnd{x}) = ECD_sweep_cuts_bkgrnd{x};
+FDS_right_filt(FDS_sample_points_bkgrnd{x}) = FDS_sweep_cuts_bkgrnd{x};
+
+end
+
+
+        %% any contraction befor TMS?
+        
+        cfg = [];
+        cfg.start_time = -0.05;
+        cfg.end_time = -0.002;
+        cfg.SR = SR;
+        cfg.sp_trigger = sp_TMS;
+        cfg.analysis = 'p2p';
+        
+        [ECD_MEP_bkgrnd] = gb_background_activity(cfg, ECD_right_filt)
+        [FDS_MEP_bkgrnd] = gb_background_activity(cfg, FDS_right_filt)
+        
+        
+        %% muscle contraction before TMS?
+        
+        MEP_excl_ECD_idx = find(ECD_MEP_bkgrnd>200) % microvolt
+        MEP_excl_FDS_idx = find(FDS_MEP_bkgrnd>200) % microvolt
+        
+        
+        %% EXCLUDE MEP<50 p2p microVolts
+        
+        MEP_excl2_ECD_idx = find(ECD_MEP_amplitudes<50) % microvolt
+       
+        
+         %% EXCLUDE out slow and out fast trials
+        
+      ECD_MEP_amplitudes(exclude) = nan;
+     
+        
+        %% exclude trials with background contraction
+        
+        ECD_MEP_amplitudes(union(MEP_excl_ECD_idx, MEP_excl_FDS_idx)) = nan;
+    
+        
+        ECD_MEP_amplitudes(MEP_excl2_ECD_idx) = nan;
+    
+ 
+        %% median amplitude based on color trigger         
+        summary.ECD_press{s,t} = nanmedian((ECD_MEP_amplitudes(CF_colors_order == 1)));
+        summary.ECD_lift{s,t} = nanmedian((ECD_MEP_amplitudes(CF_colors_order == -1)));
+        %%
+              
+cfg = [];
+cfg.start_time = -0.015;
+cfg.end_time = 0.15;
+cfg.SR = SR;
+cfg.sp_trigger = sp_TMS;
+
+[bbb.ECD_sweep_PLOT_CHECKS, bbb.ECD_sweep_time_PLOT_CHECKS, bbb.ECD_sample_points_PLOT_CHECKS] = gb_cutting_sweeps(cfg, ECD_right_filt);
+    
+        %%
+        
+        trial = (1:length(CF_colors_order))';
+        subject_rep = repmat(subjects(s), [length(CF_colors_order),1]);
+        TMS_timings_rep = repmat(TMS_timings(t), [length(CF_colors_order),1]);
+        
+        subject_table{t} = table(subject_rep,trial, ECD_MEP_amplitudes,...
+            ECD_MEP_bkgrnd,...
+            FDS_MEP_bkgrnd,...
+            TMS_timings_rep,...
+            CF_colors_order)
+        
+        
+        
+        
+        
+        
+        clearvars -except s t subjects TMS_timings threshold subject_table general_directory final_table summary
+        
+        %%
+       
+    end
+    
+    final_subject_table = vertcat(subject_table{:})
+    
+    
+     output_filename = strcat('final_subject_table_',subjects{s},'_' , '.mat');
+     save(output_filename,'final_subject_table')
+     
+     
+     final_table = vertcat(final_table, final_subject_table)
+     
+     clearvars -except s t subjects TMS_timings threshold general_directory final_table summary
+     
+     
+end
+
+
+%%
+cd(general_directory)
+filename = 'final_table_for_R_10-2500_def_o1_rep.xlsx';
+writetable(final_table,filename,'Sheet',1)
+

Pilot_analysis/gb_background_activity.m

@@ -0,0 +1,31 @@
+function [amplitudes, sweep_cut] = gb_background_activity(cfg, data)
+
+sweep_start = round(cfg.start_time*cfg.SR);
+sweep_end = round(cfg.end_time*cfg.SR);
+
+for x = 1:length(cfg.sp_trigger)
+   
+   sweep_points{x} = sweep_start+cfg.sp_trigger(x):cfg.sp_trigger(x)+sweep_end;
+    
+   sweep_cut{x} = data(sweep_points{x});
+   
+    %plot(sweep_cut{x})
+   switch  cfg.analysis
+       case 'p2p'
+   MEP_max{x} = max(sweep_cut{x});
+   MEP_min{x} = min(sweep_cut{x});
+   
+   amplitudes(x,1) = MEP_max{x} - MEP_min{x};
+   
+       case 'mean'
+         amplitudes(x,1) = mean(sweep_cut{x});  
+         
+       case 'median'
+         amplitudes(x,1) = median(sweep_cut{x});  
+         
+       case 'sum'
+         amplitudes(x,1) = sum(sweep_cut{x});   
+         
+   end
+
+end

Pilot_analysis/gb_cutting_sweeps.m

@@ -0,0 +1,38 @@
+function [sweep_cut, sweep_time,sweep_points] = gb_cutting_sweeps(cfg, data)
+
+sweep_start = round(cfg.start_time*cfg.SR);
+sweep_end = round(cfg.end_time*cfg.SR);
+
+
+for x = 1:length(cfg.sp_trigger)
+   
+   sweep_points{x} = sweep_start+cfg.sp_trigger(x):cfg.sp_trigger(x)+sweep_end;
+    
+   sweep_cut{x} = data(sweep_points{x});
+   
+   if isfield(cfg, 'baseline_rectify')
+       
+    % rectify epoch   
+    sweep_cut{x} =  abs(sweep_cut{x});
+    
+    
+      bs_start = round(cfg.baseline_rectify(1)*cfg.SR);
+    bs_end = round(cfg.baseline_rectify(2)*cfg.SR);
+    
+      bs_points{x} = bs_start+cfg.sp_trigger(x):cfg.sp_trigger(x)+bs_end;
+    
+   bs_cut{x} = abs(data(bs_points{x}));
+  
+    
+  
+    
+    % identify mean of the baselined period
+    baselined_chunk = mean(bs_cut{x},2);   
+    
+    % subtract the mean from the baselined period from the rectified epoch
+    sweep_cut{x} =  sweep_cut{x} - baselined_chunk;
+   end
+
+end
+
+ sweep_time = linspace(cfg.start_time,cfg.end_time, length(sweep_points{1})) ; 

Pilot_analysis/pilot_functions.rtf

@@ -0,0 +1,15 @@
+{\rtf1\ansi\ansicpg1252\cocoartf2636
+\cocoatextscaling0\cocoaplatform0{\fonttbl\f0\fswiss\fcharset0 Helvetica;}
+{\colortbl;\red255\green255\blue255;}
+{\*\expandedcolortbl;;}
+\paperw11900\paperh16840\margl1440\margr1440\vieww11520\viewh8400\viewkind0
+\pard\tx566\tx1133\tx1700\tx2267\tx2834\tx3401\tx3968\tx4535\tx5102\tx5669\tx6236\tx6803\pardirnatural\partightenfactor0
+
+\f0\fs24 \cf0 The folder contains the functions used to preprocess and analyse the pilot data\
+\
+JA_pilot_analysis_DEF_rep.m contains all the other functions\
+\
+JA_pilot_analysis.R was used to perform statistical analysis of pilot data\
+\
+\
+final_table_for_R_10-2500_def_o1_rep.xlsx contains the preprocessed data as a table which is used a s input fro the R script}