diff --git a/DESCRIPTION b/DESCRIPTION index 2458236..114eb2b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: jazzPanda Type: Package Title: Finding spatially relevant marker genes in image based spatial transcriptomics data -Version: 0.99.7 +Version: 0.2.3 Date: 2024-07-25 Authors@R: c(person(given = "Melody", @@ -11,7 +11,7 @@ Authors@R: comment=c(ORCID="0000-0002-2222-0958")) ) LazyData: FALSE -Depends: R (>= 4.4.0) +Depends: R (>= 4.5.0) Imports: spatstat.geom, dplyr, glmnet, caret, foreach, stats, magrittr, doParallel, BiocParallel, methods, BumpyMatrix,SpatialExperiment VignetteBuilder: knitr Suggests: diff --git a/NEWS.md b/NEWS.md index c9605ec..f660199 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,28 @@ +# jazzPanda 0.2.3 +*Updated on June 24, 2025* + +Changes to the package include: + +- **function `get_vectors()`**: + - Parameters `w_x` and `w_y` have been removed. Spatial windows are now + automatically determined **per sample** using a 5% buffer around the + coordinate range. + - Improved handling of genes that are only detected in a subset of samples. + If a gene is not present in a sample, a zero-filled vector is used for that + sample to maintain consistent vector lengths across samples. + +- **function `create_genesets()`**: + - Parameters `w_x` and `w_y` have been removed. Spatial windows are now + computed **per sample** using a 5% buffer based on input coordinates. + - Improved handling of genes that are only detected in a subset of samples. + If a gene is not present in a sample, a zero-filled vector is used for that + sample to maintain consistent vector lengths across samples. + +- **function `compute_permp()`**: + - Parameters `w_x` and `w_y` are no longer required. Spatial boundaries + are now automatically derived **per sample** using a 5% buffer. + + # jazzPanda 0.2.2 *Updated on February 23, 2025* diff --git a/R/compute_permp.R b/R/compute_permp.R index a1dcc3c..3a48055 100644 --- a/R/compute_permp.R +++ b/R/compute_permp.R @@ -27,10 +27,6 @@ #' two giving the numbers of rectangular quadrats in the x and y directions. If #' the \code{bin_type} is "hexagonal", this will be a number giving the side #' length of hexagons. Positive numbers only. -#' @param w_x a numeric vector of length two specifying the x coordinate -#' limits of enclosing box. -#' @param w_y a numeric vector of length two specifying the y coordinate -#' limits of enclosing box. #' @param use_cm A boolean value that specifies whether to create spatial #' vectors for genes using the count matrix and cell coordinates instead of #' the transcript coordinates when both types of information are available. @@ -46,39 +42,13 @@ #' \item{\code{gene_mt} }{ contains the transcript count in each grid. #' Each row refers to a grid, and each column refers to a gene.} .compute_observation<- function(x, cluster_info, correlation_method, n_cores, - test_genes,bin_type, bin_param, w_x, w_y,use_cm){ -# primary_class <- class(x)[1] -# if (primary_class == "SpatialExperiment" | -# primary_class == "SpatialFeatureExperiment"){ -# if (is.null(x$sample_id) == FALSE){ -# if (length(unique(x$sample_id))>1){ -# stop("Input x has multiple samples") -# } -# cluster_info$sample <- unique(x$sample_id) -# }else{ -# cluster_info$sample <- "sample1" -# x$sample_id <- "sample1" -# } -# }else if (primary_class == "SingleCellExperiment"){ -# if (length(names(x@assays@data@listData)) == 1){ -# cluster_info$sample <- names(x@assays@data@listData) -# }else if (length(names(x@assays@data@listData)) > 1){ -# stop("Input x has multiple samples") -# }else { -# cluster_info$sample <- "sample1" -# names(x@assays@data@listData) <- "sample1" -# } -# }else{ -# stop("The input class of 'x' is not supported. -# Please convert 'x' to one of the following supported types: -# SingleCellExperiment, SpatialExperiment, or SpatialFeatureExperiment.") -# } + test_genes,bin_type, bin_param, use_cm){ vectors_lst <- get_vectors(x=x,sample_names=unique(cluster_info$sample), cluster_info = cluster_info, bin_type=bin_type, bin_param=bin_param, test_genes=test_genes, n_cores=n_cores, - w_x=w_x, w_y=w_y,use_cm =use_cm) + use_cm =use_cm) # calculate correlation of permuted clusters and gene obs.stat <- cor(x=vectors_lst$gene_mt, y=vectors_lst$cluster_mt, @@ -106,10 +76,6 @@ #' @param n_cores A positive number specifying number of cores used for #' parallelizing permutation testing. Default is one core #' (sequential processing). -#' @param w_x a numeric vector of length two specifying the x coordinate -#' limits of enclosing box. -#' @param w_y a numeric vector of length two specifying the y coordinate -#' limits of enclosing box. #' @param gene_mt A matrix contains the transcript count in each grid. #' Each row refers to a grid, and each column refers to a gene. #' @param cluster_names A list of strings giving the name and order of the @@ -119,13 +85,24 @@ #' @importFrom stats cor #' @return A matrix with permutation statistics #' + .compute_permutation<- function(cluster_info, perm.size = 1000, correlation_method = "pearson", bin_type, - bin_param, n_cores=1, w_x,w_y, gene_mt, + bin_param, n_cores=1, gene_mt, cluster_names){ n_clusters <- length(unique(cluster_info$cluster)) t.perm.array<- array(0, dim = c(ncol(gene_mt), length(cluster_names),perm.size)) + # Compute padded ranges for this sample + x_rng <- range(cluster_info$x, na.rm = TRUE) + y_rng <- range(cluster_info$y, na.rm = TRUE) + + pad_x <- BUFFER_FRAC * diff(x_rng) + pad_y <- BUFFER_FRAC * diff(y_rng) + + w_x <- c(x_rng[1] - pad_x, x_rng[2] + pad_x) + w_y <- c(y_rng[1] - pad_y, y_rng[2] + pad_y) + if (bin_type == "hexagon"){ w <-owin(xrange=w_x, yrange=w_y) H <-hextess(W=w, bin_param[1]) @@ -209,11 +186,6 @@ #' @param n_cores A positive number specifying number of cores used for #' parallelizing permutation testing. Default is one core #' (sequential processing). -#' @param w_x a numeric vector of length two specifying the x coordinate -#' limits of enclosing box. -#' @param w_y a numeric vector of length two specifying the y coordinate -#' limits of enclosing box. -#' #' @param correction_method A character string pass to \code{\link{p.adjust}} #' specifying the correction method for multiple testing . #' @param use_cm A boolean value that specifies whether to create spatial @@ -267,14 +239,6 @@ #' row = trans_info$feature_name, col = trans_info$cell ) #' spe_sample1 <- SpatialExperiment( #' assays = list(molecules = mol),sample_id ="sample1" ) -#' w_x <- c(min(floor(min(trans_info$x)), -#' floor(min(clusters$x))), -#' max(ceiling(max(trans_info$x)), -#' ceiling(max(clusters$x)))) -#' w_y <- c(min(floor(min(trans_info$y)), -#' floor(min(clusters$y))), -#' max(ceiling(max(trans_info$y)), -#' ceiling(max(clusters$y)))) #' set.seed(100) #' corr_res <- compute_permp(x=spe_sample1, #' cluster_info=clusters, @@ -284,9 +248,7 @@ #' test_genes=unique(trans_info$feature_name), #' correlation_method = "pearson", #' n_cores=1, -#' correction_method="BH", -#' w_x=w_x , -#' w_y=w_y) +#' correction_method="BH") #' #' # raw permutation p-value #' perm_p <- get_perm_p(corr_res) @@ -298,7 +260,7 @@ compute_permp<-function(x, cluster_info, perm.size, bin_type, bin_param,test_genes, correlation_method = "pearson", n_cores=1, - correction_method="BH",w_x ,w_y,use_cm = FALSE){ + correction_method="BH",use_cm = FALSE){ message(sprintf("Correlation Method = %s", correlation_method)) tm1 <- system.time( @@ -307,7 +269,7 @@ compute_permp<-function(x, cluster_info, perm.size, bin_type, n_cores=n_cores,use_cm =use_cm, correlation_method = correlation_method, bin_type=bin_type,test_genes=test_genes, - bin_param=bin_param, w_x=w_x, w_y = w_y) + bin_param=bin_param) }) obs.stat<- obs_res$obs.stat @@ -325,7 +287,6 @@ compute_permp<-function(x, cluster_info, perm.size, bin_type, correlation_method = correlation_method, bin_type=bin_type, bin_param=bin_param, n_cores=n_cores, - w_x=w_x, w_y = w_y, gene_mt = obs_res$gene_mt, cluster_names =colnames(obs.stat) ) # permutation stats diff --git a/R/cor_mg_result.R b/R/cor_mg_result.R index ed7c356..3b00da4 100644 --- a/R/cor_mg_result.R +++ b/R/cor_mg_result.R @@ -59,14 +59,6 @@ #' row = trans_info$feature_name, col = trans_info$cell ) #' spe_sample1 <- SpatialExperiment( #' assays = list(molecules = mol),sample_id ="sample1" ) -#' w_x <- c(min(floor(min(trans_info$x)), -#' floor(min(clusters$x))), -#' max(ceiling(max(trans_info$x)), -#' ceiling(max(clusters$x)))) -#' w_y <- c(min(floor(min(trans_info$y)), -#' floor(min(clusters$y))), -#' max(ceiling(max(trans_info$y)), -#' ceiling(max(clusters$y)))) #' set.seed(100) #' corr_res <- compute_permp(x=spe_sample1, #' cluster_info=clusters, @@ -76,9 +68,7 @@ #' test_genes=unique(trans_info$feature_name), #' correlation_method = "pearson", #' n_cores=1, -#' correction_method="BH", -#' w_x=w_x , -#' w_y=w_y) +#' correction_method="BH") #' # observed correlation #' obs_corr <- get_cor(corr_res) get_cor <- function(obj) { @@ -127,14 +117,7 @@ get_cor <- function(obj) { #' row = trans_info$feature_name, col = trans_info$cell ) #' spe_sample1 <- SpatialExperiment( #' assays = list(molecules = mol),sample_id ="sample1" ) -#' w_x <- c(min(floor(min(trans_info$x)), -#' floor(min(clusters$x))), -#' max(ceiling(max(trans_info$x)), -#' ceiling(max(clusters$x)))) -#' w_y <- c(min(floor(min(trans_info$y)), -#' floor(min(clusters$y))), -#' max(ceiling(max(trans_info$y)), -#' ceiling(max(clusters$y)))) + #' set.seed(100) #' corr_res <- compute_permp(x=spe_sample1, #' cluster_info=clusters, @@ -144,9 +127,7 @@ get_cor <- function(obj) { #' test_genes=unique(trans_info$feature_name), #' correlation_method = "pearson", #' n_cores=1, -#' correction_method="BH", -#' w_x=w_x , -#' w_y=w_y) +#' correction_method="BH") #' #' # raw permutation p-value #' perm_p <- get_perm_p(corr_res) @@ -195,14 +176,7 @@ get_perm_p <- function(obj) { #' row = trans_info$feature_name, col = trans_info$cell ) #' spe_sample1 <- SpatialExperiment( #' assays = list(molecules = mol),sample_id ="sample1" ) -#' w_x <- c(min(floor(min(trans_info$x)), -#' floor(min(clusters$x))), -#' max(ceiling(max(trans_info$x)), -#' ceiling(max(clusters$x)))) -#' w_y <- c(min(floor(min(trans_info$y)), -#' floor(min(clusters$y))), -#' max(ceiling(max(trans_info$y)), -#' ceiling(max(clusters$y)))) + #' set.seed(100) #' corr_res <- compute_permp(x=spe_sample1, #' cluster_info=clusters, @@ -212,9 +186,7 @@ get_perm_p <- function(obj) { #' test_genes=unique(trans_info$feature_name), #' correlation_method = "pearson", #' n_cores=1, -#' correction_method="BH", -#' w_x=w_x , -#' w_y=w_y) +#' correction_method="BH") #' # adjusted permutation p-value #' adjusted_perm_p <- get_perm_adjp(corr_res) get_perm_adjp <- function(obj) { diff --git a/R/create_genesets.R b/R/create_genesets.R index 6e80e0f..b74ba73 100644 --- a/R/create_genesets.R +++ b/R/create_genesets.R @@ -26,10 +26,6 @@ #' two giving the numbers of rectangular quadrats in the x and y directions. If #' the \code{bin_type} is "hexagonal", this will be a number giving the side #' length of hexagons. Positive numbers only. -#' @param w_x A numeric vector of length two specifying the x coordinate -#' limits of enclosing box. -#' @param w_y A numeric vector of length two specifying the y coordinate -#' limits of enclosing box. #' @param use_cm A boolean value that specifies whether to create spatial #' vectors for genes using the count matrix and cell coordinates instead of #' the transcript coordinates when both types of information are available. @@ -52,10 +48,6 @@ #' two giving the numbers of rectangular quadrats in the x and y directions. If #' the \code{bin_type} is "hexagonal", this will be a number giving the side #' length of hexagons. Positive numbers only. -#' @param w_x A numeric vector of length two specifying the x coordinate -#' limits of enclosing box. -#' @param w_y A numeric vector of length two specifying the y coordinate -#' limits of enclosing box. #' @param cluster_info A dataframe/matrix containing the centroid coordinates, #' cluster and sample label for each cell.The column names must include #' "x" (x coordinate), "y" (y coordinate), @@ -94,25 +86,24 @@ #' name_lst=list(dummy_A=c("A","C"), #' dummy_B=c("A","B","C")), #' bin_type="square", -#' bin_param=c(2,2), -#' w_x=c(0,25), w_y=c(0,25), cluster_info=NULL) +#' bin_param=c(2,2),cluster_info=NULL) #' create_genesets<-function(x, name_lst, cluster_info, sample_names, bin_type, bin_param, - w_x, w_y, use_cm = FALSE, n_cores=1){ + use_cm = FALSE, n_cores=1){ primary_class <- class(x)[1] if (primary_class == "list"){ + uni_genes <- unique(unlist(lapply(x, function(df) df$feature_name))) for (sp in sample_names){ sub_x <- x[[sp]] req_cols <- c("feature_name","x","y") if (length(setdiff(req_cols, colnames(sub_x))) != 0){ stop("Invalid columns in input x. Every list element in x must contain columns 'feature_name','x', 'y', - for every transcipt")} - uni_genes <- unique(sub_x$feature_name) + for every transcipt")} } if (length(setdiff(as.vector(unique(unlist(name_lst))),uni_genes))>0){ stop("Invalid name_lst, - can not find genes in the input name_lst from x")}} + can not find genes in the input name_lst from x")} test_genes <- intersect(as.vector(unique(unlist(name_lst))), uni_genes) }else { if (length(setdiff(as.vector(unique(unlist(name_lst))), row.names(x)))>0){ @@ -124,7 +115,7 @@ create_genesets<-function(x, name_lst, res_lst <- get_vectors(x=x, cluster_info=cluster_info, sample_names=sample_names, bin_type=bin_type, - bin_param=bin_param,test_genes=test_genes, w_x=w_x, w_y=w_y, + bin_param=bin_param,test_genes=test_genes, use_cm = use_cm, n_cores=n_cores) gene_vecs <- as.matrix(res_lst$gene_mt) diff --git a/R/get_vectors.R b/R/get_vectors.R index d2cca47..97e51dd 100644 --- a/R/get_vectors.R +++ b/R/get_vectors.R @@ -17,42 +17,57 @@ #' the \code{bin_type} is "hexagonal", this will be a number giving the side #' length of hexagons. Positive numbers only. #' @param bin_length A positive integer giving the length of total bins -#' @param w_x A numeric vector of length two specifying the x coordinate -#' limits of enclosing box. -#' @param w_y A numeric vector of length two specifying the y coordinate -#' limits of enclosing box. +#' @param range_list A named list of spatial ranges for each sample. +#' Each element should be a list with two components: +#' \code{w_x} and \code{w_y}, which are numeric vectors of length 2 +#' specifying the x- and y-axis ranges (e.g., from cell or +#' transcript coordinates). +#' The range is calculated with 5% buffer to ensure all points fall +#' within the window. #' @importFrom BiocParallel bplapply #' @return a matrix contains the transcript count in each grid. #' Each row refers to a grid, and each column refers to a gene. .get_gene_vectors_tr<- function(trans_lst, test_genes, bin_type, bin_param, - bin_length, w_x, w_y){ + bin_length, range_list){ n_genes <- length(test_genes) n_samples <- length(trans_lst) vec_gene_mt <- as.data.frame(matrix(0, ncol=n_genes, nrow=bin_length*n_samples)) colnames(vec_gene_mt) <- test_genes - if (bin_type == "hexagon"){ - w <- owin(xrange=w_x, yrange=w_y) - H <- hextess(W=w, bin_param[1]) - } - calculate_one_gene <- function(i_gene){ + + calculate_one_gene <- function(i_gene) { vec_gene <- c() - for (rpp in trans_lst){ - curr <- rpp[rpp$feature_name==i_gene, - c("x","y")] %>% distinct() - gene_ppp <- ppp(curr$x,curr$y,w_x, w_y) - # create gene vector - if (bin_type == "hexagon"){ - vec_g <- as.vector(t(quadratcount(gene_ppp, tess=H))) - }else{ - vec_g <- as.vector(t(quadratcount(gene_ppp, - bin_param[1],bin_param[2]))) + for (sample_name in names(trans_lst)) { + rpp <- trans_lst[[sample_name]] + ranges <- range_list[[sample_name]] + if (!(i_gene %in% rpp$feature_name)) { + # Append zeros if gene not found + vec_g <- rep(0, bin_length) + } else { + curr <- rpp[rpp$feature_name == i_gene, c("x", "y")] + curr <- unique(curr) + # Create point pattern using sample-specific window + gene_ppp <- ppp(curr$x, curr$y, + window = owin(xrange = ranges$w_x, yrange = ranges$w_y) + ) + # Create binning and count + if (bin_type == "hexagon") { + H <- hextess(W = gene_ppp$window, + bin_param[1]) + vec_g <- as.vector(t(quadratcount(gene_ppp, tess = H))) + } else { + vec_g <- as.vector(t(quadratcount(gene_ppp, + nx = bin_param[1], ny = bin_param[2] + ))) + } } vec_gene <- c(vec_gene, vec_g) } + return(vec_gene) } - + + # to calculate all genes in parallel result_lst <- bplapply(test_genes, calculate_one_gene) @@ -90,10 +105,13 @@ #' @param test_genes A vector of strings giving the name of the genes you #' want to test. This will be used as column names for one of the result matrix #' \code{gene_mt}. -#' @param w_x A numeric vector of length two specifying the x coordinate -#' limits of enclosing box. -#' @param w_y A numeric vector of length two specifying the y coordinate -#' limits of enclosing box. +#' @param range_list A named list of spatial ranges for each sample. +#' Each element should be a list with two components: +#' \code{w_x} and \code{w_y}, which are numeric vectors of length 2 +#' specifying the x- and y-axis ranges (e.g., from cell or +#' transcript coordinates). +#' The range is calculated with 5% buffer to ensure all points fall +#' within the window. #' #' @importFrom BiocParallel bplapply #' @@ -112,43 +130,70 @@ #' @importFrom magrittr "%>%" #' .get_gene_vectors_cm<- function(cluster_info, cm_lst, bin_type, bin_param, - test_genes, w_x, w_y){ - # binning - bin_length <- 0 - if (bin_type == "hexagon"){ - w <- owin(xrange=w_x, yrange=w_y) - H <- hextess(W=w, bin_param[1]) - bin_length <- length(H$tiles) - }else{ bin_length <- bin_param[1] * bin_param[2] } + test_genes,range_list){ n_genes <- length(test_genes) - index_vec <- NULL - # use count matrix to build gene vector matrix - if (bin_type == "hexagon"){ - tile_indices <- tileindex(x=cluster_info$x, y=cluster_info$y, Z=H) - cluster_info$index_vec <- match(tile_indices,levels(tile_indices)) - }else{ - all_cell_intervals <- as.tess(quadratcount(ppp(cluster_info$x, - cluster_info$y, w_x, w_y),bin_param[1],bin_param[2])) - cluster_info$index_x <- findInterval(cluster_info$x, - all_cell_intervals$xgrid, all.inside = TRUE, - left.open=FALSE, rightmost.closed=TRUE ) - cluster_info$index_y <- findInterval(cluster_info$y, - all_cell_intervals$ygrid, all.inside = TRUE, - left.open=FALSE, rightmost.closed=TRUE ) - cluster_info$index_y <- bin_param[2] + 1 - cluster_info$index_y - # column major - ind_vec <-(cluster_info$index_y-1)* bin_param[2]+cluster_info$index_x - cluster_info$index_vec <- ind_vec } - n_samples <- length(cm_lst) + n_samples <- length(unique(cluster_info$sample)) + + # Use first sample as template for binning layout + example_sample <- names(range_list)[1] + ranges <- range_list[[example_sample]] + + if (bin_type == "hexagon") { + shared_H <- hextess( + W = owin(xrange = ranges$w_x, yrange = ranges$w_y), + bin_param[1] + ) + bin_length <- length(shared_H$tiles) + } else { + bin_length <- bin_param[1] * bin_param[2] + } + + cluster_info$index_vec <- NA + for (rp_nm in unique(cluster_info$sample)) { + i_info <- cluster_info[cluster_info$sample == rp_nm, ] + ranges <- range_list[[rp_nm]] + + if (bin_type == "hexagon") { + tile_indices <-tileindex(x = i_info$x, y = i_info$y, Z = shared_H) + cluster_info$index_vec[cluster_info$sample == rp_nm] <- + match(tile_indices, levels(tile_indices)) + } else { + tess <- as.tess(quadratcount(ppp(i_info$x, i_info$y, + xrange = ranges$w_x, yrange = ranges$w_y), + bin_param[1], bin_param[2] + ) + ) + + ix <- findInterval(i_info$x, tess$xgrid, + all.inside = TRUE, + left.open = FALSE, + rightmost.closed = TRUE) + + iy <- findInterval(i_info$y, tess$ygrid, + all.inside = TRUE, + left.open = FALSE, + rightmost.closed = TRUE) + # column major + iy <- bin_param[2] + 1 - iy # flip y for column-major order + ind_vec <- (iy - 1) * bin_param[2] + ix + cluster_info$index_vec[cluster_info$sample == rp_nm] <- ind_vec + } + } # create gene vector vec_gene_mt<-as.data.frame(matrix(0,ncol=n_genes,nrow=bin_length*n_samples)) colnames(vec_gene_mt) <- test_genes count_value <- NULL - calculate_one_gene <- function(i_gene){ + calculate_one_gene <- function(i_gene) { vec_gene <- c() - for (rp_nm in names(cm_lst)){ + for (rp_nm in names(cm_lst)) { cm <- cm_lst[[rp_nm]] vec_g <- rep(0, bin_length) + + # If gene is not present in this sample, use zero vector + if (!(i_gene %in% rownames(cm))) { + vec_gene <- c(vec_gene, vec_g) + next + } i_gene_mt <- cluster_info[cluster_info$sample==rp_nm, ] i_gene_mt$count_value <- as.numeric(cm[i_gene,i_gene_mt$cell_id]) i_gene_mt$index_vec <- factor(i_gene_mt$index_vec ) @@ -184,16 +229,19 @@ #' two giving the numbers of rectangular quadrats in the x and y directions. If #' the \code{bin_type} is "hexagonal", this will be a number giving the side #' length of hexagons. Positive numbers only. -#' @param w_x A numeric vector of length two specifying the x coordinate -#' limits of enclosing box. -#' @param w_y A numeric vector of length two specifying the y coordinate -#' limits of enclosing box. +#' @param range_list A named list of spatial ranges for each sample. +#' Each element should be a list with two components: +#' \code{w_x} and \code{w_y}, which are numeric vectors of length 2 +#' specifying the x- and y-axis ranges (e.g., from cell or +#' transcript coordinates). +#' The range is calculated with 5% buffer to ensure all points fall +#' within the window. #' @param sample_names a vector of strings giving the sample names #' @return a matrix contains the cell count in each grid. #' Each row refers to a grid, and each column refers to a cluster. #' .get_cluster_vectors<- function(cluster_info,bin_length,bin_type, bin_param, - w_x, w_y, sample_names){ + range_list, sample_names){ ava_sample_names <- unique(as.character(cluster_info$sample)) if (length(setdiff(sample_names,ava_sample_names))>0){ stop("Can not find input sample_names from input cluster_info") @@ -221,6 +269,8 @@ for (nm in sample_cluster_nm ){ rp <- unlist(strsplit(nm, split="--"))[1] cl <- unlist(strsplit(nm, split="--"))[2] + w_x <- range_list[[rp]]$w_x + w_y <- range_list[[rp]]$w_y cluster_ppp <- ppp(cluster_info[cluster_info$sample==rp & cluster_info$cluster==cl, "x"], @@ -253,14 +303,17 @@ #' two giving the numbers of rectangular quadrats in the x and y directions. If #' the \code{bin_type} is "hexagonal", this will be a number giving the side #' length of hexagons. Positive numbers only. -#' @param w_x A numeric vector of length two specifying the x coordinate -#' limits of enclosing box. -#' @param w_y A numeric vector of length two specifying the y coordinate -#' limits of enclosing box. +#' @param range_list A named list of spatial ranges for each sample. +#' Each element should be a list with two components: +#' \code{w_x} and \code{w_y}, which are numeric vectors of length 2 +#' specifying the x- and y-axis ranges (e.g., from cell or +#' transcript coordinates). +#' The range is calculated with 5% buffer to ensure all points fall +#' within the window. #' #' @return the length of total bins -.check_binning<- function(bin_param, bin_type, w_x, w_y){ +.check_binning<- function(bin_param, bin_type, range_list){ # binning bin_length <- 0 if (bin_type == "hexagon"){ @@ -268,7 +321,7 @@ stop("Invalid input bin_param, bin_param should be a vector of length 2 for hexagon bins") } - w <- owin(xrange=w_x, yrange=w_y) + w <- owin(xrange=range_list[[1]]$w_x, yrange=range_list[[1]]$w_y) H <- hextess(W=w, bin_param[1]) bin_length <- length(H$tiles) }else if (bin_type == "square" | bin_type == "rectangle"){ @@ -394,6 +447,7 @@ SingleCellExperiment, SpatialExperiment, or SpatialFeatureExperiment.")) } return (list(trans_lst = trans_lst, cm_lst = cm_lst)) } + #' Vectorise the spatial coordinates #' #' @description @@ -454,10 +508,6 @@ SingleCellExperiment, SpatialExperiment, or SpatialFeatureExperiment.")) } #' @param test_genes A vector of strings giving the name of the genes you want #' to create gene vector. This will be used as column names for one of the #' result matrix \code{gene_mt}. -#' @param w_x A numeric vector of length two specifying the x coordinate -#' limits of enclosing box. -#' @param w_y A numeric vector of length two specifying the y coordinate -#' limits of enclosing box. #' @param use_cm A boolean value that specifies whether to create spatial #' vectors for genes using the count matrix and cell coordinates instead of #' the transcript coordinates when both types of information are available. @@ -511,22 +561,13 @@ SingleCellExperiment, SpatialExperiment, or SpatialFeatureExperiment.")) } #' trans_info$y=as.numeric(trans_info$y) #' trans_info$cell = sample(c("cell1","cell2","cell2"),replace=TRUE, #' size=nrow(trans_info)) -#' w_x = c(min(floor(min(trans_info$x)), -#' floor(min(clusters$x))), -#' max(ceiling(max(trans_info$x)), -#' ceiling(max(clusters$x)))) -#' w_y = c(min(floor(min(trans_info$y)), -#' floor(min(clusters$y))), -#' max(ceiling(max(trans_info$y)), -#' ceiling(max(clusters$y)))) #' # use named list as input #' vecs_lst = get_vectors(x= list("sample1" = trans_info), #' sample_names=c("sample1"), #' cluster_info = clusters, #' bin_type = "square", #' bin_param = c(5,5), -#' test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2"), -#' w_x = w_x, w_y=w_y) +#' test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2")) #' # use SpatialExperiment object as input #' trans_mol <- BumpyMatrix::splitAsBumpyMatrix( #' trans_info[, c("x", "y")], @@ -537,11 +578,10 @@ SingleCellExperiment, SpatialExperiment, or SpatialFeatureExperiment.")) } #' cluster_info = clusters, #' bin_type = "square", #' bin_param = c(5,5), -#' test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2"), -#' w_x = w_x, w_y=w_y) +#' test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2")) #' get_vectors<- function(x, cluster_info, sample_names, bin_type, bin_param, - test_genes, w_x, w_y, use_cm = FALSE, n_cores=1){ + test_genes, use_cm = FALSE, n_cores=1){ # check input if ((is.null(x) ==TRUE) & (is.null(cluster_info) == TRUE) ){ stop("Invalid input, no coordinates information is specified") } @@ -579,8 +619,43 @@ get_vectors<- function(x, cluster_info, sample_names, bin_type, bin_param, # must provide cluster info if gene vectors are created from count matrix if ((is.null(cm_lst) == FALSE) & (is.null(cluster_info) == TRUE)){ stop("Missing cluster information to build gene vector matrix.")} + + if (is.null(trans_lst) == FALSE ){ + # check the range for transcript coordinates only + # For each sample, compute x and y ranges from trasncript coordinates + range_list <- lapply(trans_lst, function(d) { + x_rng <- range(d$x, na.rm = TRUE) + y_rng <- range(d$y, na.rm = TRUE) + + # 1% padding on each side + pad_x <- BUFFER_FRAC * diff(x_rng) + pad_y <- BUFFER_FRAC * diff(y_rng) + + list( + w_x = c(x_rng[1] - pad_x, x_rng[2] + pad_x), + w_y = c(y_rng[1] - pad_y, y_rng[2] + pad_y) + ) + }) + }else{ + # check the range for cell coordinates only + # For each sample, compute x and y ranges from cell coordinates + range_list <- lapply(split(cluster_info, cluster_info$sample), + function(d) { + x_rng <- range(d$x, na.rm = TRUE) + y_rng <- range(d$y, na.rm = TRUE) + + # 1% padding on each side + pad_x <- BUFFER_FRAC * diff(x_rng) + pad_y <- BUFFER_FRAC * diff(y_rng) + list(w_x = c(x_rng[1] - pad_x, x_rng[2] + pad_x), + w_y = c(y_rng[1] - pad_y, y_rng[2] + pad_y) + ) + }) + + } + bin_length<-.check_binning(bin_param=bin_param,bin_type=bin_type, - w_x=w_x,w_y=w_y) + range_list=range_list) # register for BiocParallel register(SnowParam(workers = n_cores, type = "SOCK")) # with cluster information @@ -589,22 +664,23 @@ get_vectors<- function(x, cluster_info, sample_names, bin_type, bin_param, cluster_info$cluster<- .check_valid_names(cluster_info$cluster,"cluster") cluster_info$sample<- .check_valid_names(cluster_info$sample,"sample") sample_names<- .check_valid_names(sample_names,"sample") + names(range_list) <- make.names(names(range_list), unique=TRUE) vec_cluster <- .get_cluster_vectors(cluster_info=cluster_info, bin_length=bin_length, bin_type=bin_type, - bin_param=bin_param,w_x=w_x,w_y=w_y, + bin_param=bin_param,range_list=range_list, sample_names=sample_names) } # with gene information if (is.null(trans_lst) == FALSE & is.null(cm_lst)== TRUE){ vec_gene_mt<-.get_gene_vectors_tr(trans_lst=trans_lst, test_genes=test_genes,bin_type=bin_type, bin_param=bin_param, - bin_length=bin_length,w_x=w_x,w_y=w_y)} + bin_length=bin_length,range_list=range_list)} if (is.null(trans_lst) == TRUE & is.null(cm_lst)==FALSE & is.null(cluster_info)==FALSE){ # use count matrix to build gene vector matrix vec_gene_mt<-.get_gene_vectors_cm(cluster_info=cluster_info, cm_lst=cm_lst, bin_type=bin_type, bin_param=bin_param, - test_genes=test_genes,w_x=w_x, w_y=w_y) } + test_genes=test_genes,range_list=range_list) } result <- list() if ((is.null(cluster_info) == FALSE) ){ result$cluster_mt <- as.matrix(vec_cluster)} diff --git a/R/globals.R b/R/globals.R new file mode 100644 index 0000000..84e73cc --- /dev/null +++ b/R/globals.R @@ -0,0 +1,6 @@ +#' Default padding parameter for function get_vectors() +#' +#' A numeric value giving the padding offset proportion when defining +#' the range for x and y coordinates +#' @keywords internal +BUFFER_FRAC <- 0.05 \ No newline at end of file diff --git a/R/lasso_markers.R b/R/lasso_markers.R index f32d7df..e7f62ee 100644 --- a/R/lasso_markers.R +++ b/R/lasso_markers.R @@ -264,20 +264,11 @@ #' row = trans_info$feature_name, col = trans_info$cell ) #' spe<- SpatialExperiment( #' assays = list(molecules = trans_mol),sample_id ="sample1" ) -#' w_x <- c(min(floor(min(trans_info$x)), -#' floor(min(clusters$x))), -#' max(ceiling(max(trans_info$x)), -#' ceiling(max(clusters$x)))) -#' w_y <- c(min(floor(min(trans_info$y)), -#' floor(min(clusters$y))), -#' max(ceiling(max(trans_info$y)), -#' ceiling(max(clusters$y)))) #' vecs_lst <- get_vectors(x=spe,sample_names=c("sample1"), #' cluster_info = clusters, #' bin_type = "square", #' bin_param = c(20,20), -#' test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2"), -#' w_x = w_x, w_y=w_y) +#' test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2")) #' lasso_res <- lasso_markers(gene_mt=vecs_lst$gene_mt, #' cluster_mt = vecs_lst$cluster_mt, #' sample_names=c("sample1"), diff --git a/R/lm_mg_result.R b/R/lm_mg_result.R index 31f5bb8..dd92683 100644 --- a/R/lm_mg_result.R +++ b/R/lm_mg_result.R @@ -89,20 +89,11 @@ #' row = trans_info$feature_name, col = trans_info$cell ) #' spe<- SpatialExperiment( #' assays = list(molecules = trans_mol),sample_id ="sample1" ) -#' w_x <- c(min(floor(min(trans_info$x)), -#' floor(min(clusters$x))), -#' max(ceiling(max(trans_info$x)), -#' ceiling(max(clusters$x)))) -#' w_y <- c(min(floor(min(trans_info$y)), -#' floor(min(clusters$y))), -#' max(ceiling(max(trans_info$y)), -#' ceiling(max(clusters$y)))) #' vecs_lst <- get_vectors(x=spe,sample_names=c("sample1"), #' cluster_info = clusters, #' bin_type = "square", #' bin_param = c(20,20), -#' test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2"), -#' w_x = w_x, w_y=w_y) +#' test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2")) #' lasso_res <- lasso_markers(gene_mt=vecs_lst$gene_mt, #' cluster_mt = vecs_lst$cluster_mt, #' sample_names=c("sample1"), @@ -178,20 +169,11 @@ get_top_mg <- function(obj, coef_cutoff=0.05) { #' row = trans_info$feature_name, col = trans_info$cell ) #' spe<- SpatialExperiment( #' assays = list(molecules = trans_mol),sample_id ="sample1" ) -#' w_x <- c(min(floor(min(trans_info$x)), -#' floor(min(clusters$x))), -#' max(ceiling(max(trans_info$x)), -#' ceiling(max(clusters$x)))) -#' w_y <- c(min(floor(min(trans_info$y)), -#' floor(min(clusters$y))), -#' max(ceiling(max(trans_info$y)), -#' ceiling(max(clusters$y)))) #' vecs_lst <- get_vectors(x=spe,sample_names=c("sample1"), #' cluster_info = clusters, #' bin_type = "square", #' bin_param = c(20,20), -#' test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2"), -#' w_x = w_x, w_y=w_y) +#' test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2")) #' lasso_res <- lasso_markers(gene_mt=vecs_lst$gene_mt, #' cluster_mt = vecs_lst$cluster_mt, #' sample_names=c("sample1"), diff --git a/jazzPanda.Rproj b/jazzPanda.Rproj new file mode 100644 index 0000000..cf0a065 --- /dev/null +++ b/jazzPanda.Rproj @@ -0,0 +1,18 @@ +Version: 1.0 +ProjectId: b78cc7e6-acce-4286-b91a-f34d2f9c40a7 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 4 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX + +BuildType: Package +PackageUseDevtools: Yes +PackageInstallArgs: --no-multiarch --with-keep.source diff --git a/man/BUFFER_FRAC.Rd b/man/BUFFER_FRAC.Rd new file mode 100644 index 0000000..7936da3 --- /dev/null +++ b/man/BUFFER_FRAC.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/globals.R +\docType{data} +\name{BUFFER_FRAC} +\alias{BUFFER_FRAC} +\title{Default padding parameter for function get_vectors()} +\format{ +An object of class \code{numeric} of length 1. +} +\usage{ +BUFFER_FRAC +} +\description{ +A numeric value giving the padding offset proportion when defining +the range for x and y coordinates +} +\keyword{internal} diff --git a/man/compute_permp.Rd b/man/compute_permp.Rd index dc26280..e390ca5 100644 --- a/man/compute_permp.Rd +++ b/man/compute_permp.Rd @@ -14,8 +14,6 @@ compute_permp( correlation_method = "pearson", n_cores = 1, correction_method = "BH", - w_x, - w_y, use_cm = FALSE ) } @@ -58,12 +56,6 @@ parallelizing permutation testing. Default is one core \item{correction_method}{A character string pass to \code{\link{p.adjust}} specifying the correction method for multiple testing .} -\item{w_x}{a numeric vector of length two specifying the x coordinate -limits of enclosing box.} - -\item{w_y}{a numeric vector of length two specifying the y coordinate -limits of enclosing box.} - \item{use_cm}{A boolean value that specifies whether to create spatial vectors for genes using the count matrix and cell coordinates instead of the transcript coordinates when both types of information are available. @@ -126,14 +118,6 @@ mol <- BumpyMatrix::splitAsBumpyMatrix( row = trans_info$feature_name, col = trans_info$cell ) spe_sample1 <- SpatialExperiment( assays = list(molecules = mol),sample_id ="sample1" ) -w_x <- c(min(floor(min(trans_info$x)), - floor(min(clusters$x))), - max(ceiling(max(trans_info$x)), - ceiling(max(clusters$x)))) -w_y <- c(min(floor(min(trans_info$y)), - floor(min(clusters$y))), - max(ceiling(max(trans_info$y)), - ceiling(max(clusters$y)))) set.seed(100) corr_res <- compute_permp(x=spe_sample1, cluster_info=clusters, @@ -143,9 +127,7 @@ corr_res <- compute_permp(x=spe_sample1, test_genes=unique(trans_info$feature_name), correlation_method = "pearson", n_cores=1, - correction_method="BH", - w_x=w_x , - w_y=w_y) + correction_method="BH") # raw permutation p-value perm_p <- get_perm_p(corr_res) diff --git a/man/create_genesets.Rd b/man/create_genesets.Rd index 7645f2a..987672a 100644 --- a/man/create_genesets.Rd +++ b/man/create_genesets.Rd @@ -11,8 +11,6 @@ create_genesets( sample_names, bin_type, bin_param, - w_x, - w_y, use_cm = FALSE, n_cores = 1 ) @@ -45,12 +43,6 @@ two giving the numbers of rectangular quadrats in the x and y directions. If the \code{bin_type} is "hexagonal", this will be a number giving the side length of hexagons. Positive numbers only.} -\item{w_x}{A numeric vector of length two specifying the x coordinate -limits of enclosing box.} - -\item{w_y}{A numeric vector of length two specifying the y coordinate -limits of enclosing box.} - \item{use_cm}{A boolean value that specifies whether to create spatial vectors for genes using the count matrix and cell coordinates instead of the transcript coordinates when both types of information are available. @@ -102,7 +94,6 @@ geneset_res <- create_genesets(x=rep1_spe, sample=c("sample1"), name_lst=list(dummy_A=c("A","C"), dummy_B=c("A","B","C")), bin_type="square", - bin_param=c(2,2), - w_x=c(0,25), w_y=c(0,25), cluster_info=NULL) + bin_param=c(2,2),cluster_info=NULL) } diff --git a/man/dot-check_binning.Rd b/man/dot-check_binning.Rd index dd63a0c..b4e4350 100644 --- a/man/dot-check_binning.Rd +++ b/man/dot-check_binning.Rd @@ -4,7 +4,7 @@ \alias{.check_binning} \title{helper function to check the input of binning} \usage{ -.check_binning(bin_param, bin_type, w_x, w_y) +.check_binning(bin_param, bin_type, range_list) } \arguments{ \item{bin_param}{A numeric vector indicating the size of the bin. If the @@ -16,11 +16,13 @@ length of hexagons. Positive numbers only.} \item{bin_type}{A string indicating which bin shape is to be used for vectorization. One of "square" (default), "rectangle", or "hexagon".} -\item{w_x}{A numeric vector of length two specifying the x coordinate -limits of enclosing box.} - -\item{w_y}{A numeric vector of length two specifying the y coordinate -limits of enclosing box.} +\item{range_list}{A named list of spatial ranges for each sample. +Each element should be a list with two components: +\code{w_x} and \code{w_y}, which are numeric vectors of length 2 +specifying the x- and y-axis ranges (e.g., from cell or +transcript coordinates). +The range is calculated with 5% buffer to ensure all points fall +within the window.} } \value{ the length of total bins diff --git a/man/dot-compute_observation.Rd b/man/dot-compute_observation.Rd index 603013c..7cc1965 100644 --- a/man/dot-compute_observation.Rd +++ b/man/dot-compute_observation.Rd @@ -12,8 +12,6 @@ test_genes, bin_type, bin_param, - w_x, - w_y, use_cm ) } @@ -55,12 +53,6 @@ two giving the numbers of rectangular quadrats in the x and y directions. If the \code{bin_type} is "hexagonal", this will be a number giving the side length of hexagons. Positive numbers only.} -\item{w_x}{a numeric vector of length two specifying the x coordinate -limits of enclosing box.} - -\item{w_y}{a numeric vector of length two specifying the y coordinate -limits of enclosing box.} - \item{use_cm}{A boolean value that specifies whether to create spatial vectors for genes using the count matrix and cell coordinates instead of the transcript coordinates when both types of information are available. diff --git a/man/dot-compute_permutation.Rd b/man/dot-compute_permutation.Rd index 466957e..5e1653f 100644 --- a/man/dot-compute_permutation.Rd +++ b/man/dot-compute_permutation.Rd @@ -11,8 +11,6 @@ bin_type, bin_param, n_cores = 1, - w_x, - w_y, gene_mt, cluster_names ) @@ -42,12 +40,6 @@ length of hexagons. Positive numbers only.} parallelizing permutation testing. Default is one core (sequential processing).} -\item{w_x}{a numeric vector of length two specifying the x coordinate -limits of enclosing box.} - -\item{w_y}{a numeric vector of length two specifying the y coordinate -limits of enclosing box.} - \item{gene_mt}{A matrix contains the transcript count in each grid. Each row refers to a grid, and each column refers to a gene.} diff --git a/man/dot-get_cluster_vectors.Rd b/man/dot-get_cluster_vectors.Rd index d560605..80e3852 100644 --- a/man/dot-get_cluster_vectors.Rd +++ b/man/dot-get_cluster_vectors.Rd @@ -9,8 +9,7 @@ bin_length, bin_type, bin_param, - w_x, - w_y, + range_list, sample_names ) } @@ -31,11 +30,13 @@ two giving the numbers of rectangular quadrats in the x and y directions. If the \code{bin_type} is "hexagonal", this will be a number giving the side length of hexagons. Positive numbers only.} -\item{w_x}{A numeric vector of length two specifying the x coordinate -limits of enclosing box.} - -\item{w_y}{A numeric vector of length two specifying the y coordinate -limits of enclosing box.} +\item{range_list}{A named list of spatial ranges for each sample. +Each element should be a list with two components: +\code{w_x} and \code{w_y}, which are numeric vectors of length 2 +specifying the x- and y-axis ranges (e.g., from cell or +transcript coordinates). +The range is calculated with 5% buffer to ensure all points fall +within the window.} \item{sample_names}{a vector of strings giving the sample names} } diff --git a/man/dot-get_gene_vectors_cm.Rd b/man/dot-get_gene_vectors_cm.Rd index 0d2ea0d..5b0d21e 100644 --- a/man/dot-get_gene_vectors_cm.Rd +++ b/man/dot-get_gene_vectors_cm.Rd @@ -10,8 +10,7 @@ bin_type, bin_param, test_genes, - w_x, - w_y + range_list ) } \arguments{ @@ -39,11 +38,13 @@ length of hexagons. Positive numbers only.} want to test. This will be used as column names for one of the result matrix \code{gene_mt}.} -\item{w_x}{A numeric vector of length two specifying the x coordinate -limits of enclosing box.} - -\item{w_y}{A numeric vector of length two specifying the y coordinate -limits of enclosing box.} +\item{range_list}{A named list of spatial ranges for each sample. +Each element should be a list with two components: +\code{w_x} and \code{w_y}, which are numeric vectors of length 2 +specifying the x- and y-axis ranges (e.g., from cell or +transcript coordinates). +The range is calculated with 5% buffer to ensure all points fall +within the window.} } \value{ a matrix contains the transcript count in each grid. diff --git a/man/dot-get_gene_vectors_tr.Rd b/man/dot-get_gene_vectors_tr.Rd index a2e5918..492b260 100644 --- a/man/dot-get_gene_vectors_tr.Rd +++ b/man/dot-get_gene_vectors_tr.Rd @@ -10,8 +10,7 @@ bin_type, bin_param, bin_length, - w_x, - w_y + range_list ) } \arguments{ @@ -34,11 +33,13 @@ length of hexagons. Positive numbers only.} \item{bin_length}{A positive integer giving the length of total bins} -\item{w_x}{A numeric vector of length two specifying the x coordinate -limits of enclosing box.} - -\item{w_y}{A numeric vector of length two specifying the y coordinate -limits of enclosing box.} +\item{range_list}{A named list of spatial ranges for each sample. +Each element should be a list with two components: +\code{w_x} and \code{w_y}, which are numeric vectors of length 2 +specifying the x- and y-axis ranges (e.g., from cell or +transcript coordinates). +The range is calculated with 5% buffer to ensure all points fall +within the window.} } \value{ a matrix contains the transcript count in each grid. diff --git a/man/get_cor.Rd b/man/get_cor.Rd index c47eeab..43bca55 100644 --- a/man/get_cor.Rd +++ b/man/get_cor.Rd @@ -50,14 +50,6 @@ mol <- BumpyMatrix::splitAsBumpyMatrix( row = trans_info$feature_name, col = trans_info$cell ) spe_sample1 <- SpatialExperiment( assays = list(molecules = mol),sample_id ="sample1" ) -w_x <- c(min(floor(min(trans_info$x)), - floor(min(clusters$x))), - max(ceiling(max(trans_info$x)), - ceiling(max(clusters$x)))) -w_y <- c(min(floor(min(trans_info$y)), - floor(min(clusters$y))), - max(ceiling(max(trans_info$y)), - ceiling(max(clusters$y)))) set.seed(100) corr_res <- compute_permp(x=spe_sample1, cluster_info=clusters, @@ -67,9 +59,7 @@ corr_res <- compute_permp(x=spe_sample1, test_genes=unique(trans_info$feature_name), correlation_method = "pearson", n_cores=1, - correction_method="BH", - w_x=w_x , - w_y=w_y) + correction_method="BH") # observed correlation obs_corr <- get_cor(corr_res) } diff --git a/man/get_full_mg.Rd b/man/get_full_mg.Rd index 8d5c6ab..2764559 100644 --- a/man/get_full_mg.Rd +++ b/man/get_full_mg.Rd @@ -69,20 +69,11 @@ trans_mol <- BumpyMatrix::splitAsBumpyMatrix( row = trans_info$feature_name, col = trans_info$cell ) spe<- SpatialExperiment( assays = list(molecules = trans_mol),sample_id ="sample1" ) -w_x <- c(min(floor(min(trans_info$x)), - floor(min(clusters$x))), - max(ceiling(max(trans_info$x)), - ceiling(max(clusters$x)))) -w_y <- c(min(floor(min(trans_info$y)), - floor(min(clusters$y))), - max(ceiling(max(trans_info$y)), - ceiling(max(clusters$y)))) vecs_lst <- get_vectors(x=spe,sample_names=c("sample1"), cluster_info = clusters, bin_type = "square", bin_param = c(20,20), - test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2"), - w_x = w_x, w_y=w_y) + test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2")) lasso_res <- lasso_markers(gene_mt=vecs_lst$gene_mt, cluster_mt = vecs_lst$cluster_mt, sample_names=c("sample1"), diff --git a/man/get_perm_adjp.Rd b/man/get_perm_adjp.Rd index a516025..9e2e317 100644 --- a/man/get_perm_adjp.Rd +++ b/man/get_perm_adjp.Rd @@ -49,14 +49,6 @@ mol <- BumpyMatrix::splitAsBumpyMatrix( row = trans_info$feature_name, col = trans_info$cell ) spe_sample1 <- SpatialExperiment( assays = list(molecules = mol),sample_id ="sample1" ) -w_x <- c(min(floor(min(trans_info$x)), - floor(min(clusters$x))), - max(ceiling(max(trans_info$x)), - ceiling(max(clusters$x)))) -w_y <- c(min(floor(min(trans_info$y)), - floor(min(clusters$y))), - max(ceiling(max(trans_info$y)), - ceiling(max(clusters$y)))) set.seed(100) corr_res <- compute_permp(x=spe_sample1, cluster_info=clusters, @@ -66,9 +58,7 @@ corr_res <- compute_permp(x=spe_sample1, test_genes=unique(trans_info$feature_name), correlation_method = "pearson", n_cores=1, - correction_method="BH", - w_x=w_x , - w_y=w_y) + correction_method="BH") # adjusted permutation p-value adjusted_perm_p <- get_perm_adjp(corr_res) } diff --git a/man/get_perm_p.Rd b/man/get_perm_p.Rd index 1cd6492..8cd71e0 100644 --- a/man/get_perm_p.Rd +++ b/man/get_perm_p.Rd @@ -49,14 +49,6 @@ mol <- BumpyMatrix::splitAsBumpyMatrix( row = trans_info$feature_name, col = trans_info$cell ) spe_sample1 <- SpatialExperiment( assays = list(molecules = mol),sample_id ="sample1" ) -w_x <- c(min(floor(min(trans_info$x)), - floor(min(clusters$x))), - max(ceiling(max(trans_info$x)), - ceiling(max(clusters$x)))) -w_y <- c(min(floor(min(trans_info$y)), - floor(min(clusters$y))), - max(ceiling(max(trans_info$y)), - ceiling(max(clusters$y)))) set.seed(100) corr_res <- compute_permp(x=spe_sample1, cluster_info=clusters, @@ -66,9 +58,7 @@ corr_res <- compute_permp(x=spe_sample1, test_genes=unique(trans_info$feature_name), correlation_method = "pearson", n_cores=1, - correction_method="BH", - w_x=w_x , - w_y=w_y) + correction_method="BH") # raw permutation p-value perm_p <- get_perm_p(corr_res) diff --git a/man/get_top_mg.Rd b/man/get_top_mg.Rd index c82d8f8..56c1831 100644 --- a/man/get_top_mg.Rd +++ b/man/get_top_mg.Rd @@ -69,20 +69,11 @@ trans_mol <- BumpyMatrix::splitAsBumpyMatrix( row = trans_info$feature_name, col = trans_info$cell ) spe<- SpatialExperiment( assays = list(molecules = trans_mol),sample_id ="sample1" ) -w_x <- c(min(floor(min(trans_info$x)), - floor(min(clusters$x))), - max(ceiling(max(trans_info$x)), - ceiling(max(clusters$x)))) -w_y <- c(min(floor(min(trans_info$y)), - floor(min(clusters$y))), - max(ceiling(max(trans_info$y)), - ceiling(max(clusters$y)))) vecs_lst <- get_vectors(x=spe,sample_names=c("sample1"), cluster_info = clusters, bin_type = "square", bin_param = c(20,20), - test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2"), - w_x = w_x, w_y=w_y) + test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2")) lasso_res <- lasso_markers(gene_mt=vecs_lst$gene_mt, cluster_mt = vecs_lst$cluster_mt, sample_names=c("sample1"), diff --git a/man/get_vectors.Rd b/man/get_vectors.Rd index 4269478..9b4d039 100644 --- a/man/get_vectors.Rd +++ b/man/get_vectors.Rd @@ -11,8 +11,6 @@ get_vectors( bin_type, bin_param, test_genes, - w_x, - w_y, use_cm = FALSE, n_cores = 1 ) @@ -54,12 +52,6 @@ length of hexagons. Positive numbers only.} to create gene vector. This will be used as column names for one of the result matrix \code{gene_mt}.} -\item{w_x}{A numeric vector of length two specifying the x coordinate -limits of enclosing box.} - -\item{w_y}{A numeric vector of length two specifying the y coordinate -limits of enclosing box.} - \item{use_cm}{A boolean value that specifies whether to create spatial vectors for genes using the count matrix and cell coordinates instead of the transcript coordinates when both types of information are available. @@ -135,22 +127,13 @@ trans_info$x=as.numeric(trans_info$x) trans_info$y=as.numeric(trans_info$y) trans_info$cell = sample(c("cell1","cell2","cell2"),replace=TRUE, size=nrow(trans_info)) -w_x = c(min(floor(min(trans_info$x)), - floor(min(clusters$x))), - max(ceiling(max(trans_info$x)), - ceiling(max(clusters$x)))) -w_y = c(min(floor(min(trans_info$y)), - floor(min(clusters$y))), - max(ceiling(max(trans_info$y)), - ceiling(max(clusters$y)))) # use named list as input vecs_lst = get_vectors(x= list("sample1" = trans_info), sample_names=c("sample1"), cluster_info = clusters, bin_type = "square", bin_param = c(5,5), - test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2"), - w_x = w_x, w_y=w_y) + test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2")) # use SpatialExperiment object as input trans_mol <- BumpyMatrix::splitAsBumpyMatrix( trans_info[, c("x", "y")], @@ -161,7 +144,6 @@ vecs_lst_spe = get_vectors(x=spe,sample_names=c("sample1"), cluster_info = clusters, bin_type = "square", bin_param = c(5,5), - test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2"), - w_x = w_x, w_y=w_y) + test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2")) } diff --git a/man/lasso_markers.Rd b/man/lasso_markers.Rd index a6c99af..5b30ecb 100644 --- a/man/lasso_markers.Rd +++ b/man/lasso_markers.Rd @@ -123,20 +123,11 @@ trans_mol <- BumpyMatrix::splitAsBumpyMatrix( row = trans_info$feature_name, col = trans_info$cell ) spe<- SpatialExperiment( assays = list(molecules = trans_mol),sample_id ="sample1" ) -w_x <- c(min(floor(min(trans_info$x)), - floor(min(clusters$x))), - max(ceiling(max(trans_info$x)), - ceiling(max(clusters$x)))) -w_y <- c(min(floor(min(trans_info$y)), - floor(min(clusters$y))), - max(ceiling(max(trans_info$y)), - ceiling(max(clusters$y)))) vecs_lst <- get_vectors(x=spe,sample_names=c("sample1"), cluster_info = clusters, bin_type = "square", bin_param = c(20,20), - test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2"), - w_x = w_x, w_y=w_y) + test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2")) lasso_res <- lasso_markers(gene_mt=vecs_lst$gene_mt, cluster_mt = vecs_lst$cluster_mt, sample_names=c("sample1"), diff --git a/tests/testthat/test-compute_permp.R b/tests/testthat/test-compute_permp.R index 62358ef..eaa360a 100644 --- a/tests/testthat/test-compute_permp.R +++ b/tests/testthat/test-compute_permp.R @@ -36,14 +36,6 @@ spe_rep2 <- SpatialExperiment( invalid_spe_rep1 <- SingleCellExperiment::cbind(spe_rep1,spe_rep2) -w_x = c(min(floor(min(trans_info$x)), - floor(min(clusters$x))), - max(ceiling(max(trans_info$x)), - ceiling(max(clusters$x)))) -w_y = c(min(floor(min(trans_info$y)), - floor(min(clusters$y))), - max(ceiling(max(trans_info$y)), - ceiling(max(clusters$y)))) test_that("Invaid input",{ expect_error(compute_permp(x=trans_info, @@ -54,9 +46,7 @@ test_that("Invaid input",{ test_genes=unique(trans_info$feature_name), correlation_method = "pearson", n_cores=2, - correction_method="BH", - w_x=w_x , - w_y=w_y)) + correction_method="BH")) expect_error(compute_permp(x=spe_rep1, cluster_info=clusters, @@ -66,9 +56,8 @@ test_that("Invaid input",{ test_genes=unique(trans_info$feature_name), correlation_method = "pearson", n_cores=2, - correction_method="BH", - w_x=w_x , - w_y=w_y)) + correction_method="BH" + )) expect_error(compute_permp(x=spe_rep1, cluster_info=clusters, perm.size=100, @@ -77,9 +66,8 @@ test_that("Invaid input",{ test_genes=unique(trans_info$feature_name), correlation_method = "pearson", n_cores=2, - correction_method="BH", - w_x=w_x , - w_y=w_y)) + correction_method="BH" + )) expect_error(compute_permp(x=spe_rep1, cluster_info=clusters, perm.size=100, @@ -88,9 +76,8 @@ test_that("Invaid input",{ test_genes=unique(trans_info$feature_name), correlation_method = "pearson", n_cores=2, - correction_method="BH", - w_x=w_x , - w_y=w_y)) + correction_method="BH" + )) expect_error(compute_permp(x=invalid_spe_rep1, cluster_info=clusters, @@ -100,9 +87,8 @@ test_that("Invaid input",{ test_genes=unique(trans_info$feature_name), correlation_method = "pearson", n_cores=2, - correction_method="BH", - w_x=w_x , - w_y=w_y)) + correction_method="BH" + )) }) set.seed(100) perm_p_lst_input = compute_permp(x=list("rep1"=trans_info), @@ -113,9 +99,7 @@ perm_p_lst_input = compute_permp(x=list("rep1"=trans_info), test_genes=unique(trans_info$feature_name), correlation_method = "pearson", n_cores=2, - correction_method="BH", - w_x=w_x , - w_y=w_y) + correction_method="BH") perm_p_lst = compute_permp(x=spe_rep1, cluster_info=clusters, @@ -125,9 +109,7 @@ perm_p_lst = compute_permp(x=spe_rep1, test_genes=unique(trans_info$feature_name), correlation_method = "pearson", n_cores=2, - correction_method="BH", - w_x=w_x , - w_y=w_y) + correction_method="BH") perm_p_s = compute_permp(x=spe_rep1, cluster_info=clusters, perm.size=10, @@ -136,9 +118,7 @@ perm_p_s = compute_permp(x=spe_rep1, test_genes=unique(trans_info$feature_name), correlation_method = "pearson", n_cores=1, - correction_method="BH", - w_x=w_x , - w_y=w_y) + correction_method="BH") test_that("Test permutation result - output dimension matches", { #expect_equal(length(perm_p_lst), 4) #expect_equal(dim(perm_p_lst$perm.arrays), c(4,2,10)) @@ -212,8 +192,8 @@ clusters = data.frame(x = c(1, 2,20,21,22,23,24), clusters$sample="rep1" clusters$cell_id= colnames(cm) # simulate coordiantes for genes -w_x=c(0,25) -w_y=c(0,25) +# w_x=c(0,25) +# w_y=c(0,25) # cell_1 = (1,0,0,0) # cell_2 = (1,0,0,0) # cell_3 = (0,0,0,1) @@ -231,10 +211,8 @@ perm_p_s = compute_permp(x=sce, bin_param=c(2,2), test_genes=row.names(cm), correlation_method = "pearson", - n_cores=1, - correction_method="BH", - w_x=w_x , - w_y=w_y) + n_cores=1,use_cm = TRUE, + correction_method="BH") test_that("Test permutation result - output dimension matches", { #expect_equal(length(perm_p_s), 4) @@ -257,9 +235,7 @@ test_that("Invaid input",{ test_genes=row.names(cm), correlation_method = "pearson", n_cores=1, - correction_method="BH", - w_x=w_x , - w_y=w_y)) + correction_method="BH")) expect_error(compute_permp(x=invalid_sce, cluster_info=clusters, perm.size=10, @@ -268,7 +244,5 @@ test_that("Invaid input",{ test_genes=row.names(cm), correlation_method = "pearson", n_cores=1, - correction_method="BH", - w_x=w_x , - w_y=w_y)) + correction_method="BH")) }) diff --git a/tests/testthat/test-create_genesets.R b/tests/testthat/test-create_genesets.R index 641b7a4..4940f5b 100644 --- a/tests/testthat/test-create_genesets.R +++ b/tests/testthat/test-create_genesets.R @@ -18,8 +18,7 @@ test_that("Invaid input",{ expect_error(create_genesets(x=list("rep1" = trans_invalid),sample_names = "rep1", name_lst=list(dummy_W=c("A","B")), bin_type="square", - bin_param=c(2,2),cluster_info = NULL, - w_x=c(0,25), w_y=c(0,25))) + bin_param=c(2,2),cluster_info = NULL)) }) mol <- BumpyMatrix::splitAsBumpyMatrix( trans[, c("x", "y")], @@ -43,35 +42,29 @@ spe_invalid <- SpatialExperiment( geneset_res_A_fl = create_genesets(x=list("rep1" = trans),sample_names = "rep1", name_lst=list(dummy_W=c("A")), bin_type="square", - bin_param=c(2,2),cluster_info = NULL, - w_x=c(0,25), w_y=c(0,25)) + bin_param=c(2,2),cluster_info = NULL) geneset_res_A = create_genesets(x=spe_rep1,sample_names = "rep1", name_lst=list(dummy_W=c("A")), bin_type="square", - bin_param=c(2,2),cluster_info = NULL, - w_x=c(0,25), w_y=c(0,25)) + bin_param=c(2,2),cluster_info = NULL) geneset_res_A_hex = create_genesets(x=spe_rep1,sample_names = "rep1", name_lst=list(dummy_W=c("A")), bin_type="hexagon", - bin_param=c(10),cluster_info = NULL, - w_x=c(0,25), w_y=c(0,25)) + bin_param=c(10),cluster_info = NULL) test_that("Invaid input",{ expect_error(create_genesets(x=spe_rep1,sample_names = "rep1", name_lst=list(dummy_W=c("A")), bin_type="hexagon", - bin_param=c(2,2),cluster_info = NULL, - w_x=c(0,25), w_y=c(0,25))) + bin_param=c(2,2),cluster_info = NULL)) expect_error(create_genesets(x=spe_rep1,sample_names = "rep1", name_lst=list(dummy_W=c("A")), bin_type="circle", - bin_param=c(2,2),cluster_info = NULL, - w_x=c(0,25), w_y=c(0,25))) + bin_param=c(2,2),cluster_info = NULL)) expect_error(create_genesets(x=spe_rep1,sample_names = "rep1", name_lst=list(dummy_W=c("A")), bin_type="square", - bin_param=c(2),cluster_info = NULL, - w_x=c(0,25), w_y=c(0,25))) + bin_param=c(2),cluster_info = NULL)) }) test_that("Test can create vectors for for a list input", { @@ -80,18 +73,33 @@ test_that("Test can create vectors for for a list input", { }) +# +# x_rng <- range(trans$x, na.rm = TRUE) +# y_rng <- range(trans$y, na.rm = TRUE) +# +# # 1% padding on each side +# pad_x <- 0.05 * diff(x_rng) +# pad_y <- 0.05 * diff(y_rng) +# +# w_x = c(x_rng[1] - pad_x, x_rng[2] + pad_x) +# w_y = c(y_rng[1] - pad_y, y_rng[2] + pad_y) +# +# points_ppp <- ppp(trans[trans$feature_name=="A","x"], trans[trans$feature_name=="A","y"], +# window = owin(xrange = w_x,yrange =w_y)) +# H <- hextess(W = points_ppp$window, 10) +# vec_g <- as.vector(t(quadratcount(points_ppp, tess = H))) +# vec_g test_that("Test can create vectors for a single gene set with a single element A", { expect_equal(colnames(geneset_res_A), c("dummy_W")) expect_equal(as.vector(geneset_res_A$dummy_W), c(0,0,10, 0)) - expect_equal(as.vector(geneset_res_A_hex$dummy_W), c(2,5,0,3,0,0,0)) + expect_equal(as.vector(geneset_res_A_hex$dummy_W), c(2,6,0,2,0,0,0)) }) geneset_res_C = create_genesets(x=spe_rep1,sample_names = "rep1", name_lst=list(dummy_W=c("C")), bin_type="square", - bin_param=c(2,2),cluster_info = NULL, - w_x=c(0,25), w_y=c(0,25)) + bin_param=c(2,2),cluster_info = NULL) test_that("Test can create vectors for a single gene set with a single element C", { expect_equal(colnames(geneset_res_C), c("dummy_W")) @@ -101,13 +109,11 @@ test_that("Test can create vectors for a single gene set with a single element C geneset_res1 = create_genesets(x=spe_rep1,sample_names = "rep1", name_lst=list(dummy_W=c("A","B")), bin_type="square", - bin_param=c(2,2),cluster_info = NULL, - w_x=c(0,25), w_y=c(0,25)) + bin_param=c(2,2),cluster_info = NULL) geneset_res1_lst = create_genesets(x=list("rep1" = trans),sample_names = "rep1", name_lst=list(dummy_W=c("A","B")), bin_type="square", - bin_param=c(2,2),cluster_info = NULL, - w_x=c(0,25), w_y=c(0,25)) + bin_param=c(2,2),cluster_info = NULL) test_that("Test can create vectors for a single gene set with multiple elements", { expect_equal(colnames(geneset_res1), c("dummy_W")) expect_equal(as.vector(geneset_res1$dummy_W), c(0,0,11,4)) @@ -118,8 +124,7 @@ geneset_res2 = create_genesets(x=spe_rep1,sample_names = "rep1", name_lst=list(dummy_A=c("A","C"), dummy_B=c("A","B","C")), bin_type="square", - bin_param=c(2,2),cluster_info = NULL, - w_x=c(0,25), w_y=c(0,25)) + bin_param=c(2,2),cluster_info = NULL) test_that("Test can create vectors for gene sets- output mathces", { expect_equal(colnames(geneset_res2), c("dummy_A","dummy_B")) @@ -132,8 +137,7 @@ geneset_res3 = create_genesets(x=spe_rep1,sample_names = "rep1", name_lst=list(dummy_A=c("A","B","C","C"), dummy_B=c("A","B","C")), bin_type="square",cluster_info = NULL, - bin_param=c(2,2), - w_x=c(0,25), w_y=c(0,25)) + bin_param=c(2,2)) test_that("Test can create vectors for duplicated gene sets", { expect_equal(colnames(geneset_res3), c("dummy_A","dummy_B")) @@ -144,15 +148,13 @@ test_that("Test can create vectors for duplicated gene sets", { geneset_res4 = create_genesets(x=spe_rep1,sample_names = "rep1", name_lst=list(dummy_B=c("A","B","C")), bin_type="square", - bin_param=c(2,2),cluster_info = NULL, - w_x=c(0,25), w_y=c(0,25)) + bin_param=c(2,2),cluster_info = NULL) test_that("Test can not create vectors for gene sets with non-overlapped genes", { expect_error(create_genesets(x=spe_rep1,sample_names = "rep1", name_lst=list(dummy_A=c("F","G")), bin_type="square", - bin_param=c(2,2),cluster_info = NULL, - w_x=c(0,25), w_y=c(0,25))) + bin_param=c(2,2),cluster_info = NULL)) expect_equal(colnames(geneset_res4), c("dummy_B")) expect_equal(as.vector(geneset_res4$dummy_B), c(2,7,11,5)) }) @@ -172,15 +174,14 @@ clusters$sample="sample1" clusters$cell_id= colnames(cm) sce <- SingleCellExperiment(list(sample1=cm)) # simulate coordiantes for genes -w_x=c(0,25) -w_y=c(0,25) +# w_x=c(0,25) +# w_y=c(0,25) geneset_res_cm = create_genesets(x=sce,sample_names = "sample1", name_lst=list(dummy_A=c("gene_A","gene_B"), dummy_B=c("gene_C","gene_D")), bin_type="square", - bin_param=c(2,2), - cluster_info = clusters, - w_x=w_x, w_y=w_y) + bin_param=c(2,2),use_cm = TRUE, + cluster_info = clusters) test_that("Test can define gene sets from count matrix", { expect_equal(dim(geneset_res_cm), c(4,2)) @@ -191,87 +192,64 @@ test_that("Test can define gene sets from count matrix", { dummy_B=c("gene_C","gene_D")), bin_type="square", bin_param=c(2,2), - cluster_info = NULL, - w_x=w_x, w_y=w_y)) + cluster_info = NULL)) expect_error(create_genesets(x=sce,sample_names = "sample1", name_lst=list(dummy_A=c("gene_A","gene_B"), dummy_B=c("gene_C","gene_D")), bin_type="square", bin_param=c(2,2), - cluster_info = invalid_cluster, - w_x=w_x, w_y=w_y)) + cluster_info = invalid_cluster)) }) ############################################################################ # example of two SPE objects but with different gene sets -# real world exmple: if the two samples have slightly different negative controls probes + s1 =as.data.frame(cbind(feature_name = c("A","B","C"), x= c(1,2,3), y= c(24, 22, 23), cell_id = c("cell1","cell2","cell3"))) s2 = as.data.frame(cbind(feature_name = c("A","E","F","G"), - x= c(1,2,3,4), - y= c(24, 22, 23,24), + x= c(1,2,3,3), + y= c(24,22,23,24), cell_id = c("cell4","cell5","cell6","cell7"))) -shared_nc = intersect(s1$feature_name,s2$feature_name) - -s1_added = cbind(x=NA, y=NA, - feature_name=setdiff(s2$feature_name,shared_nc), - cell_id = NA) -s2_added = cbind(x=NA, y=NA, - feature_name=setdiff(s1$feature_name,shared_nc), - cell_id = NA) - -s1 = rbind(s1,s1_added) -s2 = rbind(s2,s2_added) - -s1$category = "probe" -s1[s1$feature_name %in% c("A","E"), "category"] = "codeword" - -s2$category = "probe" -s2[s2$feature_name %in% c("C","G"), "category"] = "codeword" - -s1_mol <- BumpyMatrix::splitAsBumpyMatrix( - s1[, c("x", "y","category")], - row = s1$feature_name, col = s1$cell_id ) - -s1_spe<- SpatialExperiment( - assays = list(molecules = s1_mol),sample_id ="sample1") - -s2_mol <- BumpyMatrix::splitAsBumpyMatrix( - s2[, c("x", "y","category")], - row = s2$feature_name, col = s2$cell_id ) +s2$x = as.numeric(s2$x) +s2$y = as.numeric(s2$y) +s1$x = as.numeric(s1$x) +s1$y = as.numeric(s1$y) -s2_spe<- SpatialExperiment( - assays = list(molecules = s2_mol),sample_id ="sample2") +# library(spatstat.geom) +# w <- owin(xrange = c(0.9, 3.1), yrange = c(21.9, 24.1)) +# pp <- ppp(x = s1[s1$feature_name=="B", "x"], y = s1[s1$feature_name=="B", "y"], +# window = w) +# qc <- quadratcount(pp, nx = 2, ny = 2) +# as.vector(t(qc)) -combined_spe = SingleCellExperiment::cbind(s1_spe, s2_spe) -vecs_lst = get_vectors(x=combined_spe,sample_names=c("sample1","sample2"), +vecs_lst = get_vectors(x=list(sample1= s1, + sample2 =s2 ), + sample_names=c("sample1","sample2"), bin_type="square",test_genes = c("A","B","C","E","F","G"), bin_param=c(2,2), - w_x=w_x, w_y=w_y, cluster_info = NULL) test_that("Test can not create vectors for multi-samples with different genes", { expect_equal(as.vector(vecs_lst$gene_mt[,"A"]), c(1,0,0,0,1,0,0,0)) - expect_equal(as.vector(vecs_lst$gene_mt[,"B"]), c(1,0,0,0,0,0,0,0)) - expect_equal(as.vector(vecs_lst$gene_mt[,"C"]), c(1,0,0,0,0,0,0,0)) - expect_equal(as.vector(vecs_lst$gene_mt[,"E"]), c(0,0,0,0,1,0,0,0)) - expect_equal(as.vector(vecs_lst$gene_mt[,"F"]), c(0,0,0,0,1,0,0,0)) - expect_equal(as.vector(vecs_lst$gene_mt[,"G"]), c(0,0,0,0,1,0,0,0)) + expect_equal(as.vector(vecs_lst$gene_mt[,"B"]), c(0,0,1,0,0,0,0,0)) + expect_equal(as.vector(vecs_lst$gene_mt[,"C"]), c(0,0,0,1,0,0,0,0)) + expect_equal(as.vector(vecs_lst$gene_mt[,"E"]), c(0,0,0,0,0,0,1,0)) + expect_equal(as.vector(vecs_lst$gene_mt[,"F"]), c(0,0,0,0,0,0,0,1)) + expect_equal(as.vector(vecs_lst$gene_mt[,"G"]), c(0,0,0,0,0,1,0,0)) }) -example_vec = create_genesets(x=combined_spe, +example_vec = create_genesets(x=list(sample1= s1,sample2 =s2 ), sample_names=c("sample1","sample2"), name_lst=list(probe=c("A","E"), codeword=c("C","G")), bin_type="square", bin_param=c(2,2), - w_x=w_x, w_y=w_y, cluster_info = NULL) test_that("Test can not create vectors for multi-samples with different genes", { - expect_equal(as.vector(example_vec[,"probe"]), c(1,0,0,0,2,0,0,0)) - expect_equal(as.vector(example_vec[,"codeword"]), c(1,0,0,0,1,0,0,0)) + expect_equal(as.vector(example_vec[,"probe"]), c(1,0,0,0,1,0,1,0)) + expect_equal(as.vector(example_vec[,"codeword"]), c(0,0,0,1,0,1,0,0)) }) diff --git a/tests/testthat/test-get_vectors.R b/tests/testthat/test-get_vectors.R index 99a9d00..018f252 100644 --- a/tests/testthat/test-get_vectors.R +++ b/tests/testthat/test-get_vectors.R @@ -10,13 +10,10 @@ clusters = data.frame(x = c(1,2,20,21,22,23,24), y = c(23, 24, 1,2,3,4,5), cluster="A") clusters$sample="rep1" -w_x=c(0,25) -w_y=c(0,25) vecs_lst = get_vectors(x= NULL, cluster_info = clusters, bin_type = "square",sample_names="rep1", bin_param = c(2,2), - test_genes = NULL, - w_x = w_x, w_y=w_y) + test_genes = NULL) test_that("Test can only vectorise clusters - output length matches", { expect_equal(length(vecs_lst), 1) @@ -34,13 +31,10 @@ clusters = data.frame(x = c(1,2,20,21,22,23,24), cluster=c("erqf","qre","$+9","9-", "123","erqf","123")) clusters$sample="rep1" -w_x=c(0,25) -w_y=c(0,25) vecs_lst = get_vectors(x= NULL, cluster_info = clusters, bin_type = "square",sample_names="rep1", bin_param = c(2,2), - test_genes = NULL, - w_x = w_x, w_y=w_y) + test_genes = NULL) test_that("Test can correct invalid cluster names- output vector matches", { # expect_message(get_vectors(x= NULL, cluster_info = clusters, @@ -58,13 +52,11 @@ clusters = data.frame(x = c(1,2,20,21,22,23,24), cluster=c("erqf","qre","$+9","9-", "M+","erqf","M+")) clusters$sample=c(rep("12",times=4),rep("s1",times=3)) -w_x=c(0,25) -w_y=c(0,25) + vecs_lst = get_vectors(x= NULL, cluster_info = clusters, bin_type = "square",sample_names=c("12","s1"), bin_param = c(2,2), - test_genes = NULL, - w_x = w_x, w_y=w_y) + test_genes = NULL) test_that("Test can correct invalid sample names- output vector matches", { # expect_message(vecs_lst = get_vectors(x= NULL, cluster_info = clusters, @@ -110,16 +102,14 @@ vecs_lst_gene = get_vectors(x= spe, cluster_info = NULL, bin_type = "square", bin_param = c(2,2), - test_genes = row.names(spe), - w_x = w_x, w_y=w_y) + test_genes = row.names(spe)) test_that("Invalid colnames for transctipt with list input", { expect_error(get_vectors(x= list(sample1=df), sample_names = c("sample1"), cluster_info = NULL, bin_type = "square", bin_param = c(2,2), - test_genes = gs, - w_x = w_x, w_y=w_y)) + test_genes = gs)) }) colnames(df) =c("feature_name","cell","x","y") vecs_lst_gene_lst = get_vectors(x= list(sample1=df), @@ -127,8 +117,7 @@ vecs_lst_gene_lst = get_vectors(x= list(sample1=df), cluster_info = NULL, bin_type = "square", bin_param = c(2,2), - test_genes = gs, - w_x = w_x, w_y=w_y) + test_genes = gs) test_that("Test can only vectorise genes - output length mathces", { expect_equal(length(vecs_lst_gene), 1) @@ -170,15 +159,13 @@ spe_invalid <- SpatialExperiment( assays = list(molecules = mol_inv),sample_id ="rep1" ) -w_x=c(0,25) w_y=c(0,25) vecs_lst_cluster = get_vectors(x= NULL, cluster_info = clusters, sample_names="sample1", bin_type = "square", bin_param = c(2,2), - test_genes = NULL, - w_x = w_x, w_y=w_y) + test_genes = NULL) spe <- SpatialExperiment( assays = list(molecules = mol),sample_id ="sample1" ) @@ -186,8 +173,7 @@ vecs_lst_gene = get_vectors(x= spe,sample_names="sample1" , cluster_info = NULL, bin_type = "square", bin_param = c(2,2), - test_genes = row.names(spe), - w_x = w_x, w_y=w_y) + test_genes = row.names(spe)) test_that("Test can vectorise genes and clusters - output mathces", { expect_equal(as.vector(vecs_lst_cluster$cluster_mt), c(2,0,2,3)) @@ -200,50 +186,43 @@ test_that("Invaid input",{ cluster_info = clusters, bin_type="square", bin_param=c(2,2), - test_genes =matrix(0, nrow=5, ncol=2), - w_x=c(0,25), w_y=c(0,25))) + test_genes =matrix(0, nrow=5, ncol=2))) expect_error(get_vectors(x=trans, sample_names="sample1", cluster_info = clusters, bin_type="hexagon", bin_param=c(2,2), - test_genes = row.names(spe), - w_x=c(0,25), w_y=c(0,25))) + test_genes = row.names(spe))) expect_error(get_vectors(x=NULL, sample_names="sample1", cluster_info = clusters, bin_type="hexagon", bin_param=c(2,2), - test_genes = row.names(spe), - w_x=c(0,25), w_y=c(0,25))) + test_genes = row.names(spe))) expect_error(get_vectors(x=NULL, sample_names="sample1", cluster_info = clusters, bin_type="circle", bin_param=c(2,2), - test_genes = row.names(spe), - w_x=c(0,25), w_y=c(0,25))) + test_genes = row.names(spe))) expect_error(get_vectors(x=NULL, sample_names="sample1", cluster_info = clusters, bin_type="square", - bin_param=c(2), - w_x=c(0,25), w_y=c(0,25))) + bin_param=c(2))) expect_error(get_vectors(x=NULL, cluster_info = NULL, sample_names="sample1", bin_type="hexagon", bin_param=c(2), - test_genes = row.names(spe), - w_x=c(0,25), w_y=c(0,25))) + test_genes = row.names(spe))) expect_error(get_vectors(x=NULL, sample_names="sample1", cluster_info=invalid_cluster, bin_type="square", bin_param=c(2,2), - test_genes = row.names(spe), - w_x=c(0,25), w_y=c(0,25))) + test_genes = row.names(spe))) expect_error(get_vectors(x=spe, sample_names="sample1", cluster_info=invalid_cluster, bin_type="square", bin_param=c(2,2), test_genes = row.names(spe), - w_x=c(0,25), w_y=c(0,25), use_cm=TRUE)) + use_cm=TRUE)) # can not match test_genes from x expect_error(get_vectors(x= spe, @@ -251,24 +230,21 @@ test_that("Invaid input",{ test_genes =c("geneA","geneB"), cluster_info=clusters, bin_type="square", - bin_param=c(2,2), - w_x=c(0,25), w_y=c(0,25))) + bin_param=c(2,2))) # mismatch sample_names and clsuetr_info$sample expect_error(get_vectors(x= NULL, cluster_info = clusters, sample_names="sample2", bin_type = "square", bin_param = c(2,2), - test_genes = NULL, - w_x = w_x, w_y=w_y)) + test_genes = NULL)) # missing y coordiantes for transcript expect_error(get_vectors(x= spe_invalid, cluster_info = NULL, sample_names="rep1", bin_type = "square", bin_param = c(2,2), - test_genes = unique(invalid_trans_df$feature_name), - w_x = w_x, w_y=w_y)) + test_genes = unique(invalid_trans_df$feature_name))) }) @@ -286,8 +262,8 @@ clusters = data.frame(x = c(1, 2,20,21,22,23,24), clusters$sample="rep1" clusters$cell_id= colnames(cm) # simulate coordiantes for genes -w_x=c(0,25) -w_y=c(0,25) +# w_x=c(0,25) +# w_y=c(0,25) # cell_1 = (1,0,0,0) # cell_2 = (1,0,0,0) # cell_3 = (0,0,0,1) @@ -300,8 +276,7 @@ vecs_lst = get_vectors(x= sce, cluster_info = clusters, sample_names = "rep1", bin_type = "square", bin_param = c(2,2), - test_genes = row.names(cm), - w_x = w_x, w_y=w_y) + test_genes = row.names(cm)) test_that("Test can use count matrix (square) - output vector matches", { expect_equal(as.vector(vecs_lst$gene_mt), @@ -316,8 +291,7 @@ test_that("Invaid input",{ expect_error(get_vectors(x=sce,sample_names = "rep1", cluster_info = NULL, bin_type="square", - bin_param=c(2,2),test_genes = row.names(cm), - w_x=c(0,25), w_y=c(0,25))) + bin_param=c(2,2),test_genes = row.names(cm))) }) @@ -335,8 +309,8 @@ clusters = data.frame(x = c(1, 2,20,21,22,23,24), clusters$sample="rep1" clusters$cell_id= colnames(cm) # simulate coordinates for genes -w_x=c(0,25) -w_y=c(0,25) +# w_x=c(0,25) +# w_y=c(0,25) # cell_1 = (1,0,0,0) # cell_2 = (1,0,0,0) # cell_3 = (0,0,0,1) @@ -391,15 +365,14 @@ vecs_lst_cm = get_vectors(x= sce_m, cluster_info = clusters, sample_names = "rep1", bin_type = "square", bin_param = c(2,2), - test_genes = row.names(cm), - w_x = w_x, w_y=w_y) + test_genes = row.names(cm) + ) vecs_lst_tr = get_vectors(x= spe_m, cluster_info = clusters, sample_names = "rep1", bin_type = "square", bin_param = c(2,2), - test_genes = row.names(cm), - w_x = w_x, w_y=w_y) + test_genes = row.names(cm)) test_that("Test count matrix vectors match with transcript vectors - gene vectors", { expect_equal(as.vector(vecs_lst_tr$gene_mt), as.vector(vecs_lst_cm$gene_mt)) @@ -415,8 +388,7 @@ vecs_lst = get_vectors(x= sce_m, cluster_info = clusters, sample_names = "rep1", bin_type = "hexagon", bin_param = c(10), - test_genes = row.names(cm), - w_x = w_x, w_y=w_y) + test_genes = row.names(cm)) test_that("Test can use count matrix (hex) - output vector matches", { expect_equal(as.vector(vecs_lst$gene_mt), @@ -453,14 +425,13 @@ spe_rep1 <- SpatialExperiment( clusters = data.frame(x = c(3, 5,11,21,2,23,19), y = c(20, 24, 1,2,3,4,5), cluster="cluster_1") clusters$sample="rep1" -w_x=c(0,25) -w_y=c(0,25) + + vecs_lst_full = get_vectors(x=spe_rep1,sample_names = "rep1", cluster_info = clusters, bin_type = "square", bin_param = c(2,2), - test_genes = c("A","B","C"), - w_x = w_x, w_y=w_y,n_cores = 2) + test_genes = c("A","B","C"),n_cores = 2) test_that("Test can vectorise genes and clusters - output mathces", { @@ -508,15 +479,13 @@ vector_lst_1core = get_vectors(x= spe_rep2, sample_names = "rep2", cluster_info = NULL, bin_type="square", bin_param=c(2,2), - test_genes = feature_names, - w_x=c(0,100), w_y=c(0,100), n_cores = 1) + test_genes = feature_names,n_cores = 1) vector_lst_5core = get_vectors(x= spe_rep2, sample_names = "rep2", cluster_info = NULL, bin_type="square", bin_param=c(2,2), - test_genes = feature_names, - w_x=c(0,100), w_y=c(0,100), n_cores = 2) + test_genes = feature_names,n_cores = 2) test_that("Test can result from sequential matches with result from parallel",{ expect_equal(as.vector(vector_lst_1core$gene_mt), as.vector(vector_lst_5core$gene_mt)) @@ -556,7 +525,6 @@ vector_lst_twosample = get_vectors(x= sce_two, sample_names =c("sp1","sp2"), bin_type="square", bin_param=c(2,2), test_genes = paste("gene", 1:10, sep=""), - w_x=c(0,25), w_y=c(0,25), n_cores = 1) test_that("Test can result from sequential matches with result from parallel",{ expect_equal(as.vector(vector_lst_twosample$cluster_mt[,"A"]), @@ -589,15 +557,12 @@ clusters = data.frame(x = c(3, 5,11,21,2,23,19), y = c(20, 24, 1,2,3,4,5), cluster="cluster_1") set.seed(98) clusters$sample=sample(c("s1","s2"),size = 7,replace = TRUE) -w_x=c(0,25) -w_y=c(0,25) vecs_lst = get_vectors(x=list("s1"=trans, "s2"=trans), sample_names = c("s1","s2"), cluster_info = clusters, bin_type = "square", bin_param = c(2,2), - test_genes = c("A","B","C"), - w_x = w_x, w_y=w_y,n_cores = 1) + test_genes = c("A","B","C"),n_cores = 1) test_that("Test can work for a list input with 2 samples",{ expect_equal(as.vector(vecs_lst$cluster_mt[,"cluster_1"]), c(0,0,1,2,2,0,1,1)) diff --git a/tests/testthat/test-lasso_markers.R b/tests/testthat/test-lasso_markers.R index 79fdd52..2b12eb3 100644 --- a/tests/testthat/test-lasso_markers.R +++ b/tests/testthat/test-lasso_markers.R @@ -32,27 +32,17 @@ trans_mol <- BumpyMatrix::splitAsBumpyMatrix( spe<- SpatialExperiment( assays = list(molecules = trans_mol),sample_id ="sample1" ) -w_x = c(min(floor(min(trans_info$x)), - floor(min(clusters$x))), - max(ceiling(max(trans_info$x)), - ceiling(max(clusters$x)))) -w_y = c(min(floor(min(trans_info$y)), - floor(min(clusters$y))), - max(ceiling(max(trans_info$y)), - ceiling(max(clusters$y)))) vecs_lst = get_vectors(x=spe, cluster_info = clusters, bin_type = "square",sample_names = "sample1", bin_param = c(20,20), - test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2"), - w_x = w_x, w_y=w_y) + test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2")) #### background background_sv = create_genesets(x =spe,sample_names = "sample1", name_lst=list(dummy_W=c("gene_A1","gene_B1")), bin_type="square", - bin_param = c(20,20), - w_x = w_x, w_y=w_y,cluster_info = NULL) + bin_param = c(20,20),cluster_info = NULL) set.seed(100) lasso_res1 = lasso_markers(gene_mt=vecs_lst$gene_mt, cluster_mt = vecs_lst$cluster_mt, @@ -282,7 +272,7 @@ one_gene_res1 = lasso_markers(gene_mt=one_gene_mt, sample_names=c("rep1"), keep_positive=TRUE, background=NULL) -one_gene_top = get_top_mg(one_gene_res1, coef_cutoff = 1) +one_gene_top = get_top_mg(one_gene_res1, coef_cutoff = 2) one_gene_full =get_full_mg(one_gene_res1, coef_cutoff = 0) test_that("high resolution results in 0 significant genes", { expect_equal(unique(one_gene_top[one_gene_top$gene == "gene_A", @@ -327,30 +317,28 @@ sp2_trans_mol <- BumpyMatrix::splitAsBumpyMatrix( spe_sp2<- SpatialExperiment( assays = list(molecules = sp2_trans_mol),sample_id ="sample2" ) - -w_x = c(min(floor(min(trans_info$x)),floor(min(trans_info_sp2$x)), - floor(min(clusters$x))), - max(ceiling(max(trans_info$x)),ceiling(max(trans_info_sp2$x)), - ceiling(max(clusters$x)))) -w_y = c(min(floor(min(trans_info$y)),floor(min(trans_info_sp2$y)), - floor(min(clusters$y))), - max(ceiling(max(trans_info$y)),ceiling(max(trans_info_sp2$y)), - ceiling(max(clusters$y)))) +# +# w_x = c(min(floor(min(trans_info$x)),floor(min(trans_info_sp2$x)), +# floor(min(clusters$x))), +# max(ceiling(max(trans_info$x)),ceiling(max(trans_info_sp2$x)), +# ceiling(max(clusters$x)))) +# w_y = c(min(floor(min(trans_info$y)),floor(min(trans_info_sp2$y)), +# floor(min(clusters$y))), +# max(ceiling(max(trans_info$y)),ceiling(max(trans_info_sp2$y)), +# ceiling(max(clusters$y)))) twosample = SingleCellExperiment::cbind(spe,spe_sp2) vecs_lst = get_vectors(x=twosample, sample_names = c("sample1","sample2"), cluster_info = clusters, bin_type = "square", bin_param = c(20,20), - test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2"), - w_x = w_x, w_y=w_y) + test_genes =c("gene_A1","gene_A2","gene_B1","gene_B2")) #### background background_sv = create_genesets(x=twosample, sample_names = c("sample1","sample2"), name_lst=list(dummy_W=c("gene_A1","gene_B1")), bin_type="square", - bin_param = c(20,20), - w_x = w_x, w_y=w_y,cluster_info = NULL) + bin_param = c(20,20), cluster_info = NULL) set.seed(100) lasso_res_background = lasso_markers(gene_mt=vecs_lst$gene_mt, cluster_mt = vecs_lst$cluster_mt, @@ -381,7 +369,6 @@ background_sv_all = create_genesets(x=twosample, dummy_B2=c("gene_B2")), bin_type="square", bin_param = c(20,20), - w_x = w_x, w_y=w_y, cluster_info = NULL) lasso_nosig_bg = lasso_markers(gene_mt=vecs_lst$gene_mt, cluster_mt = vecs_lst$cluster_mt, diff --git a/vignettes/jazzPanda.Rmd b/vignettes/jazzPanda.Rmd index 0218a8e..fb5792f 100644 --- a/vignettes/jazzPanda.Rmd +++ b/vignettes/jazzPanda.Rmd @@ -342,15 +342,6 @@ specified by the `bin_type` parameter. ```{r create rep1 all vectors} seed_number<- 589 -w_x <- c(min(floor(min(rep1_transcripts$x)), - floor(min(rep1_clusters$x))), - max(ceiling(max(rep1_transcripts$x)), - ceiling(max(rep1_clusters$x)))) -w_y <- c(min(floor(min(rep1_transcripts$y)), - floor(min(rep1_clusters$y))), - max(ceiling(max(rep1_transcripts$y)), - ceiling(max(rep1_clusters$y)))) - grid_length <- 10 # get spatial vectors rep1_sq10_vectors <- get_vectors(x= rep1_sub, @@ -358,8 +349,7 @@ rep1_sq10_vectors <- get_vectors(x= rep1_sub, cluster_info = rep1_clusters, bin_type="square", bin_param=c(grid_length,grid_length), - test_genes = all_real_genes , - w_x=w_x, w_y=w_y) + test_genes = all_real_genes) ``` If your data has different scales (e.g.,` w_x` and `w_y`) for cluster coordinates and gene coordinates, you can construct separate spatial vectors @@ -368,9 +358,6 @@ for clusters and genes. Please make sure you use consistent bin parameters ```{r create rep1 different cluster-vec and gene vec} seed_number<- 589 grid_length <- 10 -# define a window for cell detections -cluster_w_x <- c(floor(min(rep1_clusters$x)), ceiling(max(rep1_clusters$x))) -cluster_w_y <- c(floor(min(rep1_clusters$y)), ceiling(max(rep1_clusters$y))) # get spatial vectors for clusters only cluster_sq10_vectors <- get_vectors(x= NULL, @@ -378,26 +365,19 @@ cluster_sq10_vectors <- get_vectors(x= NULL, cluster_info = rep1_clusters, bin_type="square", bin_param=c(grid_length,grid_length), - test_genes = all_real_genes , - w_x=cluster_w_x, w_y=cluster_w_y) + test_genes = all_real_genes) # the gene vectors can be accessed as: head(cluster_sq10_vectors$cluster_mt) # cluster_sq10_vectors$cluster_mt is equivalent to rep1_sq10_vectors$cluster_mt ####################################################### -# define a different window for transctipt detections -transcript_w_x <-c(floor(min(rep1_transcripts$x)), - ceiling(max(rep1_transcripts$x))) -transcript_w_y <- c(floor(min(rep1_transcripts$y)), - ceiling(max(rep1_transcripts$y))) # get spatial vectors for genes only gene_sq10_vectors <- get_vectors(x= rep1_sub, sample_names = "rep1", cluster_info = NULL, bin_type="square", bin_param=c(grid_length,grid_length), - test_genes = all_real_genes , - w_x=transcript_w_x, w_y=transcript_w_y) + test_genes = all_real_genes) # the gene vectors can be accessed as: head(gene_sq10_vectors$gene_mt) # gene_sq10_vectors$gene_mt is equivalent to rep1_sq10_vectors$gene_mt @@ -504,15 +484,6 @@ defining a SpatialExperiment object for your data. ```{r} seed_number<- 589 -w_x <- c(min(floor(min(rep1_transcripts$x)), - floor(min(rep1_clusters$x))), - max(ceiling(max(rep1_transcripts$x)), - ceiling(max(rep1_clusters$x)))) -w_y <- c(min(floor(min(rep1_transcripts$y)), - floor(min(rep1_clusters$y))), - max(ceiling(max(rep1_transcripts$y)), - ceiling(max(rep1_clusters$y)))) - grid_length <- 10 # get spatial vectors rep1_sq10_vectors_lst <- get_vectors(x= list("rep1" = rep1_transcripts), @@ -520,8 +491,7 @@ rep1_sq10_vectors_lst <- get_vectors(x= list("rep1" = rep1_transcripts), cluster_info = rep1_clusters, bin_type="square", bin_param=c(grid_length,grid_length), - test_genes = all_real_genes , - w_x=w_x, w_y=w_y) + test_genes = all_real_genes) # the created spatial vectors will be the same as from other input structure table(rep1_sq10_vectors$gene_mt[,"DST"]==rep1_sq10_vectors_lst$gene_mt[,"DST"]) # spatial vector for every cluster @@ -580,8 +550,6 @@ colnames(clusters_info) <- c("x","y") clusters_info$cell_id <- row.names(clusters_info) clusters_info$sample <- spe_sub$sample_id clusters_info <- merge(clusters_info, scran_clusters, by="cell_id") -w_x <- c(floor(min(clusters_info$x)), ceiling(max(clusters_info$x))) -w_y <- c(floor(min(clusters_info$y)), ceiling(max(clusters_info$y))) # ----------------------------------------------------------------------------- # build spatial vectors from count matrix and cluster coordinates example_vectors_cm <- get_vectors(x= spe_sub, @@ -590,8 +558,7 @@ example_vectors_cm <- get_vectors(x= spe_sub, bin_type="square", bin_param=c(5,5), test_genes = row.names(spe_sub)[1:5], - use_cm=TRUE, - w_x=w_x, w_y=w_y) + use_cm=TRUE) # spatial vector for every cluster head(example_vectors_cm$cluster_mt) @@ -605,8 +572,7 @@ example_vectors_tr <- get_vectors(x= spe_sub, bin_type="square", bin_param=c(5,5), test_genes = row.names(spe_sub)[1:5], - use_cm=FALSE, - w_x=w_x, w_y=w_y) + use_cm=FALSE) # spatial vector for every cluster # example_vectors_tr$cluster_mt @@ -631,7 +597,6 @@ example_vectors_cm <- get_vectors(x= sfe_example, bin_type="square", bin_param=c(5,5), test_genes = row.names(spe_sub)[1:3], - w_x=w_x, w_y=w_y, use_cm = TRUE) # spatial vector for every cluster @@ -647,7 +612,6 @@ example_vectors_tr <- get_vectors(x= sfe_example, bin_type="square", bin_param=c(5,5), test_genes = row.names(spe_sub)[1:3], - w_x=w_x, w_y=w_y, use_cm = FALSE) # spatial vector for every cluster @@ -677,7 +641,6 @@ example_vectors_cm <- get_vectors(x= sce_example, bin_type="square", bin_param=c(5,5), test_genes = row.names(spe_sub)[1:3], - w_x=w_x, w_y=w_y, use_cm=TRUE) # spatial vector for every cluster @@ -699,8 +662,7 @@ example_vectors_tr <- get_vectors(x= spe, sample_names = "sample01", cluster_info = clusters_info, bin_type="square", bin_param=c(5,5), - test_genes = row.names(spe_sub)[1:3], - w_x=w_x, w_y=w_y) + test_genes = row.names(spe_sub)[1:3]) # spatial vector for every cluster # example_vectors_tr$cluster_mt @@ -727,8 +689,6 @@ sce <- SingleCellExperiment(list(sample01 =seu_obj@assays$RNA$counts )) # cluster, x, y and sample clusters_info = clusters_info -w_x <- c(floor(min(clusters_info$x)), ceiling(max(clusters_info$x))) -w_y <- c(floor(min(clusters_info$y)), ceiling(max(clusters_info$y))) # ----------------------------------------------------------------------------- # build spatial vectors from count matrix and cluster coordinates # make sure the cluster information contains column names: @@ -738,7 +698,6 @@ example_vectors_cm <- get_vectors(x= sce, sample_names = "sample01", bin_type="square", bin_param=c(10,10), test_genes = test_genes, - w_x=w_x, w_y=w_y, use_cm=TRUE) # spatial vector for every cluster example_vectors_cm$cluster_mt @@ -760,8 +719,7 @@ example_vectors_tr <- get_vectors(x= spe, sample_names = "sample01", cluster_info = clusters_info, bin_type="square", bin_param=c(10,10), - test_genes = test_genes, - w_x=w_x, w_y=w_y) + test_genes = test_genes) # spatial vector for every cluster example_vectors_tr$cluster_mt @@ -834,14 +792,6 @@ cluster and cluster vector. ```{r rep1 permutation} -w_x <- c(min(floor(min(rep1_transcripts$x)), - floor(min(rep1_clusters$x))), - max(ceiling(max(rep1_transcripts$x)), - ceiling(max(rep1_clusters$x)))) -w_y <- c(min(floor(min(rep1_transcripts$y)), - floor(min(rep1_clusters$y))), - max(ceiling(max(rep1_transcripts$y)), - ceiling(max(rep1_clusters$y)))) set.seed(seed_number) perm_p <- compute_permp(x=rep1_sub, cluster_info=rep1_clusters, @@ -851,9 +801,7 @@ perm_p <- compute_permp(x=rep1_sub, test_genes= all_real_genes, correlation_method = "pearson", n_cores=1, - correction_method="BH", - w_x=w_x , - w_y=w_y) + correction_method="BH") # observed correlation for every pair of gene and cluster vector obs_corr <- get_cor(perm_p) @@ -993,7 +941,6 @@ rep1_nc_vectors <- create_genesets(x=rep1_neg,sample_names="rep1", codeword=codeword_nm), bin_type="square", bin_param=c(10, 10), - w_x=w_x, w_y=w_y, cluster_info = NULL) set.seed(seed_number) @@ -1193,18 +1140,6 @@ genes by providing additional sample information as the input for the function ### Linear modeling approach to detect marker genes ```{r tworep no background} -w_x <- c(min(floor(min(rep1_transcripts$x)), - floor(min(rep2_transcripts$x)), - floor(min(rep_clusters$x))), - max(ceiling(max(rep1_transcripts$x)), - ceiling(max(rep2_transcripts$x)), - ceiling(max(rep_clusters$x)))) -w_y <- c(min(floor(min(rep1_transcripts$y)), - floor(min(rep2_transcripts$y)), - floor(min(rep_clusters$y))), - max(ceiling(max(rep1_transcripts$y)), - ceiling(max(rep2_transcripts$y)), - ceiling(max(rep_clusters$y)))) grid_length<-10 twosample_spe<-cbind(rep1_sub, rep2_sub) # get spatial vectors @@ -1212,8 +1147,7 @@ two_rep_vectors<- get_vectors(x= twosample_spe, sample_names=c("rep1","rep2"), cluster_info = rep_clusters, bin_type="square", bin_param=c(grid_length, grid_length), - test_genes = all_real_genes , - w_x=w_x, w_y=w_y) + test_genes = all_real_genes) twosample_neg_spe<-cbind(rep1_neg, rep2_neg) @@ -1230,7 +1164,6 @@ two_rep_nc_vectors<-create_genesets(x=twosample_neg_spe, codeword=codeword_nm), bin_type="square", bin_param=c(10,10), - w_x=w_x, w_y=w_y, cluster_info = NULL) set.seed(seed_number) two_rep_lasso_with_nc<-lasso_markers(gene_mt=two_rep_vectors$gene_mt,