Skip to content

Commit

Permalink
fill repo
Browse files Browse the repository at this point in the history
  • Loading branch information
viko18 committed Sep 27, 2022
1 parent dd53f01 commit 172c2ca
Show file tree
Hide file tree
Showing 31 changed files with 2,749 additions and 0 deletions.
71 changes: 71 additions & 0 deletions 01_Preprocessing/GM_matrix_create.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,71 @@
#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)
# Create masked GM matrix (subj x voxels) for OPNMF analyses #

### Arguments ###

# 1: Species investigated - "chimp" or "human"
# 2: Name of meta data .csv file #
# 3: Name and relative path of GM mask Nifti image #


wd <- "~/projects/chimp_human_opnmf/"
setwd(wd)

# load in libraries needed for script #
# need to install pacman if not already installed
if (!require("pacman")) install.packages("pacman")
pacman::p_load(fs, readr, dplyr, neurobase, oro.nifti)

# sample species either human (IXI) or chimp for locating images and naming outputs #
cohort <- args[1] # "chimp" or "human

# sample meta data to gather image names #
meta_name <- args[2] #"Chimp_meta_QC_u50_n189.csv"
meta_data <- read_csv(file = paste(path_wd(), "data", cohort, meta_name,
sep = "/"))
# smoothing kernel of images - chimp 4mm & IXI 6mm #
if (cohort == "chimp") {
s_kernel <- "s4"
} else if (cohort == "human") {
s_kernel <- "s6"
} else {
print("!! ERROR don't recognise species given should be chimp or human !!")
}
# file names of images to create GM data matrix #
subj_name <- c(paste0(s_kernel, "_mwp1", meta_data$Subject, ".gz"))
# path to images#
img_path <- paste(path_wd(), "data", cohort, "mwp1", sep = "/")

# create GM index at TPM 0.3GM #
mask_name <- args[3] #"masks/Chimp_GM_TPM_03_Cortex_2mm.nii.gz"
GM_mask <- readnii(paste(path_wd(), mask_name ,sep = "/"))
# GM index mask vector #
GM_ind <- c(img_data(GM_mask > 0))

# empty lists #
imgs <- vector(mode = "list", length = nrow(meta_data))
imgs_GM <- vector(mode = "list", length = nrow(meta_data))

# fill lists with GM data from each subject #
for (i in 1:length(subj_name)) {
# read nifti images #
imgs[[i]] <- readnii(paste(img_path, subj_name[i], sep= "/"))
# mask GM data #
imgs_GM[[i]] <- c(img_data(imgs[[i]])[GM_ind])
}

# create a GM matrix subj x voxels #
GM_mat <- rbind(matrix(unlist(imgs_GM),
ncol = length(imgs_GM[[1]]), byrow = TRUE))
# Name of outpu data matrix #
if (cohort == "chimp") {
out_name <- paste0("Chimp_cort_n", nrow(meta_data), "_TPM03_mat.rds")
} else if (cohort == "human") {
out_name <- paste0("IXI_cort_n", nrow(meta_data), "_mat.rds")
} else {
print("!! ERROR don't recognise species given should be chimp or human !!")
}
# save data matrix #
saveRDS(GM_mat, file = paste(path_wd(), "data", cohort, out_name, sep="/"))

6 changes: 6 additions & 0 deletions 01_Preprocessing/README.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
# Structural Preprocessing

