Update Create_GCMs.m, MRS_QC.m, and 4 more files...

Codes/Create_GCMs.m

+%% CODE TO CREATE GCM - GIT
+
+% This code needs the usual your_paths.txt file.
+
+% Once the NMM model is run, this code
+% will create a GCM file to input in the PEB analysis. 
+
+%-----------------------------------------------------------------------
+% Written by I. Paparella
+% Cyclotron Research Centre, University of Liege, Belgium
+% 14th December 2023
+
+
+function Create_GCMs 
+
+%% Path info
+paths_txt = fileread('your_paths.txt') ;
+paths_txt = char(strsplit(paths_txt, '\n')) ;
+
+pathtodata = (strtrim(paths_txt(2,:)));
+
+%% Pre define variables
+GCM = [];
+
+%%  Select the subject to analyse
+cd(pathtodata);
+[subdir] = spm_select(Inf,'dir','Select subject_ID directory');
+IDs = str2num(subdir(:,end-2:end));
+
+%% Loop across subjects
+for isub = 1:size(subdir,1)
+    subject_ID = sprintf('sub-%03d', IDs(isub));
+    subject_Path = strcat (pathtodata, '/', subject_ID, '/TMS_EEG');
+    
+    % Go in the session folder and load DCM
+    DCMfile_sess=spm_select('List',subject_Path,'DCM_DefaultCodes_.*.mat$');
+    if isempty(DCMfile_sess)
+        warning('No DCM for subject %s, session %s was found. Please check the subjects you selected', subject_ID, session)
+        return
+    end
+    
+    Loaded_DCMfile_sess = load(strcat(subject_Path, '/', DCMfile_sess));
+    
+    % Fill the GCM according to the session with all subjects together
+    GCM = [GCM; {Loaded_DCMfile_sess.DCM}];
+
+end
+
+% Make PEB directory
+mkdir(strcat(pathtodata,'/PEB'));
+
+% Save GCM with the correct naming
+nameTosav= 'allFemale';
+
+GCM = GCM;
+save(strcat(pathtodata,'/PEB/', sprintf('GCM_%s_%s', session, nameTosav)), 'GCM')
+
+fprintf('GCM has been created. \n')
+end
\ No newline at end of file

Codes/MRS_QC.m

+%% CODE TO PERFORM ADDITIONAL QC MRS ANALYSIS REQUIRED IN 1ST ROUND REV
+
+% This code extracts QC MRS metric (linewidth, CRLB and SNR across
+% participants)
+
+% It needs;
+% a "your_paths.txt" file with all the info on your paths. An example is
+% provided on GIT + the fitted MRS data.
+%-----------------------------------------------------------------------
+% Written by I. Paparella
+% Cyclotron Research Centre, University of Liege, Belgium
+% 06th November 2024
+
+%% Path info
+paths_txt = fileread('your_paths.txt') ;
+paths_txt = char(strsplit(paths_txt, '\n')) ;
+
+pathtodata = (strtrim(paths_txt(2,:)));
+
+%%  Select the subject to analyse (all except 001, 007 and 011 which are males)
+cd(pathtodata);
+[subdir] = spm_select(Inf,'dir','Select subject_ID directory');
+IDs = str2num(subdir(:,end-2:end));
+
+%% Variables needed
+%Metrics= ['SNR', 'FWHM', '%CRLB'];
+Metabs ={'GABA', 'NAA', 'Cr'};
+%% Initalize variables
+%MRS_QC_metrics = zeros(size(subdir,1),(length(Metrics)*length(Metabs)));
+MRS_QC_metrics = [];
+%% Loop across subjects to analyse and get subjective parameters
+for isub = 1:size(subdir,1)
+    
+    % Extract path info per subject
+    subject_ID = sprintf('sub-%03d', IDs(isub));
+    subject_MRSfit = fullfile (pathtodata, subject_ID, 'MRS');
+    
+    % Load .csv summary table
+    SubMRS= readtable(fullfile(subject_MRSfit, '/summary.csv'));
+    SubMRS.Metab = categorical(SubMRS.Metab); %convert to categorical
+    
+    % Extract needed info
+    MRS_QC_metrics_Sbj = [];
+    for metabi= 1: length(Metabs)
+        MRS_QC_metrics_Sbj= [MRS_QC_metrics_Sbj table2array(SubMRS(SubMRS.Metab== Metabs(metabi), end-2:end))];
+    end
+    MRS_QC_metrics = [MRS_QC_metrics; MRS_QC_metrics_Sbj];
+end
+
+MRS_QC_metrics=array2table(MRS_QC_metrics);
+IDs_table = table(IDs);
+
+MRS_QC_metrics_full = [IDs_table MRS_QC_metrics];
+
+MRS_QC_metrics_full.Properties.VariableNames = {'Sbj_number'; 'GABA_CRLB'; 'GABA_SNR';'GABA_FWHM';...
+    'NAA_CRLB'; 'NAA_SNR';'NAA_FWHM';'Cr_CRLB'; 'Cr_SNR';'Cr_FWHM'};
+writetable(MRS_QC_metrics_full,fullfile(pathtodata,'MRS_QC_1stRev_AllFemale.csv'),'WriteRowNames',true);
\ No newline at end of file

