fix #2 - use double colon operator

CytoSpill.Rproj

@@ -18,3 +18,4 @@ StripTrailingWhitespace: Yes
 BuildType: Package
 PackageUseDevtools: Yes
 PackageInstallArgs: --no-multiarch --with-keep.source
+PackageRoxygenize: rd,collate,namespace

DESCRIPTION

@@ -10,9 +10,7 @@ Encoding: UTF-8
 LazyData: true
 biocViews:
 Imports:
-    flowCore,
-    CATALYST,
     nloptr,
     nnls,
     flexmix
-RoxygenNote: 6.1.1
+RoxygenNote: 7.1.1

NAMESPACE

@@ -1 +1,4 @@
-exportPattern("^[[:alpha:]]+")
+# Generated by roxygen2: do not edit by hand
+
+export(GetSpillMat)
+export(SpillComp)

R/DeriveCutoffs.R

@@ -1,7 +1,6 @@
 .DeriveCutoffs <- function(data, cols, n, flexrep){
   data_sample <- data[sample(1:nrow(data), n, replace = FALSE), ]
   data_work <- data_sample[,cols]
-  # cutoffs <- .DeriveCutoffsHelper(x = data_work, emt = 0.001, quantile = 0.1)
   cutoffs <- .DeriveCutoffsHelper(x = data_work, quantile = 0.1, flexrep = flexrep)
   return(cutoffs)
 }
@@ -12,8 +11,6 @@
   for (i in 1:dim(x)[2]) {
     s <- x[, i][which(x[, i] > 0)]
 
-    # s <- asinh(s/5)
-
     model_fit <- flexmix::initFlexmix(s~1, k = c(1,2,3,4,5), nrep = flexrep)
     model_chosen <- flexmix::getModel(model_fit, which = "ICL")
 
@@ -50,16 +47,6 @@
         func_list = list(d1,d2,d3,d4,d5)
       }
 
-
-      # i_min <- which.min(param[1,])
-      # sum_func <- function (x){sum(sapply(func_list, function(f) f(x)))}
-      # type1 <- 0.05
-      # min_modal <- eval(parse(text=paste("d",i_min,sep = "")))
-      # p <- function(x) {
-      #   min_modal(x)/(sum_func(x))-1+type1
-      # }
-      # interval <- c(0, min(param[1,][param[1,]!=min(param[1,])]))
-
       i_min <- which.min(param[1,])
       i_max <- which.max(param[1,])
       min_modal <- eval(parse(text=paste("d",i_min,sep = "")))
@@ -73,7 +60,6 @@
       interval <- c(0, max(param[1,]))
 
       if (p(0)>0){
-      # if (p(0) < (1-type1)){
         root <- uniroot(p,interval)$root[1]
         cutoffs[i] <- root
       } else{
@@ -84,111 +70,7 @@
       cutoffs[i] <- quantile(s, probs = quantile)
       print("using quantile")
     }
-    # cutoffs[i] <- sinh(cutoffs[i])*5
     print(cutoffs[i])
   }
   return(cutoffs)
 }
