inst/shiny-examples/MSstatsQCgui/QCMetrics.R
In eralpdogu/MSstatsQCgui: A graphical user interface for MSstatsQC package
#############################################################################################
#INPUTS : "prodata" is the data user uploads.
# "precursorSelection" is the precursor that user selects in Data Import tab. it can be either one precursor(peptide) or it can be "all peptides"
# "L" and "U" are lower and upper bound of guide set that user choose in Data Import tab.
# "metric" is one of these metrics: COL.BEST.RET,COL.FWHM, COL.TOTAL.AREA,COL.PEAK.ASS or a metric that user defines in his data set
# "normalization" is either TRUE or FALSE
#Description of function : it gets the metric column only for the precursor chosen and either return the column as it is or normalize it and then return it
getMetricData <- function(prodata, precursorSelection, L, U, metric, normalization,selectMean,selectSD, guidset_selected) {
precursor.data<-prodata[prodata$Precursor==precursorSelection,] #"Precursor" is one of the columns in data that shows the name of peptides. This gets only the part of data related to the selected peptide
metricData <- 0
if(is.null(metric)){
return(NULL)
}
metricData = precursor.data[,metric]
if(normalization == TRUE) {
if(guidset_selected) {
mu=mean(metricData[L:U]) # in-control process mean
sd=sd(metricData[L:U]) # in-control process variance
}else {
mu = selectMean
sd = selectSD
}
if(sd == 0) {sd <- 0.0001}
metricData=scale(metricData[seq_len(length(metricData))],mu,sd) # transformation for N(0,1) )
return(metricData)
} else if(normalization == FALSE){
return(metricData)
}
}
#########################################################################################################
find_custom_metrics <- function(prodata) {
prodata <- prodata[, which(colnames(prodata)=="Annotations"):ncol(prodata),drop = FALSE]
nums <- vapply(prodata, is.numeric, FUN.VALUE = logical(1))
other.metrics <- colnames(prodata[,nums])[1:ifelse(length(colnames(prodata[,nums]))<31,
length(colnames(prodata[,nums])),
30)
] # limiting custom metrics up to 30 metrics and not more
if(any(is.na(other.metrics))) {
return(c())
}
return(other.metrics)
}
################################################################
#z = metricData
#INPUT : "prodata" is the data user uploads.
# "metricData" is the column of the data related to the metric we want. Forexample if we want retention time, it gives retention time column
# "precursorSelection" is the precursor that user selects in Data Import tab. it can be either one precursor(peptide) or it can be "all peptides"
# "L" and "U" are lower and upper bound of guide set that user choose in Data Import tab.
# "type" is either 1 or 2. one is "Individual Value" plot and other "Moving Range" plot
#DESCRIPTION : returns a data frame for CUSUM that contains all the information needed to plot CUSUM
CUSUM.data.prepare <- function(prodata, metricData, precursorSelection, type) {
outRangeInRangePoz = rep(0,length(metricData))
outRangeInRangeNeg = rep(0,length(metricData))
#precursor.data only gets the data for the selected precursor
precursor.data <- prodata[prodata$Precursor==precursorSelection,] #"Precursor" is one of the columns in data that shows the name of peptides
v <- numeric(length(metricData))
Cpoz <- numeric(length(metricData))
Cneg <- numeric(length(metricData))
for(i in 2:length(metricData)) {
Cpoz[i]=max(0,(metricData[i]-0.5+Cpoz[i-1]))
Cneg[i]=max(0,((-0.5)-metricData[i]+Cneg[i-1]))
}
if(type==2) {
v[2:length(metricData)]=(sqrt(abs(metricData[2:length(metricData)]))-0.822)/0.349
for(i in 2:length(metricData)) {
Cpoz[i]=max(0,(v[i]-(0.5)+Cpoz[i-1]))
Cneg[i]=max(0,((-0.5)-v[i]+Cneg[i-1]))
}
}
QCno = seq_len(length(metricData))
for(i in seq_len(length(metricData))) {
if(Cpoz[i] >= 5 || Cpoz[i] <= -5)
outRangeInRangePoz[i] <- "OutRangeCUSUM+"
else
outRangeInRangePoz[i] <- "InRangeCUSUM+"
}
for(i in seq_len(length(metricData))) {
if(-Cneg[i] >= 5 || -Cneg[i] <= -5)
outRangeInRangeNeg[i] <- "OutRangeCUSUM-"
else
outRangeInRangeNeg[i] <- "InRangeCUSUM-"
}
plot.data =
data.frame(QCno = QCno
,CUSUM.poz = Cpoz
,CUSUM.neg = -Cneg
,Annotations=precursor.data$Annotations
,outRangeInRangePoz = outRangeInRangePoz
,outRangeInRangeNeg = outRangeInRangeNeg
)
return(plot.data)
}
###################################################################################################
#z = metricData
#INPUT : "prodata" is the data user uploads.
# "metricData" is the column of the data related to the metric we want. Forexample if we want retention time, it gives retention time column
# "type" is either 1 or 2. one is "Individual Value" plot and other "Moving Range" plot
#DESCRIPTION : returns a data frame for CP that contains all the information needed to plot Change Point
CP.data.prepare <- function(prodata, metricData, type) {
length_metricData <- length(metricData)
Et <- numeric(length_metricData-1) # this is Ct in type 1, and Dt in type 2.
SS<- numeric(length_metricData-1)
SST<- numeric(length_metricData-1)
tho.hat <- 0
if(type == 1) {
## Change point analysis for mean (Single step change model)
for(i in seq_len((length_metricData - 1))) {
Et[i]=(length_metricData-i)*(((1/(length_metricData-i))*sum(metricData[(i+1):length_metricData]))-0)^2 #change point function
}
QCno=seq_len(length_metricData - 1)
} else if(type == 2) {
## Change point analysis for variance (Single step change model)
SS[seq_len(length_metricData)]=metricData[seq_len(length_metricData)]^2
for(i in seq_len(length_metricData)) {
SST[i]=sum(SS[i:length_metricData])
Et[i]=((SST[i]/2)-((length_metricData-i+1)/2)*log(SST[i]/(length_metricData-i+1))-(length_metricData-i+1)/2) #change point function
}
QCno=seq_len(length_metricData)
}
tho.hat = which(Et==max(Et)) # change point estimate
plot.data <- data.frame(QCno,Et,tho.hat)
return(plot.data)
}
###################################################################################################
#INPUT : "prodata" is the data user uploads.
# "L" and "U" are lower and upper bound of guide set that user choose in Data Import tab.
# "data.metrics" is all the available metrics. It is defined in server.R
get_CP_tho.hat <- function(prodata, L, U, data.metrics,listMean,listSD, guidset_selected) {
tho.hat <- data.frame(tho.hat = c(), metric = c(), group = c(), y=c())
#precursors <- levels(reorder(prodata$Precursor,prodata[,COL.BEST.RET]))
precursors <- levels(prodata$Precursor)
for(metric in data.metrics) {
for (j in seq_len(nlevels(prodata$Precursor))) {
metricData <- getMetricData(prodata, precursors[j], L, U, metric = metric, normalization = TRUE,selectMean = listMean[[metric]],selectSD = listSD[[metric]], guidset_selected)
mix <- rbind(
data.frame(tho.hat = CP.data.prepare(prodata, metricData, type = 1)$tho.hat[1], metric = metric, group = "Individual Value", y=1.1),
data.frame(tho.hat = CP.data.prepare(prodata, metricData, type = 2)$tho.hat[1], metric = metric, group = "Moving Range", y=-1.1)
)
tho.hat <- rbind(tho.hat, mix)
}
}
return(tho.hat)
}
###################################################################################################
#z = metricData
#INPUT : "prodata" is the data user uploads.
# "metricData" is the column of the data related to the metric we want. Forexample if we want retention time, it gives retention time column
# "L" and "U" are lower and upper bound of guide set that user choose in Data Import tab.
# "type" is either 1 or 2. one is "Individual Value" plot and other "Moving Range" plot
#DESCRIPTION : returns a data frame for XmR that contains all the information needed to plot XmR
XmR.data.prepare <- function(prodata, metricData, L,U, type,selectMean,selectSD, guidset_selected) {
t <- numeric(length(metricData)-1)
UCL <- 0
LCL <- 0
InRangeOutRange <- rep(0,length(metricData))
t[2:length(metricData)]=abs(metricData[2:length(metricData)]-metricData[(2:length(metricData))-1]) # Compute moving range of metricData
QCno=seq_len(length(metricData))
if(type == 1) {
if(guidset_selected) {
UCL=mean(metricData[L:U])+2.66*sd(t[L:U])
LCL=mean(metricData[L:U])-2.66*sd(t[L:U])
}else {
UCL = selectMean + 2.66 * selectSD
LCL = selectMean - 2.66 * selectSD
}
t <- metricData
} else if(type == 2) {
## Calculate MR chart statistics and limits
if(guidset_selected) {
UCL=3.267*sd(t[seq_len(L)-U])
}else{
UCL = 3.267 * selectSD
}
LCL=0
}
for(i in seq_len(length(metricData))) {
if(t[i] > LCL && t[i] < UCL)
InRangeOutRange[i] <- "InRange"
else
InRangeOutRange[i] <- "OutRange"
}
plot.data=data.frame(QCno,metricData,t,UCL,LCL,InRangeOutRange)
return(plot.data)
}
############################################################################################
#INPUTS : "prodata" is the data user uploads.
# "metric" is one of these metrics: COL.BEST.RET,COL.FWHM, COL.TOTAL.AREA,COL.PEAK.ASS or a metric that user defines in his data set
# "L" and "U" are lower and upper bound of guide set that user choose in Data Import tab.
# "type" is either 1 or 2. one is "Individual Value" plot and other "Moving Range" plot
#DESCRIPTION : returns a data frame that is used in CUSUM.Summary.DataFrame function below to use for summary plot of CUSUM in Summary Tab of shiny app
CUSUM.Summary.prepare <- function(prodata, metric, L, U,type,selectMean,selectSD, guidset_selected) {
QCno <- seq_len(nrow(prodata))
y.poz <- rep(0,nrow(prodata))
y.neg <- rep(0,nrow(prodata))
counter <- rep(0,nrow(prodata))
#precursors <- levels(reorder(prodata$Precursor,prodata[,COL.BEST.RET]))
precursors <- levels(prodata$Precursor)
for(j in seq_len(length(precursors))) {
metricData <- getMetricData(prodata, precursors[j], L, U, metric = metric, normalization = T,selectMean,selectSD, guidset_selected)
counter[seq_len(length(metricData))] <- counter[seq_len(length(metricData))]+1
plot.data <- CUSUM.data.prepare(prodata, metricData, precursors[j], type)
sub.poz <- plot.data[plot.data$CUSUM.poz >= 5 | plot.data$CUSUM.poz <= -5, ]
sub.neg <- plot.data[plot.data$CUSUM.neg >= 5 | plot.data$CUSUM.neg <= -5, ]
y.poz[sub.poz$QCno] <- y.poz[sub.poz$QCno] + 1
y.neg[sub.neg$QCno] <- y.neg[sub.neg$QCno] + 1
}
max_QCno <- max(which(counter!=0))
pr.y.poz = y.poz[seq_len(max_QCno)]/counter[seq_len(max_QCno)]
pr.y.neg = y.neg[seq_len(max_QCno)]/counter[seq_len(max_QCno)]
plot.data <- data.frame(QCno = rep(seq_len(max_QCno),2),
pr.y = c(pr.y.poz, pr.y.neg),
group = ifelse(rep(type==1,2*max_QCno),
c(rep("Mean increase",max_QCno),
rep("Mean decrease",max_QCno)),
c(rep("Variability increase",max_QCno),
rep("Variability decrease",max_QCno))),
metric = rep(metric,max_QCno*2)
)
return(plot.data)
}
############################################################################################
CUSUM.Summary.DataFrame <- function(prodata, data.metrics, L, U,listMean,listSD, guidset_selected) {
dat <- data.frame(QCno = c(),
pr.y = c(),
group = c(),
metric = c())
for (metric in data.metrics) {
data.1 <- CUSUM.Summary.prepare(prodata, metric = metric, L, U,type = 1,selectMean = listMean[[metric]],selectSD = listSD[[metric]], guidset_selected)
data.2 <- CUSUM.Summary.prepare(prodata, metric = metric, L, U,type = 2,selectMean = listMean[[metric]],selectSD = listSD[[metric]], guidset_selected)
data.2$pr.y <- -(data.2$pr.y)
dat <- rbind(dat,data.1,data.2)
}
return(dat)
}
############################################################################################
#DESCRIPTION : for each metric returns a data frame of QCno, probability of out of control peptide for dispersion or mean plot
XmR.Summary.prepare <- function(prodata, metric, L, U,type,selectMean,selectSD, guidset_selected) {
QCno <- seq_len(nrow(prodata))
y.poz <- rep(0,nrow(prodata))
y.neg <- rep(0,nrow(prodata))
counter <- rep(0,nrow(prodata))
#precursors <- levels(reorder(prodata$Precursor,prodata[,COL.BEST.RET]))
precursors <- levels(prodata$Precursor)
for(j in seq_len(length(precursors))) {
metricData <- getMetricData(prodata, precursors[j], L = L, U = U, metric = metric, normalization = T,selectMean,selectSD, guidset_selected)
counter[seq_len(length(metricData))] <- counter[seq_len(length(metricData))]+1
plot.data <- XmR.data.prepare(prodata, metricData , L , U , type, selectMean,selectSD, guidset_selected)
sub.poz <- plot.data[plot.data$t >= plot.data$UCL, ]
sub.neg <- plot.data[plot.data$t <= plot.data$LCL, ]
y.poz[sub.poz$QCno] <- y.poz[sub.poz$QCno] + 1
y.neg[sub.neg$QCno] <- y.neg[sub.neg$QCno] + 1
}
max_QCno <- max(which(counter!=0))
pr.y.poz = y.poz[seq_len(max_QCno)]/counter[seq_len(max_QCno)]
pr.y.neg = y.neg[seq_len(max_QCno)]/counter[seq_len(max_QCno)]
plot.data <- data.frame(QCno = rep(seq_len(max_QCno),2),
pr.y = c(pr.y.poz, pr.y.neg),
group = ifelse(rep(type==1,2*max_QCno),
c(rep("Mean increase",max_QCno),
rep("Mean decrease",max_QCno)),
c(rep("Variability increase",max_QCno),
rep("Variability decrease",max_QCno))),
metric = rep(metric,max_QCno*2))
return(plot.data)
}
###########################################################################################
#DESCRIPTION : returns a data frame for all the metrics, probability of out of range peptides, wheter it is for metric mean or metric dispersion and i is an increase or decrease
XmR.Summary.DataFrame <- function(prodata, data.metrics, L, U,listMean, listSD, guidset_selected) {
dat <- data.frame(QCno = c(),
pr.y = c(),
group = c(),
metric = c())
for (metric in data.metrics) {
data.1 <- XmR.Summary.prepare(prodata, metric = metric, L, U,type = 1,selectMean = listMean[[metric]],selectSD = listSD[[metric]], guidset_selected)
data.2 <- XmR.Summary.prepare(prodata, metric = metric, L, U,type = 2,selectMean = listMean[[metric]],selectSD = listSD[[metric]], guidset_selected)
data.2$pr.y <- -(data.2$pr.y)
dat <- rbind(dat, data.1, data.2)
}
return(dat)
}
############################################################################################
heatmap.DataFrame <- function(prodata, data.metrics,method,peptideThresholdRed,peptideThresholdYellow,
L, U, type,listMean, listSD, guidset_selected) {
time <- c()
val <- c()
met <- c()
flag <- c()
for (metric in data.metrics) {
df <- Decision.DataFrame.prepare(prodata, metric, method, peptideThresholdRed,peptideThresholdYellow,
L, U,type,selectMean=listMean[[metric]],selectSD=listSD[[metric]], guidset_selected)
time_df <- as.character(df$AcquiredTime)
val_df <- df$pr.y
met_df <- rep(metric,length(val_df))
flag_df <- df$flag
time <- c(time,time_df)
val <- c(val,val_df)
met <- c(met,met_df)
flag <- c(flag,flag_df)
}
dataFrame <- data.frame(time = time,
value = val,
metric = met,
flag = flag,
QCno = df$QCno
)
return(dataFrame)
}
###############################################################################################################
XmR.Radar.Plot.prepare <- function(prodata,L,U, metric, type,group,XmR_type,selectMean,selectSD, guidset_selected) {
precursors <- levels(prodata$Precursor)
precursors2 <- substring(precursors, first = 1, last = 10)
QCno.length <- c()
QCno.out.range.poz <- c()
QCno.out.range.neg <- c()
for(j in seq_len(length(precursors))) {
metricData <- getMetricData(prodata, precursors[j], L = L, U = U, metric = metric, normalization = T,selectMean,selectSD, guidset_selected)
QCno.length <- c(QCno.length,length(metricData))
plot.data <- XmR.data.prepare(prodata, metricData , L = L, U = U, type ,selectMean,selectSD, guidset_selected)
QCno.out.range.poz <- c(QCno.out.range.poz,length(plot.data[plot.data$t >= plot.data$UCL, ]$QCno))
QCno.out.range.neg <- c(QCno.out.range.neg,length(plot.data[plot.data$t <= plot.data$LCL, ]$QCno))
}
#pr <- Compute.QCno.OutOfRangePeptide.XmR(prodata,L,U,metric = metric,type = type, XmR.type,selectMean,selectSD, guidset_selected)
if(XmR_type == "poz") {
dat <- data.frame(peptides = precursors2,
OutRangeQCno = QCno.out.range.poz,
group = rep(group,length(precursors)),
orderby = seq(seq_len(length(precursors))),
metric = rep(metric, length(precursors)),
tool = rep("XmR",length(precursors)),
probability = QCno.out.range.poz/QCno.length
)
} else {
dat <- data.frame(peptides = precursors2,
OutRangeQCno = QCno.out.range.neg,
group = rep(group,length(precursors)),
orderby = seq(seq_len(length(precursors))),
metric = rep(metric, length(precursors)),
tool = rep("XmR",length(precursors)),
probability = QCno.out.range.neg/QCno.length
)
}
return(dat)
#return(list(dat.poz,dat.neg))
}
################################################################################################
XmR.Radar.Plot.DataFrame <- function(prodata, data.metrics, L,U,listMean,listSD,guidset_selected) {
dat <- data.frame(peptides = c(), OutRangeQCno = c(), group = c(),
orderby = c(), metric = c(), tool = c(),
probability = c()
)
for (metric in data.metrics) {
data.1 <- XmR.Radar.Plot.prepare(prodata,L,U,metric = metric, type = 1,group = "Mean increase",XmR_type = "poz",
selectMean = listMean[[metric]],selectSD = listSD[[metric]],guidset_selected)
data.2 <- XmR.Radar.Plot.prepare(prodata,L,U,metric = metric, type = 1,group = "Mean decrease",XmR_type = "neg",
selectMean = listMean[[metric]],selectSD = listSD[[metric]],guidset_selected)
data.3 <- XmR.Radar.Plot.prepare(prodata,L,U,metric = metric, type = 2,group = "Variability increase",XmR_type = "poz",
selectMean = listMean[[metric]],selectSD = listSD[[metric]],guidset_selected)
data.4 <- XmR.Radar.Plot.prepare(prodata,L,U,metric = metric, type = 2,group = "Variability decrease",XmR_type = "neg",
selectMean = listMean[[metric]],selectSD = listSD[[metric]],guidset_selected)
dat <- rbind(dat, data.1, data.2, data.3, data.4)
}
return(dat)
}
#################################################################################################################
CUSUM.Radar.Plot.prepare <- function(prodata,L,U, metric,type,group, CUSUM.type,selectMean,selectSD, guidset_selected) {
#precursors <- levels(reorder(prodata$Precursor,prodata[,COL.BEST.RET]))
precursors <- levels(prodata$Precursor)
precursors2 <- substring(precursors, first = 1, last = 10)
QCno.length <- c()
QCno.out.range <- c()
for(j in seq_len(length(precursors))) {
metricData <- getMetricData(prodata, precursors[j], L = L, U = U, metric = metric, normalization = T,selectMean,selectSD, guidset_selected)
QCno.length <- c(QCno.length,length(metricData))
plot.data <- CUSUM.data.prepare(prodata, metricData, precursors[j], type)
if(CUSUM.type == "poz")
QCno.out.range <- c(QCno.out.range,length(plot.data[plot.data$CUSUM.poz >= 5 | plot.data$CUSUM.poz <= -5, ]$QCno))
else
QCno.out.range <- c(QCno.out.range,length(plot.data[plot.data$CUSUM.neg >= 5 | plot.data$CUSUM.neg <= -5, ]$QCno))
}
#pr <- Compute.QCno.OutOfRangePeptide.CUSUM(prodata,L,U,metric = metric,type = type, CUSUM.type,selectMean,selectSD, guidset_selected)
dat <- data.frame(peptides = precursors2,
OutRangeQCno = QCno.out.range,
group = rep(group,length(precursors)),
orderby = seq(seq_len(length(precursors))),
metric = rep(metric, length(precursors)),
tool = rep("XmR",length(precursors)),
probability = QCno.out.range/QCno.length
)
return(dat)
}
#################################################################################################
CUSUM.Radar.Plot.DataFrame <- function(prodata, data.metrics, L,U,listMean,listSD,guidset_selected) {
dat <- data.frame(peptides = c(), OutRangeQCno = c(), group = c(),
orderby = c(), metric = c(), tool = c(),
probability = c()
)
for (metric in data.metrics) {
data.1 <- CUSUM.Radar.Plot.prepare(prodata,L,U, metric = metric, type = 1, group = "Mean increase", CUSUM.type = "poz",selectMean = listMean[[metric]],selectSD = listSD[[metric]],guidset_selected)
data.2 <- CUSUM.Radar.Plot.prepare(prodata,L,U, metric = metric, type = 1, group = "Mean decrease", CUSUM.type = "neg",selectMean = listMean[[metric]],selectSD = listSD[[metric]],guidset_selected)
data.3 <- CUSUM.Radar.Plot.prepare(prodata,L,U, metric = metric, type = 2, group = "Variability increase", CUSUM.type = "poz",selectMean = listMean[[metric]],selectSD = listSD[[metric]],guidset_selected)
data.4 <- CUSUM.Radar.Plot.prepare(prodata,L,U, metric = metric, type = 2, group = "Variability decrease", CUSUM.type = "neg",selectMean = listMean[[metric]],selectSD = listSD[[metric]],guidset_selected)
dat <- rbind(dat, data.1, data.2, data.3, data.4)
}
return(dat)
}
#######################################################################################################
Decision.DataFrame.prepare <- function(prodata, metric, method, peptideThresholdRed,peptideThresholdYellow,
L, U,type,selectMean,selectSD, guidset_selected) {
AcquiredTime <- prodata$AcquiredTime
QCno <- seq_len(nrow(prodata))
y <- rep(0,nrow(prodata))
counter <- rep(0,nrow(prodata))
#precursors <- levels(reorder(prodata$Precursor,prodata[,COL.BEST.RET]))
precursors <- levels(prodata$Precursor)
if(method == "XmR") {
for(precursor in precursors) {
metricData <- getMetricData(prodata, precursor, L = L, U = U, metric = metric, normalization = T,selectMean,selectSD, guidset_selected)
counter[seq_len(length(metricData))] <- counter[seq_len(length(metricData))]+1
plot.data <- XmR.data.prepare(prodata, metricData , L , U , type,selectMean,selectSD, guidset_selected)
sub <- plot.data[plot.data$InRangeOutRange == "OutRange",]
y[sub$QCno] <- y[sub$QCno] + 1
}
} else if(method == "CUSUM") {
for(precursor in precursors) {
metricData <- getMetricData(prodata, precursor, L = L, U = U, metric = metric, normalization = T,selectMean,selectSD, guidset_selected)
counter[seq_len(length(metricData))] <- counter[seq_len(length(metricData))]+1
plot.data <- CUSUM.data.prepare(prodata, metricData, precursor, type)
sub <- plot.data[(plot.data$CUSUM.poz >= 5 | plot.data$CUSUM.poz <= -5) | (plot.data$CUSUM.neg >= 5 | plot.data$CUSUM.neg <= -5), ]
y[sub$QCno] <- y[sub$QCno] + 1
}
}
max_QCno <- max(which(counter!=0))
pr.y = y[seq_len(max_QCno)]/counter[seq_len(max_QCno)]
plot.data <- data.frame(AcquiredTime = AcquiredTime[seq_len(max_QCno)],
QCno = rep(seq_len(max_QCno),1),
pr.y = pr.y,
group = ifelse(rep(type==1,max_QCno),
rep("Metric mean",max_QCno),
rep("Metric dispersion",max_QCno)
),
metric = rep(metric,max_QCno),
flag = rep(0,max_QCno)
)
for (i in seq_len(max_QCno)) {
if(plot.data$pr.y[i] > peptideThresholdRed){
plot.data$flag[i] <- "Fail"
}
else if(plot.data$pr.y[i] > peptideThresholdYellow){
plot.data$flag[i] <- "Warning"
}
else {
plot.data$flag[i] <- "Pass"
}
}
if(type == 2) {
return(plot.data[-1,])
}
return(plot.data)
}
#######################################################################################################
number.Of.Out.Of.Range.Metrics <- function(prodata,data.metrics,method, peptideThresholdRed,
peptideThresholdYellow, L, U, type,listMean,listSD,
guidset_selected) {
metricCounterAboveRed = 0
metricCounterAboveYellowBelowRed = 0
precursors <- levels(prodata$Precursor)
for (metric in data.metrics) {
QCno <- seq_len(nrow(prodata))
y <- rep(0,nrow(prodata))
counter <- rep(0,nrow(prodata))
for(precursor in precursors) {
metricData <- getMetricData(prodata, precursor, L = L, U = U, metric = metric, normalization = T,listMean[[metric]],listSD[[metric]], guidset_selected)
counter[seq_len(length(metricData))] <- counter[seq_len(length(metricData))]+1
plot.data <- XmR.data.prepare(prodata, metricData , L , U , type,selectMean = listMean[[metric]],selectSD = listSD[[metric]], guidset_selected)
sub <- plot.data[plot.data$InRangeOutRange == "OutRange",]
y[sub$QCno] <- y[sub$QCno] + 1
}
max_QCno <- max(which(counter!=0))
pr.y = y[seq_len(max_QCno)]/counter[seq_len(max_QCno)]
if(type == 2) {
pr.y <- pr.y[-1]
}
aboveYellow <- which(pr.y > peptideThresholdYellow)
aboveYellowBelowRed <- which(pr.y > peptideThresholdYellow & pr.y <= peptideThresholdRed)
if(length(which(pr.y > peptideThresholdRed)) > 0) {
metricCounterAboveRed = metricCounterAboveRed + 1
}
if(length(aboveYellowBelowRed) > 0) {
metricCounterAboveYellowBelowRed = metricCounterAboveYellowBelowRed + 1
}
}
return(c(metricCounterAboveRed,metricCounterAboveYellowBelowRed))
}
####################################################################################################
decisionRule_warning_message_XmR <- function(prodata, data.metrics, method, peptideThresholdRed, peptideThresholdYellow,metricThresholdRed,metricThresholdYellow, L,U, type, listMean,listSD , guidset_selected) {
a1 <- number.Of.Out.Of.Range.Metrics(prodata,data.metrics,method = "XmR", peptideThresholdRed,peptideThresholdYellow,
L, U, type = 1,listMean,listSD, guidset_selected)
XmRCounterAboveRed1 <- a1[1]
XmRCounterAboveYellow1 <- a1[2]
a2 <- number.Of.Out.Of.Range.Metrics(prodata,data.metrics,method = "XmR", peptideThresholdRed,peptideThresholdYellow,
L, U, type = 2,listMean,listSD, guidset_selected)
XmRCounterAboveRed2 <- a2[1]
XmRCounterAboveYellow2 <- a2[2]
if(XmRCounterAboveRed1 > metricThresholdRed && XmRCounterAboveRed2 > metricThresholdRed) return({paste("XmR RED FLAG: System performance is UNACCEPTABLE")})
if(XmRCounterAboveRed1 > metricThresholdRed && XmRCounterAboveRed2 <= metricThresholdRed) return({paste("XmR RED FLAG: System performance is UNACCEPTABLE")})
if(XmRCounterAboveRed1 <= metricThresholdRed && XmRCounterAboveRed2 > metricThresholdRed) return({paste("XmR RED FLAG: System performance is UNACCEPTABLE")})
if(XmRCounterAboveYellow1 > metricThresholdYellow && XmRCounterAboveYellow2 > metricThresholdYellow) return({"XmR Yellow FLAG: System performance is POOR"})
if(XmRCounterAboveYellow1 <= metricThresholdYellow && XmRCounterAboveYellow2 > metricThresholdYellow) return({"XmR Yellow FLAG: System performance is POOR"})
if(XmRCounterAboveYellow1 > metricThresholdYellow && XmRCounterAboveYellow2 <= metricThresholdYellow) return({"XmR Yellow FLAG: System performance is POOR"})
return({"XmR BLUE FLAG: System performance is acceptable"})
}
############################################################################################################
decisionRule_warning_message_CUSUM <- function(prodata, data.metrics, method, peptideThresholdRed, peptideThresholdYellow,metricThresholdRed,metricThresholdYellow, L,U, type, listMean,listSD, guidset_selected ) {
a1 <- number.Of.Out.Of.Range.Metrics(prodata,data.metrics,method = "CUSUM", peptideThresholdRed,peptideThresholdYellow,
L, U, type = 1,listMean,listSD, guidset_selected)
CUSUMCounterAboveRed1 <- a1[1]
CUSUMCounterAboveYellow1 <- a1[2]
a2 <- number.Of.Out.Of.Range.Metrics(prodata,data.metrics,method = "CUSUM", peptideThresholdRed,peptideThresholdYellow,
L, U, type = 2,listMean,listSD, guidset_selected)
CUSUMCounterAboveRed2 <- a2[1]
CUSUMCounterAboveYellow2 <-a2[2]
if(CUSUMCounterAboveRed1 > metricThresholdRed && CUSUMCounterAboveRed2 > metricThresholdRed) return({paste("CUSUM RED FLAG: System performance is UNACCEPTABLE")})
if(CUSUMCounterAboveRed1 > metricThresholdRed && CUSUMCounterAboveRed2 <= metricThresholdRed) return({paste("CUSUM RED FLAG: System performance is UNACCEPTABLE")})
if(CUSUMCounterAboveRed1 <= metricThresholdRed && CUSUMCounterAboveRed2 > metricThresholdRed) return({paste("CUSUM RED FLAG: System performance is UNACCEPTABLE")})
if(CUSUMCounterAboveYellow1 > metricThresholdYellow && CUSUMCounterAboveYellow2 > metricThresholdYellow) return({"CUSUM Yellow FLAG: System performance is POOR"})
if(CUSUMCounterAboveYellow1 <= metricThresholdYellow && CUSUMCounterAboveYellow2 > metricThresholdYellow) return({"CUSUM Yellow FLAG: System performance is POOR"})
if(CUSUMCounterAboveYellow1 > metricThresholdYellow && CUSUMCounterAboveYellow2 <= metricThresholdYellow) return({"CUSUM Yellow FLAG: System performance is POOR"})
return({"CUSUM BLUE FLAG: System performance is acceptable"})
}
eralpdogu/MSstatsQCgui documentation built on May 26, 2020, 12:18 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#############################################################################################
#INPUTS : "prodata" is the data user uploads.
# "precursorSelection" is the precursor that user selects in Data Import tab. it can be either one precursor(peptide) or it can be "all peptides"
# "L" and "U" are lower and upper bound of guide set that user choose in Data Import tab.
# "metric" is one of these metrics: COL.BEST.RET,COL.FWHM, COL.TOTAL.AREA,COL.PEAK.ASS or a metric that user defines in his data set
# "normalization" is either TRUE or FALSE
#Description of function : it gets the metric column only for the precursor chosen and either return the column as it is or normalize it and then return it
getMetricData <- function(prodata, precursorSelection, L, U, metric, normalization,selectMean,selectSD, guidset_selected) {
precursor.data<-prodata[prodata$Precursor==precursorSelection,] #"Precursor" is one of the columns in data that shows the name of peptides. This gets only the part of data related to the selected peptide
metricData <- 0
if(is.null(metric)){
return(NULL)
}
metricData = precursor.data[,metric]
if(normalization == TRUE) {
if(guidset_selected) {
mu=mean(metricData[L:U]) # in-control process mean
sd=sd(metricData[L:U]) # in-control process variance
}else {
mu = selectMean
sd = selectSD
}
if(sd == 0) {sd <- 0.0001}
metricData=scale(metricData[seq_len(length(metricData))],mu,sd) # transformation for N(0,1) )
return(metricData)
} else if(normalization == FALSE){
return(metricData)
}
}
#########################################################################################################
find_custom_metrics <- function(prodata) {
prodata <- prodata[, which(colnames(prodata)=="Annotations"):ncol(prodata),drop = FALSE]
nums <- vapply(prodata, is.numeric, FUN.VALUE = logical(1))
other.metrics <- colnames(prodata[,nums])[1:ifelse(length(colnames(prodata[,nums]))<31,
length(colnames(prodata[,nums])),
30)
] # limiting custom metrics up to 30 metrics and not more
if(any(is.na(other.metrics))) {
return(c())
}
return(other.metrics)
}
################################################################
#z = metricData
#INPUT : "prodata" is the data user uploads.
# "metricData" is the column of the data related to the metric we want. Forexample if we want retention time, it gives retention time column
# "precursorSelection" is the precursor that user selects in Data Import tab. it can be either one precursor(peptide) or it can be "all peptides"
# "L" and "U" are lower and upper bound of guide set that user choose in Data Import tab.
# "type" is either 1 or 2. one is "Individual Value" plot and other "Moving Range" plot
#DESCRIPTION : returns a data frame for CUSUM that contains all the information needed to plot CUSUM
CUSUM.data.prepare <- function(prodata, metricData, precursorSelection, type) {
outRangeInRangePoz = rep(0,length(metricData))
outRangeInRangeNeg = rep(0,length(metricData))
#precursor.data only gets the data for the selected precursor
precursor.data <- prodata[prodata$Precursor==precursorSelection,] #"Precursor" is one of the columns in data that shows the name of peptides
v <- numeric(length(metricData))
Cpoz <- numeric(length(metricData))
Cneg <- numeric(length(metricData))
for(i in 2:length(metricData)) {
Cpoz[i]=max(0,(metricData[i]-0.5+Cpoz[i-1]))
Cneg[i]=max(0,((-0.5)-metricData[i]+Cneg[i-1]))
}
if(type==2) {
v[2:length(metricData)]=(sqrt(abs(metricData[2:length(metricData)]))-0.822)/0.349
for(i in 2:length(metricData)) {
Cpoz[i]=max(0,(v[i]-(0.5)+Cpoz[i-1]))
Cneg[i]=max(0,((-0.5)-v[i]+Cneg[i-1]))
}
}
QCno = seq_len(length(metricData))
for(i in seq_len(length(metricData))) {
if(Cpoz[i] >= 5 || Cpoz[i] <= -5)
outRangeInRangePoz[i] <- "OutRangeCUSUM+"
else
outRangeInRangePoz[i] <- "InRangeCUSUM+"
}
for(i in seq_len(length(metricData))) {
if(-Cneg[i] >= 5 || -Cneg[i] <= -5)
outRangeInRangeNeg[i] <- "OutRangeCUSUM-"
else
outRangeInRangeNeg[i] <- "InRangeCUSUM-"
}
plot.data =
data.frame(QCno = QCno
,CUSUM.poz = Cpoz
,CUSUM.neg = -Cneg
,Annotations=precursor.data$Annotations
,outRangeInRangePoz = outRangeInRangePoz
,outRangeInRangeNeg = outRangeInRangeNeg
)
return(plot.data)
}
###################################################################################################
#z = metricData
#INPUT : "prodata" is the data user uploads.
# "metricData" is the column of the data related to the metric we want. Forexample if we want retention time, it gives retention time column
# "type" is either 1 or 2. one is "Individual Value" plot and other "Moving Range" plot
#DESCRIPTION : returns a data frame for CP that contains all the information needed to plot Change Point
CP.data.prepare <- function(prodata, metricData, type) {
length_metricData <- length(metricData)
Et <- numeric(length_metricData-1) # this is Ct in type 1, and Dt in type 2.
SS<- numeric(length_metricData-1)
SST<- numeric(length_metricData-1)
tho.hat <- 0
if(type == 1) {
## Change point analysis for mean (Single step change model)
for(i in seq_len((length_metricData - 1))) {
Et[i]=(length_metricData-i)*(((1/(length_metricData-i))*sum(metricData[(i+1):length_metricData]))-0)^2 #change point function
}
QCno=seq_len(length_metricData - 1)
} else if(type == 2) {
## Change point analysis for variance (Single step change model)
SS[seq_len(length_metricData)]=metricData[seq_len(length_metricData)]^2
for(i in seq_len(length_metricData)) {
SST[i]=sum(SS[i:length_metricData])
Et[i]=((SST[i]/2)-((length_metricData-i+1)/2)*log(SST[i]/(length_metricData-i+1))-(length_metricData-i+1)/2) #change point function
}
QCno=seq_len(length_metricData)
}
tho.hat = which(Et==max(Et)) # change point estimate
plot.data <- data.frame(QCno,Et,tho.hat)
return(plot.data)
}
###################################################################################################
#INPUT : "prodata" is the data user uploads.
# "L" and "U" are lower and upper bound of guide set that user choose in Data Import tab.
# "data.metrics" is all the available metrics. It is defined in server.R
get_CP_tho.hat <- function(prodata, L, U, data.metrics,listMean,listSD, guidset_selected) {
tho.hat <- data.frame(tho.hat = c(), metric = c(), group = c(), y=c())
#precursors <- levels(reorder(prodata$Precursor,prodata[,COL.BEST.RET]))
precursors <- levels(prodata$Precursor)
for(metric in data.metrics) {
for (j in seq_len(nlevels(prodata$Precursor))) {
metricData <- getMetricData(prodata, precursors[j], L, U, metric = metric, normalization = TRUE,selectMean = listMean[[metric]],selectSD = listSD[[metric]], guidset_selected)
mix <- rbind(
data.frame(tho.hat = CP.data.prepare(prodata, metricData, type = 1)$tho.hat[1], metric = metric, group = "Individual Value", y=1.1),
data.frame(tho.hat = CP.data.prepare(prodata, metricData, type = 2)$tho.hat[1], metric = metric, group = "Moving Range", y=-1.1)
)
tho.hat <- rbind(tho.hat, mix)
}
}
return(tho.hat)
}
###################################################################################################
#z = metricData
#INPUT : "prodata" is the data user uploads.
# "metricData" is the column of the data related to the metric we want. Forexample if we want retention time, it gives retention time column
# "L" and "U" are lower and upper bound of guide set that user choose in Data Import tab.
# "type" is either 1 or 2. one is "Individual Value" plot and other "Moving Range" plot
#DESCRIPTION : returns a data frame for XmR that contains all the information needed to plot XmR
XmR.data.prepare <- function(prodata, metricData, L,U, type,selectMean,selectSD, guidset_selected) {
t <- numeric(length(metricData)-1)
UCL <- 0
LCL <- 0
InRangeOutRange <- rep(0,length(metricData))
t[2:length(metricData)]=abs(metricData[2:length(metricData)]-metricData[(2:length(metricData))-1]) # Compute moving range of metricData
QCno=seq_len(length(metricData))
if(type == 1) {
if(guidset_selected) {
UCL=mean(metricData[L:U])+2.66*sd(t[L:U])
LCL=mean(metricData[L:U])-2.66*sd(t[L:U])
}else {
UCL = selectMean + 2.66 * selectSD
LCL = selectMean - 2.66 * selectSD
}
t <- metricData
} else if(type == 2) {
## Calculate MR chart statistics and limits
if(guidset_selected) {
UCL=3.267*sd(t[seq_len(L)-U])
}else{
UCL = 3.267 * selectSD
}
LCL=0
}
for(i in seq_len(length(metricData))) {
if(t[i] > LCL && t[i] < UCL)
InRangeOutRange[i] <- "InRange"
else
InRangeOutRange[i] <- "OutRange"
}
plot.data=data.frame(QCno,metricData,t,UCL,LCL,InRangeOutRange)
return(plot.data)
}
############################################################################################
#INPUTS : "prodata" is the data user uploads.
# "metric" is one of these metrics: COL.BEST.RET,COL.FWHM, COL.TOTAL.AREA,COL.PEAK.ASS or a metric that user defines in his data set
# "L" and "U" are lower and upper bound of guide set that user choose in Data Import tab.
# "type" is either 1 or 2. one is "Individual Value" plot and other "Moving Range" plot
#DESCRIPTION : returns a data frame that is used in CUSUM.Summary.DataFrame function below to use for summary plot of CUSUM in Summary Tab of shiny app
CUSUM.Summary.prepare <- function(prodata, metric, L, U,type,selectMean,selectSD, guidset_selected) {
QCno <- seq_len(nrow(prodata))
y.poz <- rep(0,nrow(prodata))
y.neg <- rep(0,nrow(prodata))
counter <- rep(0,nrow(prodata))
#precursors <- levels(reorder(prodata$Precursor,prodata[,COL.BEST.RET]))
precursors <- levels(prodata$Precursor)
for(j in seq_len(length(precursors))) {
metricData <- getMetricData(prodata, precursors[j], L, U, metric = metric, normalization = T,selectMean,selectSD, guidset_selected)
counter[seq_len(length(metricData))] <- counter[seq_len(length(metricData))]+1
plot.data <- CUSUM.data.prepare(prodata, metricData, precursors[j], type)
sub.poz <- plot.data[plot.data$CUSUM.poz >= 5 | plot.data$CUSUM.poz <= -5, ]
sub.neg <- plot.data[plot.data$CUSUM.neg >= 5 | plot.data$CUSUM.neg <= -5, ]
y.poz[sub.poz$QCno] <- y.poz[sub.poz$QCno] + 1
y.neg[sub.neg$QCno] <- y.neg[sub.neg$QCno] + 1
}
max_QCno <- max(which(counter!=0))
pr.y.poz = y.poz[seq_len(max_QCno)]/counter[seq_len(max_QCno)]
pr.y.neg = y.neg[seq_len(max_QCno)]/counter[seq_len(max_QCno)]
plot.data <- data.frame(QCno = rep(seq_len(max_QCno),2),
pr.y = c(pr.y.poz, pr.y.neg),
group = ifelse(rep(type==1,2*max_QCno),
c(rep("Mean increase",max_QCno),
rep("Mean decrease",max_QCno)),
c(rep("Variability increase",max_QCno),
rep("Variability decrease",max_QCno))),
metric = rep(metric,max_QCno*2)
)
return(plot.data)
}
############################################################################################
CUSUM.Summary.DataFrame <- function(prodata, data.metrics, L, U,listMean,listSD, guidset_selected) {
dat <- data.frame(QCno = c(),
pr.y = c(),
group = c(),
metric = c())
for (metric in data.metrics) {
data.1 <- CUSUM.Summary.prepare(prodata, metric = metric, L, U,type = 1,selectMean = listMean[[metric]],selectSD = listSD[[metric]], guidset_selected)
data.2 <- CUSUM.Summary.prepare(prodata, metric = metric, L, U,type = 2,selectMean = listMean[[metric]],selectSD = listSD[[metric]], guidset_selected)
data.2$pr.y <- -(data.2$pr.y)
dat <- rbind(dat,data.1,data.2)
}
return(dat)
}
############################################################################################
#DESCRIPTION : for each metric returns a data frame of QCno, probability of out of control peptide for dispersion or mean plot
XmR.Summary.prepare <- function(prodata, metric, L, U,type,selectMean,selectSD, guidset_selected) {
QCno <- seq_len(nrow(prodata))
y.poz <- rep(0,nrow(prodata))
y.neg <- rep(0,nrow(prodata))
counter <- rep(0,nrow(prodata))
#precursors <- levels(reorder(prodata$Precursor,prodata[,COL.BEST.RET]))
precursors <- levels(prodata$Precursor)
for(j in seq_len(length(precursors))) {
metricData <- getMetricData(prodata, precursors[j], L = L, U = U, metric = metric, normalization = T,selectMean,selectSD, guidset_selected)
counter[seq_len(length(metricData))] <- counter[seq_len(length(metricData))]+1
plot.data <- XmR.data.prepare(prodata, metricData , L , U , type, selectMean,selectSD, guidset_selected)
sub.poz <- plot.data[plot.data$t >= plot.data$UCL, ]
sub.neg <- plot.data[plot.data$t <= plot.data$LCL, ]
y.poz[sub.poz$QCno] <- y.poz[sub.poz$QCno] + 1
y.neg[sub.neg$QCno] <- y.neg[sub.neg$QCno] + 1
}
max_QCno <- max(which(counter!=0))
pr.y.poz = y.poz[seq_len(max_QCno)]/counter[seq_len(max_QCno)]
pr.y.neg = y.neg[seq_len(max_QCno)]/counter[seq_len(max_QCno)]
plot.data <- data.frame(QCno = rep(seq_len(max_QCno),2),
pr.y = c(pr.y.poz, pr.y.neg),
group = ifelse(rep(type==1,2*max_QCno),
c(rep("Mean increase",max_QCno),
rep("Mean decrease",max_QCno)),
c(rep("Variability increase",max_QCno),
rep("Variability decrease",max_QCno))),
metric = rep(metric,max_QCno*2))
return(plot.data)
}
###########################################################################################
#DESCRIPTION : returns a data frame for all the metrics, probability of out of range peptides, wheter it is for metric mean or metric dispersion and i is an increase or decrease
XmR.Summary.DataFrame <- function(prodata, data.metrics, L, U,listMean, listSD, guidset_selected) {
dat <- data.frame(QCno = c(),
pr.y = c(),
group = c(),
metric = c())
for (metric in data.metrics) {
data.1 <- XmR.Summary.prepare(prodata, metric = metric, L, U,type = 1,selectMean = listMean[[metric]],selectSD = listSD[[metric]], guidset_selected)
data.2 <- XmR.Summary.prepare(prodata, metric = metric, L, U,type = 2,selectMean = listMean[[metric]],selectSD = listSD[[metric]], guidset_selected)
data.2$pr.y <- -(data.2$pr.y)
dat <- rbind(dat, data.1, data.2)
}
return(dat)
}
############################################################################################
heatmap.DataFrame <- function(prodata, data.metrics,method,peptideThresholdRed,peptideThresholdYellow,
L, U, type,listMean, listSD, guidset_selected) {
time <- c()
val <- c()
met <- c()
flag <- c()
for (metric in data.metrics) {
df <- Decision.DataFrame.prepare(prodata, metric, method, peptideThresholdRed,peptideThresholdYellow,
L, U,type,selectMean=listMean[[metric]],selectSD=listSD[[metric]], guidset_selected)
time_df <- as.character(df$AcquiredTime)
val_df <- df$pr.y
met_df <- rep(metric,length(val_df))
flag_df <- df$flag
time <- c(time,time_df)
val <- c(val,val_df)
met <- c(met,met_df)
flag <- c(flag,flag_df)
}
dataFrame <- data.frame(time = time,
value = val,
metric = met,
flag = flag,
QCno = df$QCno
)
return(dataFrame)
}
###############################################################################################################
XmR.Radar.Plot.prepare <- function(prodata,L,U, metric, type,group,XmR_type,selectMean,selectSD, guidset_selected) {
precursors <- levels(prodata$Precursor)
precursors2 <- substring(precursors, first = 1, last = 10)
QCno.length <- c()
QCno.out.range.poz <- c()
QCno.out.range.neg <- c()
for(j in seq_len(length(precursors))) {
metricData <- getMetricData(prodata, precursors[j], L = L, U = U, metric = metric, normalization = T,selectMean,selectSD, guidset_selected)
QCno.length <- c(QCno.length,length(metricData))
plot.data <- XmR.data.prepare(prodata, metricData , L = L, U = U, type ,selectMean,selectSD, guidset_selected)
QCno.out.range.poz <- c(QCno.out.range.poz,length(plot.data[plot.data$t >= plot.data$UCL, ]$QCno))
QCno.out.range.neg <- c(QCno.out.range.neg,length(plot.data[plot.data$t <= plot.data$LCL, ]$QCno))
}
#pr <- Compute.QCno.OutOfRangePeptide.XmR(prodata,L,U,metric = metric,type = type, XmR.type,selectMean,selectSD, guidset_selected)
if(XmR_type == "poz") {
dat <- data.frame(peptides = precursors2,
OutRangeQCno = QCno.out.range.poz,
group = rep(group,length(precursors)),
orderby = seq(seq_len(length(precursors))),
metric = rep(metric, length(precursors)),
tool = rep("XmR",length(precursors)),
probability = QCno.out.range.poz/QCno.length
)
} else {
dat <- data.frame(peptides = precursors2,
OutRangeQCno = QCno.out.range.neg,
group = rep(group,length(precursors)),
orderby = seq(seq_len(length(precursors))),
metric = rep(metric, length(precursors)),
tool = rep("XmR",length(precursors)),
probability = QCno.out.range.neg/QCno.length
)
}
return(dat)
#return(list(dat.poz,dat.neg))
}
################################################################################################
XmR.Radar.Plot.DataFrame <- function(prodata, data.metrics, L,U,listMean,listSD,guidset_selected) {
dat <- data.frame(peptides = c(), OutRangeQCno = c(), group = c(),
orderby = c(), metric = c(), tool = c(),
probability = c()
)
for (metric in data.metrics) {
data.1 <- XmR.Radar.Plot.prepare(prodata,L,U,metric = metric, type = 1,group = "Mean increase",XmR_type = "poz",
selectMean = listMean[[metric]],selectSD = listSD[[metric]],guidset_selected)
data.2 <- XmR.Radar.Plot.prepare(prodata,L,U,metric = metric, type = 1,group = "Mean decrease",XmR_type = "neg",
selectMean = listMean[[metric]],selectSD = listSD[[metric]],guidset_selected)
data.3 <- XmR.Radar.Plot.prepare(prodata,L,U,metric = metric, type = 2,group = "Variability increase",XmR_type = "poz",
selectMean = listMean[[metric]],selectSD = listSD[[metric]],guidset_selected)
data.4 <- XmR.Radar.Plot.prepare(prodata,L,U,metric = metric, type = 2,group = "Variability decrease",XmR_type = "neg",
selectMean = listMean[[metric]],selectSD = listSD[[metric]],guidset_selected)
dat <- rbind(dat, data.1, data.2, data.3, data.4)
}
return(dat)
}
#################################################################################################################
CUSUM.Radar.Plot.prepare <- function(prodata,L,U, metric,type,group, CUSUM.type,selectMean,selectSD, guidset_selected) {
#precursors <- levels(reorder(prodata$Precursor,prodata[,COL.BEST.RET]))
precursors <- levels(prodata$Precursor)
precursors2 <- substring(precursors, first = 1, last = 10)
QCno.length <- c()
QCno.out.range <- c()
for(j in seq_len(length(precursors))) {
metricData <- getMetricData(prodata, precursors[j], L = L, U = U, metric = metric, normalization = T,selectMean,selectSD, guidset_selected)
QCno.length <- c(QCno.length,length(metricData))
plot.data <- CUSUM.data.prepare(prodata, metricData, precursors[j], type)
if(CUSUM.type == "poz")
QCno.out.range <- c(QCno.out.range,length(plot.data[plot.data$CUSUM.poz >= 5 | plot.data$CUSUM.poz <= -5, ]$QCno))
else
QCno.out.range <- c(QCno.out.range,length(plot.data[plot.data$CUSUM.neg >= 5 | plot.data$CUSUM.neg <= -5, ]$QCno))
}
#pr <- Compute.QCno.OutOfRangePeptide.CUSUM(prodata,L,U,metric = metric,type = type, CUSUM.type,selectMean,selectSD, guidset_selected)
dat <- data.frame(peptides = precursors2,
OutRangeQCno = QCno.out.range,
group = rep(group,length(precursors)),
orderby = seq(seq_len(length(precursors))),
metric = rep(metric, length(precursors)),
tool = rep("XmR",length(precursors)),
probability = QCno.out.range/QCno.length
)
return(dat)
}
#################################################################################################
CUSUM.Radar.Plot.DataFrame <- function(prodata, data.metrics, L,U,listMean,listSD,guidset_selected) {
dat <- data.frame(peptides = c(), OutRangeQCno = c(), group = c(),
orderby = c(), metric = c(), tool = c(),
probability = c()
)
for (metric in data.metrics) {
data.1 <- CUSUM.Radar.Plot.prepare(prodata,L,U, metric = metric, type = 1, group = "Mean increase", CUSUM.type = "poz",selectMean = listMean[[metric]],selectSD = listSD[[metric]],guidset_selected)
data.2 <- CUSUM.Radar.Plot.prepare(prodata,L,U, metric = metric, type = 1, group = "Mean decrease", CUSUM.type = "neg",selectMean = listMean[[metric]],selectSD = listSD[[metric]],guidset_selected)
data.3 <- CUSUM.Radar.Plot.prepare(prodata,L,U, metric = metric, type = 2, group = "Variability increase", CUSUM.type = "poz",selectMean = listMean[[metric]],selectSD = listSD[[metric]],guidset_selected)
data.4 <- CUSUM.Radar.Plot.prepare(prodata,L,U, metric = metric, type = 2, group = "Variability decrease", CUSUM.type = "neg",selectMean = listMean[[metric]],selectSD = listSD[[metric]],guidset_selected)
dat <- rbind(dat, data.1, data.2, data.3, data.4)
}
return(dat)
}
#######################################################################################################
Decision.DataFrame.prepare <- function(prodata, metric, method, peptideThresholdRed,peptideThresholdYellow,
L, U,type,selectMean,selectSD, guidset_selected) {
AcquiredTime <- prodata$AcquiredTime
QCno <- seq_len(nrow(prodata))
y <- rep(0,nrow(prodata))
counter <- rep(0,nrow(prodata))
#precursors <- levels(reorder(prodata$Precursor,prodata[,COL.BEST.RET]))
precursors <- levels(prodata$Precursor)
if(method == "XmR") {
for(precursor in precursors) {
metricData <- getMetricData(prodata, precursor, L = L, U = U, metric = metric, normalization = T,selectMean,selectSD, guidset_selected)
counter[seq_len(length(metricData))] <- counter[seq_len(length(metricData))]+1
plot.data <- XmR.data.prepare(prodata, metricData , L , U , type,selectMean,selectSD, guidset_selected)
sub <- plot.data[plot.data$InRangeOutRange == "OutRange",]
y[sub$QCno] <- y[sub$QCno] + 1
}
} else if(method == "CUSUM") {
for(precursor in precursors) {
metricData <- getMetricData(prodata, precursor, L = L, U = U, metric = metric, normalization = T,selectMean,selectSD, guidset_selected)
counter[seq_len(length(metricData))] <- counter[seq_len(length(metricData))]+1
plot.data <- CUSUM.data.prepare(prodata, metricData, precursor, type)
sub <- plot.data[(plot.data$CUSUM.poz >= 5 | plot.data$CUSUM.poz <= -5) | (plot.data$CUSUM.neg >= 5 | plot.data$CUSUM.neg <= -5), ]
y[sub$QCno] <- y[sub$QCno] + 1
}
}
max_QCno <- max(which(counter!=0))
pr.y = y[seq_len(max_QCno)]/counter[seq_len(max_QCno)]
plot.data <- data.frame(AcquiredTime = AcquiredTime[seq_len(max_QCno)],
QCno = rep(seq_len(max_QCno),1),
pr.y = pr.y,
group = ifelse(rep(type==1,max_QCno),
rep("Metric mean",max_QCno),
rep("Metric dispersion",max_QCno)
),
metric = rep(metric,max_QCno),
flag = rep(0,max_QCno)
)
for (i in seq_len(max_QCno)) {
if(plot.data$pr.y[i] > peptideThresholdRed){
plot.data$flag[i] <- "Fail"
}
else if(plot.data$pr.y[i] > peptideThresholdYellow){
plot.data$flag[i] <- "Warning"
}
else {
plot.data$flag[i] <- "Pass"
}
}
if(type == 2) {
return(plot.data[-1,])
}
return(plot.data)
}
#######################################################################################################
number.Of.Out.Of.Range.Metrics <- function(prodata,data.metrics,method, peptideThresholdRed,
peptideThresholdYellow, L, U, type,listMean,listSD,
guidset_selected) {
metricCounterAboveRed = 0
metricCounterAboveYellowBelowRed = 0
precursors <- levels(prodata$Precursor)
for (metric in data.metrics) {
QCno <- seq_len(nrow(prodata))
y <- rep(0,nrow(prodata))
counter <- rep(0,nrow(prodata))
for(precursor in precursors) {
metricData <- getMetricData(prodata, precursor, L = L, U = U, metric = metric, normalization = T,listMean[[metric]],listSD[[metric]], guidset_selected)
counter[seq_len(length(metricData))] <- counter[seq_len(length(metricData))]+1
plot.data <- XmR.data.prepare(prodata, metricData , L , U , type,selectMean = listMean[[metric]],selectSD = listSD[[metric]], guidset_selected)
sub <- plot.data[plot.data$InRangeOutRange == "OutRange",]
y[sub$QCno] <- y[sub$QCno] + 1
}
max_QCno <- max(which(counter!=0))
pr.y = y[seq_len(max_QCno)]/counter[seq_len(max_QCno)]
if(type == 2) {
pr.y <- pr.y[-1]
}
aboveYellow <- which(pr.y > peptideThresholdYellow)
aboveYellowBelowRed <- which(pr.y > peptideThresholdYellow & pr.y <= peptideThresholdRed)
if(length(which(pr.y > peptideThresholdRed)) > 0) {
metricCounterAboveRed = metricCounterAboveRed + 1
}
if(length(aboveYellowBelowRed) > 0) {
metricCounterAboveYellowBelowRed = metricCounterAboveYellowBelowRed + 1
}
}
return(c(metricCounterAboveRed,metricCounterAboveYellowBelowRed))
}
####################################################################################################
decisionRule_warning_message_XmR <- function(prodata, data.metrics, method, peptideThresholdRed, peptideThresholdYellow,metricThresholdRed,metricThresholdYellow, L,U, type, listMean,listSD , guidset_selected) {
a1 <- number.Of.Out.Of.Range.Metrics(prodata,data.metrics,method = "XmR", peptideThresholdRed,peptideThresholdYellow,
L, U, type = 1,listMean,listSD, guidset_selected)
XmRCounterAboveRed1 <- a1[1]
XmRCounterAboveYellow1 <- a1[2]
a2 <- number.Of.Out.Of.Range.Metrics(prodata,data.metrics,method = "XmR", peptideThresholdRed,peptideThresholdYellow,
L, U, type = 2,listMean,listSD, guidset_selected)
XmRCounterAboveRed2 <- a2[1]
XmRCounterAboveYellow2 <- a2[2]
if(XmRCounterAboveRed1 > metricThresholdRed && XmRCounterAboveRed2 > metricThresholdRed) return({paste("XmR RED FLAG: System performance is UNACCEPTABLE")})
if(XmRCounterAboveRed1 > metricThresholdRed && XmRCounterAboveRed2 <= metricThresholdRed) return({paste("XmR RED FLAG: System performance is UNACCEPTABLE")})
if(XmRCounterAboveRed1 <= metricThresholdRed && XmRCounterAboveRed2 > metricThresholdRed) return({paste("XmR RED FLAG: System performance is UNACCEPTABLE")})
if(XmRCounterAboveYellow1 > metricThresholdYellow && XmRCounterAboveYellow2 > metricThresholdYellow) return({"XmR Yellow FLAG: System performance is POOR"})
if(XmRCounterAboveYellow1 <= metricThresholdYellow && XmRCounterAboveYellow2 > metricThresholdYellow) return({"XmR Yellow FLAG: System performance is POOR"})
if(XmRCounterAboveYellow1 > metricThresholdYellow && XmRCounterAboveYellow2 <= metricThresholdYellow) return({"XmR Yellow FLAG: System performance is POOR"})
return({"XmR BLUE FLAG: System performance is acceptable"})
}
############################################################################################################
decisionRule_warning_message_CUSUM <- function(prodata, data.metrics, method, peptideThresholdRed, peptideThresholdYellow,metricThresholdRed,metricThresholdYellow, L,U, type, listMean,listSD, guidset_selected ) {
a1 <- number.Of.Out.Of.Range.Metrics(prodata,data.metrics,method = "CUSUM", peptideThresholdRed,peptideThresholdYellow,
L, U, type = 1,listMean,listSD, guidset_selected)
CUSUMCounterAboveRed1 <- a1[1]
CUSUMCounterAboveYellow1 <- a1[2]
a2 <- number.Of.Out.Of.Range.Metrics(prodata,data.metrics,method = "CUSUM", peptideThresholdRed,peptideThresholdYellow,
L, U, type = 2,listMean,listSD, guidset_selected)
CUSUMCounterAboveRed2 <- a2[1]
CUSUMCounterAboveYellow2 <-a2[2]
if(CUSUMCounterAboveRed1 > metricThresholdRed && CUSUMCounterAboveRed2 > metricThresholdRed) return({paste("CUSUM RED FLAG: System performance is UNACCEPTABLE")})
if(CUSUMCounterAboveRed1 > metricThresholdRed && CUSUMCounterAboveRed2 <= metricThresholdRed) return({paste("CUSUM RED FLAG: System performance is UNACCEPTABLE")})
if(CUSUMCounterAboveRed1 <= metricThresholdRed && CUSUMCounterAboveRed2 > metricThresholdRed) return({paste("CUSUM RED FLAG: System performance is UNACCEPTABLE")})
if(CUSUMCounterAboveYellow1 > metricThresholdYellow && CUSUMCounterAboveYellow2 > metricThresholdYellow) return({"CUSUM Yellow FLAG: System performance is POOR"})
if(CUSUMCounterAboveYellow1 <= metricThresholdYellow && CUSUMCounterAboveYellow2 > metricThresholdYellow) return({"CUSUM Yellow FLAG: System performance is POOR"})
if(CUSUMCounterAboveYellow1 > metricThresholdYellow && CUSUMCounterAboveYellow2 <= metricThresholdYellow) return({"CUSUM Yellow FLAG: System performance is POOR"})
return({"CUSUM BLUE FLAG: System performance is acceptable"})
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.