Cleaned coda_new

R/cell_density.R

@@ -58,7 +58,7 @@ estimateCellDensityKde <- function(emb, sample.per.cell, sample.groups, bins, ba
 #' @param n.cores number of cores
 #' @param m numeric Maximum order of Chebyshev coeff to compute (default=50)
 #' @keywords internal
-estimateCellDensityGraph <- function(graph, sample.per.cell, sample.groups, n.cores=1, beta=30, m=50, verbose = TRUE) {
+estimateCellDensityGraph <- function(graph, sample.per.cell, sample.groups, n.cores=1, beta=30, m=50, verbose=TRUE) {
   sig.mat <- unique(sample.per.cell) %>% sapply(function(s) as.numeric(sample.per.cell == s)) %>%
     set_rownames(names(sample.per.cell)) %>% set_colnames(unique(sample.per.cell))
 
@@ -217,7 +217,7 @@ adjustZScoresByPermutations <- function(score, scores.shuffled, wins=0.01, smoot
     if (is.null(graph)) stop("graph has to be provided if smooth is TRUE")
 
     g.filt <- function(...) sccore::heatFilter(..., beta=beta)
-    score %<>% smoothSignalOnGraph(filter=g.filt, graph=graph, l.max=l.max)
+    score %<>% smoothSignalOnGraph(filter=g.filt, graph=graph, l.max=l.max, progress=FALSE)
     scores.shuffled %<>%
       smoothSignalOnGraph(filter=g.filt, graph=graph, n.cores=n.cores, l.max=l.max, progress=verbose) %>%
       as.matrix()

R/coda_new.R → R/coda.R

@@ -185,7 +185,6 @@ produceResampling <- function(cnts, groups, n.perm = 1000, seed = 239) {
 
 #' @keywords internal
 runCoda <- function(cnts, groups, n.seed=239, n.boot=1000, ref.cell.type=NULL, null.distr=FALSE) {
-
   # Create datasets as
   samples.init <- produceResampling(cnts = cnts, groups = groups, n.perm = n.boot, seed = n.seed)
   loadings <- do.call(cbind, lapply(1:length(samples.init$cnts), function(ib) {
@@ -333,7 +332,6 @@ referenceSet <- function(freqs, groups, p.thresh=0.05) {
               mx.final = mx))
 }
 
-
 #' @keywords internal
 sbpDiff <- function(freqs, groups){
   freqs[freqs == 0] <- min(freqs[freqs != 0])/2
@@ -373,99 +371,3 @@ sbpDiff <- function(freqs, groups){
 
   return(t(sbp))
 }
-
-#' @keywords internal
-pvalInLoadingsOrder <- function(cell.types.order, cnts, groups){
-  p <- c()
-  cnts[cnts <= 0] <- 0.5
-  freqs <- (cnts)/rowSums(cnts)
-  for(i in 0:(ncol(freqs) - 2)){
-
-    if(i == 0){
-      freqs1 <- freqs
-    } else {
-      freqs1 <- freqs[, -which(colnames(freqs) %in% cell.types.order[1:i] )]
-    }
-
-    psi <- coda.base::ilr_basis(ncol(freqs1), type = "default")
-    rownames(psi) <- colnames(freqs1)
-    b <- log(freqs1) %*% psi
-
-    x <- summary(lm(groups ~ b))
-    pval.lm <- pf(x$fstatistic[1],x$fstatistic[2],x$fstatistic[3],lower.tail=FALSE)
-    p <- c(p, pval.lm)
-
-  }
-  p <- c(p, 1)
-  names(p) <- cell.types.order
-  return(p)
-}
-
-#' @keywords internal
-estimateCdaSpaceNew <- function(cnts, groups){
-  # ----- Get data -----
-
-  cnts[cnts == 0] = 0.5
-  freqs <- cnts / rowSums(cnts)
-  psi <- coda.base::ilr_basis(ncol(freqs), type = "default")
-  rownames(psi) <- colnames(cnts)
-  b <- log(freqs) %*% psi
-
-  f1 = apply(exp(b %*% t(psi)),1, function(x) {x/sum(x)})
-
-  # PCA - to get new ilr basis of principal components
-  b.norm <-  apply(b, 2, function(y) y - mean(y))
-  pca.res <- prcomp(b.norm)
-  pca.loadings <- psi %*% pca.res$rotation
-
-  psi <- psi %*% pca.res$rotation  # new Psi matrix
-
-  b <- log(freqs) %*% psi  # New ilr coordinates
-  b.init <- log(freqs.init) %*% psi
-
-
-
-  # Standardization
-  b.mean <- colMeans(b)
-  b.sd <- apply(b, 2, sd)
-  b.init <- (b.init - b.mean) / b.sd
-  b <- apply(b, 2, scale)
-
-  b.plot <- rbind(b, b.init)
-
-  # ------------------------------
-  # Find the first projection
-
-  # Create dataframe
-  b.df <- data.frame(b)
-  b.df$groups <- 1*groups
-
-  # Linear model
-  res.lm <- lm(groups ~ ., data = b.df)
-  w <- res.lm$coefficients[-1]  # coefficients without intercept
-
-  w <- as.matrix(w / sqrt(sum(w ^ 2)))  # not required
-
-
-  # Projections of b.plot on w - is the first coordinate
-  b.proj = as.matrix(apply(b.plot, 1, function(x) sum(x*w)))
-  b.proj.vec = b.proj %*% t(w)
-  b.ort <- b.plot - b.proj.vec
-
-  b.ort.proj = as.matrix(apply(b.ort, 1, function(x) sum(x*w))) # <- check that all are small
-  if(sum(b.ort.proj > 10^(-10)) != 0) stop('Something is going wrong')
-
-  # ------------------------------
-  # Find the second projection - from PCA,
-  # because the all differences between groups are in the first coordinate
-
-  pca.ort <- prcomp(b.ort)
-
-
-  b.plot <- data.frame(cbind(b.proj, pca.ort$x[,1]))
-  colnames(b.plot) <- c('S1', 'S2')
-  # b.plot$group <- c(as.character(uniq.combo$site), 'init')
-  # colnames(b.plot)
-
-  return(b.plot)
-}