Skip to content
bayesrxPublic

About

No description, website, or topics provided.

Resources

Stars

0 stars

Watchers

0 watching

Forks

Latest commit

 

History

1 Commit

Folders and files

Repository files navigation

mSGR Overview

My Amazing Plot The $\texttt{mSGR}$ package provides a robust Bayesian framework for inferring spatially varying gene-gene networks from multi-resolution spatial transcriptomics (ST) data.

Installation

You can install the released version of mSGR from (https://github.com/clychenly/mSGR) with

devtools::install_github("clychenly/mSGR")

Functions

This section provides a detailed overview of the core functions within the mSGR package, including their purpose, required parameters, and return values.

1. mSGR_regression()

Performs spatially varying graph regression across multiple distinct samples (FOVs) using a Gaussian Process approximation. This function regresses a single node $Y_i$ against all other nodes $Y_j$.

Parameters:

  • Y_i: Numeric vector of length $\sum N_k$. The outcome variable (node of interest) concatenated across all samples.
  • Y_j: Numeric matrix of dimension $(\sum N_k) \times (p-1)$. The covariates (other nodes) used to regress against $Y_i$.
  • location_fov: Data frame with columns including x, y coordinates and FOV (sample) index. The column order should be c("x", "y", "fov").
  • n_k_list: Integer vector of length $K$. Contains the number of observations $N_k$ in each sample/FOV, ordered by FOV.
  • scale_Y: Logical. If TRUE, $Y_i$ and columns of $Y_j$ are standardized to mean 0 and variance 1.
  • scale_location: Logical. If TRUE, location coordinates are standardized.
  • a_gp: Numeric. Hyperparameter for the Gaussian Process (scale).
  • b_gp: Numeric. Hyperparameter for the Gaussian Process (decay).
  • poly_degree: Integer. The degree of the polynomial approximation for the Gaussian Process (determines the number of eigenfunctions $L$).
  • a_w, b_w: Numeric. Hyperparameters for the variance $\omega_{ii}^2$.
  • cor_lambda: Matrix. The correlation matrix for the latent variable $\Lambda$.
  • lambda_prior_mean: Numeric or Vector. The prior mean for $\Lambda$. Can be a scalar or a vector of length $K$.
  • init_step: Logical. If TRUE, performs an initialization step using variational Bayes.
  • max_tol: Numeric. Convergence tolerance threshold.
  • max_iter: Integer. Maximum number of iterations.
  • damped: Logical. If TRUE, applies damped fixed-point updates (also known as relaxation) to stabilize convergence. Default is TRUE.
  • step_size: Numeric. The weight $\alpha \in (0, 1]$ given to the new estimate in the damped update equation.

Output:

A list containing:

  • E_delta: Matrix of dimensions $(p-1) \times K$. Expected probability of edge inclusion.
  • E_u_i: Array of dimensions $L \times (p-1) \times K$. Expected basis coefficients.
  • E_w_i: Vector of length $K$. Expected diagonal elements for precision matrix $\Omega$.
  • Tolerance: Numeric. Final convergence tolerance value.

2. mSGR_network()

Estimates the full graph structure across multiple samples/FOVs.

Parameters:

  • Y: Numeric matrix of dimension $(\sum N_k) \times p$. The data matrix containing all nodes concatenated across all samples.
  • location_fov: Data frame with columns including x, y coordinates and FOV (sample) index. The column order should be c("x", "y", "fov").
  • n_k_list: Integer vector of length $K$. Contains the number of observations in each sample.
  • scale_Y: Logical. If TRUE, Y is column-standardized.
  • scale_location: Logical. If TRUE, location coordinates are standardized.
  • a_gp: Numeric. Hyperparameter for the Gaussian Process (scale).
  • b_gp: Numeric. Hyperparameter for the Gaussian Process (decay).
  • poly_degree: Integer. The degree of the polynomial approximation for the Gaussian Process.
  • a_w, b_w: Hyperparameters for the variance $\omega_{ii}^2$.
  • cor_lambda: Matrix. Correlation matrix for $\Lambda$.
  • lambda_prior_mean: Numeric/Vector. Prior mean for $\Lambda$.
  • init_step: Logical. If TRUE, performs an initialization step using variational Bayes.
  • max_tol: Numeric. Convergence tolerance threshold.
  • max_iter: Integer. Maximum number of iterations.
  • damped: Logical. If TRUE, applies damped fixed-point updates (relaxation).
  • step_size: Numeric. The weight $\alpha$ given to the new estimate in the damped update equation.

Output:

A list containing:

  • E_delta: Array of dimensions $p \times p \times K$. Expected inclusion probabilities.
  • E_u: List of length $K$. Each element is a $p \times p \times L$ array of expected basis coefficients.
  • E_w: Matrix of dimensions $p \times K$. Posterior sample weights.
  • Tolerance: Vector of length $p$. Final tolerance for each node’s regression.

3. selection_posthoc_multiplefov()

Applies an “And” or “Or” rule to combine edge inclusion probabilities and controls FDR for multiple FOVs/samples.

Parameters:

  • pr_array: Array of dimensions $p \times p \times K$. Expected probabilities of inclusion.
  • rule: Character. Either "and" or "or".
  • alpha: Numeric. Significance level (FDR threshold).

Output:

A list containing:

  • inclusion_array: Binary array of inclusion indicators.
  • inclusion_probability: Data frame containing detailed inclusion stats per edge/FOV.
  • cutoff: Numeric. The probability threshold used.

4. PD_rescale_factors_multiplefov()

Computes scaling factors to help ensure the estimated precision matrices are positive definite.

Parameters:

  • location: Matrix or data frame of coordinates, ordered by FOV.
  • E_u: List of expected basis coefficient arrays.
  • inclusion_array: Binary array of selected edges.
  • n_k_list: Integer vector of sample sizes for each FOV, ordered by FOV.
  • poly_degree: Integer. GP hyperparameter.
  • a_gp: Numeric. GP hyperparameter.
  • b_gp: Numeric. GP hyperparameter.

Output:

  • A list of length $K$, where each element is a numeric vector of length $p$ containing the scaling factors.

5. spatial_correlation_multiplefov()

Computes the partial correlation between two specific nodes across space, but only for FOVs where the edge is selected.

Parameters:

  • location_fov: Data frame with columns including x, y coordinates and FOV (sample) index. The column order should be c("x", "y", "fov").
  • gene_idx: Integer vector of length 2. The indices of the two nodes.
  • msgr_network_res: Output from mSGR_network.
  • selection_res: Output from selection_posthoc_multiplefov.
  • poly_degree: Integer. GP hyperparameter.
  • a_gp: Numeric. GP hyperparameter.
  • b_gp: Numeric. GP hyperparameter.
  • PD_rescale: Logical. Whether to apply PD rescaling.
  • PD_rescale_factor: Output from PD_rescale_factors_multiplefov.
  • rule: Character. Either "min" or "max". Keep $E_u$ with smaller norm or larger norm.
  • scale_location: Logical.

Output:

  • A data frame containing coordinates and the computed rho (partial correlation).

6. calculate_precision_mat_multiplefov()

Reconstructs the precision matrix for specific spatial points based on the estimated model parameters.

Parameters:

  • location_fov_scaled: Data frame of scaled coordinates and FOV ID.
  • msgr_network_res: Output from mSGR_network.
  • selection_res: Output from selection_posthoc_multiplefov.
  • poly_degree: Integer. GP hyperparameter.
  • a_gp: Numeric. GP hyperparameter.
  • b_gp: Numeric. GP hyperparameter.
  • PD_rescale: Logical. Whether to apply PD rescaling.
  • PD_rescale_factor: Output from PD_rescale_factors_multiplefov.
  • rule: Character. Either "min" or "max".

Output:

  • A list of precision matrices, one for each row in location_fov_scaled.

7. singleSGR_regression()

Performs spatially varying graph regression on a single sample.

Parameters:

  • Y_i: Numeric vector of length $N$. The outcome variable.
  • Y_j: Numeric matrix of dimension $N \times (p-1)$. Covariates.
  • location: Data frame or matrix with 2 columns (x, y).
  • scale_Y: Logical. If TRUE, Y_i and columns of Y_j are standardized to mean 0 and variance 1.
  • scale_location: Logical. If TRUE, location coordinates are standardized.
  • a_gp: Numeric. Hyperparameter for the Gaussian Process (scale).
  • b_gp: Numeric. Hyperparameter for the Gaussian Process (decay).
  • poly_degree: Integer. The degree of the polynomial approximation for the Gaussian Process.
  • a_w, b_w: Numeric. Hyperparameters for the variance $\omega_{ii}^2$.
  • pi: Numeric. Prior probability of edge inclusion (sparsity parameter).
  • init_step: Logical. If TRUE, performs an initialization step using variational Bayes.
  • max_tol: Numeric. Convergence tolerance threshold.
  • max_iter: Integer. Maximum number of iterations.
  • damped: Logical. If TRUE, applies damped fixed-point updates.
  • step_size: Numeric. The weight $\alpha \in (0, 1]$ given to the new estimate in the damped update equation.

Output:

A list containing:

  • E_delta: Vector of length $p-1$. Expected inclusion probabilities.
  • E_u: Matrix of dimension $(p-1) \times L$. Expected basis coefficients.
  • E_w: Scalar. Expected precision parameters.
  • Tolerance: Scalar. Final tolerance.

8. singleSGR_network()

Estimates the graph structure for a single sample. Parameters:

  • Y: Numeric matrix of dimension $N \times p$.
  • location: Data frame or matrix with 2 columns.
  • scale_Y: Logical. If TRUE, Y is column-standardized.
  • scale_location: Logical. If TRUE, location coordinates are standardized.
  • a_gp: Numeric. Hyperparameter for the Gaussian Process (scale).
  • b_gp: Numeric. Hyperparameter for the Gaussian Process (decay).
  • poly_degree: Integer. The degree of the polynomial approximation for the Gaussian Process.
  • a_w, b_w: Hyperparameters for the variance $\omega_{ii}^2$.
  • pi: Numeric. Prior probability of edge inclusion.
  • init_step: Logical. If TRUE, performs an initialization step using variational Bayes.
  • max_tol: Numeric. Convergence tolerance threshold.
  • max_iter: Integer. Maximum number of iterations.
  • damped: Logical. If TRUE, applies damped fixed-point updates.
  • step_size: Numeric. The weight $\alpha \in (0, 1]$ given to the new estimate in the damped update equation.

Output:

A list containing:

  • E_delta: Matrix of dimensions $p \times p$. Expected inclusion probabilities.
  • E_u: Array of dimensions $p \times p \times L$. Expected basis coefficients.
  • E_w: Vector of length $p$. Expected diagonal elements for precision matrix.
  • Tolerance: Vector of length $p$. Final tolerance.

Example

library(mSGR)
library(reshape2)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
data("example")

#Preprocess Data
Y <- example$Y
Y <- lapply(Y, function(x){
  output <- apply(x,1,scale)
  return(output)
})
Y <- do.call(rbind,Y)
cat("Structure of Y\n")
#> Structure of Y
print(str(Y))
#>  num [1:12500, 1:30] -1.35 -0.102 -1.52 -0.493 0.6 ...
#> NULL
n_k_list <- lapply(example$Y, ncol) %>% unlist
cat("\nStructure of n_k_list\n")
#> 
#> Structure of n_k_list
print(str(n_k_list))
#>  Named int [1:25] 500 500 500 500 500 500 500 500 500 500 ...
#>  - attr(*, "names")= chr [1:25] "fov1" "fov2" "fov3" "fov4" ...
#> NULL
location <- example$location_used

# Set hyperparameters
a_gp = 0.01; b_gp = 0.5; poly_degree = 10;
a_w = 10; b_w = 10;
cor_lambda = diag(25)


# Run mSGR_regression
msgr_regression_res <- mSGR_regression(Y[,1], Y[,-1],
                                       location, n_k_list,
                                       scale_Y = TRUE, scale_location = TRUE,
                                       a_gp, b_gp, poly_degree,
                                       a_w, b_w,
                                       cor_lambda, lambda_prior_mean = 0,
                                       init_step = TRUE,
                                       max_tol = 0.001, max_iter = 2000,
                                       damped = TRUE, step_size = 0.9)
#> [1] "Initialization Done!"
cat("\nStructure of results from mSGR_regression\n")
#> 
#> Structure of results from mSGR_regression
str(msgr_regression_res)
#> List of 4
#>  $ E_delta  : num [1:29, 1:25] 0.00993 0.02667 0.01542 0.04157 0.03305 ...
#>  $ E_u_i    : num [1:66, 1:29, 1:25] -0.000863 -0.000157 -0.00115 0.000314 -0.000118 ...
#>  $ E_w_i    : num [1:25, 1] 1.52 1.02 2.05 1.49 1.04 ...
#>  $ Tolerance: num 0.000967

# Run mSGR_network
msgr_network_res <- mSGR_network(Y,
                                 location, n_k_list,
                                 scale_Y = TRUE, scale_location = TRUE,
                                 a_gp = 0.01, b_gp = 0.5, poly_degree = 10,
                                 a_w = 10, b_w = 10,
                                 cor_lambda = diag(25), lambda_prior_mean = 0,
                                 init_step = TRUE,
                                 max_tol = 0.001, max_iter = 2000,
                                 damped = TRUE, step_size = 0.9)
#> [1] "Initialization Done for node 1"
#> [1] "Initialization Done for node 2"
#> [1] "Initialization Done for node 3"
#> [1] "Initialization Done for node 4"
#> [1] "Initialization Done for node 5"
#> [1] "Initialization Done for node 6"
#> [1] "Initialization Done for node 7"
#> [1] "Initialization Done for node 8"
#> [1] "Initialization Done for node 9"
#> [1] "Initialization Done for node 10"
#> [1] "Initialization Done for node 11"
#> [1] "Initialization Done for node 12"
#> [1] "Initialization Done for node 13"
#> [1] "Initialization Done for node 14"
#> [1] "Initialization Done for node 15"
#> [1] "Initialization Done for node 16"
#> [1] "Initialization Done for node 17"
#> [1] "Initialization Done for node 18"
#> [1] "Initialization Done for node 19"
#> [1] "Initialization Done for node 20"
#> [1] "Initialization Done for node 21"
#> [1] "Initialization Done for node 22"
#> [1] "Initialization Done for node 23"
#> [1] "Initialization Done for node 24"
#> [1] "Initialization Done for node 25"
#> [1] "Initialization Done for node 26"
#> [1] "Initialization Done for node 27"
#> [1] "Initialization Done for node 28"
#> [1] "Initialization Done for node 29"
#> [1] "Initialization Done for node 30"
cat("\nStructure of results from mSGR_network\n")
#> 
#> Structure of results from mSGR_network
str(msgr_network_res)
#> List of 4
#>  $ E_delta  : num [1:30, 1:30, 1:25] 0 0.011 0.0245 0.0133 0.0306 ...
#>  $ E_u      :List of 25
#>   ..$ : num [1:30, 1:30, 1:66] 0.00 -7.39e-04 -1.43e-03 -8.66e-05 -7.45e-04 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 -0.000146 -0.000037 0.00159 -0.004077 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 0.000743 0.001962 0.000139 0.361225 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 -0.000848 0.000714 0.001043 0.007218 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 -0.000443 0.001298 0.000637 -0.018415 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 -0.000151 0.001753 0.00453 0.026178 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0.00 -7.40e-06 -4.27e-04 -5.91e-04 -6.34e-04 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 -0.002281 0.00295 0.000978 0.166958 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 -0.00103 0.000882 0.000404 0.201883 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 0.000135 0.000215 0.001537 0.045177 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 -0.007683 -0.000272 -0.000262 0.001211 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 0.005532 0.000371 0.054778 -0.000506 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 0.084695 -0.001032 0.000151 0.14958 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 0.000562 0.001924 -0.000057 0.223391 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0.00 4.75e-03 7.56e-05 -1.91e-03 2.64e-02 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 -0.00011 0.26509 -0.00023 0.00183 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0.00 8.31e-04 5.34e-04 1.22e-03 -6.75e-05 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 -0.001596 0.017404 0.002448 0.000146 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0.00 -2.19e-03 -1.18e-03 -4.22e-05 -2.90e-04 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 -0.000236 0.002656 -0.000823 -0.000784 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 0.001534 -0.000502 0.003934 0.000581 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 0.003187 -0.004772 0.000054 0.000105 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 -0.001128 0.000399 -0.000983 -0.000454 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 -0.000354 0.000234 0.001381 -0.268825 ...
#>   ..$ : num [1:30, 1:30, 1:66] 0 -0.002786 -0.003235 -0.000728 -0.001008 ...
#>  $ E_w      : num [1:30, 1:25] 1.52 1.53 1.02 1.32 1.14 ...
#>  $ Tolerance: num [1:30] 0.000937 0.000992 0.000987 0.000978 0.000995 ...

#Post-hoc procedure for edge selection
selection_res <- 
  selection_posthoc_multiplefov(msgr_network_res$E_delta,"and",0.01)

cat("\nStructure of results from selection_posthoc_multiplefov\n")
#> 
#> Structure of results from selection_posthoc_multiplefov
str(selection_res)
#> List of 3
#>  $ inclusion_array      : num [1:30, 1:30, 1:25] 0 0 0 0 0 0 0 0 0 1 ...
#>  $ inclusion_probability: tibble [10,875 × 8] (S3: tbl_df/tbl/data.frame)
#>   ..$ node_idx_max: int [1:10875] 5 8 9 10 10 11 13 15 15 16 ...
#>   ..$ node_idx_min: int [1:10875] 3 5 5 3 8 1 6 10 11 4 ...
#>   ..$ FOV         : int [1:10875] 24 24 25 23 21 24 22 23 21 21 ...
#>   ..$ pr_max      : num [1:10875] 1 1 1 1 1 1 1 1 1 1 ...
#>   ..$ pr_min      : num [1:10875] 1 1 1 1 1 1 1 1 1 1 ...
#>   ..$ q           : num [1:10875] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ fdr_p       : num [1:10875] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ selected    : num [1:10875] 1 1 1 1 1 1 1 1 1 1 ...
#>  $ cutoff               : num 0.781

#To ensure PD, calcualte rescale factors
rescale_res <- PD_rescale_factors_multiplefov(location[,c(1:2)],
                                              msgr_network_res$E_u,
                                              selection_res$inclusion_array,
                                              n_k_list,
                                              poly_degree, a_gp, b_gp)
cat("\nStructure of results from PD_rescale_factors_multiplefov\n")
#> 
#> Structure of results from PD_rescale_factors_multiplefov
str(rescale_res)
#> List of 25
#>  $ : num [1:30] 4.83 4.46 1 3.46 1.68 ...
#>  $ : num [1:30] 1 1.55 1.84 2.13 3.05 ...
#>  $ : num [1:30] 5.24 2.47 2.12 1.51 9.13 ...
#>  $ : num [1:30] 2.68 1.7 5.78 5.1 1.93 ...
#>  $ : num [1:30] 0.898 5.011 2.882 5.04 6.464 ...
#>  $ : num [1:30] 1.8 9 1 1.72 1 ...
#>  $ : num [1:30] 5.53 5.34 3.31 1 9.32 ...
#>  $ : num [1:30] 7.96 3.68 5.18 4.05 2.93 ...
#>  $ : num [1:30] 6.8 6.01 1 6.65 5.61 ...
#>  $ : num [1:30] 1.94 6.13 7.88 3.32 7.15 ...
#>  $ : num [1:30] 1.96 4.53 1 9.11 5.4 ...
#>  $ : num [1:30] 1 9.72 2.5 4.48 7.1 ...
#>  $ : num [1:30] 4.71 1.5 6.78 5.27 8.04 ...
#>  $ : num [1:30] 8.59 7.34 1 3.39 6.06 ...
#>  $ : num [1:30] 4.36 5.35 4.43 4.35 7.35 ...
#>  $ : num [1:30] 6.97 5.36 3.87 8.02 2.06 ...
#>  $ : num [1:30] 2.96 1 1.1 7.03 8.83 ...
#>  $ : num [1:30] 1 4.34 2.04 4.3 7.08 ...
#>  $ : num [1:30] 10.05 1 4.99 1.89 7.04 ...
#>  $ : num [1:30] 3.8 2.17 3.28 2.49 3.8 ...
#>  $ : num [1:30] 3.05 2.76 4.78 6.94 1 ...
#>  $ : num [1:30] 1 4.87 6.44 4.68 1.2 ...
#>  $ : num [1:30] 1.37 1 6.26 1.35 1 ...
#>  $ : num [1:30] 5.08 1 4.06 2.97 7.04 ...
#>  $ : num [1:30] 2.12 1 1.82 3.81 6.32 ...

#Calculate spatial correlation pattern for selected gene pairs
spatial_cor_res <- spatial_correlation_multiplefov(location, c(24,26),
                                msgr_network_res,
                                selection_res,
                                poly_degree, a_gp, b_gp,
                                PD_rescale = TRUE,
                                PD_rescale_factor = rescale_res,
                                rule = "min")
cat("\nStructure of results from spatial_correlation_multiplefov\n")
#> 
#> Structure of results from spatial_correlation_multiplefov
str(spatial_cor_res)
#> 'data.frame':    12500 obs. of  6 variables:
#>  $ x               : num  -0.79 -0.479 -0.73 -0.479 -0.551 ...
#>  $ y               : num  0.884 0.705 0.705 0.717 0.465 ...
#>  $ fov             : num  7 7 7 7 7 7 7 7 7 7 ...
#>  $ x_original_scale: num  -0.677 -0.41 -0.626 -0.41 -0.472 ...
#>  $ y_original_scale: num  0.759 0.605 0.605 0.615 0.4 ...
#>  $ rho             : num  0.124 0.141 0.117 0.142 0.112 ...

#Calculate the precision matrix given location
location_scaled <- data.frame(x = (location$x - mean(location$x))/sd(location$x),
                              y = (location$y - mean(location$y))/sd(location$y),
                              fov = location$FOV)
pre_mat <- calculate_precision_mat_multiplefov(
  location_scaled[c(1:2), ],
  msgr_network_res,
  selection_res,
  poly_degree,
  a_gp,
  b_gp,
  PD_rescale = TRUE,
  PD_rescale_factor = rescale_res,
  rule = "min"
)
cat("\nStructure of results from calculate_precision_mat_multiplefov\n")
#> 
#> Structure of results from calculate_precision_mat_multiplefov
str(pre_mat)
#> List of 2
#>  $ : num [1:30, 1:30] 1.52 0 0 0 0 ...
#>  $ : num [1:30, 1:30] 1.52 0 0 0 0 ...

About

No description, website, or topics provided.

Resources

Stars

0 stars

Watchers

0 watching

Forks

Releases

Packages

Contributors

Languages