diff --git a/Analysis_Pipeline/AnalysisPipeline_Part1.m b/Analysis_Pipeline/AnalysisPipeline_Part1.m new file mode 100644 index 0000000..06e384b --- /dev/null +++ b/Analysis_Pipeline/AnalysisPipeline_Part1.m @@ -0,0 +1,50 @@ +%% Project Title: Functional Connectivity Fingerprints of Frontal Eye Field and Inferior Frontal Junction +%%% Update Date: 18-March-2022 + +clear +close all + +%% Part 1: Preprocessing +%% Compatible with megconnectome 3.0, fieldtrip-r10442 and MATLAB2012b + +%% Please Select! +subjectids = {fill this part}; + +trialDuration = 2; % 2, 5, or 10 seconds + + +%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% directories & paths +dir_subjects = load_path('dir_subjects'); + +addpath(genpath(load_path('megconnectome-30'))); +addpath(load_path('fieldtrip-r10442')); +ft_defaults +addpath(genpath(load_path('workingDir_code'))); + +%% +cimec = struct; +cimec.dir_subjects = dir_subjects; +cimec.scanids = {'3-Restin', '4-Restin', '5-Restin'}; +cimec.trialDuration = trialDuration; + +%% Step 1: clean_data + +f = waitbar(0, 'Looping around...'); + +for idx=1:length(subjectids) + + subjectid = subjectids{idx}; + cimec.subjectid = subjectid; + cimec.experimentid = [ subjectid, '_', 'MEG']; + cimec.pipelinedatadir = fullfile(dir_subjects, subjectid, 'MEG', 'Restin'); + + clean_data(cimec) + + waitbar(idx/length(subjectids), f, sprintf('Progress: %d %%', floor(idx/length(subjectids)*100))); + +end diff --git a/Analysis_Pipeline/AnalysisPipeline_Part2.m b/Analysis_Pipeline/AnalysisPipeline_Part2.m new file mode 100644 index 0000000..83a6309 --- /dev/null +++ b/Analysis_Pipeline/AnalysisPipeline_Part2.m @@ -0,0 +1,152 @@ +%% Project Title: Functional Connectivity Fingerprints of Frontal Eye Field and Inferior Frontal Junction +%%% Update Date: 18-March-2022 + +clear +close all + +%% Part 2: END-to-END script after Preprocessing +%% Compatible with brainstorm (Version: 28-May-2021 or later), fieldtrip-20210411 and MATLAB2020a + +%% Please Select! +% Also, please fill the load_path function inside helper_fun folder for correct directories + +seedRegions = {'L_FEF_ROI L', 'L_IFJa_ROI L', 'L_IFJp_ROI L', ... + 'R_FEF_ROI R', 'R_IFJa_ROI R', 'R_IFJp_ROI R'}; % seed-based connectivity + +trialDuration = 2; % 2, 5, or 10 seconds +connMetric = 'oPEC'; % iCOH, oPEC, dwPLI, or PDC + + +subjectids = {fill this part}; + + +%% fixed +step2 = true; %% Step 2: scouts_timeSeries +step3 = true; %% Step 3: functional_connectivity +step4 = true; %% Step 4: results_connectivity and group_average + +cleanProtocol = true; % meaning delete the protocol you created in Brainstorm +redefine = true; % to fix the floating numbers for cfg.foi +resample = true; fsample = 300; +inverse_measure = 'dspm2018'; % 'dspm2018' +top_k = 5; % save the names of highest functional couplings in Excel + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% paths +addpath(genpath(fullfile(pwd, 'helper_fun'))); +addpath(genpath(load_path('brainstorm'))); +addpath(load_path('fieldtrip-20210411')); +ft_defaults +addpath(genpath(load_path('workingDir_code'))); + +% directories +dir_subjects = load_path('dir_subjects'); +scoutsTimeSeries = load_path('scoutsTimeSeries'); +functionalConnectivity = fullfile(load_path('functionalConnectivity'), [connMetric '_' char(string(trialDuration)) 's']); +connectivityResults = fullfile(load_path('connectivityResults'), [connMetric '_' char(string(trialDuration)) 's']); +groupResults = fullfile(load_path('groupResults'), [connMetric '_' char(string(trialDuration)) 's']); + + +%% Step 2: scouts_timeSeries +if step2 == true + cimec = struct; + cimec.dir_subjects = dir_subjects; + cimec.trialDuration = trialDuration; + cimec.cleanProtocol = cleanProtocol; + cimec.fsample = fsample; + cimec.resample = resample; + cimec.inverse_measure = inverse_measure; + + % Loop over + f = waitbar(0, 'Looping around...'); + + for idx=1:length(subjectids) + + subjectid = subjectids{idx}; + cimec.subjectid = subjectid; + + scouts_timeSeries(cimec); + + waitbar(idx/length(subjectids), f, sprintf('Progress: %d %%', floor(idx/length(subjectids)*100))); + end + +end +%% Step 3: functional_connectivity +if step3 == true + cimec = struct; + cimec.scoutsTimeSeries = scoutsTimeSeries; + cimec.seedRegions = seedRegions; + cimec.trialDuration = trialDuration; + cimec.redefine = redefine; + cimec.conn_metric = connMetric; + + % Loop over + f = waitbar(0, 'Looping around...'); + + for idx=1:length(subjectids) + + subjectid = subjectids{idx}; + outputdir = fullfile(functionalConnectivity, subjectid); + mkdir(outputdir); + + cimec.subjectid = subjectid; + cimec.outputdir = outputdir; + + functional_connectivity(cimec); + + waitbar(idx/length(subjectids), f, sprintf('Progress: %d %%', floor(idx/length(subjectids)*100))); + end + +end +%% Step 4A: results_connectivity and group_average_topK_forSelectedROIs (for iCOH and dwPLI; Fieldtrip computes seed-based connectivity for these measures) +if step4 == true + + if or(strcmp(connMetric, 'iCOH'), strcmp(connMetric, 'dwPLI')) + cimec = struct; + cimec.outputdir = connectivityResults; + cimec.top_k = top_k; + cimec.trialDuration = trialDuration; + + % Loop over + f = waitbar(0, 'Looping around...'); + + for idx=1:length(subjectids) + + subjectid = subjectids{idx}; + inputdir = fullfile(functionalConnectivity, subjectid); + + cimec.subjectid = subjectid; + cimec.functionalConnectivity = inputdir; + + results_connectivity(cimec); + + waitbar(idx/length(subjectids), f, sprintf('Progress: %d %%', floor(idx/length(subjectids)*100))); + + end + + cimec = struct; + cimec.subjectids = subjectids; + cimec.top_k = top_k; + cimec.trialDuration = trialDuration; + cimec.connectivityResults = connectivityResults; + cimec.outputdir = groupResults; + + group_average_topK_forSelectedROIs(cimec); + +%% Step 4B: group_average_topK_wholeBrainConnectivity (oPEC and PDC; Fieldtrip computes the whole-brain connectivity map for these measures) + elseif or(strcmp(connMetric, 'oPEC'), strcmp(connMetric, 'PDC')) + cimec = struct; + cimec.subjectids = subjectids; + cimec.top_k = top_k; + cimec.trialDuration = trialDuration; + cimec.functionalConnectivity = functionalConnectivity; + cimec.outputdir = groupResults; + + group_average_topK_wholeConnectivity(cimec) + end + +end +%% diff --git a/Analysis_Pipeline/ReadMe.md b/Analysis_Pipeline/ReadMe.md new file mode 100644 index 0000000..8d15ec6 --- /dev/null +++ b/Analysis_Pipeline/ReadMe.md @@ -0,0 +1,14 @@ +## Analysis_Pipeline folder includes MATLAB code to implement functional connectivity analysis. + +AnalysisPipeline_Part1.m script preprocesses the MEG recordings and divides them into desired epoch lengths. It is compatible with megconnectome 3.0, fieldtrip-r10442 and MATLAB2012b. + +AnalysisPipeline_Part2.m is an END-to-END script to implement group-level functional connectivity analysis based on the preprocessed MEG data. It is compatible with brainstorm (Version: 28-May-2021), fieldtrip-20210411 and MATLAB2020a. + +Before running the code, please fill the directory information in the load_path.m script inside the helper_fun folder to run the code. + +### Ground-truth analyses + +oPEC: orthogonalized Power Envelope Correlation + +![3](https://user-images.githubusercontent.com/44211738/159386267-29a470da-98c1-4d02-bc8b-d00d20b22017.PNG) + diff --git a/Analysis_Pipeline/helper_fun/clean_data.m b/Analysis_Pipeline/helper_fun/clean_data.m new file mode 100644 index 0000000..75923a3 --- /dev/null +++ b/Analysis_Pipeline/helper_fun/clean_data.m @@ -0,0 +1,91 @@ +%% Project Title: Functional Connectivity Fingerprints of Frontal Eye Field and Inferior Frontal Junction +%%% This script for the preprocessing step. +%%% Update Date: 18-March-2022 + + +%% Compatible with megconnectome 3.0, fieldtrip-r10442 and MATLAB2012b + +%% This function was modified from the parts of megconnectome (3.0) to use my project at CIMEC, University of Trento. + +% megconnectome is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. + +% Copyright (C) 2011-2014 by the Human Connectome Project, WU-Minn Consortium (1U54MH091657) + +%% Output + +% computes the clean data. + + +%% Inputs +% %% Info about the data +% subjectid = cimec.subjectid; % e.g. '111514' +% experimentid = cimec.experimentid; % e.g. '111514_MEG' +% scanids = cimec.scanids; % e.g. {'3-Restin', '4-Restin', '5-Restin'} +% dir_subjects = cimec.dir_subjects; % e.g. 'C:\Users\ASUS\Desktop\4- Thesis\HCP\HCP_Subjects' +% pipelinedatadir = cimec.pipelinedatadir; % e.g. 'C:\Users\ASUS\Desktop\Restin' +% trialDuration = cimec.trialDuration; % e.g. 10 (seconds) + +%% Function + +function clean_data(cimec) + +% ensure that the time and date of execution are not stored in the provenance information +global ft_default +ft_default.trackcallinfo = 'no'; + +%% Info about the data +subjectid = cimec.subjectid; % e.g. '111514' +experimentid = cimec.experimentid; % e.g. '111514_MEG' +scanids = cimec.scanids; % e.g. {'3-Restin', '4-Restin', '5-Restin'} +dir_subjects = cimec.dir_subjects; % e.g. 'C:\Users\ASUS\Desktop\4- Thesis\HCP\HCP_Subjects' + + +pipelinedatadir = cimec.pipelinedatadir; % e.g. 'C:\Users\ASUS\Desktop\4- Thesis\HCP\HCP_Subjects\175237\MEG\Restin' +% change to the location of the processed data (input and output) +cd(pipelinedatadir) + +trialDuration = cimec.trialDuration; % e.g. 10 (seconds) + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% execute the pipeline +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +dir_baddata = fullfile(pipelinedatadir, 'baddata', '/'); +dir_icaclass = fullfile(pipelinedatadir, 'icaclass', '/'); +dir_cleandata = fullfile(pipelinedatadir, 'cleandata', '/'); +mkdir(dir_cleandata); + + +for idx=1:length(scanids) + + raw_data = fullfile(dir_subjects, subjectid, 'unprocessed', 'MEG', scanids{idx}, '4D', 'c,rfDC'); + + + resultprefix = sprintf('%s_%s', experimentid, scanids{idx}); + badchansuffix = 'baddata_badchannels'; % mnemonic for file that contains bad channel and segment info + badsegmsuffix = 'baddata_badsegments'; % mnemonic for file that contains bad channel and segment info + icainfosuffix = 'icaclass_vs'; % mnemonic to save results + + cfg = []; + cfg.dataset = raw_data; + cfg.trialfun = 'trialfun_Restin'; + cfg.trialdef.trialDuration = trialDuration; + cfg = ft_definetrial(cfg); + cfg.trl(:,4) = nan; % hcp_extract_allfromrun assumes a trigger code to be present + + cfg.badchanfile = [dir_baddata, resultprefix, '_', badchansuffix, '.txt']; + cfg.badsegmfile = [dir_baddata, resultprefix, '_', badsegmsuffix, '.txt']; + cfg.icainfofile = [dir_icaclass, resultprefix, '_', icainfosuffix]; + cfg.montage = hcp_exgmontage(subjectid, experimentid, scanids{idx}); + cfg.lineFreq = [60 120]; + cfg.badsegmode = 'remfull'; + cfg.outputfile = [dir_cleandata, resultprefix, '_cleanData_', sprintf('%ds', trialDuration)]; + + hcp_extract_allfromrun(cfg); + +end +end \ No newline at end of file diff --git a/Analysis_Pipeline/helper_fun/functional_connectivity.m b/Analysis_Pipeline/helper_fun/functional_connectivity.m new file mode 100644 index 0000000..4269bd0 --- /dev/null +++ b/Analysis_Pipeline/helper_fun/functional_connectivity.m @@ -0,0 +1,149 @@ +%% Project Title: Functional Connectivity Fingerprints of Frontal Eye Field and Inferior Frontal Junction +%%% Functional connectivity results per subject. +%%% Update Date: 18-March-2022 + +%% Output + +% computes the connectivity matrices per subject. + + +%% Inputs +% subjectid = cimec.subjectid; % e.g. '111514' +% seedRegions = cimec.seedRegions; % e.g. {'L_FEF_ROI L', 'L_IFJa_ROI L', 'L_IFJp_ROI L', 'R_FEF_ROI R', 'R_IFJa_ROI R', 'R_IFJp_ROI R'}; +% trialDuration = cimec.trialDuration; % e.g. 10 (seconds) +% +% redefine = cimec.redefine; % e.g. true +% +% conn_metric = cimec.conn_metric; % e.g. iCOH; oPEC; dwPLI +% +% scoutsTimeSeries = cimec.scoutsTimeSeries; % e.g. '...\HCP_Subjects\Outputs\scouts_timeSeries\' +% +% outputdir = cimec.outputdir; % e.g. '...\HCP_Subjects\Outputs\connectivity_matrices\'; + +%% Function + +function functional_connectivity(cimec) + +%% Info about the data +subjectid = cimec.subjectid; % e.g. '111514' +seedRegions = cimec.seedRegions; % e.g. {'L_FEF_ROI L', 'L_IFJa_ROI L', 'L_IFJp_ROI L', 'R_FEF_ROI R', 'R_IFJa_ROI R', 'R_IFJp_ROI R'}; +trialDuration = cimec.trialDuration; % e.g. 10 (seconds) + +redefine = cimec.redefine; % e.g. true + +scoutsTimeSeries = cimec.scoutsTimeSeries; % e.g. '...\HCP_Subjects\Outputs\scouts_timeSeries\' +inputfile = fullfile(scoutsTimeSeries, [subjectid, '_scoutsTimeSeries_', sprintf('%ds.mat', trialDuration)]); + +conn_metric = cimec.conn_metric; + +% Output file +outputdir = cimec.outputdir; % e.g. '...\HCP_Subjects\Outputs\functional_connectivity\'; + +%% Frequencies of interest +% frequency bands +band_names = { 'delta' , 'theta' , 'alpha' , 'beta' , 'gamma' }; +freq_bands = [ 1 4 ; 4 8 ; 8 13 ; 13 30 ; 30 100 ]; + +N_bands = length(band_names); + +% step +foi_step = [ 0.5 ; 0.5 ; 0.5 ; 1 ; 2 ]; +% smoothing +tapsmofrq = [ 0.5 ; 0.5 ; 0.5 ; 1 ; 5 ]; + +%% Pipeline +%% Load FT_data (scoutsTimeSeries) +load(inputfile, 'FT_data') + +% selected scouts +scouts = {}; + +for r = 1:length(seedRegions) + for i = 1:length(FT_data.label) + scouts{end+1, 1} = seedRegions{r}; + scouts{end, 2} = FT_data.label{i}; + end +end + +%% Loop over frequency bands +if redefine == true + cfg = []; + cfg.length = trialDuration; + cfg.overlap = 0; + data = ft_redefinetrial(cfg, FT_data); +else + data = FT_data; +end + +for i = 1:N_bands + + %% power spectral density (PSD) + cfg = []; + %cfg.foilim = [freq_bands(i,1) freq_bands(i,end)]; + cfg.foi = [freq_bands(i,1) : foi_step(i) : freq_bands(i,end)]; + cfg.method = 'mtmfft'; + cfg.output = 'fourier'; + cfg.tapsmofrq = tapsmofrq(i); + freq = ft_freqanalysis(cfg, data); + + % dir_freq = fullfile(outputdir, '_freqanalysis'); + % mkdir(dir_freq); + % save(fullfile(dir_freq, [subjectid, '_freqanalysis_', band_names{i}, sprintf('_%ds.mat', trialDuration)]), 'freq', '-v7.3'); + + + %% connectivity analysis + % https://www.fieldtriptoolbox.org/reference/ft_connectivityanalysis/ + + if strcmp(conn_metric, 'iCOH') + % imaginary part of coherency + cfg = []; + cfg.method = 'coh'; + cfg.complex = 'absimag'; + cfg.channelcmb = scouts; + conn = ft_connectivityanalysis(cfg, freq); + + elseif strcmp(conn_metric, 'oPEC') + % orthogonalized power envelope correlation + cfg = []; + cfg.method = 'powcorr_ortho'; + conn = ft_connectivityanalysis(cfg, freq); + + elseif strcmp(conn_metric, 'PDC') + % partial directed coherence + cfg = []; + cfg.method = 'pdc'; + conn = ft_connectivityanalysis(cfg, freq); + + elseif strcmp(conn_metric, 'dwPLI') + % https://mailman.science.ru.nl/pipermail/fieldtrip/2011-June/003892.html + % https://mailman.science.ru.nl/pipermail/fieldtrip/2019-June/039226.html + % https://mailman.science.ru.nl/pipermail/fieldtrip/2016-February/022995.html + % https://mailman.science.ru.nl/pipermail/fieldtrip/2014-January/020309.html + + % debiased weighted phase lag index + cfg = []; + cfg.method = 'wpli_debiased'; + cfg.channelcmb = scouts; + conn = ft_connectivityanalysis(cfg, freq); + conn.wpli_debiasedspctrm(conn.wpli_debiasedspctrm<0) = 0; + + end + + clear freq; + + %% average over frequencies + cfg = []; + cfg.frequency = 'all'; + cfg.avgoverfreq = 'yes'; + avgFreq_funConn = ft_selectdata(cfg, conn); + + clear conn; + + %% save funConn + avgFreq_funConn.label = FT_data.label; + avgFreq_funConn.conn_metric = conn_metric; + save(fullfile(outputdir, [subjectid, '_functionalConnectivity_', band_names{i}, sprintf('_%ds.mat', trialDuration)]), 'avgFreq_funConn', '-v7.3'); + clear avgFreq_funConn +end + +end diff --git a/Analysis_Pipeline/helper_fun/group_average_topK_forSelectedROIs.m b/Analysis_Pipeline/helper_fun/group_average_topK_forSelectedROIs.m new file mode 100644 index 0000000..1713b12 --- /dev/null +++ b/Analysis_Pipeline/helper_fun/group_average_topK_forSelectedROIs.m @@ -0,0 +1,127 @@ +%% Project Title: Functional Connectivity Fingerprints of Frontal Eye Field and Inferior Frontal Junction +%%% Group average for iCOH and dwPLI metrics. +%%% Update Date: 18-March-2022 + + +%% Compatible with fieldtrip-20210411 and MATLAB2020a + +function group_average_topK_forSelectedROIs(cimec) + +%% Info about the data +subjectids = cimec.subjectids; % e.g. subjectids = {'111514', ...}; +top_k = cimec.top_k; % e.g. 5 +trialDuration = cimec.trialDuration; % e.g. 10 (seconds) + +connectivityResults = cimec.connectivityResults; %e.g. '...\HCP_Subjects\Outputs\connectivity_results\'; + +% Output file +outputdir = cimec.outputdir; % '...\HCP_Subjects\Outputs\group_results\'; + +%% Frequency bands +band_names = { 'delta' , 'theta' , 'alpha' , 'beta' , 'gamma' }; +freq_bands = [ 1 4 ; 4 8 ; 8 13 ; 13 30 ; 30 100 ]; + +N_bands = length(band_names); + +%% Loop over frequency bands and subjects +group_results = {}; + +for band = 1:N_bands + + group_results(band).band = band_names{band}; + + for s = 1:length(subjectids) + load(fullfile(connectivityResults, [subjectids{s}, '_connectivityResults_', sprintf('%ds.mat', trialDuration)]), ... + 'conn_results'); + + tmp_data = conn_results(band).data; + group_results(band).data(:,:,s) = tmp_data; + + clear tmp_data; + end + + group_results(band).name = conn_results(band).name; + group_results(band).pair = conn_results(band).pair; + group_results(band).band_names = conn_results(band).band_names; + group_results(band).freq_bands = conn_results(band).freq_bands; + group_results(band).label = conn_results(band).label; + % mean + group_results(band).data_perSubject = group_results(band).data; + group_results(band).data = mean(group_results(band).data, 3); + + + %% Top K regions that has high connectivity with selected RIOs + N = length(group_results(band).name); + + group_results(band).top_k.name = cell(top_k, N); + group_results(band).top_k.data = zeros(top_k, N); + + for i = 1:N + + values = group_results(band).data(:, i); + [max_top, idx_top] = maxk(values, top_k); + group_results(band).top_k.data(:,i) = max_top; + group_results(band).top_k.name(:,i) = group_results(band).pair(idx_top, i); + + end + + group_results(band).subjectids = subjectids; +end + +%% Save the results + +% .mat file +id = char(string(length(subjectids))); +mkdir(outputdir); +save(fullfile(outputdir, [id, 'Subjects', '_groupResults_', sprintf('%ds.mat', trialDuration)]), 'group_results', '-v7.3'); + + +%% Excel file + +% results +for idx_band = 1:length(group_results(1).band_names) + + band_name = group_results(idx_band).band; + results(idx_band).band_name = band_name; + + for idx_RIO = 1:size(group_results(idx_band).name, 1) + + RIOs_name = group_results(idx_band).name{idx_RIO}; + top_k = group_results(idx_band).top_k.name(:, idx_RIO); + data = group_results(idx_band).top_k.data(:, idx_RIO); + + results(idx_band).top_k{idx_RIO} = top_k; + results(idx_band).data{idx_RIO} = data; + results(idx_band).RIOs_name{idx_RIO} = RIOs_name; + end +end + +% make a table +excel_name = fullfile(outputdir, [id, 'Subjects', '_topK_', sprintf('%ds.xlsx', trialDuration)]); + +for idx_band = 1:length(results) + + tmp_band = results(idx_band); + + for idx_RIO = 1:length(tmp_band.RIOs_name) + + Pair = [tmp_band.top_k{idx_RIO}; {'-------------'}]; + Icoh = [tmp_band.data{idx_RIO}; NaN]; + RIO = [repmat(tmp_band.RIOs_name(idx_RIO), length(Pair)-1, 1); {'-------------'}]; + + if idx_RIO == 1 + make_table = table(RIO, Pair, Icoh); + else + T = table(RIO, Pair, Icoh); + make_table = [make_table; T]; + end + end + + % save as an excel file + writetable(make_table, excel_name, 'Sheet', results(idx_band).band_name, 'Range', 'C3:E200'); + +end + +disp('The results are saved.') + +end \ No newline at end of file diff --git a/Analysis_Pipeline/helper_fun/group_average_topK_wholeConnectivity.m b/Analysis_Pipeline/helper_fun/group_average_topK_wholeConnectivity.m new file mode 100644 index 0000000..e56a43c --- /dev/null +++ b/Analysis_Pipeline/helper_fun/group_average_topK_wholeConnectivity.m @@ -0,0 +1,72 @@ +%% Project Title: Functional Connectivity Fingerprints of Frontal Eye Field and Inferior Frontal Junction +%%% Group average for oPEC and PDC metrics. +%%% Update Date: 18-March-2022 + + +%% Compatible with fieldtrip-20210411 and MATLAB2020a + +function group_average_topK_wholeConnectivity(cimec) + +%% Info about the data +subjectids = cimec.subjectids; % e.g. subjectids = {'111514', '140117', '162026', '164636', '191033'}; +top_k = cimec.top_k; % e.g. 5 +trialDuration = cimec.trialDuration; % e.g. 10 (seconds) + +functionalConnectivity = cimec.functionalConnectivity; %e.g. '...\HCP_Subjects\Outputs\functional_connectivity\'; + +% Output file +outputdir = cimec.outputdir; % '...\HCP_Subjects\Outputs\group_results\'; + +%% Frequency bands +band_names = { 'delta' , 'theta' , 'alpha' , 'beta' , 'gamma' }; +freq_bands = [ 1 4 ; 4 8 ; 8 13 ; 13 30 ; 30 100 ]; + +N_bands = length(band_names); + +%% Loop over frequency bands and subjects +group_results = {}; + +for band = 1:N_bands + + group_results(band).band = band_names{band}; + + for s = 1:length(subjectids) + load(fullfile(functionalConnectivity, subjectids{s}, [subjectids{s}, '_functionalConnectivity_', band_names{band}, sprintf('_%ds.mat', trialDuration)])); + conn_metric = avgFreq_funConn.conn_metric; + + if strcmp(conn_metric, 'oPEC') + % orthogonalized power envelope correlation + tmp_data = avgFreq_funConn.powcorrspctrm; + + elseif strcmp(conn_metric, 'PDC') + % partial directed coherence + tmp_data = avgFreq_funConn.pdcspctrm; + end + + group_results(band).data(:,:,s) = tmp_data; + + clear tmp_data; + end + + group_results(band).band_names = band_names; + group_results(band).freq_bands = freq_bands; + group_results(band).label = avgFreq_funConn.label; + % mean + group_results(band).data_perSubject = group_results(band).data; + group_results(band).data = mean(group_results(band).data, 3); + + + group_results(band).subjectids = subjectids; +end + +%% Save the results + +% .mat file +id = char(string(length(subjectids))); +mkdir(outputdir); +save(fullfile(outputdir, [id, 'Subjects', '_groupResults_', sprintf('%ds.mat', trialDuration)]), 'group_results', '-v7.3'); + + +disp('The results are saved.') + +end \ No newline at end of file diff --git a/Analysis_Pipeline/helper_fun/load_path.m b/Analysis_Pipeline/helper_fun/load_path.m new file mode 100644 index 0000000..57a7518 --- /dev/null +++ b/Analysis_Pipeline/helper_fun/load_path.m @@ -0,0 +1,35 @@ +%% Project Title: Functional Connectivity Fingerprints of Frontal Eye Field and Inferior Frontal Junction +%%% Loading the paths. +%%% Update Date: 18-March-2022 + +function path = load_path(directory) + +% libraries and code +if strcmp(directory, 'brainstorm') + path = '/home/orhan.soyuhos/Desktop/Orhan/MATLAB_lib/brainstorm3'; +elseif strcmp(directory, 'fieldtrip-r10442') + path = '/home/orhan.soyuhos/Desktop/Orhan/MATLAB_lib/fieldtrip-r10442'; +elseif strcmp(directory, 'fieldtrip-20210411') + path = '/home/orhan.soyuhos/Desktop/Orhan/MATLAB_lib/fieldtrip-20210411'; + +elseif strcmp(directory, 'megconnectome-30') + path = '/home/orhan.soyuhos/Desktop/Orhan/MATLAB_lib/megconnectome-3.0'; +elseif strcmp(directory, 'workingDir_code') + path = '/home/orhan.soyuhos/Desktop/Thesis_OrhanSoyuhos/Analysis_Pipeline'; + +% raw data +elseif strcmp(directory, 'dir_subjects') + path = '/home/orhan.soyuhos/Desktop/Orhan/HCP_Subjects'; + +% middle and final outputs +elseif strcmp(directory, 'scoutsTimeSeries') + path = '/home/orhan.soyuhos/Desktop/Orhan/HCP_Subjects/Outputs/scouts_timeSeries/'; +elseif strcmp(directory, 'functionalConnectivity') + path = '/home/orhan.soyuhos/Desktop/Orhan/HCP_Subjects/Outputs/functional_connectivity/'; +elseif strcmp(directory, 'connectivityResults') + path = '/home/orhan.soyuhos/Desktop/Orhan/HCP_Subjects/Outputs/connectivity_results/'; +elseif strcmp(directory, 'groupResults') + path = '/home/orhan.soyuhos/Desktop/Orhan/HCP_Subjects/Outputs/group_results/'; + + +end diff --git a/Analysis_Pipeline/helper_fun/results_connectivity.m b/Analysis_Pipeline/helper_fun/results_connectivity.m new file mode 100644 index 0000000..5ef8ff3 --- /dev/null +++ b/Analysis_Pipeline/helper_fun/results_connectivity.m @@ -0,0 +1,159 @@ +%% Project Title: Functional Connectivity Fingerprints of Frontal Eye Field and Inferior Frontal Junction +%%% Skip this part for oPEC and PDC connectivity measures! +%%% Update Date: 18-March-2022 + + +%% Compatible with fieldtrip-20210411 and MATLAB2020a + +%% Output + +% outputs per subject connectivity matrices for group level analysis +% AND finds the top-K scouts that have the best connectivity measures with the selected RIOs per subject + +%% Inputs +% subjectid = cimec.subjectid; % e.g. '111514' +% top_k = cimec.top_k; % e.g. 5 +% trialDuration = cimec.trialDuration; % e.g. 10 (seconds) +% +% functionalConnectivity = cimec.functionalConnectivity; %e.g. '...\HCP_Subjects\Outputs\functional_connectivity\' +% +% % Output file +% outputdir = cimec.outputdir; % '...\HCP_Subjects\Outputs\connectivity_results\'; + +%% Function + +function results_connectivity(cimec) + +%% Info about the data +subjectid = cimec.subjectid; % e.g. '111514' +top_k = cimec.top_k; % e.g. 5 +trialDuration = cimec.trialDuration; % e.g. 10 (seconds) + +functionalConnectivity = cimec.functionalConnectivity; %e.g. '...\HCP_Subjects\Outputs\functional_connectivity\' + +% Output file +outputdir = cimec.outputdir; % '...\HCP_Subjects\Outputs\connectivity_results\'; + +%% Frequency bands +band_names = { 'delta' , 'theta' , 'alpha' , 'beta' , 'gamma' }; +freq_bands = [ 1 4 ; 4 8 ; 8 13 ; 13 30 ; 30 100 ]; + +N_bands = length(band_names); + +%% Loop over frequency bands +conn_results = {}; + +for band = 1:N_bands + + %% Load frequency band + inputfile = fullfile(functionalConnectivity, [subjectid, '_functionalConnectivity_', band_names{band}, sprintf('_%ds.mat', trialDuration)]); + load(inputfile, 'avgFreq_funConn') + + conn_metric = avgFreq_funConn.conn_metric; + + %% Organize the pairs of scouts' connectivity value + conn_results(band).band = band_names{band}; + conn_results(band).name = unique(avgFreq_funConn.labelcmb(:,1)); + + N = length(conn_results(band).name); + max_k = size(avgFreq_funConn.labelcmb, 1); + step = max_k/N; + k = [0:step:max_k]; + + conn_results(band).data = zeros(step, N); + conn_results(band).pair = cell(step, N); + + for i = 1:N + + if strcmp(conn_metric, 'iCOH') + % for imaginary part of coherency + conn_results(band).data(:,i) = [avgFreq_funConn.cohspctrm(k(i)+1:k(i)+step)]'; + + elseif strcmp(conn_metric, 'dwPLI') + % https://mailman.science.ru.nl/pipermail/fieldtrip/2011-June/003892.html + % https://mailman.science.ru.nl/pipermail/fieldtrip/2019-June/039226.html + % https://mailman.science.ru.nl/pipermail/fieldtrip/2016-February/022995.html + % https://mailman.science.ru.nl/pipermail/fieldtrip/2014-January/020309.html + + % debiased weighted phase lag index + conn_results(band).data(:,i) = [avgFreq_funConn.wpli_debiasedspctrm(k(i)+1:k(i)+step)]'; + + + end + conn_results(band).pair(:,i) = [avgFreq_funConn.labelcmb(k(i)+1:k(i)+step, 2)]'; + + end + + %% Top K regions that has high connectivity with selected RIOs + conn_results(band).top_k.name = cell(top_k, N); + conn_results(band).top_k.data = zeros(top_k, N); + + for i = 1:N + + values = conn_results(band).data(:, i); + [max_top, idx_top] = maxk(values, top_k); + conn_results(band).top_k.data(:,i) = max_top; + conn_results(band).top_k.name(:,i) = conn_results(band).pair(idx_top, i); + + end + + conn_results(band).band_names = band_names; + conn_results(band).freq_bands = freq_bands; + conn_results(band).label = avgFreq_funConn.label; +end + +%% Save the results + +% .mat file +mkdir(outputdir); +save(fullfile(outputdir, [subjectid, '_connectivityResults_', sprintf('%ds.mat', trialDuration)]), 'conn_results', '-v7.3'); + +%% Excel file + +% results +for idx_band = 1:length(conn_results(1).band_names) + + band_name = conn_results(idx_band).band; + results(idx_band).band_name = band_name; + + for idx_RIO = 1:size(conn_results(idx_band).name, 1) + + RIOs_name = conn_results(idx_band).name{idx_RIO}; + top_k = conn_results(idx_band).top_k.name(:, idx_RIO); + data = conn_results(idx_band).top_k.data(:, idx_RIO); + + results(idx_band).top_k{idx_RIO} = top_k; + results(idx_band).data{idx_RIO} = data; + results(idx_band).RIOs_name{idx_RIO} = RIOs_name; + end +end + +% make a table +excel_name = fullfile(outputdir, [subjectid, '_topK_', sprintf('%ds.xlsx', trialDuration)]); + +for idx_band = 1:length(results) + + tmp_band = results(idx_band); + + for idx_RIO = 1:length(tmp_band.RIOs_name) + + Pair = [tmp_band.top_k{idx_RIO}; {'-------------'}]; + Icoh = [tmp_band.data{idx_RIO}; NaN]; + RIO = [repmat(tmp_band.RIOs_name(idx_RIO), length(Pair)-1, 1); {'-------------'}]; + + if idx_RIO == 1 + make_table = table(RIO, Pair, Icoh); + else + T = table(RIO, Pair, Icoh); + make_table = [make_table; T]; + end + end + + % save as an excel file + writetable(make_table, excel_name, 'Sheet', results(idx_band).band_name, 'Range', 'C3:E40'); + +end + +disp('The results are saved.') + +end \ No newline at end of file diff --git a/Analysis_Pipeline/helper_fun/scouts_timeSeries.m b/Analysis_Pipeline/helper_fun/scouts_timeSeries.m new file mode 100644 index 0000000..2767a1d --- /dev/null +++ b/Analysis_Pipeline/helper_fun/scouts_timeSeries.m @@ -0,0 +1,466 @@ +%% Project Title: Functional Connectivity Fingerprints of Frontal Eye Field and Inferior Frontal Junction +%%% This script is for source reconstruction. +%%% Update Date: 18-March-2022 + + +%% Compatible with brainstorm (Version: 28-May-2021 or later) and MATLAB2020a + +%% Output + +% computes the scout time series per subject. + +%% Inputs +% % Select the directory where the subject folder is +% dir_subjects = cimec.dir_subjects % 'C:\Users\ASUS\Desktop\4- Thesis\HCP\HCP_Subjects'; +% % Select the subject +% subjectid = cimec.subjectid % '111514'; +% trialDuration = cimec.trialDuration % 6 +% cleanProtocol = cimec.cleanProtocol % false +% fsample = cimec.fsample % 300 +% resample = cimec.resample; % true +% inverse_measure = cimec.inverse_measure; % 'dspm2018' + +%% Function + +function scouts_timeSeries(cimec) + +tic +%% ===== SUBJECT INFORMATION (*Be Sure Before Running the Script) ===== +% Select the directory where the subject folder is +dir_subjects = cimec.dir_subjects; % 'C:\Users\ASUS\Desktop\4- Thesis\HCP\HCP_Subjects' +% Select the subject +subjectid = cimec.subjectid; % '111514' + +trialDuration = cimec.trialDuration; % 6 +cleanProtocol = cimec.cleanProtocol; % false + +fsample = cimec.fsample; % 300 +resample = cimec.resample; % true + +inverse_measure = cimec.inverse_measure; % 'dspm2018' + +%% ===== FILES TO IMPORT ===== + +% You have to specify the folder in which the tutorial dataset is unzipped +if isempty(dir_subjects) || ~file_exist(dir_subjects) + error('The first argument must be the full path to the tutorial dataset folder.'); +end + +% Build the path of the files to import +extended_AnatDir = fullfile(dir_subjects, subjectid, 'T1w', subjectid); +HCP_anatomy = fullfile(dir_subjects, subjectid, 'MEG', 'anatomy'); + +Run1File = fullfile(dir_subjects, subjectid, 'MEG', 'Restin', 'cleandata', ... + [subjectid, '_MEG_3-Restin_cleanData_', sprintf('%ds.mat', trialDuration)]); +Run2File = fullfile(dir_subjects, subjectid, 'MEG', 'Restin', 'cleandata', ... + [subjectid, '_MEG_4-Restin_cleanData_', sprintf('%ds.mat', trialDuration)]); +Run3File = fullfile(dir_subjects, subjectid, 'MEG', 'Restin', 'cleandata', ... + [subjectid, '_MEG_5-Restin_cleanData_', sprintf('%ds.mat', trialDuration)]); + +NoiseFile = fullfile(dir_subjects, subjectid, 'unprocessed', 'MEG', '1-Rnoise', '4D', 'c,rfDC'); + +% Output path +OutputDir_scouts = fullfile(dir_subjects, 'Outputs', 'scouts_timeSeries'); +OutputDir_fsaverage = fullfile(dir_subjects, 'Outputs', 'fsaverage_anatomy'); +mkdir(OutputDir_scouts); +mkdir(OutputDir_fsaverage); + +% Check if the folder contains the required files +if ~file_exist(extended_AnatDir) || ~file_exist(Run1File) || ~file_exist(NoiseFile) + error(['The folder ' dir_subjects ' does not contain the selected subject from the HCP-MEG distribution.']); +end + +%% ===== CREATE PROTOCOL ===== + +% The protocol name has to be a valid folder name (no spaces, no weird characters...) +ProtocolName = ['Thesis_OrhanSoyuhos_', date]; + +% Start brainstorm without the GUI +if ~brainstorm('status') + brainstorm nogui +end + +% Delete existing protocol +gui_brainstorm('DeleteProtocol', ProtocolName); +% Create new protocol +gui_brainstorm('CreateProtocol', ProtocolName, 0, 0); +% Start a new report +bst_report('Start'); +% Reset colormaps +bst_colormaps('RestoreDefaults', 'meg'); + +%% ===== GET THE FIDUCIAL POINTS FROM THE old_anatomy ===== + +% Process: Import the old_anatomy folder to get fiducials +bst_process('CallProcess', 'process_import_anatomy', [], [], ... + 'subjectname', 'old_anatomy', ... + 'mrifile', {HCP_anatomy, 'HCPv3'}, ... + 'nvertices', 15000); + +ProtocolInfo = bst_get('ProtocolInfo'); +dirFolder = bst_fullfile(ProtocolInfo.SUBJECTS, 'old_anatomy/subjectimage_T1w_acpc_dc_restore.mat'); +sMri_oldAnatomy = load(dirFolder); + +% MRI coordinates of fiducials (for old_anatomy) +NAS = sMri_oldAnatomy.SCS.NAS; +LPA = sMri_oldAnatomy.SCS.LPA; +RPA = sMri_oldAnatomy.SCS.RPA; +AC = sMri_oldAnatomy.NCS.AC; +PC = sMri_oldAnatomy.NCS.PC; +IH = sMri_oldAnatomy.NCS.IH; + +% Convert the MRI coordinates to World coordinates (for old_anatomy) +nas = cs_convert(sMri_oldAnatomy, 'mri', 'world', NAS ./ 1000) .* 1000; +lpa = cs_convert(sMri_oldAnatomy, 'mri', 'world', LPA ./ 1000) .* 1000; +rpa = cs_convert(sMri_oldAnatomy, 'mri', 'world', RPA ./ 1000) .* 1000; +ac = cs_convert(sMri_oldAnatomy, 'mri', 'world', AC ./ 1000) .* 1000; +pc = cs_convert(sMri_oldAnatomy, 'mri', 'world', PC ./ 1000) .* 1000; +ih = cs_convert(sMri_oldAnatomy, 'mri', 'world', IH ./ 1000) .* 1000; + +%% ===== IMPORT THE DATA FREESURFER ANATOMY (new_anatomy) ===== + +% find the MRI file to import +MriFile = file_find(extended_AnatDir, 'T1.mgz', [], 1); +% Process: Import MRI (Because we need sMri for the Freesurfer anatomy (new_anatomy)) +bst_process('CallProcess', 'process_import_mri', [], [], ... + 'subjectname', 'Tmp_Anat', ... + 'mrifile', {MriFile, 'MGH'}, ... + 'nas', [0, 0, 0], ... + 'lpa', [0, 0, 0], ... + 'rpa', [0, 0, 0], ... + 'ac', [0, 0, 0], ... + 'pc', [0, 0, 0], ... + 'ih', [0, 0, 0]); + +dirFolder = bst_fullfile(ProtocolInfo.SUBJECTS, 'Tmp_Anat/subjectimage_T1.mat'); +sMri_newAnatomy = load(dirFolder); + +% Convert the World coordinates of fidual points from old_anatomy to the MRI coordinates of new_anatomy +final_NAS = cs_convert(sMri_newAnatomy, 'world', 'mri', nas ./ 1000) .* 1000; +final_LPA = cs_convert(sMri_newAnatomy, 'world', 'mri', lpa ./ 1000) .* 1000; +final_RPA = cs_convert(sMri_newAnatomy, 'world', 'mri', rpa ./ 1000) .* 1000; +final_AC = cs_convert(sMri_newAnatomy, 'world', 'mri', ac ./ 1000) .* 1000; +final_PC = cs_convert(sMri_newAnatomy, 'world', 'mri', pc ./ 1000) .* 1000; +final_IH =cs_convert(sMri_newAnatomy, 'world', 'mri', ih ./ 1000) .* 1000; + +% Process: Import FreeSurfer folder +bst_process('CallProcess', 'process_import_anatomy', [], [], ... + 'subjectname', subjectid, ... + 'mrifile', {extended_AnatDir, 'FreeSurfer'}, ... + 'nvertices', 15000, ... + 'nas', final_NAS, ... + 'lpa', final_LPA, ... + 'rpa', final_RPA, ... + 'ac', final_AC, ... + 'pc', final_PC, ... + 'ih', final_IH); + + +% Get all the subjects in the protocol +sProtocolSubjects = bst_get('ProtocolSubjects'); +% Choose the subject +del1_iSubject = find(strcmpi('Tmp_Anat', {sProtocolSubjects.Subject.Name})); +del2_iSubject = find(strcmpi('old_anatomy', {sProtocolSubjects.Subject.Name})); +% Delete the subject +db_delete_subjects(del1_iSubject); +db_delete_subjects(del2_iSubject); +% Reload the Data +db_reload_database('current'); + +%% ===== IMPORT the CLEAN and EPOCHED DATA ===== + +% Process: Import MEG/EEG: Existing epochs +sFileRun1 = bst_process('CallProcess', 'process_import_data_epoch', [], [], ... + 'subjectname', subjectid, ... + 'condition', sprintf('3-Restin_cleandata_%ds', trialDuration), ... + 'datafile', {Run1File, 'FT-TIMELOCK'}, ... + 'iepochs', [], ... + 'eventtypes', '', ... + 'createcond', 0, ... + 'channelalign', 1, ... + 'usectfcomp', 1, ... + 'usessp', 1, ... + 'freq', [], ... + 'baseline', []); +% Process: DC offset correction: [All file] +sFileRun1 = bst_process('CallProcess', 'process_baseline', sFileRun1, [], ... + 'baseline', [], ... + 'sensortypes', 'MEG, EEG', ... + 'method', 'bl', ... % DC offset correction: x_std = x - μ + 'overwrite', 1); +if resample == true + % Process: Resample: 300Hz + sFileRun1 = bst_process('CallProcess', 'process_resample', sFileRun1, [], ... + 'freq', fsample, ... + 'overwrite', 1); +else + dirFolder = bst_fullfile(ProtocolInfo.STUDIES, sFileRun1(1).FileName); + tmp_data = load(dirFolder); + oldRate = 1 ./ (tmp_data.Time(2)-tmp_data.Time(1)); + fsample = oldRate; +end + +% Process: Import MEG/EEG: Existing epochs +sFileRun2 = bst_process('CallProcess', 'process_import_data_epoch', [], [], ... + 'subjectname', subjectid, ... + 'condition', sprintf('4-Restin_cleandata_%ds', trialDuration), ... + 'datafile', {Run2File, 'FT-TIMELOCK'}, ... + 'iepochs', [], ... + 'eventtypes', '', ... + 'createcond', 0, ... + 'channelalign', 1, ... + 'usectfcomp', 1, ... + 'usessp', 1, ... + 'freq', [], ... + 'baseline', []); +% Process: DC offset correction: [All file] +sFileRun2 = bst_process('CallProcess', 'process_baseline', sFileRun2, [], ... + 'baseline', [], ... + 'sensortypes', 'MEG, EEG', ... + 'method', 'bl', ... % DC offset correction: x_std = x - μ + 'overwrite', 1); +if resample == true + % Process: Resample: 300Hz + sFileRun2 = bst_process('CallProcess', 'process_resample', sFileRun2, [], ... + 'freq', fsample, ... + 'overwrite', 1); +end + +% Process: Import MEG/EEG: Existing epochs +sFileRun3 = bst_process('CallProcess', 'process_import_data_epoch', [], [], ... + 'subjectname', subjectid, ... + 'condition', sprintf('5-Restin_cleandata_%ds', trialDuration), ... + 'datafile', {Run3File, 'FT-TIMELOCK'}, ... + 'iepochs', [], ... + 'eventtypes', '', ... + 'createcond', 0, ... + 'channelalign', 1, ... + 'usectfcomp', 1, ... + 'usessp', 1, ... + 'freq', [], ... + 'baseline', []); +% Process: DC offset correction: [All file] +sFileRun3 = bst_process('CallProcess', 'process_baseline', sFileRun3, [], ... + 'baseline', [], ... + 'sensortypes', 'MEG, EEG', ... + 'method', 'bl', ... % DC offset correction: x_std = x - μ + 'overwrite', 1); +if resample == true + % Process: Resample: 300Hz + sFileRun3 = bst_process('CallProcess', 'process_resample', sFileRun3, [], ... + 'freq', fsample, ... + 'overwrite', 1); +end + +%% !!! check the freq sampling and line freq later!!! +sFileNoise = bst_process('CallProcess', 'process_import_data_raw', [], [], ... + 'subjectname', subjectid, ... + 'datafile', {NoiseFile, '4D'}, ... + 'channelalign', 1); + + +sFilesClean = [sFileRun1, sFileRun2, sFileRun3, sFileNoise]; +sFilesRest = [sFileRun1, sFileRun2, sFileRun3]; + +%% ===== SOURCE ESTIMATION ===== + +% Process: Select file names with tag: task-rest +sFilesNoise = bst_process('CallProcess', 'process_select_tag', sFilesClean, [], ... + 'tag', '1-Rnoise', ... + 'search', 1, ... % Search the file names + 'select', 1); % Select only the files with the tag + +% Process: Compute covariance (noise or data) +bst_process('CallProcess', 'process_noisecov', sFilesNoise, [], ... + 'baseline', [], ... + 'sensortypes', 'MEG', ... + 'target', 1, ... % Noise covariance (covariance over baseline time window) + 'dcoffset', 1, ... % Block by block, to avoid effects of slow shifts in data + 'identity', 0, ... + 'copycond', 1, ... + 'copysubj', 0, ... + 'replacefile', 1); % Replace + +% Process: Compute head model +bst_process('CallProcess', 'process_headmodel', sFilesRest, [], ... + 'sourcespace', 1, ... % Cortex surface + 'meg', 3); % Overlapping spheres + +% Process: Compute sources [2018] +sSrcRest = bst_process('CallProcess', 'process_inverse_2018', sFilesRest, [], ... + 'output', 1, ... % Kernel only: shared + 'inverse', struct(... + 'Comment', 'dSPM: MEG', ... + 'InverseMethod', 'minnorm', ... + 'InverseMeasure', inverse_measure, ... + 'SourceOrient', {{'fixed'}}, ... + 'Loose', 0.2, ... + 'UseDepth', 1, ... + 'WeightExp', 0.5, ... + 'WeightLimit', 10, ... + 'NoiseMethod', 'reg', ... + 'NoiseReg', 0.1, ... + 'SnrMethod', 'fixed', ... + 'SnrRms', 1e-06, ... + 'SnrFixed', 3, ... + 'ComputeKernel', 1, ... + 'DataTypes', {{'MEG'}})); + +% Reload the Data +db_reload_database('current'); + +%% ===== IMPORT FSAverage ===== + +% DO NOT FORGET to put the HCP-MMP1.0-fsaverage inside brainstorm directory + +% Add fsaverage subject +[sSubject, iSubject] = db_add_subject('fsaverage'); + +% Add HCP-MMP1.0-fsaverage template +sTemplates = bst_get('AnatomyDefaults'); +iTemplate = find(strcmpi('HCP-MMP1.0-fsaverage', {sTemplates.Name})); +db_set_template(iSubject, sTemplates(iTemplate), 0); + +%% ===== PROJECTING SOURCES on FSAverage and SAVING SCOUTS TIME SERIES [360 scouts] as Fieldtrip Structure ===== + +% Get all the subjects in the protocol +sProtocolSubjects = bst_get('ProtocolSubjects'); +% Choose the subject +iSubject = find(strcmpi('fsaverage', {sProtocolSubjects.Subject.Name})); +% Get all the cortex surfaces of selected subject +sCortex = bst_get('SurfaceFileByType', iSubject, 'Cortex', 0); +% Choose the cortex +iCortex = find(strcmpi('cortex_15002V', {sCortex.Comment})); + +% Choose the scouts for later +ScoutSel = {'HCP-MMP1.0', {'L_10d_ROI L', 'L_10pp_ROI L', 'L_10r_ROI L', 'L_10v_ROI L', 'L_11l_ROI L', 'L_13l_ROI L', 'L_1_ROI L', 'L_23c_ROI L', 'L_23d_ROI L', 'L_24dd_ROI L', 'L_24dv_ROI L', 'L_25_ROI L', 'L_2_ROI L', 'L_31a_ROI L', 'L_31pd_ROI L', 'L_31pv_ROI L', 'L_33pr_ROI L', 'L_3a_ROI L', 'L_3b_ROI L', 'L_43_ROI L', 'L_44_ROI L', 'L_45_ROI L', 'L_46_ROI L', 'L_47l_ROI L', 'L_47m_ROI L', 'L_47s_ROI L', 'L_4_ROI L', 'L_52_ROI L', 'L_55b_ROI L', 'L_5L_ROI L', 'L_5m_ROI L', 'L_5mv_ROI L', 'L_6a_ROI L', 'L_6d_ROI L', 'L_6ma_ROI L', 'L_6mp_ROI L', 'L_6r_ROI L', 'L_6v_ROI L', 'L_7AL_ROI L', 'L_7Am_ROI L', 'L_7PC_ROI L', 'L_7PL_ROI L', 'L_7Pm_ROI L', 'L_7m_ROI L', 'L_8Ad_ROI L', 'L_8Av_ROI L', 'L_8BL_ROI L', 'L_8BM_ROI L', 'L_8C_ROI L', 'L_9-46d_ROI L', 'L_9a_ROI L', 'L_9m_ROI L', 'L_9p_ROI L', 'L_A1_ROI L', 'L_A4_ROI L', 'L_A5_ROI L', 'L_AAIC_ROI L', 'L_AIP_ROI L', 'L_AVI_ROI L', 'L_DVT_ROI L', 'L_EC_ROI L', 'L_FEF_ROI L', 'L_FFC_ROI L', 'L_FOP1_ROI L', 'L_FOP2_ROI L', 'L_FOP3_ROI L', 'L_FOP4_ROI L', 'L_FOP5_ROI L', 'L_FST_ROI L', 'L_H_ROI L', 'L_IFJa_ROI L', 'L_IFJp_ROI L', 'L_IFSa_ROI L', 'L_IFSp_ROI L', 'L_IP0_ROI L', 'L_IP1_ROI L', 'L_IP2_ROI L', 'L_IPS1_ROI L', 'L_Ig_ROI L', 'L_LBelt_ROI L', 'L_LIPd_ROI L', 'L_LIPv_ROI L', 'L_LO1_ROI L', 'L_LO2_ROI L', 'L_LO3_ROI L', 'L_MBelt_ROI L', 'L_MIP_ROI L', 'L_MI_ROI L', 'L_MST_ROI L', 'L_MT_ROI L', 'L_OFC_ROI L', 'L_OP1_ROI L', 'L_OP2-3_ROI L', 'L_OP4_ROI L', 'L_PBelt_ROI L', 'L_PCV_ROI L', 'L_PEF_ROI L', 'L_PF_ROI L', 'L_PFcm_ROI L', 'L_PFm_ROI L', 'L_PFop_ROI L', 'L_PFt_ROI L', 'L_PGi_ROI L', 'L_PGp_ROI L', 'L_PGs_ROI L', 'L_PHA1_ROI L', 'L_PHA2_ROI L', 'L_PHA3_ROI L', 'L_PHT_ROI L', 'L_PH_ROI L', 'L_PIT_ROI L', 'L_PI_ROI L', 'L_POS1_ROI L', 'L_POS2_ROI L', 'L_PSL_ROI L', 'L_PeEc_ROI L', 'L_Pir_ROI L', 'L_PoI1_ROI L', 'L_PoI2_ROI L', 'L_PreS_ROI L', 'L_ProS_ROI L', 'L_RI_ROI L', 'L_RSC_ROI L', 'L_SCEF_ROI L', 'L_SFL_ROI L', 'L_STGa_ROI L', 'L_STSda_ROI L', 'L_STSdp_ROI L', 'L_STSva_ROI L', 'L_STSvp_ROI L', 'L_STV_ROI L', 'L_TA2_ROI L', 'L_TE1a_ROI L', 'L_TE1m_ROI L', 'L_TE1p_ROI L', 'L_TE2a_ROI L', 'L_TE2p_ROI L', 'L_TF_ROI L', 'L_TGd_ROI L', 'L_TGv_ROI L', 'L_TPOJ1_ROI L', 'L_TPOJ2_ROI L', 'L_TPOJ3_ROI L', 'L_V1_ROI L', 'L_V2_ROI L', 'L_V3A_ROI L', 'L_V3B_ROI L', 'L_V3CD_ROI L', 'L_V3_ROI L', 'L_V4_ROI L', 'L_V4t_ROI L', 'L_V6A_ROI L', 'L_V6_ROI L', 'L_V7_ROI L', 'L_V8_ROI L', 'L_VIP_ROI L', 'L_VMV1_ROI L', 'L_VMV2_ROI L', 'L_VMV3_ROI L', 'L_VVC_ROI L', 'L_a10p_ROI L', 'L_a24_ROI L', 'L_a24pr_ROI L', 'L_a32pr_ROI L', 'L_a47r_ROI L', 'L_a9-46v_ROI L', 'L_d23ab_ROI L', 'L_d32_ROI L', 'L_i6-8_ROI L', 'L_p10p_ROI L', 'L_p24_ROI L', 'L_p24pr_ROI L', 'L_p32_ROI L', 'L_p32pr_ROI L', 'L_p47r_ROI L', 'L_p9-46v_ROI L', 'L_pOFC_ROI L', 'L_s32_ROI L', 'L_s6-8_ROI L', 'L_v23ab_ROI L', 'R_10d_ROI R', 'R_10pp_ROI R', 'R_10r_ROI R', 'R_10v_ROI R', 'R_11l_ROI R', 'R_13l_ROI R', 'R_1_ROI R', 'R_23c_ROI R', 'R_23d_ROI R', 'R_24dd_ROI R', 'R_24dv_ROI R', 'R_25_ROI R', 'R_2_ROI R', 'R_31a_ROI R', 'R_31pd_ROI R', 'R_31pv_ROI R', 'R_33pr_ROI R', 'R_3a_ROI R', 'R_3b_ROI R', 'R_43_ROI R', 'R_44_ROI R', 'R_45_ROI R', 'R_46_ROI R', 'R_47l_ROI R', 'R_47m_ROI R', 'R_47s_ROI R', 'R_4_ROI R', 'R_52_ROI R', 'R_55b_ROI R', 'R_5L_ROI R', 'R_5m_ROI R', 'R_5mv_ROI R', 'R_6a_ROI R', 'R_6d_ROI R', 'R_6ma_ROI R', 'R_6mp_ROI R', 'R_6r_ROI R', 'R_6v_ROI R', 'R_7AL_ROI R', 'R_7Am_ROI R', 'R_7PC_ROI R', 'R_7PL_ROI R', 'R_7Pm_ROI R', 'R_7m_ROI R', 'R_8Ad_ROI R', 'R_8Av_ROI R', 'R_8BL_ROI R', 'R_8BM_ROI R', 'R_8C_ROI R', 'R_9-46d_ROI R', 'R_9a_ROI R', 'R_9m_ROI R', 'R_9p_ROI R', 'R_A1_ROI R', 'R_A4_ROI R', 'R_A5_ROI R', 'R_AAIC_ROI R', 'R_AIP_ROI R', 'R_AVI_ROI R', 'R_DVT_ROI R', 'R_EC_ROI R', 'R_FEF_ROI R', 'R_FFC_ROI R', 'R_FOP1_ROI R', 'R_FOP2_ROI R', 'R_FOP3_ROI R', 'R_FOP4_ROI R', 'R_FOP5_ROI R', 'R_FST_ROI R', 'R_H_ROI R', 'R_IFJa_ROI R', 'R_IFJp_ROI R', 'R_IFSa_ROI R', 'R_IFSp_ROI R', 'R_IP0_ROI R', 'R_IP1_ROI R', 'R_IP2_ROI R', 'R_IPS1_ROI R', 'R_Ig_ROI R', 'R_LBelt_ROI R', 'R_LIPd_ROI R', 'R_LIPv_ROI R', 'R_LO1_ROI R', 'R_LO2_ROI R', 'R_LO3_ROI R', 'R_MBelt_ROI R', 'R_MIP_ROI R', 'R_MI_ROI R', 'R_MST_ROI R', 'R_MT_ROI R', 'R_OFC_ROI R', 'R_OP1_ROI R', 'R_OP2-3_ROI R', 'R_OP4_ROI R', 'R_PBelt_ROI R', 'R_PCV_ROI R', 'R_PEF_ROI R', 'R_PF_ROI R', 'R_PFcm_ROI R', 'R_PFm_ROI R', 'R_PFop_ROI R', 'R_PFt_ROI R', 'R_PGi_ROI R', 'R_PGp_ROI R', 'R_PGs_ROI R', 'R_PHA1_ROI R', 'R_PHA2_ROI R', 'R_PHA3_ROI R', 'R_PHT_ROI R', 'R_PH_ROI R', 'R_PIT_ROI R', 'R_PI_ROI R', 'R_POS1_ROI R', 'R_POS2_ROI R', 'R_PSL_ROI R', 'R_PeEc_ROI R', 'R_Pir_ROI R', 'R_PoI1_ROI R', 'R_PoI2_ROI R', 'R_PreS_ROI R', 'R_ProS_ROI R', 'R_RI_ROI R', 'R_RSC_ROI R', 'R_SCEF_ROI R', 'R_SFL_ROI R', 'R_STGa_ROI R', 'R_STSda_ROI R', 'R_STSdp_ROI R', 'R_STSva_ROI R', 'R_STSvp_ROI R', 'R_STV_ROI R', 'R_TA2_ROI R', 'R_TE1a_ROI R', 'R_TE1m_ROI R', 'R_TE1p_ROI R', 'R_TE2a_ROI R', 'R_TE2p_ROI R', 'R_TF_ROI R', 'R_TGd_ROI R', 'R_TGv_ROI R', 'R_TPOJ1_ROI R', 'R_TPOJ2_ROI R', 'R_TPOJ3_ROI R', 'R_V1_ROI R', 'R_V2_ROI R', 'R_V3A_ROI R', 'R_V3B_ROI R', 'R_V3CD_ROI R', 'R_V3_ROI R', 'R_V4_ROI R', 'R_V4t_ROI R', 'R_V6A_ROI R', 'R_V6_ROI R', 'R_V7_ROI R', 'R_V8_ROI R', 'R_VIP_ROI R', 'R_VMV1_ROI R', 'R_VMV2_ROI R', 'R_VMV3_ROI R', 'R_VVC_ROI R', 'R_a10p_ROI R', 'R_a24_ROI R', 'R_a24pr_ROI R', 'R_a32pr_ROI R', 'R_a47r_ROI R', 'R_a9-46v_ROI R', 'R_d23ab_ROI R', 'R_d32_ROI R', 'R_i6-8_ROI R', 'R_p10p_ROI R', 'R_p24_ROI R', 'R_p24pr_ROI R', 'R_p32_ROI R', 'R_p32pr_ROI R', 'R_p47r_ROI R', 'R_p9-46v_ROI R', 'R_pOFC_ROI R', 'R_s32_ROI R', 'R_s6-8_ROI R', 'R_v23ab_ROI R'}}; + +% empty Fieldtrip data to export later +FT_data = struct; + + +%%% 3-Restin +% Process: Select Sources in 3-Restin_cleandata_ +Sources_1 = bst_process('CallProcess', 'process_select_files_results', [], [], ... + 'subjectname', subjectid, ... + 'condition', sprintf('3-Restin_cleandata_%ds', trialDuration), ... + 'includebad', 0); + +length_s_1 = sum(~cellfun(@isempty,{Sources_1.FileName})); +for idx = 1 : length_s_1 + % Projecting sources on FSAverage + ResultsFile = cellstr(Sources_1(idx).FileName); + Projection = bst_project_sources(ResultsFile, sCortex(iCortex).FileName); + + % Converting to Fieldtrip structure + ScoutFunc = 'Max'; + isNorm = 0; + [ftData, sInput, VertConn] = out_fieldtrip_results( Projection{1}, ScoutSel, ScoutFunc, [], isNorm); + + FT_data.trial{idx} = ftData.pow; + FT_data.trialinfo(idx) = {sInput.Comment}; + + display(FT_data.trialinfo(idx)) + + % Process: Delete selected files + bst_process('CallProcess', 'process_delete', Projection, [], ... + 'target', 1); % Delete selected files +end + + +%%% 4-Restin +% Process: Select Sources in 4-Restin_cleandata_ +Sources_2 = bst_process('CallProcess', 'process_select_files_results', [], [], ... + 'subjectname', subjectid, ... + 'condition', sprintf('4-Restin_cleandata_%ds', trialDuration), ... + 'includebad', 0); + +length_s_2 = sum(~cellfun(@isempty,{Sources_2.FileName})); +for idx = 1+length_s_1 : length_s_1+length_s_2 + % Projecting sources on FSAverage + ResultsFile = cellstr(Sources_2(idx-length_s_1).FileName); + Projection = bst_project_sources(ResultsFile, sCortex(iCortex).FileName); + + % Converting to Fieldtrip structure + ScoutFunc = 2; % Max + isNorm = 0; + [ftData, sInput, VertConn] = out_fieldtrip_results( Projection{1}, ScoutSel, ScoutFunc, [], isNorm); + + FT_data.trial{idx} = ftData.pow; + FT_data.trialinfo(idx) = {sInput.Comment}; + + display(FT_data.trialinfo(idx)) + + + % Process: Delete selected files + bst_process('CallProcess', 'process_delete', Projection, [], ... + 'target', 1); % Delete selected files +end + + +%%% 5-Restin +% Process: Select Sources in 5-Restin_cleandata_ +Sources_3 = bst_process('CallProcess', 'process_select_files_results', [], [], ... + 'subjectname', subjectid, ... + 'condition', sprintf('5-Restin_cleandata_%ds', trialDuration), ... + 'includebad', 0); + +length_s_3 = sum(~cellfun(@isempty,{Sources_3.FileName})); +for idx = 1+length_s_1+length_s_2 : length_s_1+length_s_2+length_s_3 + % Projecting sources on FSAverage + ResultsFile = cellstr(Sources_3(idx-length_s_1-length_s_2).FileName); + Projection = bst_project_sources(ResultsFile, sCortex(iCortex).FileName); + + % Converting to Fieldtrip structure + ScoutFunc = 2; % Max + isNorm = 0; + [ftData, sInput, VertConn] = out_fieldtrip_results( Projection{1}, ScoutSel, ScoutFunc, [], isNorm); + + FT_data.trial{idx} = ftData.pow; + FT_data.trialinfo(idx) = {sInput.Comment}; + + display(FT_data.trialinfo(idx)) + + % Process: Delete selected files + bst_process('CallProcess', 'process_delete', Projection, [], ... + 'target', 1); % Delete selected files +end + + +time = cell(1, length(FT_data.trial)); +time(:)={ftData.time}; +FT_data.time = time; + +FT_data.label = sInput.RowNames; +FT_data.fsample = fsample; + +FT_data.trialDuration = trialDuration; + + +% Reload the Data +db_reload_database('current'); + +%% ===== WRITE THE SUBJECT STRUCTURE TO DISK ===== + +% % Make it memory efficent +% FT_data = ft_struct2single(FT_data); +% fsaverage = ft_struct2single(fsaverage); + +save(fullfile(OutputDir_scouts, [subjectid, '_scoutsTimeSeries_', sprintf('%ds.mat', trialDuration)]), 'FT_data', '-v7.3'); + +filename = fullfile(OutputDir_fsaverage, ['fsaverage_', sprintf('%ds.mat', trialDuration)]); +if ~isfile(filename) + % Save information related to anatomy + fsaverage = in_tess_bst( bst_fullfile(ProtocolInfo.SUBJECTS, sCortex(iCortex).FileName), 0); + save(filename, 'fsaverage', '-v7.3'); +end + +%% ===== THE END ===== +disp(['THE END for ', subjectid, '!']) + +if cleanProtocol == true + gui_brainstorm('DeleteProtocol', ProtocolName); +end + +toc +end diff --git a/README.md b/README.md new file mode 100644 index 0000000..92d81c8 --- /dev/null +++ b/README.md @@ -0,0 +1,24 @@ +# Functional connectivity fingerprints of the frontal eye fields and inferior frontal junction in the dorsal vs. ventral prefrontal cortex +Inside this repository, you can find the analysis pipeline for our research on the contrasting functional connectivity profiles of FEF and IFJ. + +![Figure_1](https://user-images.githubusercontent.com/44211738/172026736-e0b69fc6-d48d-4988-baf7-a6447d9288fd.png) +FIGURE 1 Predominant functional connectivity maps of FEF and IFJa across frequency bands and hemispheres for (A) oPEC and (B) iCOH metrics (two-sided Wilcoxon signed-rank test, p < 0.05, FDR-corrected for 33 ROIs). The results shown here are based on 2 s epoch segmentation. + +## Abstract +#### Orhan Soyuhos | Daniel Baldauf +#### Centre for Mind/Brain Sciences (CIMeC), University of Trento, 38068 Rovereto, Italy. +#### E-mail addresses: orhan.soyuhos@studenti.unitn.it (O. Soyuhos), daniel.baldauf@unitn.it (D. Baldauf). + +Neuroimaging evidence suggests that the frontal eye field (FEF) and inferior frontal junction (IFJ) govern the encoding of spatial and non-spatial (such as feature- or object-based) representations, respectively, both during visual attention and working memory tasks. However, it is still unclear whether such contrasting functional segregation is also reflected in their underlying functional connectivity patterns. Here, we hypothesized that FEF has predominant functional coupling with spatiotopically organized regions in the dorsal ('where') visual stream, whereas IFJ has predominant functional connectivity with the ventral ('what') visual stream. We applied seed-based functional connectivity analyses to temporally high-resolving resting-state magnetoencephalography (MEG) recordings. We parcellated the brain according to the multimodal Glasser atlas and tested, for various frequency bands, whether the spontaneous activity of each parcel in the ventral and dorsal visual pathway has predominant functional connectivity with FEF or IFJ. The results show that FEF has a robust power correlation with the dorsal visual pathway in beta and gamma bands. In contrast, anterior IFJ (IFJa) has a strong power coupling with the ventral visual stream in delta, beta, and gamma oscillations. Moreover, while FEF is directly phase-coupled with the superior parietal lobe in the beta band, IFJa is directly phase-coupled with the middle and inferior temporal cortex in delta and gamma oscillations. We argue that these intrinsic connectivity fingerprints are congruent with each brain region's function. Therefore, we conclude that FEF and IFJ have dissociable connectivity patterns that fit their respective functional roles in spatial vs. non-spatial top-down attention and working memory control. + +#### KEYWORDS +frontal eye field, inferior frontal junction, functional connectivity, prefrontal cortex, visual pathways, visual attention, working memory. + +### Acknowledgement +This research is supported by Ministero degli Affari Esteri e della Cooperazione Internazionale (to O.S.) and Fondazione Cassa Di Risparmio Di Trento E Rovereto (to D.B). + +### Data Availability Statement +The data used in the current study provided and publicly made available by the 1200 Subjects Release of Human Connectome Project (Larson-Prior et al., 2013; Van Essen et al., 2013) in ConnectomeDB (db.humanconnectome.org; Hodge et al., 2016). + +### ORCID +Orhan Soyuhos https://orcid.org/0000-0002-2858-3979; Daniel Baldauf https://orcid.org/0000-0001-5764-506X diff --git a/Visualize_Results/ReadMe.md b/Visualize_Results/ReadMe.md new file mode 100644 index 0000000..4604e7c --- /dev/null +++ b/Visualize_Results/ReadMe.md @@ -0,0 +1,8 @@ +## Visualize_Results folder includes MATLAB code to apply statistical masks and visualize the resulting functional connectivity matrices from the AnalysisPipeline step. + +The fsaverage folder contains the MATLAB files for the fsaverage surface-based template with several smoothing values. + +The results folder has group-level connectivity matrices from the AnalysisPipeline step. We could not upload all connectivity matrices due to limits on the file size but provided some examples to run the code. + +The scripts folder has the code for this part. A detailed explanation is provided in its ReadMe.md file. + diff --git a/Visualize_Results/fsaverage/fsaverage.mat b/Visualize_Results/fsaverage/fsaverage.mat new file mode 100644 index 0000000..9cbe904 Binary files /dev/null and b/Visualize_Results/fsaverage/fsaverage.mat differ diff --git a/Visualize_Results/fsaverage/fsaverage_inflated_100.mat b/Visualize_Results/fsaverage/fsaverage_inflated_100.mat new file mode 100644 index 0000000..b964e2e Binary files /dev/null and b/Visualize_Results/fsaverage/fsaverage_inflated_100.mat differ diff --git a/Visualize_Results/fsaverage/fsaverage_inflated_70.mat b/Visualize_Results/fsaverage/fsaverage_inflated_70.mat new file mode 100644 index 0000000..10bc51e Binary files /dev/null and b/Visualize_Results/fsaverage/fsaverage_inflated_70.mat differ diff --git a/Visualize_Results/results/ReadMe.md b/Visualize_Results/results/ReadMe.md new file mode 100644 index 0000000..499a0eb --- /dev/null +++ b/Visualize_Results/results/ReadMe.md @@ -0,0 +1,3 @@ +## Example functional connectivity matrix +The results for power envelope correlation based on oPEC metric. +- To run the code please download from Google Drive and locate into the current directory: https://drive.google.com/drive/folders/1XbMTGpCcEaH_bcEQUyZkYxxOKqOy0r5O?usp=sharing diff --git a/Visualize_Results/scripts/RUNME_visualizeResults_effectiveConn.m b/Visualize_Results/scripts/RUNME_visualizeResults_effectiveConn.m new file mode 100644 index 0000000..bd7c819 --- /dev/null +++ b/Visualize_Results/scripts/RUNME_visualizeResults_effectiveConn.m @@ -0,0 +1,70 @@ +close all +clear + +addpath(genpath(fullfile(pwd, 'helper_fun'))); + +%% MEANING of colors +% y: target region +% +% GREEN : seed --> y direction is significant. A seed region explains the varience in a target better. +% ~ sending signal +% MAGENTA : y--> seed direction is significant. A target region explains the varience in a seed better. +% ~ receiving signal + +%% Select +statistics.analysis_type = 'ROI'; % 'ROI' 'exploratory' +% 'ROI' : ROI analysis per hemisphere. FDR correction for 33+3 seed regions. +% 'exploratory' : exploratory analysis per hemisphere. FDR correction for 180 regions. + +duration = '2s'; %fixed % only '2s' +conn_metric = 'pdc'; % 'pdc' +band_name = 'beta'; % 'delta' 'theta' 'alpha' 'beta' 'gamma' +seed_target = 'right-right'; % 'right-right' 'right-left' 'left-left' 'left-right' + +%% +% figures +fsaverage_fig = true; +violinplot = false; + +% Paired Wilcoxon paired test against zero +statistics.alpha = 0.001; % 0.05 0.01 0.001 +% FDR correction +statistics.corrected = true; % true false +statistics.correction = 'FDR'; + +% to show the contours of parcels; +%%keep it false for a faster processing +flag.contours = false; + +% save the figures automatic? +save_fig = false; +dir_fig = load_path('save_fig'); + +%% +colorlim_value = [-4 4]; +%%% fixed +circular_graph = false; +font_size = 15; +flag.fsaverage = fsaverage_fig; % to show the figure. +flag.circular_graph = circular_graph; +flag.font_size = font_size; +flag.save_fig = save_fig; +flag.duration = duration; +flag.dir_fig = dir_fig; +inputfile_conn = fullfile(load_path(['Outputs_' duration '_55subjs_' conn_metric]), ['55Subjects_groupResults_' duration '.mat']); + +%% Run the function +cimec.inputfile_conn = inputfile_conn; +cimec.band_name = band_name; +cimec.seed_target = seed_target; +cimec.conn_metric = conn_metric; +cimec.statistics = statistics; +cimec.colorlim_value = colorlim_value; +cimec.flag = flag; + +outputs = helper_effectiveConn_fsaverage(cimec); + +%% violinplot +if violinplot == true + helper_violinplot(outputs, conn_metric, band_name, seed_target, duration) +end \ No newline at end of file diff --git a/Visualize_Results/scripts/RUNME_visualizeResults_functionConn.m b/Visualize_Results/scripts/RUNME_visualizeResults_functionConn.m new file mode 100644 index 0000000..33fddaa --- /dev/null +++ b/Visualize_Results/scripts/RUNME_visualizeResults_functionConn.m @@ -0,0 +1,62 @@ +clear; +close all; + +addpath(genpath(fullfile(pwd, 'helper_fun'))); +% colormapeditor +set(0, 'DefaultFigureColormap', cmap_rbw); + +%% Select +statistics.analysis_type = 'ROI'; % 'ROI' 'exploratory' +% 'ROI' : ROI analysis per hemisphere. FDR correction for 33 regions. +% 'exploratory' : exploratory analysis per hemisphere. FDR correction for 180 regions. + +duration = '2s'; % '2s' '5s' '10s' +conn_metric = 'powcorrOrtho'; % 'powcorrOrtho' 'icoh' 'dwPLI' +band_name = 'delta'; % 'delta' 'theta' 'alpha' 'beta' 'gamma' 'collapseFreq' +seed_target = 'right-right' ; % 'right-right' 'right-left' 'left-left' 'left-right' + +% contrast +contrast = true; + +% figures +fsaverage_fig = true; +circular_graph = false; +font_size = 15; + +% Wilcoxon signed rank test +statistics.alpha = 0.05; % 0.05 0.01 0.001 +% FDR correction +statistics.corrected = true; % true false +statistics.correction = 'FDR'; + +% to show the contours of parcels; +%%keep it false for a faster processing +flag.contours = false; + +% save the figures automatic? +save_fig = false; +dir_fig = load_path('save_fig'); + +%% fixed +flag.fsaverage = fsaverage_fig; % to show the figure. +flag.circular_graph = circular_graph; +flag.font_size = font_size; +flag.save_fig = save_fig; +flag.dir_fig = dir_fig; +flag.duration = duration; +inputfile_conn = fullfile(load_path(['Outputs_' duration '_55subjs_' conn_metric]), ['55Subjects_groupResults_' duration '.mat']); + +%% Run the function +cimec.inputfile_conn = inputfile_conn; +cimec.band_name = band_name; +cimec.seed_target = seed_target; +cimec.conn_metric = conn_metric; +cimec.statistics = statistics; +cimec.contrast = contrast; +cimec.flag = flag; + +if strcmp(band_name, 'collapseFreq') + outputs = helper_functionConn_fsaverage_freqCollapsed(cimec); +else + outputs = helper_functionConn_fsaverage(cimec); +end diff --git a/Visualize_Results/scripts/RUNME_visualizeResults_functionConn_collapseTargets.m b/Visualize_Results/scripts/RUNME_visualizeResults_functionConn_collapseTargets.m new file mode 100644 index 0000000..f8fca8a --- /dev/null +++ b/Visualize_Results/scripts/RUNME_visualizeResults_functionConn_collapseTargets.m @@ -0,0 +1,294 @@ +clear +close all + +addpath(genpath(fullfile(pwd, 'helper_fun'))); + +%% Select +analysis_type = 'ROI'; % 'ROI' 'exploratory' +% 'ROI' : ROI analysis per hemisphere. FDR correction for 33 regions. +% 'exploratory' : exploratory analysis per hemisphere. FDR correction for 180 regions. + +duration = '2s'; % '2s' '5s' '10s' +conn_metric = 'powcorrOrtho'; % 'powcorrOrtho' 'icoh' 'dwPLI' +band_name = 'collapseFreq'; % 'delta' 'theta' 'alpha' 'beta' 'gamma' 'collapseFreq' +seed_targets = {'right-right' 'left-left'}; % 'right-right' 'right-left' 'left-left' 'left-right' + +% contrast +contrast = true; + +% figures +fsaverage_fig = true; +circular_graph = true; + +font_size = 18; %18 +flag.contours = true; + +% save? +flag.save_fig = false; +flag.dir_fig = load_path('save_fig'); + +%% fixed +flag.script = 'collapseTargets'; +flag.seed_targets = seed_targets; + +outputs_first = helper_visualizeResults_functionConn_leftORright(analysis_type, contrast, duration, conn_metric, band_name, seed_targets{1}); +outputs_second = helper_visualizeResults_functionConn_leftORright(analysis_type, contrast, duration, conn_metric, band_name, seed_targets{2}); + +alpha = 0.05; +corrected = true; +correction ='FDR'; +flag.duration = duration; + +%% +seed_target = 'right-right'; +idx_IFJa = outputs_first.helper.idx_IFJa; +idx_IFJp = outputs_first.helper.idx_IFJp; +idx_FEF = outputs_first.helper.idx_FEF; +ROIs = outputs_first.helper.ROIs; +label = outputs_first.label; + +if strcmp(band_name, 'collapseFreq') + h_IFJa_all = logical(outputs_first.stat.h_IFJa_all) + logical(outputs_second.stat.h_IFJa_all); + h_IFJp_all = logical(outputs_first.stat.h_IFJp_all) + logical(outputs_second.stat.h_IFJp_all); + if contrast == false + h_FEF_all = logical(outputs_first.stat.h_FEF_all) + logical(outputs_second.stat.h_FEF_all); + end +else + h_IFJa_all = logical(outputs_first.stat.h_adj_IFJa) + logical(outputs_second.stat.h_adj_IFJa); + h_IFJp_all = logical(outputs_first.stat.h_adj_IFJp) + logical(outputs_second.stat.h_adj_IFJp); + if contrast == false + h_FEF_all = logical(outputs_first.stat.h_adj_FEF) + logical(outputs_second.stat.h_adj_FEF); + end +end + +tmp_cohspctrm = outputs_first.conn_matrix(181:360,181:360) + outputs_second.conn_matrix(1:180,1:180); +conn_matrix = zeros(360,360); +conn_matrix(181:360,181:360) = tmp_cohspctrm; +conn_matrix(idx_IFJa, ROIs.idx) = conn_matrix(idx_IFJa, ROIs.idx)./h_IFJa_all; +conn_matrix(idx_IFJp, ROIs.idx) = conn_matrix(idx_IFJp, ROIs.idx)./h_IFJp_all; +if contrast == false + conn_matrix(idx_FEF, ROIs.idx) = conn_matrix(idx_FEF, ROIs.idx)./h_FEF_all; +end +conn_matrix(isinf(conn_matrix)) = 0; +conn_matrix(isnan(conn_matrix)) = 0; +conn_matrix(ROIs.idx, idx_IFJa) = conn_matrix(idx_IFJa, ROIs.idx); +conn_matrix(ROIs.idx, idx_IFJp) = conn_matrix(idx_IFJp, ROIs.idx); +if contrast == false + conn_matrix(ROIs.idx, idx_FEF) = conn_matrix(idx_FEF, ROIs.idx); +end + +conn_mat = struct; +conn_mat.cohspctrm = conn_matrix; +conn_mat.label = outputs_first.label; +conn_mat.brainordinate = outputs_first.helper.brainordinate; + +if contrast == true + colorlim_value = [-4 4]; +else + colorlim_value = [0 6]; +end + +%% fsaverage +if fsaverage_fig == true + + % to make black the areas outside ROI + if and(contrast == true, strcmp(analysis_type, 'ROI')) + cbar = cmap_rbw; + n = size(cbar, 1); + up = cbar(1:n/2,:); + middle = [0.5 0.5 0.5]; + down = cbar(n/2+1:n,:); + set(0, 'DefaultFigureColormap', [up;middle;down]); + + conn_mat.cohspctrm(ROIs.idx(~h_IFJa_all),idx_IFJa) = 0.04; + conn_mat.cohspctrm(ROIs.idx(~h_IFJp_all),idx_IFJp) = 0.04; + conn_mat.cohspctrm(idx_IFJa,ROIs.idx(~h_IFJa_all)) = 0.04; + conn_mat.cohspctrm(idx_IFJp,ROIs.idx(~h_IFJp_all)) = 0.04; + + elseif and(contrast == false, strcmp(analysis_type, 'ROI')) + cbar = hot; + base1 = [1 1 1]; + base2 = [0.25 0.25 0.25]; + set(0, 'DefaultFigureColormap', [base1; base2; cbar]); + + conn_mat.cohspctrm(ROIs.idx(~h_IFJa_all),idx_IFJa) = -0.06; + conn_mat.cohspctrm(ROIs.idx(~h_IFJp_all),idx_IFJp) = -0.06; + conn_mat.cohspctrm(ROIs.idx(~h_FEF_all),idx_FEF) = -0.06; + conn_mat.cohspctrm(idx_IFJa,ROIs.idx(~h_IFJa_all)) = -0.06; + conn_mat.cohspctrm(idx_IFJp,ROIs.idx(~h_IFJp_all)) = -0.06; + conn_mat.cohspctrm(idx_FEF,ROIs.idx(~h_FEF_all)) = -0.06; + + minv = -0.06; + maxv = 6; + colorlim_value = [minv maxv]; + elseif contrast == false + set(0, 'DefaultFigureColormap', hot); + end + + % in order to run the function + conn_mat.brainordinate.pos = conn_mat.brainordinate.pos*1000; + + %% save the figure + if flag.save_fig == true + + cimec = struct; + cimec.conn_metric = conn_metric; + cimec.band_name = band_name; + cimec.seed_target = seed_target; + cimec.statistics.analysis_type = analysis_type; + cimec.statistics.alpha = alpha; + cimec.statistics.corrected = corrected; + cimec.statistics.correction = correction; + cimec.flag = flag; + + if seed_target(1) == 'r' + seed = 'IFJa'; + cimec.seed = seed; + pos2d_IFJa_R = [52.5323 -69.3638 71.0650; 52.5323 69.4183 71.0650]; + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_IFJa_R); + helper_save_figure(cimec); + close all; + + seed = 'IFJp'; + cimec.seed = seed; + pos2d_IFJp_R = [45.6139 -69.3638 76.8706 ; 45.6139 69.4183 76.8706]; + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_IFJp_R); + helper_save_figure(cimec); + close all; + + if contrast == false + seed = 'FEF'; + cimec.seed = seed; + pos2d_FEF_R = [25.8639 -69.3638 96.3111; 25.8639 69.4183 96.3111]; + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_FEF_R); + helper_save_figure(cimec); + close all; + end + + elseif seed_target(1) == 'l' + seed = 'IFJa'; + cimec.seed = seed; + pos2d_IFJa_L = [50.5105 69.4183 69.5634 ; 53.1465 -69.3638 68.2860]; + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_IFJa_L); + helper_save_figure(cimec); + close all; + + seed = 'IFJp'; + cimec.seed = seed; + pos2d_IFJp_L = [43.6193 69.4183 80.0619 ; 43.6193 -69.3638 80.0619]; + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_IFJp_L); + helper_save_figure(cimec); + close all; + + if contrast == false + seed = 'FEF'; + cimec.seed = seed; + pos2d_FEF_L = [30.9891 69.4183 105.2148 ; 34.2627 -69.3638 98.8225]; + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_FEF_L); + helper_save_figure(cimec); + close all; + end + end + else + tutorial_nwa_connectivityviewer(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours); + end +end + +%% circularGraph +if circular_graph == true + results = struct; + results.stat.analysis_type = analysis_type; + results.helper.band_name = band_name; + if strcmp(band_name, 'collapseFreq') + results.stat.h_IFJa_all = outputs_first.stat.h_IFJa_all + outputs_second.stat.h_IFJa_all; + results.stat.h_IFJp_all = outputs_first.stat.h_IFJp_all + outputs_second.stat.h_IFJp_all; + if contrast == false + results.stat.h_FEF_all = outputs_first.stat.h_FEF_all + outputs_second.stat.h_FEF_all; + end + else + results.stat.h_adj_IFJa = outputs_first.stat.h_adj_IFJa + outputs_second.stat.h_adj_IFJa; + results.stat.h_adj_IFJp = outputs_first.stat.h_adj_IFJp + outputs_second.stat.h_adj_IFJp; + if contrast == false + results.stat.h_adj_FEF = outputs_first.stat.h_adj_FEF + outputs_second.stat.h_adj_FEF; + end + end + + tmp_data = conn_matrix; + tmp_data(tmp_data == 1) = NaN; + tmp_data(idx_IFJa, idx_IFJa) = 0; + tmp_data(idx_IFJp, idx_IFJp) = 0; + tmp_data(idx_FEF, idx_FEF) = 0; + results.conn_matrix = tmp_data; + + results.helper.idx_FEF = idx_FEF; + results.helper.idx_IFJa = idx_IFJa; + results.helper.idx_IFJp = idx_IFJp; + results.helper.ROIs = ROIs; + results.label = label; + + %% + flag.collapseTargets = true; + flag.contrast = contrast; + results.cimec = outputs_first.cimec; + results.cimec.band_name = band_name; + results.cimec.seed_target = char(string(join(seed_targets))); + fontSize = font_size; + + helper_functionConn_circularGraph(results, fontSize, flag); +end + +%% helper function +function outputs = helper_visualizeResults_functionConn_leftORright(analysis_type, contrast, duration, conn_metric, band_name, seed_target) + +addpath(genpath(fullfile(pwd, 'helper_fun'))); +% colormapeditor +set(0, 'DefaultFigureColormap', cmap_rbw); + +statistics.analysis_type = analysis_type; % 'ROI' 'exploratory' +% 'ROI' : ROI analysis per hemisphere. FDR correction for 33 regions. +% 'exploratory' : exploratory analysis per hemisphere. FDR correction for 180 regions. + +% figures +fsaverage_fig = false; +circular_graph = false; +font_size = 15; + +% Wilcoxon signed rank test +statistics.alpha = 0.05; % 0.05 0.01 0.001 +% FDR correction +statistics.corrected = true; % true false +statistics.correction = 'FDR'; + +% to show the contours of parcels; +%%keep it false for a faster processing +flag.contours = false; + +% save the figures automatic? +save_fig = false; +dir_fig = load_path('save_fig'); + +%% fixed +flag.fsaverage = fsaverage_fig; % to show the figure. +flag.circular_graph = circular_graph; +flag.font_size = font_size; +flag.save_fig = save_fig; +flag.dir_fig = dir_fig; +flag.duration = duration; +inputfile_conn = fullfile(load_path(['Outputs_' duration '_55subjs_' conn_metric]), ['55Subjects_groupResults_' duration '.mat']); + +%% Run the function +cimec.inputfile_conn = inputfile_conn; +cimec.band_name = band_name; +cimec.seed_target = seed_target; +cimec.conn_metric = conn_metric; +cimec.statistics = statistics; +cimec.contrast = contrast; +cimec.flag = flag; + +if strcmp(band_name, 'collapseFreq') + outputs = helper_functionConn_fsaverage_freqCollapsed(cimec); +else + outputs = helper_functionConn_fsaverage(cimec); +end + +end diff --git a/Visualize_Results/scripts/ReadMe.md b/Visualize_Results/scripts/ReadMe.md new file mode 100644 index 0000000..6d27c10 --- /dev/null +++ b/Visualize_Results/scripts/ReadMe.md @@ -0,0 +1,15 @@ +## The code to visualize the results for predominant functional connectivities. + +- The helper_fun folder includes custom functions used in the scripts below. Please fill the directory information inside load_path.m file. + +### Main scripts + +- RUNME_visualizeResults_functionConn.m: To apply statistical analysis and visualize FEF and IFJ's predominant functional connectivities in the ipsilateral or contralateral sides. It produces the results on fsaverage cortex and circular graphs. + +- RUNME_visualizeResults_functionConn_collapseTargets.m: Same as the above script but irrespective of laterilazation. + +- RUNME_visualizeResults_effectiveConn.m: To apply statistical analysis and visualize FEF and IFJ's direction of interaction with the visual streams or the rest of the brain. + +### Mass production of the figures + +- helper_saveFig_effectiveConn.m and helper_saveFig_functionalConn.m is for the mass production of the figures outputted from the main scripts above. It helps to save figures across frequency bands and lateralization. diff --git a/Visualize_Results/scripts/helper_fun/ReadMe.md b/Visualize_Results/scripts/helper_fun/ReadMe.md new file mode 100644 index 0000000..33ddece --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/ReadMe.md @@ -0,0 +1,18 @@ +- The codes inside fdr_bh.m, Violinplot-Matlab-master, and circularGraph folders belong to their respective owners: + - fdr_bh.m: https://www.mathworks.com/matlabcentral/fileexchange/27418-fdr_bh + - Violinplot-Matlab-master: https://github.com/bastibe/Violinplot-Matlab + - circularGraph: https://it.mathworks.com/matlabcentral/fileexchange/48576-circulargraph/ +- Additionally, the modified tutorial_nwa_connectivityviewer.m and tutorial_nwa_connectivityviewer.m files are from the Fieldtrip library. +- Finally, cmap_rbw.m file is taken from the Brainstorm library. + +### Functions + +- load_path: To load the paths. +- cmap_rbw: Color scale for the contrast map. +- tutorial_nwa_connectivityviewer: For 3d visualization. +- fdr_bh: For FDR correction. +- circularGraph folder: Toolbox for circular graphs. +- helper_effectiveConn_fsaverage: Helper function for visualizing on fsaverage; for effective connectivity (pdc). +- helper_functionConn_fsaverage: Helper function for visualizing on fsaverage; for functional connectivity (amplitude- or phase-based). +- helper_functionConn_fsaverage_freqCollapsed: Helper function for collapsing across freqs; for functional connectivity (amplitude- or phase-based). +- helper_functionConn_circularGraph: Helper function for circularGraph; for functional conn and ROI analyses. diff --git a/Visualize_Results/scripts/helper_fun/Violinplot-Matlab-master/README.md b/Visualize_Results/scripts/helper_fun/Violinplot-Matlab-master/README.md new file mode 100644 index 0000000..2ff14d1 --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/Violinplot-Matlab-master/README.md @@ -0,0 +1,41 @@ +# Violin Plots for Matlab + +A violin plot is an easy to read substitute for a box plot that +replaces the box shape with a kernel density estimate of the data, and +optionally overlays the data points itself. The original boxplot shape +is still included as a grey box/line in the center of the violin. + +Violin plots are a superset of box plots, and give a much richer +understanding of the data distribution, while not taking more space. +You will be able to instantly spot too-sparse data, or multi-modal +distributions, which could go unnoticed in boxplots. + +`violinplot` is meant as a direct substitute for `boxplot` (excluding +named arguments). Additional constructor parameters include the width +of the plot, the bandwidth of the kernel density estimation, and the +X-axis position of the violin plot. + +For more information about violin plots, read "[Violin plots: a box +plot-density trace synergism](http://www.stat.cmu.edu/~rnugent/PCMI2016/papers/ViolinPlots.pdf)" +by J. L. Hintze and R. D. Nelson in The American Statistician, vol. +52, no. 2, pp. 181-184, 1998 (DOI: 10.2307/2685478). + +```matlab +load carbig MPG Origin +Origin = cellstr(Origin); +figure +vs = violinplot(MPG, Origin); +ylabel('Fuel Economy in MPG'); +xlim([0.5, 7.5]); +``` + +![example image](example.png) + +## Citation + +[![DOI](https://zenodo.org/badge/60771923.svg)](https://zenodo.org/badge/latestdoi/60771923) + +If you want to cite this repository, use + +> Bechtold, Bastian, 2016. Violin Plots for Matlab, Github Project +> https://github.com/bastibe/Violinplot-Matlab, DOI: 10.5281/zenodo.4559847 diff --git a/Visualize_Results/scripts/helper_fun/Violinplot-Matlab-master/Violin.m b/Visualize_Results/scripts/helper_fun/Violinplot-Matlab-master/Violin.m new file mode 100644 index 0000000..e1c74ca --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/Violinplot-Matlab-master/Violin.m @@ -0,0 +1,346 @@ +classdef Violin < handle + % Violin creates violin plots for some data + % A violin plot is an easy to read substitute for a box plot + % that replaces the box shape with a kernel density estimate of + % the data, and optionally overlays the data points itself. + % + % Additional constructor parameters include the width of the + % plot, the bandwidth of the kernel density estimation, and the + % X-axis position of the violin plot. + % + % Use violinplot for a + % boxplot-like wrapper for + % interactive plotting. + % + % See for more information on Violin Plots: + % J. L. Hintze and R. D. Nelson, "Violin plots: a box + % plot-density trace synergism," The American Statistician, vol. + % 52, no. 2, pp. 181-184, 1998. + % + % Violin Properties: + % ViolinColor - Fill color of the violin area and data points. + % Defaults to the next default color cycle. + % ViolinAlpha - Transparency of the ciolin area and data points. + % Defaults to 0.3. + % EdgeColor - Color of the violin area outline. + % Defaults to [0.5 0.5 0.5] + % BoxColor - Color of the box, whiskers, and the outlines of + % the median point and the notch indicators. + % Defaults to [0.5 0.5 0.5] + % MedianColor - Fill color of the median and notch indicators. + % Defaults to [1 1 1] + % ShowData - Whether to show data points. + % Defaults to true + % ShowNotches - Whether to show notch indicators. + % Defaults to false + % ShowMean - Whether to show mean indicator. + % Defaults to false + % + % Violin Children: + % ScatterPlot - scatter plot of the data points + % ViolinPlot - fill plot of the kernel density estimate + % BoxPlot - fill plot of the box between the quartiles + % WhiskerPlot - line plot between the whisker ends + % MedianPlot - scatter plot of the median (one point) + % NotchPlots - scatter plots for the notch indicators + % MeanPlot - line plot at mean value + + % Copyright (c) 2016, Bastian Bechtold + % This code is released under the terms of the BSD 3-clause license + + properties + ScatterPlot % scatter plot of the data points + ViolinPlot % fill plot of the kernel density estimate + BoxPlot % fill plot of the box between the quartiles + WhiskerPlot % line plot between the whisker ends + MedianPlot % scatter plot of the median (one point) + NotchPlots % scatter plots for the notch indicators + MeanPlot % line plot of the mean (horizontal line) + end + + properties (Dependent=true) + ViolinColor % fill color of the violin area and data points + ViolinAlpha % transparency of the violin area and data points + EdgeColor % color of the violin area outline + BoxColor % color of box, whiskers, and median/notch edges + BoxWidth % width of box between the quartiles in axis space (default 10% of Violin plot width, 0.03) + MedianColor % fill color of median and notches + ShowData % whether to show data points + ShowNotches % whether to show notch indicators + ShowMean % whether to show mean indicator + end + + methods + function obj = Violin(data, pos, varargin) + %Violin plots a violin plot of some data at pos + % VIOLIN(DATA, POS) plots a violin at x-position POS for + % a vector of DATA points. + % + % VIOLIN(..., 'PARAM1', val1, 'PARAM2', val2, ...) + % specifies optional name/value pairs: + % 'Width' Width of the violin in axis space. + % Defaults to 0.3 + % 'Bandwidth' Bandwidth of the kernel density + % estimate. Should be between 10% and + % 40% of the data range. + % 'ViolinColor' Fill color of the violin area and + % data points. Defaults to the next + % default color cycle. + % 'ViolinAlpha' Transparency of the violin area and + % data points. Defaults to 0.3. + % 'EdgeColor' Color of the violin area outline. + % Defaults to [0.5 0.5 0.5] + % 'BoxColor' Color of the box, whiskers, and the + % outlines of the median point and the + % notch indicators. Defaults to + % [0.5 0.5 0.5] + % 'MedianColor' Fill color of the median and notch + % indicators. Defaults to [1 1 1] + % 'ShowData' Whether to show data points. + % Defaults to true + % 'ShowNotches' Whether to show notch indicators. + % Defaults to false + % 'ShowMean' Whether to show mean indicator. + % Defaults to false + + args = obj.checkInputs(data, pos, varargin{:}); + data = data(not(isnan(data))); + if numel(data) == 1 + obj.MedianPlot = scatter(pos, data, 'filled'); + obj.MedianColor = args.MedianColor; + obj.MedianPlot.MarkerEdgeColor = args.EdgeColor; + return + end + + hold('on'); + + % calculate kernel density estimation for the violin + if isempty(data) + return + end + [density, value] = ksdensity(data, 'bandwidth', args.Bandwidth); + density = density(value >= min(data) & value <= max(data)); + value = value(value >= min(data) & value <= max(data)); + value(1) = min(data); + value(end) = max(data); + + % all data is identical + if min(data) == max(data) + density = 1; + end + + width = args.Width/max(density); + + % plot the data points within the violin area + if length(density) > 1 + jitterstrength = interp1(value, density*width, data); + else % all data is identical: + jitterstrength = density*width; + end + jitter = 2*(rand(size(data))-0.5); + obj.ScatterPlot = ... + scatter(pos + jitter.*jitterstrength, data, 'filled'); + + % plot the violin + obj.ViolinPlot = ... % plot color will be overwritten later + fill([pos+density*width pos-density(end:-1:1)*width], ... + [value value(end:-1:1)], [1 1 1]); + + % plot the mini-boxplot within the violin + quartiles = quantile(data, [0.25, 0.5, 0.75]); + obj.BoxPlot = ... % plot color will be overwritten later + fill(pos+[-1,1,1,-1]*args.BoxWidth, ... + [quartiles(1) quartiles(1) quartiles(3) quartiles(3)], ... + [1 1 1]); + + % plot the data mean + meanValue = mean(data); + if length(density) > 1 + meanDensityWidth = interp1(value, density, meanValue)*width; + else % all data is identical: + meanDensityWidth = density*width; + end + if meanDensityWidth lowhisker))); + hiwhisker = quartiles(3) + 1.5*IQR; + hiwhisker = min(hiwhisker, max(data(data < hiwhisker))); + if ~isempty(lowhisker) && ~isempty(hiwhisker) + obj.WhiskerPlot = plot([pos pos], [lowhisker hiwhisker]); + end + obj.MedianPlot = scatter(pos, quartiles(2), [], [1 1 1], 'filled'); + + obj.NotchPlots = ... + scatter(pos, quartiles(2)-1.57*IQR/sqrt(length(data)), ... + [], [1 1 1], 'filled', '^'); + obj.NotchPlots(2) = ... + scatter(pos, quartiles(2)+1.57*IQR/sqrt(length(data)), ... + [], [1 1 1], 'filled', 'v'); + + obj.EdgeColor = args.EdgeColor; + obj.BoxColor = args.BoxColor; + obj.BoxWidth = args.BoxWidth; + obj.MedianColor = args.MedianColor; + if not(isempty(args.ViolinColor)) + obj.ViolinColor = args.ViolinColor; + else + obj.ViolinColor = obj.ScatterPlot.CData; + end + obj.ViolinAlpha = args.ViolinAlpha; + obj.ShowData = args.ShowData; + obj.ShowNotches = args.ShowNotches; + obj.ShowMean = args.ShowMean; + end + + function set.EdgeColor(obj, color) + if ~isempty(obj.ViolinPlot) + obj.ViolinPlot.EdgeColor = color; + end + end + + function color = get.EdgeColor(obj) + if ~isempty(obj.ViolinPlot) + color = obj.ViolinPlot.EdgeColor; + end + end + + function set.MedianColor(obj, color) + obj.MedianPlot.MarkerFaceColor = color; + if ~isempty(obj.NotchPlots) + obj.NotchPlots(1).MarkerFaceColor = color; + obj.NotchPlots(2).MarkerFaceColor = color; + end + end + + function color = get.MedianColor(obj) + color = obj.MedianPlot.MarkerFaceColor; + end + + function set.BoxColor(obj, color) + if ~isempty(obj.BoxPlot) + obj.BoxPlot.FaceColor = color; + obj.BoxPlot.EdgeColor = color; + obj.WhiskerPlot.Color = color; + obj.MedianPlot.MarkerEdgeColor = color; + obj.NotchPlots(1).MarkerFaceColor = color; + obj.NotchPlots(2).MarkerFaceColor = color; + end + end + + function color = get.BoxColor(obj) + if ~isempty(obj.BoxPlot) + color = obj.BoxPlot.FaceColor; + end + end + + function set.BoxWidth(obj,width) + if ~isempty(obj.BoxPlot) + pos=mean(obj.BoxPlot.XData); + obj.BoxPlot.XData=pos+[-1,1,1,-1]*width; + end + end + + function width = get.BoxWidth(obj) + width=max(obj.BoxPlot.XData)-min(obj.BoxPlot.XData); + end + + function set.ViolinColor(obj, color) + obj.ViolinPlot.FaceColor = color; + obj.ScatterPlot.MarkerFaceColor = color; + obj.MeanPlot.Color = color; + end + + function color = get.ViolinColor(obj) + color = obj.ViolinPlot.FaceColor; + end + + function set.ViolinAlpha(obj, alpha) + obj.ScatterPlot.MarkerFaceAlpha = alpha; + obj.ViolinPlot.FaceAlpha = alpha; + end + + function alpha = get.ViolinAlpha(obj) + alpha = obj.ViolinPlot.FaceAlpha; + end + + function set.ShowData(obj, yesno) + if yesno + obj.ScatterPlot.Visible = 'on'; + else + obj.ScatterPlot.Visible = 'off'; + end + end + + function yesno = get.ShowData(obj) + if ~isempty(obj.ScatterPlot) + yesno = strcmp(obj.ScatterPlot.Visible, 'on'); + end + end + + function set.ShowNotches(obj, yesno) + if ~isempty(obj.NotchPlots) + if yesno + obj.NotchPlots(1).Visible = 'on'; + obj.NotchPlots(2).Visible = 'on'; + else + obj.NotchPlots(1).Visible = 'off'; + obj.NotchPlots(2).Visible = 'off'; + end + end + end + + function yesno = get.ShowNotches(obj) + if ~isempty(obj.NotchPlots) + yesno = strcmp(obj.NotchPlots(1).Visible, 'on'); + end + end + + function set.ShowMean(obj, yesno) + if ~isempty(obj.MeanPlot) + if yesno + obj.MeanPlot.Visible = 'on'; + else + obj.MeanPlot.Visible = 'off'; + end + end + end + + function yesno = get.ShowMean(obj) + if ~isempty(obj.MeanPlot) + yesno = strcmp(obj.MeanPlot.Visible, 'on'); + end + end + end + + methods (Access=private) + function results = checkInputs(obj, data, pos, varargin) + isscalarnumber = @(x) (isnumeric(x) & isscalar(x)); + p = inputParser(); + p.addRequired('Data', @isnumeric); + p.addRequired('Pos', isscalarnumber); + p.addParameter('Width', 0.3, isscalarnumber); + p.addParameter('Bandwidth', [], isscalarnumber); + iscolor = @(x) (isnumeric(x) & length(x) == 3); + p.addParameter('ViolinColor', [], iscolor); + p.addParameter('BoxColor', [0.5 0.5 0.5], iscolor); + p.addParameter('BoxWidth', 0.01, isscalarnumber); + p.addParameter('EdgeColor', [0.5 0.5 0.5], iscolor); + p.addParameter('MedianColor', [1 1 1], iscolor); + p.addParameter('ViolinAlpha', 0.3, isscalarnumber); + isscalarlogical = @(x) (islogical(x) & isscalar(x)); + p.addParameter('ShowData', true, isscalarlogical); + p.addParameter('ShowNotches', false, isscalarlogical); + p.addParameter('ShowMean', false, isscalarlogical); + + p.parse(data, pos, varargin{:}); + results = p.Results; + end + end +end diff --git a/Visualize_Results/scripts/helper_fun/Violinplot-Matlab-master/example.png b/Visualize_Results/scripts/helper_fun/Violinplot-Matlab-master/example.png new file mode 100644 index 0000000..240e509 Binary files /dev/null and b/Visualize_Results/scripts/helper_fun/Violinplot-Matlab-master/example.png differ diff --git a/Visualize_Results/scripts/helper_fun/Violinplot-Matlab-master/testviolinplot.m b/Visualize_Results/scripts/helper_fun/Violinplot-Matlab-master/testviolinplot.m new file mode 100644 index 0000000..6047002 --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/Violinplot-Matlab-master/testviolinplot.m @@ -0,0 +1,21 @@ + +% TEST CASE 1 +display('Test 1: Violin plot default options'); +load carbig MPG Origin +Origin = cellstr(Origin); +figure +vs = violinplot(MPG, Origin); +ylabel('Fuel Economy in MPG'); +xlim([0.5, 7.5]); + +display('Test 1 passed ok'); + +% TEST CASE 2 +display('Test 2: Test the plot ordering option'); +grouporder={'USA','Sweden','Japan','Italy','Germany','France','England'}; + +figure; +vs2 = violinplot(MPG,Origin,'GroupOrder',grouporder); +display('Test 2 passed ok'); + +%other test cases could be added here diff --git a/Visualize_Results/scripts/helper_fun/Violinplot-Matlab-master/violinplot.m b/Visualize_Results/scripts/helper_fun/Violinplot-Matlab-master/violinplot.m new file mode 100644 index 0000000..57264f8 --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/Violinplot-Matlab-master/violinplot.m @@ -0,0 +1,127 @@ +function violins = violinplot(data, cats, varargin) +%Violinplots plots violin plots of some data and categories +% VIOLINPLOT(DATA) plots a violin of a double vector DATA +% +% VIOLINPLOT(DATAMATRIX) plots violins for each column in +% DATAMATRIX. +% +% VIOLINPLOT(TABLE), VIOLINPLOT(STRUCT), VIOLINPLOT(DATASET) +% plots violins for each column in TABLE, each field in STRUCT, and +% each variable in DATASET. The violins are labeled according to +% the table/dataset variable name or the struct field name. +% +% VIOLINPLOT(DATAMATRIX, CATEGORYNAMES) plots violins for each +% column in DATAMATRIX and labels them according to the names in the +% cell-of-strings CATEGORYNAMES. +% +% VIOLINPLOT(DATA, CATEGORIES) where double vector DATA and vector +% CATEGORIES are of equal length; plots violins for each category in +% DATA. +% +% violins = VIOLINPLOT(...) returns an object array of +% Violin objects. +% +% VIOLINPLOT(..., 'PARAM1', val1, 'PARAM2', val2, ...) +% specifies optional name/value pairs for all violins: +% 'Width' Width of the violin in axis space. +% Defaults to 0.3 +% 'Bandwidth' Bandwidth of the kernel density estimate. +% Should be between 10% and 40% of the data range. +% 'ViolinColor' Fill color of the violin area and data points. +% Defaults to the next default color cycle. +% 'ViolinAlpha' Transparency of the violin area and data points. +% Defaults to 0.3. +% 'EdgeColor' Color of the violin area outline. +% Defaults to [0.5 0.5 0.5] +% 'BoxColor' Color of the box, whiskers, and the outlines of +% the median point and the notch indicators. +% Defaults to [0.5 0.5 0.5] +% 'MedianColor' Fill color of the median and notch indicators. +% Defaults to [1 1 1] +% 'ShowData' Whether to show data points. +% Defaults to true +% 'ShowNotches' Whether to show notch indicators. +% Defaults to false +% 'ShowMean' Whether to show mean indicator +% Defaults to false +% 'GroupOrder' Cell of category names in order to be plotted. +% Defaults to alphabetical ordering + +% Copyright (c) 2016, Bastian Bechtold +% This code is released under the terms of the BSD 3-clause license + + hascategories = exist('cats','var') && not(isempty(cats)); + + %parse the optional grouporder argument + %if it exists parse the categories order + % but also delete it from the arguments passed to Violin + grouporder = {}; + idx=find(strcmp(varargin, 'GroupOrder')); + if ~isempty(idx) && numel(varargin)>idx + if iscell(varargin{idx+1}) + grouporder = varargin{idx+1}; + varargin(idx:idx+1)=[]; + else + error('Second argument of ''GroupOrder'' optional arg must be a cell of category names') + end + end + + % tabular data + if isa(data, 'dataset') || isstruct(data) || istable(data) + if isa(data, 'dataset') + colnames = data.Properties.VarNames; + elseif istable(data) + colnames = data.Properties.VariableNames; + elseif isstruct(data) + colnames = fieldnames(data); + end + catnames = {}; + for n=1:length(colnames) + if isnumeric(data.(colnames{n})) + catnames = [catnames colnames{n}]; + end + end + for n=1:length(catnames) + thisData = data.(catnames{n}); + violins(n) = Violin(thisData, n, varargin{:}); + end + set(gca, 'XTick', 1:length(catnames), 'XTickLabels', catnames); + + % 1D data, one category for each data point + elseif hascategories && numel(data) == numel(cats) + + if isempty(grouporder) + cats = categorical(cats); + else + cats = categorical(cats, grouporder); + end + + catnames = (unique(cats)); % this ignores categories without any data + catnames_labels = {}; + for n = 1:length(catnames) + thisCat = catnames(n); + catnames_labels{n} = char(thisCat); + thisData = data(cats == thisCat); + violins(n) = Violin(thisData, n, varargin{:}); + end + set(gca, 'XTick', 1:length(catnames), 'XTickLabels', catnames_labels); + + % 1D data, no categories + elseif not(hascategories) && isvector(data) + violins = Violin(data, 1, varargin{:}); + set(gca, 'XTick', 1); + + % 2D data with or without categories + elseif ismatrix(data) + for n=1:size(data, 2) + thisData = data(:, n); + violins(n) = Violin(thisData, n, varargin{:}); + end + set(gca, 'XTick', 1:size(data, 2)); + if hascategories && length(cats) == size(data, 2) + set(gca, 'XTickLabels', cats); + end + + end + +end diff --git a/Visualize_Results/scripts/helper_fun/circularGraph/README.md b/Visualize_Results/scripts/helper_fun/circularGraph/README.md new file mode 100644 index 0000000..a548a67 --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/circularGraph/README.md @@ -0,0 +1,10 @@ +# circularGraph + +## Description +A **circular graph** is a visualization of a network of nodes and their connections. The nodes are laid out along a circle, and the connections are drawn within the circle. Click on a node to make the connections that emanate from it more visible or less visible. Click on the **Show All** button to make all nodes and their connections visible. Click on the **Hide All** button to make all nodes and their connections less visible. + +## Required Products +* MATLAB 8.4 (R2014b) + +## Tags +adjacency, adjacency matrix, circle, circle graph, circle network, circular, circular graph, circular network, connection, connections, connectivity, directed graph, graph, graph theory, network, network plot, node, nodes, undirected graph diff --git a/Visualize_Results/scripts/helper_fun/circularGraph/circularGraph.PNG b/Visualize_Results/scripts/helper_fun/circularGraph/circularGraph.PNG new file mode 100644 index 0000000..01c4209 Binary files /dev/null and b/Visualize_Results/scripts/helper_fun/circularGraph/circularGraph.PNG differ diff --git a/Visualize_Results/scripts/helper_fun/circularGraph/circularGraph.m b/Visualize_Results/scripts/helper_fun/circularGraph/circularGraph.m new file mode 100644 index 0000000..2a39db7 --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/circularGraph/circularGraph.m @@ -0,0 +1,197 @@ +classdef circularGraph < handle +% CIRCULARGRAPH Plot an interactive circular graph to illustrate connections in a network. +% +%% Syntax +% circularGraph(X) +% circularGraph(X,'PropertyName',propertyvalue,...) +% h = circularGraph(...) +% +%% Description +% A 'circular graph' is a visualization of a network of nodes and their +% connections. The nodes are laid out along a circle, and the connections +% are drawn within the circle. Click on a node to make the connections that +% emanate from it more visible or less visible. Click on the 'Show All' +% button to make all nodes and their connections visible. Click on the +% 'Hide All' button to make all nodes and their connections less visible. +% +% Required input arguments. +% X : A symmetric matrix of numeric or logical values. +% +% Optional properties. +% Colormap : A N by 3 matrix of [r g b] triples, where N is the +% length(adjacenyMatrix). +% Label : A cell array of N strings. +%% +% Copyright 2016 The MathWorks, Inc. + properties + Node = node(0,0); % Array of nodes + ColorMap; % Colormap + Label; % Cell array of strings + ShowButton; % Turn all nodes on + HideButton; % Turn all nodes off + end + + methods + function this = circularGraph(adjacencyMatrix,varargin) + % Constructor + p = inputParser; + + defaultColorMap = parula(length(adjacencyMatrix)); + defaultLabel = cell(length(adjacencyMatrix)); + for i = 1:length(defaultLabel) + defaultLabel{i} = num2str(i); + end + + addRequired(p,'adjacencyMatrix',@(x)(isnumeric(x) || islogical(x))); + addParameter(p,'ColorMap',defaultColorMap,@(colormap)length(colormap) == length(adjacencyMatrix)); + addParameter(p,'Label' ,defaultLabel ,@iscell); + + parse(p,adjacencyMatrix,varargin{:}); + this.ColorMap = p.Results.ColorMap; + this.Label = p.Results.Label; + + this.ShowButton = uicontrol(... + 'Style','pushbutton',... + 'Position',[0 40 80 40],... + 'String','Show All',... + 'Callback',@circularGraph.showNodes,... + 'UserData',this); + + this.HideButton = uicontrol(... + 'Style','pushbutton',... + 'Position',[0 0 80 40],... + 'String','Hide All',... + 'Callback',@circularGraph.hideNodes,... + 'UserData',this); + + fig = gcf; + set(fig,... + 'UserData',this,... + 'CloseRequestFcn',@circularGraph.CloseRequestFcn); + + % Draw the nodes + delete(this.Node); + t = linspace(-pi,pi,length(adjacencyMatrix) + 1).'; % theta for each node + extent = zeros(length(adjacencyMatrix),1); + for i = 1:length(adjacencyMatrix) + this.Node(i) = node(cos(t(i)),sin(t(i))); + this.Node(i).Color = this.ColorMap(i,:); + this.Node(i).Label = this.Label{i}; + end + + % Find non-zero values of s and their indices + [row,col,v] = find(adjacencyMatrix); + + % Calculate line widths based on values of s (stored in v). + minLineWidth = 0.1; + lineWidthCoef = 1.5; +% Orhan: + lineWidth = v.^2; %v./max(v); + if sum(lineWidth) == numel(lineWidth) % all lines are the same width. + lineWidth = repmat(minLineWidth,numel(lineWidth),1); + else % lines of variable width. + lineWidth = lineWidthCoef*lineWidth + minLineWidth; + end + + % Draw connections on the Poincare hyperbolic disk. + % + % Equation of the circles on the disk: + % x^2 + y^2 + % + 2*(u(2)-v(2))/(u(1)*v(2)-u(2)*v(1))*x + % - 2*(u(1)-v(1))/(u(1)*v(2)-u(2)*v(1))*y + 1 = 0, + % where u and v are points on the boundary. + % + % Standard form of equation of a circle + % (x - x0)^2 + (y - y0)^2 = r^2 + % + % Therefore we can identify + % x0 = -(u(2)-v(2))/(u(1)*v(2)-u(2)*v(1)); + % y0 = (u(1)-v(1))/(u(1)*v(2)-u(2)*v(1)); + % r^2 = x0^2 + y0^2 - 1 + + for i = 1:length(v) + if row(i) ~= col(i) + if abs(row(i) - col(i)) - length(adjacencyMatrix)/2 == 0 + % points are diametric, so draw a straight line + u = [cos(t(row(i)));sin(t(row(i)))]; + v = [cos(t(col(i)));sin(t(col(i)))]; + this.Node(row(i)).Connection(end+1) = line(... + [u(1);v(1)],... + [u(2);v(2)],... + 'LineWidth', lineWidth(i),... + 'Color', this.ColorMap(row(i),:),... + 'PickableParts','none'); + else % points are not diametric, so draw an arc + u = [cos(t(row(i)));sin(t(row(i)))]; + v = [cos(t(col(i)));sin(t(col(i)))]; + x0 = -(u(2)-v(2))/(u(1)*v(2)-u(2)*v(1)); + y0 = (u(1)-v(1))/(u(1)*v(2)-u(2)*v(1)); + r = sqrt(x0^2 + y0^2 - 1); + thetaLim(1) = atan2(u(2)-y0,u(1)-x0); + thetaLim(2) = atan2(v(2)-y0,v(1)-x0); + + if u(1) >= 0 && v(1) >= 0 + % ensure the arc is within the unit disk + theta = [linspace(max(thetaLim),pi,50),... + linspace(-pi,min(thetaLim),50)].'; + else + theta = linspace(thetaLim(1),thetaLim(2)).'; + end + + this.Node(row(i)).Connection(end+1) = line(... + r*cos(theta)+x0,... + r*sin(theta)+y0,... + 'LineWidth', lineWidth(i),... + 'Color', this.ColorMap(row(i),:),... + 'PickableParts','none'); + end + end + end + + axis image; + ax = gca; + for i = 1:length(adjacencyMatrix) + extent(i) = this.Node(i).Extent; + end + extent = max(extent(:)); + ax.XLim = ax.XLim + extent*[-1 1]; + fudgeFactor = 1.75; % Not sure why this is necessary. Eyeballed it. + ax.YLim = ax.YLim + fudgeFactor*extent*[-1 1]; + ax.Visible = 'off'; + ax.SortMethod = 'depth'; + + fig = gcf; + fig.Color = [1 1 1]; + end + + end + + methods (Static = true) + function showNodes(this,~) + % Callback for 'Show All' button + n = this.UserData.Node; + for i = 1:length(n) + n(i).Visible = true; + end + end + + function hideNodes(this,~) + % Callback for 'Hide All' button + n = this.UserData.Node; + for i = 1:length(n) + n(i).Visible = false; + end + end + + function CloseRequestFcn(this,~) + % Callback for figure CloseRequestFcn + c = this.UserData; + for i = 1:length(c.Node) + delete(c.Node(i)); + end + delete(gcf); + end + + end + +end \ No newline at end of file diff --git a/Visualize_Results/scripts/helper_fun/circularGraph/circularGraph.mltbx b/Visualize_Results/scripts/helper_fun/circularGraph/circularGraph.mltbx new file mode 100644 index 0000000..be649b6 Binary files /dev/null and b/Visualize_Results/scripts/helper_fun/circularGraph/circularGraph.mltbx differ diff --git a/Visualize_Results/scripts/helper_fun/circularGraph/circularGraph.prj b/Visualize_Results/scripts/helper_fun/circularGraph/circularGraph.prj new file mode 100644 index 0000000..f910a9a --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/circularGraph/circularGraph.prj @@ -0,0 +1,131 @@ + + + circularGraph + Paul Kassebaum + paul.kassebaum@mathworks.com + MathWorks + Plot an interactive circular graph to illustrate connections in a network. + A 'circular graph' is a visualization of a network of nodes and their connections. The nodes are laid out along a circle, and the connections are drawn within the circle. Click on a node to make the connections that emanate from it more visible or less visible. Click on the 'Show All' button to make all nodes and their connections visible. Click on the 'Hide All' button to make all nodes and their connections less visible. + ${PROJECT_ROOT}\circularGraph.PNG + 2.0 + ${PROJECT_ROOT}\circularGraph.mltbx + + + + + 04027212-1e6d-49d4-84fa-f66cc29418de + % List files contained in your toolbox folder that you would like to exclude +% from packaging. Excludes should be listed relative to the toolbox folder. +% Some examples of how to specify excludes are provided below: +% +% A single file in the toolbox folder: +% .svn +% +% A single file in a subfolder of the toolbox folder: +% example/.svn +% +% All files in a subfolder of the toolbox folder: +% example/* +% +% All files of a certain name in all subfolders of the toolbox folder: +% **/.svn +% +% All files matching a pattern in all subfolders of the toolbox folder: +% **/*.bak +% +.git + <?xml version="1.0" encoding="utf-8"?> +<examples> + <exampleCategory name="circularGraph"> + <example name="example" type="html"> + <file type="source">/html/example.html</file> + <file type="main">/example.m</file> + <file type="thumbnail">/html/example.png</file> + <file type="image">/html/example_01.png</file> + <file type="image">/html/example_02.png</file> + </example> + </exampleCategory> +</examples> + + + + + + + + + + + + + + + + + + + C:\Users\pkasseba\Documents\MATLAB\circularGraph + + + ${PROJECT_ROOT}\circularGraph.m + ${PROJECT_ROOT}\circularGraph.PNG + ${PROJECT_ROOT}\demos.xml + ${PROJECT_ROOT}\example.m + ${PROJECT_ROOT}\html + ${PROJECT_ROOT}\license.txt + ${PROJECT_ROOT}\node.m + ${PROJECT_ROOT}\README.md + + + + + + C:\Users\pkasseba\Documents\MATLAB\circularGraph\circularGraph.mltbx + + + + C:\Program Files\MATLAB\R2016b + + + + + + + + + true + + + + + true + + + + + true + + + + + true + + + + + false + false + true + false + false + false + false + false + 6.2 + false + true + win64 + true + + + \ No newline at end of file diff --git a/Visualize_Results/scripts/helper_fun/circularGraph/demos.xml b/Visualize_Results/scripts/helper_fun/circularGraph/demos.xml new file mode 100644 index 0000000..497e4d6 --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/circularGraph/demos.xml @@ -0,0 +1,17 @@ + + + circularGraph + toolbox + HelpIcon.DEMOS + + A 'circular graph' is a visualization of a network of nodes and their connections. The nodes are laid out along a circle, and the connections are drawn within the circle. Click on a node to make the connections that emanate from it more visible or less visible. Click on the 'Show All' button to make all nodes and their connections visible. Click on the 'Hide All' button to make all nodes and their connections less visible. + + + + + other + example + html/example.html + + + \ No newline at end of file diff --git a/Visualize_Results/scripts/helper_fun/circularGraph/example.m b/Visualize_Results/scripts/helper_fun/circularGraph/example.m new file mode 100644 index 0000000..2ee55c0 --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/circularGraph/example.m @@ -0,0 +1,61 @@ +%% Circular Graph Examples +% Copyright 2016 The MathWorks, Inc. + +%% 1. Adjacency matrix of 1s and 0s +% Create an example adjacency matrix made up of ones and zeros. +rng(0); +x = rand(50); +thresh = 0.93; +x(x > thresh) = 1; +x(x <= thresh) = 0; + +%% +% Call CIRCULARGRAPH with only the adjacency matrix as an argument. +circularGraph(x); + +%% +% Click on a node to make the connections that emanate from it more visible +% or less visible. Click on the 'Show All' button to make all nodes and +% their connections visible. Click on the 'Hide All' button to make all +% nodes and their connections less visible. + +%% 2. Supply custom properties +% Create an example adjacency matrix made up of various values and supply +% custom properties. +rng(0); +x = rand(20); +thresh = 0.93; +x(x > thresh) = 1; +x(x <= thresh) = 0; +for i = 1:numel(x) + if x(i) > 0 + x(i) = rand(1,1); + end +end + +%% +% Create custom node labels +myLabel = cell(length(x)); +for i = 1:length(x) + myLabel{i} = num2str(round(1000000*rand(1,1))); +end + +%% +% Create custom colormap +figure; +myColorMap = lines(length(x)); + +circularGraph(x,'Colormap',myColorMap,'Label',myLabel); + +fg1 = gca; +fg2 = get(fg1, 'children'); +fg2(end-2).FontSize = 15 + +for ii = 1:length(fg2) + try + tmp_fg = fg2(ii); + tmp_fg.FontSize = 15; + catch + fprintf('No field.\n', i); + end +end diff --git a/Visualize_Results/scripts/helper_fun/circularGraph/html/example.html b/Visualize_Results/scripts/helper_fun/circularGraph/html/example.html new file mode 100644 index 0000000..253d05e --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/circularGraph/html/example.html @@ -0,0 +1,144 @@ + + + + + Circular Graph Examples

