Running python PLS on voxelwise data and normalizing Jacobians

Running a PLS on voxelwise data in python uses, for the most part, the instructions here: https://github.com/CoBrALab/documentation/wiki/Running-PLS-in-Python:-pyls, but this wiki will add some instructions for several specific alterations.

Prepare Jacobians for vertexwise PLS

Thus far, we have used helper scripts in R to prepare the Jacobians for the PLS. The goal of the script is to read in the Jacobians in minc files, write them to a matrix restricted to the mask, and save this as a csv.

#Load necessary packages
library(RMINC)
library(lme4)
library(lmerTest)

# Read in the csv with the listed jacobians that you want to include in the analysis, as you would for a classic twolevel analysis 
# They must be mincs, in this case in a column named "file" 
braintrix = read.csv("t2_full_smooth_jacs.csv")

# Read in the minc of the mask 
mask = mincGetVolume("../../make_mask/adult_mask.mnc")

# Create a column filled with the jacobians, but restricted to the mask
braintrix$jacdata = t(sapply(braintrix$file, function(f) mincGetVolume(f)[mask > 0]))

# making the jacobian data its own object
braintrix2 = braintrix$jacdata

# Write out the matrix 
# This can then become the brain inputs to your PLS.
write.csv(braintrix2, file ="./t2_brain_matrix.csv", row.names=FALSE)

Write out bootstrapped brain variables into mincs

After you have run the PLS and want to visualize the bootstrapped brain results from a significant LV, do the following in R.

# Load necessary libraries
library(tidyverse)       # for convenient syntax
library(RMINC)           # for conducting neuroanatomy statistics
library(MRIcrotome)      # for visualisation

# Conjure the dataframe of bootstraps from a csv file. Note that you have to tell r that there's no headers on this file otherwise, it will think the first row of data is a header.
booties = read.csv("bootres_x_weights_normed.csv", header=FALSE)

# Load in the same mask object
mask = mincGetVolume("../../make_mask/adult_mask.mnc")

# Subset down to the LV that you are interested in
# In this case, I used LV2, aka second column. Change the "X" in booties[,X] to the LV you want. 
# then set equal to where mask is 1 (in the brain)
mask[mask == 1] = booties[,2]

# Create the mnc file. Note that a BSR threshold of 3.29 corresponds to p < 0.001, 2.57 to p < 0.01 and 1.95 to p < 0.05.
mincWriteVolume(mask, "./LV2_bootstrap_map.mnc")

################
# Visualize with MRIcrotome images

# Program the average and mask in:
anatVol = mincArray(mincGetVolume("../../dbm_outputs/template_sharpen_shapeupdate.mnc"))
averagemask = mincArray(mincGetVolume("../../make_mask/adult_mask.mnc"))

#Here we figure out the extent of the mask, so we can use it to limit the FOV
# You may want to adjust these slighty buy adding/subtracting, so you don't quite go to the edge of the mask
# dimension 1 is sagittal, dim 2 is coronal dim 3 is axial
dim1_begin <- min(which(averagemask == 1, arr.ind=TRUE)[,"dim1"]) + 10
dim1_end <-max(which(averagemask == 1, arr.ind=TRUE)[,"dim1"]) - 10
dim2_begin <- min(which(averagemask == 1, arr.ind=TRUE)[,"dim2"]) + 15
dim2_end <- max(which(averagemask == 1, arr.ind=TRUE)[,"dim2"]) - 15
dim3_begin <- min(which(averagemask == 1, arr.ind=TRUE)[,"dim3"]) + 15
dim3_end <- max(which(averagemask == 1, arr.ind=TRUE)[,"dim3"]) - 5

