Introduction

BioGSP provides Spectral Graph Wavelet Transform (SGWT) analysis for spatial biological data. This package enables multi-scale analysis of spatial patterns using graph signal processing techniques.

Key Capabilities

Object-oriented SGWT workflow with clear separation of concerns
Multi-scale spatial analysis using spectral graph wavelets
Flexible data input with custom column naming
Energy-normalized similarity analysis with comprehensive metrics
Multiple wavelet kernels (Mexican Hat, Meyer, Heat)
Comprehensive visualization tools
Detailed inverse transform results (low-pass, band-pass approximations)

New Workflow Overview

The new workflow consists of five main steps:

initSGWT(): Initialize SGWT object with data and parameters
checkKband(): Analyze k-band limited property of signals (optional)
runSpecGraph(): Build spectral graph structure
runSGWT(): Perform forward and inverse SGWT transforms
runSGCC(): Calculate energy-normalized weighted similarity

Setup

# Load required packages
library(BioGSP)
library(ggplot2)
library(patchwork)
library(viridis)

# Set random seed for reproducibility
set.seed(123)

New SGWT Workflow Demonstration

Generate Example Pattern

# Create a spatial pattern with concentric circles
demo_pattern <- simulate_multiscale(
  grid_size = 20,        # Moderate size for fast computation
  n_centers = 1,         # Single center pattern
  Ra_seq = 5,            # Inner circle radius
  n_steps = 1,           # Single step (one pattern)
  outer_start = 10,      # Outer ring radius
  seed = 123
)[[1]]

## Generating 1 patterns...
##   |                                                                              |                                                                      |   0%  |                                                                              |======================================================================| 100%
## 
## Multicenter multiscale simulation completed!

# Display pattern structure
cat("Generated pattern dimensions:", nrow(demo_pattern), "x", ncol(demo_pattern), "\n")

## Generated pattern dimensions: 400 x 4

cat("Column names:", paste(colnames(demo_pattern), collapse = ", "), "\n")

## Column names: X, Y, signal_1, signal_2

cat("Signals available:", sum(demo_pattern$signal_1), "inner pixels,", sum(demo_pattern$signal_2), "outer pixels\n")

## Signals available: 81 inner pixels, 234 outer pixels

Visualize Input Pattern

# Single combined visualization (one plot for this section)
demo_pattern$pattern_type <- ifelse(
  demo_pattern$signal_1 == 1,
  "Inner Circle",
  ifelse(demo_pattern$signal_2 == 1, "Outer Ring", "Background")
)

p_input <- ggplot(demo_pattern, aes(X, Y, fill = pattern_type)) +
  geom_tile() +
  scale_fill_manual(
    values = c("Background" = "white", "Inner Circle" = "#e31a1c", "Outer Ring" = "#1f78b4"),
    name = "Pattern"
  ) +
  coord_fixed() +
  theme_void() +
  ggtitle("Input Pattern (Signal 1 + Signal 2)")

print(p_input)

Step 1: Initialize SGWT Object

# Initialize SGWT object with custom column names
SG <- initSGWT(
  data.in = demo_pattern,
  x_col = "X",           # Custom X coordinate column name
  y_col = "Y",           # Custom Y coordinate column name
  signals = c("signal_1", "signal_2"),  # Analyze both signals
  J = 2,                 # Number of wavelet scales
  scaling_factor = 1,    # Scale progression factor
  kernel_type = "mexican_hat"
)

# Display initialized object
print(SG)

## SGWT Object
## ===========
## Data:
##   Dimensions: 400 x 5 
##   Coordinates: X , Y 
##   Signals: signal_1, signal_2 
## 
## Parameters:
##   k (neighbors): mexican_hat 
##   J (scales): 2 
##   Kernel type: mexican_hat 
##   Laplacian type: 
## 
## Status:
##   Graph computed: FALSE 
##   Forward computed: FALSE 
##   Inverse computed: FALSE

Step 2: Build Spectral Graph

# Build spectral graph structure
SG <- runSpecGraph(SG, k = 12, laplacian_type = "normalized", length_eigenvalue = 200, verbose = TRUE)

