From 73243feaa39f6d9f9ffafb948c0f4a8ac6b3f768 Mon Sep 17 00:00:00 2001 From: Dmec5133 Date: Tue, 31 Oct 2023 11:39:19 +1100 Subject: [PATCH 1/5] CRAN Dep remove rgeos/maptools depdencies --- DESCRIPTION | 5 +- R/cellshape_metrics.R | 211 ++++++++++++++++++++++++++---------------- vignettes/CellSPA.Rmd | 22 +++-- 3 files changed, 147 insertions(+), 91 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 4e47abd..0a9722a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -25,17 +25,18 @@ Imports: BiocParallel, grDevices, stats, + lwgeom, Matrix, Metrics, proxyC, RcppParallel, + raster, alphahull, igraph, - rgeos, shotGroups, sp, + spatstat, spatstat.geom, - maptools, entropy, data.table, grr diff --git a/R/cellshape_metrics.R b/R/cellshape_metrics.R index dc7d719..4cba630 100644 --- a/R/cellshape_metrics.R +++ b/R/cellshape_metrics.R @@ -81,11 +81,6 @@ chull_to_poly2 <- function(coord){ return(sp_poly) } - - - - - # for functions requiring the bounding box we just need lines rather than sp_polygon, # you can also get the bounding box from polygon but this is more computationally expensive get_ashape_lines <- function(coord){ @@ -103,45 +98,23 @@ get_ashape_lines <- function(coord){ } } -# -# ashape_to_poly <- function(x,y){ -# points <- cbind(x,y) -# -# for (alpha_var in 1:100){ -# alpha_shape <- alphahull::ashape(points, alpha = alpha_var) -# pth <- tryCatch(ashape2poly(alpha_shape),error=function(e) e) -# -# if ( !inherits(pth, "error")){ -# -# hull_points <- points[pth,] -# poly <- sp::Polygon(hull_points) -# poly_list <- list(poly) -# -# # Create spatial polygon object -# sp_poly <- sp::SpatialPolygons(list(sp::Polygons(poly_list, ID = 1))) -# -# -# return(sp_poly) -# } -# } -# } -# + get_ashape_chull_poly <- function(x) { - polygon <- try(ashape_to_poly2(x), silent = TRUE) - if (methods::is(polygon, "try-error")) { - polygon <- list(NA, NA) + poly_res <- try(ashape_to_poly2(x), silent = TRUE) + if (methods::is(poly_res, "try-error")) { + poly_res <- list(NA, NA) } - ch_polygon <- try(chull_to_poly2(x), silent = TRUE) - if (methods::is(ch_polygon, "try-error")) { - ch_polygon <- NA + ch_poly_res <- try(chull_to_poly2(x), silent = TRUE) + if (methods::is(ch_poly_res, "try-error")) { + ch_poly_res <- NA } - return(list(ashape_poly = polygon[[1]], - chull_poly = ch_polygon, - ashape_alpha = polygon[[2]])) + return(list(ashape_poly = poly_res[[1]], + chull_poly = ch_poly_res, + ashape_alpha = poly_res[[2]])) } elongation <- function(x) { @@ -164,35 +137,27 @@ elongation <- function(x) { } - eccentricity <- function(x){ if (methods::is(x, "matrix") | methods::is(x, "data.frame")) { # Calculate the alpha shape - lines <- get_ashape_lines(x) + sp_poly <- ashape_to_poly2(x)[[1]] + sf_poly <- st_as_sf(sp_poly) } if (methods::is(x, "SpatialPolygons")) { - lines <- x@polygons[[1]]@Polygons[[1]]@coords + sf_poly <- st_as_sf(x) } + major_axis <- find_major_axis(sf_poly) + minor_axis <- find_minor_axis(sf_poly,major_axis$major_axis_endpoints) - #nice function to extract bbox that is not a flat square - bb <- shotGroups::getMinBBox(lines) - - - #this is super similar to elongation but basically we take the max and min, for many cells we get the same result as elongation - major_axis <- max(bb$height, bb$width) - minor_axis <- min(bb$height, bb$width) - - return(minor_axis/major_axis) + return(as.numeric(minor_axis$minor_axis_length/major_axis$major_axis_length)) } - - sphericity <- function(x){ if (methods::is(x, "matrix") | methods::is(x, "data.frame")) { @@ -204,41 +169,34 @@ sphericity <- function(x){ res <- x } - - #spatstat needs a different type of spatial object - pol <- maptools::as.owin.SpatialPolygons(res) - #if (is(pol, "try-error")) { return(NULL) } + bbox <- res@bbox + pol <- spatstat.geom::owin(xrange = c(bbox["x", "min"], bbox["x", "max"]), yrange = c(bbox["y", "min"], bbox["y", "max"])) #calculate outer circle outer_radius <- spatstat.geom::boundingradius(pol)[[1]] #get inner circle inner_radius <- spatstat.geom::incircle(pol) + #perform calculation of radius of inner circle divided by radius of outer circle return(inner_radius$r/outer_radius) } - -solidity <- function(x, chull_polygon = NULL){ - +solidity <- function(x, chull_polygon = NULL) { if (methods::is(x, "matrix") | methods::is(x, "data.frame")) { - polygon <- ashape_to_poly2(x)[[1]] + poly_res <- ashape_to_poly2(x)[[1]] ch_polygon <- chull_to_poly2(x) } if (methods::is(x, "SpatialPolygons") & methods::is(chull_polygon, "SpatialPolygons")) { - polygon <- x + poly_res <- x ch_polygon <- chull_polygon } + area_p <- as.numeric(raster::area(poly_res)) + ch_area <- as.numeric(raster::area(ch_polygon)) - area <- as.numeric(rgeos::gArea(polygon)) - - ch_area <- as.numeric(rgeos::gArea(ch_polygon)) - - return(area/ch_area) - + return(area_p / ch_area) } - convexity <- function(x, chull_polygon = NULL){ @@ -255,16 +213,15 @@ convexity <- function(x, chull_polygon = NULL){ # Calculate area of polygon object - perimeter <- as.numeric(rgeos::gLength(polygon)) + perimeter <- as.numeric(lwgeom::st_perimeter_2d(st_as_sf(polygon))) # Calculate perimeter of convex hull object - ch_perimeter <- as.numeric(rgeos::gLength(ch_polygon)) + ch_perimeter <- as.numeric(lwgeom::st_perimeter_2d(st_as_sf(ch_polygon))) return(ch_perimeter/perimeter) } - circularity <- function(x, chull_polygon = NULL){ #This one kinda sucks for us because we dont really have a ground truth polygon, its basically just alpha convex hull(tighter polygon) # versus a normal convex hull @@ -282,40 +239,132 @@ circularity <- function(x, chull_polygon = NULL){ # Calculate area of polygon object - area <- as.numeric(rgeos::gArea(polygon)) + area <- as.numeric(raster::area(polygon)) # Calculate perimeter of convex hull object - convex_perimeter <- as.numeric(rgeos::gLength(ch_polygon)) + ch_perimeter <- as.numeric(lwgeom::st_perimeter_2d(st_as_sf(ch_polygon))) #circularity calculation - return((4*pi*area)/ ((convex_perimeter)^2)) + return((4*pi*area)/ ((ch_perimeter)^2)) } - - compactness <- function(x) { if (methods::is(x, "matrix") | methods::is(x, "data.frame")) { # Calculate the alpha shape - polygon <- ashape_to_poly2(x)[[1]] + poly_res <- ashape_to_poly2(x)[[1]] } if (methods::is(x, "SpatialPolygons")) { - polygon <- x + poly_res <- x } # Calculate area of polygon object - area <- as.numeric(rgeos::gArea(polygon)) + area <- as.numeric(raster::area(poly_res)) + + # Calculate perimeter of hull object + perimeter <- as.numeric(lwgeom::st_perimeter_2d(st_as_sf(poly_res))) - # Calculate perimeter of convex hull object - perimeter <- as.numeric(rgeos::gLength(polygon)) - #Compactness calculation return((4*pi*area)/ ((perimeter)^2)) } +# support functions for eccentricity to find major and minor axis, this function essentially breaks a line into #num samples number of points to find the longest line through the polygon that is within the polygon +# the more samples the more accurate but also computational expensive, angle threshhold allows for slight leeway in minor not being perfectly perpendicular to major, again a computationally expensive trade off for better accuracy +find_major_axis <- function(poly_res, num_samples = 25) { + + #we get the boundary poly_res and then create lots of samples along its perimeter + boundary <- st_boundary(poly_res) + boundary_length <- st_length(boundary)[[1]] + segment_length <- boundary_length / num_samples + boundary_points <- st_coordinates(st_segmentize(boundary, segment_length)) + + max_distance <- -Inf + max_indices <- c(0, 0) + num_boundary_points <- nrow(boundary_points) + #loop through sampled x,y boundary points to find lines, keep biggest + for (i in 1:(num_boundary_points - 1)) { + for (j in (i + 1):num_boundary_points) { + line_segment <- st_linestring(matrix(c(boundary_points[i,], boundary_points[j,]), nrow = 2, byrow = TRUE)) + contains_result <- st_contains(poly_res, st_sfc(line_segment)) + + if (length(contains_result[[1]]) > 0) { # line segment is contained inside the poly_res or part of its perimeter + current_distance <- sqrt((boundary_points[i, 1] - boundary_points[j, 1])^2 + + (boundary_points[i, 2] - boundary_points[j, 2])^2) + + if (current_distance > max_distance) { + max_distance <- current_distance + max_indices <- c(i, j) + } + } + } + } + + major_axis_endpoints <- boundary_points[max_indices, ] + major_axis_length <- max_distance + #here we return the endpoints and length, ultimately eccentricity just is minor axis/major axis but this way we can do plots + result <- list( + major_axis_endpoints = major_axis_endpoints, + major_axis_length = major_axis_length + ) + + return(result) +} +find_minor_axis <- function(poly_res, major_axis_endpoints, num_samples = 25, angle_threshold = 2) { + boundary <- st_boundary(poly_res) + boundary_length <- st_length(boundary)[[1]] + segment_length <- boundary_length / num_samples + boundary_points <- st_coordinates(st_segmentize(boundary, segment_length)) + + major_axis_vector <- major_axis_endpoints[2,] - major_axis_endpoints[1,] + major_axis_angle <- atan2(major_axis_vector[2], major_axis_vector[1]) + major_axis_line <- st_linestring(major_axis_endpoints) + + max_distance <- -Inf + max_indices <- c(0, 0) + num_boundary_points <- nrow(boundary_points) + + for (i in 1:(num_boundary_points - 1)) { + for (j in (i + 1):num_boundary_points) { + line_segment <- st_linestring(matrix(c(boundary_points[i,], boundary_points[j,]), nrow = 2, byrow = TRUE)) + contains_result <- st_contains(poly_res, st_sfc(line_segment)) + + if (length(contains_result[[1]]) > 0) { # line segment is contained inside the poly_res or part of its perimeter + # Check if the line segment intersects the major axis + intersects_result <- st_intersects(major_axis_line, st_sfc(line_segment)) + + if (length(intersects_result[[1]]) > 0) { + line_segment_vector <- boundary_points[j,] - boundary_points[i,] + line_segment_angle <- atan2(line_segment_vector[2], line_segment_vector[1]) + + angle_difference <- abs(major_axis_angle - line_segment_angle) * (180 / pi) + angle_difference <- min(angle_difference, 360 - angle_difference) + + if (angle_difference >= (90 - angle_threshold) && angle_difference <= (90 + angle_threshold)) { + current_distance <- sqrt((boundary_points[i, 1] - boundary_points[j, 1])^2 + + (boundary_points[i, 2] - boundary_points[j, 2])^2) + + if (current_distance > max_distance) { + max_distance <- current_distance + max_indices <- c(i, j) + } + } + } + } + } + } + + minor_axis_endpoints <- boundary_points[max_indices, ] + minor_axis_length <- max_distance + result <- list( + minor_axis_endpoints = minor_axis_endpoints, + minor_axis_length = minor_axis_length + ) + + return(result) +} cell_shape_fn <- function(x, method) { diff --git a/vignettes/CellSPA.Rmd b/vignettes/CellSPA.Rmd index d8545f8..8dc4f1b 100644 --- a/vignettes/CellSPA.Rmd +++ b/vignettes/CellSPA.Rmd @@ -27,30 +27,35 @@ library(BiocStyle) -```{r} +```{r, include=FALSE} + library(CellSPA) library(ggplot2) library(ggthemes) library(scater) library(SingleCellExperiment) library(SpatialExperiment) +#Not seemingly starting +library(sf) +library(spatstat) +library(lwgeom) theme_set(theme_bw()) ``` # Read BIDCell output -```{r} -data_dir <- system.file("extdata/BIDCell_csv_output", package = "CellSPA") -data_dir -tiff_path <- system.file("extdata/BIDCell_output_subset.tif", package = "CellSPA") +```{r, cache=TRUE} + +# data_dir <- system.file("extdata/BIDCell_cells_breast1/", package = "CellSPA") +# tiff_path <- system.file("extdata/BIDCell_cells_breast1.tif", package = "CellSPA") spe <- readBIDCell(data_dir, tiff_path = tiff_path, method_name = "BIDCell") ``` -```{r} +```{r, cache=TRUE} spe <- processingSPE(spe, qc_range = list(total_transciprts = c(20, 2000), total_genes = c(20, Inf))) @@ -59,9 +64,9 @@ spe <- processingSPE(spe, ```{r} -# subset a set of cells for illustration +#subset a set of cells for illustration spe <- CellSPA::subset(spe, 1:500) -spe + ``` @@ -69,6 +74,7 @@ spe ```{r} + spe <- generatePolygon(spe) spe <- calBaselineAllMetrics(spe, verbose = TRUE) head(rowData(spe)) From 1d7d673b51befe50a1485568152aaec449c9bdf1 Mon Sep 17 00:00:00 2001 From: Dmec5133 Date: Tue, 31 Oct 2023 12:07:50 +1100 Subject: [PATCH 2/5] TestRun --- vignettes/CellSPA.Rmd | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/vignettes/CellSPA.Rmd b/vignettes/CellSPA.Rmd index 8dc4f1b..934efdd 100644 --- a/vignettes/CellSPA.Rmd +++ b/vignettes/CellSPA.Rmd @@ -35,7 +35,7 @@ library(ggthemes) library(scater) library(SingleCellExperiment) library(SpatialExperiment) -#Not seemingly starting +#Why do I need to call these in the rmd,the dependent functions fail if I dont..? library(sf) library(spatstat) library(lwgeom) @@ -46,9 +46,11 @@ theme_set(theme_bw()) # Read BIDCell output ```{r, cache=TRUE} - +source("C:\\Users\\dmech\\OneDrive - The University of Sydney (Students)\\PhD\\Metrics\\CellSPA\\R\\cellshape_metrics.R") # data_dir <- system.file("extdata/BIDCell_cells_breast1/", package = "CellSPA") # tiff_path <- system.file("extdata/BIDCell_cells_breast1.tif", package = "CellSPA") +data_dir <- "data/BIDCell_cells_breast1" +tiff_path <- "data/BidCell_cells_breast1.tif" spe <- readBIDCell(data_dir, tiff_path = tiff_path, method_name = "BIDCell") @@ -65,7 +67,7 @@ spe <- processingSPE(spe, ```{r} #subset a set of cells for illustration -spe <- CellSPA::subset(spe, 1:500) +spe <- CellSPA::subset(spe, 1:50) ``` From 2f243888b02149e564014c68b9bdf44d10a48440 Mon Sep 17 00:00:00 2001 From: Dmec5133 Date: Tue, 31 Oct 2023 13:18:44 +1100 Subject: [PATCH 3/5] devtools check --- DESCRIPTION | 1 + R/cellshape_metrics.R | 36 ++++++++++++++++++------------------ 2 files changed, 19 insertions(+), 18 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 0a9722a..c07868d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -35,6 +35,7 @@ Imports: igraph, shotGroups, sp, + sf, spatstat, spatstat.geom, entropy, diff --git a/R/cellshape_metrics.R b/R/cellshape_metrics.R index 4cba630..aafa2e1 100644 --- a/R/cellshape_metrics.R +++ b/R/cellshape_metrics.R @@ -144,11 +144,11 @@ eccentricity <- function(x){ if (methods::is(x, "matrix") | methods::is(x, "data.frame")) { # Calculate the alpha shape sp_poly <- ashape_to_poly2(x)[[1]] - sf_poly <- st_as_sf(sp_poly) + sf_poly <- sf::st_as_sf(sp_poly) } if (methods::is(x, "SpatialPolygons")) { - sf_poly <- st_as_sf(x) + sf_poly <- sf::st_as_sf(x) } major_axis <- find_major_axis(sf_poly) @@ -213,10 +213,10 @@ convexity <- function(x, chull_polygon = NULL){ # Calculate area of polygon object - perimeter <- as.numeric(lwgeom::st_perimeter_2d(st_as_sf(polygon))) + perimeter <- as.numeric(lwgeom::st_perimeter_2d(sf::st_as_sf(polygon))) # Calculate perimeter of convex hull object - ch_perimeter <- as.numeric(lwgeom::st_perimeter_2d(st_as_sf(ch_polygon))) + ch_perimeter <- as.numeric(lwgeom::st_perimeter_2d(sf::st_as_sf(ch_polygon))) return(ch_perimeter/perimeter) @@ -242,7 +242,7 @@ circularity <- function(x, chull_polygon = NULL){ area <- as.numeric(raster::area(polygon)) # Calculate perimeter of convex hull object - ch_perimeter <- as.numeric(lwgeom::st_perimeter_2d(st_as_sf(ch_polygon))) + ch_perimeter <- as.numeric(lwgeom::st_perimeter_2d(sf::st_as_sf(ch_polygon))) #circularity calculation return((4*pi*area)/ ((ch_perimeter)^2)) @@ -264,7 +264,7 @@ compactness <- function(x) { area <- as.numeric(raster::area(poly_res)) # Calculate perimeter of hull object - perimeter <- as.numeric(lwgeom::st_perimeter_2d(st_as_sf(poly_res))) + perimeter <- as.numeric(lwgeom::st_perimeter_2d(sf::st_as_sf(poly_res))) return((4*pi*area)/ ((perimeter)^2)) @@ -275,10 +275,10 @@ compactness <- function(x) { find_major_axis <- function(poly_res, num_samples = 25) { #we get the boundary poly_res and then create lots of samples along its perimeter - boundary <- st_boundary(poly_res) - boundary_length <- st_length(boundary)[[1]] + boundary <- sf::st_boundary(poly_res) + boundary_length <- sf::st_length(boundary)[[1]] segment_length <- boundary_length / num_samples - boundary_points <- st_coordinates(st_segmentize(boundary, segment_length)) + boundary_points <- sf::st_coordinates(sf::st_segmentize(boundary, segment_length)) max_distance <- -Inf max_indices <- c(0, 0) @@ -286,8 +286,8 @@ find_major_axis <- function(poly_res, num_samples = 25) { #loop through sampled x,y boundary points to find lines, keep biggest for (i in 1:(num_boundary_points - 1)) { for (j in (i + 1):num_boundary_points) { - line_segment <- st_linestring(matrix(c(boundary_points[i,], boundary_points[j,]), nrow = 2, byrow = TRUE)) - contains_result <- st_contains(poly_res, st_sfc(line_segment)) + line_segment <- sf::st_linestring(matrix(c(boundary_points[i,], boundary_points[j,]), nrow = 2, byrow = TRUE)) + contains_result <- sf::st_contains(poly_res, sf::st_sfc(line_segment)) if (length(contains_result[[1]]) > 0) { # line segment is contained inside the poly_res or part of its perimeter current_distance <- sqrt((boundary_points[i, 1] - boundary_points[j, 1])^2 + @@ -312,14 +312,14 @@ find_major_axis <- function(poly_res, num_samples = 25) { return(result) } find_minor_axis <- function(poly_res, major_axis_endpoints, num_samples = 25, angle_threshold = 2) { - boundary <- st_boundary(poly_res) - boundary_length <- st_length(boundary)[[1]] + boundary <- sf::st_boundary(poly_res) + boundary_length <- sf::st_length(boundary)[[1]] segment_length <- boundary_length / num_samples - boundary_points <- st_coordinates(st_segmentize(boundary, segment_length)) + boundary_points <- sf::st_coordinates(sf::st_segmentize(boundary, segment_length)) major_axis_vector <- major_axis_endpoints[2,] - major_axis_endpoints[1,] major_axis_angle <- atan2(major_axis_vector[2], major_axis_vector[1]) - major_axis_line <- st_linestring(major_axis_endpoints) + major_axis_line <- sf::st_linestring(major_axis_endpoints) max_distance <- -Inf max_indices <- c(0, 0) @@ -327,12 +327,12 @@ find_minor_axis <- function(poly_res, major_axis_endpoints, num_samples = 25, an for (i in 1:(num_boundary_points - 1)) { for (j in (i + 1):num_boundary_points) { - line_segment <- st_linestring(matrix(c(boundary_points[i,], boundary_points[j,]), nrow = 2, byrow = TRUE)) - contains_result <- st_contains(poly_res, st_sfc(line_segment)) + line_segment <- sf::st_linestring(matrix(c(boundary_points[i,], boundary_points[j,]), nrow = 2, byrow = TRUE)) + contains_result <- sf::st_contains(poly_res, sf::st_sfc(line_segment)) if (length(contains_result[[1]]) > 0) { # line segment is contained inside the poly_res or part of its perimeter # Check if the line segment intersects the major axis - intersects_result <- st_intersects(major_axis_line, st_sfc(line_segment)) + intersects_result <- sf::st_intersects(major_axis_line, sf::st_sfc(line_segment)) if (length(intersects_result[[1]]) > 0) { line_segment_vector <- boundary_points[j,] - boundary_points[i,] From 6c0b94275749d3944903410ee2b66740cfb60fa4 Mon Sep 17 00:00:00 2001 From: Dmec5133 Date: Tue, 31 Oct 2023 13:42:16 +1100 Subject: [PATCH 4/5] devtools check --- vignettes/CellSPA.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/CellSPA.Rmd b/vignettes/CellSPA.Rmd index 934efdd..e412805 100644 --- a/vignettes/CellSPA.Rmd +++ b/vignettes/CellSPA.Rmd @@ -46,7 +46,7 @@ theme_set(theme_bw()) # Read BIDCell output ```{r, cache=TRUE} -source("C:\\Users\\dmech\\OneDrive - The University of Sydney (Students)\\PhD\\Metrics\\CellSPA\\R\\cellshape_metrics.R") +#source("C:\\Users\\dmech\\OneDrive - The University of Sydney (Students)\\PhD\\Metrics\\CellSPA\\R\\cellshape_metrics.R") # data_dir <- system.file("extdata/BIDCell_cells_breast1/", package = "CellSPA") # tiff_path <- system.file("extdata/BIDCell_cells_breast1.tif", package = "CellSPA") data_dir <- "data/BIDCell_cells_breast1" From d3cc7e7e774d22321d69f0b7e6b61f49706c8c64 Mon Sep 17 00:00:00 2001 From: Nick Date: Wed, 1 Nov 2023 15:23:46 +1100 Subject: [PATCH 5/5] remove spatstat doep --- DESCRIPTION | 1 - 1 file changed, 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index c07868d..fe47eae 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -36,7 +36,6 @@ Imports: shotGroups, sp, sf, - spatstat, spatstat.geom, entropy, data.table,