#We use a tryCatch here to keep going if there's some kind of error in a an individual plot
tryCatch({
  # Set the thresholds of the bootstraps
  lowerthreshold = 1.95
  upperthreshold = 3.29
  # make the figure
  sliceSeries(nrow = 6, ncol= 3, begin = dim2_begin, end = dim2_end, dimension = 2) %>%  # slice sequence to display stats
    anatomy(anatVol, low=1, high=5.9) %>%
    overlay(mincArray(mask),     # specify the overlay stuff
            low=lowerthreshold, high=upperthreshold, symmetric = T, alpha=0.6) %>%       # stat range; symmetric = T means show positive and negative values
    addtitle("Coronal") %>%
    sliceSeries(nrow = 6, ncol= 2, begin = dim1_begin, end = dim1_end, dimension = 1) %>%  # slice sequence to display stats
    anatomy(anatVol, low=1, high=5.9) %>%
    overlay(mincArray(mask),     # specify the overlay stuff
            low=lowerthreshold, high=upperthreshold, symmetric = T, alpha=0.6) %>%       # stat range; symmetric = T means show positive and negative values
    addtitle("Sagittal") %>%
    sliceSeries(nrow = 3, ncol= 3, begin = dim3_begin, end = dim3_end, dimension = 3) %>%  # slice sequence to display stats
    anatomy(anatVol, low=1, high=5.9) %>%
    overlay(mincArray(mask),     # specify the overlay stuff
            low=lowerthreshold, high=upperthreshold, symmetric = T, alpha=0.6) %>%       # stat range; symmetric = T means show positive and negative values
    addtitle("Axial") %>%
    legend("Bootstrap") %>%
    draw()
  #dev.off()
}, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})

Normalize Jacobians by residualizing by a variable before PLS

Finally, you may find you want to normalize your Jacobian matrix by a variable in your demographics file. For example, I want to normalize by the weight of the mouse here. To do this, we will use R functions built by Yohan Yee.

#First, until it is installed throughout the CIC, you have to install the packages the first time you want to use it:
# install packages one time only
devtools::install_github("yohanyee/normtools")

# Load necessary libraries
library(RMINC)
library(lme4)
library(lmerTest)
library(normtools)

The first steps are identical to preparing the Jacobians, identified above:

# Conjure the dataframe in question from a csv file
braintrix = read.csv("t2_full_smooth_jacs.csv")
dems <- read.csv(file = "/data/chamal/projects/lani/adult_thc/dems_with_weight_resids.csv") #, header = TRUE)

# Create an object that embodies the mask file I created
mask = mincGetVolume("../../make_mask/adult_mask.mnc")

# make a new column in my data that only contains the masked jacobians
braintrix$jacdata = t(sapply(braintrix$file, function(f) mincGetVolume(f)[mask > 0]))

# making the jacobian data its own thing
braintrix2 = braintrix$jacdata

Now, before writing out the braintrix2, we want to residualize it.

# Read in your dems
dems <- read.csv(file = "/data/chamal/projects/lani/adult_thc/dems_with_weight_resids.csv") #, header = TRUE)

# At this point, make sure that the nrow(dems) matches the number of subjects included in the braintrix2 (nrow(braintrix2)). Also, make sure that the order of your subjects in dems is identical to the order of the Jacobians. 
nrow(dems)==nrow(braintrix2) # Check if TRUE

# Then, use Yohan's function to normalize the Jacobians:
# Here, X is the matrix with the Jacobians (nrow is number of subjects, ncol is the voxels). The formula leaves out the left half of the equation, but includes the variables to be normalized for on the right. df is the demographics dataframe which must include the variable (here weight). 
# Generally takes a few minutes to run
normalized_jacs <- norm_by_model(X = braintrix2, formula = ~ weight, df = dems)

# Then I unnamed the Jacobians matrix (but I'm not sure if this is entirely necessary)
unnamed_norm_jacs <- unname(normalized_jacs)

# Write it out into a csv, and this becomes the input for your PLS!
write.csv(unnamed_norm_jacs, file ="./t2_normalized_jacs.csv", row.names=FALSE)