## Building graph from spatial coordinates...
## Computing Laplacian and eigendecomposition...
## Auto-generated scales: 1.0188, 1.0188 
## Graph construction completed.

# Check updated object
cat("Graph construction completed!\n")

## Graph construction completed!

cat("Adjacency matrix dimensions:", dim(SG$Graph$adjacency_matrix), "\n")

## Adjacency matrix dimensions: 400 400

cat("Number of eigenvalues computed:", length(SG$Graph$eigenvalues), "\n")

## Number of eigenvalues computed: 200

Visualize Fourier Modes (Eigenvectors)

# Visualize Fourier modes (eigenvectors) - 5 low and 5 high frequency modes
fourier_modes <- plot_FM(SG, mode_type = "both", n_modes = 2, ncol = 2)

cat("Fourier Modes Visualization:\n")

## Fourier Modes Visualization:

cat("- Low-frequency modes (2-6): Smooth spatial patterns, excluding DC component\n")

## - Low-frequency modes (2-6): Smooth spatial patterns, excluding DC component

cat("- High-frequency modes (96-100): Fine-detailed spatial patterns\n")

## - High-frequency modes (96-100): Fine-detailed spatial patterns

cat("- Each mode shows its eigenvalue (λ) indicating frequency content\n")

## - Each mode shows its eigenvalue (λ) indicating frequency content

Step 3: Run SGWT Analysis

# Perform SGWT forward and inverse transforms
SG <- runSGWT(SG, verbose = TRUE)

## Performing SGWT analysis for 2 signals...
## Using batch processing for efficiency...
## SGWT analysis completed.

# Display final object with all results
print(SG)

## SGWT Object
## ===========
## Data:
##   Dimensions: 400 x 5 
##   Coordinates: X , Y 
##   Signals: signal_1, signal_2 
## 
## Parameters:
##   k (neighbors): mexican_hat 
##   J (scales): 2 
##   Kernel type: mexican_hat 
##   Laplacian type: 
##   Scales: 1.0188, 1.0188 
## 
## Status:
##   Graph computed: TRUE 
##   Forward computed: TRUE 
##   Inverse computed: TRUE 
## 
## Reconstruction Errors:
##    signal_1 : 0.129674 
##    signal_2 : 0.234629

Visualize SGWT Decomposition

# Plot SGWT decomposition results for signal_1
decomp_plots <- plot_sgwt_decomposition(
  SG = SG,
  signal_name = "signal_1",
  plot_scales = 1,       # One plot for this section
  ncol = 1
)

print(decomp_plots)

## TableGrob (4 x 1) "arrange": 4 grobs
##               z     cells    name           grob
## original      1 (1-1,1-1) arrange gtable[layout]
## scaling       2 (2-2,1-1) arrange gtable[layout]
## wavelet_1     3 (3-3,1-1) arrange gtable[layout]
## reconstructed 4 (4-4,1-1) arrange gtable[layout]

Energy Analysis

# Analyze energy distribution across scales for signal_1
energy_analysis <- sgwt_energy_analysis(SG, "signal_1")
print(energy_analysis)

##       scale    energy energy_ratio scale_value   signal
## 1  low_pass 57.784192   0.95227414    1.018772 signal_1
## 2 wavelet_1  1.448008   0.02386293    1.018772 signal_1
## 3 wavelet_2  1.448008   0.02386293    1.018772 signal_1

# Create energy distribution plot
energy_plot <- ggplot(energy_analysis, aes(x = scale, y = energy_ratio)) +
  geom_bar(stat = "identity", fill = "steelblue", alpha = 0.7) +
  geom_text(aes(label = paste0(round(energy_ratio * 100, 1), "%")), 
            vjust = -0.5, size = 3) +
  labs(title = "Energy Distribution Across SGWT Scales (Signal 1)",
       x = "Scale", y = "Energy Ratio") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

print(energy_plot)

Advanced Features Demonstration

Different Wavelet Kernels

# Compare different kernel types
kernels <- c("mexican_hat", "meyer", "heat")
kernel_results <- list()