-#
-# .DeriveCutoffsHelper <- function(x, quantile){
-#   cutoffs <- rep(NA, dim(x)[2])
-#
-#   for (i in 1:dim(x)[2]) {
-#     s <- x[, i][which(x[, i] > 0)]
-#     s_asinh <- asinh(s / 5)
-#     q = quantile(s_asinh, 0.999)
-#     s_asinh[s_asinh > q] = q
-#     model_fit <- flexmix::initFlexmix(s_asinh~1, k = c(1,2,3,4,5))
-#     model_chosen <- flexmix::getModel(model_fit, which = "ICL")
-#
-#     if (length(model_chosen@prior)>1){
-#       param_list <- list()
-#       param <- flexmix::parameters(model_chosen)
-#       func_list <- list()
-#       #get parameters for each component, prior probability, mean, sd
-#       k <- length(model_chosen@prior)
-#       if (k == 2){
-#         d1 <- function(x) model_chosen@prior[1]*dnorm(x, param[1,1], param[2,1])
-#         d2 <- function(x) model_chosen@prior[2]*dnorm(x, param[1,2], param[2,2])
-#         func_list = list(d1,d2)
-#       }
-#       if (k == 3){
-#         d1 <- function(x) model_chosen@prior[1]*dnorm(x, param[1,1], param[2,1])
-#         d2 <- function(x) model_chosen@prior[2]*dnorm(x, param[1,2], param[2,2])
-#         d3 <- function(x) model_chosen@prior[3]*dnorm(x, param[1,3], param[2,3])
-#         func_list = list(d1,d2,d3)
-#       }
-#       if (k == 4){
-#         d1 <- function(x) model_chosen@prior[1]*dnorm(x, param[1,1], param[2,1])
-#         d2 <- function(x) model_chosen@prior[2]*dnorm(x, param[1,2], param[2,2])
-#         d3 <- function(x) model_chosen@prior[3]*dnorm(x, param[1,3], param[2,3])
-#         d4 <- function(x) model_chosen@prior[4]*dnorm(x, param[1,4], param[2,4])
-#         func_list = list(d1,d2,d3,d4)
-#       }
-#       if (k == 5){
-#         d1 <- function(x) model_chosen@prior[1]*dnorm(x, param[1,1], param[2,1])
-#         d2 <- function(x) model_chosen@prior[2]*dnorm(x, param[1,2], param[2,2])
-#         d3 <- function(x) model_chosen@prior[3]*dnorm(x, param[1,3], param[2,3])
-#         d4 <- function(x) model_chosen@prior[4]*dnorm(x, param[1,4], param[2,4])
-#         d5 <- function(x) model_chosen@prior[5]*dnorm(x, param[1,5], param[2,5])
-#         func_list = list(d1,d2,d3,d4,d5)
-#       }
-#       # for (i in seq_along(model_chosen@prior)){
-#       #   #prior, mean, sd
-#       #   param_list[[i]] <- c(model_chosen@prior[i], param[1,i], param[2,i])
-#       #   var <- paste("d", i, sep = "")
-#       #   assign(var, (function(x) model_chosen@prior[i]*dnorm(x, param[1,i], param[2,i])))
-#       #   # assign(var, (function(x) substitute(model_chosen@prior[i]*dnorm(x, param[1,i], param[2,i]), list(i=i))))
-#       #   func_list[[i]] <- eval(parse(text=var))
-#       # }
-#       i_min <- which.min(param[1,])
-#       sum_func <- function (x){sum(sapply(func_list, function(f) f(x)))}
-#       type1 <- 0.05
-#       min_modal <- eval(parse(text=paste("d",i_min,sep = "")))
-#       p <- function(x) {
-#         min_modal(x)/(sum_func(x))-1+type1
-#       }
-#       interval <- c(0, min(param[1,][param[1,]!=min(param[1,])]))
-#       if (p(0)>0){
-#         root <- uniroot(p,interval)$root
-#         cutoffs[i] <- (sinh(root)) * 5
-#       } else{
-#         cutoffs[i] <- quantile(s, probs = quantile)
-#       }
-#     } else{
-#       cutoffs[i] <- quantile(s, probs = quantile)
-#     }
-#   }
-#   return(cutoffs)
-# }
-
-# .DeriveCutoffsHelper <- function(x, emt, quantile) {
-#   # arguments: x: data matrix, emt: EM criterion
-#   # em, cutoff are the functions from package cutoff which are edited to suit the current R version
-#   cutoffs <- rep(NA, dim(x)[2])
-#   for (i in 1:dim(x)[2]) {
-#     a <- x[, i][which(x[, i] > 0)]
-#     a_asinh <- asinh(a / 5)
-#     if (diptest::dip.test(a_asinh)[2] > 0.05) {
-#       cutoff_channeli <- quantile(a_asinh, probs = quantile)
-#     } else {
-#       fit_channeli <- try(em(a_asinh, "normal", "log-normal", t = 0.001), silent=FALSE)
-#       if (isTRUE(class(fit_channeli) == "try-error")) {
-#         cutoff_channeli <- quantile(a_asinh, probs = quantile)
-#       } else {
-#         cutoff_channeli <- try(cutoff(fit_channeli, t = 0.001, nb = 10, distr = 1, type1 = 0.05, level = 0.95))
-#         if (isTRUE(class(cutoff_channeli) == "try-error")) {
-#           cutoff_channeli <- fit_channeli$param[1] + fit_channeli$param[2]
-#         } else{
-#           cutoff_channeli <- cutoff_channeli[1]
-#         }
-#       }
-#     }
-#     cutoff_channeli <- (sinh(cutoff_channeli)) * 5
-#     if (cutoff_channeli <= min(a)){
-#       cutoff_channeli <- quantile(a[which(a > 0)], probs = quantile)
-#     }
-#     cutoffs[i] <- cutoff_channeli
-#   }
-#   return(cutoffs)
-# }

