fMRI gICA & DR Analysis

Before preprocessing your scans please read the RABIES documentation (https://github.com/CoBrALab/RABIES_tutorial)

Disclaimer: This wiki was created my Marina Broomfield. Please cite this page and give credit when using or copying the following code.

Group ICA - Classifying Networks

Once scans have been preprocessed following the RABIES documentation, one can run gICA. gICA is a data driven approach which allows you to obtain patters of brain activity simultaneously firing together. Usually, gICA is initially run with 30 independent components (ICs), however, this can change based on the number of networks observed and whether these networks are broken down into different ICs. To validate whether each IC is a reliable network or an artifact, you can visualize the ICs by opening the .nii.gz file created by RABIES. Once the nii.gz file is open, to go over each IC, press the keys "1" and "2".

module load minc-toolkit-v2
Display /analysis_outputs/analysis_datasink/group_ICA_dir/group_melodic.ica/melodic_IC.nii.gz 

Make sure to also go over the .png images of each IC component. You can find these in:

/analysis_outputs/analysis_datasink/group_ICA_dir/group_melodic.ica/report/

After classifying each IC into a network or an artifact, make sure you correctly lable which network the IC belongs to. You can find reference networks in RABIES and the following paper (https://pubmed.ncbi.nlm.nih.gov/26296501/).

Classifying Outliers

After labeling each network and running dual regression in RABIES, it is important to make sure your data is clean. to do this follow the RABIES documentation and go over the following png images created:

#For each network:
/diagnosis_outputs/data_diagnosis_datasink/analysis_QC/analysis_QC/non_parametric_stats/ 
/diagnosis_outputs/data_diagnosis_datasink/analysis_QC/analysis_QC/parametric_stats/
#For each subject:
/diagnosis_outputs/data_diagnosis_datasink/figure_spatial_diagnosis/
/diagnosis_outputs/data_diagnosis_datasink/figure_temporal_diagnosis/

You should also look at the DICE values for each subject across each network. Usually a threshold of 0.4 should be placed for the somatomotor/somatosensory networks. This threshold however might be susceptible to change depending on the network being evaluated. The following code runs in R and allows you to visualize each subject's DICE for a specific network:

# Load the required libraries
library(ggplot2)
library(dplyr)
library(tidyr)
library(car)

# Read the CSV file
data <- read.csv(file = "/data/scratch2/marina/experiment1_oct_2023/scanning/functional/stats/outlier_detection.csv", head = TRUE, sep = ",")
data$sub_cohort <- factor(data$sub_cohort)
data$cohort <- factor(data$cohort)

#Replace 'subject_ID', 'diceDR0', and 'sex' with the actual column names in your data frame
DR0 <- ggplot(data, aes(x = scan_ID, y = diceDR0, color = sub_cohort)) +
  geom_point(size = 3) +
  geom_vline(aes(xintercept = as.numeric(scan_ID)), linetype = "dashed", alpha = 0.5) +
  labs(title = "Scatter plot of somatosensory by subject",
       x = "Subject ID",
       y = "diceDR1",
       color = "Sex")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),  # Adjust the size parameter for x-axis labels
        axis.text.y = element_text(hjust = 1, size = 12),
        plot.title = element_text(hjust = 0.5, size = 16))  # Adjust the size and centering for the title

# Print the plot
print(DR0)

Isolating each Network by Subject

You will first have to isolate each of the networks for each subject. To do this you can run the following script in python:

# Here we are importing all the python packages that will be used in the notebook
import numpy as np
import pandas as pd
import nibabel as nb
import nilearn
import matplotlib.pyplot as plt
import glob
import os

# Choose where to save output nifti images
output_dir='/data/scratch2/marina/experiment1_oct_2023/scanning/functional/isolated_IC_components/'
# Select which ICs are brain networks
network_idx=np.array([0,1,2,8])

# Path to where your dual regression outputs are saved
file_list=glob.glob("/data/scratch2/marina/experiment1_oct_2023/scanning/functional/diagnosis_outputs/analysis_datasink/dual_regression_nii/*/*.nii.gz")

