###############################################################################
# code to load, preprocess, analyze the dataset associated with the paper     #
# "Repeated double cross validation applied to the PCA-LDA classification of  #
# SERS spectra: a case study with serum samples from hepatocellular carcinoma #
# patients" by Elisa Gurian , Alessia Di Silvestre, Elisa Mitri, Devis Pascut #
# Claudio Tiribelli, Mauro Giuffre', Lory Saveria Croce', Valter Sergo and    #
# Alois Bonifacio (abonifacio@units.it).                                      #
#                                                                             #
# Code by Alois Bonifacio using R ver 4.0.3 (2020-10-10)                      #
# platform x86_64-w64-mingw32/x64,                                            #   
# run on a Intel(R) Core(TM) i7-6500U CPU @ 2.50GHz, 2592 Mhz, 2 cores        #
# with a Microsoft Windows 10 Home and RStudio Version 1.3.1093               #  
#                                                                             #
# the code allows to analyze the data preduce the figures and tables as shown #
# in the paper, although since each run is different, there might be slight   #
# differences.                                                                #
###############################################################################

###############################################################################
# ATTENTION!!! The code requires the following libraries for R installed:     #
# hyperSpec, baseline, caret, ROCR, MASS, e1071, MALDIquant, ggplot2, binom,  #
# cvAUC, ggsignif.                                                            # 
###############################################################################

# PART 1: load dataset ########################################################
# to import data from 144 ASCII files containing spectra

require (hyperSpec) # call hyperSpec library

# specify the folder with the 144 txt spectra 
# NOTE: ACTION NEEDED!!! YOU NEED TO CHANGE THIS WITH YOUR PATH!!!
mydir <- "c:/users/Alois Bonifacio/mydata_folder"
setwd(mydir)

txt_files <- list.files(pattern=".*txt") #get filenames of files containing spectra
spec_table <- read.table(txt_files[1]) #import first spectrum as table
spc <- new ("hyperSpec", wavelength = spec_table[,1],spc = spec_table[,2]) #create new hyperSpec object spc 
spc$date <- substr(txt_files[1],1,8) #extract metadata from filename: date of acquisition
spc$batch <- substr(txt_files[1],16,16) #substrate batch: A,B or C
spc$class <- substr(txt_files[1],18,20) #class: CTR (control) or H0T (cancer)
spc$sample_code <- substr(txt_files[1],22,24) #sample code
rm(spec_table)

for (i in c(2:length(txt_files))){ #for cycle to load all files one by one, merging them in spc
  spec_table <- read.table(txt_files[i]) #import i-th spectrum
  spc_temp <- new ("hyperSpec", wavelength = spec_table[,1],spc = spec_table[,2]) #create temporary hyperSpec obj
  spc_temp$date <- substr(txt_files[i],1,8) #extract metadata from i-th filename
  spc_temp$batch <- substr(txt_files[i],16,16)
  spc_temp$class <- substr(txt_files[i],18,20)
  spc_temp$sample_code <- substr(txt_files[i],22,24)
  spc <- rbind(spc, spc_temp) # merge temporary object to spc hyperSpec object
  rm(spec_table) # remove temporary objects
  rm(spc_temp)
}

rm(i); rm(txt_files) #remove un-necessary temporary objects

pos.class <- "H0T" #specify positive class
neg.class <- "CTR" #specify negative class


# PART 2: data preprocessing ##################################################
# to crop spectral region, interpolate data, fit and subtract a baseline 
# and normalize spectra

require(baseline)

# 2.a) crop spectral region and loess interpolate (smoothing)
spc <- spc.loess(spc, seq(400,1800,2)) 

# 2.b) baseline fitting and removal
bl <- baseline (spc [[]], method = "modpolyfit", degree = 4) #fit baselines
spc$spc <- getCorrected(bl)  #subtract baselines
rm(bl) #remove object with fitted baselines
spc <- spc[,,430~1730] # crop data to leave "baseline artifacts" out

# 2.c) normalization (vector normalization)
factor <- apply (spc, 1,function(x){sqrt(sum(x^2))}) #calculate factors for normalization
spc <- sweep (spc, 1, factor, "/") #normalize intensity values
rm(factor) #remove temporary factors


# PART 3: Repeated Double Cross Validation (RDCV) of PCA-LDA models ###########
# NOTE: the strucuture of the RDCV is that of a triple nested loop, as from
# the RDCV paper from Filzmoser, Liebmann, Varmuza, J.Chemometrics 2009, 23, 160-171

X <- data.frame(spc$spc) #create a dataframe "X" with all intensities and Raman shifts 
colnames(X) <- spc@wavelength #name columns as wavenumbers
grp <- as.factor(spc$class) #create a separate vector "grp" with class labels

# set the parameters for the RDCV loops
# note: segments (for both inner and outer loops) are automatically stratified
repetitions <- 100 # number of repetitions of the cross-validation (typically 20-100)
out.segments <- 3 # number of segments of the OUTER LOOP (i.e. the cross-validation) (typically 3-10 segments)
inn.segments <- 7 # number of segments of the INNER LOOP (i.e. the parameter optimization loop) (typically 7-10)
param <- c(1:7) # set the paramter range (i.e. number of principal components, PCs) to consider

