diff --git a/.Rhistory b/.Rhistory
deleted file mode 100644
index 43ad527..0000000
--- a/.Rhistory
+++ /dev/null
@@ -1,512 +0,0 @@
-alpha.LM
-alpha.ML
-degree(g)
-sort(degree(g))
-sort(degree(g),decreasing=FALSE)
-sort(degree(g),decreasing=F)
-sort(degree(g),decreasing=false)
-sort(degree(g), decreasing = TRUE)
-head(sort(degree(g), decreasing = TRUE))
-stddev(degree(g))
-sd(degree(g))
-tail(sort(degree(g), decreasing = TRUE))
-plot(log(g.breaks.clean), log(g.probs.clean))
-# Homework 4 for the University of Tulsa' s CS-7863 Network Theory Course
-# Degree Distribution
-# Professor: Dr. McKinney, Spring 2022
-# Noah Schrick - 1492657
-library(igraph)
-library(igraphdata)
-data(yeast)
-g <- yeast
-g.netname <- "Yeast"
-################# Set up Work #################
-g.vec <- degree(g)
-g.hist <- hist(g.vec, freq=FALSE, main=paste("Histogram of the", g.netname,
-" Network"))
-legend("topright", c("Guess", "Poisson", "Least-Squares Fit",
-"Max Log-Likelihood"), lty=c(1,2,3,4), col=c("#40B0A6",
-"#006CD1", "#E66100", "#D35FB7"))
-g.mean <- mean(g.vec)
-g.seq <- 0:max(g.vec) # x-axis
-################# Guessing Alpha #################
-alpha.guess <- 1.5
-lines(g.seq, g.seq^(-alpha.guess), col="#40B0A6", lty=1, lwd=3)
-################# Poisson #################
-g.pois <- dpois(g.seq, g.mean, log=F)
-lines(g.seq, g.pois, col="#006CD1", lty=2, lwd=3)
-################# Linear model: Least-Squares Fit #################
-g.breaks <- g.hist$breaks[-c(1)] # remove 0
-g.probs <- g.hist$density[-1] # make lengths match
-# Need to clean up probabilities that are 0
-nz.probs.mask <- g.probs!=0
-g.breaks.clean <- g.breaks[nz.probs.mask]
-g.probs.clean <- g.probs[nz.probs.mask]
-plot(log(g.breaks.clean), log(g.probs.clean))
-g.fit <- lm(log(g.probs.clean)~log(g.breaks.clean))
-summary(g.fit)
-alpha.LM <- coef(g.fit)[2]
-lines(g.seq, g.seq^(-alpha.LM), col="#E66100", lty=3, lwd=3)
-################# Max-Log-Likelihood #################
-n <- length(g.breaks.clean)
-kmin <- g.breaks.clean[1]
-alpha.ML <- 1 + n/sum(log(g.breaks.clean/kmin))
-alpha.ML
-lines(g.seq, g.seq^(-alpha.ML), col="#D35FB7", lty=4, lwd=3)
-plot(log(g.breaks.clean), log(g.probs.clean))
-g.breaks.clean <- g.breaks[nz.probs.mask]
-g.probs.clean <- g.probs[nz.probs.mask]
-plot(log(g.breaks.clean), log(g.probs.clean))
-BiocManager::install()
-simDats <- createSimulation2(data.type = "discrete", avg.maf = 0.2, sim.type = "mainEffect",
-pct.train = 0.5, pct.holdout = 0.5, pct.validation = 0,
-main.bias = 0.4, pct.signals = 0.2)
-# npdro::createSimulation2() example for gwas simulation
-library(npdro)
-# npdro::createSimulation2() example for gwas simulation
-if (!require("npdro")) install.packages("npdro")
-library(npdro)
-# npdro::createSimulation2() example for gwas simulation
-install_github("insilico/npdro")
-# npdro::createSimulation2() example for gwas simulation
-if (!require("devtools")) install.packages("devtools")
-library(devtools)
-install_github("insilico/npdro")
-# npdro::createSimulation2() example for gwas simulation
-if (!require("devtools")) install.packages("devtools")
-#### 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(1,20)
-time
-# use seq to create vector time and use hq to create hrf vectors
-time <- seq(0,20)
-hq
-# 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(blocks.model, type="l")
-lines(time,hrf1,col="red")
-pwd
-ls
-voxel_time <- scan("059_069_025.1D", character(), quote = "\n")
-## Set Working Directory to file directory - RStudio approach
-setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
-voxel_time <- scan("059_069_025.1D", character(), quote = "\n")
-voxel_time
-plot(blocks.model, type="l")
-lines(voxel_time,hrf1,col="red")
-voxel_time <- scan("059_069_025.1D", character(), quote = "\n")
-t=seq(0,max(voxel_time))
-blocks.model = double()
-for (curr_t in voxel_time){
-summed_hrf=0.0
-for (start in onsets){
-summed_hrf=summed_hrf+HRF(curr_t-start,20)
-}
-blocks.model = c(blocks.model,summed_hrf)
-}
-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")
-blocks.model = double()
-for (curr_t in t){
-summed_hrf=0.0
-for (start in voxel_time){
-summed_hrf=summed_hrf+HRF(curr_t-start,20)
-}
-blocks.model = c(blocks.model,summed_hrf)
-}
-voxel_time
-onsets
-class(onsets)
-class(voxel_time)
-voxel_time <- as.numeric(scan("059_069_025.1D", character(), quote = "\n"))
-t=seq(0,max(voxel_time))
-blocks.model = double()
-for (curr_t in t){
-summed_hrf=0.0
-for (start in voxel_time){
-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)
-plot(blocks.norm.subset,voxel.data.vec,
-xlab="block model",ylab="voxel data",main="regression fit")
-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)
-voxel.data.vec <- matrix(unlist(voxel.data),ncol=1)
-length(blocks.norm.subset)
-voxel.data.vec <- matrix(unlist(voxel.data),ncol=1)
-length(voxel.data.vec)
-plot(blocks.norm.subset,voxel.data.vec,
-xlab="block model",ylab="voxel data",main="regression fit")
-plot(blocks.norm.subset,voxel.data.vec,
-xlab="block model",ylab="voxel data",main="regression fit")
-blocks.norm.subset
-voxel.data.vec
-voxel.fit <- lm(blocks.norm.subet ~ voxel.data.vec)
-voxel.fit <- lm(blocks.norm.subset ~ voxel.data.vec)
-# use abline(voxel.fit) to overlay a line with fit coefficients
-abline(voxel.fit)
-## 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")
-# 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)
-voxel.fitr
-voxel.fit
-voxel.data
-blocks.normal
-blocks.model
-max(blocks.model)
-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
-# 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)
-# 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)
-blocks.model
-count(blocks.model > 1)
-which(blocks.model > 1)
-sum(blocks.model < 1)
-length(blocks.model)
-100* sum(blocks.model < 1)/length(blocks.model)
-#### 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
-dim(multArray)
-size(multArray)
-length(multArray)
-multArray
-multArray[,,1] # matrix 1
-multArray[,,2] # matrix 2
-image(multArray[,,1]) # plot slice 1
-image(multArray[,,2]) # plot slice 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("AnalyzeFMRI")) install.packages("AnalyzeFMRI")
-if (!require("fmri")) install.packages("fmri")
-#    If this fails, try library(devtools) (install.packages("devtools"))
-#    and then install_github(https://github.com/cran/fmri.git)
-library(fmri)
-## Neuroimaging Informatics Technology Initiative
-if (!require("fastICA")) install.packages("fastICA")
-if (!require("AnalyzeFMRI")) install.packages("AnalyzeFMRI")
-install.packages("R.matlab")
-install.packages("Biostrings")
-BiocManager::install("Biostrings")
-BiocManager::install("Biostrings")
-BiocManager::install("AnalyzeFMRI")
-BiocManager::install()
-if (!require("fmri")) install.packages("fmri")
-if (!require("fmri")) install.packages("aws")
-if (!require("fmri")) install.packages("gsl")
-install.packages("gsl")
-install.packages("GSL")
-install_github(https://github.com/cran/fmri.git)
-install_github("https://github.com/cran/fmri.git")
-#    If this fails, try library(devtools) (install.packages("devtools"))
-#    and then install_github(https://github.com/cran/fmri.git)
-if (!require("devtools")) install.packages("devtools")
-library(devtools)
-install_github("https://github.com/cran/fmri.git")
-install_github("https://github.com/cran/aws.git")
-install_github("https://github.com/cran/gsl.git")
-#    If this fails, try library(devtools) (insta
-install.packages("gsl")
-install.packages("fmri")
-if (!require("AnalyzeFMRI")) install.packages("AnalyzeFMRI")
-if (!require("fastICA")) install.packages("fastICA")
-library(fastICA)
-# read in 4d nifti
-# uses fmri library, takes about 3min to load
-img <- read.NIFTI("rest_res2standard.nii")
-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")
-length(ttt)
-ttt
-dim(ttt)
-length(ttt)
-numScans
-## 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")
-mask
-slice.mask
-## 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)
-dataMat
-dataMat[1,]
-dataMat[99,]
-## 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 $
-S
-dim(S)
-dim(f$K)
-dim(f$W)
-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)")