Updated vignettes

vignettes/bivariate.Rmd

@@ -91,10 +91,8 @@ if (!file.exists("tissue_lowres_5a.jpeg")) {
 ```
 
 ```{r}
-sfe <- addImg(sfe, file = "tissue_lowres_5a.jpeg", sample_id = "Vis5A", 
-              image_id = "lowres", 
-              scale_fct = 1024/22208)
-sfe <- mirrorImg(sfe, sample_id = "Vis5A", image_id = "lowres")
+sfe <- addImg(sfe, imageSource = "tissue_lowres_5a.jpeg", sample_id = "Vis5A", 
+              image_id = "lowres", scale_fct = 1024/22208)
 ```
 
 ```{r}

vignettes/localc.Rmd

@@ -94,10 +94,9 @@ if (!file.exists("tissue_lowres_5a.jpeg")) {
 ```
 
 ```{r}
-sfe <- addImg(sfe, file = "tissue_lowres_5a.jpeg", sample_id = "Vis5A", 
+sfe <- addImg(sfe, imageSource = "tissue_lowres_5a.jpeg", sample_id = "Vis5A", 
               image_id = "lowres", 
               scale_fct = 1024/22208)
-sfe <- mirrorImg(sfe, sample_id = "Vis5A", image_id = "lowres")
 ```
 
 ```{r}

vignettes/multispati.Rmd

@@ -90,6 +90,8 @@ library(BiocSingular)
 library(BiocParallel)
 library(sf)
 library(patchwork)
+library(spdep)
+set.ZeroPolicyOption(TRUE)
 theme_set(theme_bw())
 ```
 
@@ -143,7 +145,7 @@ Next we compute Moran's I for QC metrics, which requires a spatial neighborhood
 system.time(
     colGraph(sfe, "knn5") <- findSpatialNeighbors(sfe, method = "knearneigh", 
                                                   dist_type = "idw", k = 5, 
-                                                  style = "W")
+                                                  style = "W", zero.policy = TRUE)
 )
 ```
 
@@ -167,14 +169,33 @@ inds2 <- st_intersects(cellSeg(sfe), st_as_sfc(st_bbox(bbox_use)),
 sfe_sub <- sfe[, inds2]
 ```
 
+Note that since SpatialFeatureExperiment v1.8, the sptial graphs are subsetted rather than reconstructed when the SFE object is subsetted because reconstruction tends to be more time consuming and the BPPARAM and BNPARAM arguments can't be stored and saved by `alabaster.sfe`. Subsetting removes some of the neighbors of cells near the boundary of this bounding box but the vast majority of cells still have all 5 nearest neighbors.
+
+```{r}
+table(card(colGraph(sfe_sub, "knn5")$neighbours))
+```
+
 ```{r}
 (mb <- moranBounds(colGraph(sfe_sub, "knn5")))
 ```
 
-So considering the bounds, the MOran's I values of the QC metrics are more like
+The lower bound is quite different from if we reconstruct the knn graph on the subset
+```{r}
+colGraph(sfe_sub, "knn5_reconst") <- findSpatialNeighbors(sfe_sub, method = "knearneigh", 
+                                                          k = 5L, dist_type = "idw",
+                                                          style = "W")
+```
+
+```{r}
+(mb2 <- moranBounds(colGraph(sfe_sub, "knn5_reconst")))
+```
+
+It would be cool to systematically investigate the effects of perturbations to the spatial neighborhood graph on Moran's I and other spatial statistics.
+
+So considering the bounds, the Moran's I values of the QC metrics are more like
 ```{r}
 setNames(colFeatureData(sfe)[c("nCounts", "nGenes", "volume"), 
-                             "moran.mc_statistic_sample01"] / mb["Imin"],
+                             "moran.mc_statistic_sample01"] / mb2["Imin"],
          features_use)
 ```
 
@@ -358,7 +379,7 @@ df_moran |>
     geom_line() +
     scale_color_manual(values = ditto_colors) +
     geom_hline(yintercept = 0, color = "gray") +
-    geom_hline(yintercept = mb, linetype = 2, color = "gray") +
+    geom_hline(yintercept = mb2, linetype = 2, color = "gray") +
     scale_y_continuous(breaks = scales::breaks_pretty()) +
     scale_x_continuous(breaks = scales::breaks_width(5)) +
     labs(y = "Moran's I", color = "Type", x = "Component")
@@ -367,7 +388,7 @@ df_moran |>
 In MULTISPATI, Moran's I is high in PC1 and PC2, but then sharply drops. Moran's I for the PC with the most negative eigenvalues is not very negative, which means the large magnitude of that eigenvalue comes from explaining more variance. However, considering the lower bound of Moran's I that is around -0.6 instead of -1, the magnitude of Moran's I for the PC with the most negative eigenvalue is not trivial.
 
 ```{r}
-min(df_moran$MULTISPATI_neg) / mb[1]
+min(df_moran$MULTISPATI_neg) / mb2[1]
 ```
 
 Non-spatial PCs are not sorted by Moran's I; PC5 has surprising large Moran's I. 

vignettes/vig1_visium_basic.Rmd

@@ -82,6 +82,7 @@ library(dplyr)
 library(reticulate)
 library(concordexR)
 library(BiocNeighbors)
+library(pheatmap)
 theme_set(theme_bw(10))
 ```
 
@@ -125,10 +126,11 @@ if (!file.exists("tissue_lowres_5a.jpeg")) {
 ```
 
 ```{r}
-sfe <- addImg(sfe, file = "tissue_lowres_5a.jpeg", sample_id = "Vis5A", 
+sfe <- addImg(sfe, imageSource = "tissue_lowres_5a.jpeg", sample_id = "Vis5A", 
               image_id = "lowres", 
               scale_fct = 1024/22208)
-sfe <- mirrorImg(sfe, sample_id = "Vis5A", image_id = "lowres")
+# Don't need to mirror for terra >= 1.7.83
+#sfe <- mirrorImg(sfe, sample_id = "Vis5A", image_id = "lowres")
 ```
 
 
@@ -313,24 +315,29 @@ plotUMAP(sfe_tissue, colour_by = "cluster")
 UMAP is often used to visualize clusters. An alternative to UMAP is [`concordex`](https://bioconductor.org/packages/devel/bioc/html/concordexR.html), which quantitatively shows the proportion of neighbors on the k nearest neighbor graph with the same cluster label. To be consistent with the default in igraph Leiden clustering, we use k = 10.
 
 ```{r}
-g <- findKNN(reducedDim(sfe_tissue, "PCA")[,1:3], k = 10)
-res <- calculateConcordex(g$index, labels = sfe_tissue$cluster, k = 10, 
-                          return.map = TRUE)
+res <- calculateConcordex(sfe_tissue, labels = sfe_tissue$cluster, 
+                          n_neighbors = 18,
+                          BLUSPARAM = SNNGraphParam(
+                              cluster.fun = "leiden",
+                              cluster.args = list(
+                                  resolution_parameter = 0.5,
+                                  objective_function = "modularity")))
 ```
 
-The cluster labels are permuted to estimate a null distribution, and the observed values can be compared to the simulated values:
+The result is a neighborhood consolidation matrix which can help finding spatial regions, see the [concordex paper](https://www.biorxiv.org/content/10.1101/2023.06.28.546949v2). The heatmap here shows the proportion of the spatial neighborhood (k nearest neighbors with k = 18 here for 2nd order neighbors on the hexagonal array) of each spot that belongs to a given gene expression based cluster; this matrix shown in this heatmap is then clustered 
 ```{r}
-plotConcordexSim(res)
+pheatmap(res)
 ```
 
-The observed value is much higher than in the simulated values, indicating good clustering. The single number is the average of all clusters. Values for different clusters can be plotted in a heatmap:
-```{r, fig.width=5, fig.height=4}
-heatConcordex(res, angle_col = 0, cluster_rows = FALSE, cluster_cols = FALSE)
+```{r}
+sfe_tissue$concordex <- attr(res, "shrs")
 ```
 
-The diagonal represents the proportion of neighbors of cells from each cluster that are from the same cluster. That the off diagonal entries are very low indicate good clustering.
+```{r}
+plotSpatialFeature(sfe_tissue, "concordex", image = "lowres", maxcell = 5e4)
+```
 
-More interesting in spatial transcriptomics, is to locate the clusters in space, and this can be done as follows:
+Also plot the clusters that did not use spatial information here
 ```{r}
 plotSpatialFeature(sfe_tissue, "cluster", colGeometryName = "spotPoly",
                    image = "lowres")

vignettes/vig2_visium.Rmd

@@ -109,10 +109,9 @@ if (!file.exists("tissue_lowres_5a.jpeg")) {
 
 The image can be added to the SFE object and plotted behind the geometries, and needs to be flipped to align to the spots because the origin is at the top left for the image but bottom left for geometries.
 ```{r}
-sfe <- addImg(sfe, file = "tissue_lowres_5a.jpeg", sample_id = "Vis5A", 
+sfe <- addImg(sfe, imageSource = "tissue_lowres_5a.jpeg", sample_id = "Vis5A", 
               image_id = "lowres", 
               scale_fct = 1024/22208)
-sfe <- mirrorImg(sfe, sample_id = "Vis5A", image_id = "lowres")
 ```
 
 # Exploratory data analysis