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
}
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent