---
title: "A triple-network organization for the mouse brain"
author: "Joanes Grandjean"
output: github_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval=FALSE)
```

## Summary
In this section, we perform an optognetic analysis in mice implanted in the anterior insula cortex. We further compare this to optogenetics in the [dorsal raphe](https://doi.org/10.1038/s41467-018-08256-w) 


### Inputs
Preprocessed optoINS data
Preprocessed optoDR data


### Outputs
**Figure 3B, S7B**   
Group optogenetics map (assets/opto/thresh_thr_optoINS.nii.gz)  

**Figure 3DEF**   
Group optogenetics map as a function of behaviour (assets/opto/thresh_thr_behav.nii.gz)  
Contrast parameter estimates extracted from the posterior cingulate area (assets/opto/Behav_RSP.txt)

**Figure 3G**   
Neuroquerry for positive valance (assets/opto/human/positive_valance.nii.gz)   


## Detailed code 

### Loading the R main variables      

```{r load main variables and pacakges, eval=FALSE}
library('stringr')
library('ggplot2')
require('svglite')
library('nat')
library('oro.nifti')
library('reshape2')
library('extrafont')
library('multcomp')
library('parameters')
library('effectsize')
library('reshape2')

#set the same variables as for bash_env.sh

init_folder <- "/home/traaffneu/joagra/code/insula"
data_folder <- "/project/4180000.18"

study_folder <- paste(data_folder, "/rest", sep='')
roi_folder <- paste(data_folder, "/SBA_roi", sep='')
SBA.folder <- paste(study_folder, '/SBA/', sep = '')
viral.folder <- paste(study_folder, "tracer", sep='')

template_mask <- "assets/template/mask_200um.nii"
template <- "assets/template/ABI_template_200um.nii"
atlas <- "assets/template/ABI_atlas_reduced_TIDO_200um.nii.gz"
atlas.label <- "assets/template/ABI_atlas_reduced_TIDO.csv"

