# Lab 11 for the University of Tulsa's CS-6643 Bioinformatics Course
# Introduction to fMRI Analysis and ICA
# Professor: Dr. McKinney, Fall 2022
# Noah L. Schrick - 1492657

## Set Working Directory to file directory - RStudio approach
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#### Part A: Haemodynamic response functions (HRF) and block design
## Plot Basic HRF from 0-20s
hq <- function(t,q=4){
  # q=4 or 5, where 5 has more of a delay
  return (t^q * exp(-t)/(q^q * exp(-q)))
}

# use seq to create vector time and use hq to create hrf vectors
time <- seq(0,20)
hrf1 <- hq(time)
hrf2 <- hq(time, 5)

# plot      
plot(time,hrf1,type="l") 
lines(time,hrf2,col="red")

## Deconvolve with task onset times
# grabbed from afni c code
# basis_block_hrf4 from 3dDeconvolve.c
HRF <- function(t, d){
  if (t<0){
    y=0.0
  }else{
    y = 1/256*exp(4-t)*(-24-24*t-12*t^2-4*t^3-t^4 + exp(min(d,t))*(24+24*(t-min(d,t)) + 12*(t-min(d,t))^2+4*(t-min(d,t))^3+(t-min(d,t))^4))
  }
  return(y)
}

t=seq(0,360,len=360)
onsets=c(14,174,254)
blocks.model = double()
for (curr_t in t){
  summed_hrf=0.0
  for (start in onsets){
    summed_hrf=summed_hrf+HRF(curr_t-start,20)
  }
  blocks.model = c(blocks.model,summed_hrf)
}

plot(blocks.model,type="l")


## Plot voxel time series data and the block design curve
voxel.data <- read.delim("059_069_025.1D")
plot(seq(1,2*dim(voxel.data)[1],by=2),t(voxel.data),
     type="l", xlab="time",ylab="intensity")
# normalize the height of the blocks model
blocks.normal <- max(voxel.data)*blocks.model/max(blocks.model)
lines(blocks.normal,type="l",col="red")

# regression
length(blocks.normal) # too long
dim(voxel.data)[1]
# grab elements from blocks.normal to make a vector same as data
blocks.norm.subset <- blocks.normal[seq(1,length(blocks.normal),len=dim(voxel.data)[1])]
length(blocks.norm.subset)

voxel.data.vec <- matrix(unlist(voxel.data),ncol=1)

# use lm to create voxel.fit <- lm...
voxel.fit <- lm(voxel.data.vec ~ blocks.norm.subset)


plot(blocks.norm.subset,voxel.data.vec,
     xlab="block model",ylab="voxel data",main="regression fit")

# use abline(voxel.fit) to overlay a line with fit coefficients 
abline(voxel.fit)


#### Part B: Resting state fMRI visualization, multidimensional arrays and independent component analysis (ICA). 
## example: multi-dimensional array 
#                        2 3x4 matrices
multArray <- array(1:24,dim=c(3,4,2))
dim(multArray)
multArray[,,1] # matrix 1
multArray[,,2] # matrix 2
image(multArray[,,1]) # plot slice 1
image(multArray[,,2]) # plot slice 2

## Neuroimaging Informatics Technology Initiative
if (!require("fastICA")) install.packages("fastICA")
if (!require("fmri")) install.packages("fmri")
library(fastICA)
library(fmri)

# read in 4d nifti 
# uses fmri library, takes about 3min to load 
img <- read.NIFTI("rest_res2standard.nii") 

mask <- img$mask  # Boolean mask for brain voxels
dim(mask)
ttt <- extractData(img) # extract 4d data cube
numScans <- dim(ttt)[4]

# plot a voxel's time series
plot(ttt[30,30,30,],type="l",xlab="time",ylab="activity")

## Plot a 2D Slice
yslice <- 35
scan2dslice <- ttt[,yslice,,50] # grab 2d slice at t=50
image(scan2dslice,main="no masking") # no mask
# mask it off
slice.mask <- mask[,yslice,]
scan2dslice[slice.mask] <- NA # NA's become white
image(scan2dslice,main="masked")

## ICA
t1 <- Sys.time() # for timing purposes
dataMat <- NULL
for(t in seq(1,numScans)){
  scan <- ttt[,,,t]
  # stretched out the 61x73x61 3d matrix into one row
  # apply mask and stack 
  dataMat <- rbind(dataMat,scan[mask])
}
t2 <- Sys.time()
difftime(t2,t1) # 3 minutes
dim(dataMat)

## ICA analysis
# input X: rows observations (voxels) and cols variables (time)
X <- t(dataMat)
m <- 20 # specify number of ICA components
t1<-Sys.time()
f<-fastICA(X,n.comp=m,method="C")
t2<-Sys.time()
difftime(t2,t1) # 1.23min
# S=XKW,
# K is a pre-whitening PCA matrix (components by time)
# S is the matrix of m ICAs (columns of S are spatial signals)
# S has dimensions voxel x components
S<-f$S  # you can find K and W with $

ica.comp <- 5 # look at 5th component 
# plot the 5th time ICA
plot(f$K[,ica.comp],type="l",xlab="time",ylab="signal",main="ICA component")

# threshold S matrix 
theta <- 2
S[S<=theta] <- NA

# turn S back into 4d multidimensional array
xdim<-dim(ttt)[1]
ydim<-dim(ttt)[2]
zdim<-dim(ttt)[3]
ica.4dArray <- array(matrix(S,ncol=1),dim=c(xdim,ydim,zdim,m))
dim(ica.4dArray) # 61x73x61x10
yslice <- 35
# grab 2d slice at y=yslice and ica 5
ica2dslice <- ica.4dArray[,yslice,,ica.comp] 
image(ica2dslice,main="ica component (spatial locations)")