kBET doesn't give any output?

[getanid123]

Hi I am trying out kBET for my counts data from alignment for three single cell RNA seq batches with [6,8,17cells] respectively. I used the following Code and it executed perfectly without any errors [ after I troubleshooted previous errors I faced using some answered issues]. But I am not able to generate any output file. Can you please help as I want to know how the batch corrected files can be obtained? Thanks so much for your time and in advance!!

library(kBET)
scdata<-read.csv("batch_merge.txt", skip=1,row.names = 1, header=FALSE, sep='\t')
scdata
df<- t(scdata)
df
batch<-factor (c(1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3))

kBET <- function(
  df, batch, k0 = 8, knn = NULL,
  testSize = NULL, do.pca = TRUE, dim.pca = 50,
  heuristic = TRUE, n_repeat = 100,
  alpha = 0.05, addTest = FALSE,
  verbose = TRUE, plot = TRUE, adapt = TRUE
) {
  #create a subsetting mode:
  #if (is.data.frame(batch) && dim(batch)[2]>1) {
  #  cat('Complex study design detected.\n')
  #  design.names <- colnames(batch)
  #  cat(paste0('Subset by ', design.names[[2]], '\n'))
  #  bio <- unlist(unique(batch[[design.names[2]]]))
  #drop subset batch
  #  new.batch <- base::subset(batch,
  # batch[[2]]== bio[1])[,!design.names %in% design.names[[2]]]
  #  sapply(df[batch[[2]]==bio[1],], kBET, new.batch,
  # k0, knn, testSize, heuristic,n_repeat, alpha, addTest, verbose)

  #}

  #preliminaries:
  dof <- length(unique(batch)) - 1 #degrees of freedom
  if (is.factor(batch)) {
    batch <- droplevels(batch)
  }

  frequencies <- table(batch)/length(batch)
  #get 3 different permutations of the batch label
  batch.shuff <- replicate(3, batch[sample.int(length(batch))])

  class.frequency <- data.frame(class = names(frequencies),
                                freq = as.numeric(frequencies))
  dataset <- df
  dim.dataset <- dim(dataset)
  #check the feasibility of data input
  if (dim.dataset[1] != length(batch) && dim.dataset[2] != length(batch)) {
    stop("Input matrix and batch information do not match. Execution halted.")
  }

  if (dim.dataset[2] == length(batch) && dim.dataset[1] != length(batch)) {
    if (verbose) {
      cat('Input matrix has samples as columns. kBET needs samples as rows. Transposing...\n')
    }
    dataset <- t(dataset)
    dim.dataset <- dim(dataset)
  }
  #check if the dataset is too small per se
  if (dim.dataset[1]<=10){
    if (verbose){
      cat("Your dataset has less than 10 samples. Abort and return NA.\n")
    }
    return(NA)
  }

  stopifnot(class(n_repeat) == 'numeric', n_repeat > 0)

  do_heuristic <- FALSE
  if (is.null(k0) || k0 >= dim.dataset[1]) {
    do_heuristic <- heuristic
    if (!heuristic) {
      #default environment size: quarter the size of the largest batch
      k0 <- floor(mean(class.frequency$freq)*dim.dataset[1]/4)
    } else {
      #default environment size: three quarter the size of the largest batch
      k0 <- floor(mean(class.frequency$freq)*dim.dataset[1]*0.75)
      if (k0 < 10) {
        if (verbose){
          warning(
            "Your dataset has too few samples to run a heuristic.\n",
            "Return NA.\n",
            "Please assign k0 and set heuristic=FALSE."
          )
        }

        return(NA)
      }
    }
    if (verbose) {
      cat(paste0('Initial neighbourhood size is set to ', k0, '.\n'))
    }
  }

  #if k0 was set by the user and is too small & we do not operate on a
  #knn graph, abort
  #the reason is that if we want to test kBET on knn graph data integration methods,
  #we usually face small numbers of nearest neighbours.
  if (k0 < 10 & is.null(knn)) {
    if (verbose){
      warning(
        "Your dataset has too few samples to run a heuristic.\n",
        "Return NA.\n",
        "Please assign k0 and set heuristic=FALSE."
      )
    }
    return(NA)
  }
  # find KNNs
  if (is.null(knn)) {
    if (!do.pca) {
      if (verbose) {
        cat('finding knns...')
        tic <- proc.time()
      }
      #use the nearest neighbour index directly for further use in the package
      knn <- get.knn(dataset, k = k0, algorithm = 'cover_tree')$nn.index
    } else {
      dim.comp <- min(dim.pca, dim.dataset[2])
      if (verbose) {
        cat('reducing dimensions with svd first...\n')
      }
      data.pca <- svd(x = dataset, nu = dim.comp, nv = 0)
      if (verbose) {
        cat('finding knns...')
        tic <- proc.time()
      }
      knn <- get.knn(data.pca$u,  k = k0, algorithm = 'cover_tree')
    }
    if (verbose) {
      cat('done. Time:\n')
      print(proc.time() - tic)
    }
  }

  #backward compatibility for knn-graph
  if (class(knn) =='list'){
    knn <- knn$nn.index
    if (verbose){
      cat('KNN input is a list, extracting nearest neighbour index.\n')
    }
  }

  #set number of tests
  if (is.null(testSize) || (floor(testSize) < 1 || dim.dataset[1] < testSize)) {
    test.frac <- 0.1
    testSize <- ceiling(dim.dataset[1]*test.frac)
    if (testSize < 25 && dim.dataset[1] > 25) {
      testSize <- 25
    }
    if (verbose) {
      cat('Number of kBET tests is set to ')
      cat(paste0(testSize, '.\n'))
    }
  }
  #decide to adapt general frequencies
  if (adapt) {
    # idx.run <- sample.int(dim.dataset[1], size = min(2*testSize, dim.dataset[1]))
    outsider <- which(!(seq_len(dim.dataset[1]) %in% knn[, seq_len(k0 - 1)]))
    is.imbalanced <- FALSE #initialisation
    p.out <- 1
    #avoid unwanted things happen if length(outsider) == 0
    if (length(outsider) > 0) {
      p.out <- chi_batch_test(outsider, class.frequency, batch,  dof)
      if (!is.na(p.out)){
        is.imbalanced <- p.out < alpha
        if (is.imbalanced) {
          new.frequencies <- table(batch[-outsider])/length(batch[-outsider])
          new.class.frequency <- data.frame(class = names(new.frequencies),
                                            freq = as.numeric(new.frequencies))
          if (verbose) {
            outs_percent <- round(length(outsider) / length(batch) * 100, 3)
            cat(paste(
              sprintf(
                'There are %s cells (%s%%) that do not appear in any neighbourhood.',
                length(outsider), outs_percent
              ),
              'The expected frequencies for each category have been adapted.',
              'Cell indexes are saved to result list.',
              '', sep = '\n'
            ))
          }
        } else {
          if (verbose) {
            cat(paste0('No outsiders found.'))
          }
        }
      } else{
        if (verbose) {
          cat(paste0('No outsiders found.'))
        }
      }
    }
  }

  if (do_heuristic) {
    #btw, when we bisect here that returns some interval with
    #the optimal neihbourhood size
    if (verbose) {
      cat('Determining optimal neighbourhood size ...')
    }
    opt.k <- bisect(scan_nb, bounds = c(10, k0), known = NULL, dataset, batch, knn)
    #result
    if (length(opt.k) > 1) {
      k0 <- opt.k[2]
      if (verbose) {
        cat(paste0('done.\nNew size of neighbourhood is set to ', k0, '.\n'))
      }
    } else {
      if (verbose) {
        cat(paste(
          'done.',
          'Heuristic did not change the neighbourhood.',
          sprintf('If results appear inconclusive, change k0=%s.', k0),
          '', sep = '\n'
        ))
      }
    }
  }

  #initialise result list
  rejection <- list()
  rejection$summary <- data.frame(
    kBET.expected = numeric(4),
    kBET.observed = numeric(4),
    kBET.signif = numeric(4)
  )

  rejection$results <- data.frame(
    tested = numeric(dim.dataset[1]),
    kBET.pvalue.test = rep(0,dim.dataset[1]),
    kBET.pvalue.null = rep(0, dim.dataset[1])
  )
  #get average residual score
  env <- as.vector(cbind(knn[, seq_len(k0 - 1)], seq_len(dim.dataset[1])))
  cf <- if (adapt && is.imbalanced) new.class.frequency else class.frequency
  rejection$average.pval <- 1 - pchisq(k0 * residual_score_batch(env, cf, batch), dof)

  #initialise intermediates
  kBET.expected <- numeric(n_repeat)
  kBET.observed <- numeric(n_repeat)
  kBET.signif <- numeric(n_repeat)

  if (addTest) {
    #initialize result list
    rejection$summary$lrt.expected <- numeric(4)
    rejection$summary$lrt.observed <- numeric(4)

    rejection$results$lrt.pvalue.test <- rep(0,dim.dataset[1])
    rejection$results$lrt.pvalue.null <- rep(0, dim.dataset[1])

    lrt.expected <- numeric(n_repeat)
    lrt.observed <- numeric(n_repeat)
    lrt.signif <- numeric(n_repeat)
    #decide to perform exact test or not
    if (choose(k0 + dof,dof) < 5e5 && k0 <= min(table(batch))) {
      exact.expected <- numeric(n_repeat)
      exact.observed <- numeric(n_repeat)
      exact.signif <- numeric(n_repeat)

      rejection$summary$exact.expected <- numeric(4)
      rejection$summary$exact.observed <- numeric(4)
      rejection$results$exact.pvalue.test <- rep(0, dim.dataset[1])
      rejection$results$exact.pvalue.null <- rep(0, dim.dataset[1])
    }

    for (i in seq_len(n_repeat)) {
      # choose a random sample from dataset (rows: samples, columns: features)
      idx.runs <- sample.int(dim.dataset[1], size = testSize)
      env <- cbind(knn[idx.runs, seq_len(k0 - 1)], idx.runs)
      #env.rand <- t(sapply(rep(dim.dataset[1],testSize),  sample.int, k0))

      #perform test
      if (adapt && is.imbalanced) {
        p.val.test <- apply(env, 1, FUN = chi_batch_test,
                            new.class.frequency, batch,  dof)
      } else {
        p.val.test <- apply(env, 1, FUN = chi_batch_test,
                            class.frequency, batch,  dof)
      }

      is.rejected <- p.val.test < alpha

      #p.val.test <- apply(env, 1, FUN = chi_batch_test, class.frequency,
      #batch,  dof)
      p.val.test.null <-  apply(apply(batch.shuff, 2,
                                      function(x, freq, dof, envir) {
                                        apply(envir, 1, FUN = chi_batch_test, freq, x, dof)},
                                      class.frequency, dof, env), 1, mean, na.rm=TRUE)
      #p.val.test.null <- apply(env.rand, 1, FUN = chi_batch_test,
      #class.frequency, batch, dof)

      #summarise test results
      kBET.expected[i] <- sum(p.val.test.null < alpha,
                              na.rm=TRUE) / sum(!is.na(p.val.test.null))
      kBET.observed[i] <- sum(is.rejected,na.rm=TRUE) / sum(!is.na(p.val.test))

      #compute significance
      kBET.signif[i] <-
        1 - ptnorm(kBET.observed[i],
                   mu = kBET.expected[i],
                   sd = sqrt(kBET.expected[i] * (1 - kBET.expected[i]) / testSize),
                   alpha = alpha)

      #assign results to result table
      rejection$results$tested[idx.runs] <- 1
      rejection$results$kBET.pvalue.test[idx.runs] <- p.val.test
      rejection$results$kBET.pvalue.null[idx.runs] <- p.val.test.null

      #compute likelihood-ratio test (approximation for multinomial exact test)
      cf <- if (adapt && is.imbalanced) new.class.frequency else class.frequency
      p.val.test.lrt <- apply(env, 1, FUN = lrt_approximation, cf, batch, dof)
      p.val.test.lrt.null <- apply(apply(batch.shuff, 2,
                                         function(x, freq, dof, envir) {
                                           apply(envir, 1, FUN = lrt_approximation, freq, x, dof)},
                                         class.frequency, dof, env), 1, mean, na.rm=TRUE)

      lrt.expected[i] <-
        sum(p.val.test.lrt.null < alpha, na.rm=TRUE) / sum(!is.na(p.val.test.lrt.null))
      lrt.observed[i] <-
        sum(p.val.test.lrt < alpha, na.rm=TRUE) / sum(!is.na(p.val.test.lrt))

      lrt.signif[i] <-
        1 - ptnorm(lrt.observed[i],
                   mu = lrt.expected[i],
                   sd = sqrt(lrt.expected[i] * (1 - lrt.expected[i]) / testSize),
                   alpha = alpha)

      rejection$results$lrt.pvalue.test[idx.runs] <- p.val.test.lrt
      rejection$results$lrt.pvalue.null[idx.runs] <- p.val.test.lrt.null

      #exact result test consume a fairly high amount of computation time,
      #as the exact test computes the probability of ALL
      #possible configurations (under the assumption that all batches
      #are large enough to 'imitate' sampling with replacement)
      #For example: k0=33 and dof=5 yields 501942 possible
      #choices and a computation time of several seconds (on a 16GB RAM machine)
      if (exists(x = 'exact.observed')) {
        if (adapt && is.imbalanced) {
          p.val.test.exact <- apply(env, 1, multiNom,
                                    new.class.frequency$freq, batch)
        } else {
          p.val.test.exact <- apply(env, 1, multiNom,
                                    class.frequency$freq, batch)
        }
        p.val.test.exact.null <-  apply(apply(batch.shuff, 2,
                                              function(x, freq, envir) {
                                                apply(envir, 1, FUN = multiNom, freq, x)},
                                              class.frequency$freq, env), 1, mean, na.rm=TRUE)
        # apply(env, 1, multiNom, class.frequency$freq, batch.shuff)

        exact.expected[i] <- sum(p.val.test.exact.null < alpha, na.rm=TRUE)/testSize
        exact.observed[i] <- sum(p.val.test.exact < alpha, na.rm=TRUE)/testSize
        #compute the significance level for the number of rejected data points
        exact.signif[i] <-
          1 - ptnorm(exact.observed[i],
                     mu = exact.expected[i],
                     sd = sqrt(exact.expected[i] * (1 - exact.expected[i]) / testSize),
                     alpha = alpha)
        #p-value distribution
        rejection$results$exact.pvalue.test[idx.runs] <- p.val.test.exact
        rejection$results$exact.pvalue.null[idx.runs] <- p.val.test.exact.null
      }
    }

    if (n_repeat > 1) {
      #summarize chi2-results
      CI95 <- c(0.025,0.5,0.975)
      rejection$summary$kBET.expected <-  c(mean(kBET.expected, na.rm=TRUE) ,
                                            quantile(kBET.expected, CI95,
                                                     na.rm=TRUE))
      rownames(rejection$summary) <- c('mean', '2.5%', '50%', '97.5%')
      rejection$summary$kBET.observed <-  c(mean(kBET.observed,na.rm=TRUE) ,
                                            quantile(kBET.observed, CI95,
                                                     na.rm=TRUE))
      rejection$summary$kBET.signif <- c(mean(kBET.signif,na.rm=TRUE) ,
                                         quantile(kBET.signif, CI95,
                                                  na.rm=TRUE))
      #summarize lrt-results
      rejection$summary$lrt.expected <-  c(mean(lrt.expected,na.rm=TRUE) ,
                                           quantile(lrt.expected, CI95,
                                                    na.rm=TRUE))
      rejection$summary$lrt.observed <-  c(mean(lrt.observed,na.rm=TRUE) ,
                                           quantile(lrt.observed, CI95,
                                                    na.rm=TRUE))
      rejection$summary$lrt.signif <- c(mean(lrt.signif,na.rm=TRUE) ,
                                        quantile(lrt.signif, CI95,
                                                 na.rm=TRUE))
      #summarize exact test results
      if (exists(x = 'exact.observed')) {
        rejection$summary$exact.expected <-  c(mean(exact.expected,na.rm=TRUE) ,
                                               quantile(exact.expected, CI95,
                                                        na.rm=TRUE))
        rejection$summary$exact.observed <-  c(mean(exact.observed,na.rm=TRUE) ,
                                               quantile(exact.observed, CI95,
                                                        na.rm=TRUE))
        rejection$summary$exact.signif <- c(mean(exact.signif,na.rm=TRUE) ,
                                            quantile(exact.signif, CI95,
                                                     na.rm=TRUE))
      }

      if (n_repeat < 10) {
        cat('Warning: The quantile computation for ')
        cat(paste0(n_repeat))
        cat(' subset results is not meaningful.')
      }

      if (plot && exists(x = 'exact.observed')) {
        plot.data <- data.frame(class = rep(c('kBET', 'kBET (random)',
                                              'lrt', 'lrt (random)',
                                              'exact', 'exact (random)'),
                                            each = n_repeat),
                                data =  c(kBET.observed, kBET.expected,
                                          lrt.observed, lrt.expected,
                                          exact.observed, exact.expected))
        g <- ggplot(plot.data, aes(class, data)) + geom_boxplot() +
          theme_bw() + labs(x = 'Test', y = 'Rejection rate')  +
          theme(axis.text.x = element_text(angle = 45, hjust = 1))
        print(g)
      }
      if (plot && !exists(x = 'exact.observed')) {
        plot.data <- data.frame(class = rep(c('kBET', 'kBET (random)',
                                              'lrt', 'lrt (random)'),
                                            each = n_repeat),
                                data =  c(kBET.observed, kBET.expected,
                                          lrt.observed, lrt.expected))
        g <- ggplot(plot.data, aes(class, data)) + geom_boxplot() +
          theme_bw() + labs(x = 'Test', y  = 'Rejection rate') +
          scale_y_continuous(limits = c(0,1))  +
          theme(axis.text.x = element_text(angle = 45, hjust = 1))
        print(g)
      }

    } else {#i.e. no n_repeat
      rejection$summary$kBET.expected <- kBET.expected
      rejection$summary$kBET.observed <- kBET.observed
      rejection$summary$kBET.signif <- kBET.signif

      if (addTest) {
        rejection$summary$lrt.expected <- lrt.expected
        rejection$summary$lrt.observed <- lrt.observed
        rejection$summary$lrt.signif <- lrt.signif
        if (exists(x = 'exact.observed')) {
          rejection$summary$exact.expected <-  exact.expected
          rejection$summary$exact.observed <-  exact.observed
          rejection$summary$exact.signif <- exact.signif
        }
      }

    }
  } else {#kBET only
    for (i in seq_len(n_repeat)) {
      # choose a random sample from dataset
      #(rows: samples, columns: parameters)
      idx.runs <- sample.int(dim.dataset[1], size = testSize)
      env <- cbind(knn[idx.runs,seq_len(k0 - 1)], idx.runs)

      #perform test
      cf <- if (adapt && is.imbalanced) new.class.frequency else class.frequency
      p.val.test <- apply(env, 1, chi_batch_test, cf, batch, dof)

      #print(dim(env))
      is.rejected <- p.val.test < alpha

      p.val.test.null <- apply(
        batch.shuff, 2,
        function(x) apply(env, 1, chi_batch_test, class.frequency, x, dof)
      )
      # p.val.test.null <- apply(env, 1, FUN = chi_batch_test,
      #class.frequency, batch.shuff, dof)

      #summarise test results
      #kBET.expected[i] <- sum(p.val.test.null < alpha) / length(p.val.test.null)
      kBET.expected[i] <- mean(apply(
        p.val.test.null, 2,
        function(x) sum(x < alpha, na.rm =TRUE) / sum(!is.na(x))
      ))

      kBET.observed[i] <- sum(is.rejected,na.rm=TRUE) / sum(!is.na(p.val.test))

      #compute significance
      kBET.signif[i] <- 1 - ptnorm(
        kBET.observed[i],
        mu = kBET.expected[i],
        sd = sqrt(kBET.expected[i] * (1 - kBET.expected[i]) / testSize),
        alpha = alpha
      )
      #assign results to result table
      rejection$results$tested[idx.runs] <- 1
      rejection$results$kBET.pvalue.test[idx.runs] <- p.val.test
      rejection$results$kBET.pvalue.null[idx.runs] <- rowMeans(p.val.test.null, na.rm=TRUE)
    }

    if (n_repeat > 1) {
      #summarize chi2-results
      CI95 <- c(0.025,0.5,0.975)
      rejection$summary$kBET.expected <- c(mean(kBET.expected,na.rm=TRUE),
                                           quantile(kBET.expected,
                                                    CI95,na.rm=TRUE))
      rownames(rejection$summary) <- c('mean', '2.5%', '50%', '97.5%')
      rejection$summary$kBET.observed <- c(mean(kBET.observed,na.rm=TRUE),
                                           quantile(kBET.observed, CI95,
                                                    na.rm=TRUE))
      rejection$summary$kBET.signif <- c(mean(kBET.signif,na.rm=TRUE),
                                         quantile(kBET.signif, CI95,
                                                  na.rm=TRUE))

      #return also n_repeat
      rejection$stats$kBET.expected <- kBET.expected
      rejection$stats$kBET.observed <- kBET.observed
      rejection$stats$kBET.signif <- kBET.signif

      if (n_repeat < 10) {
        cat('Warning: The quantile computation for ')
        cat(paste0(n_repeat))
        cat(' subset results is not meaningful.')
      }
      if (plot) {
        plot.data <- data.frame(class = rep(c('observed(kBET)',
                                              'expected(random)'),
                                            each = n_repeat),
                                data =  c( kBET.observed,kBET.expected))
        g <- ggplot(plot.data, aes(class, data)) +
          geom_boxplot() + theme_bw() +
          labs(x = 'Test', y = 'Rejection rate') +
          scale_y_continuous(limits = c(0,1))
        print(g)
      }
    } else {
      rejection$summary$kBET.expected <- kBET.expected[1]
      rejection$summary$kBET.observed <- kBET.observed[1]
      rejection$summary$kBET.signif <- kBET.signif[1]
    }
  }
  #collect parameters
  rejection$params <- list()
  rejection$params$k0 <- k0
  rejection$params$testSize <- testSize
  rejection$params$do.pca <- do.pca
  rejection$params$dim.pca <- dim.pca
  rejection$params$heuristic <- heuristic
  rejection$params$n_repeat <- n_repeat
  rejection$params$alpha <- alpha
  rejection$params$addTest <- addTest
  rejection$params$verbose <- verbose
  rejection$params$plot <- plot

  #add outsiders
  if (adapt) {
    rejection$outsider <- list()
    rejection$outsider$index <- outsider
    rejection$outsider$categories <- table(batch[outsider])
    rejection$outsider$p.val <- p.out
  }
  rejection
}