study <- read.csv(paste(study_folder,'/participants.tsv',sep='')
study <- study[!study$dataset %in% c('aes1', 'MEPSI', 'AD1', 'ME1'), ]
```


## Optogenetics fMRI
Datasets are publicly available in raw BIDs format from the openneuro.org resource. 
[Insula opto-fMRI](10.18112/openneuro.ds003464.v1.0.0)
[Dorsal raphe opto-fMRI](10.18112/openneuro.ds001541.v1.1.2)

Both datasets were preprocessed using the following [code](https://github.com/grandjeanlab/MouseMRIPrep).  
The code below assumes also that a GLM was fitted using the task files. 

This first set of chunks tests for the main effect of insula photostimulation. 
### Figure 3B

```{r opto GLM, eval=FALSE}


######### 1) Folder containing the glm images#######
# Do not add the final '/' to the directory
glm_dir <- ('/project/4180000.15/opto_ins/preprocessing')


########## 2) folder list containing the experiment information########
folder_list <- ('/project/4180000.15/opto_ins/study.csv')


########### 5) where you want to export the statistics #########

export_folder <-
  '/project/4180000.15/opto_ins/statistics/opto_select_20200609'

dir.create(export_folder, recursive = T, showWarnings = F)
file_prefix <- 'task-opto_acq-EPI'
file_sufix <- 'avg_cope.nii.gz'

# 5) exclude data
exclude <- T

#read the folder list
study_list <- read.csv(folder_list, header = T, sep = ',')

if (exclude) {
  study_list <- study_list[study_list$exclude == 'no', ]
}
#if(exclude){study_list<-study_list[study_list$response=='yes',]}

#study_list$condition_virus<-paste(study_list$condition,study_list$virus,sep='_')

levels(study_list$group) <- c('ctl', 'opto', 'opto', 'opto')


ID_factor <- levels(factor(study_list$ID))
Group_factor <- levels(factor(study_list$group))




ID_design <- c()
Group_design <- c()

ID_group <- c()

file_list <- c()

ii <- 1
for (w in 1:dim(study_list)[1])
{
  file_dir <-
    paste(
      glm_dir,
      '/sub-',
      study_list$ID[w],
      '/ses-',
      study_list$session[w],
      '/func/',
      file_sufix,
      sep = ''
    )
  
  
  if (file.exists(file_dir)) {
    file_list[ii] <- file_dir
    ID_vec <- rep(0, length(ID_factor))
    ID_vec[which(ID_factor == study_list$ID[w])] <- 1
    ID_design <- rbind(ID_design, ID_vec)
    
    file_list[ii] <- file_dir
    ID_vec <- rep(0, length(ID_factor))
    ID_vec[which(ID_factor == study_list$ID[w])] <- 1
    ID_design <- rbind(ID_design, ID_vec)
    
    Group_vec <- rep(0, length(Group_factor))
    Group_vec[which(Group_factor == study_list$group[w])] <- 1
    Group_design <- rbind(Group_design, Group_vec)
    
    ID_group[ii] <- which(ID_factor == study_list$ID[w])
    ii <- ii + 1
  } else{
    print(file_dir)
  }
  
  
}

final_design <- cbind(Group_design)
write.table(
  file_list,
  file.path(export_folder, 'scan_list.txt'),
  row.names = F,
  col.names = F,
  quote = F
)



K <- rbind(c(-1, 1),
           c(1, -1))

con_name <- c("ctl<opto", "ctl>opto")

desgin_mat <- c()
ii <- 1
desgin_mat[ii] <- paste('/NumWaves',  dim(final_design)[2], sep = "\t")
ii <- ii + 1
desgin_mat[ii] <- paste('/NumPoints',  dim(final_design)[1], sep = "\t")
ii <- ii + 1
desgin_mat[ii] <-
  paste('/PPheights\t', paste(rep('1', dim(final_design)[1]), collapse = '\t'), sep =
          "\t")
ii <- ii + 1
desgin_mat[ii] <- ''
ii <- ii + 1
desgin_mat[ii] <- paste('/Matrix', sep = "\t")
for (i in 1:dim(final_design)[1]) {
  ii <- ii + 1
  desgin_mat[ii] <-
    paste(as.character(final_design[i, ]), collapse = '\t')
  
}


write.table(
  desgin_mat,
  file.path(export_folder, 'design.mat'),
  quote = F,
  row.names = F,
  col.names = F
)



desgin_grp <- c()
ii <- 1
desgin_grp[ii] <- paste('/NumWaves',  '1', sep = "\t")
ii <- ii + 1
desgin_grp[ii] <- paste('/NumPoints',  dim(final_design)[1], sep = "\t")
ii <- ii + 1
desgin_grp[ii] <- ''
ii <- ii + 1
desgin_grp[ii] <- paste('/Matrix', sep = "\t")
for (i in 1:length(ID_group)) {
  ii <- ii + 1
  desgin_grp[ii] <- paste(as.character(ID_group[i]), collapse = '\t')
  
}


write.table(
  desgin_grp,
  file.path(export_folder, 'design.grp'),
  quote = F,
  row.names = F,
  col.names = F
)






desgin_con <- c()
ii <- 1
for (i in 1:dim(K)[1]) {
  desgin_con[ii] <-
    paste(paste('/ContrastName', i, sep = ''),  con_name[i], sep = "\t")
  ii <- ii + 1
}
desgin_con[ii] <- paste('/NumWaves',  dim(K)[2], sep = "\t")
ii <- ii + 1
desgin_con[ii] <- paste('/NumContrasts',  dim(K)[1], sep = "\t")
ii <- ii + 1
desgin_con[ii] <-
  paste('/PPheights\t', paste(rep('1', dim(K)[1]), collapse = '\t'), sep =
          "\t")
ii <- ii + 1
desgin_con[ii] <-
  paste('/RequiredEffect\t', paste(rep('1', dim(K)[1]), collapse = '\t'), sep =
          "\t")
ii <- ii + 1
desgin_con[ii] <- ''
ii <- ii + 1
desgin_con[ii] <- paste('/Matrix', sep = "\t")
for (i in 1:dim(K)[1]) {
  ii <- ii + 1
  desgin_con[ii] <- paste(as.character(K[i, ]), collapse = '\t')
  
}

write.table(
  desgin_con,
  file.path(export_folder, 'design.con'),
  quote = F,
  row.names = F,
  col.names = F
)

ii <- 1
script.randomize <- c()
script.randomize[ii] <- paste('#!/bin/bash', sep = "\t")
ii <- ii + 1
script.randomize[ii] <- paste('cd', export_folder, sep = ' ')
ii <- ii + 1
script.randomize[ii] <-
  paste('fslmerge -t merge $(cat scan_list.txt)', sep = "\t")
ii <- ii + 1
script.randomize[ii] <-
  paste('fslmeants -i merge.nii.gz -o COPE -m ',
        template_mask,
        ' --label=\"',
        atlas,
        '\"',
        sep = "")
ii <- ii + 1
script.randomize[ii] <-
  paste(
    'fsl_glm -i merge.nii.gz --out_z=glm_z --out_p=glm_p -d design.mat -c design.con -m',
    template_mask,
    sep = " "
  )
ii <- ii + 1
script.randomize[ii] <- paste('fslsplit glm_z', sep = " ")
ii <- ii + 1
script.randomize[ii] <-
  paste('easythresh vol0000.nii.gz',
        template_mask,
        '1.92 0.05',
        template ,
        'thr_pos',
        sep = " ")
ii <- ii + 1
script.randomize[ii] <-
  paste('easythresh vol0001.nii.gz',
        template_mask,
        '1.92 0.05',
        template ,
        'thr_neg',
        sep = " ")
ii <- ii + 1
script.randomize[ii] <-
  paste(
    'fslmaths thresh_thr_neg.nii.gz -mul -1 -add thresh_thr_pos.nii.gz thresh_thr.nii.gz',
    sep = " "
  )



write.table(
  script.randomize,
  file.path(export_folder, 'randomize.sh'),
  quote = F,
  row.names = F,
  col.names = F
)
```

```{bash evalutate opto GLM test, eval=FALSE}
source bash_env.sh

cd /project/4180000.15/opto_ins/statistics/opto_select_20200609
chmod +x ${PWD}/randomize.sh

#because we have access to a cluster system, we submit jobs. 
qsub -l 'procs=1,mem=16gb,walltime=04:00:00' ${PWD}/randomize.sh
#if a cluster is not available, comment the line above and use the following instead. Might be a tad slow to run. 
#bash ${PWD}/randomize.sh

fslmeants -i glm_z.nii.gz -o $template_mask -o $init_folder'/assets/opto/opto_z' --label=$atlas

```

This set of chunks test for insula photostimulation X behavioural response interactions. 
### Figure 3E
```{r opto behaviour glm, eval=FALSE}
glm_dir <- ('/project/4180000.15/opto_ins/preprocessing')


########## 2) folder list containing the experiment information########
folder_list <- ('/project/4180000.15/opto_ins/study.csv')


########### 5) where you want to export the statistics #########

export_folder <-
  '/project/4180000.15/opto_ins/statistics/opto_behav_select_20200609'

dir.create(export_folder, recursive = T, showWarnings = F)
file_prefix <- 'task-opto_acq-EPI'
file_sufix <- 'avg_cope.nii.gz'

# 5) exclude data
exclude <- T

#read the folder list
study_list <- read.csv(folder_list, header = T, sep = ',')

if (exclude) {
  study_list <- study_list[study_list$exclude == 'no', ]
}
if (exclude) {
  study_list <- study_list[!is.na(study_list$delta.pref), ]
}
#if(exclude){study_list<-study_list[study_list$exclude=='no',]}
#if(exclude){study_list<-study_list[study_list$group %in% c('Ins_front','ctl'),]}
#study_list$condition_virus<-paste(study_list$condition,study_list$virus,sep='_')


levels(study_list$group) <- c('ctl', 'opto', 'opto', 'opto')

ID_factor <- levels(factor(study_list$ID))
Group_factor <- levels(factor(study_list$group))

ID_design <- c()

Behaviour_design <- c()

ID_group <- c()

file_list <- c()

ii <- 1
for (w in 1:dim(study_list)[1])
{
  file_dir <-
    paste(
      glm_dir,
      '/sub-',
      study_list$ID[w],
      '/ses-',
      study_list$session[w],
      '/func/',
      file_sufix,
      sep = ''
    )
  
  
  if (file.exists(file_dir) &
      !is.na(study_list$delta.pref[w]) & study_list$group[w] == 'opto') {
    file_list[ii] <- file_dir
    ID_vec <- rep(0, length(ID_factor))
    ID_vec[which(ID_factor == study_list$ID[w])] <- 1
    ID_design <- rbind(ID_design, ID_vec)
    
    
    Behaviour_design <- rbind(Behaviour_design, study_list$delta.pref[w])
    ii <- ii + 1
    ID_group[ii] <-
      which(ID_factor == study_list$ID[w])
  } else{
    print(file_dir)
  }
  
  
}

final_design <- cbind(Behaviour_design)
write.table(
  file_list,
  paste(export_folder, '/', 'scan_list.txt', sep = ''),
  row.names = F,
  col.names = F,
  quote = F
)

##set design now just for  Group_design
#Group_factor
#


K <- rbind(c(1), c(-1))

con_name <- c("Behav+", "Behav-")


desgin_mat <- c()
ii <- 1
desgin_mat[ii] <- paste('/NumWaves',  dim(final_design)[2], sep = "\t")
ii <- ii + 1
desgin_mat[ii] <- paste('/NumPoints',  dim(final_design)[1], sep = "\t")
ii <- ii + 1
desgin_mat[ii] <-
  paste('/PPheights\t', paste(rep('1', dim(final_design)[1]), collapse = '\t'), sep =
          "\t")
ii <- ii + 1
desgin_mat[ii] <- ''
ii <- ii + 1
desgin_mat[ii] <- paste('/Matrix', sep = "\t")
for (i in 1:dim(final_design)[1]) {
  ii <- ii + 1
  desgin_mat[ii] <-
    paste(as.character(final_design[i, ]), collapse = '\t')
  
}


write.table(
  desgin_mat,
  file.path(export_folder, 'design.mat'),
  quote = F,
  row.names = F,
  col.names = F
)



desgin_grp <- c()
ii <- 1
desgin_grp[ii] <- paste('/NumWaves',  '1', sep = "\t")
ii <- ii + 1
desgin_grp[ii] <- paste('/NumPoints',  dim(final_design)[1], sep = "\t")
ii <- ii + 1
desgin_grp[ii] <- ''
ii <- ii + 1
desgin_grp[ii] <- paste('/Matrix', sep = "\t")
for (i in 1:length(ID_group)) {
  ii <- ii + 1
  desgin_grp[ii] <- paste(as.character(ID_group[i]), collapse = '\t')
  
}


write.table(
  desgin_grp,
  file.path(export_folder, 'design.grp'),
  quote = F,
  row.names = F,
  col.names = F
)






desgin_con <- c()
ii <- 1
for (i in 1:dim(K)[1]) {
  desgin_con[ii] <-
    paste(paste('/ContrastName', i, sep = ''),  con_name[i], sep = "\t")
  ii <- ii + 1
}
desgin_con[ii] <- paste('/NumWaves',  dim(K)[2], sep = "\t")
ii <- ii + 1
desgin_con[ii] <- paste('/NumContrasts',  dim(K)[1], sep = "\t")
ii <- ii + 1
desgin_con[ii] <-
  paste('/PPheights\t', paste(rep('1', dim(K)[1]), collapse = '\t'), sep =
          "\t")
ii <- ii + 1
desgin_con[ii] <-
  paste('/RequiredEffect\t', paste(rep('1', dim(K)[1]), collapse = '\t'), sep =
          "\t")
ii <- ii + 1
desgin_con[ii] <- ''
ii <- ii + 1
desgin_con[ii] <- paste('/Matrix', sep = "\t")
for (i in 1:dim(K)[1]) {
  ii <- ii + 1
  desgin_con[ii] <- paste(as.character(K[i, ]), collapse = '\t')
  
}

write.table(
  desgin_con,
  file.path(export_folder, 'design.con'),
  quote = F,
  row.names = F,
  col.names = F
)


ii <- 1
script.randomize <- c()
script.randomize[ii] <- paste('#!/bin/bash', sep = "\t")
ii <- ii + 1
script.randomize[ii] <- paste('cd', export_folder, sep = ' ')
ii <- ii + 1
script.randomize[ii] <-
  paste('fslmerge -t merge $(cat scan_list.txt)', sep = "\t")
ii <- ii + 1
script.randomize[ii] <-
  paste('fslmeants -i merge.nii.gz -o COPE -m ',
        template_mask,
        ' --label=\"',
        atlas,
        '\"',
        sep = "")
ii <- ii + 1
script.randomize[ii] <-
  paste(
    'fsl_glm -i merge.nii.gz --out_z=glm_z --out_p=glm_p -d design.mat -c design.con -m',
    template_mask,
    sep = " "
  )
ii <- ii + 1
script.randomize[ii] <- paste('fslsplit glm_z', sep = " ")
ii <- ii + 1
script.randomize[ii] <-
  paste('easythresh vol0000.nii.gz',
        template_mask,
        '1.92 0.05',
        template ,
        'thr_pos',
        sep = " ")
ii <- ii + 1
script.randomize[ii] <-
  paste('easythresh vol0001.nii.gz',
        template_mask,
        '1.92 0.05',
        template ,
        'thr_neg',
        sep = " ")
ii <- ii + 1
script.randomize[ii] <-
  paste(
    'fslmaths thresh_thr_neg.nii.gz -mul -1 -add thresh_thr_pos.nii.gz thresh_thr.nii.gz',
    sep = " "
  )

write.table(
  script.randomize,
  file.path(export_folder, 'randomize.sh'),
  quote = F,
  row.names = F,
  col.names = F
)
```
```{bash evalutate opto GLM, eval=FALSE}
source bash_env.sh

cd /project/4180000.15/opto_ins/statistics/opto_behav_select_20200609
chmod +x ${PWD}/randomize.sh

#because we have access to a cluster system, we submit jobs. 
qsub -l 'procs=1,mem=16gb,walltime=04:00:00' ${PWD}/randomize.sh
#if a cluster is not available, comment the line above and use the following instead. Might be a tad slow to run. 
#bash ${PWD}/randomize.sh

fslmeants -i merge.nii.gz -m /groupshare/traaffneu/preclinimg/templates/my_templates/mask_200um.nii -o COPE --label=cluster_mask_thr_pos.nii.gz

```

### Figure 3G
get the map from human (querry on 09/06/2020). We use the default threshold of 4-10, as per default on neuroquery database.
https://neuroquery.org/query?text=positive+valance

This plots the optogenetics evoked response
### Figure 3C
```{bash get time series, eval=FALSE}
bids='/project/4180000.15/opto_ins/preprocessing'
output='/project/4180000.15/opto_ins/TS_AIdv'


template='/home/traaffneu/joagra/my_templates/MRI_exvivo_template_200um_skullstripped.nii'
template_mask='/home/traaffneu/joagra/my_templates/mask_200um.nii'
#template_atlas='/home/traaffneu/joagra/my_templates/atlas/ABI_atlas_reduced_200um.nii.gz'
template_atlas='/project/4180000.15/opto_ins/l_AIdv_200um.nii.gz'

affine=reg/func2anat0GenericAffine.mat
warp2=reg/func2anat1Warp.nii.gz
warp1=reg/anat2abi.nii.gz

cd $bids
mkdir -p $output


ls -d */ | while read subject
do

cd $subject


ls -d */ | while read session
do

cd $session'func'

echo 'now doing '$subject' and '$session


ls -d *task*/ | while read func
do
cd $func
func_noext=${func%/}

outputfile=$output'/'${func%/}'_ts' 

outdir=$(pwd)
echo "cd "$outdir > $outdir/run.sh
echo "outputfile="${outputfile} >> $outdir/run.sh
echo "template="${template} >> $outdir/run.sh
echo "template_mask="${template_mask} >> $outdir/run.sh
echo "template_atlas="${template_atlas} >> $outdir/run.sh
echo "WarpTimeSeriesImageMultiTransform 4 filtered_func_data_clean.nii.gz filtered_func_data_clean_reg.nii.gz -R $template $warp1 $warp2 $affine" >> $outdir/run.sh
echo "fslmeants -i filtered_func_data_clean_reg.nii.gz -m $template_mask --label=$template_atlas -o $outputfile" >> $outdir/run.sh
echo "rm filtered_func_data_clean_reg.nii.gz" >> $outdir/run.sh

chmod +x $outdir/run.sh
qsub -l 'procs=1,mem=4gb,walltime=00:30:00' $outdir/run.sh

cd ..
done
cd ../..
done
cd ..
done

```


```{r plot trace, eval=FALSE}

study <- read.csv('/project/4180000.15/opto_ins/study.csv', header = T)
study.sub <- study[study$exclude == 'no', ]
data.dir <- '/project/4180000.15/opto_ins/TS_AIdv/'

dir.list <- dir(data.dir, full.names = T)

TS.ctl <- c()
TS.exp <- c()

for (i in 1:length(dir.list)) {
  ID <-
    strsplit(strsplit(dir.list[i], 'jgroptoINS')[[1]][2], '_ses')[[1]][1]
  ses <-
    strsplit(strsplit(dir.list[i], '_ses-')[[1]][2], '_task-opto')[[1]][1]
  
  TS.select <-
    intersect(which(study.sub$ID.original == ID),
              which(study.sub$session == ses))
  if (length(TS.select) == 1) {
    group <- study.sub$group[TS.select]
    TS <- read.table(dir.list[i])$V1
    if (group == 'ctl') {
      TS.ctl <- rbind(TS.ctl, t(TS))
    }
    if (group != 'ctl') {
      TS.exp <- rbind(TS.exp, t(TS))
    }
  }
}



TS.ctl.block <- array(0, dim = c(dim(TS.ctl)[1], 61))
TS.exp.block <- array(0, dim = c(dim(TS.exp)[1], 61))


TS.merge <- rbind(
  c(45:(45 + 60)),
  c(105:(105 + 60)),
  c(165:(165 + 60)),
  c(225:(225 + 60)),
  c(285:(285 + 60)),
  c(345:(345 + 60)),
  c(405:(405 + 60)),
  c(465:(465 + 60)),
  c(525:(525 + 60)),
  c(585:(525 + 60))
)

for (j in 1:dim(TS.merge)[1]) {
  TS.ctl.block <- TS.ctl.block + TS.ctl[, TS.merge[j, ]]
  TS.exp.block <- TS.exp.block + TS.exp[, TS.merge[j, ]]
}

TS.ctl.block <- TS.ctl.block / 10
TS.exp.block <- TS.exp.block / 10

TS.ctl.block <-
  (TS.ctl.block / apply(TS.ctl.block[, 1:5], 1, mean) * 100) - 100
TS.exp.block <-
  (TS.exp.block / apply(TS.exp.block[, 1:5], 1, mean) * 100) - 100

TS.ctl.block.mean <- colMeans(TS.ctl.block)
TS.exp.block.mean <- colMeans(TS.exp.block)

TS.ctl.block.uppper <-
  colMeans(TS.ctl.block)#+apply(TS.ctl.block,2,sd)
TS.exp.block.upper <- colMeans(TS.exp.block) + apply(TS.exp.block, 2, sd)

TS.ctl.block.lower <- colMeans(TS.ctl.block) - apply(TS.ctl.block, 2, sd)
TS.exp.block.lower <- colMeans(TS.exp.block)#-apply(TS.exp.block,2,sd)

time <- c(-5:55)

dd <-
  data.frame(
    c(rep('ctl', length(time)), rep('exp', length(time))),
    c(time, time),
    c(TS.ctl.block.mean, TS.exp.block.mean),
    c(TS.ctl.block.uppper, TS.exp.block.upper),
    c(TS.ctl.block.lower, TS.exp.block.lower)
  )
names(dd) <- c('group', 'time', 'BOLD', 'upper', 'lower')

p <-
  ggplot(dd, aes(
    x = time,
    y = BOLD,
    group = group,
    colour = group,
    fill = group
  ))
p + geom_hline(yintercept = 0) + geom_rect(
  xmin = 0,
  xmax = 10,
  ymin = -1,
  ymax = -0.1,
  colour = '#add8e6',
  alpha = 0.005
) +
  geom_line(size = 0.5) +
  geom_errorbar(aes(ymax = upper, ymin = lower),
                width = 0,
                size = 0.2) +
  theme_classic() + theme(legend.position = 'none', text = element_text(size =
                                                                          6)) +
  scale_color_manual(values = c('#44708F', '#CB650B'))
ggsave(
  'assets/opto/trace/trace_plot.svg',
  device = 'svg',
  units = 'mm',
  dpi = 300,
  width = 41,
  height = 31
)
ggsave(
  'assets/opto/trace/trace_plot.png',
  device = 'png',
  units = 'mm',
  dpi = 300,
  width = 41,
  height = 31
)

```
![Figure 3c](assets/opto/trace/trace_plot.svg)

This chunks exports the behaviour parameters for plotting with DABEST and the linear model with the parameter estimates from the posterior cingulate area 
### Figure 3F
```{r plot opto COPE, eval=FALSE}
folder_list <- ('/project/4180000.15/opto_ins/study.csv')
ROI_list <- read.csv2(atlas.label)

study_list <- read.csv(folder_list, header = T, sep = ',')
levels(study_list$group) <- c('ctl', 'opto', 'opto', 'opto')

study_list.allCOPE <- study_list[study_list$exclude == 'no', ]



study_list.allCOPE$insula <-
  read.table('/project/4180000.15/opto_ins/statistics/opto_select_20200609/COPE')$V19
COPE_export_dabest <-
  dcast(study_list.allCOPE, insula ~ group, value.var = "insula")
COPE_export_dabest <- COPE_export_dabest[, -1]
write.csv(COPE_export_dabest,
          'assets/opto/COPE_insula.csv',
          row.names = FALSE)




study_list.behavCOPE <-
  study_list.allCOPE[!is.na(study_list.allCOPE$delta.pref), ]
Behav_export_dabest <-
  dcast(study_list.behavCOPE, delta.pref ~ group, value.var = "delta.pref")
Behav_export_dabest <- Behav_export_dabest[, -1]
write.csv(Behav_export_dabest, 'assets/opto/Behav.csv', row.names = FALSE)

study_list.behavCOPE.opto <-
  study_list.behavCOPE[study_list.behavCOPE$group == 'opto', ]
study_list.behavCOPE.opto$RSP <-
  read.table('/project/4180000.15/opto_ins/statistics/opto_behav_select_20200609/COPE')$V78
study_list.behavCOPE.opto$ACA <-
  read.table('/project/4180000.15/opto_ins/statistics/opto_behav_select_20200609/COPE')$V76


p <- ggplot(study_list.behavCOPE.opto, aes(x = ACA, y = delta.pref))
p + geom_point(colour = '#CB650B') +
  geom_smooth(method = "lm",
              formula = y ~ x,
              colour = '#CB650B') + theme_classic() + labs(title = 'Anterior cingulate area', y =
                                                             'Conditioned place preference', x = 'COPE')
s <- parameters(lm(delta.pref ~ ACA, study_list.behavCOPE.opto),summary  = TRUE, standardize = 'refit')
capture.output(s, file = "assets/opto/Behav_ACA.txt")

ggsave(
  file = "assets/opto/Behav_ACA.svg",
  device = 'svg',
  width = 75,
  height = 75,
  units = 'mm',
  dpi = 300
)
ggsave(
  file = "assets/opto/Behav_ACA.png",
  device = 'png',
  width = 75,
  height = 75,
  units = 'mm',
  dpi = 300
)

p <- ggplot(study_list.behavCOPE.opto, aes(x = RSP, y = delta.pref))
p + geom_point(colour = '#CB650B') +
  geom_smooth(method = "lm",
              formula = y ~ x,
              colour = '#CB650B') + theme_classic() + labs(title = 'Posterior cingulate area', y =
                                                             'Conditioned place preference', x = 'COPE')
s <- parameters(lm(delta.pref ~ RSP, study_list.behavCOPE.opto),summary  = TRUE, standardize = 'refit')
capture.output(s, file = "assets/opto/Behav_RSP.txt")
s

ggsave(
  file = "assets/opto/Behav_RSP.svg",
  device = 'svg',
  width = 75,
  height = 75,
  units = 'mm',
  dpi = 300
)
ggsave(
  file = "assets/opto/Behav_RSP.png",
  device = 'png',
  width = 75,
  height = 75,
  units = 'mm',
  dpi = 300
)

```
![Figure 3f](assets/opto/Behav_RSP.svg)

### optogenetics dorsal raphe
This is the code for the preprocessing of optogenetics in the dorsal raphe nucleus. 
### Figure S7A

```{r make haemodynamic response model, eval=FALSE}
setwd('/project/4180000.18/optoDR/')
library('neuRosim')

contrast <- 1

subject.list <- dir('bids')
for (i in subject.list) {
  session.list <- (dir(file.path('bids', i)))
  for (j in session.list) {
    func.list <-
      (dir(
        file.path('bids', i, j, 'func'),
        pattern = 'events.tsv',
        full.names = TRUE
      ))
    bold.id <- grep(x = func.list, pattern = 'run-4')
    if (length(bold.id) >= 1) {
      func.list <- func.list[-c(bold.id)]
    }
    for (k in func.list) {
      print(k)
      design.table <- read.table(k, header = TRUE)
      
      
      model.convolve <-
        specifydesign(
          totaltime = 360 * 2,
          onsets = list(design.table$onset),
          durations = list(design.table$duration),
          effectsize = 1,
          TR = 2,
          conv = "Balloon"
        )
      out.dir <-
        file.path('preprocessing',
                  i,
                  j,
                  'func',
                  paste(sub('events.tsv', '', basename(k)), 'cbv', sep = ''),
                  'glm')
      dir.create(out.dir, recursive = TRUE)
      write.table(
        model.convolve,
        file.path(out.dir, 'design.txt'),
        quote = FALSE,
        row.names = FALSE,
        col.names = FALSE
      )
      write.table(
        contrast,
        file.path(out.dir, 'contrast.txt'),
        quote = FALSE,
        row.names = FALSE,
        col.names = FALSE
      )
      
    }
  }
}
```

After running the fsl_glm functions, generate a second-level analysis
```{r opto dorsal raphe glm, eval=FALSE }

setwd('/project/4180000.18/optoDR')

export_folder <- 'statistics'
dir.create(export_folder, recursive = TRUE, showWarnings = FALSE)


study <- read.table('participants.tsv', sep = '\t', header = TRUE)
study <- study[study$run1 != 'l100' & study$run1 != 'ul100', ]



file_list <- c()
design_ctl <- c()
design_exp <- c()
design_fluo <- c()
design_restrain <- c()


ii <- 1

for (i in 1:dim(study)[1]) {
  subject <- paste('sub-', study$participant_id[i], sep = '')
  session <- paste('ses-', study$session[i], sep = '')
  func1 <-
    paste(subject, session, 'task-opto', 'acq-EPI', 'run-1', 'cbv', sep = '_')
  func2 <-
    paste(subject, session, 'task-opto', 'acq-EPI', 'run-2', 'cbv', sep = '_')
  func3 <-
    paste(subject, session, 'task-opto', 'acq-EPI', 'run-3', 'cbv', sep = '_')
  
  cope1 <-
    file.path(
      getwd(),
      'preprocessing',
      subject,
      session,
      'func',
      func1,
      paste(func1, '_COPE_reg.nii.gz', sep = '')
    )
  cope2 <-
    file.path(
      getwd(),
      'preprocessing',
      subject,
      session,
      'func',
      func2,
      paste(func2, '_COPE_reg.nii.gz', sep = '')
    )
  cope3 <-
    file.path(
      getwd(),
      'preprocessing',
      subject,
      session,
      'func',
      func3,
      paste(func3, '_COPE_reg.nii.gz', sep = '')
    )
  
  
  
  if (file.exists(cope1)) {
    if (study$group[i] == 'con') {
      run <- 'ctl'
      design_ctl[ii] <- 1
      design_exp[ii] <- 0
      design_fluo[ii] <- 0
      design_restrain[ii] <- 0
    } else {
      if (study$run1[i] %in% c('33', '66', '100', 'g100', 'o100')) {
        run <-
          '100'
        design_ctl[ii] <-
          0
        design_exp[ii] <- 1
        design_fluo[ii] <- 0
        design_restrain[ii] <- 0
      }
      if (study$run1[i] %in% c('r100')) {
        run <-
          'r100'
        design_ctl[ii] <-
          0
        design_exp[ii] <- 0
        design_fluo[ii] <- 0
        design_restrain[ii] <- 1
      }
      if (study$run1[i] %in% c('f100')) {
        run <-
          'f100'
        design_ctl[ii] <-
          0
        design_exp[ii] <- 0
        design_fluo[ii] <- 1
        design_restrain[ii] <- 0
      }
    }
    
    
    file_list[ii] <- cope1
    ii <- ii + 1
    
  }
  if (file.exists(cope2)) {
    if (study$group[i] == 'con') {
      run <- 'ctl'
      design_ctl[ii] <- 1
      design_exp[ii] <- 0
      design_fluo[ii] <- 0
      design_restrain[ii] <- 0
    } else {
      if (study$run2[i] %in% c('33', '66', '100', 'g100', 'o100')) {
        run <-
          '100'
        design_ctl[ii] <-
          0
        design_exp[ii] <- 1
        design_fluo[ii] <- 0
        design_restrain[ii] <- 0
      }
      if (study$run2[i] %in% c('r100')) {
        run <-
          'r100'
        design_ctl[ii] <-
          0
        design_exp[ii] <- 0
        design_fluo[ii] <- 0
        design_restrain[ii] <- 1
      }
      if (study$run2[i] %in% c('f100')) {
        run <-
          'f100'
        design_ctl[ii] <-
          0
        design_exp[ii] <- 0
        design_fluo[ii] <- 1
        design_restrain[ii] <- 0
      }
    }
    
    
    file_list[ii] <- cope2
    ii <- ii + 1
    
  }
  if (file.exists(cope3)) {
    if (study$group[i] == 'con') {
      run <- 'ctl'
      design_ctl[ii] <- 1
      design_exp[ii] <- 0
      design_fluo[ii] <- 0
      design_restrain[ii] <- 0
    } else {
      if (study$run3[i] %in% c('33', '66', '100', 'g100', 'o100')) {
        run <-
          '100'
        design_ctl[ii] <-
          0
        design_exp[ii] <- 1
        design_fluo[ii] <- 0
        design_restrain[ii] <- 0
      }
      if (study$run3[i] %in% c('r100')) {
        run <-
          'r100'
        design_ctl[ii] <-
          0
        design_exp[ii] <- 0
        design_fluo[ii] <- 0
        design_restrain[ii] <- 1
      }
      if (study$run3[i] %in% c('f100')) {
        run <-
          'f100'
        design_ctl[ii] <-
          0
        design_exp[ii] <- 0
        design_fluo[ii] <- 1
        design_restrain[ii] <- 0
      }
    }
    
    
    file_list[ii] <- cope3
    ii <- ii + 1
    
  }
}
final_design <-
  cbind(design_ctl, design_exp, design_fluo, design_restrain)
write.table(
  file_list,
  paste(export_folder, '/', 'scan_list.txt', sep = ''),
  row.names = F,
  col.names = F,
  quote = F
)



K <- rbind(c(-1, 1, 0, 0),
           c(1, -1, 0, 0))


con_name <- c("exp>ctl", "exp<ctl")

desgin_mat <- c()
ii <- 1
desgin_mat[ii] <- paste('/NumWaves',  dim(final_design)[2], sep = "\t")
ii <- ii + 1
desgin_mat[ii] <-
  paste('/NumPoints',  dim(final_design)[1], sep = "\t")
ii <- ii + 1
desgin_mat[ii] <-
  paste('/PPheights\t', paste(rep('1', dim(final_design)[1]), collapse = '\t'), sep =
          "\t")
ii <- ii + 1
desgin_mat[ii] <- ''
ii <- ii + 1
desgin_mat[ii] <- paste('/Matrix', sep = "\t")
for (i in 1:dim(final_design)[1]) {
  ii <- ii + 1
  desgin_mat[ii] <-
    paste(as.character(final_design[i, ]), collapse = '\t')
  
}
write.table(
  desgin_mat,
  file.path(export_folder, 'design.mat'),
  quote = F,
  row.names = F,
  col.names = F
)


desgin_con <- c()
ii <- 1
for (i in 1:dim(K)[1]) {
  desgin_con[ii] <-
    paste(paste('/ContrastName', i, sep = ''),  con_name[i], sep = "\t")
  ii <- ii + 1
}
desgin_con[ii] <- paste('/NumWaves',  dim(K)[2], sep = "\t")
ii <- ii + 1
desgin_con[ii] <- paste('/NumContrasts',  dim(K)[1], sep = "\t")
ii <- ii + 1
desgin_con[ii] <-
  paste('/PPheights\t', paste(rep('1', dim(K)[1]), collapse = '\t'), sep =
          "\t")
ii <- ii + 1
desgin_con[ii] <-
  paste('/RequiredEffect\t', paste(rep('1', dim(K)[1]), collapse = '\t'), sep =
          "\t")
ii <- ii + 1
desgin_con[ii] <- ''
ii <- ii + 1
desgin_con[ii] <- paste('/Matrix', sep = "\t")
for (i in 1:dim(K)[1]) {
  ii <- ii + 1
  desgin_con[ii] <- paste(as.character(K[i, ]), collapse = '\t')
  
}
write.table(
  desgin_con,
  file.path(export_folder, 'design.con'),
  quote = F,
  row.names = F,
  col.names = F
)


ii <- 1
script.randomize <- c()
script.randomize[ii] <- paste('#!/bin/bash', sep = "\t")
ii <- ii + 1
script.randomize[ii] <- paste('cd', export_folder, sep = ' ')
ii <- ii + 1
script.randomize[ii] <-
  paste('fslmerge -t merge $(cat scan_list.txt)', sep = "\t")
ii <- ii + 1
script.randomize[ii] <-
  paste(
    'fsl_glm -i merge.nii.gz --out_z=glm_z --out_p=glm_p -d design.mat -c design.con -m',
    template_mask,
    sep = " "
  )
ii <- ii + 1
script.randomize[ii] <- paste('fslsplit glm_z', sep = " ")
ii <- ii + 1
script.randomize[ii] <-
  paste('easythresh vol0000.nii.gz',
        template_mask,
        '2.3 0.05',
        template ,
        'thr_pos',
        sep = " ")
ii <- ii + 1
script.randomize[ii] <-
  paste('easythresh vol0001.nii.gz',
        template_mask,
        '2.3 0.05',
        template ,
        'thr_neg',
        sep = " ")
ii <- ii + 1
script.randomize[ii] <-
  paste(
    'fslmaths thresh_thr_neg.nii.gz -mul -1 -add thresh_thr_pos.nii.gz thresh_thr.nii.gz',
    sep = " "
  )


write.table(
  script.randomize,
  file.path(export_folder, 'randomize.sh'),
  quote = F,
  row.names = F,
  col.names = F
)
```
run the group level GLM, and extract z-statistics from the statistical maps
```{bash run opto dorsal raphe glm, eval=FALSE}
cd /project/4180000.18/optoDR/statistics
bash randomize.sh
fslmeants -i vol0000.nii.gz -m $template_mask -o $init_folder'/assets/opto/optoDR_z' --label=$atlas
```



### plotting optogenetics response as a function of viral tracer maps. 
### Figure S7B and C
```{r plotting opto, eval=FALSE}
library(ggplot2)
library(gridExtra)
library(parameters)


atlas.roi <- read.csv2(atlas.label)

atlas.roi$aAI.z <-
  as.numeric(unlist(read.table('assets/img/SBA_aAI_z')[1, ]))
atlas.roi$ACA.z <-
  as.numeric(unlist(read.table('assets/img/SBA_ACA_z')[1, ]))
atlas.roi$MOp.z <-
  as.numeric(unlist(read.table('assets/img/SBA_MOp_z')[1, ]))

atlas.roi$aAI.tracer.z <-
  as.numeric(unlist(read.table('assets/img/SBA_tracer_aAI_z')[1, ]))
atlas.roi$ACA.tracer.z <-
  as.numeric(unlist(read.table('assets/img/SBA_tracer_ACA_z')[1, ]))
atlas.roi$MOp.tracer.z <-
  as.numeric(unlist(read.table('assets/img/SBA_tracer_MOp_z')[1, ]))

atlas.roi$opto.AI.z <-
  as.numeric(unlist(read.table('assets/opto/opto_z')[1, ]))
atlas.roi$opto.DR.z <-
  as.numeric(unlist(read.table('assets/opto/optoDR_z')[1, ]))

#remove 0's from the FC rows, due to incomplete FOV coverage in the FC but not SC experiment.
atlas.roi$aAI.z[which(atlas.roi$aAI.z == 0)] <- NaN
atlas.roi <- atlas.roi[complete.cases(atlas.roi), ]

atlas.roi$hex <- paste('#', as.character(atlas.roi$hex), sep = '')
atlas.roi$hex[atlas.roi$hex == '#19399'] <- '#019399'
atlas.roi <- atlas.roi[-c(which(atlas.roi$hex == '#0')), ]

write.csv2(atlas.roi, 'assets/table/zstat_roi.csv', row.names = FALSE)

rsn.list <- factor(c('SN', 'DMN', 'LCN'), levels = c('DMN', 'SN', 'LCN'))

mod.full <-
  lm(opto.DR.z ~ aAI.tracer.z + ACA.tracer.z + MOp.tracer.z, atlas.roi)
s <- parameters(mod.full, standardize = 'refit')

capture.output(s, file = "assets/img/optoDR2tracer.txt")

p1<- qplot(opto.DR.z,ACA.tracer.z,data=atlas.roi) +   geom_smooth(method = "lm", formula = y ~ x, color='#e41a1c',fill='#e41a1c') + theme_classic() + labs(y ='ACA Tracer [z-score]', x = 'CBV response [z-score]')+theme(text = element_text(size=12))
p2<- qplot(opto.DR.z,aAI.tracer.z,data=atlas.roi) +   geom_smooth(method = "lm", formula = y ~ x, color='#4daf4a',fill='#4daf4a') + theme_classic() + labs(y ='aAI Tracer [z-score]', x = 'CBV response [z-score]')+theme(text = element_text(size=12))
p3<- qplot(opto.DR.z,MOp.tracer.z,data=atlas.roi) +   geom_smooth(method = "lm", formula = y ~ x, color='#377eb8',fill='#377eb8') + theme_classic() + labs(y ='MOp Tracer [z-score]', x = 'CBV response [z-score]')+theme(text = element_text(size=12))

p4<-grid.arrange(p1, p2, p3, nrow = 2)


ggsave(
  'assets/img/optoDR2tracer.png',
  plot = p4,
  device = 'png',
  units = 'mm',
  dpi = 300,
  width = 166,
  height = 120
)
ggsave(
  'assets/img/optoDR2tracer.svg',
  plot = p4,
  device = 'svg',
  units = 'mm',
  dpi = 300,
  width = 166,
  height = 120
)


mod.full <-
  lm(opto.AI.z ~ aAI.tracer.z + ACA.tracer.z + MOp.tracer.z, atlas.roi)
s <- parameters(mod.full, standardize = 'refit')

capture.output(s, file = "assets/img/optoAI2tracer.txt")

p1<- qplot(opto.AI.z,ACA.tracer.z,data=atlas.roi) +   geom_smooth(method = "lm", formula = y ~ x, color='#e41a1c',fill='#e41a1c') + theme_classic() + labs(y ='ACA Tracer [z-score]', x = 'BOLD response [z-score]')+theme(text = element_text(size=12))
p2<- qplot(opto.AI.z,aAI.tracer.z,data=atlas.roi) +   geom_smooth(method = "lm", formula = y ~ x, color='#4daf4a',fill='#4daf4a') + theme_classic() + labs(y ='aAI Tracer [z-score]', x = 'BOLD response [z-score]')+theme(text = element_text(size=12))
p3<- qplot(opto.AI.z,MOp.tracer.z,data=atlas.roi) +   geom_smooth(method = "lm", formula = y ~ x, color='#377eb8',fill='#377eb8') + theme_classic() + labs(y ='MOp Tracer [z-score]', x = 'BOLD response [z-score]')+theme(text = element_text(size=12))

p4<-grid.arrange(p1, p2, p3, nrow = 2)


ggsave(
  'assets/img/optoAI2tracer.png',
  plot = p4,
  device = 'png',
  units = 'mm',
  dpi = 300,
  width = 166,
  height = 120
)
ggsave(
  'assets/img/optoAI2tracer.svg',
  plot = p4,
  device = 'svg',
  units = 'mm',
  dpi = 300,
  width = 166,
  height = 120
)


```

![Figure S7b](assets/img/optoDR2tracer.svg)

![Figure S7b](assets/img/optoAI2tracer.svg)
