The
You can install the released version of mSGR from (https://github.com/clychenly/mSGR) with
devtools::install_github("clychenly/mSGR")This section provides a detailed overview of the core functions within
the mSGR package, including their purpose, required parameters, and
return values.
Performs spatially varying graph regression across multiple distinct
samples (FOVs) using a Gaussian Process approximation. This function
regresses a single node
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 bec("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. IfTRUE,$Y_i$ and columns of$Y_j$ are standardized to mean 0 and variance 1. -
scale_location: Logical. IfTRUE, 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. IfTRUE, performs an initialization step using variational Bayes. -
max_tol: Numeric. Convergence tolerance threshold. -
max_iter: Integer. Maximum number of iterations. -
damped: Logical. IfTRUE, applies damped fixed-point updates (also known as relaxation) to stabilize convergence. Default isTRUE. -
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.
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 bec("x", "y", "fov"). -
n_k_list: Integer vector of length$K$ . Contains the number of observations in each sample. -
scale_Y: Logical. IfTRUE,Yis column-standardized. -
scale_location: Logical. IfTRUE, 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. IfTRUE, performs an initialization step using variational Bayes. -
max_tol: Numeric. Convergence tolerance threshold. -
max_iter: Integer. Maximum number of iterations. -
damped: Logical. IfTRUE, 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.
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.
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.
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 bec("x", "y", "fov"). -
gene_idx: Integer vector of length 2. The indices of the two nodes. -
msgr_network_res: Output frommSGR_network. -
selection_res: Output fromselection_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 fromPD_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).
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 frommSGR_network.selection_res: Output fromselection_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 fromPD_rescale_factors_multiplefov.rule: Character. Either"min"or"max".
Output:
- A list of precision matrices, one for each row in
location_fov_scaled.
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. IfTRUE,Y_iand columns ofY_jare standardized to mean 0 and variance 1. -
scale_location: Logical. IfTRUE, 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. IfTRUE, performs an initialization step using variational Bayes. -
max_tol: Numeric. Convergence tolerance threshold. -
max_iter: Integer. Maximum number of iterations. -
damped: Logical. IfTRUE, 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.
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. IfTRUE,Yis column-standardized. -
scale_location: Logical. IfTRUE, 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. IfTRUE, performs an initialization step using variational Bayes. -
max_tol: Numeric. Convergence tolerance threshold. -
max_iter: Integer. Maximum number of iterations. -
damped: Logical. IfTRUE, 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.
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 ...