Codes/NMM_QCfit.m

+%% CODE TO PERFORM ADDITIONAL QC MRS ANALYSIS REQUIRED IN 1ST ROUND REV
+
+% This code extracts QC NMM fit metrics (R2 and correlations between
+% observed and modeled data)
+
+% It needs
+% a "your_paths.txt" file with all the info on your paths. An example is
+% provided on GIT + the DCM estimated.
+%-----------------------------------------------------------------------
+% Written by I. Paparella
+% Cyclotron Research Centre, University of Liege, Belgium
+% 07th November 2024
+
+%% Path info
+paths_txt = fileread('your_paths.txt') ;
+paths_txt = char(strsplit(paths_txt, '\n')) ;
+
+pathtodata = (strtrim(paths_txt(2,:)));
+
+%%  Select the subject to analyse (all except 001, 007 and 011 which are males)
+cd(pathtodata);
+[subdir] = spm_select(Inf,'dir','Select subject_ID directory');
+IDs = str2num(subdir(:,end-2:end));
+
+%% Initalize variables
+%MRS_QC_metrics = zeros(size(subdir,1),(length(Metrics)*length(Metabs)));
+R2 = [];
+r = [];
+modes= [1 2 3];
+%% Loop across subjects to analyse and get subjective parameters
+for isub = 1:size(subdir,1)
+    
+    % Extract path info per subject
+    subject_ID = sprintf('sub-%03d', IDs(isub));
+    subject_DCMest = fullfile (pathtodata, subject_ID, 'TMS_EEG');
+    
+    % Load DCM file
+    DCMfile = spm_select('List',subject_DCMest,'DCM_DefaultCodes_.*.mat$');
+    SubDCM=load(fullfile(subject_DCMest, '/', DCMfile(1,:)));
+    
+% Compute R2    
+    for m = 1:length(modes)
+        % compute and store R2 for ERP
+        MSE_1 = sum(SubDCM.DCM.R{1}(:,m).^2); % for condition 1
+        SS_1  = sum((SubDCM.DCM.H{1}(:,m)+SubDCM.DCM.R{1}(:,m)).^2);
+        R2(isub,m,1) = 1-(MSE_1/SS_1);
+        clear MSE_1 SS_1
+    end
+    
+% Compute correlations between observed and modeled data    
+    for m = 1:length(modes)
+        % correlation
+        y1 = SubDCM.DCM.H{1}(:,m)+SubDCM.DCM.R{1}(:,m);        
+        yhat1 = SubDCM.DCM.H{1}(:,m);
+        r(isub,m,1) = corr(y1(:,1),yhat1(:,1),'type','Pearson');
+    end    
+end
+
+% Convert to table and save
+r=array2table(r);
+R2=array2table(R2);
+IDs_table = table(IDs);
+
+QC_metrics_full = [IDs_table r R2];
+
+QC_metrics_full.Properties.VariableNames = {'Sbj_number'; 'Mode1_r'; 'Mode2_r'; 'Mode3_r';
+    'Mode1_R2'; 'Mode2_R2'; 'Mode3_R2'};
+writetable(QC_metrics_full,fullfile(pathtodata,'DCM_QC_1stRev_AllFemale.csv'),'WriteRowNames',true);
\ No newline at end of file

Codes/PEB_analysis.m