# creation of a series of empty objects (vectors,arrays) to store results generated during the RDCV
optpar <- NA # to store optimized parameters
PCAscores <- array(NA, c(nrow(X),20,out.segments*repetitions)) # to store PCA scores
PCAload <- array(NA, c(ncol(X),20,out.segments*repetitions)) # to store PCA loadings
res.auc <- NA # to store AUC values
res.mat <- array(NA, c(length(levels(grp)),length(levels(grp)), # res.mat to store confusion matrices
                       out.segments*repetitions), dimnames = list(levels(grp),levels(grp))) 
# res.pred to store reference values and predictions
res.pred <-  array(NA, c(ceiling(length(grp)/(out.segments-1)), 2, out.segments*repetitions))
colnames(res.pred) <- c("ref","pred")
# res.roc to store ROC curves values
res.roc <-  array(NA, c(ceiling(length(grp)/(out.segments-1)), 2, out.segments*repetitions))
colnames(res.roc) <- c("FPR","TPR")
# create dataframe to store results from the (repeated) INNER LOOP CV
res.cv <- array(NA, c(length(param), 3, out.segments*repetitions))
colnames(res.cv) <- c("par", "meanCVerr", "sdCVerr")
# create data frame with all the prediction probabilities for each repetition
pred.prob <- array(NA, c(nrow(X), repetitions))
# create data frame to store the LD scores
LD_scores <- data.frame(array(NA, c(nrow(spc), repetitions*out.segments)))
colnames(LD_scores) <- c(1:(repetitions*out.segments))
  
# call libraries required for the RDCV loops
require(caret) 
require(ROCR) 
require(MASS)
require(e1071)

i <- 0 #create incremental index (for each fold, outer loop)
j <- 0 #create a second incremental index (for each iteration, inner loop)