for (kn in kernels) {
  # Initialize with different kernel
  SG_temp <- initSGWT(
    data.in = demo_pattern,
    x_col = "X", y_col = "Y",
    signals = "signal_1",
    J = 3,                    # Reduced for faster computation
    kernel_type = kn
  )
  
  # Run full workflow
  SG_temp <- runSpecGraph(SG_temp, k = 12, laplacian_type = "normalized", length_eigenvalue = 100,verbose = FALSE)
  SG_temp <- runSGWT(SG_temp, verbose = FALSE)
  
  kernel_results[[kn]] <- SG_temp
}

# Compare reconstruction errors
comparison_df <- data.frame(
  Kernel = kernels,
  RMSE = sapply(kernel_results, function(x) x$Inverse$signal_1$reconstruction_error),
  stringsAsFactors = FALSE
)

print("Kernel Performance Comparison:")

## [1] "Kernel Performance Comparison:"

print(comparison_df)

##                  Kernel      RMSE
## mexican_hat mexican_hat 0.2113782
## meyer             meyer 0.2019181
## heat               heat 0.1928501

# Plot comparison
ggplot(comparison_df, aes(x = Kernel, y = RMSE, fill = Kernel)) +
  geom_bar(stat = "identity", alpha = 0.7) +
  geom_text(aes(label = round(RMSE, 6)), vjust = -0.5) +
  scale_fill_viridis_d() +
  labs(title = "SGWT Reconstruction Error by Kernel Type",
       x = "Kernel Type", y = "RMSE") +
  theme_minimal()

Step 4: Pattern Similarity Analysis

# Calculate similarity between signal_1 and signal_2 in same object
similarity_within <- runSGCC("signal_1", "signal_2", SG = SG, return_parts = TRUE)

cat("Pattern Similarity Analysis (within same graph):\n")

## Pattern Similarity Analysis (within same graph):

cat(sprintf("Overall similarity: %.4f\n", similarity_within$S))

## Overall similarity: -0.5930

cat(sprintf("Low-frequency similarity: %.4f\n", similarity_within$c_low))

## Low-frequency similarity: -0.5996

cat(sprintf("Non-low-frequency similarity: %.4f\n", similarity_within$c_nonlow))

## Non-low-frequency similarity: -0.5040

cat(sprintf("Energy weights - Low: %.3f, Non-low: %.3f\n", 
           similarity_within$w_low, similarity_within$w_NL))

## Energy weights - Low: 0.932, Non-low: 0.068

# Generate a second pattern for cross-comparison
pattern_2 <- simulate_multiscale(
  grid_size = 40,
  n_centers = 1,
  Ra_seq = 8,            # Different inner radius
  n_steps = 1,           # Single step
  outer_start = 15,      # Different outer radius
  seed = 456
)[[1]]

## Generating 1 patterns...
##   |                                                                              |                                                                      |   0%  |                                                                              |======================================================================| 100%
## 
## Multicenter multiscale simulation completed!

# Create second SGWT object
SG2 <- initSGWT(pattern_2, x_col = "X", y_col = "Y", signals = "signal_1", 
                J = 4)
SG2 <- runSpecGraph(SG2, k = 12, laplacian_type = "normalized",length_eigenvalue = 100, verbose = FALSE)
SG2 <- runSGWT(SG2, verbose = FALSE)

# Note: Cross-object similarity comparison removed due to different eigenvalue counts
# between SG (200 eigenvalues) and SG2 (100 eigenvalues) causing dimension mismatch

# Visualize both patterns for comparison (single faceted plot for this section)
pattern_compare <- rbind(
  data.frame(
    X = demo_pattern$X,
    Y = demo_pattern$Y,
    signal_1 = demo_pattern$signal_1,
    Pattern = "Pattern A (R_in=5, R_out=10)"
  ),
  data.frame(
    X = pattern_2$X,
    Y = pattern_2$Y,
    signal_1 = pattern_2$signal_1,
    Pattern = "Pattern B (R_in=8, R_out=15)"
  )
)

pattern_compare$Signal <- ifelse(pattern_compare$signal_1 == 1, "Signal", "Background")