+%% CODE TO RUN PEB MODEL AND TEST MRS-GABA EFFECTS
+
+% This code needs the usual your_paths.txt file with all the info on your paths. 
+% An example is provided on GIT.
+
+% This code is an authomatized script to run PEB analysis with a regressor.
+% In this case we regress the MRS-derived metrics to the NMM output. We
+% will just focus on the T and H output of the NMM because they store
+% receptors and connections info respectively. 
+
+%-----------------------------------------------------------------------
+% Written by I. Paparella
+% Cyclotron Research Centre, University of Liege, Belgium
+% 14th December 2023
+
+
+function PEB_analysis
+
+
+%% Path info
+paths_txt = fileread('your_paths.txt') ;
+paths_txt = char(strsplit(paths_txt, '\n')) ;
+
+paths = [];
+paths.data = (strtrim(paths_txt(2,:)));
+paths.spm = (strtrim(paths_txt(4,:)));
+paths.codes = (strtrim(paths_txt(8,:)));
+
+addpath(paths.spm)
+%% Define folder to store analysis into
+PEB_path= strcat(paths.data, '/PEB');
+
+%% Start the PEB analysis
+% To test whether TMS-EEG and MRS GABA metrics are related we will
+% conduct a PEB analysis putting the MRS derived GABA metric as regressor
+% as described here: (https://en.wikibooks.org/wiki/SPM/Parametric_Empirical_Bayes_(PEB))
+
+
+% First let's create a folder to store the output of this PEB analysis
+PEB_MRS_NMM = strcat(PEB_path, '/PEB_MRS_NMM_FHK_AllFemale_GABA_withoutCycle');
+mkdir(PEB_MRS_NMM);
+
+% We load the GCMs 
+GCM_all = load(strcat(PEB_path, '/GCM_FHK_allFemale.mat'));
+
+% Load the MRS info
+MRS = load(strcat(paths.data, '/MRS_GABA_FHK_AllSubjects.txt'));
+
+% We define a design matrix X1 with regressors:
+% Mean
+% MRS-derived GABA
+% Interpret: as correlation (if positive = higher MRS, higher NMM; if
+% negative the other way around)
+    
+% Then we create the design matrix
+X1= [ones(size(GCM_all.GCM,1),1) MRS(:,2)];
+X1(:,2) = X1(:,2) - mean(X1(:,2)); % mean center the 2nd column so that the
+%first column can be taken as mean connectivity across subjects. If we
+%don't mean center the first column will be taken as baseline
+
+X1_labels = {'Mean_connectivity', 'MRS'};
+
+% Fit to all DCMs from old in a single column vector:
+GCM = GCM_all.GCM;
+
+% PEB setting and Run PEB
+M = struct();
+M.Q = 'all';
+M.X= X1;
+M.Xnames = X1_labels;
+
+[PEB ,RCM] = spm_dcm_peb(GCM, M, {'T', 'H'}); % need to set the field otherwise the
+% default will be matrix A and B (set on the MRI) which we don't have in
+% NMM
+
+% Save the output in the respective folder & run BMA
+save(strcat(PEB_MRS_NMM, '/PEB_MRS_NMM_FHK_allFemale_GABA.mat'),'PEB','RCM');
+
+BMA = spm_dcm_peb_bmc(PEB);
+save(strcat(PEB_MRS_NMM,'/BMA_PEB_MRS_NMM_FHK_allFemale_GABA.mat'),'BMA');
+
+% LOO cross-validation
+spm_dcm_loo(GCM,M,{'H(3,3)'}); % on IIfeedbackLoop
+spm_dcm_loo(GCM,M,{'H(3,4)'}); % on II-DP connection
+    
+end
+
+% spm_dcm_peb_review(BMA); to review results
\ No newline at end of file

Codes/PEB_analysis_withCycle.m

+%% CODE TO RUN PEB MODEL AND TEST MRS-GABA EFFECTS
+
+% This code needs the usual your_paths.txt file with all the info on your paths. 
+% An example is provided on GIT.
+
+% This code is an authomatized script to run PEB analysis with a regressor.
+% In this case we regress the MRS-derived metrics to the NMM output. We
+% will just focus on the T and H output of the NMM because they store
+% receptors and connections info respectively. 
+
+%-----------------------------------------------------------------------
+% Written by I. Paparella
+% Cyclotron Research Centre, University of Liege, Belgium
+% 26th February 2025
+
+
+function PEB_analysis_withCycle
+
+
+%% Path info
+paths_txt = fileread('your_paths.txt') ;
+paths_txt = char(strsplit(paths_txt, '\n')) ;
+
+paths = [];
+paths.data = (strtrim(paths_txt(2,:)));
+paths.spm = (strtrim(paths_txt(4,:)));
+paths.codes = (strtrim(paths_txt(8,:)));
+
+addpath(paths.spm)
+%% Define folder to store analysis into
+PEB_path= strcat(paths.data, '/PEB');
+
+%% Start the PEB analysis
+% To test whether TMS-EEG and MRS GABA metrics are related we will
+% conduct a PEB analysis putting the MRS derived GABA metric as regressor
+% as described here: (https://en.wikibooks.org/wiki/SPM/Parametric_Empirical_Bayes_(PEB))
+
+
+% First let's create a folder to store the output of this PEB analysis
+PEB_MRS_NMM = strcat(PEB_path, '/PEB_MRS_NMM_FHK_AllFemale_GABA_withCycle');
+mkdir(PEB_MRS_NMM);
+
+% We load the GCMs 
+GCM_all = load(strcat(PEB_path, '/GCM_FHK_allFemale.mat'));
+GCM_all.GCM{11,1}= []; % remove STR7T_HS_015 -> no menstrual cycle due to contracceptive pill
+
+% Load the MRS info
+subjToExclude = 15; %015 who has no cycle due to contracceptive pill
+MRS = load(strcat(paths.data, '/MRS_GABA_FHK_AllSubjects.txt'));
+MRS(ismember(MRS(:,1), subjToExclude),:) = [];
+
+% Load table with cycle info
+MensCycle = load(strcat(paths.data, '/MensCycle.txt'));
+
+% We define a design matrix X1 with regressors:
+% Mean
+% MRS-derived GABA
+% Interpret: as correlation (if positive = higher MRS, higher NMM; if
+% negative the other way around)
+    
+% Then we create the design matrix
+X1= [ones(size(GCM_all.GCM,1),1) MRS(:,2) MensCycle];
+X1(:,2) = X1(:,2) - mean(X1(:,2)); % mean center the 2nd column so that the
+%first column can be taken as mean connectivity across subjects. If we
+%don't mean center the first column will be taken as baseline
+X1(:,3) = X1(:,3) - mean(X1(:,3));
+
+X1_labels = {'Mean_connectivity', 'MRS', 'Cycle'};
+
+% Fit to all DCMs from old in a single column vector:
+GCM = GCM_all.GCM;
+
+% PEB setting and Run PEB
+M = struct();
+M.Q = 'all';
+M.X= X1;
+M.Xnames = X1_labels;
+
+[PEB ,RCM] = spm_dcm_peb(GCM, M, {'T', 'H'}); % need to set the field otherwise the
+% default will be matrix A and B (set on the MRI) which we don't have in
+% NMM
+
+% Save the output in the respective folder & run BMA
+save(strcat(PEB_MRS_NMM, '/PEB_MRS_NMM_FHK_allFemale_GABA.mat'),'PEB','RCM');
+
+BMA = spm_dcm_peb_bmc(PEB);
+save(strcat(PEB_MRS_NMM,'/BMA_PEB_MRS_NMM_FHK_allFemale_GABA.mat'),'BMA');
+
+% LOO cross-validation
+spm_dcm_loo(GCM,M,{'H(3,3)'}); % on IIfeedbackLoop
+spm_dcm_loo(GCM,M,{'H(3,4)'}); % on II-DP connection
+    
+end
+
+% spm_dcm_peb_review(BMA); to review results
\ No newline at end of file

Codes/RegressionPlots.R

+#### STR7T_HS -- CORRELATION ANALYSIS TMS-EEG DERIVED GABA MEASURES AND MRS ONES
+## -------------------------------------------------------------------------------------------------------------------
+## Created: Thursday Dec 21 14:37 2023 -- author: ipaparella@uliege.be
+# To check citation (for referencing) for any of the package used in this script just run citation("package")
+## -------------------------------------------------------------------------------------------------------------------
+
+
+### WHAT IS THIS SCRIPT DOING?
+## 1. Regression plots of the relationships shown by PEB analysis.
+## --------------------------------------------------------------------------------------------------------------------
+## --------------------------------------------------------------------------------------------------------------------
+## --------------------------------------------------------------------------------------------------------------------
+## --------------------------------------------------------------------------------------------------------------------
+
+# Import Libraries
+library("ggpubr")
+library("tidyverse")
+library("ggplot2")
+library(viridis)
+
+# Import data
+GABAmetrix_onlyFem <- read.csv("/mnt/data/STR7T_HS/Prepared/STR7T_NMM_perSess_DefaultCodes_AllSubjects.txt")
+
+### 1. 
+
+EI_MRS = GABAmetrix_onlyFem$MRS_GABA/GABAmetrix_onlyFem$MRS_GLU
+data.frame(EI_MRS)
+GABAmetrix_onlyFem = cbind(GABAmetrix_onlyFem, EI_MRS)
+
+# Average of tonic and phasic inhib (just the significant values)
+TonicInhib_sign = rowMeans(subset(GABAmetrix_onlyFem, select = c(SC, II, DP)), na.rm = TRUE)
+PhasicInhib_sign = rowMeans(subset(GABAmetrix_onlyFem, select = c(II.DP)), na.rm = TRUE)
+# Add this to table
+GABAmetrix_onlyFem = cbind(GABAmetrix_onlyFem, TonicInhib_sign, PhasicInhib_sign)
+
+# Check colors in virdis palette
+# library(scales)
+# v_colors =  viridis(30, option = 'D') -> number means how many colors you want to see
+# show_col(v_colors)
+# Plots
+ggplot(GABAmetrix_onlyFem, aes(x=MRS_GABA, y=GABAa)) +
+  geom_point(alpha=0.7, color="#3C508BFF", size = 4) +
+  geom_smooth(method=lm , color="#3C508BFF", fill="#3C508BFF", se=TRUE)+
+  theme_classic()+
+  xlab("GABA:Cr") +
+  ylab("GABAa")+
+  font("xy.text", size = 20, color = "black")+
+  scale_y_continuous(labels = scales::label_number(accuracy = 0.1))
+
+# 
+GABA_MRS_Tonic = cor.test(GABAmetrix_onlyFem$TonicInhib_sign, GABAmetrix_onlyFem$MRS_GABA, 
+                          method="spearman",exact=FALSE)
+GABA_MRS_Tonic
+ggplot(GABAmetrix_onlyFem, aes(x=MRS_GABA, y=TonicInhib_sign)) +
+  geom_point(alpha=0.7, color="#1F948CFF", size = 4) +
+  geom_smooth(method=lm , color="#1F948CFF", fill="#1F948CFF", se=TRUE)+
+  theme_classic()+
+  xlab("GABA:Cr") +
+  ylab("Tonic Inhibition")+
+  font("xy.text", size = 20, color = "black")
+
+#### -> 3rd round revision menstrual cycle effect. Plot linear regression with menstrual cycle category
+# Create extra column with cycle factor
+Cycle=c('M', 'M', 'M', 'M', 'L', 'L', 'M', 'F', 'L', 'F', 'None', 'M', 'F', 'F', 'M', 'F', 'F', 'M', 'L', 'L')
+data.frame(Cycle)
+GABAmetrix_onlyFem_cycle = cbind(GABAmetrix_onlyFem, Cycle)
+
+# Plots
+ggplot(GABAmetrix_onlyFem_cycle, aes(x=MRS_GABA, y=GABAa, col=Cycle)) +
+  geom_point(alpha=0.7, size = 4) +
+  geom_smooth(method=lm , color="#3C508BFF", fill="#3C508BFF", se=TRUE)+
+  theme_classic()+
+  xlab("GABA:Cr") +
+  ylab("GABAa")+
+  font("xy.text", size = 20, color = "black")+
+  scale_y_continuous(labels = scales::label_number(accuracy = 0.1))
+
+# 
+ggplot(GABAmetrix_onlyFem_cycle, aes(x=MRS_GABA, y=TonicInhib_sign, col=Cycle)) +
+  geom_point(alpha=0.7, size = 4) +
+  geom_smooth(method=lm , color="#1F948CFF", fill="#1F948CFF", se=TRUE)+
+  theme_classic()+
+  xlab("GABA:Cr") +
+  ylab("Tonic Inhibition")+
+  font("xy.text", size = 20, color = "black")