# start REPETITION loop 
for (r in seq(1:repetitions)){
  tfolds <- createFolds(grp, k=out.segments, list=TRUE) # create data folds
  
  # start OUTER loop 
  # create optimization (t-1 segments) and test (1 segment) partitions
  for (t in seq(1:out.segments)){
    i <- i + 1 # increment overall index
    test <- sort(unlist(tfolds[t], use.names=F)) # test segment
    opt <- sort(unlist(tfolds[-t], use.names=F)) # optimization segment
    # split optimization set into k segments (automatically stratified)
    kfolds <- createFolds(grp[opt], k=inn.segments, list=TRUE) # create data folds
    # create empty array to store CV misclassification errors for each parameter
    cvERR <-  array(NA, c(inn.segments, 1, length(param)))
    
    # start INNER loop
    for (p in c(1:length(param))){ #loop for the different parameters 
      for (k in seq(1:inn.segments)){
        j <- j+1 # increment index j
        # create validation (1 segment) and training set (k-1 segments)
        val <- unlist(kfolds[k], use.names=F) # validation segment
        train <- unlist(kfolds[-k], use.names=F) # train segment
        # create model with traning set and apply to validation set
        # set the preprocessing model (PCA) using only the training data
        PCAmodel <- prcomp(X[opt,][train,], scale=FALSE, center=TRUE) #data are NOT scaled
        # process/scale the training data (PCA)
        PCAscores_train <- PCAmodel$x
        # process/scale the test data (PCA)
        PCAscores_val <- predict(PCAmodel, X[opt,][val,])
        # create the LDA model with the training set
        train_df <- data.frame(PCAscores_train[,1:param[p]])
        colnames(train_df) <- c(1:param[p])
        rule <- lda(train_df, grouping = grp[opt][train],  CV = FALSE)
        # predict the test set using the LDA model
        val_df <- data.frame(PCAscores_val[,1:param[p]])
        colnames(val_df) <- c(1:param[p])
        pred.k <- predict(rule, newdata = val_df) # predictions
        # calculate misclassification error
        cm <- as.matrix(table(pred.k$class,grp[opt][val]))
        delta <- row(cm) - col(cm)
        misclas <- cm[delta < 0 | delta > 0]
        cvERR[k,,p] <- sum(misclas)/sum(cm) #store errors in a temporary object cvERR
      } # end loop
    } # end INNER loop
    # remove un-necessary object
    rm(val);rm(train); rm(PCAmodel); rm(PCAscores_train);rm(PCAscores_val); rm(p);rm(k)
    rm(train_df); rm(val_df); rm(rule); rm(pred.k); rm(cm); rm(delta); rm(misclas)
    
    # store results from the (repeated) INNER loop 
    res.cv[,1,i] <- param
    res.cv[,2,i] <- apply(cvERR,3, mean)
    res.cv[,3,i] <- apply(cvERR,3, sd)
    
    # calculate threshold according to the "one-standard-error-rule"
    CVmin <- res.cv[res.cv[,2,i] == min(res.cv[,2,i]),,i]
    if(class(CVmin)=="numeric"){CVmin <- t(as.matrix(CVmin))}
    threshold <- (CVmin[1,2] + CVmin[1,3]) # set threshold (1 std. error)
    # identify/store optimum paramter
    optparam <- res.cv[res.cv[,2,i] < threshold,,i]
    if (class(optparam)=="numeric"){optparam <- t(as.matrix(optparam))}
    optpar[i] <- optparam[1,1] #optimum paramter (according to "one standard error rule")
    
    # build model and predict test segment with optimized parameter
    # set the preprocessing model (PCA) using only the training data 
    PCAmodel <- prcomp(X[opt,], scale=FALSE, center=TRUE) #data are NOT scaled
    # process/scale the training data (PCA) - store scores for the opt (training) set
    PCAscores_opt <- PCAmodel$x
    # store the PC scores (first 20 PCs)
    PCAscores[,,i][1:nrow(PCAmodel$x),] <- PCAmodel$x[,1:20]
    # store the PC loadings (first 20 PCs)
    PCAload[,,i] <- PCAmodel$rotation[, 1:20]
    # build a dataframe with the classes and the scores for the train set
    PCAscores_train_df <- data.frame(grp[opt], PCAmodel$x[,1:20])
    # process/scale the test data (PCA) 
    # NOTE: test data scores are predicted using the train PCA model, so that no
    # information from the test is used to get the scores
    PCAscores_test <- predict(PCAmodel, X[test,])
    # make predictions with optimized model
    rule <- lda(data.frame(PCAscores_opt[,1:optpar[i]]), grouping = grp[opt],  CV = FALSE)
    pred <- predict(rule, data.frame(PCAscores_test[,1:optpar[i]])) # predictions
    # store LD scores (LD1)
    LD_scores[test,i] <- pred$x
    # scores or probabilities (for ROC)
    pred.prob[test,r] <- pred$posterior[,colnames(pred$posterior)==pos.class]
    # store results 
    cm <- confusionMatrix(data=as.factor(pred$class), reference=grp[test], positive=pos.class) # extract confusion matrix
    res.mat[,,i] <- cm$table # store confusion matrix
    res.pred[,,i][1:length(grp[test]),1] <- grp[test] # store reference class labels of TEST
    res.pred[,,i][1:length(pred$class),2] <- as.factor(pred$class) # store predictions for TEST
    # FOR ROC calculation (using ROCR package)
    prob.roc <- prediction(pred.prob[test,r], grp[test])
    perf.roc <- performance(prob.roc, "tpr","fpr")
    res.auc[i] <- performance(prob.roc, "auc")@y.values[[1]] # store AUC values
    # store ROC values
    res.roc[,,i][1:length(perf.roc@x.values[[1]]),1] <- perf.roc@x.values[[1]] 
    res.roc[,,i][1:length(perf.roc@y.values[[1]]),2] <- perf.roc@y.values[[1]]
    # monitor progress by printing loop indexes
    print (paste("repetition =",r,"/", repetitions, "  fold =",t, "/", out.segments, "  model",i)) 
  } # end OUTER loop
  # remove un-necessary objects
  rm(test); rm(opt); rm(kfolds); rm(cvERR); rm(CVmin); rm(threshold); rm(optparam)
  rm(PCAmodel); rm(PCAscores_opt); rm(PCAscores_test); rm(PCAscores_train_df)
  rm(rule); rm(cm); rm(prob.roc); rm(perf.roc); rm(pred); rm(t)
} # end REPETION loop 
# remove un-necessary objects
rm(tfolds)
rm(i);rm(j);rm(r)

# POST PROCESS the results in new objects
# merge confusion matrices over folds in res.mat
res.mat_folds <- res.mat
res.matFold <- array(NA, c(2, 2, repetitions)) 
nfold <- out.segments
j <- 1
for (i in seq(1,(nfold*repetitions), by =nfold)){
  res.matFold[,,j] <- apply(res.mat[,,c(i:(i+(nfold-1)))], c(1,2), sum)
  j <- j+1
}
res.mat <- res.matFold
rm(i); rm(j); rm(res.matFold); rm(res.mat_folds)
# average ROCs over folds (func = mean)
res.rocFold <- array(NA, c(nrow(res.roc), 2, repetitions)) 
nfold <- out.segments
j <- 1
for (i in seq(1,(nfold*repetitions), by =nfold)){
  res.rocFold[,,j] <- apply(res.roc[,,c(i:(i+(nfold-1)))], c(1,2), mean)
  j <- j+1
}
res.roc <- res.rocFold
rm(res.rocFold)


# PART 4: CREATE FIGURES and TABLES FROM RDCV RESULTS #########################

# Figure 1 ####################################################################
# Comparison between the medians of SERS spectra of serum from H0T (n = 72) and 
# CTR groups (n=72). Interquartile ranges of the SERS intensity for the two groups 
# are shown as shaded areas. Medians and interquartile ranges were calculated from 
# intensity normalized spectra. The intensity difference between H0T and CTR is 
# reported in the lower part of the figure.