ggplot(pattern_compare, aes(X, Y, fill = Signal)) +
  geom_tile() +
  scale_fill_manual(values = c("Background" = "white", "Signal" = "#377eb8")) +
  coord_fixed() +
  facet_wrap(~ Pattern) +
  theme_void() +
  theme(legend.position = "none") +
  ggtitle("Pattern Comparison (signal_1)")

Low-frequency Only Analysis

# Compare using only low-frequency components
similarity_low <- runSGCC("signal_1", "signal_2", SG = SG, low_only = TRUE, return_parts = TRUE)

cat("Low-frequency Only Similarity Analysis:\n")

## Low-frequency Only Similarity Analysis:

cat(sprintf("Low-frequency similarity: %.4f\n", similarity_low$c_low))

## Low-frequency similarity: -0.5996

cat(sprintf("Overall score (same as low-freq): %.4f\n", similarity_low$S))

## Overall score (same as low-freq): -0.5996

cat("Note: Non-low-frequency components are ignored (c_nonlow = NA)\n")

## Note: Non-low-frequency components are ignored (c_nonlow = NA)

# Compare with full analysis
cat("\nComparison:\n")

## 
## Comparison:

cat(sprintf("Full analysis similarity: %.4f\n", similarity_within$S))

## Full analysis similarity: -0.5930

cat(sprintf("Low-only similarity: %.4f\n", similarity_low$S))

## Low-only similarity: -0.5996

cat(sprintf("Difference: %.4f\n", abs(similarity_within$S - similarity_low$S)))

## Difference: 0.0065

Multiple Signal Analysis

# Analyze energy distribution for both signals
energy_signal1 <- sgwt_energy_analysis(SG, "signal_1")
energy_signal2 <- sgwt_energy_analysis(SG, "signal_2")

# Combine for comparison
energy_comparison <- rbind(energy_signal1, energy_signal2)

# Plot comparison
ggplot(energy_comparison, aes(x = scale, y = energy_ratio, fill = signal)) +
  geom_bar(stat = "identity", position = "dodge", alpha = 0.7) +
  geom_text(aes(label = paste0(round(energy_ratio * 100, 1), "%")), 
            position = position_dodge(width = 0.9), vjust = -0.5, size = 3) +
  scale_fill_viridis_d() +
  labs(title = "Energy Distribution Comparison: Signal 1 vs Signal 2",
       x = "Scale", y = "Energy Ratio", fill = "Signal") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Similarity Space Visualization with Multiple Patterns

# Generate 9 paired patterns using simulate_multiscale
cat("Generating 9 paired patterns for similarity analysis...\n")

## Generating 9 paired patterns for similarity analysis...

patterns_9 <- simulate_multiscale(
  grid_size = 40,
  n_centers = 1,
  Ra_seq = c(3, 6, 9),     # Inner circle radii
  n_steps = 3,             # 3 shrinkage steps
  outer_start = 20,        # Fixed starting outer radius
  seed = 123
)

## Generating 9 patterns...
##   |                                                                              |                                                                      |   0%  |                                                                              |========                                                              |  11%  |                                                                              |================                                                      |  22%  |                                                                              |=======================                                               |  33%  |                                                                              |===============================                                       |  44%  |                                                                              |=======================================                               |  56%  |                                                                              |===============================================                       |  67%  |                                                                              |======================================================                |  78%  |                                                                              |==============================================================        |  89%  |                                                                              |======================================================================| 100%
## 
## Multicenter multiscale simulation completed!

cat("Generated", length(patterns_9), "patterns\n")

## Generated 9 patterns

cat("Pattern names:", names(patterns_9), "\n")

## Pattern names: Ra_3_Step_1 Ra_3_Step_2 Ra_3_Step_3 Ra_6_Step_1 Ra_6_Step_2 Ra_6_Step_3 Ra_9_Step_1 Ra_9_Step_2 Ra_9_Step_3

# Create SGWT objects for all 9 patterns and compute similarities
similarity_results <- list()
sgwt_objects <- list()

cat("\nProcessing patterns and computing SGWT analysis...\n")

