From 0ffcc55deb476a2235e5c422d673674939c346db Mon Sep 17 00:00:00 2001 From: PratibhaPanwar Date: Tue, 27 Aug 2024 11:02:40 +1000 Subject: [PATCH] Modified adative smoothing. Entropy values are used to generate a matrix of weights, where each column contains weights corresponding to a unique entropy value. --- R/adaptiveSmoothing.R | 72 +++++++++++++++++++++++++++++++++---------- R/expKernel.R | 8 +++-- R/gaussKernel.R | 8 +++-- 3 files changed, 68 insertions(+), 20 deletions(-) diff --git a/R/adaptiveSmoothing.R b/R/adaptiveSmoothing.R index 3c3ccdd..5c20671 100644 --- a/R/adaptiveSmoothing.R +++ b/R/adaptiveSmoothing.R @@ -22,34 +22,38 @@ #### Smoothing adaptiveSmoothing <- function(spe, nnCells, NN, kernel, spread) { - e <- unique(spe$entropy) - gXc <- as(logcounts(spe), "sparseMatrix") + ed = unique(spe$entropy) + gXc = as(logcounts(spe), "sparseMatrix") + if (kernel == "G") { + # normal distribution weights + # for each unique entropy value + weights = .gauss_kernel(ed, NN + 1, spread) + } else if (kernel == "E") { + # exponential distribution weights + # for each unique entropy value + weights = .exp_kernel(ed, NN + 1, spread) + } + outMats = matrix(nrow = nrow(spe), ncol = 0) - for (val in c(1:length(e))) { + for (val in 1:length(ed)) { # cells with same entropy value will have same weights - entCells <- colnames(spe[, spe$entropy == e[val]]) - inMatList = list() - if (kernel == "G") { - # normal distribution-weights - weight <- .gauss_kernel(e[val], NN + 1, spread) - } else if (kernel == "E") { - # exponential distribution-weights - weight <- .exp_kernel(e[val], NN + 1, spread) - } + entCells <- colnames(spe[, spe$entropy == ed[val]]) + e = paste0("E", ed[val]) if (length(entCells) == 1) { region <- as.vector(nnCells[entCells, ]) inMat <- as.matrix(gXc[, region]) - tmpMat <- .smoothedData(inMat, weight) + tmpMat <- .smoothedData(inMat, weights[, e]) colnames(tmpMat) <- entCells } else { + inMatList = list() for (x in entCells) { # names of central cell + NN cells region <- as.vector(nnCells[x, ]) # gene expression matrix of NN cells inMatList[[x]] <- as.matrix(gXc[, region]) } - out <- lapply(inMatList, .smoothedData, weight = weight) + out <- lapply(inMatList, .smoothedData, weight = weights[, e]) tmpMat <- matrix(unlist(out), ncol = length(out), dimnames = list(rownames(out[[1]]), sapply(out, names))) @@ -57,7 +61,6 @@ adaptiveSmoothing <- function(spe, nnCells, NN, kernel, spread) { } outMats <- cbind(outMats, tmpMat) } - # rearrange the cells in smoothed data matrix smoothMat <- outMats[, as.vector(colnames(spe))] @@ -70,8 +73,45 @@ adaptiveSmoothing <- function(spe, nnCells, NN, kernel, spread) { stop("Missing values in smoothed data.") } else { # add data to spatial experiment - assay(spe, i = "smoothed") <- as(smoothMat, "sparseMatrix") + assay(spe, "smoothed") <- as(smoothMat, "sparseMatrix") print(paste("Smoothing performed. NN =", NN, "Kernel =", kernel, "Spread =", spread, Sys.time())) } return(spe) } + +# cl <- makeCluster(threads) +# doParallel::registerDoParallel(cl) +# y <- foreach (x = c(1:ncol(gXc))) %dopar% { +# # index cell name +# cell = colnames(gXc)[x] +# # cell entropy +# e = paste0("E", spe$entropy[x]) +# # names of index cell + NN neighbors +# region = as.vector(nnCells[x, ]) +# # gene expression matrix of NN cells +# inMat = as.matrix(gXc[, region]) +# smat <- .smoothedData(inMat, weights[, e]) +# colnames(smat) <- cell +# smat +# } +# stopCluster(cl) +# +# # QC check +# check.cells = identical(sapply(y, colnames), as.vector(spe[[cells]])) +# check.genes = identical(rownames(logcounts(spe)), rownames(y[[1]])) +# check.NA = sum(is.na(unlist(y))) +# if (check.cells == FALSE) { +# stop("Order of cells in smoothed data does not match cell order in spatial experiment object.") +# } else if (check.genes == FALSE) { +# stop("Order of genes in smoothed data does not match gene order in spatial experiment object.") +# } else if (check.NA != 0) { +# stop("Missing values in smoothed data.") +# } else { +# smoothMat = matrix(unlist(y), ncol = length(y), +# dimnames = list(rownames(y[[1]]), sapply(y, colnames))) +# # add data to spatial experiment +# assay(spe, i = "smoothed") <- as(smoothMat, "sparseMatrix") +# print(paste("Smoothing performed. NN =", NN, "Kernel =", kernel, "Spread =", spread, Sys.time())) +# } + + diff --git a/R/expKernel.R b/R/expKernel.R index 5bdee04..fae2856 100644 --- a/R/expKernel.R +++ b/R/expKernel.R @@ -14,6 +14,10 @@ # function to retrieve exponential distribution weights for neighbors .exp_kernel <- function(ed, NN, rate) { # distribution from 0 to entropy, with cells in smoothing radius as cut points - cutpoints <- seq(0, ed, length.out = NN) - return(as.matrix(dexp(cutpoints, rate = rate))) + + # cutpoints <- seq(0, ed, length.out = NN) + # return(as.matrix(dexp(cutpoints, rate = rate))) + weights = dexp(sapply(ed, function(x) seq(0, x, length.out = NN)), rate = rate) + colnames(weights) <- paste0("E", ed) + return (weights) } diff --git a/R/gaussKernel.R b/R/gaussKernel.R index 1630209..37d1081 100644 --- a/R/gaussKernel.R +++ b/R/gaussKernel.R @@ -14,6 +14,10 @@ # function to retrieve normal distribution weights for neighbors .gauss_kernel <- function(ed, NN, sd) { # distribution from 0 to entropy, with cells in smoothing radius as cut points - cutpoints <- seq(0, ed, length.out = NN) - return(as.matrix(dnorm(cutpoints, sd = sd))) + + # cutpoints <- seq(0, ed, length.out = NN) + # return(as.matrix(dnorm(cutpoints, sd = sd))) + weights = dnorm(sapply(ed, function(x) seq(0, x, length.out = NN)), sd = sd) + colnames(weights) <- paste0("E", ed) + return (weights) }