# calculate interquartiles and median for neg.class
med.neg <- apply(spc[spc$class==neg.class]$spc, 2, median)
qt3n <- apply(spc[spc$class==neg.class]$spc, 2, function(x) quantile(x, 3/4))
qt1n <- apply(spc[spc$class==neg.class]$spc, 2, function(x) quantile(x, 1/4))

# calculate interquartiles and median for pos.class
med.pos <- apply(spc[spc$class==pos.class]$spc, 2, median)
qt3p <- apply(spc[spc$class==pos.class]$spc, 2, function(x) quantile(x, 3/4))
qt1p <- apply(spc[spc$class==pos.class]$spc, 2, function(x) quantile(x, 1/4))

# calculate all difference spectra (pos.class - neg.class) - might take some time
spc.diff <- matrix(data = NA, 
                   nrow = nrow(spc[spc$class == pos.class])*nrow(spc[spc$class == neg.class]),
                   ncol = ncol(spc$spc))
k <- 1
for(i in (1:nrow(spc[spc$class == pos.class]))){
  for(j in (1:nrow(spc[spc$class == neg.class]))){
    spc.diff[k,] <- as.vector(spc[spc$class == pos.class][i]$spc - spc[spc$class == neg.class][j]$spc)
    k <- k+1}}

# calculate median and interquartile for difference spectra
diff.med <- apply(spc.diff, 2, median)
diff.qt3 <- apply(spc.diff, 2, function(x) quantile(x, 3/4))
diff.qt1 <- apply(spc.diff, 2, function(x) quantile(x, 1/4))

# generate PLOT
require(MALDIquant) # call library for peak labeling
X11(width = 6, height = 6)
# plot median and interquartiles for pos and neg class
plot(spc@wavelength, med.neg, type="l", ylim=c(-0.1,0.18), yaxt="n",
     xlab=expression(paste("Raman shift (", cm^-1, ")", sep = "")), ylab="")
segments(spc@wavelength,qt1n,spc@wavelength,qt3n, col=rgb(0,0,1,0.1), lwd=2)
lines(spc@wavelength, med.neg, col="blue", lwd=1)
segments(spc@wavelength,qt1p,spc@wavelength,qt3p, col=rgb(1,0,0,0.1), lwd=2)
lines(spc@wavelength, med.pos, col="red", lwd=1)
# label peaks
mspec1 <- createMassSpectrum(as.vector(spc@wavelength),as.vector(med.neg))
peaks1 <- detectPeaks(mspec1,halfWindowSize=11, method=c("MAD", "SuperSmoother"), SNR=0.8)
labelPeaks(peaks1, digits=0, underline=FALSE,  srt=90,avoidOverlap=FALSE , adj=c(-1.0,0.4)) 
# plot median and interquartile for difference spectra
diff.offset <- 0.05 # set offset for difference spectra
segments(spc@wavelength,diff.qt1-diff.offset,spc@wavelength,diff.qt3-diff.offset, col="darkgrey", lwd=2)
lines(spc@wavelength, diff.med-diff.offset)
# label peaks
mspec1 <- createMassSpectrum(as.vector(spc@wavelength),as.vector(diff.med))
peaks1 <- detectPeaks(mspec1,halfWindowSize=11, method=c("MAD", "SuperSmoother"), SNR=0.8)
labelPeaks(peaks1, absoluteVerticalPos = peaks1@intensity - diff.offset, 
           srt=90,avoidOverlap=FALSE, digits=0, underline=FALSE, adj=c(-0.6,0.4)) 
mspec1 <- createMassSpectrum(as.vector(spc@wavelength),as.vector(-diff.med))
peaks1 <- detectPeaks(mspec1,halfWindowSize=11, method=c("MAD", "SuperSmoother"), SNR=0.8)
labelPeaks(peaks1, absoluteVerticalPos = -peaks1@intensity - diff.offset,
           digits=0, underline=FALSE, srt=90,avoidOverlap=FALSE , adj=c(1.6,0.4)) 
# draw lines 
abline(h=0)
abline(h=0-diff.offset, lty=2)
# draw legends and labels
title(ylab="Normalized intensity (a.u.)", line=1, cex.lab=1)
text(1640, -0.1, label = paste(pos.class,"-", neg.class), cex=0.9)
legend(1450, 0.18, legend = c(neg.class, pos.class), lty=1, col=c("blue", "red"),
       border="white", box.col="white", cex=0.9, y.intersp=1.4)
# remove un-necessary temporary objects
rm(med.neg); rm(med.pos); rm(qt3n); rm(qt3p); rm(qt1n);rm(qt1p); rm(spc.diff)
rm(diff.med); rm(diff.qt1); rm(diff.qt3); rm(mspec1); rm(peaks1); rm(diff.offset)
rm(k); rm(i); rm(j)


# Figure 2A ###################################################################
# Characterization of the PCA-LDA models produced by the RDCV: (A) curves for the 
# inner cycle of the RDCV, showing the cross-validation error (CVerr) when using 
# a different number of PC