# Loop through each subject
for j in file_list: 
        print(j)
        dual_reg_file=nb.load(j) # Load dual regression file
        sub=os.path.basename(j)[:-7] # Name of your nifti file (subject,session)
        im_ar=np.asarray(dual_reg_file.dataobj) # Create matrix (array) for the dual regression file

        # Loop through all ICs
        for i in network_idx:  
                print(i)
                ic = im_ar[:,:,:,i] 
                img= nb.Nifti1Image(ic, dual_reg_file.affine, dual_reg_file.header) # Create nifit image for a single IC
                nb.save(img,f'{output_dir}/{sub}_IC{i}.nii.gz') # Save nifti image

Comparing Network Amplitude Across the Brain

After you obtain each subject's network transform the nifi files to .mnc files. Create an input matrix .csv file with the following information: subject_ID, group, Network1, Netwrok2 ..etc. Under each Network, put the path to the .mnc files for each subject.

Once the input.csv has been created you can visualize and compare each network across the brain, before comparing functional MRI outputs, make sure you choose the appropriate statistical model to run you analysis, this can be done by running AIC. Once a model has been chosen, spatial maps comparing network amplitude across the brain can be visualized with the following code:

#Load necessary libraries
library(ggplot2) # Library for plotting
library(RMINC) # Library to work with images
library(lme4) # Library for running linear models
library(tidyverse)


# Read csv with input data. Columns are variables. Rows are subjects.
data<-read.csv(file="/data/scratch2/marina/experiment1_oct_2023/scanning/functional/stats/spatialmaps_sex/input.csv",head=TRUE,sep=",")
# Create variable with a string that contains the path to the mask to use for the analysis
mask = "/data/scratch2/marina/experiment1_oct_2023/scanning/functional/stats/spatialmaps_sex/spatialMaps_masksub-C01S015_ses-1_task-rest_acq-EPI_bold_RAS_EPI_brain_mask.mnc"
# Order sex
data$sex <- factor(data$sex, ordered = FALSE, levels=c("M","F"))
# Set the control group as reference. This is done for interpretability of the model
data$sex = relevel(data$sex, ref = "M")

# Now you can run the linear model using the mincLm function
SMM_data <- mincLm(SMM ~ sex + Mean_FD, data, mask = mask)

# Create a minc file with the t-values for the specific predictor in the model 'cheese'- this will allow us to visualize the t-values on top of out brain template later on
mincWriteVolume(SMM_data, '/data/scratch2/marina/experiment1_oct_2023/scanning/functional/stats/spatialmaps_sex/figures/lmSMM_SexFD.mnc',column = 'tvalue-sexF')
# Save mincLm model "cheese" for future reference
#save.image("/data/scratch2/marina/experiment1_oct_2023/scanning/functional/stats/spatialmaps_sex.RData")

SMM_data = mincFDR(SMM_data, mask= mask)
SMM_data

mincWriteVolume(FDR, '/data/scratch2/marina/experiment1_oct_2023/scanning/functional/stats/spatialmaps_sex/figures/lm_SMM_sexFD_FDR.mnc',column = 'tvalue-sexF')


# Load your population template in which the t-stats will be overlapped. It needs to be the same resolution as the t-stats minc volume created above
functemp <- mincArray(mincGetVolume("/data/scratch2/marina/experiment1_oct_2023/scanning/functional/stats/spatialmaps_sex/template_resampled.mnc"))


#PLOTTING STATS MAP TO TEMPLATE

# Limits for the most posterior to most anterior slices that we want to plot
ylow<-10
yhigh<-30
# Plot slices with the t-values overlapped
mincPlotSliceSeries(functemp, dimension=2,mincArray(SMM_data,'tvalue-sexF')
                    ,plottitle='SMM Amplitude Across Sex -FD Regressed'
                    ,mfrow = c(4,4),                   # layout
                    ,low=2.02,high =3.55,symmetric=T
                    ,col =  heat.colors(255)
                    ,anatCol = gray.colors(255, start=0.0)
                    ,begin=ylow,end=yhigh
                    ,legend="t-statistics")
save.image("/data/scratch2/marina/experiment1_oct_2023/scanning/functional/stats/spatialmaps_sex/figures/lmSMMspatialmap_SexFD.png")