## ----global_options, include=FALSE---------------------------------------
knitr::opts_chunk$set(fig.width=6, fig.height=6, fig.path='figures/')

## ----loading, echo=TRUE, warning=FALSE, message=FALSE--------------------
# loading the packages
library(proFIA)
library(plasFIA)

## ----paths---------------------------------------------------------------
# finding the directory of the raw files
path <- system.file(package="plasFIA", "mzML")
list.files(path)

## ----profiaset,fig.show="hide",message=FALSE-----------------------------
# defining the ppm parameter adapted to the Orbitrap Fusion
ppm <- 3

# performing the first step of the workflow
plasSet <- proFIAset(path, ppm=ppm, parallel=FALSE)

## ----raw_plot------------------------------------------------------------
# loading the spiked molecules data frame
data("plasMols")

# plotting the raw region aroung the Diphenhydramine mass signal
plasMols[7,]
mzrange <- c(plasMols[7,"mass_M+H"]-0.1,plasMols[7,"mass_M+H"]+0.1)
plotRaw(plasSet, type="r", sample=3, ylim=mzrange, size=0.6)

## ----peaks_plot----------------------------------------------------------
# plotting the filter Dipehnhydramine region.
plotRaw(plasSet, type="p", sample=3, ylim=mzrange, size=0.6)

## ----plot_injection------------------------------------------------------
# plotting the injection peak
plotInjectionPeaks(plasSet)

## ----group,message=FALSE-------------------------------------------------
# selecting the parameters
ppmgroup <- 1

# due to the experimental design, sample fraction was set to 0.2
fracGroup <- 0.2

# grouping
plasSet <- group.FIA(plasSet, ppm=ppmgroup, fracGroup=fracGroup)

## ----plotEICs------------------------------------------------------------
# plotting the EICs of the parameters.
plotEICs(plasSet,mz=plasMols[4,"mass_M+H"])

## ----find_group----------------------------------------------------------
# searching for match group with 2 ppm tolerance
lMatch <- findMzGroup(plasSet,plasMols[,"mass_M+H"],tol=3)

# index of the 40 molecules which may be used with plotEICs
molFound <- data.frame(names=plasMols[,"names"],found=lMatch)
head(molFound)

# getting the molecules which are not detected
plasMols[which(is.na(lMatch)),]

## ----datamatrix----------------------------------------------------------
# building the data matrix
plasSet <- makeDataMatrix(plasSet, maxo=FALSE)

## ----fillpeaks, warning=FALSE, eval=FALSE--------------------------------
## # since there is only 1 replicate, k is set to 1, which should not be the case in a standard experiment.
## k <- 1
## 
## # imputation
## plasSet <- fillPeaks.WKNN(plasSet, k=k)

## ----analyzeAcquisitionFIA, eval=FALSE-----------------------------------
## # selecting the parameters
## ppm <- 2
## ppmgroup <- 1
## fracGroup <- 0.2
## k <- 1
## 
## # running the whole workflow in a single step
## plasSet <- analyzeAcquisitionFIA(path, ppm=ppm, ppmgroup=ppmgroup, k=k,fracGroup = fracGroup,parallel=FALSE)

## ----export--------------------------------------------------------------
# Expression Set.
eset <- exportExpressionSet(plasSet)
eset

# Peak Table.
peakTable <- exportPeakTable(plasSet)


# 3 W4M Tables:
dataMatrix <- exportDataMatrix(plasSet)
sampleMetadata <- exportSampleMetadata(plasSet)
variableMetadata <- exportVariableMetadata(plasSet)

## ----conc----------------------------------------------------------------
# getting the concentrations of the spiked compounds (log10 scale)
data("plasSamples")
vconc <- log10(plasSamples[,"concentration_ng_ml"])

## ----opls, eval=FALSE----------------------------------------------------
## # building the OPLS model
## library(ropls)
## plasSet.opls <- opls(dataMatrix, vconc, predI = 1, orthoI = NA, log10L = TRUE, crossvalI=5)

## ----opls_eval, echo=FALSE-----------------------------------------------
# building the OPLS model
library(ropls)
plasSet.opls <- opls(t(dataMatrix), vconc, predI = 1, orthoI = NA, log10L = TRUE, crossvalI=5, plotL = FALSE)
layout(matrix(1:4, nrow = 2, byrow = TRUE))
for(typeC in c("permutation", "overview", "outlier", "x-score"))
plot(plasSet.opls, typeVc = typeC, parDevNewL = FALSE)

## ----matrix_effect_plot--------------------------------------------------
# matrix effec metric
MEmetricVn <- 1 - variableMetadata[, "corMean"]
nnaVl <- !is.na(MEmetricVn)
MEmetricVn <- MEmetricVn[nnaVl]
ordVi <- order(MEmetricVn)
MEmetricVn <- MEmetricVn[ordVi]

# predictive and orthogonal VIP
vipVn <- getVipVn(plasSet.opls)[nnaVl][ordVi]
orthoVipVn <- getVipVn(plasSet.opls, orthoL = TRUE)[nnaVl][ordVi]

# plot
colVc <- rev(rainbow(sum(nnaVl), end = 4/6))
plot(vipVn, orthoVipVn, pch = 16, col = colVc,
     xlab = "VIP_pred", ylab = "VIP_ortho", main = "VIP_ortho vs VIP_pred",lwd=3)

# adding black circles around spiked compounds
matchVl <- (1:nrow(variableMetadata) %in% lMatch[!is.na(lMatch)])[nnaVl][ordVi]
points(vipVn[matchVl], orthoVipVn[matchVl], cex=1.2,pch=1,col="black",lwd=2)
legend("topright", legend = c(paste0(c("ME max: ", "ME min: "), round(rev(range(MEmetricVn)), 2)),"Spiked molecules"), pch=c(16,16,1),col = c(rev(colVc[c(1, length(colVc))]),1))

## ----sessionInfo, echo=FALSE---------------------------------------------
sessionInfo()