Dear Gitlab users, due to maintenance reasons, Gitlab will not be available on Thursday 30.09.2021 from 5:00 pm to approximately 5:30 pm.

Commit 7190b4b6 authored by joerg.buescher's avatar joerg.buescher
Browse files

bugfix in case of NA peak candidate

parent 48a3d2c4
......@@ -28,148 +28,150 @@ get_quality_scores <- function(x1, x2, x3, peakstart, peaktop, peakend, nowusetr
'QS_1vs2' = 0 ,
'QS_dRTs' = 0)
if ((peakstart < peaktop) && (peakend > peaktop) ) {
if (!is.na(peaktop)) {
if ((peakstart < peaktop) && (peakend > peaktop) ) {
# QS1: average correlation (quant, qual12c, qual13c)
nowcor <- suppressWarnings(cor(cbind(x1[peakstart:peakend],x2[peakstart:peakend],x3[peakstart:peakend]),use="everything") ) #use="complete.obs"
cor_value <- 0.5*(mean(c(nowcor[1, 2:3], nowcor[2,3]), na.rm = TRUE)+1)
qs_list$QS_cor123 <- check_qs_value(cor_value)
# QS1: average correlation (quant, qual12c, qual13c)
nowcor <- suppressWarnings(cor(cbind(x1[peakstart:peakend],x2[peakstart:peakend],x3[peakstart:peakend]),use="everything") ) #use="complete.obs"
cor_value <- 0.5*(mean(c(nowcor[1, 2:3], nowcor[2,3]), na.rm = TRUE)+1)
qs_list$QS_cor123 <- check_qs_value(cor_value)
# QS2 ratio peak top to highest point outside of series
t_o_value <- max(c(x1[peaktop] / (x1[peaktop] + max(x1[c(1:(max(c(peakstart-1,1))),(min(c(peakend+1,length(x1)))):length(x1))], na.rm = TRUE)),0))
qs_list$QS_T.o <- check_qs_value(t_o_value)
# QS2 ratio peak top to highest point outside of series
t_o_value <- max(c(x1[peaktop] / (x1[peaktop] + max(x1[c(1:(max(c(peakstart-1,1))),(min(c(peakend+1,length(x1)))):length(x1))], na.rm = TRUE)),0))
qs_list$QS_T.o <- check_qs_value(t_o_value)
# QS3 ratio peak top/max(peakend,peakstart)
t_b_h_value <- ((2 * x1[peaktop])/ (x1[peaktop] + max(x1[peakstart],x1[peakend])) ) -1
qs_list$QS_T.h <- check_qs_value(t_b_h_value)
# QS3 ratio peak top/max(peakend,peakstart)
t_b_h_value <- ((2 * x1[peaktop])/ (x1[peaktop] + max(x1[peakstart],x1[peakend])) ) -1
qs_list$QS_T.h <- check_qs_value(t_b_h_value)
# QS4 ratio peak top/min(peakend,peakstart)
t_b_l_value <- ((2 * x1[peaktop])/ (x1[peaktop] + min(x1[peakstart],x1[peakend])) ) -1
qs_list$QS_T.l <- check_qs_value(t_b_l_value)
# QS4 ratio peak top/min(peakend,peakstart)
t_b_l_value <- ((2 * x1[peaktop])/ (x1[peaktop] + min(x1[peakstart],x1[peakend])) ) -1
qs_list$QS_T.l <- check_qs_value(t_b_l_value)
# QS5 deviation from expected RT
# algo in mrm_integrate: max(c(1/(1 + (0.02*abs(peaktops[i] - metparam$peaktopa) / (1+metparam$peaktopadelta) ) )^2,0))
peaktopa <- round(metab[[im]]$RT * 60 * prm$samplingfrequency)
dRT_value <- (1/(1+abs(prm$unirt[peaktop + rangestart] - prm$unirt[peaktopa]) ))^2
qs_list$QS_dRT <- check_qs_value(dRT_value)
# QS5 deviation from expected RT
# algo in mrm_integrate: max(c(1/(1 + (0.02*abs(peaktops[i] - metparam$peaktopa) / (1+metparam$peaktopadelta) ) )^2,0))
peaktopa <- round(metab[[im]]$RT * 60 * prm$samplingfrequency)
dRT_value <- (1/(1+abs(prm$unirt[peaktop + rangestart] - prm$unirt[peaktopa]) ))^2
qs_list$QS_dRT <- check_qs_value(dRT_value)
# deviation from predicted RT
if (prm$ml_type != 'mltrain_initial') {
peaktopa <- round(metab[[im]]$RT * 60 * prm$samplingfrequency) + rtshift
peaktopa <- min(max(c(1, peaktopa + rangestart)) , prm$nscan)
dRT_value <- (1/(1+abs(prm$unirt[peaktop] - prm$unirt[peaktopa]) ))^2
qs_list$QS_dRTs <- check_qs_value(dRT_value)
}
# deviation from predicted RT
if (prm$ml_type != 'mltrain_initial') {
peaktopa <- round(metab[[im]]$RT * 60 * prm$samplingfrequency) + rtshift
peaktopa <- min(max(c(1, peaktopa + rangestart)) , prm$nscan)
dRT_value <- (1/(1+abs(prm$unirt[peaktop] - prm$unirt[peaktopa]) ))^2
qs_list$QS_dRTs <- check_qs_value(dRT_value)
}
# QS6 widthscore
if (peakstart != peakend){
widthscore <- 1/ (1+(prm$unirt[peakend] - prm$unirt[peakstart]))
}else{
widthscore <- 0
}
qs_list$QS_w <- check_qs_value(widthscore)
# QS6 widthscore
if (peakstart != peakend){
widthscore <- 1/ (1+(prm$unirt[peakend] - prm$unirt[peakstart]))
}else{
widthscore <- 0
}
qs_list$QS_w <- check_qs_value(widthscore)
# QS7 | QS8 | QS9 average correlation (quant, qual12c, qual13c)
#Check if nowusetrace is 0
if (nowusetrace[1] == FALSE | nowusetrace[2] == FALSE) {
if(prm$ml_type == "mltrain_initial"){
cor_value <- 0
# QS7 | QS8 | QS9 average correlation (quant, qual12c, qual13c)
#Check if nowusetrace is 0
if (nowusetrace[1] == FALSE | nowusetrace[2] == FALSE) {
if(prm$ml_type == "mltrain_initial"){
cor_value <- 0
}else{
cor_value <- prm$median_values$QS_cor12
}
}else{
cor_value <- prm$median_values$QS_cor12
cor_value <- cor(x1[peakstart:peakend],x2[peakstart:peakend],use="everything") * 0.5 + 0.5
}
}else{
cor_value <- cor(x1[peakstart:peakend],x2[peakstart:peakend],use="everything") * 0.5 + 0.5
}
qs_list$QS_cor12 <- check_qs_value(cor_value)
qs_list$QS_cor12 <- check_qs_value(cor_value)
#Check if nowusetrace is 0
if (nowusetrace[1] == FALSE || nowusetrace[3] == FALSE) {
if(prm$ml_type == "mltrain_initial"){
cor_value <- 0
#Check if nowusetrace is 0
if (nowusetrace[1] == FALSE || nowusetrace[3] == FALSE) {
if(prm$ml_type == "mltrain_initial"){
cor_value <- 0
}else{
cor_value <- prm$
median_values$QS_cor13
}
}else{
cor_value <- prm$
median_values$QS_cor13
cor_value <- cor(x1[peakstart:peakend],x3[peakstart:peakend],use="everything") * 0.5 + 0.5
}
}else{
cor_value <- cor(x1[peakstart:peakend],x3[peakstart:peakend],use="everything") * 0.5 + 0.5
}
qs_list$QS_cor13 <- check_qs_value(cor_value)
#Check if nowusetrace is 0
if (nowusetrace[2] == FALSE || nowusetrace[3] == FALSE) {
if(prm$ml_type == "mltrain_initial"){
cor_value <- 0
qs_list$QS_cor13 <- check_qs_value(cor_value)
#Check if nowusetrace is 0
if (nowusetrace[2] == FALSE || nowusetrace[3] == FALSE) {
if(prm$ml_type == "mltrain_initial"){
cor_value <- 0
}else{
cor_value <- prm$median_values$QS_cor23
}
}else{
cor_value <- prm$median_values$QS_cor23
cor_value <- cor(x2[peakstart:peakend],x3[peakstart:peakend],use="everything") * 0.5 + 0.5
}
}else{
cor_value <- cor(x2[peakstart:peakend],x3[peakstart:peakend],use="everything") * 0.5 + 0.5
}
qs_list$QS_cor23 <- check_qs_value(cor_value)
qs_list$QS_cor23 <- check_qs_value(cor_value)
# QS10 Peak intensity
peak_int_value <- log10(x1[peaktop])
qs_list$QS_height <- check_qs_value(peak_int_value)
# QS10 Peak intensity
peak_int_value <- log10(x1[peaktop])
qs_list$QS_height <- check_qs_value(peak_int_value)
# QS11 #retention time relative to method length
RT_value <- prm$unirt[peaktop] / max(prm$unirt, na.rm = TRUE)
qs_list$QS_relRT <- check_qs_value(RT_value)
# QS11 #retention time relative to method length
RT_value <- prm$unirt[peaktop] / max(prm$unirt, na.rm = TRUE)
qs_list$QS_relRT <- check_qs_value(RT_value)
# QS12 | QS13 (12c & 13c qualifier): ratio peak top to highest point outside
if (is.null(x2) || nowusetrace[2] == FALSE || (peakstart == peakend)){
if(prm$ml_type == "mltrain_initial"){
qual_t_o_value <- 0
# QS12 | QS13 (12c & 13c qualifier): ratio peak top to highest point outside
if (is.null(x2) || nowusetrace[2] == FALSE || (peakstart == peakend)){
if(prm$ml_type == "mltrain_initial"){
qual_t_o_value <- 0
}else{
qual_t_o_value <- prm$median_values$QS_T.o2
}
}else{
qual_t_o_value <- prm$median_values$QS_T.o2
qual_t_o_value <- max((max(x2[peakstart:peakend],na.rm=TRUE) / ( (max(x2[peakstart:peakend],na.rm=TRUE) + max(c(x2[1:max(c(peakstart-1),1)],x2[min(c(peakend+1,length(x2))):length(x2)]),na.rm=TRUE)) )),0)
}
}else{
qual_t_o_value <- max((max(x2[peakstart:peakend],na.rm=TRUE) / ( (max(x2[peakstart:peakend],na.rm=TRUE) + max(c(x2[1:max(c(peakstart-1),1)],x2[min(c(peakend+1,length(x2))):length(x2)]),na.rm=TRUE)) )),0)
}
qs_list$QS_T.o2 <- check_qs_value(qual_t_o_value)
# gaussean shape of quantifier
nowoutsidepoints <- which(x1 < 0.5*x1[peaktop]) - peaktop
nowfwhm <- min(nowoutsidepoints[nowoutsidepoints > 0]) - max(nowoutsidepoints[nowoutsidepoints < 0])
qs_list$QS_gauss <- check_qs_value(max(c(0,cor(x1[peakstart:peakend], dnorm(peakstart:peakend,mean=(peaktop+1),sd=(0.5*nowfwhm)) )),na.rm=TRUE) )# correlation to gaussian peak
# ratio quant/12C_qual
# Check if nowusetrace is 0
if (nowusetrace[1] == FALSE || nowusetrace[2] == FALSE ||
is.na(metab[[im]]$quant$intens) || is.na(metab[[im]]$qual12c$intens) ) {
if(prm$ml_type == "mltrain_initial"){
ratio_value <- 0
}else{
ratio_value <- prm$median_values$QS_1vs2
qs_list$QS_T.o2 <- check_qs_value(qual_t_o_value)
# gaussean shape of quantifier
nowoutsidepoints <- which(x1 < 0.5*x1[peaktop]) - peaktop
nowfwhm <- min(nowoutsidepoints[nowoutsidepoints > 0]) - max(nowoutsidepoints[nowoutsidepoints < 0])
qs_list$QS_gauss <- check_qs_value(max(c(0,cor(x1[peakstart:peakend], dnorm(peakstart:peakend,mean=(peaktop+1),sd=(0.5*nowfwhm)) )),na.rm=TRUE) )# correlation to gaussian peak
# ratio quant/12C_qual
# Check if nowusetrace is 0
if (nowusetrace[1] == FALSE || nowusetrace[2] == FALSE ||
is.na(metab[[im]]$quant$intens) || is.na(metab[[im]]$qual12c$intens) ) {
if(prm$ml_type == "mltrain_initial"){
ratio_value <- 0
}else{
ratio_value <- prm$median_values$QS_1vs2
}
} else {
ratio_value <- (1 / (1 + abs(log(x1[peaktop] / x2[peaktop], 2) - log(metab[[im]]$quant$intens / metab[[im]]$qual12c$intens,2))^2 ) )^2
}
} else {
ratio_value <- (1 / (1 + abs(log(x1[peaktop] / x2[peaktop], 2) - log(metab[[im]]$quant$intens / metab[[im]]$qual12c$intens,2))^2 ) )^2
}
qs_list$QS_1vs2 <- check_qs_value(ratio_value)
qs_list$QS_1vs2 <- check_qs_value(ratio_value)
if (is.null(x3) || nowusetrace[3] == FALSE || (peakstart == peakend)){
if(prm$ml_type == "mltrain_initial"){
qual_t_o_value <- 0
if (is.null(x3) || nowusetrace[3] == FALSE || (peakstart == peakend)){
if(prm$ml_type == "mltrain_initial"){
qual_t_o_value <- 0
}else{
qual_t_o_value <- prm$median_values$QS_T.o3
}
}else{
qual_t_o_value <- prm$median_values$QS_T.o3
qual_t_o_value <- max((max(x3[peakstart:peakend],na.rm=TRUE) / ( (max(x3[peakstart:peakend],na.rm=TRUE) + max(c(x3[1:max(c(peakstart-1),1)],x3[min(c(peakend+1,length(x3))):length(x3)]),na.rm=TRUE)) )),0)
}
}else{
qual_t_o_value <- max((max(x3[peakstart:peakend],na.rm=TRUE) / ( (max(x3[peakstart:peakend],na.rm=TRUE) + max(c(x3[1:max(c(peakstart-1),1)],x3[min(c(peakend+1,length(x3))):length(x3)]),na.rm=TRUE)) )),0)
}
qs_list$QS_T.o3 <- check_qs_value(qual_t_o_value)
qs_list$QS_T.o3 <- check_qs_value(qual_t_o_value)
} # endif peaktop is between peakstart and peakend
} # endif peaktop is between peakstart and peakend
} # endif !is.na(peaktop)
# return qs_list
qs_list
......
......@@ -5,13 +5,13 @@ read_metabdb <- function(prm){
mrms <- openxlsx::read.xlsx(prm$excelfile)
cat(paste0('\nExcelfile: ',prm$excelfile,'\n'), file=prm$log_con)
cat(paste0('\nSize of excel table: ', as.character(nrow(mrms)),' by ', as.character(ncol(mrms))),file=prm$log_con)
# check for beak line
breakpos <- which(mrms[,1] == 'break')
if (length(breakpos) > 0) {
mrms <- mrms[1:(breakpos[1]-1), ]
}
# remove old '_13C' appendices to metabolite names
for (im in 1:nrow(mrms)) {
nownchar <- nchar(mrms$Name[im])
......@@ -21,18 +21,18 @@ read_metabdb <- function(prm){
}
}
}
# generate list of metabolites
mets <- unique(mrms[ ,'Name'])
# mets <- mets[which(regexpr(pattern='_13C',mets)<0)] # remove 13C from list
# pre-allocate output
metab <- list(length(mets))
chrominfo <- list("MRM" = character(),
"intens" = numeric(),
"use" = numeric())
minfo <- list("name" = character(),
"formula" = character(),
"HMDB" = character(),
......@@ -49,31 +49,29 @@ read_metabdb <- function(prm){
rawareasig = NA,
rawheightsig = NA
)
metcount <- 0
cat('\nNow reading metabolites from excel file. \n Counting metabolites... ', file=prm$log_con)
rows12c <- which(mrms$Isotope == '12C')
rows13c <- which(mrms$Isotope == '13C')
for (im in 1:length(mets)){
# log progress
cat(paste0(as.character(im),'-'), file=prm$log_con)
cat(paste0(as.character(im),'-'), file=prm$log_con)
# find rows in excel list
nowrows <- intersect( which(mrms[ ,'Name'] == mets[im]) , rows12c)
now13crow <- intersect( which(mrms[ ,'Name'] == mets[im]) , rows13c)
if (length(nowrows) >= 1) {
# Get indices of quantifier and qualifiers and metabolites to ignore
quantrow <- nowrows[which(mrms[nowrows,'Quantifier'] == 1)][1]
if (is.na(quantrow)) quantrow <- nowrows[1]
qualrow <- setdiff(nowrows, quantrow)
if (mrms[quantrow,'UseTrace'] != 0) {
metcount <- metcount + 1
# pre-allocate all entries
mi <- minfo
mi$name <- mrms[nowrows[1],'Name']
......@@ -82,30 +80,35 @@ read_metabdb <- function(prm){
mi$CAS <- mrms[nowrows[1],'CAS']
mi$KEGG <- mrms[nowrows[1],'KEGG']
mi$pubchem <- mrms[nowrows[1],'pubchem']
# Read retention time
if (!is.null(prm$RTcolumn)) {
if (prm$RTcolumn %in% colnames(mrms)) {
mi$RT <- mrms[nowrows[1],prm$RTcolumn]
mi$RT <- as.numeric(mrms[nowrows[1],prm$RTcolumn])
}
} else if ('RT' %in% colnames(mrms)) {
mi$RT <- mrms[nowrows[1],'RT']
mi$RT <- as.numeric(mrms[nowrows[1],'RT'])
} else {
print(paste0('ERROR in ', prm$excelfile, ': cannot find RT column'))
stop()
}
if (is.na(mi$RT)) {
# skip this entry
next
}
# preprocess data for subsequent steps
mi$RTnum <- round(mi$RT * 60 * prm$samplingfrequency)
mi$RTnum <- round(mi$RT * 60 * prm$samplingfrequency)
# Get quantifier info
mi$quant$MRM <- paste0(as.character(prm$polarities[mrms[quantrow,'Polarity']]),
mi$quant$MRM <- paste0(as.character(prm$polarities[mrms[quantrow,'Polarity']]),
as.character(round(as.numeric(mrms[quantrow,'Q1mz']),1)),
'>',
as.character(round(as.numeric(mrms[quantrow,'Q3mz']),1)) )
mi$quant$intens <- as.numeric(mrms[quantrow,'Abundance'])
mi$quant$use <- 1
# Get 12C qualifier info (can be multiple qualifiers in list)
if (length(qualrow) >= 1) {
for (iq in 1:length(qualrow)) {
......@@ -116,7 +119,7 @@ read_metabdb <- function(prm){
mi$qual12c$use[iq] <- as.numeric( mrms[qualrow[iq],'UseTrace'])
}
}
# Get 13C qualifier info (can be multiple qualifiers in list)
if (length(now13crow) >= 1) {
for (iq in 1:length(now13crow)) {
......@@ -127,13 +130,14 @@ read_metabdb <- function(prm){
mi$qual13c$use[iq] <- as.numeric( mrms[now13crow[iq],'UseTrace'])
}
}
metcount <- metcount + 1
metab[[metcount]] <- mi
} # endif 12C row has been found
} # endif use this metabolite
} # endfor read all MRM database entries
metab
} # endfunction
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment