Part A: Haemodynamic response functions (HRF) and block design

Schrick-Noah_CS-6643_Lab-11.R

@@ -0,0 +1,78 @@
+# 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)
+
+