## 
## Processing patterns and computing SGWT analysis...

# Process each pattern
for (i in seq_along(patterns_9)) {
  pattern_name <- names(patterns_9)[i]
  pattern_data <- patterns_9[[i]]
  
  cat("Processing", pattern_name, "...\n")
  
  # Create SGWT object
  SG_temp <- initSGWT(
    data.in = pattern_data,
    x_col = "X", 
    y_col = "Y", 
    signals = c("signal_1", "signal_2"),
    J = 4,
    kernel_type = "heat"
  )
  
  # Build graph and run SGWT
  SG_temp <- runSpecGraph(SG_temp, k = 12, laplacian_type = "normalized", 
                         length_eigenvalue = 50, verbose = FALSE)
  SG_temp <- runSGWT(SG_temp, verbose = FALSE)
  
  sgwt_objects[[pattern_name]] <- SG_temp
  
  # Compute within-pattern similarity (signal_1 vs signal_2)
  sim_within <- runSGCC("signal_1", "signal_2", SG = SG_temp, return_parts = TRUE)
  similarity_results[[pattern_name]] <- sim_within
}

## Processing Ra_3_Step_1 ...
## Processing Ra_3_Step_2 ...
## Processing Ra_3_Step_3 ...
## Processing Ra_6_Step_1 ...
## Processing Ra_6_Step_2 ...
## Processing Ra_6_Step_3 ...
## Processing Ra_9_Step_1 ...
## Processing Ra_9_Step_2 ...
## Processing Ra_9_Step_3 ...

# Print similarity results summary
cat("\nSimilarity Results Summary (signal_1 vs signal_2 within each pattern):\n")

## 
## Similarity Results Summary (signal_1 vs signal_2 within each pattern):

similarity_df <- data.frame(
  Pattern = names(similarity_results),
  S = sapply(similarity_results, function(x) x$S),
  c_low = sapply(similarity_results, function(x) x$c_low),
  c_nonlow = sapply(similarity_results, function(x) x$c_nonlow),
  w_low = sapply(similarity_results, function(x) x$w_low),
  w_NL = sapply(similarity_results, function(x) x$w_NL),
  stringsAsFactors = FALSE
)

# Extract Ra and Step values for better labeling
similarity_df$Ra <- as.numeric(gsub(".*Ra_([0-9.]+)_Step.*", "\\1", similarity_df$Pattern))
similarity_df$Step <- as.numeric(gsub(".*Step_([0-9]+)$", "\\1", similarity_df$Pattern))
similarity_df$Label <- paste0("Ra=", similarity_df$Ra, ",Step=", similarity_df$Step)

# Create similarity space visualization
similarity_plot <- visualize_similarity_xy(
  similarity_results,
  point_size = 4,
  point_color = "steelblue",
  add_diagonal = TRUE,
  add_axes_lines = TRUE,
  title = "SGWT Similarity Space: 9 Paired Patterns (signal_1 vs signal_2)",
  show_labels = TRUE
)

print(similarity_plot)

# Analysis of patterns
cat("\nPattern Analysis:\n")

## 
## Pattern Analysis:

cat("Ra (Inner Radius) Effect:\n")

## Ra (Inner Radius) Effect:

for (ra in unique(similarity_df$Ra)) {
  subset_data <- similarity_df[similarity_df$Ra == ra, ]
  cat(sprintf("  Ra=%g: Mean c_low=%.3f, Mean c_nonlow=%.3f, Mean S=%.3f\n", 
              ra, mean(subset_data$c_low), mean(subset_data$c_nonlow), mean(subset_data$S)))
}

##   Ra=3: Mean c_low=0.459, Mean c_nonlow=0.216, Mean S=0.382
##   Ra=6: Mean c_low=0.205, Mean c_nonlow=-0.153, Mean S=0.097
##   Ra=9: Mean c_low=-0.018, Mean c_nonlow=-0.300, Mean S=-0.102

cat("\nStep (Shrinkage Step) Effect:\n")

## 
## Step (Shrinkage Step) Effect:

for (step in unique(similarity_df$Step)) {
  subset_data <- similarity_df[similarity_df$Step == step, ]
  cat(sprintf("  Step=%g: Mean c_low=%.3f, Mean c_nonlow=%.3f, Mean S=%.3f\n", 
              step, mean(subset_data$c_low), mean(subset_data$c_nonlow), mean(subset_data$S)))
}

##   Step=1: Mean c_low=-0.266, Mean c_nonlow=-0.489, Mean S=-0.334
##   Step=2: Mean c_low=0.185, Mean c_nonlow=-0.256, Mean S=0.050
##   Step=3: Mean c_low=0.726, Mean c_nonlow=0.508, Mean S=0.661

cat("\nSimilarity Space Interpretation:\n")

## 
## Similarity Space Interpretation:

cat("- Each point represents similarity between signal_1 (inner circle) and signal_2 (outer ring)\n")

## - Each point represents similarity between signal_1 (inner circle) and signal_2 (outer ring)

cat("- Color indicates inner radius (Ra): different colors for different inner radii\n")

## - Color indicates inner radius (Ra): different colors for different inner radii

cat("- Shape indicates shrinkage step: circle/triangle/square for different steps\n")

## - Shape indicates shrinkage step: circle/triangle/square for different steps

cat("- Position shows frequency-domain similarity characteristics\n")

## - Position shows frequency-domain similarity characteristics

Advanced Workflow Features

Custom Parameters

# Demonstrate advanced parameter customization
SG_advanced <- initSGWT(
  data.in = demo_pattern,
  x_col = "X", y_col = "Y",
  signals = "signal_1",
  J = 6,                     # More scales for finer analysis
  scaling_factor = 1.5,      # Closer scales
  kernel_type = "heat"       # Heat kernel
)

SG_advanced <- runSpecGraph(SG_advanced, k = 15, laplacian_type = "randomwalk",length_eigenvalue = 30, verbose = FALSE)
SG_advanced <- runSGWT(SG_advanced, verbose = FALSE)

cat("Advanced Parameters Results:\n")

## Advanced Parameters Results:

cat("Number of scales:", length(SG_advanced$Parameters$scales), "\n")

## Number of scales: 6

cat("Scales:", paste(round(SG_advanced$Parameters$scales, 4), collapse = ", "), "\n")

## Scales: 0.5399, 0.36, 0.24, 0.16, 0.1067, 0.0711

cat("Reconstruction error:", round(SG_advanced$Inverse$signal_1$reconstruction_error, 6), "\n")

## Reconstruction error: 0.538056

Summary

Key Takeaways

New Object-Oriented Workflow: Clear separation of initialization, graph building, analysis, and similarity computation
Flexible Input: Supports custom column naming for coordinates and signals
Multiple Kernels: Choose from Mexican Hat, Meyer, or Heat kernel families
Comprehensive Similarity: Energy-normalized weighted similarity with detailed diagnostics
Advanced Analysis: Low-frequency only analysis, cross-object comparison, multiple signal support

New Workflow Benefits

# Complete workflow in 5 clear steps:
SG <- initSGWT(data, signals = c("signal1", "signal2"))           # 1. Initialize
kband <- checkKband(SG, signals = c("signal1", "signal2"))       # 2. Check k-band (optional)
SG <- runSpecGraph(SG, k = 12, laplacian_type = "normalized")    # 3. Build graph  
SG <- runSGWT(SG)                                                 # 4. Run SGWT
similarity <- runSGCC("signal1", "signal2", SG = SG)             # 5. Compare patterns

Parameter Guidelines

Small datasets (<200 points): k=8-12, J=3-4
Medium datasets (200-1000 points): k=12-20, J=4-6
Large datasets (>1000 points): k=15-25, J=5-8

Next Steps

Explore your own spatial datasets using the demonstrated workflow
Experiment with different kernel types and parameters
Use comprehensive similarity analysis for comparative studies
Refer to function documentation for detailed parameter specifications

BioGSP Package - Biological Graph Signal Processing for Spatial Data Analysis

BioGSP: Spectral Graph Wavelet Transform for Spatial Data (New Workflow)

Yuzhou

2026-01-28