Structural preprocessing (T1-weighted images) for chimpanzee and human cohorts utilizing [SPM12](https://www.fil.ion.ucl.ac.uk/spm/software/spm12/) and [CAT12](http://www.neuro.uni-jena.de/cat/). The Matlab scripts are to be used as batch inputs to then adapted in the SPM batch editor.

Following preprocessing move the ```mwp1*``` and smoothed modulated images into ```$project_dir/data/$cohort/wmp1/``` directory and then run the ```run_GM_matrix_create.sh``` to create the data matrix for OPNMF parcellation.

162 changes: 162 additions & 0 deletions 01_Preprocessing/chimp_preproc_batch.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,162 @@
clc, clear

% ------------------------------------------------------------------------
% Chimpanzee batch options that can be changed in the GUI or here directly
% ------------------------------------------------------------------------
% open the chimpanzee preprocessing from cat12
cat12('chimpanzee');
% set the number of cores for parallel processing
try
numcores = max(cat_get_defaults('extopts.nproc'),1);
catch
numcores = 0;
end
%% CAT preprocessing
% ------------------------------------------------------------------------
% Dir containing chimp original T1w images %
chimps_dir = ...
'C:\Users\svickery\Documents\preproc_test';
% chimp images %
chimp_img = dir(fullfile(chimps_dir,'*.nii'));
% IMG need to be in a cell to be read by batch %
for subj = 1:numel(chimp_img)
matlabbatch{1}.spm.tools.cat.estwrite.data{subj,1} = ...
fullfile(chimps_dir, chimp_img(subj).name);
end
% Check name and location of template dir in SPM dir %
% This should be the dir with cat v1725
cattempdir = fullfile(spm('dir'),'toolbox','cat12','animal_templates');
% CAT preprocessing expert options
% SPM parameter
matlabbatch{1}.spm.tools.cat.estwrite.data_wmh = {''};
matlabbatch{1}.spm.tools.cat.estwrite.nproc = numcores;
matlabbatch{1}.spm.tools.cat.estwrite.opts.tpm = {fullfile(cattempdir,'chimpanzee_TPM.nii')}; % Juna chimp TPM
matlabbatch{1}.spm.tools.cat.estwrite.opts.affreg = 'none';
matlabbatch{1}.spm.tools.cat.estwrite.opts.ngaus = [1 1 2 3 4 2];
matlabbatch{1}.spm.tools.cat.estwrite.opts.warpreg = [0 0.001 0.5 0.05 0.2];
matlabbatch{1}.spm.tools.cat.estwrite.opts.bias.spm.biasfwhm = 30; % small values are important to remove the bias but 30 mm is more or less the limit
matlabbatch{1}.spm.tools.cat.estwrite.opts.bias.spm.biasreg = 0.001; %
matlabbatch{1}.spm.tools.cat.estwrite.opts.acc.spm.samp = 1.5; % ############# higher resolutions helps but takes much more time (e.g., 1.5 about 4 hours), so 1.0 to 1.5 mm seems to be adequate
matlabbatch{1}.spm.tools.cat.estwrite.opts.acc.spm.tol = 1e-06; % smaller values are better and important to remove the bias in some image
matlabbatch{1}.spm.tools.cat.estwrite.opts.redspmres = 0;

% segmentation options
matlabbatch{1}.spm.tools.cat.estwrite.extopts.segmentation.APP = 1070;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.segmentation.NCstr = -Inf;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.segmentation.spm_kamap = 0;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.segmentation.LASstr = 1.0;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.segmentation.gcutstr = 2;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.segmentation.cleanupstr = 0.5;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.segmentation.BVCstr = 0.5;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.segmentation.WMHC = 0;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.segmentation.SLC = 0;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.segmentation.mrf = 1;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.segmentation.restypes.best = [0.5 0.3];

% registration options
matlabbatch{1}.spm.tools.cat.estwrite.extopts.registration.T1 = {fullfile(cattempdir,'chimpanzee_T1.nii')}; % Juna chimp T1
matlabbatch{1}.spm.tools.cat.estwrite.extopts.registration.brainmask = {fullfile(cattempdir,'chimpanzee_brainmask.nii')}; % Juna chimp brainmask
matlabbatch{1}.spm.tools.cat.estwrite.extopts.registration.cat12atlas = {fullfile(cattempdir,'chimpanzee_cat.nii')}; % Juna chimp cat atlas
matlabbatch{1}.spm.tools.cat.estwrite.extopts.registration.darteltpm = {fullfile(cattempdir,'chimpanzee_Template_1.nii')}; % there is no Juna chimp Dartel template as shooting is much better
matlabbatch{1}.spm.tools.cat.estwrite.extopts.registration.shootingtpm = {fullfile(cattempdir,'chimpanzee_Template_0_GS.nii')}; % Juna chimp shooting template
matlabbatch{1}.spm.tools.cat.estwrite.extopts.registration.regstr = 0.5;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.registration.vox = 2; % Downsample to 2mm for faster OPNMF processing

% surface options
matlabbatch{1}.spm.tools.cat.estwrite.extopts.surface.pbtres = 0.5;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.surface.pbtmethod = 'pbt2x';
matlabbatch{1}.spm.tools.cat.estwrite.extopts.surface.pbtlas = 0;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.surface.collcorr = 0;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.surface.reduce_mesh = 1;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.surface.vdist = 1.33333333333333;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.surface.scale_cortex = 0.7;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.surface.add_parahipp = 0.1;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.surface.close_parahipp = 0;

% admin options
matlabbatch{1}.spm.tools.cat.estwrite.extopts.admin.experimental = 0;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.admin.new_release = 0;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.admin.lazy = 0; % ############## avoid reprocessing
matlabbatch{1}.spm.tools.cat.estwrite.extopts.admin.ignoreErrors = 1;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.admin.verb = 2;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.admin.print = 2;

% output options
matlabbatch{1}.spm.tools.cat.estwrite.output.surface = 0; % surface reconstruction - not yet optimised for non-human primates
matlabbatch{1}.spm.tools.cat.estwrite.output.surf_measures = 3;

% volume atlas maps
matlabbatch{1}.spm.tools.cat.estwrite.output.ROImenu.atlases.chimpanzee_atlas_davi = 1;
matlabbatch{1}.spm.tools.cat.estwrite.output.ROImenu.atlases.ownatlas = {''}; % you can add own atlas maps but they have to be in the same orientation as the other template files especially the final GS template

% surface atlas maps
matlabbatch{1}.spm.tools.cat.estwrite.output.sROImenu.satlases.Desikan = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.sROImenu.satlases.Destrieux = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.sROImenu.satlases.HCP = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.sROImenu.satlases.Schaefer2018_100P_17N = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.sROImenu.satlases.Schaefer2018_200P_17N = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.sROImenu.satlases.Schaefer2018_400P_17N = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.sROImenu.satlases.Schaefer2018_600P_17N = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.sROImenu.satlases.ownatlas = {''};

% volume output
matlabbatch{1}.spm.tools.cat.estwrite.output.GM.native = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.GM.warped = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.GM.mod = 1; % needed for VBM
matlabbatch{1}.spm.tools.cat.estwrite.output.GM.dartel = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.WM.native = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.WM.warped = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.WM.mod = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.WM.dartel = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.CSF.native = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.CSF.warped = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.CSF.mod = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.CSF.dartel = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.bias.native = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.bias.warped = 1;
matlabbatch{1}.spm.tools.cat.estwrite.output.bias.dartel = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.jacobianwarped = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.warps = [0 0];

% special maps
matlabbatch{1}.spm.tools.cat.estwrite.output.ct.native = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.ct.warped = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.ct.dartel = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.pp.native = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.pp.warped = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.pp.dartel = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.WMH.native = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.WMH.warped = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.WMH.mod = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.WMH.dartel = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.SL.native = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.SL.warped = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.SL.mod = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.SL.dartel = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.TPMC.native = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.TPMC.warped = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.TPMC.mod = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.TPMC.dartel = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.atlas.native = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.atlas.warped = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.atlas.dartel = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.label.native = 1;
matlabbatch{1}.spm.tools.cat.estwrite.output.label.warped = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.label.dartel = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.las.native = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.las.warped = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.las.dartel = 0;

% ------------------------------------------------------------------------

%% smoothing
% ------------------------------------------------------------------------
matlabbatch{2}.spm.spatial.smooth.fwhm = repmat(4,1,3); % smoothing filter size 4x4x4
matlabbatch{2}.spm.spatial.smooth.dtype = 0;
matlabbatch{2}.spm.spatial.smooth.im = 0;
matlabbatch{2}.spm.spatial.smooth.prefix = 's4_';
% ------------------------------------------------------------------------

%%
% Use to directly run script keep commented out if want to adapt using GUI
%spm_jobman('run',matlabbatch);
Loading

0 comments on commit 172c2ca

Please sign in to comment.