Descriptions and example scripts for running network preprocessing and analysis functions contained in this repository. See paper for details when it comes out.
Because of the known effect of motion on measures of connectivity, we followed up standard preprocessing in FSL with an extended nuisance regression. Affine transformation parameters from motion correction, CSF, white matter, and whole-brain signals are regressed against preprocessed 4D data, along with the squares, derivatives, and squared derivatives of these confounds. See Satterthwaite et al 2013 for details. Bash code:
for i in /data/engine/rgerraty/learn_dyncon/4*/Learn*/filtered_func_data.nii.gz
do
subdir=$(dirname $i)
#extract confound timecourses from preprocessed data
#need to provide feat directory as well as anatomical directory
#can also provide z-score cut-off for high-motion timepoints (spikes)
~/GitHub/rl_flexibility/fsl_extract_confts.sh $subdir\
$subdir/../structural/mprage.anat 3
#run 1st level confound regression, using template .fsf and confound files
~/GitHub/rl_flexibility/1st_level_conf.sh $i $subdir/36par+spikes.txt
doneAfter nuisance regression has been run, the residual timeseries needs to be transformed into standard space (in this case, MNI). Make sure fsl_anat has been run on each structural image first. The following bash code was used to perform these transformations:
#fnirt has already been run, just applying transformation
for i in /data/engine/rgerraty/learn_dyncon/4*/Learn*;
do
#run linear registration on example functional image
flirt -ref $i/../structural/mprage.anat/T1_biascorr_brain.nii.gz\
-in $i/reg/example_func.nii.gz\
-omat $i/reg/example_func2highres.mat;
echo warping $i;
#apply warp from FNIRT to preprocessed 4D data
applywarp --ref=$FSLDIR/data/standard/MNI152_T1_2mm.nii.gz\
--in=$i/36par+spikes.feat/stats/res4d.nii.gz\
--out=$i/36par+spikes.feat/stats/res4d_std.nii.gz\
--warp=$i/../structural/mprage.anat/T1_to_MNI_nonlin_field.nii.gz\
--premat=$i/reg/example_func2highres.mat;
doneOnce the preprocessed images have been registered, we extract mean timecourses for each Harvard-Oxford ROI, using the function extract_ROIs.sh. The output of this function is a timecourse for each ROI in the specified input folder, as well as a .txt file containing all of the ROIs. The bash code used to run this function on each learning block for each subject is below:
for i in /data/engine/rgerraty/learn_dyncon/4*/Learn?_PEprior.feat/36par+spikes.feat/;
do
#extract timeseries (mean or 1st eigenvector, see function) data from each ROI in ~/Harvard-Oxford_ROIs/
~/GitHub/rl_flexibility/extract_ROIs.sh $i/stats/res4d_std.nii.gz ~/Harvard-Oxford_ROIs/ $i/H-O_rois/;
doneConnectivity between pairs of ROIs was measured by average magnitude squared coherence in the .06-.12 Hz band, computed in MATLAB. The code below calls a function for creating a coherence matrices in a specified frequency range for specified time windows (in this case 25 TRs, or 50 s). These are saved as a .mat file for multi-slice community detection.
addpath ~/GitHub/rl_flexibility
%read in all subject/run ROI timeseries directories
[a,b]=system('ls -d /data/engine/rgerraty/learn_dyncon/4*/Learn?_PEprior.feat/36par+spikes.feat/H-O_rois');
c=strread(b,'%s');
for i=1:size(c,1)
%calculate coherence per time window from concatenated ROI file
filename=char(strcat(c(i),'/all_rois.txt'))
%need to specify filename, window length in TR, sampling rate, bandpass
conn_cell=coherence_by_block(filename,25,.5,.06,.12);
save(char(strcat(c(i),'/conn_cells')),'conn_cell')
endInput coherence matrix for each block. Also need number of blocks, resolution and coupling parameters. In Matlab
%need multi-slice, flexibility codes not yet on GitHub for network_diags to run
addpath ~/GitHub/rl_flexibility
addpath ~/GitHub/rl_flexibility/GenLouvain/
addpath ~/GitHub/rl_flexibility/Bassett_Code/
%read in data
[a,b]=system('ls -d /data/engine/rgerraty/learn_dyncon/4*/Learn?_PEprior.feat/36par+spikes.feat/H-O_rois/');
c=strread(b,'%s');
%concatenate runs for each subject
numruns=4
k=1;
for j=1:size(c,1)/numruns
c(k)
conn_cell_cat=[];
for i=1:numruns
load(strcat(char(c(k-1+i)),'/conn_cells'))
conn_cell_cat=cat(3,conn_cell_cat,conn_cell)
end
%network_diags code:
%runs multi-slice community detection
%gives flexibility for each run
%also allegiance matrix (not using yet)
%need to specify number of blocks, simulations, coupling, resolution
[a_mat,flex]=network_diags(conn_cell_cat,4,500,1,1.1813)
save(char(strcat(c(k),'/../../../a_mat')),'a_mat')
save(char(strcat(c(k),'/../../../flex')),'flex')
k=k+numruns;
endFor plotting and preparing for heirarchical models. Matlab.
%load data and concatenate flexibility statistics
[a,b]=system('ls -d /data/engine/rgerraty/learn_dyncon/4*/flex.mat');
c=strread(b,'%s');
flex_cat=[];
for j=1:size(c,1)
load(char(c(j)))
flex_cat=cat(3,flex_cat,flex)
end
plot(squeeze(mean(flex_cat)))
block=repmat([1:4]',22,1);
sub=repmat([1:22]',1,4)'
sub=sub(:);
%reshape whole-brain average flexibility
meanflex=squeeze(mean(flex_cat));
meanflex=meanflex(:);
%get striatal average flexibility
%check to make sure indices are correct
[trash,roi_names]=system('ls ~/Harvard-Oxford_ROIs/*nii.gz | xargs -n1 basename');
roi_names=strread(roi_names,'%s');
str_ind=[49,51,54,104,106,109];
roi_names(str_ind)
strflex=squeeze(mean(flex_cat(str_ind,:,:)));
strflex=strflex(:);
plot(squeeze(mean(flex_cat(str_ind,:,:))))
%write out csv for modeling in R
flexdata=[sub block meanflex strflex]
dlmwrite('/data/engine/rgerraty/learn_dyncon/flexdata.csv',flexdata)
%get flexibility scores for each ROI for each run for whole-brain search
%can prob do this more effeciently but this is easier to see, harder to botch
flex_allrois=[];
for i=1:size(flex_cat,3)
flex_allrois=[flex_allrois;flex_cat(:,:,i)'];
end
dlmwrite('/data/engine/rgerraty/learn_dyncon/flex_allrois.csv',flex_allrois) Estimate the effect of striatal and whole-brain flexibility on reinforcement learning. See models.Rmd and models.pdf for more details. R
library(reshape2)
library(lme4)
library(brms)
#prepare data for binomial logistic modelling
library(lme4)
library(reshape2)
library(ggplot2)
flexdat<-read.csv("/data/engine/rgerraty/learn_dyncon/flex_allrois.csv",header=F)
#read in trial-by-by trial behavioral data
data<-read.delim('/data/engine/rgerraty/learn_dyncon/behav_data.tsv',header=1)
data$block<-rep(rep(seq(1,4,1),each=30),22)
flexdat$subject<-rep(seq(1,22,1),each=4)
flexdat$Block<-rep(seq(1,4,1),times=22)
flexdat$Corr<-melt(tapply(data$optCor,list(data$block,data$subjectNum),mean,na.rm=1))$value
nacount<-is.na(data$optCor)
weights<-melt(tapply(nacount,list(data$block,data$subjectNum),sum))
flexdat$weights<-30-weights[,3]
str_ind=c(49,51,54,104,106,109);
flexdat<-melt(flexdat,id=c("subject","Block","Corr","weights"))
names(flexdat)[c(5,6)]<-c("ROI","flex")
flexdat$ROI<-as.factor(as.numeric(flexdat$ROI))
flexdat_str<-subset(flexdat,ROI %in% str_ind)
flexdat_str$ROI<-as.factor(flexdat_str$ROI)
flexdat_str_avg<-melt(tapply(flexdat_str$flex,list(flexdat_str$subject,flexdat_str$ROI),mean,na.rm=1))
names(flexdat_str_avg)<-c("subject","ROI","flex")
flexdat_str_avg<-subset(flexdat_str_avg,ROI %in% str_ind)
flexdat_str$meanflex<-rep(flexdat_str_avg$flex,each=4)
flexdat_str_melt<-melt(tapply(flexdat_str$flex,list(flexdat_str$subject,flexdat_str$Block),mean))
names(flexdat_str_melt)<-c("subject","Block","flex")
flexdat_str_melt$correct<-melt(tapply(flexdat_str$Corr,list(flexdat_str$subject,flexdat_str$Block),mean))$value
flexdat_str_melt$weights<-melt(tapply(flexdat_str$weights,list(flexdat_str$subject,flexdat_str$Block),mean))$value
m_str<-glmer(correct~flex+(flex||subject),
family=binomial,weights=weights,
data=flexdat_str_melt)
write.csv(flexdat_str_melt,'/data/engine/rgerraty/learn_dyncon/flex_behav.csv')
#for posterior inference, run bayesian model using brms wrapper for stan
options(mc.cores = parallel::detectCores())
flexdat_str_melt$numcorr<-as.integer(flexdat_str_melt$correct*flexdat_str_melt$weights)
mlearn_stan<-brm(numcorr~flex+(flex|subject),data=flexdat_str_melt,family=binomial)
R
#load in data
library(lme4)
roi_data<-read.csv('/data/engine/rgerraty/learn_dyncon/flex_allrois.csv',header=0)
behav<-read.csv('/data/engine/rgerraty/learn_dyncon/flex_behav.csv')
p<-NULL
b<-NULL
for(i in 1:dim(roi_data)[2]){
mtmp<-summary(glmer(behav$correct~roi_data[,i]+(roi_data[,i]||behav$subject),family=binomial,weights=behav$weights))
if(dim(summary(mtmp)[10]$coefficients)[2]>=4){
p[i]<-summary(mtmp)[10]$coefficients[2,4]
b[i]<-summary(mtmp)$coefficients[2,1]
} else{
p[i]<-NA
b[i]<-NA
}
}
write(p,'/data/engine/rgerraty/learn_dyncon/p_vals_learning_glm.csv')
write(b,'/data/engine/rgerraty/learn_dyncon/b_vals_learning_glm.csv')Matlab
str_ind=[49,51,54,104,106,109];
[a,b]=system('ls -d /data/engine/rgerraty/learn_dyncon/4*/a_mat.mat');
c=strread(b,'%s');
h=1
for s=1:length(c)
clear a_mat clear a_mat_str
load(c{s})
a_mat(a_mat==1)=NaN;
a_mat_str=a_mat(:,str_ind,:);
for i=1:size(a_mat_str,1)
for j=1:size(a_mat_str,2)
for k=1:size(a_mat_str,3)
a_mat_long(h,:)=[a_mat_str(i,j,k),i,j,ceil(k/8),s];
h=h+1;
end
end
end
end
header={'allegiance','ROI','str_roi','block','sub'};
header2=sprintf('%s,',header{:});header2(end)=[];
dlmwrite('/data/engine/rgerraty/learn_dyncon/alleg_long.csv',...
header2,'')
dlmwrite('/data/engine/rgerraty/learn_dyncon/alleg_long.csv',...
a_mat_long,'-append','delimiter',',')