forked from james-cole/brainageR
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpredict_new_data_gm_wm_csf.R
executable file
·74 lines (65 loc) · 3.25 KB
/
predict_new_data_gm_wm_csf.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
#!/usr/bin/env Rscript
## kernlab regression using nifti files
## James Cole 24/09/2019
rm(list = ls())
args <- commandArgs(trailingOnly = TRUE)
## test if there are two argument: if not, return an error
if (length(args) < 6) {
stop("Six arguments must be supplied (brainageR directory full path; list of GM nifti files; list of WM nifti files; list of CSF nifti files, kernlab PCA model; output filename).n", call. = FALSE)
}
## set remote or local for testing. Uncomment accordingly
brainageR_dir <- args[1]
## libraries
library(kernlab)
library(RNifti)
library(stringr)
## get new data and load masks
gm.list <- read.table(file = args[2], header = FALSE, colClasses = "character")$V1
wm.list <- read.table(file = args[3], header = FALSE, colClasses = "character")$V1
csf.list <- read.table(file = args[4], header = FALSE, colClasses = "character")$V1
## read in average smwc files for GM, WM and CSF and convert to vectors
gm_average <- as.vector(readNifti(paste(brainageR_dir, "/software/templates/average_smwc1.nii", sep = "")))
wm_average <- as.vector(readNifti(paste(brainageR_dir, "/software/templates/average_smwc2.nii", sep = "")))
csf_average <- as.vector(readNifti(paste(brainageR_dir, "/software/templates/average_smwc3.nii", sep = "")))
## create nifti vector matrices
## set mask threshold
threshold <- 0.3
## function to read and mask nifti files from list above.
read_mask_nii <- function(arg1){
gm <- as.vector(readNifti(gm.list[arg1]))
gm <- gm[gm_average > threshold]
wm <- as.vector(readNifti(wm.list[arg1]))
wm <- wm[wm_average > threshold]
csf <- as.vector(readNifti(csf.list[arg1]))
csf <- csf[csf_average > threshold]
gm.wm.csf <- c(gm, wm, csf)
return(gm.wm.csf)
}
paste("loading nifti data", date(),sep = " " )
new_data_mat <- matrix(unlist(lapply(1:length(gm.list), function(x) read_mask_nii(x))), nrow = length(gm.list), byrow = TRUE)
dim(new_data_mat)
## load and then apply PCA parameters
rotation <- readRDS(file = paste(brainageR_dir, "/software/pca_rotation.rds", sep = ""))
center <- readRDS(file = paste(brainageR_dir, "/software/pca_center.rds", sep = ""))
scale <- readRDS(file = paste(brainageR_dir, "/software/pca_scale.rds", sep = ""))
newx <- scale(new_data_mat, center, scale) %*% rotation
## load previously trained regression model
paste("loading regression model", date(),sep = " ")
load(args[5])
## generate predictions
brain.ages <- read.csv(paste(brainageR_dir, "/software/brains.train_labels.csv", sep = ""), header = FALSE)
test.res.gpr <- as.data.frame(predict(model.gpr, newx) + mean(brain.ages$V1))
names(test.res.gpr) <- "gpr.brain.age"
## generate prediction confidence intervals
test.sd.gpr <- predict(model.gpr, newx, type = "sdev")
test.res.gpr$lower.CI <- as.numeric(test.res.gpr$gpr.brain.age - (1.96 * test.sd.gpr))
test.res.gpr$upper.CI <- as.numeric(test.res.gpr$gpr.brain.age + (1.96 * test.sd.gpr))
## save predictions to text file
paste("saving new results", date(),sep = " ")
str_sub(gm.list, 1, str_locate(gm.list, "smwc1")[,2]) <- ""
str_sub(gm.list, str_locate(gm.list, ".nii")[,1], str_length(gm.list)) <- ""
write.table(cbind(gm.list, round(test.res.gpr,4)),
file = args[6],
row.names = F,
quote = F,
col.names = c("File", "brain.predicted_age", "lower.CI", "upper.CI"), sep = ",")