diff --git a/Schrick-Noah_CS-6643_Lab-11.R b/Schrick-Noah_CS-6643_Lab-11.R new file mode 100644 index 0000000..45e98ea --- /dev/null +++ b/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) + + diff --git a/Schrick-Noah_CS-6643_Lab-11.docx b/Schrick-Noah_CS-6643_Lab-11.docx index c2da1df..6fe8202 100644 Binary files a/Schrick-Noah_CS-6643_Lab-11.docx and b/Schrick-Noah_CS-6643_Lab-11.docx differ