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

data dependent determination of smoothing factors in consensus_peak_update

parent a5fe5271
......@@ -20,6 +20,20 @@ for (id in 1:prm$nsmpl) {
}
refsmpl <- which.max(bestqs)[1]
# find typical fwhm to determine smoothing factor
if ('fwhm' %in% names(msd[[1]][[1]]) ) {
allfwhm <- numeric()
allqs <- numeric()
for (id in 1:prm$nsmpl) {
allfwhm <- c(allfwhm, unlist(sapply(msd, '[[', id)['fwhm', ]) )
allqs <- c(allqs, unlist(sapply(msd, '[[', id)['qs', ]) )
}
typicalfwhm <- median(allfwhm[allqs > 0.9], na.rm = TRUE)
} else {
typicalfwhm <- 15
}
# find shift relative to refsmpl
corweight <- list()
for (id in 1:prm$nsmpl) {
......@@ -191,9 +205,10 @@ for (im in 1:prm$nmet) {
p2[is.na(p2)] <- 0
p3[is.na(p3)] <- 0
# peak detection
# peak detection on weighted sum of single traces
p00 <- rowSums(cbind(p1/max(p1,na.rm = TRUE), p2/max(p2,na.rm = TRUE), p3/max(p3,na.rm = TRUE)), na.rm = TRUE)
# determine candidates and their qs scores
p0 <- zoo::rollmean(rowSums(cbind(p1/max(p1,na.rm = TRUE), p2/max(p2,na.rm = TRUE), p3/max(p3,na.rm = TRUE)), na.rm = TRUE), 3, fill=0, align='center')
p0 <- zoo::rollmean(p00, round(typicalfwhm/5), fill=0, align='center')
sp0 <- diff(p0)
maxsp0 <- max(abs(sp0))
local_max <- which(diff(sign(sp0))==-2)
......@@ -207,7 +222,7 @@ for (im in 1:prm$nmet) {
}
}
# second pass with more smoothing
p0 <- zoo::rollmean(rowSums(cbind(p1/max(p1,na.rm = TRUE), p2/max(p2,na.rm = TRUE), p3/max(p3,na.rm = TRUE)), na.rm = TRUE), 15, fill=0, align='center')
p0 <- zoo::rollmean(p00, typicalfwhm, fill=0, align='center')
sp0 <- diff(p0)
maxsp0 <- max(abs(sp0))
local_max <- which(diff(sign(sp0))==-2)
......@@ -233,7 +248,7 @@ for (im in 1:prm$nmet) {
nowpend <- prm$unirt[max(c(3, min(c(prm$nscan, peakend[ip] + tr[1])) )) ]
peakoldqs[ip] <- sum(originalqs[intersect(which((originalrt - metabshift) > nowpstart),
which((originalrt - metabshift) < nowpend))]) /
sqrt(peakend[ip] - peakstart[ip])
sqrt(peakend[ip] - peakstart[ip])
}
}
peakoldqs[is.na(peakoldqs)] <- 0
......
......@@ -34,7 +34,7 @@ peak_detection <- function(metab, smpl, prm) {
'label1' = character(), # simple 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
'fwhm' = NA # full width half maximum
)
msd2 <- list()
......@@ -217,7 +217,7 @@ peak_detection <- function(metab, smpl, prm) {
qs_list <- list()
peak_list <- list()
peakstarts <- peakends <- peak_min_start <- peak_max_end <- rep(NA,prm$initial_candidates)
qs <- area1 <- area2 <- area13c <- height1 <- height2 <- height13c <- RT <- rep(NA,prm$initial_candidates)
fwhm <- qs <- area1 <- area2 <- area13c <- height1 <- height2 <- height13c <- RT <- rep(NA,prm$initial_candidates)
# peakparameter adjustment
for (i in 1:length(peaktops)){
......@@ -294,6 +294,11 @@ peak_detection <- function(metab, smpl, prm) {
area2[i] <- if(!is.null(x2)) sum(x2[peakstarts[i]:peakends[i]])/prm$samplingfrequency else NA
area13c[i] <- if(!is.null(x3)) sum(x3[peakstarts[i]:peakends[i]])/prm$samplingfrequency else NA
height1[i] <- x1[peaktops[i]]
if (!is.na(height1[i])) {
fwhm[i] <- min(c(peakends[i], intersect(which(x1 < 0.5*height1[i]), peaktops[i]:prm$nscan)), na.rm = TRUE) -
max(c(peakstarts[i], intersect(which(x1 < 0.5*height1[i]), 1:peaktops[i])), na.rm = TRUE)
}
# save peak-heights for met
if (length(x2) == 0){
......@@ -354,6 +359,7 @@ peak_detection <- function(metab, smpl, prm) {
peak$peakstart <- peakstarts[i]
peak$peaktop <- peaktops[i]
peak$peakend <- peakends[i]
peak$fwhm <- fwhm[i]
peak$ylim <- c(0,1.1*max(c(100,height1[i],nowusetrace[2]*height2[i],nowusetrace[3]*height13c[i]),na.rm = TRUE))
peak$label1 <- paste0(metab[[im]]$name, '\n',
smpl$samplenames[[id]], '\n',
......@@ -442,6 +448,7 @@ peak_detection <- function(metab, smpl, prm) {
msd[[im]][[id]]$peakstart <- peak_list[[bestpeak]]$peakstart
msd[[im]][[id]]$peaktop <- peak_list[[bestpeak]]$peaktop
msd[[im]][[id]]$peakend <- peak_list[[bestpeak]]$peakend
msd[[im]][[id]]$fwhm <- peak_list[[bestpeak]]$fwhm
msd[[im]][[id]]$ylim <- peak_list[[bestpeak]]$ylim
msd[[im]][[id]]$label1 <- peak_list[[bestpeak]]$label1
msd[[im]][[id]]$label2 <- peak_list[[bestpeak]]$label2
......
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