diff --git a/FUB1outlier.numbers b/FUB1outlier.numbers new file mode 100644 index 0000000..2ec7e92 Binary files /dev/null and b/FUB1outlier.numbers differ diff --git a/FUB1outlier.xlsx b/FUB1outlier.xlsx new file mode 100644 index 0000000..7c5e8ae Binary files /dev/null and b/FUB1outlier.xlsx differ diff --git a/dev/ROUT.R b/dev/ROUT.R new file mode 100644 index 0000000..a23babb --- /dev/null +++ b/dev/ROUT.R @@ -0,0 +1,267 @@ +################################################################################ +# F.Innerbichler +# Jun 2026 +# Robust Outlier detection +################################################################################ + +# Code acc. to Mutulsky and Brown hpptps://pubmed.ncbi.nlm.nih.gov/16526949/ +# except line 175 p <- 2*(1-pt(z, length(residuals)-5)) was exchanged with p <- 2*(1-pcauchy(z)) to be robust + +library(minpack.lm) +library(ggplot2) +library(rstudioapi) + +# Getting path for current open file +current_path = rstudioapi::getActiveDocumentContext()$path +setwd(dirname(current_path)) + +outsPlot_FUN <- function(all_l,OUTs_, TS, PROC = "not specified", PROBE="", Q=0.01, par) { +#browser() + #all_l_rout <- all_l[KEEP,] + klein <- min(log(all_l$conc)) + gross <- max(log(all_l$conc)) + x_seq <- seq(klein, gross, (gross-klein)/100) + samTRUE <- f4pl(x=x_seq, bottom=par["a"], top=par["d"],hill=par["b"], logEC50=par["cs"]-par["r"]) + refTRUE <- f4pl(x=x_seq, bottom=par["a"], top=par["d"],hill=par["b"], logEC50=par["cs"]) + + pl_T <- data.frame(cbind(x_seq, refTRUE, samTRUE)) + p <- ggplot(all_l) + + geom_point(aes(x=log(conc), y=y, shape=factor(isRef)), size=2) + + theme_bw() + + p2 <- p + geom_point(OUTs_, mapping=aes(x=log(conc), y=y, shape=factor(isRef)), + size=5, color="violetred", stroke=2, inherit.aes = FALSE) + + scale_shape_manual(label=c("R",TS), values = c(21,25)) + + ggtitle(paste(PROC, PROBE, TS, "Threshold=", Q)) + + theme(legend.title=element_blank()) + p3 <- p2 + geom_line(data=pl_T, aes(x=x_seq, y=refTRUE), color="blue", inherit.aes = F) + + geom_line(data=pl_T, aes(x=x_seq, y=samTRUE), color="red", inherit.aes = F) + p3 + +} + +#### 4PL model ---- +f4pl <- function(x, bottom, top, logEC50, hill) { + bottom+(top-bottom)/(1+exp((logEC50-x)*hill)) +} + +#### residuals for joint REF/TEST fit ---- +resid_4pl_joint <- function(par,x,y, is_ref) { + bottom <- par["bottom"] + top <- par["top"] + hill <- par["hill"] + le50_r <- par["logEC50_ref"] + le50_t <- par["logEC50_test"] + + le50 <- ifelse(is_ref, le50_r, le50_t) + yhat <- f4pl(x, bottom, top, le50, hill) + y - yhat +} + + +jac_4pl_joint <- function(par, x,y,is_ref) { + bottom <- par["bottom"] + top <- par["top"] + hill <- par["hill"] + le50_r <- par["logEC50_ref"] + le50_t <- par["logEC50_test"] + + le50 <- ifelse(is_ref, le50_r, le50_t) + + # useful intermediates + + texp <- exp((le50-x)*hill) # t = 10^((le50-x)*hill)) + den <- (1-texp) + frac <- (top-bottom)/den + + dyhat_dbottom <- 1-1/den + dyhat_dtop <- 1/den + + d_invden_dle50 <- -(texp*log(2.718282)*hill)/(den^2) + + d_invden_dhill <- -(texp*log(2.718282)*(le50-x))/(den^2) + + dyhat_dle50 <- (top-bottom)* d_invden_dle50 + dyhat_dhill <- (top-bottom)* d_invden_dhill + + J <- cbind( + bottom = -dyhat_dbottom, + topm = -dyhat_dtop, + hill = -dyhat_dhill, + logEC50_ref = ifelse(is_ref, -dyhat_dle50,0), + logEC50_test = ifelse(!is_ref, -dyhat_dle50,0) + ) + J +} + + +robust_fit_irls <- function(par_start, x,y, is_ref, + max_outer = 30, control = nls.lm.control(maxiter=200, ftol=1e-10, ptol=1e-10)) { + par <- par_start + #browser() + for (k in seq_len(max_outer)) { + r <- resid_4pl_joint(par,x,y,is_ref) + + # Mot Brown + rSD <- quantile(r, 0.6827)*length(is_ref)/(length(is_ref)-5) + # Robust scale (MAD); fallback if MAD ~0 + # s <- median(abs(r))/0.6745 + + # Lorentzian / Cauchy weights acc. Mot/Brown + w <- 1/log(1+(abs(r)/rSD)^2) + # weighted residual function lor LM: sqrt(w)*r + fn_w <- function(p) { + rr <- resid_4pl_joint(p,x,y,is_ref) + sqrt(w)*rr + # w*rr + } + jac_w <- function(p) { + JJ <- jac_4pl_joint(p,x,y,is_ref) + JJ*sqrt(w) # row-wise scaling + # JJ*w + } + + #browser() + fit <- tryCatch(nls.lm(par=par, fn=fn_w, jac=jac_w, control=control), + error=function(e) { + paste("error at robust fit irls", k) + # tau <<- tau *10 + #Return zero delta to avoid NaN updates + return(rep(0, length(par))) + }) + + if (max(abs(fit$par -par)) < 1e-8) { + par <- fit$par + break + } + par <- fit$par + } + r_final <- resid_4pl_joint(par,x,y,is_ref) + # s_final <- median(abs(r_final))/0.6745 + rSD <- quantile(r, 0.6827)*length(r)/(length(r)-5) + # Robust scale (MAD); fallback if MAD ~0 + # s <- median(abs(r))/0.6745 + list(par=par, residuals=r_final, scale=rSD) + +} + + # ROUT outlier detection via BH-FDR [1] https://cran.r-project.org/web//packages/minpack.lm/refman/minpack.lm.html)[6] + +rout_detect <- function(residuals, scale, Q=0.01) { + z <- abs(residuals/scale) + p <- 2*(1 - pcauchy(z)) + # p <- 2*(1-pt(z, length(residuals)-5)) + + o <- order(p) + p_sorted <- p[o] + m <- length(p) + # thresh <- (seq_len(m)/m) # *Q + keep <- p_sorted <= Q # thresh + (out_idx <- o[keep]) + sort(out_idx) + +} + +# Final OLS LM refit after removing outliers + +ols_refit <- function(par_start, x,y,is_ref, keep_idx, + control = nls.lm.control(maxiter=400, ftol=1e-12, ptol=1e-12)) { + xk <- x[keep_idx]; yk <- y[keep_idx]; rk <- is_ref[keep_idx] + + fn_w <- function(p) resid_4pl_joint(p,xk,yk,rk) + jac_w <- function(p) jac_4pl_joint(p,xk,yk,rk) +#browser() + fit <- tryCatch(nls.lm(par=par_start, fn=fn_w, jac=jac_w, control=control), + error=function(e) { + paste("error at ols_refits") + # tau <<- tau *10 + #Return zero delta to avoid NaN updates + return(rep(0, length(par))) + }) + + list(par = fit$par, residuals = fn_w(fit$par), keep=keep_idx) +} + +#. high-level ROUT 4PL analysis + +rout_4pl_potency <- function(df, Q=0.01, par_start=NULL, max_outer = 30) { + stopifnot(all(c("sample","conc","y") %in% names(df))) # "rep", + #browser() + x <- log(df$conc) + y <- df$y + is_ref <- df$isRef == 1 + + # start values if not provided + if (is.null(par_start)) { + bottom0 = min(y, na.rm = T) + top0 = max(y, na.rm = T) + if (df$y[1]>df$y[8]) hill0 <- 1 else hill0 <- -1 + le50_r0 <- median(x[is_ref], na.rm = T) + le50_t0 <- median(x[!is_ref], na.rm = T) + par_start <- c(bottom=bottom0, top=top0, hill=hill0, logEC50_ref = le50_r0 , logEC50_test = le50_t0) + } + # Step1 : robust fit + robust <- robust_fit_irls(par_start, x,y,is_ref, max_outer = max_outer) + # Step 2: outlier detection + out_idx <- rout_detect(robust$residuals, robust$scale, Q=Q) + keep_idx <- setdiff(seq_len(nrow(df)), out_idx) + + # Step 3: OLS refit + ols <- ols_refit(robust$par, x,y,is_ref, keep_idx) + + # Relative potency + rp <- exp(ols$par["logEC50_ref"] - ols$par["logEC50_test"])*100 + + list(Q=Q, outliers=out_idx, + kept=keep_idx, + par_robust = robust$par, + par_final = ols$par, + rp_percent = rp, + df=df) + +} + +# sheets <- openxlsx::getSheetNames("~/plateflow/outlierFUB.xlsx") +# all_dat <- lapply(sheets, openxlsx::read.xlsx, xlsxFile="~/plateflow/outlierFUB.xlsx") +# names(all_dat) <- sheets +# +# Plate <- 1 +# +# PlanteN <- sheets[Plate] +# NoDils <- 8 +# Nreps <- 3 +# DAT <- all_dat[[Plate]] +# colnames(DAT) <- DAT[2,] +# REF <- DAT[3:10,1:4] +# SAM1 <- DAT[13:20,1:4] +# SAM2 <- DAT[23:30,1:4] +# SAM3 <- DAT[33:40,1:4] +# +# REF_ <- sapply(REF, as.numeric) +# colnames(REF_) <- colnames(REF) +# SAM1_ <- sapply(SAM1, as.numeric) +# colnames(SAM1_) <- colnames(SAM1) +# SAM2_ <- sapply(SAM2, as.numeric) +# colnames(SAM2_) <- colnames(SAM2) +# SAM3_ <- sapply(SAM3, as.numeric) +# colnames(SAM3_) <- colnames(SAM3) +# +# plot_df <- rbind(REF_, SAM1_, SAM2_,SAM3_) +# +# plot_df_ <- cbind(as.data.frame(plot_df), sample=c(rep("R", NoDils),rep("S1R", NoDils),rep("S2", NoDils),rep("S3", NoDils))) +# test_df <- plot_df_[1:16,] +# test_df2 <- plot_df_[c(1:8, 17:24),] +# test_df3 <- plot_df_[c(1:8, 25:32),] +# +# all_l <- melt(test_df, id.vars = c("sample","Dose"), variable.name="replname", value.name="readout") +# colnames(all_l) <- c("sample","conc","rep","y") +# +# for (Q in c(0.01,0.015,0.02)) { +# res <- rout_4pl_potency(all_l, Q) +# OUTs_ <- all_l[res$outliers,] +# print(outsPlot_FUN(all_l, OUTs_, TS=all_l$sample[9], PROC="ROUT",PROBE="Sheet 1",Q,par=res$par_final)) +# } + + + + diff --git a/dev/app.R b/dev/app.R index fbebd28..6512adb 100644 --- a/dev/app.R +++ b/dev/app.R @@ -31,7 +31,7 @@ library(scales) library(tolerance) source("../R/Global.R") - +source("ROUT.R") #### ui ---- @@ -177,7 +177,6 @@ server <- function(input, output, session) { ), uiOutput(outputId = "sheetName"), "For data format in the EXCEL file see Data template", - "If no data are uploaded, the settings to the right are used for calculations.", tags$head(tags$style(HTML("label {font-size:80%;margin-bottom: 3px;margin-top: 3px;}"))), div(checkboxInput("PureErr", "Should pure error be used for calculation of CIs?", FALSE), style = "font-size: 24px !important;color: #C2173F" @@ -300,7 +299,7 @@ server <- function(input, output, session) { ) ), tabPanel( - "Tests and ANOVAA", + "Tests and ANOVA", column( 12, h3("Tests for linear PLA:"), @@ -325,6 +324,16 @@ server <- function(input, output, session) { ) ) ), + tabPanel("robust outlier testing", + sliderInput("Qslider", "adjust Q-value in %",min=0.1, max = 20, value = 1, step=0.1), + plotOutput("OutlierPlot"), + tableOutput("OutlierDF"), + "GUIDANCE: The procedure of Motulsky & Brown allows for robust outlier testing.", + "Adjust the slider to mark the suspected outliers. Mostly, a Q-value of 2% is sufficient.", + "If the general variability of the data is high, many datapoints will be flagged, also ones that are not deemed to be outliers", + "An indicator for high variability is, when the Q-value needs to be increased above 5%, to flag suspected outliers.", + "Then, please re-consider, if the suspected outlier is not 'just' normal variability." + ), tabPanel( "parameter estimates", htmlOutput("PureErrWParEst"), @@ -729,7 +738,7 @@ server <- function(input, output, session) { reset(id = "") # from shinyjs package }) -#### input optim XL file ---- +#### input Wizard XL file ---- observe({ if (!is.null(input$MiFile)) { MinFile <- input$MiFile @@ -804,7 +813,50 @@ server <- function(input, output, session) { # all_l$readout[all_l$readout < 0] <- 0.01 REP$all_l <- all_l - #### XLSX eval ---- + ##### ROUT outlier testing ---- + #browser() + if(!is.null(input$Qslider)) { + all_lROUT <- all_l + colnames(all_lROUT) <- c("log_dose","sample","y","isRef","isSample","conc") + res <- rout_4pl_potency(all_lROUT, Q=input$Qslider/100) + OUTs_ <- all_lROUT[res$outliers,] + + all_l_rout <- all_lROUT[res$kept,] + # all_l_rout$log_dose <- log(all_l_rout$Conc) + # colnames(all_l_rout) <- c("log_dose","sample","y","isRef","isSample","conc") + if (all_l_rout$conc[1]>all_l_rout$conc[6]) { + if (all_l_rout$y[1]>all_l_rout$y[6]) SLOPE <- 1 else SLOPE<- -1 + } else { + if (all_l_rout$y[1]>all_l_rout$y[6]) SLOPE <- -1 else SLOPE<- 1 + } + + startlist <- list(a = min(all_l_rout$y), b = SLOPE, d = max(all_l_rout$y), cs = mean(log(all_l_rout$conc)), r = 0) + + mr <- tryCatch( + { + gsl_nls( + fn = y ~ a + (d - a) / (1 + exp(b * ((cs - r * isSample) - log_dose))), + data = all_l_rout, + start = startlist, # race=T, + control = gsl_nls_control(xtol = 1e-6, ftol = 1e-6, gtol = 1e-6) + ) + }, + warning = function(e) { + mr <<- "In nlsModel singular gradient matrix" + + }) + PAR <- summary(mr)$coefficients[,1] + ROUTplot <- outsPlot_FUN(all_l_rout, OUTs_, TS=all_l_rout$sample[13], PROC="ROUT",PROBE=input$sheet,Q=input$Qslider/100,par=PAR) + output$OutlierDF <- renderTable({ + OUTs_ + }) + output$OutlierPlot <- renderPlot({ + print(ROUTplot) + }) + } + + + ##### XLSX eval ---- #if (CORro < 0) SLOPE <- -1 else SLOPE <- 1 FITs <- Fitting_FUNC(XLdat2, TransFlag = FALSE, nameWS="")