# plot CV curves ALL TOGETHER plus median and IQR
X11(width = 6, height = 6)
cols <- rgb(0,0,0,0.02)
plot(res.cv[,,1], type="l", col=cols, xlab="Principal Components", 
       ylab="CVerr", lwd=2, ylim=c(min(res.cv[,2,]),max(res.cv[,2,])))
for (p in c(2:dim(res.cv)[3])){
  lines(res.cv[,,p], type="l", col=cols, lwd=2)
}
lines(param, apply(res.cv[,2,],1,median), lwd=2)
lines(param, apply(res.cv[,2,],c(1), FUN=function(x) quantile(x, 1/4)), lty=2)
lines(param, apply(res.cv[,2,],c(1), FUN=function(x) quantile(x, 3/4)), lty=2)
rm(cols); rm(p)
  

# Figure 2B ###################################################################
# Characterization of the PCA-LDA models produced by the RDCV: (B) frequency plot
# for optimized models, showing the number of models generated (i.e. Frequency) 
# using a specific number of PC, as a consequence of model optimization.

X11(width = 6, height = 6)
parfreq <- data.frame(param, 0)
colnames(parfreq) <- c("optpar", "Freq")
parfreq[match(data.frame(table(optpar))[,1], parfreq[,1]),2] <- data.frame(table(optpar))[,2]
par(lend=2)
plot(parfreq[,2], type="h", lwd=10, xlab="Optimal number of Principal Components", 
     col="darkgrey", ylab="Frequency", xaxt="n")
axis(side=1, at = c(1:nrow(parfreq)), las=1, labels=parfreq[,1])
rm(parfreq)


# Figure 3 ####################################################################
# Statistics for the Confusion Matrices resulting from the predictions of the 
# RDCV optimized models. Median values are shown in red.
  
# name the values in the res.mat
TN <- res.mat[1,1,]
TP <- res.mat[2,2,]
FP <- res.mat[2,1,]
FN <- res.mat[1,2,]
# create dataframe for ggplot
CM.df <- data.frame(c(TP, FP, FN, TN), 
                    c(rep("TP",length(TP)), rep("FP",length(FP)),
                      rep("FN",length(FN)), rep("TN",length(TN))))
colnames(CM.df) <- c("counts", "type")
CM.df$true <- NA
CM.df$pred <- NA
CM.df[CM.df$type=="TN" | CM.df$type=="FP",]$true <- "neg"
CM.df[CM.df$type=="FN" | CM.df$type=="TP",]$true <- "pos"
CM.df[CM.df$type=="TN" | CM.df$type=="FN",]$pred <- "neg"
CM.df[CM.df$type=="FP" | CM.df$type=="TP",]$pred <- "pos"
# to set the order
CM.df$true_f <- factor(CM.df$true, levels=c("pos", "neg"))
CM.df$pred_f <- factor(CM.df$pred, levels=c("pos", "neg"))
# new facet label names 
classnames <- c("H0T", "CTR")
names(classnames) <- c("pos", "neg")
# set horizontal coordinate for medians labels
y0 <- max(hist(TP, breaks=seq(0,unique(TP + FN), 2), include.lowest=TRUE, plot=FALSE)$counts)/2
  
