-
Notifications
You must be signed in to change notification settings - Fork 0
/
BLDI_Main_Script_voxel_wise_ttest.R
110 lines (82 loc) · 4.28 KB
/
BLDI_Main_Script_voxel_wise_ttest.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
# Bayes Factor Lesion-Deficit Mapping
# Main script
####################################################
#
# Within this script, the main analysis parameters and the
# paths to the data have to be defined
####################################################
# empty working space
rm(list = ls())
# load packages
print("Loading dependencies")
library(RNifti)
library(BayesFactor)
####################################################
# start of input required by user
####################################################
# provide the path to the folder that contains the lesions in NIFTI format
# (Attention: the symbol - R uses "/" (NOT: "\") in path names!)
path = "D:/Folder/Data"
# provide the path to the csv file that contains the scores
# (you can change the separator by changing the sep argument)
dat_scores <- read.csv("D:/Folder/Deficit_Scores.csv",
header = TRUE,sep = ";")
# define the minimum number of lesions in a voxel for the voxel to be included
min_vox_threshold <- 5
# provide the path to folder of BLDI main script
script_path <- "D:/Folder/BLDI"
####################################################
# end of input required by user
####################################################
# flag for binary lesion data - 1 yes, 0 no
binary_data <- 1
# set working directory to path of the scripts
setwd(script_path)
source(paste(script_path,"/bldi_create_analysis_mask.R",sep = ""))
source(paste(script_path,"/bldi_load_lesions.R",sep = ""))
source(paste(script_path,"/bldi_create_ttest_bf_map.R",sep = ""))
# create a lesion overlap and get the index of all voxels to be included in the analysis
output_overlap <- bldi_create_analysis_mask(dat_scores,path,min_vox_threshold)
idx_voxels_included<-unlist(output_overlap[1])
# report number of voxels
writeLines(sprintf("%d voxles are damaged in at least %d patients \n and will be included in the analysis",
length(idx_voxels_included),min_vox_threshold))
####################################################
# load the lesions
lesions_full <- bldi_load_lesions(dat_scores,path,idx_included_voxels,binary_data)
print("Done - Loading lesion data")
print(sprintf("Starting Bayes Factor mapping for score %s", colnames(dat_scores)[2]))
print("The procedure may take up to several hours! Please wait!")
BF_result <- bldi_create_ttest_bf_map(dat_scores,idx_voxels_included,lesions_full)
####################################################
# re-load the first lesion as a reference for image image
lesion_path <- paste(path, "/", dat_scores[1,1], ".nii", sep="")
lesion <- readNifti(lesion_path)
dimension_1st_lesion <- dim(lesion)
lesion_vect_length <- length(c(lesion))
# rebuild a 3D image of the BF
BF_full_vect <- matrix(0, lesion_vect_length)
BF_full_vect[idx_voxels_included] <- BF_result
# generate a NIFTI of the BF map
BF_full_image <- array(BF_full_vect, dimension_1st_lesion)
# take the first lesion image as reference for header info
BF_image <- asNifti(BF_full_image, reference = lesion_path)
filename <- paste("Output_BLDI_BF_map_ttest_",as.character(colnames(dat_scores)[2]),"_",format(Sys.time(), "%Y-%m-%d_%H_%M"),sep = "")
writeNifti(BF_image, filename)
# additionally, create a log(BF) map
filename <- paste("Output_BLDI_logBF_map_ttest_",as.character(colnames(dat_scores)[2]),"_",format(Sys.time(), "%Y-%m-%d_%H_%M"),sep = "")
writeNifti(log10(BF_image), filename)
## create a results txt with some additional information
filename <- paste("Output_BLDI_info_ttest_",as.character(colnames(dat_scores)[2]),"_",format(Sys.time(), "%Y-%m-%d_%H_%M"), ".txt",sep = "")
writeLines("Results of voxel-wise BLDI with t-tests", filename)
line <- paste("For deficit",as.character(colnames(dat_scores)[2]))
write(line,file = filename,append = TRUE)
write("",file = filename,append = TRUE)
line <- paste("Voxels with at least ",as.character(min_vox_threshold),"lesions were analysed")
write(line,file = filename,append = TRUE)
line <- paste(as.character(length(idx_voxels_included)),"voxels were analysed")
write(line,file = filename,append = TRUE)
line <- paste("The maximum BF is ", as.character(max(BF_result)), ", the minimum BF is ", as.character(min(BF_result)))
write(line,file = filename,append = TRUE)
## end of txt file creation
print("Done - Analyses finished!")