session_info_kbet.docx

[mbuttner]

Hi, thanks for your interest in kBET. For troubleshooting, you can set verbose=TRUE. Your problem is that you have very few cells and thus, the neighbourhood size k0 is lower than 10. I introduced a check of the input neighbourhood size:

if (k0 < 10) {
        if (verbose){
          warning(
            "Your dataset has too few samples to run a heuristic.\n",
            "Return NA.\n",
            "Please assign k0 and set heuristic=FALSE."
          )
        }
 return(NA)

Actually, I test for it twice:

if (k0 < 10 & is.null(knn)) {
    if (verbose){
      warning(
        "Your dataset has too few samples to run a heuristic.\n",
        "Return NA.\n",
        "Please assign k0 and set heuristic=FALSE."
      )
    }
    return(NA)
  }

You may adapt the function such that it ignores the small k0, but I am not sure if the results are reasonable.

[getanid123]

Hi mbuttner,

Thank you so much for your help and insight.

As per your suggestion, I changed k0 to 10 and the program ran like before without errors, but there was no output. The verbose=TRUE is already set, but I do not see any messages when I run in R?

Thanks for your time in advance!!

[mbuttner]

Hi, what does the call test_result <- kBET(df=df, batch=batch) return?

[getanid123]

Hi

Thank you so much for getting back and sorry for the delay in responding. I tried the call you mentioned and another call using FNN library, Please find the details below:

library(devtools) install_github('theislab/kBET',force = TRUE) library(kBET) scdata<-read.csv("counts.txt", skip=1,row.names = 1, header=FALSE, sep='\t') df<- t(scdata) batch<-factor (c(1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3)) test_result <- kBET(df=df, batch=batch) test_result --->[the result was NA]

When I tried along with FNN library I got the result below: library(devtools) install_github('theislab/kBET',force = TRUE) scdata<-read.csv("counts.txt", skip=1,row.names = 1, header=FALSE, sep='\t') df<- t(scdata) batch<-factor (c(1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3)) library('FNN') k0=6 knn <- get.knn(df, k=k0, algorithm = 'cover_tree') batch.estimate <- kBET(df, batch, k = k0, knn = knn, heuristic = FALSE) ---> I got a Rejection rate plot. Rplot.pdf

[mbuttner]

OK, the batch.estimate variable contains the result of kBET that you are looking for.

[getanid123]

Yes, got it. Thank you!!