Circular Graph Examples

Contents

1. Adjacency matrix of 1s and 0s

Create an example adjacency matrix made up of ones and zeros.

rng(0);
+x = rand(50);
+thresh = 0.93;
+x(x >  thresh) = 1;
+x(x <= thresh) = 0;
+

Call CIRCULARGRAPH with only the adjacency matrix as an argument.

circularGraph(x);
+

Click on a node to make the connections that emanate from it more visible or less visible. Click on the 'Show All' button to make all nodes and their connections visible. Click on the 'Hide All' button to make all nodes and their connections less visible.

2. Supply custom properties

Create an example adjacency matrix made up of various values and supply custom properties.

rng(0);
+x = rand(20);
+thresh = 0.93;
+x(x >  thresh) = 1;
+x(x <= thresh) = 0;
+for i = 1:numel(x)
+  if x(i) > 0
+    x(i) = rand(1,1);
+  end
+end
+

Create custom node labels

myLabel = cell(length(x));
+for i = 1:length(x)
+  myLabel{i} = num2str(round(1000000*rand(1,1)));
+end
+

Create custom colormap

figure;
+myColorMap = lines(length(x));
+
+circularGraph(x,'Colormap',myColorMap,'Label',myLabel);
+
\ No newline at end of file diff --git a/Visualize_Results/scripts/helper_fun/circularGraph/html/example.png b/Visualize_Results/scripts/helper_fun/circularGraph/html/example.png new file mode 100644 index 0000000..67ebd07 Binary files /dev/null and b/Visualize_Results/scripts/helper_fun/circularGraph/html/example.png differ diff --git a/Visualize_Results/scripts/helper_fun/circularGraph/html/example_01.png b/Visualize_Results/scripts/helper_fun/circularGraph/html/example_01.png new file mode 100644 index 0000000..d03b7ce Binary files /dev/null and b/Visualize_Results/scripts/helper_fun/circularGraph/html/example_01.png differ diff --git a/Visualize_Results/scripts/helper_fun/circularGraph/html/example_02.png b/Visualize_Results/scripts/helper_fun/circularGraph/html/example_02.png new file mode 100644 index 0000000..c3732c4 Binary files /dev/null and b/Visualize_Results/scripts/helper_fun/circularGraph/html/example_02.png differ diff --git a/Visualize_Results/scripts/helper_fun/circularGraph/license.txt b/Visualize_Results/scripts/helper_fun/circularGraph/license.txt new file mode 100644 index 0000000..599af41 --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/circularGraph/license.txt @@ -0,0 +1,27 @@ +Copyright (c) 2016, The MathWorks, Inc. +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the distribution + * Neither the name of the The MathWorks, Inc. nor the names + of its contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. diff --git a/Visualize_Results/scripts/helper_fun/circularGraph/node.m b/Visualize_Results/scripts/helper_fun/circularGraph/node.m new file mode 100644 index 0000000..f362b2e --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/circularGraph/node.m @@ -0,0 +1,127 @@ +classdef node < handle +% NODE Helper class for circularGraph. Not intended for direct user manipulation. +%% +% Copyright 2016 The MathWorks, Inc. + properties (Access = public) + Label = ''; % String + Connection = line(0,0); % Array of lines + Position; % [x,y] coordinates + Color = [0 0 0]; % [r g b] + Visible = true; % Logical true or false + end + + properties (Access = public, Dependent = true) + Extent; % Width of text label + end + + properties (Access = private) + TextLabel; % Text graphics object + NodeMarker; % Line that makes the node visible + Marker = 'o'; % Marker symbol when the node is 'on' + end + + properties (Access = private, Constant) + labelOffsetFactor = 1.1; + end + + methods + function this = node(x,y) + % Constructor + this.Position = [x,y]; + this.Connection = line(0,0); + makeLine(this); + end + + function makeLine(this) + % Make the node's line graphic object + this.NodeMarker = line(... + this.Position(1),... + this.Position(2),... + 2,... + 'Color',this.Color,... + 'Marker',this.Marker,... + 'LineStyle','none',... + 'PickableParts','all',... + 'ButtonDownFcn',@node.ButtonDownFcn,... + 'UserData',this); + end + + function set.Visible(this,value) + this.Visible = value; + updateVisible(this); + end + + function set.Color(this,value) + this.Color = value; + updateColor(this); + end + + function set.Label(this,value) + this.Label = value; + updateTextLabel(this); + end + + function value = get.Extent(this) + value = this.TextLabel.Extent(3); + end + + function updateVisible(this) + if this.Visible + this.NodeMarker.Marker = 'o'; + set(this.Connection,'Color',this.Color); + for i = 1:length(this.Connection) + this.Connection(i).ZData = ones(size(this.Connection(i).XData)); + end + else + this.NodeMarker.Marker = 'x'; + set(this.Connection,'Color',0.9*[1 1 1]); + for i = 1:length(this.Connection) + this.Connection(i).ZData = zeros(size(this.Connection(i).XData)); + end + end + end + + function updateColor(this) + this.NodeMarker.Color = this.Color; + set(this.Connection,'Color',this.Color); + end + + function updateTextLabel(this) + delete(this.TextLabel); + + x = this.Position(1); + y = this.Position(2); + t = atan2(y,x); + + this.TextLabel = text(0,0,this.Label); + + this.TextLabel.Position = node.labelOffsetFactor*this.Position; + if abs(t) > pi/2 + this.TextLabel.Rotation = 180*(t/pi + 1); + this.TextLabel.HorizontalAlignment = 'right'; + else + this.TextLabel.Rotation = t*180/pi; + end + end + + function delete(this) + % Destructor + delete(this.Connection(:)) + delete(this.TextLabel); + delete(this.NodeMarker); + delete(this); + end + + end + + methods (Static = true) + function ButtonDownFcn(this,~) + n = this.UserData; + if n.Visible + n.Visible = false; + else + n.Visible = true; + end + end + end +end \ No newline at end of file diff --git a/Visualize_Results/scripts/helper_fun/cmap_rbw.m b/Visualize_Results/scripts/helper_fun/cmap_rbw.m new file mode 100644 index 0000000..1168afe --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/cmap_rbw.m @@ -0,0 +1,42 @@ +function CMap = cmap_rbw(N) +% CMAP_RBW: Red-blue-white colormap. + +% @============================================================================= +% This function is part of the Brainstorm software: +% https://neuroimage.usc.edu/brainstorm +% +% Copyright (c)2000-2020 University of Southern California & McGill University +% This software is distributed under the terms of the GNU General Public License +% as published by the Free Software Foundation. Further details on the GPLv3 +% license can be found at http://www.gnu.org/copyleft/gpl.html. +% +% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE +% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY +% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF +% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY +% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. +% +% For more information type "brainstorm license" at command prompt. +% =============================================================================@ +% +% Authors: Sylvain Baillet + +if (nargin < 1) + N = 128; % Number of color levels +end + +Wwidth = 0; % Width of the white zero level +Cmin = 0; +Cmax = 1; + +Half = floor(N/2); + +R1 = linspace(Cmin,Cmax,Half-Wwidth); +R2 = ones(1,Wwidth*2); +R3 = ones(1,Half-Wwidth)*Cmax; +R = [R1 R2 R3]; +B = fliplr(R); +G3 = fliplr(R1); +G = [R1 R2 G3]; +CMap = [R' G' B']; + diff --git a/Visualize_Results/scripts/helper_fun/customMap_effectiveConn.m b/Visualize_Results/scripts/helper_fun/customMap_effectiveConn.m new file mode 100644 index 0000000..6effdd7 --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/customMap_effectiveConn.m @@ -0,0 +1,52 @@ +function CMap = customMap_effectiveConn(N) + +% % Orhan: +% https://jdherman.github.io/colormap/ + +% @============================================================================= +% This function is part of the Brainstorm software: +% https://neuroimage.usc.edu/brainstorm +% +% Copyright (c)2000-2020 University of Southern California & McGill University +% This software is distributed under the terms of the GNU General Public License +% as published by the Free Software Foundation. Further details on the GPLv3 +% license can be found at http://www.gnu.org/copyleft/gpl.html. +% +% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE +% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY +% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF +% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY +% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. +% +% For more information type "brainstorm license" at command prompt. +% =============================================================================@ +% +% Authors: Sylvain Baillet + +if (nargin < 1) + N = 128; % Number of color levels +end + +%% + +Wwidth = 0; % Width of the white zero level +Half = floor(N/2); + +R1 = linspace(0.1569, 1, Half-Wwidth); +R2 = ones(1,Wwidth*2); +R3 = flip(linspace(0, 1, Half-Wwidth)); +R = [R1 R2 R3]'; + +G1 = linspace(0, 1, Half-Wwidth); +G2 = ones(1,Wwidth*2); +G3 = flip(linspace(0.3137, 1, Half-Wwidth)); +G = [G1 G2 G3]'; + +B1 = linspace(0.3137, 1, Half-Wwidth); +B2 = ones(1,Wwidth*2); +B3 = flip(linspace(0.1569, 1, Half-Wwidth)); +B = [B1 B2 B3]'; + +CMap = [R G B]; + + diff --git a/Visualize_Results/scripts/helper_fun/fdr_bh.m b/Visualize_Results/scripts/helper_fun/fdr_bh.m new file mode 100644 index 0000000..ef12d8b --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/fdr_bh.m @@ -0,0 +1,226 @@ +% fdr_bh() - Executes the Benjamini & Hochberg (1995) and the Benjamini & +% Yekutieli (2001) procedure for controlling the false discovery +% rate (FDR) of a family of hypothesis tests. FDR is the expected +% proportion of rejected hypotheses that are mistakenly rejected +% (i.e., the null hypothesis is actually true for those tests). +% FDR is a somewhat less conservative/more powerful method for +% correcting for multiple comparisons than procedures like Bonferroni +% correction that provide strong control of the family-wise +% error rate (i.e., the probability that one or more null +% hypotheses are mistakenly rejected). +% +% This function also returns the false coverage-statement rate +% (FCR)-adjusted selected confidence interval coverage (i.e., +% the coverage needed to construct multiple comparison corrected +% confidence intervals that correspond to the FDR-adjusted p-values). +% +% +% Usage: +% >> [h, crit_p, adj_ci_cvrg, adj_p]=fdr_bh(pvals,q,method,report); +% +% Required Input: +% pvals - A vector or matrix (two dimensions or more) containing the +% p-value of each individual test in a family of tests. +% +% Optional Inputs: +% q - The desired false discovery rate. {default: 0.05} +% method - ['pdep' or 'dep'] If 'pdep,' the original Bejnamini & Hochberg +% FDR procedure is used, which is guaranteed to be accurate if +% the individual tests are independent or positively dependent +% (e.g., Gaussian variables that are positively correlated or +% independent). If 'dep,' the FDR procedure +% described in Benjamini & Yekutieli (2001) that is guaranteed +% to be accurate for any test dependency structure (e.g., +% Gaussian variables with any covariance matrix) is used. 'dep' +% is always appropriate to use but is less powerful than 'pdep.' +% {default: 'pdep'} +% report - ['yes' or 'no'] If 'yes', a brief summary of FDR results are +% output to the MATLAB command line {default: 'no'} +% +% +% Outputs: +% h - A binary vector or matrix of the same size as the input "pvals." +% If the ith element of h is 1, then the test that produced the +% ith p-value in pvals is significant (i.e., the null hypothesis +% of the test is rejected). +% crit_p - All uncorrected p-values less than or equal to crit_p are +% significant (i.e., their null hypotheses are rejected). If +% no p-values are significant, crit_p=0. +% adj_ci_cvrg - The FCR-adjusted BH- or BY-selected +% confidence interval coverage. For any p-values that +% are significant after FDR adjustment, this gives you the +% proportion of coverage (e.g., 0.99) you should use when generating +% confidence intervals for those parameters. In other words, +% this allows you to correct your confidence intervals for +% multiple comparisons. You can NOT obtain confidence intervals +% for non-significant p-values. The adjusted confidence intervals +% guarantee that the expected FCR is less than or equal to q +% if using the appropriate FDR control algorithm for the +% dependency structure of your data (Benjamini & Yekutieli, 2005). +% FCR (i.e., false coverage-statement rate) is the proportion +% of confidence intervals you construct +% that miss the true value of the parameter. adj_ci=NaN if no +% p-values are significant after adjustment. +% adj_p - All adjusted p-values less than or equal to q are significant +% (i.e., their null hypotheses are rejected). Note, adjusted +% p-values can be greater than 1. +% +% +% References: +% Benjamini, Y. & Hochberg, Y. (1995) Controlling the false discovery +% rate: A practical and powerful approach to multiple testing. Journal +% of the Royal Statistical Society, Series B (Methodological). 57(1), +% 289-300. +% +% Benjamini, Y. & Yekutieli, D. (2001) The control of the false discovery +% rate in multiple testing under dependency. The Annals of Statistics. +% 29(4), 1165-1188. +% +% Benjamini, Y., & Yekutieli, D. (2005). False discovery rate?adjusted +% multiple confidence intervals for selected parameters. Journal of the +% American Statistical Association, 100(469), 71?81. doi:10.1198/016214504000001907 +% +% +% Example: +% nullVars=randn(12,15); +% [~, p_null]=ttest(nullVars); %15 tests where the null hypothesis +% %is true +% effectVars=randn(12,5)+1; +% [~, p_effect]=ttest(effectVars); %5 tests where the null +% %hypothesis is false +% [h, crit_p, adj_ci_cvrg, adj_p]=fdr_bh([p_null p_effect],.05,'pdep','yes'); +% data=[nullVars effectVars]; +% fcr_adj_cis=NaN*zeros(2,20); %initialize confidence interval bounds to NaN +% if ~isnan(adj_ci_cvrg), +% sigIds=find(h); +% fcr_adj_cis(:,sigIds)=tCIs(data(:,sigIds),adj_ci_cvrg); % tCIs.m is available on the +% %Mathworks File Exchagne +% end +% +% +% For a review of false discovery rate control and other contemporary +% techniques for correcting for multiple comparisons see: +% +% Groppe, D.M., Urbach, T.P., & Kutas, M. (2011) Mass univariate analysis +% of event-related brain potentials/fields I: A critical tutorial review. +% Psychophysiology, 48(12) pp. 1711-1725, DOI: 10.1111/j.1469-8986.2011.01273.x +% http://www.cogsci.ucsd.edu/~dgroppe/PUBLICATIONS/mass_uni_preprint1.pdf +% +% +% For a review of FCR-adjusted confidence intervals (CIs) and other techniques +% for adjusting CIs for multiple comparisons see: +% +% Groppe, D.M. (in press) Combating the scientific decline effect with +% confidence (intervals). Psychophysiology. +% http://biorxiv.org/content/biorxiv/early/2015/12/10/034074.full.pdf +% +% +% Author: +% David M. Groppe +% Kutaslab +% Dept. of Cognitive Science +% University of California, San Diego +% March 24, 2010 + +%%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%% +% +% 5/7/2010-Added FDR adjusted p-values +% 5/14/2013- D.H.J. Poot, Erasmus MC, improved run-time complexity +% 10/2015- Now returns FCR adjusted confidence intervals + +function [h, crit_p, adj_ci_cvrg, adj_p]=fdr_bh(pvals,q,method,report) + +if nargin<1, + error('You need to provide a vector or matrix of p-values.'); +else + if ~isempty(find(pvals<0,1)), + error('Some p-values are less than 0.'); + elseif ~isempty(find(pvals>1,1)), + error('Some p-values are greater than 1.'); + end +end + +if nargin<2, + q=.05; +end + +if nargin<3, + method='pdep'; +end + +if nargin<4, + report='no'; +end + +s=size(pvals); +if (length(s)>2) || s(1)>1, + [p_sorted, sort_ids]=sort(reshape(pvals,1,prod(s))); +else + %p-values are already a row vector + [p_sorted, sort_ids]=sort(pvals); +end +[dummy, unsort_ids]=sort(sort_ids); %indexes to return p_sorted to pvals order +m=length(p_sorted); %number of tests + +if strcmpi(method,'pdep'), + %BH procedure for independence or positive dependence + thresh=(1:m)*q/m; + wtd_p=m*p_sorted./(1:m); + +elseif strcmpi(method,'dep') + %BH procedure for any dependency structure + denom=m*sum(1./(1:m)); + thresh=(1:m)*q/denom; + wtd_p=denom*p_sorted./[1:m]; + %Note, it can produce adjusted p-values greater than 1! + %compute adjusted p-values +else + error('Argument ''method'' needs to be ''pdep'' or ''dep''.'); +end + +if nargout>3, + %compute adjusted p-values; This can be a bit computationally intensive + adj_p=zeros(1,m)*NaN; + [wtd_p_sorted, wtd_p_sindex] = sort( wtd_p ); + nextfill = 1; + for k = 1 : m + if wtd_p_sindex(k)>=nextfill + adj_p(nextfill:wtd_p_sindex(k)) = wtd_p_sorted(k); + nextfill = wtd_p_sindex(k)+1; + if nextfill>m + break; + end; + end; + end; + adj_p=reshape(adj_p(unsort_ids),s); +end + +rej=p_sorted<=thresh; +max_id=find(rej,1,'last'); %find greatest significant pvalue +if isempty(max_id), + crit_p=0; + h=pvals*0; + adj_ci_cvrg=NaN; +else + crit_p=p_sorted(max_id); + h=pvals<=crit_p; + adj_ci_cvrg=1-thresh(max_id); +end + +if strcmpi(report,'yes'), + n_sig=sum(p_sorted<=crit_p); + if n_sig==1, + fprintf('Out of %d tests, %d is significant using a false discovery rate of %f.\n',m,n_sig,q); + else + fprintf('Out of %d tests, %d are significant using a false discovery rate of %f.\n',m,n_sig,q); + end + if strcmpi(method,'pdep'), + fprintf('FDR/FCR procedure used is guaranteed valid for independent or positively dependent tests.\n'); + else + fprintf('FDR/FCR procedure used is guaranteed valid for independent or dependent tests.\n'); + end +end + + + + diff --git a/Visualize_Results/scripts/helper_fun/helper_boxplot.m b/Visualize_Results/scripts/helper_fun/helper_boxplot.m new file mode 100644 index 0000000..80c1a04 --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/helper_boxplot.m @@ -0,0 +1,91 @@ +function helper_boxplot(seed1, seed2, seeds, ROIs) + +%% Run with MATLAB 2021b and later. + +n_seeds = length(seeds); + +% labels +label = ROIs.name; +for i = 1:length(label) + myLabel{i} = strrep(label{i},'_','-'); + myLabel{i} = myLabel{i}(1:end-6); +end + +%% dorsal stream +idx_dorsal = 1:17; +dorsal = myLabel(idx_dorsal); +n_targets = length(dorsal); +ROI_FEF = seed1(ROIs.idx(idx_dorsal),:)'; +ROI_IFJ = seed2(ROIs.idx(idx_dorsal),:)'; + +tmp_Seeds = repmat(repmat(seeds,1,1), 1, 55*n_targets)'; +Seeds = tmp_Seeds(:); +Targets = repmat(repmat(dorsal',55,1), n_seeds, 1); +Subjects = repmat(repmat(1:55,1,n_targets)', n_seeds, 1); +targetsOrder = dorsal; +Conn_values = [ROI_FEF(:);ROI_IFJ(:)]; +results_tab = table(Seeds,Targets,Conn_values); +results_tab.Targets = categorical(results_tab.Targets, targetsOrder); +figure; h_one = boxchart(results_tab.Targets, results_tab.Conn_values, 'GroupByColor', results_tab.Seeds); +legend +title('Dorsam Visual Stream', 'fontweight','bold', 'FontSize', 15) +set(gcf, 'Position', get(0, 'Screensize')); + +clear results_tab + +h_one(1).SeriesIndex = 2; +h_one(2).SeriesIndex = 1; +h_one(1).WhiskerLineStyle = ':'; +h_one(2).WhiskerLineStyle = ':'; +h_one(1).MarkerSize = 6; +h_one(2).MarkerSize = 6; + + +xlabel('ROIs', 'fontweight','bold', 'FontSize', 12) +ylabel('Connectivity Values', 'fontweight','bold', 'FontSize', 12) +ax=gca; +ax.FontSize = 12; + +print(gcf,'C:\Users\ASUS\Desktop\dorsal_boxchart.png','-dpng','-r300'); + + +%% ventral stream +idx_ventral = 18:33; +ventral = myLabel(idx_ventral); +n_targets = length(ventral); +ROI_FEF = seed1(ROIs.idx(idx_ventral),:)'; +ROI_IFJ = seed2(ROIs.idx(idx_ventral),:)'; + +tmp_Seeds = repmat(repmat(seeds,1,1), 1, 55*n_targets)'; +Seeds = tmp_Seeds(:); +Targets = repmat(repmat(ventral',55,1), n_seeds, 1); +Subjects = repmat(repmat(1:55,1,n_targets)', n_seeds, 1); +targetsOrder = ventral; +Conn_values = [ROI_FEF(:);ROI_IFJ(:)]; +results_tab = table(Seeds,Targets,Conn_values); +results_tab.Targets = categorical(results_tab.Targets, targetsOrder); +figure; h_two = boxchart(results_tab.Targets, results_tab.Conn_values, 'GroupByColor', results_tab.Seeds); +legend +title('Ventral Visual Stream', 'fontweight','bold', 'FontSize', 15) +set(gcf, 'Position', get(0, 'Screensize')); + +clear results_tab + +h_two(1).SeriesIndex = 2; +h_two(2).SeriesIndex = 1; +h_two(1).WhiskerLineStyle = ':'; +h_two(2).WhiskerLineStyle = ':'; +h_two(1).MarkerSize = 6; +h_two(2).MarkerSize = 6; + + +xlabel('ROIs', 'fontweight','bold', 'FontSize', 12) +ylabel('Connectivity Values', 'fontweight','bold', 'FontSize', 12) +ax=gca; +ax.FontSize = 12; + + +print(gcf,'C:\Users\ASUS\Desktop\ventral_boxchart.png','-dpng','-r300'); + + +end \ No newline at end of file diff --git a/Visualize_Results/scripts/helper_fun/helper_effectiveConn_fsaverage.m b/Visualize_Results/scripts/helper_fun/helper_effectiveConn_fsaverage.m new file mode 100644 index 0000000..e69dff6 --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/helper_effectiveConn_fsaverage.m @@ -0,0 +1,508 @@ +function results = helper_effectiveConn_fsaverage(cimec) + +inputfile_conn = cimec.inputfile_conn; +band_name = cimec.band_name; +seed_target = cimec.seed_target; +conn_metric = cimec.conn_metric; +statistics = cimec.statistics; +colorlim_value = cimec.colorlim_value; +flag = cimec.flag; + +%% +addpath(genpath(fullfile(fileparts(fileparts(pwd)), 'helper_fun'))); + +statistics.contrast = true; +statistics.stat = true; +statistics.f_stats = @signrank; +statistics.f_correction = @fdr_bh; + +IFJ_merged = false; +ROIs.do = true; % fixed + +colorlim = true; + +%% +% fsaverage +inflated = true; +if strcmp(statistics.analysis_type, 'exploratory') + smoothing = '70'; % percent +else + smoothing = '100'; % percent +end +inputfile_fs = fullfile(load_path('inputfile_fs'), 'fsaverage.mat'); +inputfile_fs_inflated = fullfile(load_path('inputfile_fs'), sprintf('fsaverage_inflated_%s.mat', smoothing)); + +% seed to target +dorsal_R = { 'R_V6_ROI R', 'R_V6A_ROI R', 'R_V7_ROI R', 'R_IPS1_ROI R', 'R_IP1_ROI R', 'R_MIP_ROI R', ... + 'R_VIP_ROI R', 'R_LIPd_ROI R', 'R_LIPv_ROI R', 'R_7AL_ROI R', 'R_7PC_ROI R', 'R_7PL_ROI R', ... + 'R_V3A_ROI R', 'R_V3B_ROI R', 'R_V3CD_ROI R', 'R_LO3_ROI R', 'R_MT_ROI R'}; +ventral_R = {'R_V8_ROI R', 'R_VMV1_ROI R', 'R_VMV2_ROI R', 'R_VMV3_ROI R', 'R_PHT_ROI R', 'R_PH_ROI R', ... + 'R_TE1a_ROI R', 'R_TE1m_ROI R', 'R_TE1p_ROI R', 'R_TE2a_ROI R', 'R_TE2p_ROI R', ... + 'R_TF_ROI R', 'R_TGv_ROI R', 'R_FFC_ROI R', 'R_PIT_ROI R', 'R_VVC_ROI R'}; + +dorsal_L = { 'L_V6_ROI L', 'L_V6A_ROI L', 'L_V7_ROI L', 'L_IPS1_ROI L', 'L_IP1_ROI L', 'L_MIP_ROI L', ... + 'L_VIP_ROI L', 'L_LIPd_ROI L', 'L_LIPv_ROI L', 'L_7AL_ROI L', 'L_7PC_ROI L', 'L_7PL_ROI L', ... + 'L_V3A_ROI L', 'L_V3B_ROI L', 'L_V3CD_ROI L', 'L_LO3_ROI L', 'L_MT_ROI L'}; +ventral_L = {'L_V8_ROI L', 'L_VMV1_ROI L', 'L_VMV2_ROI L', 'L_VMV3_ROI L', 'L_PHT_ROI L', 'L_PH_ROI L', ... + 'L_TE1a_ROI L', 'L_TE1m_ROI L', 'L_TE1p_ROI L', 'L_TE2a_ROI L', 'L_TE2p_ROI L', ... + 'L_TF_ROI L', 'L_TGv_ROI L', 'L_FFC_ROI L', 'L_PIT_ROI L', 'L_VVC_ROI L'}; + + +if strcmp(seed_target, 'right-right') + between = {'R_IFJa_ROI R', 'R_IFJp_ROI R', 'R_FEF_ROI R'}; % 1-3 & 2-3 + ROIs.name = ['R_FEF_ROI R', 'R_IFJa_ROI R', 'R_IFJp_ROI R', dorsal_R, ventral_R]; +elseif strcmp(seed_target, 'right-left') + between = {'R_IFJa_ROI R', 'R_IFJp_ROI R', 'R_FEF_ROI R'}; % 1-3 & 2-3 + ROIs.name = ['R_FEF_ROI R', 'R_IFJa_ROI R', 'R_IFJp_ROI R', dorsal_L ventral_L]; +elseif strcmp(seed_target, 'left-left') + between = {'L_IFJa_ROI L', 'L_IFJp_ROI L', 'L_FEF_ROI L'}; % 1-3 & 2-3 + ROIs.name = ['L_FEF_ROI L', 'L_IFJa_ROI L', 'L_IFJp_ROI L', dorsal_L ventral_L]; +elseif strcmp(seed_target, 'left-right') + between = {'L_IFJa_ROI L', 'L_IFJp_ROI L', 'L_FEF_ROI L'}; % 1-3 & 2-3 + ROIs.name = ['L_FEF_ROI L', 'L_IFJa_ROI L', 'L_IFJp_ROI L', dorsal_R ventral_R]; +end + + +%% Inputs +cimec = struct; +cimec.between = between; +cimec.statistics = statistics; +cimec.band_name = band_name; +cimec.IFJ_merged = IFJ_merged; +cimec.ROIs = ROIs; +cimec.colorlim_value = colorlim_value; +cimec.colorlim = colorlim; +cimec.inputfile_conn = inputfile_conn; +cimec.inputfile_fs = inputfile_fs; +cimec.inflated = inflated; +cimec.inputfile_fs_inflated = inputfile_fs_inflated; +cimec.flag = flag; +cimec.conn_metric = conn_metric; +cimec.seed_target = seed_target; + +results = visualize_connectivity(cimec); + +%% Function + function results = visualize_connectivity(cimec) + + %% Info about the data + inputfile_conn = cimec.inputfile_conn; % e.g. '...\HCP_Subjects\Outputs\connectivity_results\175237_connectivityResults__2s'; + inputfile_fs = cimec.inputfile_fs; % e.g. '...\HCP_Subjects\Outputs\fsaverage_anatomy\175237_fsaverage'; + + contrast = cimec.statistics.contrast; % contrast map or not + between = cimec.between; % {'R_IFJa_ROI R', 'R_IFJp_ROI R', 'R_FEF_ROI R'}; % 1-2 & 1-3 + + stat = cimec.statistics.stat; + f_stats = cimec.statistics.f_stats; % e.g. @signrank + alpha = cimec.statistics.alpha; % e.g. 0.05 + f_correction = cimec.statistics.f_correction; % correction method e.g. @fdr_bh + corrected = cimec.statistics.corrected; % e.g. true + + IFJ_merged = cimec.IFJ_merged; + + ROIs = cimec.ROIs; + + band_name = cimec.band_name; % e.g. 'alpha'; + colorlim_value = cimec.colorlim_value; % e.g. [0 0.10]; + colorlim = cimec.colorlim; % e.g. true; + + inflated = cimec.inflated; % e.g. true + inputfile_fs_inflated = cimec.inputfile_fs_inflated; + + flag = cimec.flag; + conn_metric = cimec.conn_metric; + seed_target = cimec.seed_target; + + %% Frequency bands + band_names = { 'delta' , 'theta' , 'alpha' , 'beta' , 'gamma' }; + freq_bands = [ 1 4 ; 4 8 ; 8 13 ; 13 30 ; 30 100 ]; + + %% Load data + % avgFreq_icoh + load(inputfile_conn, 'group_results'); + conn_results = group_results; + + % fsaverage + load(inputfile_fs, 'fsaverage') + + %% Prepare the anatomy + % Parcellation info (Brainstorm -> Fieldtrip) + tess_cortex_pial = fsaverage; + + if inflated == true + fs_inflated = load(inputfile_fs_inflated); + tess_cortex_pial.Vertices = fs_inflated.Vertices; + tess_cortex_pial.Faces = fs_inflated.Faces; + end + + Scouts = tess_cortex_pial.Atlas(end-1).Scouts; + Scouts_Vertices = {Scouts.Vertices}'; + Scouts_Label = {Scouts.Label}'; + + parcellation = zeros(length([tess_cortex_pial.Vertices]), 1); + for idx = 1:length(Scouts_Vertices) + per_parcel = Scouts_Vertices{idx}; + + for ii = 1:length(per_parcel) + parcellation(per_parcel(ii)) = idx; + end + end + + brainordinate.parcellation = parcellation; + brainordinate.parcellationlabel = Scouts_Label; + brainordinate.pos = tess_cortex_pial.Vertices; + brainordinate.tri = tess_cortex_pial.Faces; + brainordinate.brainstructure = tess_cortex_pial.SulciMap+1; + brainordinate.brainstructurelabel = {'CORTEX_LEFT', 'CORTEX_RIGHT'}'; + + %% Prepare the connectivity matrix + conn_mat = struct; + conn_mat.label = Scouts_Label; + conn_mat.band_names = band_name; + conn_mat.brainordinate = brainordinate; + + idx_band = find(strcmp(band_names, band_name)); + tmp_data = group_results(idx_band).data_perSubject; + tmp_data(isnan(tmp_data)) = 1; + + %% seeds + idx_IFJa = find(strcmp(conn_mat.label, between{1})); + idx_IFJp = find(strcmp(conn_mat.label, between{2})); + idx_FEF = find(strcmp(conn_mat.label, between{3})); + + IFJa_12(:, :) = tmp_data(idx_IFJa,:,:); + IFJp_12(:, :) = tmp_data(idx_IFJp,:,:); + FEF_12(:, :) = tmp_data(idx_FEF,:,:); + IFJa_21(:, :) = tmp_data(:,idx_IFJa,:); + IFJp_21(:, :) = tmp_data(:,idx_IFJp,:); + FEF_21(:, :) = tmp_data(:,idx_FEF,:); + + + % subtract two directions + % positive values mean seed region predicts the target + % negative values mean target is predicted by the seed region + IFJa(:, :) = IFJa_12 - IFJa_21; + IFJp(:, :) = IFJp_12 - IFJp_21; + FEF(:, :) = FEF_12 - FEF_21; + + %% ROIs/exploratory + if strcmp(statistics.analysis_type, 'exploratory') + all_labels = conn_mat.label; + if strcmp(between{1}(1), 'R') + ROIs.name = all_labels(181:360); + elseif strcmp(between{1}(1), 'L') + ROIs.name = all_labels(1:180); + end + end + + nROIs = length(ROIs.name); + for ii = 1:nROIs + ROIs.idx(ii) = find(strcmp(conn_mat.label, ROIs.name{ii})); + end + + + mask = zeros(360,1); + mask(ROIs.idx) = 1; + + %% normality test + % stat_IFJa = IFJa; + % stat_IFJa(stat_IFJa(:,1)==1) = NaN; + % stat_FEF = FEF; + % stat_FEF(stat_FEF(:,1)==1) = NaN; + % + % + % for ii = 1:nROIs + % figure + % % histogram(atanh(stat_IFJa(ROIs.idx(ii),:)-stat_FEF(ROIs.idx(ii),:))); + % % histogram(stat_IFJa(ROIs.idx(ii),:)-stat_FEF(ROIs.idx(ii),:)); + % % histogram(stat_IFJa(ROIs.idx(ii),:)); + % histogram(atanh(stat_IFJa(ROIs.idx(ii),:))); + % end + % + % figure + % hist(atanh(stat_IFJa(ROIs.idx(33),:)-stat_FEF(ROIs.idx(33),:))); + % figure + % hist(stat_IFJa(ROIs.idx(33),:)-stat_FEF(ROIs.idx(33),:)); + + + %% + if contrast == true + %% Inferential statistics + %% Wilcoxon signed rank test + if IFJ_merged == false + % R_IFJa + for ii = 1:nROIs + pair_1 = IFJa(ROIs.idx(ii),:)'; + pair_2 = 0; + + [P, H, ~] = f_stats(pair_1, pair_2, 'Alpha', alpha, 'Tail', 'both'); + + H_Wilcoxon_IFJa(ii) = H; + P_Wilcoxon_IFJa(ii) = P; + end + + % R_IFJp + for ii = 1:nROIs + pair_1 = IFJp(ROIs.idx(ii),:)'; + pair_2 = 0; + + [P, H, ~] = f_stats(pair_1, pair_2, 'Alpha', alpha, 'Tail', 'both'); + + H_Wilcoxon_IFJp(ii) = H; + P_Wilcoxon_IFJp(ii) = P; + end + + % FEF + for ii = 1:nROIs + pair_1 = FEF(ROIs.idx(ii),:)'; + pair_2 = 0; + + [P, H, ~] = f_stats(pair_1, pair_2, 'Alpha', alpha, 'Tail', 'both'); + + H_Wilcoxon_FEF(ii) = H; + P_Wilcoxon_FEF(ii) = P; + end + end + %% FDR correction + if corrected == false + % clear NaNs + H_Wilcoxon_IFJa(isnan(H_Wilcoxon_IFJa)) = 0; + H_Wilcoxon_IFJp(isnan(H_Wilcoxon_IFJp)) = 0; + + % mask connectivity values + group_IFJa_contrast = mean(IFJa, 2); + right_group_IFJa = group_IFJa_contrast(ROIs.idx); + right_group_IFJa(~H_Wilcoxon_IFJa) = 0; + group_IFJa = zeros(360,1); + group_IFJa(ROIs.idx) = right_group_IFJa; + stats_IFJa = group_IFJa; + + group_IFJp_contrast = mean(IFJp, 2); + right_group_IFJp = group_IFJp_contrast(ROIs.idx); + right_group_IFJp(~H_Wilcoxon_IFJp) = 0; + group_IFJp = zeros(360,1); + group_IFJp(ROIs.idx) = right_group_IFJp; + stats_IFJp = group_IFJp; + + % clear NaNs + H_Wilcoxon_FEF(isnan(H_Wilcoxon_FEF)) = 0; + + % mask connectivity values + group_FEF_contrast = mean(FEF, 2); + right_group_FEF = group_FEF_contrast(ROIs.idx); + right_group_FEF(~H_Wilcoxon_FEF) = 0; + group_FEF = zeros(360,1); + group_FEF(ROIs.idx) = right_group_FEF; + stats_FEF = group_FEF; + + elseif corrected == true + % FDR correction + P_Wilcoxon_IFJa(isnan(P_Wilcoxon_IFJa)) = 1; + P_Wilcoxon_IFJp(isnan(P_Wilcoxon_IFJp)) = 1; + P_Wilcoxon_FEF(isnan(P_Wilcoxon_FEF)) = 1; + [h_adj_IFJa, ~, ~, p_adj_IFJa] = f_correction(P_Wilcoxon_IFJa(~isnan(P_Wilcoxon_IFJa)), alpha, 'pdep', 'yes'); + [h_adj_IFJp, ~, ~, p_adj_IFJp] = f_correction(P_Wilcoxon_IFJp(~isnan(P_Wilcoxon_IFJp)), alpha, 'pdep', 'yes'); + [h_adj_FEF, ~, ~, p_adj_FEF] = f_correction(P_Wilcoxon_FEF(~isnan(P_Wilcoxon_FEF)), alpha, 'pdep', 'yes'); + + + % z-score + % for two-tailed: abs(norminv(p/2)); + z_adj_IFJa = abs(norminv(p_adj_IFJa/2)); + z_adj_IFJp = abs(norminv(p_adj_IFJp/2)); + z_adj_FEF = abs(norminv(p_adj_FEF/2)); + + + group_IFJa_contrast = mean(IFJa, 2); + right_group_IFJa = group_IFJa_contrast(ROIs.idx); + right_group_IFJa(~h_adj_IFJa) = 0; + group_IFJa = zeros(360,1); + group_IFJa(ROIs.idx) = z_adj_IFJa'.*sign(right_group_IFJa); + stats_IFJa = group_IFJa; + + group_IFJp_contrast = mean(IFJp, 2); + right_group_IFJp = group_IFJp_contrast(ROIs.idx); + right_group_IFJp(~h_adj_IFJp) = 0; + group_IFJp = zeros(360,1); + group_IFJp(ROIs.idx) = z_adj_IFJp'.*sign(right_group_IFJp); + stats_IFJp = group_IFJp; + + group_FEF_contrast = mean(FEF, 2); + right_group_FEF = group_FEF_contrast(ROIs.idx); + right_group_FEF(~h_adj_FEF) = 0; + group_FEF = zeros(360,1); + group_FEF(ROIs.idx) = z_adj_FEF'.*sign(right_group_FEF); + stats_FEF = group_FEF; + + % % mask connectivity values + % group_IFJa_contrast = mean(IFJa, 2); + % right_group_IFJa = group_IFJa_contrast(ROIs.idx); + % right_group_IFJa(~h_adj_IFJa) = 0; + % group_IFJa = zeros(360,1); + % group_IFJa(ROIs.idx) = right_group_IFJa; + % stats_IFJa = group_IFJa; + % + % group_IFJp_contrast = mean(IFJp, 2); + % right_group_IFJp = group_IFJp_contrast(ROIs.idx); + % right_group_IFJp(~h_adj_IFJp) = 0; + % group_IFJp = zeros(360,1); + % group_IFJp(ROIs.idx) = right_group_IFJp; + % stats_IFJp = group_IFJp; + +% % mask connectivity values +% group_FEF_contrast = mean(FEF, 2); +% right_group_FEF = group_FEF_contrast(ROIs.idx); +% right_group_FEF(~h_adj_FEF) = 0; +% group_FEF = zeros(360,1); +% group_FEF(ROIs.idx) = right_group_FEF; +% stats_FEF = group_FEF; + end + + %% Final matrix (360x360) + if IFJ_merged == false + tmp_data = diag(ones(360,1)); + + tmp_data(:,idx_IFJa) = stats_IFJa'; + tmp_data(:,idx_IFJp) = stats_IFJp'; + tmp_data(:,idx_FEF) = stats_FEF'; + end + + end + + %% To make black the areas outside ROI + cbar = customMap_effectiveConn; + if strcmp(cimec.statistics.analysis_type, 'ROI') + n = size(cbar, 1); + up = cbar(1:n/2,:); + middle = [0.5 0.5 0.5]; + down = cbar(n/2+1:n,:); + set(0, 'DefaultFigureColormap', [up;middle;down]); + else + set(0, 'DefaultFigureColormap', cbar); + end + + if strcmp(statistics.analysis_type, 'ROI') + tmp_data(ROIs.idx(~h_adj_IFJa),idx_IFJa) = 0.04; + tmp_data(ROIs.idx(~h_adj_IFJp),idx_IFJp) = 0.04; + tmp_data(ROIs.idx(~h_adj_FEF),idx_FEF) = 0.04; + tmp_data(idx_IFJa,ROIs.idx(~h_adj_IFJa)) = 0.04; + tmp_data(idx_IFJp,ROIs.idx(~h_adj_IFJp)) = 0.04; + tmp_data(idx_FEF,ROIs.idx(~h_adj_FEF)) = 0.04; + end + + + %% + conn_mat.cohspctrm = tmp_data; + + % output + results.helper.idx_IFJa = idx_IFJa; + results.helper.idx_IFJp = idx_IFJp; + results.helper.idx_FEF = idx_FEF; + results.helper.ROIs = ROIs; + results.helper.brainordinate = brainordinate; + results.conn_matrix = tmp_data; + results.label = conn_mat.label; + results.stat.h_adj_IFJa = h_adj_IFJa; + results.stat.h_adj_IFJp = h_adj_IFJp; + results.stat.h_adj_FEF = h_adj_FEF; + + results.perSubject.IFJa_12 = IFJa_12; + results.perSubject.IFJp_12 = IFJp_12; + results.perSubject.FEF_12 = FEF_12; + + results.perSubject.IFJa_21 = IFJa_21; + results.perSubject.IFJp_21 = IFJp_21; + results.perSubject.FEF_21 = FEF_21; + + %% Visualize using the helper function + if flag.fsaverage == true + % in order to run the function + conn_mat.brainordinate.pos = conn_mat.brainordinate.pos*1000; + + %% save the figure + if and(colorlim == true, flag.save_fig == true) + + if seed_target(1) == 'r' + seed = 'IFJa'; + cimec.seed = seed; + pos2d_IFJa_R = [52.5323 -69.3638 71.0650; 52.5323 69.4183 71.0650]; + if strcmp(statistics.analysis_type, 'ROI') + pos2d_IFJa_R = [52.5323 -69.3638 71.0650; 52.5323 69.4183 71.0650]; + else + pos2d_IFJa_R = [49.7454 -69.6590 70.1540; 49.7454 69.1230 70.1540]; + end + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_IFJa_R); + helper_save_figure(cimec); + close all; + + seed = 'IFJp'; + cimec.seed = seed; + if strcmp(statistics.analysis_type, 'ROI') + pos2d_IFJp_R = [45.6139 -69.3638 76.8706 ; 45.6139 69.4183 76.8706]; + else + pos2d_IFJp_R = [44.7673 -69.6590 76.5544 ; 44.7673 69.1230 76.5544]; + end + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_IFJp_R); + helper_save_figure(cimec); + close all; + + seed = 'FEF'; + cimec.seed = seed; + if strcmp(statistics.analysis_type, 'ROI') + pos2d_FEF_R = [30.1309 -69.3638 103.7783 ; 30.1309 69.4183 103.7783]; + else + pos2d_FEF_R = [30.8997 -69.6590 101.8004 ; 30.8997 69.1230 101.8004]; + end + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_FEF_R); + helper_save_figure(cimec); + close all; + + elseif seed_target(1) == 'l' + seed = 'IFJa'; + cimec.seed = seed; + if strcmp(statistics.analysis_type, 'ROI') + pos2d_IFJa_L = [50.5105 69.4183 69.5634 ; 53.1465 -69.3638 68.2860]; + else + pos2d_IFJa_L = [47.6119 69.1230 68.3761 ; 47.6119 -69.6590 68.3761]; + end + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_IFJa_L); + helper_save_figure(cimec); + close all; + + seed = 'IFJp'; + cimec.seed = seed; + cimec.seed = seed; + if strcmp(statistics.analysis_type, 'ROI') + pos2d_IFJp_L = [43.6193 69.4183 80.0619 ; 43.6193 -69.3638 80.0619]; + else + pos2d_IFJp_L = [40.8559 69.1230 76.9099 ; 40.8559 -69.6590 76.9099]; + end + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_IFJp_L); + helper_save_figure(cimec); + close all; + + seed = 'FEF'; + cimec.seed = seed; + if strcmp(statistics.analysis_type, 'ROI') + pos2d_FEF_L = [30.6545 69.4183 105.7922 ; 30.6545 -69.3638 105.7922]; + else + pos2d_FEF_L = [33.7443 69.1230 101.0893 ; 33.7443 -69.6590 101.0893]; + end + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_FEF_L); + helper_save_figure(cimec); + close all; + end + + %% show the figure + elseif colorlim == true + tutorial_nwa_connectivityviewer(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours); + elseif colorlim == false + tutorial_nwa_connectivityviewer(conn_mat, 'cohspctrm', seed_target, flag.contours); + end + + if contrast == false + colormap('parula') + end + end + + end +end diff --git a/Visualize_Results/scripts/helper_fun/helper_functionConn_circularGraph.m b/Visualize_Results/scripts/helper_fun/helper_functionConn_circularGraph.m new file mode 100644 index 0000000..23eab99 --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/helper_functionConn_circularGraph.m @@ -0,0 +1,360 @@ +function helper_functionConn_circularGraph(results, fontSize, flag) + +if flag.contrast == true + figure; + contrast_pair = 'FEF-IFJp'; + helper_circularGraph_contrast(results, fontSize, flag, contrast_pair) + contrast_pair = 'FEF-IFJa'; + helper_circularGraph_contrast(results, fontSize, flag, contrast_pair) +else + figure; + seed = 'FEF'; + helper_circularGraph_notContrast(results, fontSize, flag, seed) + figure; + seed = 'IFJa'; + helper_circularGraph_notContrast(results, fontSize, flag, seed) + figure; + seed = 'IFJp'; + helper_circularGraph_notContrast(results, fontSize, flag, seed) +end + +%% + function helper_circularGraph_contrast(results, fontSize, flag, pair) + + tmp_mat = results.conn_matrix; + idx_FEF = results.helper.idx_FEF; + idx_IFJa = results.helper.idx_IFJa; + idx_IFJp = results.helper.idx_IFJp; + ROIs = results.helper.ROIs; + analysis_type = results.stat.analysis_type; + + band_name = results.cimec.band_name ; + conn_metric = results.cimec.conn_metric; + seed_target = results.cimec.seed_target; + alpha = char(string(results.cimec.statistics.alpha)); + correction = results.cimec.statistics.correction; + duration = results.cimec.flag.duration; + + % https://sashamaps.net/docs/resources/20-colors/ + % https://html-color.codes/blue + FEF_color_1 = [255,0,0]/255; + FEF_color_2 = [255,160,122]/255; + IFJa_color = [0,0,255]/255; %[0, 92, 171]/255; + IFJp_color = [173,216,230]/255; %[70, 240, 240]/255; + + if strcmp(pair, 'FEF-IFJa') + + %% ROI analysis + if strcmp(analysis_type, 'ROI') + + idxs = [idx_FEF, idx_IFJa, idx_IFJp, flip(ROIs.idx)]; + data = tmp_mat(idxs, idxs); + data(isnan(data)) = 0; + + % to make FEF results explicit + data = data*-1; + + data(1,data(2,:)<0) = abs(data(2,data(2,:)<0)); + data(2,data(2,:)<0) = 0; + data(:,1) = data(1,:); + data(:,2) = data(2,:); + data(:,3) = 0; + data(3,:) = 0; + + label = results.label(idxs); + + if flag.collapseTargets == true + % labels + myLabel = cell(length(label)); + for i = 1:length(label) + myLabel{i} = strrep(label{i},'_','-'); + myLabel{i} = myLabel{i}(3:end-6); + end + else + % labels + myLabel = cell(length(label)); + for i = 1:length(label) + myLabel{i} = strrep(label{i},'_','-'); + myLabel{i} = myLabel{i}(1:end-6); + end + end + + % color + % https://sashamaps.net/docs/resources/20-colors/ + myColorMap = repmat([1 1 1], length(idxs), 1); + myColorMap(1,:) = FEF_color_1; + myColorMap(2,:) = IFJa_color; + + %% Exploratory analysis + elseif strcmp(analysis_type, 'exploratory') + + figure; + + idxs = ROIs.idx; + data = tmp_mat(idxs, idxs); + data(isnan(data)) = 0; + + % to make FEF results explicit + data = data*-1; + + data(1,data(2,:)<0) = abs(data(2,data(2,:)<0)); + data(2,data(2,:)<0) = 0; + data(:,1) = data(1,:); + data(:,2) = data(2,:); + data(:,3) = 0; + data(3,:) = 0; + + label = results.label(idxs); + + if flag.collapseTargets == true + % labels + myLabel = cell(length(label)); + for i = 1:length(label) + myLabel{i} = strrep(label{i},'_','-'); + myLabel{i} = myLabel{i}(3:end-6); + end + else + % labels + myLabel = cell(length(label)); + for i = 1:length(label) + myLabel{i} = strrep(label{i},'_','-'); + myLabel{i} = myLabel{i}(1:end-6); + end + end + + % color + % https://sashamaps.net/docs/resources/20-colors/ + myColorMap = repmat([1 1 1], length(idxs), 1); + myColorMap(idx_FEF-180,:) = FEF_color_1; + myColorMap(idx_IFJa-180,:) = IFJa_color; + end + + %% Circular graph for all together + circularGraph(data', 'Colormap', myColorMap, 'Label', myLabel); + title('All Frequency Bands') + set(gcf, 'Position', get(0, 'Screensize')); + + fg1 = gca; + fg2 = get(fg1, 'children'); + + for ii = 1:length(fg2) + try + tmp_fg = fg2(ii); + tmp_fg.FontSize = fontSize; + catch + fprintf('Almost ready...\n', i); + end + end + + if flag.save_fig == true + name = [ contrast_pair '_circularGraph_' conn_metric '_' band_name '_' seed_target '_' analysis_type '_' alpha '_' correction '_' duration]; + print(gcf, [flag.dir_fig name '.png'], '-dpng', '-r300'); + close all + end + + elseif strcmp(pair, 'FEF-IFJp') + %% ROI analysis + if strcmp(analysis_type, 'ROI') + + idxs = [idx_FEF, idx_IFJa, idx_IFJp, flip(ROIs.idx)]; + data = tmp_mat(idxs, idxs); + data(isnan(data)) = 0; + + % to make FEF results explicit + data = data*-1; + + data(1,data(3,:)<0) = abs(data(3,data(3,:)<0)); + data(3,data(3,:)<0) = 0; + data(:,1) = data(1,:); + data(:,3) = data(3,:); + data(:,2) = 0; + data(2,:) = 0; + + label = results.label(idxs); + + if flag.collapseTargets == true + % labels + myLabel = cell(length(label)); + for i = 1:length(label) + myLabel{i} = strrep(label{i},'_','-'); + myLabel{i} = myLabel{i}(3:end-6); + end + else + % labels + myLabel = cell(length(label)); + for i = 1:length(label) + myLabel{i} = strrep(label{i},'_','-'); + myLabel{i} = myLabel{i}(1:end-6); + end + end + + % color + % https://sashamaps.net/docs/resources/20-colors/ + myColorMap = repmat([1 1 1], length(idxs), 1); + myColorMap(1,:) = FEF_color_2; + myColorMap(3,:) = IFJp_color; + + %% Exploratory analysis + elseif strcmp(analysis_type, 'exploratory') + + idxs = ROIs.idx; + data = tmp_mat(idxs, idxs); + data(isnan(data)) = 0; + + % to make FEF results explicit + data = data*-1; + + data(1,data(3,:)<0) = abs(data(3,data(3,:)<0)); + data(3,data(3,:)<0) = 0; + data(:,1) = data(1,:); + data(:,3) = data(3,:); + data(:,2) = 0; + data(2,:) = 0; + + label = results.label(idxs); + + if flag.collapseTargets == true + % labels + myLabel = cell(length(label)); + for i = 1:length(label) + myLabel{i} = strrep(label{i},'_','-'); + myLabel{i} = myLabel{i}(3:end-6); + end + else + % labels + myLabel = cell(length(label)); + for i = 1:length(label) + myLabel{i} = strrep(label{i},'_','-'); + myLabel{i} = myLabel{i}(1:end-6); + end + end + + % color + % https://sashamaps.net/docs/resources/20-colors/ + myColorMap = repmat([1 1 1], length(idxs), 1); + myColorMap(idx_FEF-180,:) = FEF_color_2; + myColorMap(idx_IFJp-180,:) = IFJp_color; + end + %% Circular graph for all together + circularGraph(data', 'Colormap', myColorMap, 'Label', myLabel); + title('All Frequency Bands') + set(gcf, 'Position', get(0, 'Screensize')); + + fg1 = gca; + fg2 = get(fg1, 'children'); + + for ii = 1:length(fg2) + try + tmp_fg = fg2(ii); + tmp_fg.FontSize = fontSize; + catch + fprintf('Almost ready...\n', i); + end + end + + if flag.save_fig == true + name = [contrast_pair '_circularGraph_' conn_metric '_' band_name '_' seed_target '_' analysis_type '_' alpha '_' correction '_' duration]; + print(gcf, [flag.dir_fig name '.png'], '-dpng', '-r300'); + close all + end + + + end + end + + function helper_circularGraph_notContrast(results, fontSize, flag, seed) + + tmp_mat = results.conn_matrix; + idx_FEF = results.helper.idx_FEF; + idx_IFJa = results.helper.idx_IFJa; + idx_IFJp = results.helper.idx_IFJp; + ROIs = results.helper.ROIs; + analysis_type = results.stat.analysis_type; + + band_name = results.cimec.band_name ; + conn_metric = results.cimec.conn_metric; + seed_target = results.cimec.seed_target; + alpha = char(string(results.cimec.statistics.alpha)); + correction = results.cimec.statistics.correction; + duration = results.cimec.flag.duration; + + % https://sashamaps.net/docs/resources/20-colors/ + % https://html-color.codes/blue + FEF_color = [255,0,0]/255; + IFJa_color = [0,0,255]/255; %[0, 92, 171]/255; + IFJp_color = [173,216,230]/255; %[70, 240, 240]/255; + + %% ROI analysis + if strcmp(analysis_type, 'ROI') + + if strcmp(seed, 'FEF') + idxs = [idx_FEF, flip(ROIs.idx)]; + data = tmp_mat(idxs, idxs); + data(isnan(data)) = 0; + data(:,1) = data(1,:); + elseif strcmp(seed, 'IFJa') + idxs = [idx_IFJa, flip(ROIs.idx)]; + data = tmp_mat(idxs, idxs); + data(isnan(data)) = 0; + data(:,1) = data(1,:); + elseif strcmp(seed, 'IFJp') + idxs = [idx_IFJp, flip(ROIs.idx)]; + data = tmp_mat(idxs, idxs); + data(isnan(data)) = 0; + data(:,1) = data(1,:); + end + + label = results.label(idxs); + + if flag.collapseTargets == true + % labels + myLabel = cell(length(label)); + for i = 1:length(label) + myLabel{i} = strrep(label{i},'_','-'); + myLabel{i} = myLabel{i}(3:end-6); + end + else + % labels + myLabel = cell(length(label)); + for i = 1:length(label) + myLabel{i} = strrep(label{i},'_','-'); + myLabel{i} = myLabel{i}(1:end-6); + end + end + + % color + % https://sashamaps.net/docs/resources/20-colors/ + myColorMap = repmat([1 1 1], length(idxs), 1); + if strcmp(seed, 'FEF') + myColorMap(1,:) = FEF_color; + elseif strcmp(seed, 'IFJa') + myColorMap(1,:) = IFJa_color; + elseif strcmp(seed, 'IFJp') + myColorMap(1,:) = IFJp_color; + end + end + + %% Circular graph for all together + circularGraph(data', 'Colormap', myColorMap, 'Label', myLabel); + title('All Frequency Bands') + set(gcf, 'Position', get(0, 'Screensize')); + + fg1 = gca; + fg2 = get(fg1, 'children'); + + for ii = 1:length(fg2) + try + tmp_fg = fg2(ii); + tmp_fg.FontSize = fontSize; + catch + fprintf('Almost ready...\n', i); + end + end + + if flag.save_fig == true + name = [ seed '_circularGraph_' conn_metric '_' band_name '_' seed_target '_' analysis_type '_' alpha '_' correction '_' duration]; + print(gcf, [flag.dir_fig name '.png'], '-dpng', '-r300'); + close all + end + end +end diff --git a/Visualize_Results/scripts/helper_fun/helper_functionConn_fsaverage.m b/Visualize_Results/scripts/helper_fun/helper_functionConn_fsaverage.m new file mode 100644 index 0000000..642914c --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/helper_functionConn_fsaverage.m @@ -0,0 +1,639 @@ +function results = helper_functionConn_fsaverage(cimec) + +inputfile_conn = cimec.inputfile_conn; +band_name = cimec.band_name; +seed_target = cimec.seed_target; +conn_metric = cimec.conn_metric; +statistics = cimec.statistics; +contrast = cimec.contrast; +flag = cimec.flag; + +%% +addpath(genpath(fullfile(fileparts(fileparts(pwd)), 'helper_fun'))); + +statistics.contrast = contrast; +statistics.stat = true; +statistics.f_stats = @signrank; +statistics.f_correction = @fdr_bh; + +ROIs.do = true; % fixed + +colorlim = true; + +%% +% fsaverage +inflated = true; +if strcmp(statistics.analysis_type, 'exploratory') + smoothing = '70'; % percent +else + smoothing = '100'; % percent +end +inputfile_fs = fullfile(load_path('inputfile_fs'), 'fsaverage.mat'); +inputfile_fs_inflated = fullfile(load_path('inputfile_fs'), sprintf('fsaverage_inflated_%s.mat', smoothing)); + +% seed to target +dorsal_R = { 'R_V6_ROI R', 'R_V6A_ROI R', 'R_V7_ROI R', 'R_IPS1_ROI R', 'R_IP1_ROI R', 'R_MIP_ROI R', ... + 'R_VIP_ROI R', 'R_LIPd_ROI R', 'R_LIPv_ROI R', 'R_7AL_ROI R', 'R_7PC_ROI R', 'R_7PL_ROI R', ... + 'R_V3A_ROI R', 'R_V3B_ROI R', 'R_V3CD_ROI R', 'R_LO3_ROI R', 'R_MT_ROI R'}; +ventral_R = {'R_V8_ROI R', 'R_VMV1_ROI R', 'R_VMV2_ROI R', 'R_VMV3_ROI R', 'R_PHT_ROI R', 'R_PH_ROI R', ... + 'R_TE1a_ROI R', 'R_TE1m_ROI R', 'R_TE1p_ROI R', 'R_TE2a_ROI R', 'R_TE2p_ROI R', ... + 'R_TF_ROI R', 'R_TGv_ROI R', 'R_FFC_ROI R', 'R_PIT_ROI R', 'R_VVC_ROI R'}; + +dorsal_L = { 'L_V6_ROI L', 'L_V6A_ROI L', 'L_V7_ROI L', 'L_IPS1_ROI L', 'L_IP1_ROI L', 'L_MIP_ROI L', ... + 'L_VIP_ROI L', 'L_LIPd_ROI L', 'L_LIPv_ROI L', 'L_7AL_ROI L', 'L_7PC_ROI L', 'L_7PL_ROI L', ... + 'L_V3A_ROI L', 'L_V3B_ROI L', 'L_V3CD_ROI L', 'L_LO3_ROI L', 'L_MT_ROI L'}; +ventral_L = {'L_V8_ROI L', 'L_VMV1_ROI L', 'L_VMV2_ROI L', 'L_VMV3_ROI L', 'L_PHT_ROI L', 'L_PH_ROI L', ... + 'L_TE1a_ROI L', 'L_TE1m_ROI L', 'L_TE1p_ROI L', 'L_TE2a_ROI L', 'L_TE2p_ROI L', ... + 'L_TF_ROI L', 'L_TGv_ROI L', 'L_FFC_ROI L', 'L_PIT_ROI L', 'L_VVC_ROI L'}; + + +if strcmp(seed_target, 'right-right') + between = {'R_IFJa_ROI R', 'R_IFJp_ROI R', 'R_FEF_ROI R'}; % 1-3 & 2-3 + ROIs.name = [dorsal_R ventral_R]; +elseif strcmp(seed_target, 'right-left') + between = {'R_IFJa_ROI R', 'R_IFJp_ROI R', 'R_FEF_ROI R'}; % 1-3 & 2-3 + ROIs.name = [dorsal_L ventral_L]; +elseif strcmp(seed_target, 'left-left') + between = {'L_IFJa_ROI L', 'L_IFJp_ROI L', 'L_FEF_ROI L'}; % 1-3 & 2-3 + ROIs.name = [dorsal_L ventral_L]; +elseif strcmp(seed_target, 'left-right') + between = {'L_IFJa_ROI L', 'L_IFJp_ROI L', 'L_FEF_ROI L'}; % 1-3 & 2-3 + ROIs.name = [dorsal_R ventral_R]; +end + + + +%% Inputs +cimec = struct; +cimec.between = between; +cimec.statistics = statistics; +cimec.band_name = band_name; +cimec.ROIs = ROIs; +cimec.colorlim = colorlim; +cimec.inputfile_conn = inputfile_conn; +cimec.inputfile_fs = inputfile_fs; +cimec.inflated = inflated; +cimec.inputfile_fs_inflated = inputfile_fs_inflated; +cimec.flag = flag; +cimec.conn_metric = conn_metric; +cimec.seed_target = seed_target; + +results = visualize_connectivity(cimec); + +%% Function + function results = visualize_connectivity(cimec) + + %% Info about the data + inputfile_conn = cimec.inputfile_conn; % e.g. '...\HCP_Subjects\Outputs\connectivity_results\175237_connectivityResults__2s'; + inputfile_fs = cimec.inputfile_fs; % e.g. '...\HCP_Subjects\Outputs\fsaverage_anatomy\175237_fsaverage'; + + contrast = cimec.statistics.contrast; % contrast map or not + between = cimec.between; % {'R_IFJa_ROI R', 'R_IFJp_ROI R', 'R_FEF_ROI R'}; % 1-2 & 1-3 + + stat = cimec.statistics.stat; + f_stats = cimec.statistics.f_stats; % e.g. @signrank + alpha = cimec.statistics.alpha; % e.g. 0.05 + f_correction = cimec.statistics.f_correction; % correction method e.g. @fdr_bh + corrected = cimec.statistics.corrected; % e.g. true + + ROIs = cimec.ROIs; + + band_name = cimec.band_name; % e.g. 'alpha'; + colorlim = cimec.colorlim; % e.g. true; + + inflated = cimec.inflated; % e.g. true + inputfile_fs_inflated = cimec.inputfile_fs_inflated; + + flag = cimec.flag; + conn_metric = cimec.conn_metric; + seed_target = cimec.seed_target; + + %% Frequency bands + band_names = { 'delta' , 'theta' , 'alpha' , 'beta' , 'gamma' }; + freq_bands = [ 1 4 ; 4 8 ; 8 13 ; 13 30 ; 30 100 ]; + + %% Load data + % avgFreq_icoh + load(inputfile_conn, 'group_results'); + conn_results = group_results; + + % fsaverage + load(inputfile_fs, 'fsaverage') + + %% Prepare the anatomy + % Parcellation info (Brainstorm -> Fieldtrip) + tess_cortex_pial = fsaverage; + + if inflated == true + fs_inflated = load(inputfile_fs_inflated); + tess_cortex_pial.Vertices = fs_inflated.Vertices; + tess_cortex_pial.Faces = fs_inflated.Faces; + end + + Scouts = tess_cortex_pial.Atlas(end-1).Scouts; + Scouts_Vertices = {Scouts.Vertices}'; + Scouts_Label = {Scouts.Label}'; + + parcellation = zeros(length([tess_cortex_pial.Vertices]), 1); + for idx = 1:length(Scouts_Vertices) + per_parcel = Scouts_Vertices{idx}; + + for ii = 1:length(per_parcel) + parcellation(per_parcel(ii)) = idx; + end + end + + brainordinate.parcellation = parcellation; + brainordinate.parcellationlabel = Scouts_Label; + brainordinate.pos = tess_cortex_pial.Vertices; + brainordinate.tri = tess_cortex_pial.Faces; + brainordinate.brainstructure = tess_cortex_pial.SulciMap+1; + brainordinate.brainstructurelabel = {'CORTEX_LEFT', 'CORTEX_RIGHT'}'; + + %% Prepare the connectivity matrix + + conn_mat = struct; + conn_mat.label = Scouts_Label; + conn_mat.band_names = band_name; + conn_mat.brainordinate = brainordinate; + + idx_band = find(strcmp(band_names, band_name)); + tmp_data = group_results(idx_band).data_perSubject; + tmp_data(isnan(tmp_data)) = 1; + + %% seeds + idx_IFJa = find(strcmp(conn_mat.label, between{1})); + idx_IFJp = find(strcmp(conn_mat.label, between{2})); + idx_FEF = find(strcmp(conn_mat.label, between{3})); + + IFJa(:, :) = tmp_data(idx_IFJa,:,:); + IFJp(:, :) = tmp_data(idx_IFJp,:,:); + FEF(:, :) = tmp_data(idx_FEF,:,:); + + %% mean values for the statistical test for 'not contrast' + stat_against = squeeze(mean(mean(tmp_data))); + + %% ROIs/exploratory + if strcmp(statistics.analysis_type, 'exploratory') + all_labels = conn_mat.label'; + if or(strcmp(seed_target, 'right-left'), strcmp(seed_target, 'left-left')) + ROIs.name = all_labels(1:180); + elseif or(strcmp(seed_target, 'right-right'), strcmp(seed_target, 'left-right')) + ROIs.name = all_labels(181:360); + end + end + + nROIs = length(ROIs.name); + for ii = 1:nROIs + ROIs.idx(ii) = find(strcmp(conn_mat.label, ROIs.name{ii})); + end + + mask = zeros(360,1); + mask(ROIs.idx) = 1; + + %% normality test + % stat_IFJa = IFJa; + % stat_IFJa(stat_IFJa(:,1)==1) = NaN; + % stat_FEF = FEF; + % stat_FEF(stat_FEF(:,1)==1) = NaN; + % + % + % for ii = 1:nROIs + % figure + % % histogram(atanh(stat_IFJa(ROIs.idx(ii),:)-stat_FEF(ROIs.idx(ii),:))); + % % histogram(stat_IFJa(ROIs.idx(ii),:)-stat_FEF(ROIs.idx(ii),:)); + % % histogram(stat_IFJa(ROIs.idx(ii),:)); + % histogram(atanh(stat_IFJa(ROIs.idx(ii),:))); + % end + % + % figure + % hist(atanh(stat_IFJa(ROIs.idx(33),:)-stat_FEF(ROIs.idx(33),:))); + % figure + % hist(stat_IFJa(ROIs.idx(33),:)-stat_FEF(ROIs.idx(33),:)); + + + %% + if contrast == true + if stat == true + + %% Inferential statistics + %% Wilcoxon signed rank test + % IFJa vs FEF + for ii = 1:nROIs + pair_1 = IFJa(ROIs.idx(ii),:)'; + pair_2 = FEF(ROIs.idx(ii),:)'; + + [P, H, ~] = f_stats(pair_1, pair_2, 'Alpha', alpha, 'Tail', 'both'); + + H_Wilcoxon_IFJa(ii) = H; + P_Wilcoxon_IFJa(ii) = P; + end + + % IFJp vs FEF + for ii = 1:nROIs + pair_1 = IFJp(ROIs.idx(ii),:)'; + pair_2 = FEF(ROIs.idx(ii),:)'; + + [P, H, ~] = f_stats(pair_1, pair_2, 'Alpha', alpha, 'Tail', 'both'); + + H_Wilcoxon_IFJp(ii) = H; + P_Wilcoxon_IFJp(ii) = P; + end + + %% FDR correction + if corrected == false + % clear NaNs + H_Wilcoxon_IFJa(isnan(H_Wilcoxon_IFJa)) = 0; + P_Wilcoxon_IFJa(isnan(P_Wilcoxon_IFJa)) = 0; + H_Wilcoxon_IFJp(isnan(H_Wilcoxon_IFJp)) = 0; + P_Wilcoxon_IFJp(isnan(P_Wilcoxon_IFJp)) = 0; + + % z-score + % for two-tailed: abs(norminv(p/2)); + Z_Wilcoxon_IFJa = abs(norminv(P_Wilcoxon_IFJa/2)); + Z_Wilcoxon_IFJp = abs(norminv(P_Wilcoxon_IFJp/2)); + + group_IFJa_contrast = mean(FEF-IFJa, 2); + right_group_IFJa = group_IFJa_contrast(ROIs.idx); + right_group_IFJa(~H_Wilcoxon_IFJa) = 0; + group_IFJa = zeros(360,1); + group_IFJa(ROIs.idx) = Z_Wilcoxon_IFJa'.*sign(right_group_IFJa); + stats_IFJa = group_IFJa; + + group_IFJp_contrast = mean(FEF-IFJp, 2); + right_group_IFJp = group_IFJp_contrast(ROIs.idx); + right_group_IFJp(~H_Wilcoxon_IFJp) = 0; + group_IFJp = zeros(360,1); + group_IFJp(ROIs.idx) = Z_Wilcoxon_IFJp'.*sign(right_group_IFJp); + stats_IFJp = group_IFJp; + + elseif corrected == true + % FDR correction + [h_adj_IFJa, ~, ~, p_adj_IFJa] = f_correction(P_Wilcoxon_IFJa(~isnan(P_Wilcoxon_IFJa)), alpha, 'pdep', 'yes'); + [h_adj_IFJp, ~, ~, p_adj_IFJp] = f_correction(P_Wilcoxon_IFJp(~isnan(P_Wilcoxon_IFJp)), alpha, 'pdep', 'yes'); + + % z-score + % for two-tailed: abs(norminv(p/2)); + z_adj_IFJa = abs(norminv(p_adj_IFJa/2)); + z_adj_IFJp = abs(norminv(p_adj_IFJp/2)); + + % mask connectivity values + group_IFJa_contrast = mean(FEF-IFJa, 2); + right_group_IFJa = group_IFJa_contrast(ROIs.idx); + right_group_IFJa(~h_adj_IFJa) = 0; + group_IFJa = zeros(360,1); + group_IFJa(ROIs.idx) = z_adj_IFJa'.*sign(right_group_IFJa); + stats_IFJa = group_IFJa; + + group_IFJp_contrast = mean(FEF-IFJp, 2); + right_group_IFJp = group_IFJp_contrast(ROIs.idx); + right_group_IFJp(~h_adj_IFJp) = 0; + group_IFJp = zeros(360,1); + group_IFJp(ROIs.idx) = z_adj_IFJp'.*sign(right_group_IFJp); + stats_IFJp = group_IFJp; + end + + %% Final matrix (360x360) + tmp_data = diag(ones(360,1)); + tmp_data(idx_IFJa,:) = stats_IFJa; + tmp_data(idx_IFJp,:) = stats_IFJp; + tmp_data(:,idx_IFJa) = stats_IFJa'; + tmp_data(:,idx_IFJp) = stats_IFJp'; + + else + tmp_data = diag(ones(360,1)); + tmp_IFJa = mean(IFJa-FEF, 2); + tmp_IFJp = mean(IFJp-FEF, 2); + + tmp_IFJa(~mask) = 0; + tmp_IFJp(~mask) = 0; + + tmp_data(idx_IFJa,:) = tmp_IFJa; + tmp_data(idx_IFJp,:) = tmp_IFJp; + tmp_data(:,idx_IFJa) = tmp_IFJa'; + tmp_data(:,idx_IFJp) = tmp_IFJp'; + end + + else % not a contrast map + if stat == true + + %% Inferential statistics + %% Wilcoxon signed rank test + % IFJa + for ii = 1:nROIs + pair_1 = IFJa(ROIs.idx(ii),:)'; + + [P, H, ~] = f_stats(pair_1, stat_against, 'Alpha', alpha, 'Tail', 'right'); + + H_Wilcoxon_IFJa(ii) = H; + P_Wilcoxon_IFJa(ii) = P; + end + % IFJp + for ii = 1:nROIs + pair_1 = IFJp(ROIs.idx(ii),:)'; + + [P, H, ~] = f_stats(pair_1, stat_against, 'Alpha', alpha, 'Tail', 'right'); + + H_Wilcoxon_IFJp(ii) = H; + P_Wilcoxon_IFJp(ii) = P; + end + % FEF + for ii = 1:nROIs + pair_1 = FEF(ROIs.idx(ii),:)'; + + [P, H, ~] = f_stats(pair_1, stat_against, 'Alpha', alpha, 'Tail', 'right'); + + H_Wilcoxon_FEF(ii) = H; + P_Wilcoxon_FEF(ii) = P; + end + %% FDR correction + if corrected == false + % clear NaNs + H_Wilcoxon_IFJa(isnan(H_Wilcoxon_IFJa)) = 0; + H_Wilcoxon_IFJp(isnan(H_Wilcoxon_IFJp)) = 0; + H_Wilcoxon_FEF(isnan(H_Wilcoxon_FEF)) = 0; + P_Wilcoxon_IFJa(isnan(P_Wilcoxon_IFJa)) = 0; + P_Wilcoxon_IFJp(isnan(P_Wilcoxon_IFJp)) = 0; + P_Wilcoxon_FEF(isnan(P_Wilcoxon_FEF)) = 0; + + % z-score + % for one-tailed: abs(norminv(p)); + Z_Wilcoxon_IFJa = abs(norminv(P_Wilcoxon_IFJa)); + Z_Wilcoxon_IFJp = abs(norminv(P_Wilcoxon_IFJp)); + Z_Wilcoxon_FEF = abs(norminv(P_Wilcoxon_FEF)); + + % mask connectivity values + group_IFJa_notContrast = mean(IFJa, 2); + right_group_IFJa = group_IFJa_notContrast(ROIs.idx); + right_group_IFJa(~H_Wilcoxon_IFJa) = 0; + group_IFJa = zeros(360,1); + group_IFJa(ROIs.idx) = Z_Wilcoxon_IFJa'.*sign(right_group_IFJa); + stats_IFJa = group_IFJa; + + group_IFJp_notContrast = mean(IFJp, 2); + right_group_IFJp = group_IFJp_notContrast(ROIs.idx); + right_group_IFJp(~H_Wilcoxon_IFJp) = 0; + group_IFJp = zeros(360,1); + group_IFJp(ROIs.idx) = Z_Wilcoxon_IFJp'.*sign(right_group_IFJp); + stats_IFJp = group_IFJp; + + group_FEF_notContrast = mean(FEF, 2); + right_group_FEF = group_FEF_notContrast(ROIs.idx); + right_group_FEF(~H_Wilcoxon_FEF) = 0; + group_FEF = zeros(360,1); + group_FEF(ROIs.idx) = Z_Wilcoxon_FEF'.*sign(right_group_FEF); + stats_FEF = group_FEF; + + elseif corrected == true + % FDR correction + [h_adj_IFJa, ~, ~, p_adj_IFJa] = f_correction(P_Wilcoxon_IFJa(~isnan(P_Wilcoxon_IFJa)), alpha, 'pdep', 'yes'); + [h_adj_IFJp, ~, ~, p_adj_IFJp] = f_correction(P_Wilcoxon_IFJp(~isnan(P_Wilcoxon_IFJp)), alpha, 'pdep', 'yes'); + [h_adj_FEF, ~, ~, p_adj_FEF] = f_correction(P_Wilcoxon_FEF(~isnan(P_Wilcoxon_FEF)) , alpha, 'pdep', 'yes'); + + % z-score + % for one-tailed: abs(norminv(p)); + z_adj_IFJa = abs(norminv(p_adj_IFJa)); + z_adj_IFJp = abs(norminv(p_adj_IFJp)); + z_adj_FEF = abs(norminv(p_adj_FEF)); + + % mask connectivity values + group_IFJa_notContrast = mean(IFJa, 2); + right_group_IFJa = group_IFJa_notContrast(ROIs.idx); + right_group_IFJa(~h_adj_IFJa) = 0; + group_IFJa = zeros(360,1); + group_IFJa(ROIs.idx) = z_adj_IFJa'.*sign(right_group_IFJa); + stats_IFJa = group_IFJa; + + group_IFJp_notContrast = mean(IFJp, 2); + right_group_IFJp = group_IFJp_notContrast(ROIs.idx); + right_group_IFJp(~h_adj_IFJp) = 0; + group_IFJp = zeros(360,1); + group_IFJp(ROIs.idx) = z_adj_IFJp'.*sign(right_group_IFJp); + stats_IFJp = group_IFJp; + + group_FEF_notContrast = mean(FEF, 2); + right_group_FEF = group_FEF_notContrast(ROIs.idx); + right_group_FEF(~h_adj_FEF) = 0; + group_FEF = zeros(360,1); + group_FEF(ROIs.idx) = z_adj_FEF'.*sign(right_group_FEF); + stats_FEF = group_FEF; + end + + %% Final matrix (360x360) + tmp_data = zeros(360,360); + tmp_data(idx_IFJa,:) = stats_IFJa; + tmp_data(idx_IFJp,:) = stats_IFJp; + tmp_data(idx_FEF,:) = stats_FEF; + tmp_data(:,idx_IFJa) = stats_IFJa'; + tmp_data(:,idx_IFJp) = stats_IFJp'; + tmp_data(:,idx_FEF) = stats_FEF'; + + else + tmp_data = zeros(360,360); + tmp_IFJa = mean(IFJa, 2); + tmp_IFJp = mean(IFJp, 2); + tmp_FEF = mean(FEF, 2); + + tmp_IFJa(~mask) = 0; + tmp_IFJp(~mask) = 0; + tmp_FEF(~mask) = 0; + + tmp_data(idx_IFJa,:) = tmp_IFJa; + tmp_data(idx_IFJp,:) = tmp_IFJp; + tmp_data(idx_FEF,:) = tmp_FEF; + tmp_data(:,idx_IFJa) = tmp_IFJa'; + tmp_data(:,idx_IFJp) = tmp_IFJp'; + tmp_data(:,idx_FEF) = tmp_FEF'; + end + end + + conn_mat.cohspctrm = tmp_data; + + % output + results.conn_matrix = tmp_data; + results.contrast = contrast; + results.helper.idx_IFJa = idx_IFJa; + results.helper.idx_IFJp = idx_IFJp; + results.helper.idx_FEF = idx_FEF; + results.helper.ROIs = ROIs; + results.helper.brainordinate = brainordinate; + results.label = conn_mat.label; + if exist('h_adj_IFJa','var'); results.stat.h_adj_IFJa = h_adj_IFJa; end + if exist('h_adj_IFJp','var'); results.stat.h_adj_IFJp = h_adj_IFJp; end + if exist('h_adj_FEF','var'); results.stat.h_adj_FEF = h_adj_FEF; end + results.stat.analysis_type = cimec.statistics.analysis_type; + + results.perSubject.IFJa = IFJa; + results.perSubject.IFJp = IFJp; + results.perSubject.FEF = FEF; + + results.cimec = cimec; + + %% Visualize using the helper function + % for colormap scale + if contrast == false && stat == false + tmp_data(logical(eye(360))) = NaN; + tmp_data = tmp_data(idx_seed1,:); + minv = min(tmp_data); + maxv = max(tmp_data); + elseif contrast == true && stat == false + tmp_data(logical(eye(360))) = NaN; + tmp_data = tmp_data(idx_seed1,:); + + tmp_data(tmp_data > 0.85) = NaN; + tmp_data(tmp_data < -0.85) = NaN; + tmp_minv = min(tmp_data); + tmp_maxv = max(tmp_data); + + minv = -1*round((abs(tmp_minv) + tmp_maxv)/2, 3); + maxv = round((abs(tmp_minv) + tmp_maxv)/2, 3); + elseif contrast == true && stat == true + minv = -4; + maxv = 4; + elseif contrast == false && stat == true + minv = 0; + maxv = 6; + end + colorlim_value = [minv maxv]; + + if flag.fsaverage == true + % in order to run the function + conn_mat.brainordinate.pos = conn_mat.brainordinate.pos*1000; + + % To make black the areas outside ROI + if and(contrast == true, strcmp(statistics.analysis_type, 'ROI')) + cbar = cmap_rbw; + n = size(cbar, 1); + up = cbar(1:n/2,:); + middle = [0.5 0.5 0.5]; + down = cbar(n/2+1:n,:); + set(0, 'DefaultFigureColormap', [up;middle;down]); + + conn_mat.cohspctrm(ROIs.idx(~h_adj_IFJa),idx_IFJa) = 0.04; + conn_mat.cohspctrm(ROIs.idx(~h_adj_IFJp),idx_IFJp) = 0.04; + conn_mat.cohspctrm(idx_IFJa,ROIs.idx(~h_adj_IFJa)) = 0.04; + conn_mat.cohspctrm(idx_IFJp,ROIs.idx(~h_adj_IFJp)) = 0.04; + + elseif and(contrast == false, strcmp(statistics.analysis_type, 'ROI')) + cbar = hot; + base1 = [1 1 1]; + base2 = [0.25 0.25 0.25]; + set(0, 'DefaultFigureColormap', [base1; base2; cbar]); + + conn_mat.cohspctrm(ROIs.idx(~h_adj_IFJa),idx_IFJa) = -0.06; + conn_mat.cohspctrm(ROIs.idx(~h_adj_IFJp),idx_IFJp) = -0.06; + conn_mat.cohspctrm(ROIs.idx(~h_adj_FEF),idx_FEF) = -0.06; + conn_mat.cohspctrm(idx_IFJa,ROIs.idx(~h_adj_IFJa)) = -0.06; + conn_mat.cohspctrm(idx_IFJp,ROIs.idx(~h_adj_IFJp)) = -0.06; + conn_mat.cohspctrm(idx_FEF,ROIs.idx(~h_adj_FEF)) = -0.06; + + minv = -0.06; + maxv = 6; + colorlim_value = [minv maxv]; + elseif contrast == false + set(0, 'DefaultFigureColormap', hot); + end + + %% Save the figure + if flag.save_fig == true + + if seed_target(1) == 'r' + seed = 'IFJa'; + cimec.seed = seed; + if strcmp(statistics.analysis_type, 'ROI') + pos2d_IFJa_R = [52.5323 -69.3638 71.0650; 52.5323 69.4183 71.0650]; + else + pos2d_IFJa_R = [49.7454 -69.6590 70.1540; 49.7454 69.1230 70.1540]; + end + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_IFJa_R); + helper_save_figure(cimec); + close all; + + seed = 'IFJp'; + cimec.seed = seed; + if strcmp(statistics.analysis_type, 'ROI') + pos2d_IFJp_R = [45.6139 -69.3638 76.8706 ; 45.6139 69.4183 76.8706]; + else + pos2d_IFJp_R = [44.7673 -69.6590 76.5544 ; 44.7673 69.1230 76.5544]; + end + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_IFJp_R); + helper_save_figure(cimec); + close all; + + if contrast == false + seed = 'FEF'; + cimec.seed = seed; + if strcmp(statistics.analysis_type, 'ROI') + pos2d_FEF_R = [25.8639 -69.3638 96.3111; 25.8639 69.4183 96.3111]; + else + pos2d_FEF_R = [28.0551 -69.6590 101.8004 ; 28.0551 69.1230 101.8004]; + end + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_FEF_R); + helper_save_figure(cimec); + close all; + end + + elseif seed_target(1) == 'l' + seed = 'IFJa'; + cimec.seed = seed; + if strcmp(statistics.analysis_type, 'ROI') + pos2d_IFJa_L = [50.5105 69.4183 69.5634 ; 53.1465 -69.3638 68.2860]; + else + pos2d_IFJa_L = [47.6119 69.1230 68.3761 ; 47.6119 -69.6590 68.3761]; + end + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_IFJa_L); + helper_save_figure(cimec); + close all; + + seed = 'IFJp'; + cimec.seed = seed; + if strcmp(statistics.analysis_type, 'ROI') + pos2d_IFJp_L = [43.6193 69.4183 80.0619 ; 43.6193 -69.3638 80.0619]; + else + pos2d_IFJp_L = [40.8559 69.1230 76.9099 ; 40.8559 -69.6590 76.9099]; + end + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_IFJp_L); + helper_save_figure(cimec); + close all; + + if contrast == false + seed = 'FEF'; + cimec.seed = seed; + if strcmp(statistics.analysis_type, 'ROI') + pos2d_FEF_L = [30.9891 69.4183 105.2148 ; 34.2627 -69.3638 98.8225]; + else + pos2d_FEF_L = [26.2772 69.1230 101.4449 ; 26.2772 -69.6590 101.4449]; + end + tutorial_nwa_connectivityviewer_save_figures(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_FEF_L); + helper_save_figure(cimec); + close all; + end + end + + %% Show the figure + elseif colorlim == true + tutorial_nwa_connectivityviewer(conn_mat, 'cohspctrm', colorlim_value, seed_target, flag.contours); + end + end + + %% circularGraph + if flag.circular_graph == true + + results.contrast = contrast; + results.helper.stat.analysis_type = cimec.statistics.analysis_type; + results.helper.stat.h_IFJa_all = h_adj_IFJa; + results.helper.stat.h_IFJp_all = h_adj_IFJp; + if exist('h_adj_FEF','var'); results.helper.stat.h_FEF_all = h_adj_FEF; end + + fontSize = cimec.flag.font_size; + save_fig = cimec.flag.save_fig; + results.cimec = cimec; + + flag.save_fig = save_fig; + flag.collapseTargets = false; + flag.contrast = contrast; + helper_functionConn_circularGraph(results, fontSize, flag); + end + end + +end diff --git a/Visualize_Results/scripts/helper_fun/helper_functionConn_fsaverage_freqCollapsed.m b/Visualize_Results/scripts/helper_fun/helper_functionConn_fsaverage_freqCollapsed.m new file mode 100644 index 0000000..d30c03d --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/helper_functionConn_fsaverage_freqCollapsed.m @@ -0,0 +1,211 @@ +function results = helper_functionConn_fsaverage_freqCollapsed(cimec) + +band_name = cimec.band_name; +flag = cimec.flag; +conn_metric = cimec.conn_metric; +contrast = cimec.contrast; +analysis_type = cimec.statistics.analysis_type; +seed_target = cimec.seed_target; +alpha_value = cimec.statistics.alpha; +duration = cimec.flag.duration; +correction = cimec.statistics.correction; +corrected = cimec.statistics.corrected; + +%% All freq bands +cimec.flag.circular_graph = false; +cimec.flag.fsaverage = false; +cimec.band_name = 'delta'; +delta = helper_functionConn_fsaverage(cimec); +delta.conn_matrix(delta.conn_matrix == 1) = NaN; +cimec.band_name = 'theta'; +theta = helper_functionConn_fsaverage(cimec); +theta.conn_matrix(theta.conn_matrix == 1) = NaN; +cimec.band_name = 'alpha'; +alpha = helper_functionConn_fsaverage(cimec); +alpha.conn_matrix(alpha.conn_matrix == 1) = NaN; +cimec.band_name = 'beta'; +beta = helper_functionConn_fsaverage(cimec); +beta.conn_matrix(beta.conn_matrix == 1) = NaN; +cimec.band_name = 'gamma'; +gamma = helper_functionConn_fsaverage(cimec); +gamma.conn_matrix(gamma.conn_matrix == 1) = NaN; + +idx_IFJa = delta.helper.idx_IFJa; +idx_IFJp = delta.helper.idx_IFJp; +idx_FEF = delta.helper.idx_FEF; +ROIs.idx = delta.helper.ROIs.idx; + +h_IFJa_all = delta.stat.h_adj_IFJa + theta.stat.h_adj_IFJa + alpha.stat.h_adj_IFJa + beta.stat.h_adj_IFJa + gamma.stat.h_adj_IFJa; +h_IFJp_all = delta.stat.h_adj_IFJp + theta.stat.h_adj_IFJp + alpha.stat.h_adj_IFJp + beta.stat.h_adj_IFJp + gamma.stat.h_adj_IFJp; +if contrast == false + h_FEF_all = delta.stat.h_adj_FEF + theta.stat.h_adj_FEF + alpha.stat.h_adj_FEF + beta.stat.h_adj_FEF + gamma.stat.h_adj_FEF; +end + + +conn_matrix = zeros(360,360); +tmp_mat = (delta.conn_matrix + theta.conn_matrix + alpha.conn_matrix + beta.conn_matrix + gamma.conn_matrix); +tmp_mat(idx_IFJa, ROIs.idx) = tmp_mat(idx_IFJa, ROIs.idx)./h_IFJa_all; +tmp_mat(idx_IFJp, ROIs.idx) = tmp_mat(idx_IFJp, ROIs.idx)./h_IFJp_all; +if contrast == false + tmp_mat(idx_FEF, ROIs.idx) = tmp_mat(idx_FEF, ROIs.idx)./h_FEF_all; +end +tmp_mat(isinf(tmp_mat)) = 0; +tmp_mat(isnan(tmp_mat)) = 0; + +conn_matrix(:,idx_IFJa) = tmp_mat(idx_IFJa,:); +conn_matrix(:,idx_IFJp) = tmp_mat(idx_IFJp,:); +conn_matrix(idx_IFJa,:) = tmp_mat(idx_IFJa,:); +conn_matrix(idx_IFJp,:) = tmp_mat(idx_IFJp,:); +if contrast == false + conn_matrix(:,idx_FEF) = tmp_mat(idx_FEF,:); + conn_matrix(idx_FEF,:) = tmp_mat(idx_FEF,:); +end +conn_matrix(idx_IFJa, idx_IFJa) = 0; +conn_matrix(idx_IFJp, idx_IFJp) = 0; +conn_matrix(idx_FEF, idx_FEF) = 0; + + +%% Visualize using the helper function +if flag.fsaverage == true + + allFreq.cohspctrm = conn_matrix; + allFreq.brainordinate = delta.helper.brainordinate; + % in order to run the function + allFreq.brainordinate.pos = allFreq.brainordinate.pos*1000; + if contrast == true + colorlim_value = [-4 4]; + else + colorlim_value = [0 6]; + end + + %% To make black the areas outside ROI + if and(contrast == true, strcmp(analysis_type, 'ROI')) + cbar = cmap_rbw; + n = size(cbar, 1); + up = cbar(1:n/2,:); + middle = [0.5 0.5 0.5]; + down = cbar(n/2+1:n,:); + set(0, 'DefaultFigureColormap', [up;middle;down]); + + allFreq.cohspctrm(ROIs.idx(~h_IFJa_all),idx_IFJa) = 0.04; + allFreq.cohspctrm(ROIs.idx(~h_IFJp_all),idx_IFJp) = 0.04; + allFreq.cohspctrm(idx_IFJa,ROIs.idx(~h_IFJa_all)) = 0.04; + allFreq.cohspctrm(idx_IFJp,ROIs.idx(~h_IFJp_all)) = 0.04; + + elseif and(contrast == false, strcmp(analysis_type, 'ROI')) + cbar = hot; + base1 = [1 1 1]; + base2 = [0.25 0.25 0.25]; + set(0, 'DefaultFigureColormap', [base1; base2; cbar]); + + allFreq.cohspctrm(ROIs.idx(~h_IFJa_all),idx_IFJa) = -0.06; + allFreq.cohspctrm(ROIs.idx(~h_IFJp_all),idx_IFJp) = -0.06; + allFreq.cohspctrm(ROIs.idx(~h_FEF_all),idx_FEF) = -0.06; + allFreq.cohspctrm(idx_IFJa,ROIs.idx(~h_IFJa_all)) = -0.06; + allFreq.cohspctrm(idx_IFJp,ROIs.idx(~h_IFJp_all)) = -0.06; + allFreq.cohspctrm(idx_FEF,ROIs.idx(~h_FEF_all)) = -0.06; + + minv = -0.06; + maxv = 6; + colorlim_value = [minv maxv]; + elseif contrast == false + set(0, 'DefaultFigureColormap', hot); + end + + %% save the figure + if flag.save_fig == true + cimec = struct; + cimec.conn_metric = conn_metric; + cimec.band_name = band_name; + cimec.seed_target = seed_target; + cimec.statistics.analysis_type = analysis_type; + cimec.statistics.alpha = alpha_value; + cimec.statistics.corrected = corrected; + cimec.statistics.correction = correction; + cimec.flag.duration = duration; + cimec.flag.dir_fig = flag.dir_fig; + + if seed_target(1) == 'r' + seed = 'IFJa'; + cimec.seed = seed; + pos2d_IFJa_R = [52.5323 -69.3638 71.0650; 52.5323 69.4183 71.0650]; + tutorial_nwa_connectivityviewer_save_figures(allFreq, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_IFJa_R); + helper_save_figure(cimec); + close all; + + seed = 'IFJp'; + cimec.seed = seed; + pos2d_IFJp_R = [45.6139 -69.3638 76.8706 ; 45.6139 69.4183 76.8706]; + tutorial_nwa_connectivityviewer_save_figures(allFreq, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_IFJp_R); + helper_save_figure(cimec); + close all; + + if contrast == false + seed = 'FEF'; + cimec.seed = seed; + pos2d_FEF_R = [25.8639 -69.3638 96.3111; 25.8639 69.4183 96.3111]; + tutorial_nwa_connectivityviewer_save_figures(allFreq, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_FEF_R); + helper_save_figure(cimec); + close all; + end + + elseif seed_target(1) == 'l' + seed = 'IFJa'; + cimec.seed = seed; + pos2d_IFJa_L = [50.5105 69.4183 69.5634 ; 53.1465 -69.3638 68.2860]; + tutorial_nwa_connectivityviewer_save_figures(allFreq, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_IFJa_L); + helper_save_figure(cimec); + close all; + + seed = 'IFJp'; + cimec.seed = seed; + pos2d_IFJp_L = [43.6193 69.4183 80.0619 ; 43.6193 -69.3638 80.0619]; + tutorial_nwa_connectivityviewer_save_figures(allFreq, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_IFJp_L); + helper_save_figure(cimec); + close all; + + if contrast == false + seed = 'FEF'; + cimec.seed = seed; + pos2d_FEF_L = [30.9891 69.4183 105.2148 ; 34.2627 -69.3638 98.8225]; + tutorial_nwa_connectivityviewer_save_figures(allFreq, 'cohspctrm', colorlim_value, seed_target, flag.contours, pos2d_FEF_L); + helper_save_figure(cimec); + close all; + end + end + else + tutorial_nwa_connectivityviewer(allFreq, 'cohspctrm', colorlim_value, cimec.seed_target, cimec.flag.contours); + end +end + +results.conn_matrix = conn_matrix; +results.contrast = contrast; +results.label = delta.label; +results.helper.brainordinate = delta.helper.brainordinate; +results.helper.idx_FEF = delta.helper.idx_FEF; +results.helper.idx_IFJa = delta.helper.idx_IFJa; +results.helper.idx_IFJp = delta.helper.idx_IFJp; +results.helper.ROIs = delta.helper.ROIs; +results.stat.analysis_type = analysis_type; +results.helper.band_name = band_name; +results.stat.h_IFJa_all = h_IFJa_all; +results.stat.h_IFJp_all = h_IFJp_all; +if exist('h_FEF_all','var'); results.stat.h_FEF_all = h_FEF_all; end + +results.perSubject.IFJa = (delta.perSubject.IFJa + theta.perSubject.IFJa + alpha.perSubject.IFJa + beta.perSubject.IFJa + gamma.perSubject.IFJa)/5; +results.perSubject.IFJp = (delta.perSubject.IFJp + theta.perSubject.IFJp + alpha.perSubject.IFJp + beta.perSubject.IFJp + gamma.perSubject.IFJp)/5; +results.perSubject.FEF = (delta.perSubject.FEF + theta.perSubject.FEF + alpha.perSubject.FEF + beta.perSubject.FEF + gamma.perSubject.FEF)/5; + +cimec.band_name = band_name; +results.cimec = cimec; + +%% circularGraph +if flag.circular_graph == true + fontSize = flag.font_size; + flag.collapseTargets = false; + flag.contrast = contrast; + + helper_functionConn_circularGraph(results, fontSize, flag); +end + +end diff --git a/Visualize_Results/scripts/helper_fun/helper_save_figure.m b/Visualize_Results/scripts/helper_fun/helper_save_figure.m new file mode 100644 index 0000000..9136322 --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/helper_save_figure.m @@ -0,0 +1,165 @@ +function helper_save_figure(cimec) +%% 'target right hemisphere' + +conn_metric = cimec.conn_metric; +band_name = cimec.band_name; +seed = cimec.seed; +seed_target = cimec.seed_target; +analysis_type = cimec.statistics.analysis_type; +alpha = char(string(cimec.statistics.alpha)); +correction = cimec.statistics.correction; +flag = cimec.flag; +duration = flag.duration; +dir_fig = flag.dir_fig; + +if or(strcmp(seed_target, 'right-right'), strcmp(seed_target, 'left-right')) + + try + if strcmp(flag.script, 'collapseTargets') + seed_target = char(string(join(flag.seed_targets))); + end + catch + fprintf('Almost ready...\n', i); + end + + %% 0 + view(270,360) + xlim([-100 100]) + ylim([-100 0]) + zlim([-20 150]) + set(gcf, 'Position', get(0, 'Screensize')); + + %% 1 + % save the figure + lighting none + material dull + lighting phong + cam1 = camlight('left'); + + % for light + view(360,270) + cam2 = camlight('headlight'); + view(360,360) + cam3 = camlight('left'); + + if strcmp(analysis_type, 'exploratory') + cam4 = camlight('right'); + + %% 1 + number = '1'; + name = [number '_' conn_metric '_' band_name '_' seed '_' seed_target '_' analysis_type '_' alpha '_' correction '_' duration]; + print(gcf,[dir_fig name '.png'],'-dpng','-r250'); + + %% 3 + view(180,360) + cam3 =camlight('headlight'); + number = '3'; + name = [number '_' conn_metric '_' band_name '_' seed '_' seed_target '_' analysis_type '_' alpha '_' correction '_' duration]; + print(gcf,[dir_fig name '.png'],'-dpng','-r250'); + else + %% 1 + number = '1'; + name = [number '_' conn_metric '_' band_name '_' seed '_' seed_target '_' analysis_type '_' alpha '_' correction '_' duration]; + print(gcf,[dir_fig name '.png'],'-dpng','-r250'); + + %% 2 + view(315,5) + number = '2'; + name = [number '_' conn_metric '_' band_name '_' seed '_' seed_target '_' analysis_type '_' alpha '_' correction '_' duration]; + print(gcf,[dir_fig name '.png'],'-dpng','-r250'); + + %% 3 + view(180,360) + cam3 =camlight('headlight'); + number = '3'; + name = [number '_' conn_metric '_' band_name '_' seed '_' seed_target '_' analysis_type '_' alpha '_' correction '_' duration]; + print(gcf,[dir_fig name '.png'],'-dpng','-r250'); + + %% 4 + if strcmp(conn_metric, 'pdc') + view(360,270) + number = '4'; + name = [number '_' conn_metric '_' band_name '_' seed '_' seed_target '_' analysis_type '_' alpha '_' correction '_' duration]; + print(gcf,[dir_fig name '.png'],'-dpng','-r250'); + end + end + +elseif or(strcmp(seed_target, 'right-left'), strcmp(seed_target, 'left-left')) + + try + if strcmp(flag.script, 'collapseTargets') + seed_target = char(string(join(flag.seed_targets))); + end + catch + fprintf('Almost ready...\n', i); + end + + + + %% RUN + %% 0 + view(270,360) + xlim([-100 100]) + ylim([0 100]) + zlim([-20 150]) + set(gcf, 'Position', get(0, 'Screensize')); + + %% 1 + % save the figure + lighting none + material dull + lighting phong + cam1 = camlight('right'); + + % for light + view(360,270) + cam2 = camlight('headlight'); + view(180,360) + cam3 = camlight('right'); + + if strcmp(analysis_type, 'exploratory') + %% 1 + cam4 = camlight('left'); + + number = '1'; + name = [number '_' conn_metric '_' band_name '_' seed '_' seed_target '_' analysis_type '_' alpha '_' correction '_' duration]; + print(gcf,[dir_fig name '.png'],'-dpng','-r250'); + + %% 3 + view(360,360) + cam3 =camlight('headlight'); + number = '3'; + name = [number '_' conn_metric '_' band_name '_' seed '_' seed_target '_' analysis_type '_' alpha '_' correction '_' duration]; + print(gcf,[dir_fig name '.png'],'-dpng','-r250'); + + else + %% 1 + number = '1'; + name = [number '_' conn_metric '_' band_name '_' seed '_' seed_target '_' analysis_type '_' alpha '_' correction '_' duration]; + print(gcf,[dir_fig name '.png'],'-dpng','-r250'); + + %% 2 + view(225,5) + number = '2'; + name = [number '_' conn_metric '_' band_name '_' seed '_' seed_target '_' analysis_type '_' alpha '_' correction '_' duration]; + print(gcf,[dir_fig name '.png'],'-dpng','-r250'); + + %% 3 + view(360,360) + cam3 =camlight('headlight'); + number = '3'; + name = [number '_' conn_metric '_' band_name '_' seed '_' seed_target '_' analysis_type '_' alpha '_' correction '_' duration]; + print(gcf,[dir_fig name '.png'],'-dpng','-r250'); + + %% 4 + if strcmp(conn_metric, 'pdc') + view(360,270) + number = '4'; + name = [number '_' conn_metric '_' band_name '_' seed '_' seed_target '_' analysis_type '_' alpha '_' correction '_' duration]; + print(gcf,[dir_fig name '.png'],'-dpng','-r250'); + end + end +end + +end +%% \ No newline at end of file diff --git a/Visualize_Results/scripts/helper_fun/helper_violinplot.m b/Visualize_Results/scripts/helper_fun/helper_violinplot.m new file mode 100644 index 0000000..5bb02ce --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/helper_violinplot.m @@ -0,0 +1,370 @@ +function helper_violinplot(outputs, conn_metric, band_name, seed_target, duration) + +%% Cite: https://github.com/bastibe/Violinplot-Matlab +idx_FEF_R = 242; +idx_IFJa_R = 251; +idx_IFJp_R = 252; + +idx_FEF_L = 62; +idx_IFJa_L = 71; +idx_IFJp_L = 72; + +if strcmp(band_name, 'delta') + color = [240 50 230]/255; +elseif strcmp(band_name, 'theta') + color = [255 225 25]/255; +elseif strcmp(band_name, 'alpha') + color = [60 180 75]/255; +elseif strcmp(band_name, 'beta') + color = [245 130 48]/255; +elseif strcmp(band_name, 'gamma') + color = [128 128 128]/255; +end + +if strcmp(band_name, 'gamma') + yvalues = [0 0.24]; +else + yvalues = [0 0.14]; +end + +%% 'right-right' +if strcmp(seed_target, 'right-right') + + seed_12 = outputs.perSubject.FEF_12; + seed_21 = outputs.perSubject.FEF_21; + + name = [conn_metric '_' band_name '_' seed_target '_' duration '_' 'R-FEF and R-IFJa']; + legend = {'R-FEF to R-IFJa' 'R-IFJa to R-FEF'}; + labels = repmat(repmat(legend,55,1), 1, 1); + labels = labels(:); + + values = [seed_12(idx_IFJa_R,:)'; seed_21(idx_IFJa_R,:)']; + + figure; vs = violinplot(values, labels, 'ShowMean', true); + + ylabel('Partial Directed Coherence'); ylim(yvalues); + title('Granger Causal Influences') + vs(1).ViolinColor = color; + vs(2).ViolinColor = color; + + print(gcf,[load_path('save_fig') name '.png'],'-dpng','-r300'); + close all; + + + %% + seed_12 = outputs.perSubject.FEF_12; + seed_21 = outputs.perSubject.FEF_21; + + name = [conn_metric '_' band_name '_' seed_target '_' duration '_' 'R-FEF and R-IFJp']; + legend = {'R-FEF to R-IFJp' 'R-IFJp to R-FEF'}; + labels = repmat(repmat(legend,55,1), 1, 1); + labels = labels(:); + + values = [seed_12(idx_IFJp_R,:)'; seed_21(idx_IFJp_R,:)']; + + figure; vs = violinplot(values, labels, 'ShowMean', true); + + ylabel('Partial Directed Coherence'); ylim(yvalues); + title('Granger Causal Influences') + vs(1).ViolinColor = color; + vs(2).ViolinColor = color; + + print(gcf,[load_path('save_fig') name '.png'],'-dpng','-r300'); + close all; + + %% + seed_12 = outputs.perSubject.IFJa_12; + seed_21 = outputs.perSubject.IFJa_21; + + name = [conn_metric '_' band_name '_' seed_target '_' duration '_' 'R-IFJa and R-IFJp']; + legend = {'R-IFJa to R-IFJp' 'R-IFJp to R-IFJa'}; + labels = repmat(repmat(legend,55,1), 1, 1); + labels = labels(:); + + values = [seed_12(idx_IFJp_R,:)'; seed_21(idx_IFJp_R,:)']; + + figure; vs = violinplot(values, labels, 'ShowMean', true); + + ylabel('Partial Directed Coherence'); ylim(yvalues); + title('Granger Causal Influences') + vs(1).ViolinColor = color; + vs(2).ViolinColor = color; + + print(gcf,[load_path('save_fig') name '.png'],'-dpng','-r300'); + close all; + + + + %% 'left-left' +elseif strcmp(seed_target, 'left-left') + + seed_12 = outputs.perSubject.FEF_12; + seed_21 = outputs.perSubject.FEF_21; + + name = [conn_metric '_' band_name '_' seed_target '_' duration '_' 'L-FEF and L-IFJa']; + legend = {'L-FEF to L-IFJa' 'L-IFJa to L-FEF'}; + labels = repmat(repmat(legend,55,1), 1, 1); + labels = labels(:); + + values = [seed_12(idx_IFJa_L,:)'; seed_21(idx_IFJa_L,:)']; + + figure; vs = violinplot(values, labels, 'ShowMean', true); + + ylabel('Partial Directed Coherence'); ylim(yvalues); + title('Granger Causal Influences') + vs(1).ViolinColor = color; + vs(2).ViolinColor = color; + + print(gcf,[load_path('save_fig') name '.png'],'-dpng','-r300'); + close all; + + + %% + seed_12 = outputs.perSubject.FEF_12; + seed_21 = outputs.perSubject.FEF_21; + + name = [conn_metric '_' band_name '_' seed_target '_' duration '_' 'L-FEF and L-IFJp']; + legend = {'L-FEF to L-IFJp' 'L-IFJp to L-FEF'}; + labels = repmat(repmat(legend,55,1), 1, 1); + labels = labels(:); + + values = [seed_12(idx_IFJp_L,:)'; seed_21(idx_IFJp_L,:)']; + + figure; vs = violinplot(values, labels, 'ShowMean', true); + + ylabel('Partial Directed Coherence'); ylim(yvalues); + title('Granger Causal Influences') + vs(1).ViolinColor = color; + vs(2).ViolinColor = color; + + print(gcf,[load_path('save_fig') name '.png'],'-dpng','-r300'); + close all; + + %% + seed_12 = outputs.perSubject.IFJa_12; + seed_21 = outputs.perSubject.IFJa_21; + + name = [conn_metric '_' band_name '_' seed_target '_' duration '_' 'L-IFJa and L-IFJp']; + legend = {'L-IFJa to L-IFJp' 'L-IFJp to L-IFJa'}; + labels = repmat(repmat(legend,55,1), 1, 1); + labels = labels(:); + + values = [seed_12(idx_IFJp_L,:)'; seed_21(idx_IFJp_L,:)']; + + figure; vs = violinplot(values, labels, 'ShowMean', true); + + ylabel('Partial Directed Coherence'); ylim(yvalues); + title('Granger Causal Influences') + vs(1).ViolinColor = color; + vs(2).ViolinColor = color; + + print(gcf,[load_path('save_fig') name '.png'],'-dpng','-r300'); + close all; + + + + + %% 'right-left' +elseif strcmp(seed_target, 'right-left') + + seed_12 = outputs.perSubject.FEF_12; + seed_21 = outputs.perSubject.FEF_21; + + name = [conn_metric '_' band_name '_' seed_target '_' duration '_' 'R-FEF and L-IFJa']; + legend = {'R-FEF to L-IFJa' 'L-IFJa to R-FEF'}; + labels = repmat(repmat(legend,55,1), 1, 1); + labels = labels(:); + + values = [seed_12(idx_IFJa_L,:)'; seed_21(idx_IFJa_L,:)']; + + figure; vs = violinplot(values, labels, 'ShowMean', true); + + ylabel('Partial Directed Coherence'); ylim(yvalues); + title('Granger Causal Influences') + vs(1).ViolinColor = color; + vs(2).ViolinColor = color; + + print(gcf,[load_path('save_fig') name '.png'],'-dpng','-r300'); + close all; + + + %% + seed_12 = outputs.perSubject.FEF_12; + seed_21 = outputs.perSubject.FEF_21; + + name = [conn_metric '_' band_name '_' seed_target '_' duration '_' 'R-FEF and L-IFJp']; + legend = {'R-FEF to L-IFJp' 'L-IFJp to R-FEF'}; + labels = repmat(repmat(legend,55,1), 1, 1); + labels = labels(:); + + values = [seed_12(idx_IFJp_L,:)'; seed_21(idx_IFJp_L,:)']; + + figure; vs = violinplot(values, labels, 'ShowMean', true); + + ylabel('Partial Directed Coherence'); ylim(yvalues); + title('Granger Causal Influences') + vs(1).ViolinColor = color; + vs(2).ViolinColor = color; + + print(gcf,[load_path('save_fig') name '.png'],'-dpng','-r300'); + close all; + + + %% + seed_12 = outputs.perSubject.IFJa_12; + seed_21 = outputs.perSubject.IFJa_21; + + name = [conn_metric '_' band_name '_' seed_target '_' duration '_' 'R-IFJa and L-IFJp']; + legend = {'R-IFJa to L-IFJp' 'L-IFJp to R-IFJa'}; + labels = repmat(repmat(legend,55,1), 1, 1); + labels = labels(:); + + values = [seed_12(idx_IFJp_L,:)'; seed_21(idx_IFJp_L,:)']; + + figure; vs = violinplot(values, labels, 'ShowMean', true); + + ylabel('Partial Directed Coherence'); ylim(yvalues); + title('Granger Causal Influences') + vs(1).ViolinColor = color; + vs(2).ViolinColor = color; + + print(gcf,[load_path('save_fig') name '.png'],'-dpng','-r300'); + close all; + + + %% 'left-right' +elseif strcmp(seed_target, 'left-right') + + seed_12 = outputs.perSubject.FEF_12; + seed_21 = outputs.perSubject.FEF_21; + + name = [conn_metric '_' band_name '_' seed_target '_' duration '_' 'L-FEF and R-IFJa']; + legend = {'L-FEF to R-IFJa' 'R-IFJa to L-FEF'}; + labels = repmat(repmat(legend,55,1), 1, 1); + labels = labels(:); + + values = [seed_12(idx_IFJa_R,:)'; seed_21(idx_IFJa_R,:)']; + + figure; vs = violinplot(values, labels, 'ShowMean', true); + + ylabel('Partial Directed Coherence'); ylim(yvalues); + title('Granger Causal Influences') + vs(1).ViolinColor = color; + vs(2).ViolinColor = color; + + print(gcf,[load_path('save_fig') name '.png'],'-dpng','-r300'); + close all; + + + %% + seed_12 = outputs.perSubject.FEF_12; + seed_21 = outputs.perSubject.FEF_21; + + name = [conn_metric '_' band_name '_' seed_target '_' duration '_' 'L-FEF and R-IFJp']; + legend = {'L-FEF to R-IFJp' 'R-IFJp to L-FEF'}; + labels = repmat(repmat(legend,55,1), 1, 1); + labels = labels(:); + + values = [seed_12(idx_IFJp_R,:)'; seed_21(idx_IFJp_R,:)']; + + figure; vs = violinplot(values, labels, 'ShowMean', true); + + ylabel('Partial Directed Coherence'); ylim(yvalues); + title('Granger Causal Influences') + vs(1).ViolinColor = color; + vs(2).ViolinColor = color; + + print(gcf,[load_path('save_fig') name '.png'],'-dpng','-r300'); + close all; + + %% + seed_12 = outputs.perSubject.IFJa_12; + seed_21 = outputs.perSubject.IFJa_21; + + name = [conn_metric '_' band_name '_' seed_target '_' duration '_' 'L-IFJa and R-IFJp']; + legend = {'L-IFJa to R-IFJp' 'R-IFJp to L-IFJa'}; + labels = repmat(repmat(legend,55,1), 1, 1); + labels = labels(:); + + values = [seed_12(idx_IFJp_R,:)'; seed_21(idx_IFJp_R,:)']; + + figure; vs = violinplot(values, labels, 'ShowMean', true); + + ylabel('Partial Directed Coherence'); ylim(yvalues); + title('Granger Causal Influences') + vs(1).ViolinColor = color; + vs(2).ViolinColor = color; + + print(gcf,[load_path('save_fig') name '.png'],'-dpng','-r300'); + close all; + + %% same regions L and R + %% + seed_12 = outputs.perSubject.FEF_12; + seed_21 = outputs.perSubject.FEF_21; + + name = [conn_metric '_' band_name '_' seed_target '_' duration '_' 'L-FEF and R-FEF']; + legend = {'L-FEF to R-FEF' 'R-FEF to L-FEF'}; + labels = repmat(repmat(legend,55,1), 1, 1); + labels = labels(:); + + values = [seed_12(idx_FEF_R,:)'; seed_21(idx_FEF_R,:)']; + + figure; vs = violinplot(values, labels, 'ShowMean', true); + + ylabel('Partial Directed Coherence'); ylim(yvalues); + title('Granger Causal Influences') + vs(1).ViolinColor = color; + vs(2).ViolinColor = color; + + print(gcf,[load_path('save_fig') name '.png'],'-dpng','-r300'); + close all; + + + %% + seed_12 = outputs.perSubject.IFJa_12; + seed_21 = outputs.perSubject.IFJa_21; + + name = [conn_metric '_' band_name '_' seed_target '_' duration '_' 'L-IFJa and R-IFJa']; + legend = {'L-IFJa to R-IFJa' 'R-IFJa to L-IFJa'}; + labels = repmat(repmat(legend,55,1), 1, 1); + labels = labels(:); + + values = [seed_12(idx_IFJa_R,:)'; seed_21(idx_IFJa_R,:)']; + + figure; vs = violinplot(values, labels, 'ShowMean', true); + + ylabel('Partial Directed Coherence'); ylim(yvalues); + title('Granger Causal Influences') + vs(1).ViolinColor = color; + vs(2).ViolinColor = color; + + print(gcf,[load_path('save_fig') name '.png'],'-dpng','-r300'); + close all; + + %% + seed_12 = outputs.perSubject.IFJp_12; + seed_21 = outputs.perSubject.IFJp_21; + + name = [conn_metric '_' band_name '_' seed_target '_' duration '_' 'L-IFJp and R-IFJp']; + legend = {'L-IFJp to R-IFJp' 'R-IFJp to L-IFJp'}; + labels = repmat(repmat(legend,55,1), 1, 1); + labels = labels(:); + + values = [seed_12(idx_IFJp_R,:)'; seed_21(idx_IFJp_R,:)']; + + figure; vs = violinplot(values, labels, 'ShowMean', true); + + ylabel('Partial Directed Coherence'); ylim(yvalues); + title('Granger Causal Influences') + vs(1).ViolinColor = color; + vs(2).ViolinColor = color; + + print(gcf,[load_path('save_fig') name '.png'],'-dpng','-r300'); + close all; + + +end + + +end \ No newline at end of file diff --git a/Visualize_Results/scripts/helper_fun/load_path.m b/Visualize_Results/scripts/helper_fun/load_path.m new file mode 100644 index 0000000..7e54ab6 --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/load_path.m @@ -0,0 +1,52 @@ +%% Please set the working directory first before running the scripts. + +function path = load_path(directory) + +workingDir = 'C:\Users\ASUS\Desktop\Thesis Project\'; +fieldTrip_dir = 'C:\Users\ASUS\Desktop\MATLAB\fieldtrip-20210411'; +save_fig = 'C:\Users\ASUS\Desktop\Thesis Project\Main Figures\'; + +% restoredefaultpath % if needed +addpath(fieldTrip_dir); +ft_defaults + +if strcmp(directory, 'inputfile_fs') + path = fullfile(workingDir, 'Visualize results', 'fsaverage'); +elseif strcmp(directory, 'Outputs_Colclough') + path = fullfile(workingDir, 'Visualize results', 'results', 'Outputs_Colclough'); +elseif strcmp(directory, 'Outputs_2s_55subjs_icoh') + path = fullfile(workingDir, 'Visualize results', 'results', 'Outputs_2s_55subjs_icoh'); +elseif strcmp(directory, 'Outputs_2s_55subjs_wPLI') + path = fullfile(workingDir, 'Visualize results', 'results', 'Outputs_2s_55subjs_wPLI'); +elseif strcmp(directory, 'Outputs_2s_55subjs_dwPLI') + path = fullfile(workingDir, 'Visualize results', 'results', 'Outputs_2s_55subjs_dwPLI'); +elseif strcmp(directory, 'Outputs_2s_55subjs_powcorrOrtho') + path = fullfile(workingDir, 'Visualize results', 'results', 'Outputs_2s_55subjs_powcorrOrtho'); +elseif strcmp(directory, 'Outputs_2s_55subjs_granger') + path = fullfile(workingDir, 'Visualize results', 'results', 'Outputs_2s_55subjs_granger'); +elseif strcmp(directory, 'Outputs_2s_55subjs_pdc') + path = fullfile(workingDir, 'Visualize results', 'results', 'Outputs_2s_55subjs_pdc'); + +elseif strcmp(directory, 'Outputs_5s_55subjs_icoh') + path = fullfile(workingDir, 'Visualize results', 'results', 'Outputs_5s_55subjs_icoh'); +elseif strcmp(directory, 'Outputs_5s_55subjs_wPLI') + path = fullfile(workingDir, 'Visualize results', 'results', 'Outputs_5s_55subjs_wPLI'); +elseif strcmp(directory, 'Outputs_5s_55subjs_powcorrOrtho') + path = fullfile(workingDir, 'Visualize results', 'results', 'Outputs_5s_55subjs_powcorrOrtho'); + +elseif strcmp(directory, 'Outputs_10s_55subjs_icoh') + path = fullfile(workingDir, 'Visualize results', 'results', 'Outputs_10s_55subjs_icoh'); +elseif strcmp(directory, 'Outputs_10s_55subjs_powcorrOrtho') + path = fullfile(workingDir, 'Visualize results', 'results', 'Outputs_10s_55subjs_powcorrOrtho'); + +elseif strcmp(directory, 'workingDir') + path = workingDir; +elseif strcmp(directory, 'helper_fun') + path = fullfile(workingDir, 'Visualize results', 'scripts', 'helper_fun'); +elseif strcmp(directory, 'save_fig') + path = save_fig; + + +end + + \ No newline at end of file diff --git a/Visualize_Results/scripts/helper_fun/tutorial_nwa_connectivityviewer.m b/Visualize_Results/scripts/helper_fun/tutorial_nwa_connectivityviewer.m new file mode 100644 index 0000000..0466dfd --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/tutorial_nwa_connectivityviewer.m @@ -0,0 +1,325 @@ +function tutorial_nwa_connectivityviewer(data, param, colorlim, seed_target, flag) + +if nargin<2, + param = 'cohspctrm'; +end +if nargin<3, + tmp = data.(param); + if all(tmp(:)>=0), + colorlim = [0 max(tmp(:))]; + else + colorlim = [-1 1].*max(abs(tmp(:))); + end +end + +if isfield(data, 'brainordinate') && ~isfield(data, 'pos') + % this seems to be parcellated data, 'unparcellate' on the fly + tmpdat = data.(param); + data = data.brainordinate; + atlas = data; + + dat = zeros(size(data.pos,1)); + for k = 1:numel(data.parcellationlabel) + for m = 1:numel(data.parcellationlabel) + dat(data.parcellation==k,data.parcellation==m) = tmpdat(k,m); + end + end +else + dat = data.(param); + atlas = []; +end + +% create sphere for location highlighting +[x,y,z] = sphere(10); + +figure; hold on; +ft_plot_mesh(data,'vertexcolor', mean(dat,2)) +h1 = gca; +h2 = get(h1, 'children'); +if flag == false + h2.EdgeColor = [0.65,0.65,0.65]; %[0.80,0.80,0.80];%[0.94,0.94,0.94]; + h2.LineStyle = ':'; +end +h3 = surf(h1, 2*x,2*y,2*z+100, 'facecolor', 'w', 'edgecolor', 'none'); + +% % Orhan: +% change the view angle +if seed_target(1) == 'r' + view(360,0) +elseif seed_target(1) == 'l' + view(180,0) +end + +set(gcf,'windowbuttondownfcn',{@fancyplotting, [h1 h2 h3], data.pos, data.tri, dat, atlas, seed_target, flag}) +% Orhan: +cbar = colorbar; +cbar.FontSize = 14; +set(cbar.XLabel,{'String','Rotation','Position','FontWeight'},{'Z',0,[0.5 colorlim(1)-0.1],'bold'}) + +caxis(colorlim); + +function fancyplotting(hObject, ~, axh, pos, tri, dat, atlas, seed_target, flag) + +% get the clicked position in axis coordinates +pos2d = get(axh(1),'CurrentPoint'); + +% Orhan: Delete later +% nn = 0; +% if mod(nn,2) == 0 +% pos2d = [ 53.9546 -69.3638 74.6208 ; 53.9546 69.4183 74.6208]; +% nn = nn+1; +% else +% pos2d = [51.6434 -52.0650 132.1000 ; 51.6434 -52.0650 -6.6821]; +% nn = nn+1; +% end + +% get the intersection with the mesh +[ipos, d] = intersect_line(pos, tri, pos2d(1,:), pos2d(2,:)); +[md, ix] = min(abs(d)); + +% get the closest mesh point, and the corresponding index +dpos = pos - ipos(ix*ones(size(pos,1),1),:); +indx = nearest(sum(dpos.^2,2),0); + +x = get(axh(3),'xdata'); +x = x-mean(x(:)); +y = get(axh(3),'ydata'); +y = y-mean(y(:)); +z = get(axh(3),'zdata'); +z = z-mean(z(:)); + +set(axh(3),'xdata', x+pos(indx,1)); +set(axh(3),'ydata', y+pos(indx,2)); +set(axh(3),'zdata', z+pos(indx,3)); + +% update the functional data +set(axh(2), 'facevertexcdata', dat(:, indx)); + +% Orhan: Not used anymore, but keep it for later +% for i = 1:numel(atlas.parcellationlabel) +% scout_idx = (atlas.parcellation == i); +% P = atlas.pos(scout_idx,:); +% k = boundary(P); +% hold on; +% trisurf(k,P(:,1),P(:,2),P(:,3),'Facecolor','black','FaceAlpha', 0.1,'EdgeColor', 'none') +% end + +% Orhan: Only for ROIs' contours; you need add 'ROI_idx' parameter to func +% for i = 1:length(ROI_idx) +% scout_idx = (atlas.parcellation == ROI_idx(i)); +% ft_plot_mesh(atlas, 'vertexcolor', dat(:, indx), 'contour', scout_idx, 'contourcolor', 'k', 'contourlinewidth', 0.5); +% end + +% % Orhan: for contour of the parcels +if flag == true + if or(strcmp(seed_target, 'right-right'), strcmp(seed_target, 'left-right')) + for i = 1:180 + mask = atlas.parcellation == i+180; + scout_idx{i,:} = double(mask); + end + ft_plot_mesh(atlas, 'vertexcolor', dat(:, indx), 'contour', scout_idx, 'contourcolor', [185, 185, 185]/255, 'contourlinewidth', 0.5); + else + for i = 1:180 + mask = atlas.parcellation == i; + scout_idx{i,:} = double(mask); + end + ft_plot_mesh(atlas, 'vertexcolor', dat(:, indx), 'contour', scout_idx, 'contourcolor', [185, 185, 185]/255, 'contourlinewidth', 0.5); + end +end + +if isempty(atlas) + fprintf('you clicked on position %d with coordinates %3.2f %3.2f %3.2f\n', [indx pos(indx,:)]); +else + fprintf('you clicked on parcel %s\n', atlas.parcellationlabel{atlas.parcellation(indx)}); + title(atlas.parcellationlabel{atlas.parcellation(indx)}, 'Interpreter', 'none'); +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% the code below is copied from fieldtrip/private + +function [points, pos, indx] = intersect_line(pnt, tri, pnt1, pnt2) + +% INTERSECT_LINE finds the intersection points between a mesh and a line. +% +% Use as: +% [points, pos, indx] = intersect_line(pnt, tri, pnt1, pnt2) +% +% Where pnt (Nx3) and tri (Mx3) define the mesh, and pnt1 (1x3) and pnt2 +% (1x3) define the line. The output argument points (Px3) are the +% intersection points, pos (Px1) the location on the line (relative to +% pnt1) and indx is the index to the triangles of the mesh that are +% intersected. +% +% This code is based from a function from the geom3d toolbox, that can be +% found on matlab's file exchange. The original help is pasted below. The +% original function was released under the BSD-license. +% +% Adapted to FieldTrip by Jan-Mathijs Schoffelen 2012 + +%INTERSECTLINEMESH3D Intersection points of a 3D line with a mesh +% +% INTERS = intersectLineMesh3d(LINE, VERTICES, FACES) +% Compute the intersection points between a 3D line and a 3D mesh defined +% by vertices and faces. +% +% [INTERS POS INDS] = intersectLineMesh3d(LINE, VERTICES, FACES) +% Also returns the position of each intersection point on the input line, +% and the index of the intersected faces. +% If POS > 0, the point is also on the ray corresponding to the line. +% +% Example +% intersectLineMesh3d +% +% See also +% meshes3d, triangulateFaces +% +% ------ +% Author: David Legland +% e-mail: david.legland@grignon.inra.fr +% Created: 2011-12-20, using Matlab 7.9.0.529 (R2009b) +% Copyright 2011 INRA - Cepia Software Platform. + + +tol = 1e-12; + +ntri = size(tri,1); + +% normals to the triangles +n = normals(pnt, tri, 'triangle'); + +% vectors describing the edges +t0 = pnt(tri(:,1),:); +u = pnt(tri(:,2),:) - t0; +v = pnt(tri(:,3),:) - t0; + +% get the direction of the line +d = pnt2 - pnt1; +d = d/norm(d); + +% vector between triangle origin and line origin +w0 = pnt1(ones(ntri,1),:) - t0; + +a = -sum(n.*w0, 2); +b = n*d'; + +valid = abs(b)>tol & sqrt(sum(n.^2,2))>tol; + +% compute intersection point of line with supporting plane +% If pos < 0: point before ray +% IF pos > |dir|: point after edge +pos = a ./ b; + +% coordinates of intersection point +points = pnt1(ones(ntri,1),:) + pos*d; + +% normalize direction vectors of triangle edges +uu = sum(u.^2, 2); +uv = sum(u.*v, 2); +vv = sum(v.^2, 2); + +% coordinates of vector v in triangle basis +w = points - t0; +wu = sum(w.*u, 2); +wv = sum(w.*v, 2); + +% normalization constant +D = uv.^2 - uu .* vv; + +% test first coordinate +s = (uv .* wv - vv .* wu) ./ D; +ind1 = s < 0.0 | s > 1.0; +points(ind1, :) = NaN; +pos(ind1) = NaN; + +% test second coordinate, and third triangle edge +t = (uv .* wu - uu .* wv) ./ D; +ind2 = t < 0.0 | (s + t) > 1.0; +points(ind2, :) = NaN; +pos(ind2) = NaN; + +% keep only interesting points +inds = ~ind1 & ~ind2 & valid; +points = points(inds, :); + +pos = pos(inds); +indx = find(inds); + +function [nrm] = normals(pnt, tri, opt) + +% NORMALS compute the surface normals of a triangular mesh +% for each triangle or for each vertex +% +% [nrm] = normals(pnt, tri, opt) +% where opt is either 'vertex' or 'triangle' + +% Copyright (C) 2002-2007, Robert Oostenveld +% +% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org +% for the documentation and details. +% +% FieldTrip is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% FieldTrip is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with FieldTrip. If not, see . +% +% $Id$ + +if nargin<3 + opt='vertex'; +elseif (opt(1)=='v' | opt(1)=='V') + opt='vertex'; +elseif (opt(1)=='t' | opt(1)=='T') + opt='triangle'; +else + error('invalid optional argument'); +end + +npnt = size(pnt,1); +ntri = size(tri,1); + +% shift to center +pnt(:,1) = pnt(:,1)-mean(pnt(:,1),1); +pnt(:,2) = pnt(:,2)-mean(pnt(:,2),1); +pnt(:,3) = pnt(:,3)-mean(pnt(:,3),1); + +% compute triangle normals +% nrm_tri = zeros(ntri, 3); +% for i=1:ntri +% v2 = pnt(tri(i,2),:) - pnt(tri(i,1),:); +% v3 = pnt(tri(i,3),:) - pnt(tri(i,1),:); +% nrm_tri(i,:) = cross(v2, v3); +% end + +% vectorized version of the previous part +v2 = pnt(tri(:,2),:) - pnt(tri(:,1),:); +v3 = pnt(tri(:,3),:) - pnt(tri(:,1),:); +nrm_tri = cross(v2, v3); + + +if strcmp(opt, 'vertex') + % compute vertex normals + nrm_pnt = zeros(npnt, 3); + for i=1:ntri + nrm_pnt(tri(i,1),:) = nrm_pnt(tri(i,1),:) + nrm_tri(i,:); + nrm_pnt(tri(i,2),:) = nrm_pnt(tri(i,2),:) + nrm_tri(i,:); + nrm_pnt(tri(i,3),:) = nrm_pnt(tri(i,3),:) + nrm_tri(i,:); + end + % normalise the direction vectors to have length one + nrm = nrm_pnt ./ (sqrt(sum(nrm_pnt.^2, 2)) * ones(1,3)); +else + % normalise the direction vectors to have length one + nrm = nrm_tri ./ (sqrt(sum(nrm_tri.^2, 2)) * ones(1,3)); +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% fast cross product to replace the MATLAB standard version +function [c] = cross(a,b) +c = [a(:,2).*b(:,3)-a(:,3).*b(:,2) a(:,3).*b(:,1)-a(:,1).*b(:,3) a(:,1).*b(:,2)-a(:,2).*b(:,1)]; diff --git a/Visualize_Results/scripts/helper_fun/tutorial_nwa_connectivityviewer_save_figures.m b/Visualize_Results/scripts/helper_fun/tutorial_nwa_connectivityviewer_save_figures.m new file mode 100644 index 0000000..24ee6e9 --- /dev/null +++ b/Visualize_Results/scripts/helper_fun/tutorial_nwa_connectivityviewer_save_figures.m @@ -0,0 +1,328 @@ +function tutorial_nwa_connectivityviewer_save_figures(data, param, colorlim, seed_target, flag, pos2d) + +if nargin<2, + param = 'cohspctrm'; +end +if nargin<3, + tmp = data.(param); + if all(tmp(:)>=0), + colorlim = [0 max(tmp(:))]; + else + colorlim = [-1 1].*max(abs(tmp(:))); + end +end + +if isfield(data, 'brainordinate') && ~isfield(data, 'pos') + % this seems to be parcellated data, 'unparcellate' on the fly + tmpdat = data.(param); + data = data.brainordinate; + atlas = data; + + dat = zeros(size(data.pos,1)); + for k = 1:numel(data.parcellationlabel) + for m = 1:numel(data.parcellationlabel) + dat(data.parcellation==k,data.parcellation==m) = tmpdat(k,m); + end + end +else + dat = data.(param); + atlas = []; +end + +% create sphere for location highlighting +[x,y,z] = sphere(10); + +figure; hold on; +ft_plot_mesh(data,'vertexcolor', mean(dat,2)) +h1 = gca; +h2 = get(h1, 'children'); +if flag == false + h2.EdgeColor = [0.65,0.65,0.65]; %[0.80,0.80,0.80];%[0.94,0.94,0.94]; + h2.LineStyle = ':'; +end +h3 = surf(h1, 2*x,2*y,2*z+100, 'facecolor', 'w', 'edgecolor', 'none'); + +% % Orhan: +% change the view angle +if seed_target(1) == 'r' + view(360,0) +elseif seed_target(1) == 'l' + view(180,0) +end + +% Orhan: +cbar = colorbar; +cbar.FontSize = 14; +set(cbar.XLabel,{'String','Rotation','Position','FontWeight'},{'Z',0,[0.5 colorlim(1)-0.1],'bold'}) + +% % Orhan: +fancyplotting([h1 h2 h3], data.pos, data.tri, dat, atlas, seed_target, flag, pos2d); + +caxis(colorlim); +end + +function fancyplotting(axh, pos, tri, dat, atlas, seed_target, flag, pos2d) + +% get the clicked position in axis coordinates +% pos2d = get(axh(1),'CurrentPoint'); + +% get the intersection with the mesh +[ipos, d] = intersect_line(pos, tri, pos2d(1,:), pos2d(2,:)); +[md, ix] = min(abs(d)); + +% get the closest mesh point, and the corresponding index +dpos = pos - ipos(ix*ones(size(pos,1),1),:); +indx = nearest(sum(dpos.^2,2),0); + +x = get(axh(3),'xdata'); +x = x-mean(x(:)); +y = get(axh(3),'ydata'); +y = y-mean(y(:)); +z = get(axh(3),'zdata'); +z = z-mean(z(:)); + +set(axh(3),'xdata', x+pos(indx,1)); +set(axh(3),'ydata', y+pos(indx,2)); +set(axh(3),'zdata', z+pos(indx,3)); + +% update the functional data +set(axh(2), 'facevertexcdata', dat(:, indx)); + +% Orhan: Not used anymore, but keep it for later +% for i = 1:numel(atlas.parcellationlabel) +% scout_idx = (atlas.parcellation == i); +% P = atlas.pos(scout_idx,:); +% k = boundary(P); +% hold on; +% trisurf(k,P(:,1),P(:,2),P(:,3),'Facecolor','black','FaceAlpha', 0.1,'EdgeColor', 'none') +% end + +% Orhan: Only for ROIs' contours; you need add 'ROI_idx' parameter to func +% for i = 1:length(ROI_idx) +% scout_idx = (atlas.parcellation == ROI_idx(i)); +% ft_plot_mesh(atlas, 'vertexcolor', dat(:, indx), 'contour', scout_idx, 'contourcolor', 'k', 'contourlinewidth', 0.5); +% end + +% % Orhan: for contour of the parcels +if flag == true + if or(strcmp(seed_target, 'right-right'), strcmp(seed_target, 'left-right')) + for i = 1:180 + mask = atlas.parcellation == i+180; + scout_idx{i,:} = double(mask); + end + ft_plot_mesh(atlas, 'vertexcolor', dat(:, indx), 'contour', scout_idx, 'contourcolor', [185, 185, 185]/255, 'contourlinewidth', 0.5); + else + for i = 1:180 + mask = atlas.parcellation == i; + scout_idx{i,:} = double(mask); + end + ft_plot_mesh(atlas, 'vertexcolor', dat(:, indx), 'contour', scout_idx, 'contourcolor', [185, 185, 185]/255, 'contourlinewidth', 0.5); + end +end + +if isempty(atlas) + fprintf('you clicked on position %d with coordinates %3.2f %3.2f %3.2f\n', [indx pos(indx,:)]); +else + fprintf('you clicked on parcel %s\n', atlas.parcellationlabel{atlas.parcellation(indx)}); + title(atlas.parcellationlabel{atlas.parcellation(indx)}, 'Interpreter', 'none'); +end + + + + + +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% the code below is copied from fieldtrip/private + +function [points, pos, indx] = intersect_line(pnt, tri, pnt1, pnt2) + +% INTERSECT_LINE finds the intersection points between a mesh and a line. +% +% Use as: +% [points, pos, indx] = intersect_line(pnt, tri, pnt1, pnt2) +% +% Where pnt (Nx3) and tri (Mx3) define the mesh, and pnt1 (1x3) and pnt2 +% (1x3) define the line. The output argument points (Px3) are the +% intersection points, pos (Px1) the location on the line (relative to +% pnt1) and indx is the index to the triangles of the mesh that are +% intersected. +% +% This code is based from a function from the geom3d toolbox, that can be +% found on matlab's file exchange. The original help is pasted below. The +% original function was released under the BSD-license. +% +% Adapted to FieldTrip by Jan-Mathijs Schoffelen 2012 + +%INTERSECTLINEMESH3D Intersection points of a 3D line with a mesh +% +% INTERS = intersectLineMesh3d(LINE, VERTICES, FACES) +% Compute the intersection points between a 3D line and a 3D mesh defined +% by vertices and faces. +% +% [INTERS POS INDS] = intersectLineMesh3d(LINE, VERTICES, FACES) +% Also returns the position of each intersection point on the input line, +% and the index of the intersected faces. +% If POS > 0, the point is also on the ray corresponding to the line. +% +% Example +% intersectLineMesh3d +% +% See also +% meshes3d, triangulateFaces +% +% ------ +% Author: David Legland +% e-mail: david.legland@grignon.inra.fr +% Created: 2011-12-20, using Matlab 7.9.0.529 (R2009b) +% Copyright 2011 INRA - Cepia Software Platform. + + +tol = 1e-12; + +ntri = size(tri,1); + +% normals to the triangles +n = normals(pnt, tri, 'triangle'); + +% vectors describing the edges +t0 = pnt(tri(:,1),:); +u = pnt(tri(:,2),:) - t0; +v = pnt(tri(:,3),:) - t0; + +% get the direction of the line +d = pnt2 - pnt1; +d = d/norm(d); + +% vector between triangle origin and line origin +w0 = pnt1(ones(ntri,1),:) - t0; + +a = -sum(n.*w0, 2); +b = n*d'; + +valid = abs(b)>tol & sqrt(sum(n.^2,2))>tol; + +% compute intersection point of line with supporting plane +% If pos < 0: point before ray +% IF pos > |dir|: point after edge +pos = a ./ b; + +% coordinates of intersection point +points = pnt1(ones(ntri,1),:) + pos*d; + +% normalize direction vectors of triangle edges +uu = sum(u.^2, 2); +uv = sum(u.*v, 2); +vv = sum(v.^2, 2); + +% coordinates of vector v in triangle basis +w = points - t0; +wu = sum(w.*u, 2); +wv = sum(w.*v, 2); + +% normalization constant +D = uv.^2 - uu .* vv; + +% test first coordinate +s = (uv .* wv - vv .* wu) ./ D; +ind1 = s < 0.0 | s > 1.0; +points(ind1, :) = NaN; +pos(ind1) = NaN; + +% test second coordinate, and third triangle edge +t = (uv .* wu - uu .* wv) ./ D; +ind2 = t < 0.0 | (s + t) > 1.0; +points(ind2, :) = NaN; +pos(ind2) = NaN; + +% keep only interesting points +inds = ~ind1 & ~ind2 & valid; +points = points(inds, :); + +pos = pos(inds); +indx = find(inds); + +end + +function [nrm] = normals(pnt, tri, opt) + +% NORMALS compute the surface normals of a triangular mesh +% for each triangle or for each vertex +% +% [nrm] = normals(pnt, tri, opt) +% where opt is either 'vertex' or 'triangle' + +% Copyright (C) 2002-2007, Robert Oostenveld +% +% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org +% for the documentation and details. +% +% FieldTrip is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% FieldTrip is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with FieldTrip. If not, see . +% +% $Id$ + +if nargin<3 + opt='vertex'; +elseif (opt(1)=='v' | opt(1)=='V') + opt='vertex'; +elseif (opt(1)=='t' | opt(1)=='T') + opt='triangle'; +else + error('invalid optional argument'); +end + +npnt = size(pnt,1); +ntri = size(tri,1); + +% shift to center +pnt(:,1) = pnt(:,1)-mean(pnt(:,1),1); +pnt(:,2) = pnt(:,2)-mean(pnt(:,2),1); +pnt(:,3) = pnt(:,3)-mean(pnt(:,3),1); + +% compute triangle normals +% nrm_tri = zeros(ntri, 3); +% for i=1:ntri +% v2 = pnt(tri(i,2),:) - pnt(tri(i,1),:); +% v3 = pnt(tri(i,3),:) - pnt(tri(i,1),:); +% nrm_tri(i,:) = cross(v2, v3); +% end + +% vectorized version of the previous part +v2 = pnt(tri(:,2),:) - pnt(tri(:,1),:); +v3 = pnt(tri(:,3),:) - pnt(tri(:,1),:); +nrm_tri = cross(v2, v3); + + +if strcmp(opt, 'vertex') + % compute vertex normals + nrm_pnt = zeros(npnt, 3); + for i=1:ntri + nrm_pnt(tri(i,1),:) = nrm_pnt(tri(i,1),:) + nrm_tri(i,:); + nrm_pnt(tri(i,2),:) = nrm_pnt(tri(i,2),:) + nrm_tri(i,:); + nrm_pnt(tri(i,3),:) = nrm_pnt(tri(i,3),:) + nrm_tri(i,:); + end + % normalise the direction vectors to have length one + nrm = nrm_pnt ./ (sqrt(sum(nrm_pnt.^2, 2)) * ones(1,3)); +else + % normalise the direction vectors to have length one + nrm = nrm_tri ./ (sqrt(sum(nrm_tri.^2, 2)) * ones(1,3)); +end +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% fast cross product to replace the MATLAB standard version +function [c] = cross(a,b) +c = [a(:,2).*b(:,3)-a(:,3).*b(:,2) a(:,3).*b(:,1)-a(:,1).*b(:,3) a(:,1).*b(:,2)-a(:,2).*b(:,1)]; +end diff --git a/Visualize_Results/scripts/helper_saveFig_effectiveConn.m b/Visualize_Results/scripts/helper_saveFig_effectiveConn.m new file mode 100644 index 0000000..6f7758e --- /dev/null +++ b/Visualize_Results/scripts/helper_saveFig_effectiveConn.m @@ -0,0 +1,81 @@ +band_name = ["delta" "theta" "alpha" "beta" "gamma"]; +seed_target = [ "right-left" "left-right"]; %"left-left" "right-right" + +for ii = 1:5 + for kk = 1:2 + helper_automaticSave_effective(char(band_name(ii)), char(seed_target(kk))) + end +end + + +function helper_automaticSave_effective(band_name, seed_target) + +addpath(genpath(fullfile(pwd, 'helper_fun'))); + +%% MEANING of colors +% y: target region +% +% GREEN : seed --> y direction is significant. A seed region explains the varience in a target better. +% ~ sending signal +% MAGENTA : y--> seed direction is significant. A target region explains the varience in a seed better. +% ~ receiving signal + +%% Select +statistics.analysis_type = 'ROI'; % 'ROI' 'exploratory' +% 'ROI' : ROI analysis per hemisphere. FDR correction for 33+3 seed regions. +% 'exploratory' : exploratory analysis per hemisphere. FDR correction for 180 regions. + +duration = '2s'; %fixed % only '2s' +conn_metric = 'pdc'; % 'pdc' +% band_name = 'beta'; % 'delta' 'theta' 'alpha' 'beta' 'gamma' +% seed_target = 'right-right'; % 'right-right' 'right-left' 'left-left' 'left-right' + +%% +% figures +fsaverage_fig = true; +violinplot = false; + +% Paired Wilcoxon paired test against zero +statistics.alpha = 0.001; % 0.05 0.01 0.001 +% FDR correction +statistics.corrected = true; % true false +statistics.correction = 'FDR'; + +% to show the contours of parcels; +%%keep it false for a faster processing +flag.contours = false; + +% save the figures automatic? +save_fig = true; +dir_fig = load_path('save_fig'); + +%% +colorlim_value = [-4 4]; +%%% fixed +circular_graph = false; +font_size = 15; +flag.fsaverage = fsaverage_fig; % to show the figure. +flag.circular_graph = circular_graph; +flag.font_size = font_size; +flag.save_fig = save_fig; +flag.duration = duration; +flag.dir_fig = dir_fig; +inputfile_conn = fullfile(load_path(['Outputs_' duration '_55subjs_' conn_metric]), ['55Subjects_groupResults_' duration '.mat']); + +%% Run the function +cimec.inputfile_conn = inputfile_conn; +cimec.band_name = band_name; +cimec.seed_target = seed_target; +cimec.conn_metric = conn_metric; +cimec.statistics = statistics; +cimec.colorlim_value = colorlim_value; +cimec.flag = flag; + +outputs = helper_effectiveConn_fsaverage(cimec); + +%% violinplot +if violinplot == true + helper_violinplot(outputs, conn_metric, band_name, seed_target, duration) +end + +end \ No newline at end of file diff --git a/Visualize_Results/scripts/helper_saveFig_functionalConn.m b/Visualize_Results/scripts/helper_saveFig_functionalConn.m new file mode 100644 index 0000000..a7d78b4 --- /dev/null +++ b/Visualize_Results/scripts/helper_saveFig_functionalConn.m @@ -0,0 +1,73 @@ +band_name = ["delta" "theta" "alpha" "beta" "gamma"]; +seed_target = ["left-left" "right-right"]; + +for ii = 1:5 + for kk = 1:2 + helper_automaticSave_functional(char(band_name(ii)), char(seed_target(kk))) + end +end + + +function helper_automaticSave_functional(band_name, seed_target) + +addpath(genpath(fullfile(pwd, 'helper_fun'))); +% colormapeditor +set(0, 'DefaultFigureColormap', cmap_rbw); + +%% Select +statistics.analysis_type = 'ROI'; % 'ROI' 'exploratory' +% 'ROI' : ROI analysis per hemisphere. FDR correction for 33 regions. +% 'exploratory' : exploratory analysis per hemisphere. FDR correction for 180 regions. + +duration = '2s'; % '2s' '5s' '10s' +conn_metric = 'powcorrOrtho'; % 'powcorrOrtho' 'icoh' 'dwPLI' +% band_name = 'delta'; % 'delta' 'theta' 'alpha' 'beta' 'gamma' 'collapseFreq' +% seed_target = 'left-left'; % 'right-right' 'right-left' 'left-left' 'left-right' + +% contrast +contrast = true; + +% figures +fsaverage_fig = true; +circular_graph = false; +font_size = 15; + +% Wilcoxon signed rank test +statistics.alpha = 0.05; % 0.05 0.01 0.001 +% FDR correction +statistics.corrected = true; % true false +statistics.correction = 'FDR'; + +% to show the contours of parcels; +%%keep it false for a faster processing +flag.contours = false; + +% save the figures automatic? +save_fig = true; +dir_fig = load_path('save_fig'); + +%% fixed +flag.fsaverage = fsaverage_fig; % to show the figure. +flag.circular_graph = circular_graph; +flag.font_size = font_size; +flag.save_fig = save_fig; +flag.dir_fig = dir_fig; +flag.duration = duration; +inputfile_conn = fullfile(load_path(['Outputs_' duration '_55subjs_' conn_metric]), ['55Subjects_groupResults_' duration '.mat']); + +%% Run the function +cimec.inputfile_conn = inputfile_conn; +cimec.band_name = band_name; +cimec.seed_target = seed_target; +cimec.conn_metric = conn_metric; +cimec.statistics = statistics; +cimec.contrast = contrast; +cimec.flag = flag; + +if strcmp(band_name, 'collapseFreq') + outputs = helper_functionConn_fsaverage_freqCollapsed(cimec); +else + outputs = helper_functionConn_fsaverage(cimec); +end + +end \ No newline at end of file