Commit a5fe5271 authored by joerg.buescher's avatar joerg.buescher
Browse files

little bugfixes here and there

parent b0780c82
*.log
*.DS_Store
......@@ -348,6 +348,7 @@ if (prm$verbose >= 1) {
# write log of quality scores for training to tsv file
write.table(qslog, file = file.path(prm$batchdir, 'qslog.tsv'), sep = '\t', row.names = FALSE)
# write manual_peakcheck template
nowscore <- matrix(nrow = prm$nsmpl, ncol = prm$nmet, data = NA)
# use old scores if available and peaks match
......@@ -369,6 +370,16 @@ if (prm$verbose >= 1) {
}
}
}
} else { # if there is no previous solution, at least fill in the clear cases
for (id in 1:prm$nsmpl) {
for (im in 1:prm$nmet) {
if (msd[[im]][[id]]$qs < 0.2) {
nowscore[id,im] <- 0
} else if (msd[[im]][[id]]$qs > 0.95) {
nowscore[id,im] <- 2
}
}
}
}
# write template for training solution to xlsx file
......@@ -389,6 +400,7 @@ if (prm$verbose >= 1) {
openxlsx::saveWorkbook(wb, file.path(prm$batchdir, "manual_peakcheck_template.xlsx"), overwrite = TRUE)
}
msd # return msd
} # endfunction
......
......@@ -38,6 +38,7 @@ evaluate_peaks <- function(msd, prm){
final_data <- as.data.frame(c(msd[[im]][[id]]$qs_list, nowquantile))
final_data_prep <- predict(prm$train_preprocessing_final,final_data)
final_pred <- predict(prm$model_final,final_data_prep,type="prob")$H
msd[[im]][[id]]$report_p <- final_pred
#update plotdata label1
nowstartpos <- gregexpr('\n', msd[[im]][[id]]$label1)[[1]][2]
......
......@@ -153,7 +153,8 @@ get_quality_scores <- function(x1, x2, x3, peakstart, peaktop, peakend, nowusetr
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
ratio_value <- (1 / (1 + abs( log(x1[peaktop] / x2[peaktop], 2) -
log(metab[[im]]$quant$intens / metab[[im]]$qual12c$intens[1],2))^2 ) )^2
}
qs_list$QS_1vs2 <- check_qs_value(ratio_value)
......
......@@ -32,7 +32,9 @@ peak_detection <- function(metab, smpl, prm) {
'peakend' = NA, # end of peak (scan number)
'ylim' = c(0, 100), # y axis limits for plotting
'label1' = character(), # simple label for plot
'label2' = character() # debugging label for plot
'label2' = character(), # debugging label for plot
'report' = TRUE, # report in excel output only if quality is sufficient
'report_p' = 0 # probability that peak should be reported, will be filled based on finalp model in evaluate_peaks
)
msd2 <- list()
......@@ -145,10 +147,10 @@ peak_detection <- function(metab, smpl, prm) {
msd[[im]][[id]]$combidata <- nowdata
if(TRUE){
x1 <- msd[[im]][[id]]$x1
x2 <- msd[[im]][[id]]$x2
x3 <- msd[[im]][[id]]$x3
nowdata <- msd[[im]][[id]]$coimidata
# x1 <- msd[[im]][[id]]$x1
# x2 <- msd[[im]][[id]]$x2
# x3 <- msd[[im]][[id]]$x3
# nowdata <- msd[[im]][[id]]$coimidata
nowusetrace <- as.logical( c( 1, c(metab[[im]]$qual12c$use, 0)[1] , c(metab[[im]]$qual13c$use,0)[1] ) )
nowusetrace[is.na(nowusetrace)] <- FALSE
......@@ -447,7 +449,7 @@ peak_detection <- function(metab, smpl, prm) {
qsnames <- names( peak_list[[bestpeak]]$qs_list )
# append qslog info for training (last 2 cols are x12_13c_flags for qs median values)
for(p in 1:length(peak_list)){
for (p in 1:length(peak_list)) {
metab_name <- metab[[im]]$name
sample_name <- smpl$samplenames[id]
peak_id <- p
......
......@@ -6,14 +6,14 @@ plot_peaks <- function(metab, smpl, msd, prm) {
cat('Creating peakoverview.pdf ', file=prm$log_con)
print(' ')
pdf_file <- "peakoverview.pdf"
pdf_file <- file.path(prm$batchdir,"peakoverview.pdf")
pdf(file = pdf_file, height = 2*prm$nsmpl, width = 5* prm$nmet, family = "Helvetica")
par(mai = c(0.5, 0.5, 0.8, 0.5))
layout(matrix(c(1:(prm$nmet*prm$nsmpl)), prm$nsmpl, prm$nmet, byrow = TRUE)) # initiate subplots
# Reverse looping and indexing (samples first, then metabolites)
for (id in 1:prm$nsmpl){
for (id in prm$typeorder){
for (im in 1:prm$nmet){
if (prm$verbose >=2) {
cat('\r', paste('plot_peaks: metab', im, 'sample', id, ' '))
......
......@@ -36,7 +36,7 @@ process_batch <- function(nowfolder = "", parameters = list(), return_msd=FALSE)
prm <- initialize_prm()
# create new log file
prm$log_con <- file(file.path(nowfolder, "R_messages.log"),open="a")
prm$log_con <- file(file.path(nowfolder, 'R_messages.log'),open="a")
# Update prm from tsv config file and arguments
prm <- update_prm(nowfolder, parameters, prm)
......@@ -88,7 +88,7 @@ process_batch <- function(nowfolder = "", parameters = list(), return_msd=FALSE)
prm$model_final <- model_finalp
prm$train_preprocessing_final <- train_preprocessing_finalp
# only keep what we need
rm('model_finalp', 'train_preprocessing_finalp', 'train_dummy_vars_finalp')
rm('model_finalp', 'train_preprocessing_finalp')
if (prm$verbose >= 2) {
print(paste0('Successfully loaded 2. model from ', prm$model_final_path))
}
......
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