---
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 preprocess the data using the Co-activation pattern method, or [CAP](https://doi.org/10.1016/j.neuroimage.2020.116621). We use a combination of human, macaque, and mouse resting-state fMRI data. 

### Inputs
awake and CPS mouse datasets scans 
10 preprocessed macaque fMRI scans. 
15 randomly selected HCP 7T rsfMRI scans. (see 7T scan selection chunk)
seeds (see [asset](assest.md))

### Outputs
**Figure 1A, S1, S2**   
Anesthetized mouse CAPs (asset/CAP/mouse/CAP_NIFTI_*.nii)   
Awake mouse CAPs (asset/CAP/awake/CAP_NIFTI_*.nii)    
Anesthetized macaque CAPs (assets/CAP/monkey/CAP_NIFTI_*.nii)   
Awake human CAPs (assets/CAP/human/CAP_NIFTI_*.nii)   

**Figure 1BC, S3**   
3D reconstructed synthetics CAP 3 (assets/gene_translation/weighted_dk_atlas_CAP3.nii)   
Atlas-extracted synthetics CAPs (assets/gene_translation/weighted_array.csv)   
Atlas-extracted human CAPs (assets/gene_translation/human_array.csv)   
CAP correlation (assets/gene_translation/Dabest_match.csv)

**Figure 1DE**   
Anesthetized mouse CAPs transition (assets/CAP/mouse_transition.csv)   
Anesthetized macaque CAPs transition (assets/CAP/monkey_transition.csv)   
Awake human CAPs transition (assets/CAP/human_transition.csv)   

**Figure 1F**   
Anesthetized mouse CAPs entries (assets/CAP/mouse/CSD_CAP_20200511.csv)  

**Figure S3**   
https://neuroquery.org/query?text=depression


## 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'), ]
```

## CAPs
### CAPS in human
The chunk below prepare 7T S1200 data prior to co-activation pattern analysis with the [CAP toolbox](https://c4science.ch/source/CAP_Toolbox.git)

### Figure 1A
```{bash 7T scan selection, eval=FALSE}
source bash_env.sh
module unload fsl
module load fsl/5.0.9 #newer versions of FSL don't allow file change to non-GZ NIFTI
template='/opt/fsl/6.0.0/data/standard/MNI152_T1_2mm.nii.gz'
template_mask='/opt/fsl/6.0.0/data/standard/MNI152_T1_2mm_brain_mask.nii.gz'
root_dir='/project/4180000.18/Human'
data_export=$root_dir'/CAP_data'
data_dir='/project_freenas/3022017.01/S1200/'

mkdir -p $data_export

cd $root_dir

counter=0
scan_name='rfMRI_REST1_7T_PA'

#loop through the processed scans, for each scans, check if a 7T scan is present. Select the first 15 due to 64Gb memory restrictions on high performance cluster matlab node  
ls $data_dir | while read subject
do
echo "now doing scan "$subject
echo "counters "$counter

file_path=$data_dir'/'$subject'/MNINonLinear/Results/'$scan_name'/'$scan_name'_hp2000_clean.nii.gz'
file_path_mc=$data_dir'/'$subject'/MNINonLinear/Results/'$scan_name'/'$scan_name'/Movement_Regressors.txt'


if [[ "$counter" == '15' ]]; then
    break
fi
if [ -f "$file_path" ]; then
counter=$((counter+1))


qscript='fslsplit_job_'$counter'.sh'

output_sba=${data_export}/${subject}

  echo "cd "$root_dir > $root_dir/script/$qscript
  #echo "module unload fsl" >> $root_dir/script/$qscript
  #echo "module load fsl/5.0.9" >> $root_dir/script/$qscript
  echo "mkdir -p " $output_sba >> $root_dir/script/$qscript
  echo "3dresample -inset "$file_path" -prefix " $output_sba'/tmp.nii.gz -rmode Cu -master '$template >> $root_dir/script/$qscript
  echo "3dBlurInMask -input "$output_sba'/tmp.nii.gz -prefix ' $output_sba'/tmp_blur.nii.gz -FWHM 5 -mask '$template_mask >> $root_dir/script/$qscript
  echo "3dresample -inset "$output_sba'/tmp_blur.nii.gz -prefix ' $output_sba'/tmp_rs.nii.gz -rmode Cu -dxyz 3 3 3' >> $root_dir/script/$qscript
  echo "rm "$output_sba'/tmp.nii.gz' >> $root_dir/script/$qscript
  echo "rm "$output_sba'/tmp_blur.nii.gz' >> $root_dir/script/$qscript
  echo "fslsplit "$output_sba'/tmp_rs.nii.gz' $output_sba'/fmri' >> $root_dir/script/$qscript
  echo "rm "$output_sba'/tmp_rs.nii.gz' >> $root_dir/script/$qscript
  echo "ls "$output_sba"/fmri* | while read line" >> $root_dir/script/$qscript
  echo "do" >> $root_dir/script/$qscript
  echo 'fslchfiletype NIFTI ${line}' >> $root_dir/script/$qscript
  echo "done" >> $root_dir/script/$qscript
  echo "cp "$file_path $output_sba'/mc.txt' >> $root_dir/script/$qscript


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


 
fi

done
rm script/fslsplit_job_*
```

### Mouse <-> Human translation
mouse CAPs were transformed into human synthetic brain maps using [this software](https://gitlab.socsci.ru.nl/preclinical-neuroimaging/homologous-brain-expression-map.git).

### Figure 1C, Figure S4
```{r plotting CAP translation, eval=FALSE}



true.CAP <-
  t(read.table(
    'assets/gene_translation/human_array.csv',
    sep = ',',
    header = F
  ))


synthetic.CAP <-
  t(read.table(
    'assets/gene_translation/weighted_array.csv',
    sep = ',',
    header = F
  ))

merge.CAP <- cbind(true.CAP, synthetic.CAP)

merge.CAP <- merge.CAP[-c(3, 28), ]

merge.CAP.cor <- cor(merge.CAP)


merge.CAP.sub <- merge.CAP.cor[1:6, 7:12]



CAP_export_dabest <-
  data.frame(merge.CAP.sub[which(upper.tri(merge.CAP.sub, diag = FALSE)) ||
                             lower.tri(merge.CAP.sub, diag = FALSE)], c(diag(merge.CAP.sub), rep(NaN, 30)))
names(CAP_export_dabest) <- c('mismatched', 'matched')

write.csv(CAP_export_dabest,
          'assets/gene_translation/Dabest_match.csv',
          row.names = FALSE)

merge.CAP.cor.sub.melt <- melt(merge.CAP.cor[1:6, 7:12])
merge.CAP.cor.sub.melt$value[merge.CAP.cor.sub.melt$value < -0.3] <-
  -0.3
merge.CAP.cor.sub.melt$value[merge.CAP.cor.sub.melt$value > 0.3] <-
  0.3

p <- ggplot(merge.CAP.cor.sub.melt, aes(x = Var2, y = Var1, fill = value))
p + geom_tile() +
  #scale_fill_gradient(low = "blue", high = "yellow", na.value = NA)+
  theme_classic() + scale_fill_gradient2(low = "blue", mid = "white", high = "red")
ggsave(
  'assets/gene_translation/true_synthetic_cor.svg',
  device = 'svg',
  width = 70,
  height = 70,
  units = 'mm',
  dpi = 300
)
ggsave(
  'assets/gene_translation/true_synthetic_cor.png',
  device = 'png',
  width = 70,
  height = 70,
  units = 'mm',
  dpi = 300
)


merge.CAP.cor.true.melt <- melt(merge.CAP.cor[1:6, 1:6])
merge.CAP.cor.true.melt$value[merge.CAP.cor.true.melt$value < -0.7] <-
  -0.7
merge.CAP.cor.true.melt$value[merge.CAP.cor.true.melt$value > 0.7] <-
  0.7

p <- ggplot(merge.CAP.cor.true.melt, aes(x = Var2, y = Var1, fill = value))
p + geom_tile() +
  #scale_fill_gradient(low = "blue", high = "yellow", na.value = NA)+
  theme_classic() + scale_fill_gradient2(low = "blue", mid = "white", high = "red")
ggsave(
  'assets/gene_translation/true_true_cor.svg',
  device = 'svg',
  width = 70,
  height = 70,
  units = 'mm',
  dpi = 300
)
ggsave(
  'assets/gene_translation/true_true_cor.png',
  device = 'png',
  width = 70,
  height = 70,
  units = 'mm',
  dpi = 300
)


df <- data.frame(true.CAP[, 3], synthetic.CAP[, 3])
names(df) <- c('true', 'synthetic')


p <- ggplot(df, aes(x = true, y = synthetic))
p + geom_point(colour = '#CB650B') +
  geom_smooth(method = "lm",
              formula = y ~ x,
              colour = '#CB650B') + theme_classic() + labs(title = 'CAP #3', y = 'synthetic CAP', x =
                                                             'true CAP')

mod<-lm(synthetic ~ true, df)
s <- parameters(mod,standardize = 'refit',summary=TRUE)
capture.output(s, file = "assets/gene_translation/CAP3.txt")

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

```

![Figure 1C](assets/gene_translation/CAP3.svg)
![Figure S4C](assets/gene_translation/output_Dabest_match.svg)


### Compare CAP transitions across species

### Figure 1DE
```{r compare caps species, eval=FALSE}
library(parameters)

cap.human <-
  as.matrix(read.table(
    'assets/CAP/human_transition.csv',
    sep = ',',
    header = FALSE
  ))
cap.monkey <-
  as.matrix(read.table(
    'assets/CAP/monkey_transition.csv',
    sep = ',',
    header = FALSE
  ))
cap.mouse <-
  as.matrix(read.table(
    'assets/CAP/mouse_transition.csv',
    sep = ',',
    header = FALSE
  ))

#resize matrix to keep CAP transitions
cap.human <- cap.human[3:8, 3:8]
cap.monkey <- cap.monkey[3:8, 3:8]
cap.mouse <- cap.mouse[3:8, 3:8]


cap.human.order <- c(3, 1, 4, 5, 6, 2)
cap.monkey.order <- c(1, 3, 5, 4, 2, 6)
cap.mouse.order <- c(1, 4, 3, 2, 5, 6)

cap.human <- cap.human[cap.human.order, cap.human.order]
cap.monkey <- cap.monkey[cap.monkey.order, cap.monkey.order]
cap.mouse <- cap.mouse[cap.mouse.order, cap.mouse.order]

mod<-lm(c(cap.mouse) ~ c(cap.human))
parameters(mod,standardize = 'refit',summary=TRUE)
mod<-lm(c(cap.mouse) ~ c(cap.monkey))
parameters(mod,standardize = 'refit',summary=TRUE)

plot(c(cap.mouse), c(cap.human))

colnames(cap.human) <- 1:6
colnames(cap.monkey) <- 1:6
colnames(cap.mouse) <- 1:6

m.cap.mouse <- melt(cap.mouse)
p <- ggplot(m.cap.mouse, aes(x = Var2, y = Var1, fill = value))
p + geom_tile() +
  scale_fill_gradient(low = "red",
                      high = "yellow",
                      na.value = NA) +
  theme_classic()
ggsave(
  'assets/CAP/mouse_transition.svg',
  device = 'svg',
  width = 70,
  height = 70,
  units = 'mm',
  dpi = 300
)
ggsave(
  'assets/CAP/mouse_transition.png',
  device = 'png',
  width = 70,
  height = 70,
  units = 'mm',
  dpi = 300
)

m.cap.monkey <- melt(cap.monkey)
p <- ggplot(m.cap.monkey, aes(x = Var2, y = Var1, fill = value))
p + geom_tile() +
  scale_fill_gradient(low = "red",
                      high = "yellow",
                      na.value = NA) +
  theme_classic()
ggsave(
  'assets/CAP/monkey_transition.svg',
  device = 'svg',
  width = 70,
  height = 70,
  units = 'mm',
  dpi = 300
)
ggsave(
  'assets/CAP/monkey_transition.png',
  device = 'png',
  width = 70,
  height = 70,
  units = 'mm',
  dpi = 300
)


m.cap.human <- melt(cap.human)
p <- ggplot(m.cap.human, aes(x = Var2, y = Var1, fill = value))
p + geom_tile() +
  scale_fill_gradient(low = "red",
                      high = "yellow",
                      na.value = NA) +
  theme_classic()
ggsave(
  'assets/CAP/human_transition.svg',
  device = 'svg',
  width = 70,
  height = 70,
  units = 'mm',
  dpi = 300
)
ggsave(
  'assets/CAP/human_transition.png',
  device = 'png',
  width = 70,
  height = 70,
  units = 'mm',
  dpi = 300
)


m.cap.combined <-
  data.frame(m.cap.mouse$value, m.cap.monkey$value, m.cap.human$value)
names(m.cap.combined) <- c('mouse', 'monkey', 'human')

p <- ggplot(m.cap.combined, aes(x = mouse, y = monkey))
p + geom_point(colour = '#CB650B') +
  geom_smooth(method = "lm",
              formula = y ~ x,
              colour = '#CB650B') + theme_classic() + labs(title = 'CAP transition probability', y =

                                                                                                                       'Macaque', x = 'Mouse')
mod<-lm(monkey ~ mouse, m.cap.combined)
s <- parameters(mod,standardize = 'refit',summary=TRUE)
capture.output(s, file = "assets/CAP/mouse_monkey.txt")


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

p <- ggplot(m.cap.combined, aes(x = mouse, y = human))
p + geom_point(colour = '#CB650B') +
  geom_smooth(method = "lm",
              formula = y ~ x,
              colour = '#CB650B') + theme_classic() + labs(title = 'CAP transition probability', y =
                                                       'Human', x = 'Mouse')
mod <- lm(human ~ mouse, m.cap.combined)
s <- parameters(mod,standardize = 'refit',summary=TRUE)
capture.output(s, file = "assets/CAP/mouse_human.txt")

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

```
![figure 1D](assets/CAP/mouse_transition.svg)
![figure 1E left](assets/CAP/mouse_human.svg)
![figure 1E right](assets/CAP/mouse_monkey.svg)

### Comparing CAPs projected in the Chronic Social Stress dataset
This is the preparation for DABEST plots. I used the following [resource](https://www.estimationstats.com/#/) to make the plots. I thus reformatted the list of CAPs entries to make it comparable with this resource. 

### Figure 1F

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


study <- read.csv('assets/CAP/mouse/CSD_CAP_20200511.csv')

m <-
  melt(
    study,
    id.vars = c('ID', 'group'),
    measure.vars = c(
      'entrie.cap1',
      'entrie.cap2',
      'entrie.cap3',
      'entrie.cap4',
      'entrie.cap5',
      'entrie.cap6'
    )
  )
tmp <- dcast(m, ID + value ~ paste(variable, group, sep = '_'), fill = NaN)
tmp <- tmp[, -c(1, 2)]
write.csv(tmp, 'assets/CAP/mouse_entries.csv', row.names = FALSE)

```

![figure 1E right](assets/CAP/output_mouse_entries.svg)