X11(width = 6, height = 6)   
require(ggplot2)
ggplot(CM.df, aes(x = counts)) +
  theme_bw() +
  theme(plot.margin = unit(c(3,3,1,1), "lines"), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  geom_histogram(binwidth = 2, color = "darkgrey", fill = "grey") +
  xlim(0,unique(TP + FN)) +
  labs (x="Number of elements", y="Frequency") + 
  geom_text(x=unique(TP + FN)/2, y=y0, label=median(TP), size=5, 
            color="red",data = subset(CM.df, type =="TP")) +
  geom_text(x=unique(TP + FN)/2, y=y0, label=median(TN), size=5, 
            color="red",data = subset(CM.df, type =="TN")) +
  geom_text(x=unique(TP + FN)/2, y=y0, label=median(FP), size=5, 
            color="red",data = subset(CM.df, type =="FP")) +
  geom_text(x=unique(TP + FN)/2, y=y0, label=median(FN), size=5, 
            color="red",data = subset(CM.df, type =="FN")) +
  geom_vline(aes(xintercept = median(TP)), data = subset(CM.df, type =="TP"),  colour="red") +
  geom_vline(aes(xintercept = median(TN)), data = subset(CM.df, type =="TN"),  colour="red") +
  geom_vline(aes(xintercept = median(FP)), data = subset(CM.df, type =="FP"),  colour="red") +
  geom_vline(aes(xintercept = median(FN)), data = subset(CM.df, type =="FN"),  colour="red") +
  facet_grid(cols=vars(true_f), rows=vars(pred_f), labeller = labeller(true_f=classnames, pred_f=classnames)) 
grid.text("TRUE", x = 0.475, y = 0.93)
grid.text("PREDICTED", x = 0.93, y = 0.475, rot=270)
# remove un-necessary temporary objects
rm(TN); rm(TP); rm(FN); rm(FP); rm(y0); rm(CM.df); rm(classnames)
           
  
# Table 2 #####################################################################
# Figures of merit calculated from the optimized models generated from the RDCV

# create empty variables to store results
acc <- NA
sen <- NA
spe <- NA
ppv <- NA
npv <- NA
# name elements according to caret convention for confusionMatrix
#  A B
#  C D
A <- res.mat[1,1,]
B <- res.mat[1,2,]
C <- res.mat[2,1,]
D <- res.mat[2,2,]
# calculate FOMs
spe <- A/(A+C) # for gr B positive
sen <- D/(B+D)
acc <- (A+D)/(A+B+C+D)
ppv <- D/(C+D)
npv <- A/(A+B)
# store results in a single data.frame (the last is for ROC values, calculated later)
CVresults <- cbind(acc, sen, spe, ppv, npv, NA)
# calculate expected values (means) and 95% CI for FOMs
require(binom)
# sensitivity
p <- mean(sen)
n <- unique(B+D)
CIsen <- binom.confint(x=p*n, n=n, conf.level = 0.95)
# specifity
p <- mean(spe)
n <- unique(A+C)
CIspe <- binom.confint(x=p*n, n=n, conf.level = 0.95)
# accuracy
p <- mean(acc)
n <- unique(A+B+C+D)
CIacc <- binom.confint(x=p*n, n=n, conf.level = 0.95)
# PPV
p <- mean(ppv)
n <- unique(C+D)
CIppv <- binom.confint(x=p*n, n=n, conf.level = 0.95)
# NPV
p <- mean(npv)
n <- unique(A+B)
CInpv <- binom.confint(x=p*n, n=n, conf.level = 0.95)
# calculate confidence intervals for AUC
require(cvAUC)
# create matrix with labels (reference) with same dimensions as glm.probs
cilab <- matrix(rep(grp,repetitions),nrow(pred.prob),repetitions)
# calculate confidence intervals (and estimate value) for AUC
CIauc <- ci.cvAUC(predictions = pred.prob, labels =cilab, 
                      label.ordering = NULL, folds = NULL, confidence = 0.95)
# create data.frame with all means and CI
CImethod <- "bayes" # specify the method to use for CI
cvCI <- data.frame(array(NA,c(6,3)))
colnames(cvCI) <- c("average", "lowCI", "highCI")
rownames(cvCI) <- c("acc", "sen", "spe", "ppv", "npv", "auc")
cvCI[1,1:3] <- c(CIacc[CIacc$method==CImethod,4], min(CIacc[CIacc$method==CImethod,5]), max(CIacc[CIacc$method==CImethod,6]))
cvCI[2,1:3] <- c(min(CIsen[CIsen$method==CImethod,4]), min(CIsen[CIsen$method==CImethod,5]), max(CIsen[CIsen$method==CImethod,6]))
cvCI[3,1:3] <- c(min(CIspe[CIspe$method==CImethod,4]), min(CIspe[CIspe$method==CImethod,5]), max(CIspe[CIspe$method==CImethod,6]))
cvCI[4,1:3] <- c(min(CIppv[CIppv$method==CImethod,4]), min(CIppv[CIppv$method==CImethod,5]), max(CIppv[CIppv$method==CImethod,6]))
cvCI[5,1:3] <- c(min(CInpv[CInpv$method==CImethod,4]), min(CInpv[CInpv$method==CImethod,5]), max(CInpv[CInpv$method==CImethod,6]))
cvCI[6,1:3] <- c(CIauc$cvAUC, CIauc$ci[1],CIauc$ci[2])
# TABLE WITH RESULTS FROM AUC
round(cvCI*100,1)

# remove un-necessary temporary objects
rm(acc); rm(sen); rm(spe); rm(ppv); rm(npv); rm(A); rm(B); rm(C); rm(D); rm(CVresults)
rm(p); rm(n); rm(CIsen); rm(CIspe); rm(CIacc); rm(CIppv); rm(CInpv); rm(cilab); rm(CIauc)
rm(CImethod); rm(cvCI)
  

# Figure 4A ###################################################################
# Medians of the LD scores for each sample, calculated over the optimized models from the RDCV

X11(width = 6, height = 6)
LDdf <- data.frame(grp, apply(LD_scores, 1 , function(x) median(x,na.rm=T)))
colnames(LDdf) <- c("class", "LDscore")
levels(LDdf$class) <- c(neg.class, pos.class)
boxplot(LDdf$LDscore ~ LDdf$class, outline=F, boxwex=0.5,col=rgb(0,0,0,0),
          ylab="Median Linear Discriminant scores", xlab="Class",
          ylim = c(min(LDdf$LDscore), max(LDdf$LDscore)))
stripchart(LDdf[LDdf$class == neg.class,]$LDscore, add=T, vertical = T,
           pch=19, col=rgb(0,0,1,0.2), cex=1.5,
           method="jitter", jitter=0.1)
stripchart(LDdf[LDdf$class == pos.class,]$LDscore, add=T,vertical = T,
           pch=17, col=rgb(1,0,0,0.2), cex=1.5, at=2,
           method="jitter", jitter=0.1)
boxplot(LDdf$LDscore ~ LDdf$class, outline=F, boxwex=0.5, add=T, col=rgb(0,0,0,0))
abline(h=0, lty=2, col="darkgrey")
rm(LDdf)  
  
  
# Figure 4B ###################################################################
# ROC curves (B) of the optimized models from the RDCV. The average ROC is shown
# as non-transparent, black trace.
  
X11(width = 6, height = 6)
cols <- rgb(0,0,0,0.02)
plot(res.roc[,,1], type="l", col=cols, xlab="FPR", ylab="TPR", lwd=3)
for (p in c(2:dim(res.roc)[3])){
  lines(res.roc[,,p], type="l", col=rgb(0,0,0,0.05), lwd=2)
}
roc.median <- apply(res.roc,c(1,2),median, na.rm=TRUE)
lines(roc.median, lwd=3, lty=1, col="black")
segments(0,0,1,1, lty=2, col="black")
rm(cols); rm(p); rm(roc.median)


# Figure 5 ####################################################################
# Medians of the PCA scores for the first 4 principal components, grouped according
# to class, calculated over the optimized models from the RDCV; the significance 
# with respect to the Mann-Whitney U test for the 2 classes is reported for each component.
# NOTE: the scores of some PCs could be "reversed" with respect to the figure in the paper

# ALIGN VALUES OF PCA SCORES AND LOADINGS FROM DIFFERENT ITERATIONS (Folds/repetitions)
# create a vector with the sum of intensities (absolute) of 
# difference with respect to PC loadings from optimized model #1
diff.load <- t(array(data=NA, dim(PCAload[1,,])))
for (j in (1:20)){
  # j = PC
  for (i in (1:(repetitions*out.segments))){
    # i = model
    diff.load[i,j] <- sum(sqrt((PCAload[,j,i] - PCAload[,j,1])^4))
  }
}
# create a +1/-1 coefficient to multiply
PCAcoeff <- array(data=NA, c(repetitions*out.segments, 20) )
for (i in c(1:20)){
  PCAcoeff[match(diff.load[,i][diff.load[,i] > max(diff.load[,i])/2], diff.load[,i]),i] <- -1
  PCAcoeff[match(diff.load[,i][diff.load[,i] < max(diff.load[,i])/2], diff.load[,i]),i] <- 1
}
# align PCA scores
PCAscores_aligned <- array(data=NA, dim(PCAscores) )
for (i in c(1:20)){
  PCAscores_aligned[,i,] <- sweep (PCAscores[,i,], 2, PCAcoeff[,i], "*")
}
# calculate median scores
rdcv_summary <- LD_scores
rdcv_summary[!is.na(rdcv_summary)] <- 1
rdcv_summary[is.na(rdcv_summary)] <- 0

# calculate medians of PCA scores
PCAscores_temp <- array(NA, c(nrow(spc), 20, out.segments*repetitions))
for (i in (1:(out.segments*repetitions))){
  PCAscores_temp[rdcv_summary[,i] == 0,,i] <- na.omit(PCAscores_aligned[,,i])
}

PCAscores_med <- apply(PCAscores_temp, c(1,2), median, na.rm=TRUE)
maxoptPC <- 4 # specify maximum PC to include in the graph
PCAscores_df <- data.frame(grp, PCAscores_med[,1], rep( paste("PC",1, sep="") ,nrow(spc)))
colnames(PCAscores_df) <- c("class", "scores", "PC")
for (i in c(2:maxoptPC)){
  PCAscores_df_temp <- data.frame(grp, PCAscores_med[,i], rep(paste("PC",i, sep=""),nrow(spc)))
  colnames(PCAscores_df_temp) <- c("class", "scores", "PC")
  PCAscores_df <- rbind(PCAscores_df, PCAscores_df_temp)
}
levels(PCAscores_df$class) <- c(neg.class, pos.class)

require(ggsignif)
require(ggplot2)

# generate plot
X11(width = 12, height = 6)
p <- ggplot(PCAscores_df, aes(x=class, y=scores, color=factor(class), shape = factor(class))) + 
  geom_boxplot(fill="white", outlier.alpha = 0, color="black") +
  geom_jitter(position=position_jitter(0.09), cex=3) +
  geom_hline(yintercept=0, linetype="dashed", color = "darkgrey") +
  geom_signif(comparisons = list(c(neg.class, pos.class)), color = "black", test = wilcox.test,
              map_signif_level=c("***"=0.001, "**"=0.01, "*"=0.05, " "=2), 
              size=0.5, textsize = 4, margin_top = 0.08, vjust=0.4)
p + facet_wrap(. ~ PC , scales = 'free_y', nrow = 1) + 
  scale_color_manual(values=c(rgb(0,0,1,0.1), rgb(1,0,0,0.1))) + 
  scale_shape_manual(values=c(19, 17)) + 
  labs(y = "Median PCA scores") +
  theme(legend.position="none",
        panel.border = element_rect(colour = "black", fill=NA, size=0.5),
        panel.background = element_blank(),
        axis.text.x = element_text(color="black", size = 10),
        axis.text.y = element_text(color="black"),
        axis.title.y = element_text(vjust = 2.5),
        strip.background = element_rect(colour="black", fill="white", 
                                        size=0, linetype="solid"))

# remove un-necessary temporary objects
rm(diff.load); rm(PCAcoeff); rm(i); rm(j); rm(PCAscores_aligned); rm(PCAscores_temp)
rm(PCAscores_med); rm(maxoptPC); rm(PCAscores_df); rm(PCAscores_df_temp)


# Figure 6 ####################################################################
# Medians of the loadings for principal components 1,3 and 4, calculated over 
# the optimized models from the RDCV; interquartile ranges are reported in grey
# NOTE: the loadings of some PCs could be "reversed" with respect to the figure in the paper

# ALIGN VALUES OF PCA SCORES AND LOADINGS FROM DIFFERENT ITERATIONS (Folds/repetitions)
# create a vector with the sum of intensities (absolute) of 
# difference with respect to PC loadings from optimized model #1
diff.load <- t(array(data=NA, dim(PCAload[1,,])))
for (j in (1:20)){
  # j = PC
  for (i in (1:(repetitions*out.segments))){
    # i = model
    diff.load[i,j] <- sum(sqrt((PCAload[,j,i] - PCAload[,j,1])^4))
  }
}
# create a +1/-1 coefficient to multiply
PCAcoeff <- array(data=NA, c(repetitions*out.segments, 20) )
for (i in c(1:20)){
  PCAcoeff[match(diff.load[,i][diff.load[,i] > max(diff.load[,i])/2], diff.load[,i]),i] <- -1
  PCAcoeff[match(diff.load[,i][diff.load[,i] < max(diff.load[,i])/2], diff.load[,i]),i] <- 1
}
# align PCA loadlings 
PCAload_aligned <- array(data=NA, dim(PCAload) )
for (i in c(1:20)){
  PCAload_aligned[,i,] <- sweep (PCAload[,i,], 2, PCAcoeff[,i], "*")
}

X11(width = 6, height = 9)
nPC <- c(1,3,4) # specify number of PC to be plotted
stackoff <- 0.4 # stacking between two loadings
margins <- 0.08 # top and bottom margins
zerolines <- c(0:(length(nPC)-1))*stackoff
myloads <- apply(PCAload_aligned, c(1,2), median)[,nPC]
myloads_q3 <- apply(PCAload_aligned, c(1,2), function(x) quantile(x, 3/4))[,nPC] #quartile
myloads_q1 <- apply(PCAload_aligned, c(1,2), function(x) quantile(x, 1/4))[,nPC] #quartile
plot(spc@wavelength, myloads[,1], type="l", ylim=c(min(myloads[,1])- margins,
     zerolines[length(nPC)] + max(myloads[,length(nPC)] + margins)),
     xlab = expression(paste("Raman shift (", cm^-1, ")", sep = "")), 
     ylab = "Loadings of the RDCV models (a.u.)", yaxt="n")
segments(spc@wavelength,myloads_q1[,1], spc@wavelength,myloads_q3[,1], col="darkgrey", lwd=2)
lines(spc@wavelength, myloads[,1])
for (i in 2:length(nPC)){
  segments(spc@wavelength,myloads_q1[,i]+ zerolines[i], 
           spc@wavelength,myloads_q3[,i]+ zerolines[i], col="darkgrey", lwd=2)
  lines(spc@wavelength, myloads[,i]+ zerolines[i])
}
abline(h=zerolines, lty=2)
axis(side = 2, at = zerolines, labels = paste(rep("PC",length(nPC)), nPC, sep=""), las=2, cex.axis=0.8)
# label peaks
require(MALDIquant)
for (i in (1:length(nPC))){
  mspec1 <- createMassSpectrum(as.vector(spc@wavelength),as.vector(myloads[,i]))
  peaks1 <- detectPeaks(mspec1,halfWindowSize=11, method=c("MAD", "SuperSmoother"), SNR=0.8)
  peaks1@intensity <- peaks1@intensity + zerolines[i] 
  labelPeaks(peaks1, digits=0, underline=FALSE, srt=90,avoidOverlap=FALSE , adj=c(-0.2,0.4)) 
  mspec1 <- createMassSpectrum(as.vector(spc@wavelength),as.vector(-myloads[,i]))
  peaks1 <- detectPeaks(mspec1,halfWindowSize=11, method=c("MAD", "SuperSmoother"), SNR=0.8)
  labelPeaks(peaks1, absoluteVerticalPos = -peaks1@intensity + zerolines [i],
             digits=0, underline=FALSE, srt=90,avoidOverlap=FALSE , adj=c(1.2,0.4)) 
}

# remove un-necessary temporary objects
rm(i); rm(nPC); rm(stackoff); rm(margins); rm(zerolines); rm(myloads); rm(myloads_q1)
rm(myloads_q3); rm(diff.load); rm(PCAcoeff); rm(PCAload_aligned); rm(j)