## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( eval = TRUE, collapse = TRUE, comment = "#>", out.width = "100%", dev = "png", dpi = 50, fig.height = 4.2, fig.width = 5.6 ) ## ----eval=FALSE--------------------------------------------------------------- # devtools::install_github("jinming-cheng/SpNeigh") ## ----loadPackage-------------------------------------------------------------- library(SpNeigh) library(ggplot2) ## ----------------------------------------------------------------------------- coords <- readRDS(system.file("extdata", "MouseBrainCoords.rds", package = "SpNeigh" )) head(coords) ## ----------------------------------------------------------------------------- rownames(coords) <- coords$cell ## ----------------------------------------------------------------------------- logNorm_expr <- readRDS(system.file("extdata", "LogNormExpr.rds", package = "SpNeigh" )) class(logNorm_expr) ## ----------------------------------------------------------------------------- new_cell_types <- c( "0" = "Meninges", "1" = "Cerebral_cortex", "2" = "White_matter", "3" = "Cerebral_cortex", "4" = "Cerebral_cortex", "5" = "Thalamus", "6" = "Cerebral_cortex", "7" = "Hippocampus", "8" = "Hippocampus", "9" = "Hippocampus", "10" = "White_matter", "11" = "Cerebral_cortex" ) ## ----------------------------------------------------------------------------- coords_sub <- subset(coords, cluster %in% c("0", "2")) coords_sub <- as.matrix(coords_sub[, c("x", "y")]) ## ----------------------------------------------------------------------------- metadata_sub <- subset(coords[, c("cell", "cluster")], cluster %in% c("0", "2")) ## ----------------------------------------------------------------------------- spe <- SpatialExperiment::SpatialExperiment( assay = list("logcounts" = logNorm_expr), colData = metadata_sub, spatialCoords = coords_sub ) spe ## ----------------------------------------------------------------------------- spe_coords <- extractCoords(spe) ## ----------------------------------------------------------------------------- spe_logcounts <- SingleCellExperiment::logcounts(spe) ## ----------------------------------------------------------------------------- library(Seurat) ## ----------------------------------------------------------------------------- seu_sp <- CreateSeuratObject( assay = "Spatial", counts = logNorm_expr, meta.data = metadata_sub ) seu_sp ## ----------------------------------------------------------------------------- LayerData(seu_sp, assay = "Spatial", layer = "data") <- logNorm_expr ## ----------------------------------------------------------------------------- cents <- CreateCentroids(coords_sub[, c("x", "y")]) fov <- CreateFOV( coords = list("centroids" = cents), type = c("centroids"), assay = "Spatial" ) seu_sp[["fov"]] <- fov ## ----------------------------------------------------------------------------- seu_sp$seurat_clusters <- seu_sp$cluster ## ----------------------------------------------------------------------------- seu_sp_coords <- extractCoords(seu_sp) ## ----------------------------------------------------------------------------- seu_sp_log_expr <- GetAssayData(seu_sp) ## ----PlotBoundary_spe, fig.width=6.6, fig.height=4.2-------------------------- plotBoundary(spe, one_cluster = "2") ## ----PlotBoundary_seu, fig.width=6.6, fig.height=4.2-------------------------- plotBoundary(seu_sp, one_cluster = "2", eps = 120) ## ----------------------------------------------------------------------------- bon_points_c2 <- getBoundary( data = coords, one_cluster = 2, eps = 120, minPts = 10 ) table(bon_points_c2$region_id) ## ----PlotWithBoundary_C2, fig.width=6.6, fig.height=4.2----------------------- plotBoundary(coords, boundary = bon_points_c2) ## ----------------------------------------------------------------------------- # plotBoundary(coords) + addBoundary(boundary = bon_points_c2) ## ----PlotRegion_C2,fig.height = 2.8, fig.width = 5.6------------------------- bon_polys_c2 <- buildBoundaryPoly(bon_points_c2) plotRegion(bon_polys_c2) ## ----------------------------------------------------------------------------- ring_regions <- getRingRegion(boundary = bon_points_c2) ## ----PlotRing_C2,fig.height = 2.8, fig.width = 5.6--------------------------- plotRegion(ring_regions) ## ----------------------------------------------------------------------------- cells_ring <- getCellsInside(data = coords, boundary = ring_regions) cells_ring ## ----PlotCellsInsideRing_C2,fig.height = 3.6, fig.width = 5.6---------------- plotCellsInside(cells_ring, point_size = 0.2) ## ----------------------------------------------------------------------------- stats_ring <- statsCellsInside(cells_ring) ## ----PlotStatsBar_Proportion_C2_Ring, fig.height = 4.2, fig.width = 6.6------ plotStatsBar(stats_ring, stat_column = "proportion") ## ----PlotStats_Donut_C2_Ring, fig.height = 3.5, fig.width = 5.6-------------- plotStatsPie(stats_ring, plot_donut = TRUE) ## ----------------------------------------------------------------------------- interaction_matrix <- computeSpatialInteractionMatrix( subset(coords, cell %in% cells_ring$cell) ) interaction_matrix ## ----Heatmap_InteractionMatrix_C2_Ring, fig.height = 4.2, fig.width = 5.6---- plotInteractionMatrix(interaction_matrix) ## ----PlotCellsInside_SplitByCluster_C2_Ring, fig.width=8.4, fig.height=4.2---- plotCellsInside(cells_ring) + facet_wrap(~cluster) + Seurat::RotatedAxis() + addBoundaryPoly(ring_regions, linewidth_boundary = 0.1) ## ----------------------------------------------------------------------------- cells_inside <- getCellsInside(data = coords, boundary = bon_points_c2) cells_inside ## ----------------------------------------------------------------------------- cells_in <- subset(cells_inside, region_id == 1 & cluster == 2)[["cell"]] cells_out <- subset(cells_ring, region_id == 1 & cluster == 2)[["cell"]] ## ----------------------------------------------------------------------------- tab <- runLimmaDE( exp_mat = logNorm_expr, min.pct = 0.25, cells_reference = cells_in, cells_target = cells_out ) ## ----------------------------------------------------------------------------- head(tab[, c("gene", "logFC", "adj.P.Val", "pct.reference", "pct.target")]) ## ----------------------------------------------------------------------------- set.seed(123) ## ----Expression_Slc17a7_c2, fig.width = 8, fig.height = 2.5------------------- plotExpression( data = coords[colnames(logNorm_expr), ], exp_mat = logNorm_expr, genes = tab$gene[1:2], sub_plot = TRUE, one_cluster = 2, shuffle = TRUE, point_size = 0.1, angle_x_label = 45 ) ## ----Expression_Slc17a7_InOutBoundary,fig.height = 2.8, fig.width = 5.6------- plotExpression( data = coords[colnames(logNorm_expr), ], exp_mat = logNorm_expr, genes = "Slc17a7", sub_plot = TRUE, shuffle = TRUE, sub_cells = c(cells_in, cells_out), point_size = 0.3 ) + addBoundaryPoly(bon_polys_c2[1, ], linewidth_boundary = 0.5) ## ----PlotExpression_Sox10,fig.height = 2.8, fig.width = 5.6------------------ plotExpression( data = coords[colnames(logNorm_expr), ], exp_mat = logNorm_expr, sub_cells = c(cells_in, cells_out), sub_plot = TRUE, genes = "Sox10", shuffle = TRUE, point_size = 0.3 ) + addBoundaryPoly(bon_polys_c2[1, ], linewidth_boundary = 0.5) ## ----------------------------------------------------------------------------- cells_c0 <- subset(coords, cluster == 0)[, "cell"] ## ----Plot_NoBoundary_C0, fig.width=8*0.7+1, fig.height=6*0.7------------------ plotBoundary(coords[cells_c0, ]) ## ----------------------------------------------------------------------------- bon_points_c0 <- getBoundary(data = coords, one_cluster = 0) bon_polys_c0 <- buildBoundaryPoly(bon_points_c0) ## ----PlotBoundary_C0, fig.width=6.6, fig.height=4.2--------------------------- plotEdge(boundary_poly = bon_polys_c0, linewidth_boundary = 0.6) ## ----------------------------------------------------------------------------- weights_bon <- computeBoundaryWeights( data = coords, cell_ids = cells_c0, boundary = bon_points_c0 ) ## ----------------------------------------------------------------------------- weights_bon[1:3] ## ----------------------------------------------------------------------------- weights_cen <- computeCentroidWeights(data = coords, cell_ids = cells_c0) ## ----fig.width=6.6, fig.height=4.2-------------------------------------------- plotWeights(data = coords, weights = weights_bon, point_size = 0.8) + labs(title = "Boundary weights") ## ----------------------------------------------------------------------------- tab_sp <- runSpatialDE( exp_mat = logNorm_expr[, cells_c0], spatial_distance = weights_bon, cell_ids = cells_c0 ) ## ----------------------------------------------------------------------------- head(tab_sp) ## ----Expression_top_genes_boundary, fig.width = 8, fig.height = 2.5----------- plotExpression( data = coords[colnames(logNorm_expr), ], exp_mat = logNorm_expr, genes = tab_sp$gene[1:2], sub_plot = TRUE, one_cluster = 0, shuffle = TRUE, point_size = 0.1, angle_x_label = 45 ) ## ----Expression_along_boundary_weights_c0, fig.width=6.4, fig.height=4-------- plotSpatialExpression( exp_mat = logNorm_expr[, cells_c0], spatial_distance = weights_bon, scale_method = "minmax", genes = tab_sp$gene[1:10], label_x = "Boundary weights" ) ## ----------------------------------------------------------------------------- SEI_bon <- computeSpatialEnrichmentIndex( exp_mat = logNorm_expr[, cells_c0], weights = weights_bon ) head(SEI_bon) ## ----Expression_top_genes_SEI_boundary, fig.width = 8, fig.height = 2.5------- plotExpression( data = coords[colnames(logNorm_expr), ], exp_mat = logNorm_expr, genes = SEI_bon$gene[1:2], sub_plot = TRUE, one_cluster = 0, shuffle = TRUE, point_size = 0.1, angle_x_label = 45 ) ## ----------------------------------------------------------------------------- sessionInfo()