R/EstimateSpill.R

@@ -1,37 +1,11 @@
 .EstimateSpill <- function(data, cutoffs, cols, upperbound, neighbor) {
   results <- list()
   data <- data[,cols]
-  spill_cols <- .SpillColsData(data, l= CATALYST::isotope_list, nneighbor=neighbor)
+  spill_cols <- .SpillColsData(data, l= CytoSpill::isotope_list, nneighbor=neighbor)
   xcols <- .GetFmla(data, spill_cols = spill_cols)
 
-  # #old method#
-  # for (i in 1:ncol(data)) {
-  #   if (!is.na(xcols[[i]][1])) {
-  #     A = as.matrix(data[which(data[,i] < cutoffs[i]), xcols[[i]]])
-  #     b = data[which(data[,i] < cutoffs[i]), i]
-  #     x0 = runif(ncol(A), min = 0, max = upperbound)
-  #     fn = function(x) {
-  #       vec = A%*%x-b
-  #       norm(vec, type = "2")
-  #     }
-  #     result = try(nloptr::slsqp(x0, fn, lower=rep(0, length(x0)), upper = rep(upperbound, length(x0))))
-  #     if (isTRUE(class(result) == "try-error")) {
-  #       result <- NULL
-  #       xcols[[i]] <- NA
-  #     } else{
-  #       result <- result$par
-  #     }
-  #     results[[i]] <- result
-  #   } else {
-  #     results[[i]] <- NULL
-  #   }
-  # }
-
-
-
   for (i in 1:ncol(data)) {
     if (!is.na(xcols[[i]][1])) {
-      # A = as.matrix(data[which(data[,i] < cutoffs[i]), xcols[[i]]])
 
       A = as.matrix(data[which(data[,i] < cutoffs[i]), c(i,xcols[[i]])])
       for (j in seq_along(xcols[[i]])){
@@ -65,28 +39,8 @@
           result <- result$par
         }
         results[[i]] <- result
-
-        # ##
-        # print(dim(A))
-        # ##
-        # model <- nnls::nnls(A=A,b=b)
-        # ##
         print("sqp modelling done")
-        # ##
-        # print(model$x)
-        # results[[i]] <- model$x
       }
-      # ##
-      # print(dim(A))
-      # ##
-      #
-      # # b = data[which(data[,i] < cutoffs[i]), i]
-      # model <- nnls::nnls(A=A,b=b)
-      # ##
-      # print("nnls modelling done")
-      # ##
-      # print(model$x)
-      # results[[i]] <- model$x
     } else {
       results[[i]] <- NA
       print("no spillover cols for this channel")
@@ -95,84 +49,6 @@
   return(list(results, xcols))
 }
 
-#
-# .EstimateSpill <- function(data, cutoffs, cols, upperbound = 0.1) {
-#   results <- list()
-#   data <- data[,cols]
-#   xcols <- .GetFmla(data, spill_cols = .SpillColsData(data))
-#   # for (i in 1:ncol(data)) {
-#   #   if (!is.na(xcols[[i]][1])) {
-#   #     A = as.matrix(data[which(data[,i] < cutoffs[i]), xcols[[i]]])
-#   #     b = data[which(data[,i] < cutoffs[i]), i]
-#   #     x0 = runif(ncol(A), min = 0, max = upperbound)
-#   #     fn = function(x) {
-#   #       vec = A%*%x-b
-#   #       norm(vec, type = "2")
-#   #     }
-#   #     result = try(nloptr::slsqp(x0, fn, lower=rep(0, length(x0)), upper = rep(upperbound, length(x0))))
-#   #     if (isTRUE(class(result) == "try-error")) {
-#   #       result <- NULL
-#   #       xcols[[i]] <- NA
-#   #     } else{
-#   #       result <- result$par
-#   #     }
-#   #     results[[i]] <- result
-#   #   } else {
-#   #     results[[i]] <- NULL
-#   #   }
-#
-#   for (i in 1:ncol(data)) {
-#     if (!is.na(xcols[[i]][1])) {
-#       # A = as.matrix(data[which(data[,i] < cutoffs[i]), xcols[[i]]])
-#
-#       A = as.matrix(data[which(data[,i] < cutoffs[i]), c(i,xcols[[i]])])
-#       for (j in seq_along(xcols[[i]])){
-#         A = A[which(A[,j+1]>=cutoffs[xcols[[i]][j]]),]
-#         if (is.null(dim(A))) {A = matrix(A, nrow = 1)}
-#       }
-#       ##
-#       print(dim(A))
-#       ##
-#       if (nrow(A) == 0) { #where there is only 1 row remains
-#         results[[i]] <- NULL
-#         print("no datapoint available for this channel")
-#       } else {
-#         if (nrow(A) == 1) {
-#           b = A[,1]
-#           A = matrix(A[,-1], nrow = 1)
-#         } else {
-#           b = A[,1]
-#           A = as.matrix(A[,-1])
-#         }
-#         ##
-#         print(dim(A))
-#         ##
-#         model <- nnls::nnls(A=A,b=b)
-#         ##
-#         print("nnls modelling done")
-#         ##
-#         print(model$x)
-#         results[[i]] <- model$x
-#       }
-#       # ##
-#       # print(dim(A))
-#       # ##
-#       #
-#       # # b = data[which(data[,i] < cutoffs[i]), i]
-#       # model <- nnls::nnls(A=A,b=b)
-#       # ##
-#       # print("nnls modelling done")
-#       # ##
-#       # print(model$x)
-#       # results[[i]] <- model$x
-#     } else {
-#       results[[i]] <- NULL
-#       print("no spillover cols for this channel")
-#       }
-#   }
-#   return(list(results, xcols))
-# }
-
 .GetFmla <- function(data, spill_cols) {
   fmlacols <- list()
   cs <- c(1:ncol(data))
@@ -226,15 +102,6 @@
       spill_cols[[i]] <- unique(c(m1, m2, p1, p2, iso, ox))
     }
   }
-  # #only consider +-1 neighboring channel
-  # for (i in seq_along(ms)) {
-  #   p1  <- m1  <- ox <- iso <- NULL
-  #   if ((ms[i] + 1)  %in% ms) p1 <- which(ms == (ms[i] + 1))
-  #   if ((ms[i] - 1)  %in% ms) m1 <- which(ms == (ms[i] - 1))
-  #   if ((ms[i] + 16) %in% ms) ox <- which(ms == (ms[i] + 16))
-  #   iso <- l[[mets[i]]]
-  #   iso <- which(ms %in% iso[iso != ms[i]])
-  #   spill_cols[[i]] <- unique(c(m1, p1, iso, ox))
-  # }
+
   return(spill_cols)
 }

R/SpillComp.R

@@ -11,5 +11,5 @@ SpillComp <- function(file = NULL, data = NULL, cols, n, output = NULL, threshol
   data[,cols] <- data_compensated
   colnames(data) <- data_colnames
   if (!is.null(output)) {flowCore::write.FCS(flowCore::flowFrame(data), filename = output)}
-  return(list(flowFrame(data), spillmat, cutoffs))
+  return(list(flowCore::flowFrame(data), spillmat, cutoffs))
 }