SecAct Algorithm: Ridge Regression with Permutation Testing

Problem Formulation

Given: - Y: Gene expression matrix (genes × samples) - X: Signature matrix (genes × secreted proteins)

We aim to infer the activity matrix β (secreted proteins × samples) that best explains the observed gene expression changes.

Ridge Regression Model

SecAct employs ridge regression to solve this inverse problem:

\[\hat{\beta} = (X^TX + \lambda I)^{-1} X^T Y\]

Where: - \(\hat{\beta}\): Estimated activity coefficients - \(\lambda\): Regularization parameter (default: 5×10⁵) - \(I\): Identity matrix

The ridge penalty \(\lambda I\) provides: 1. Numerical stability when \(X^TX\) is ill-conditioned 2. Regularization to prevent overfitting 3. Unique solution even when features > samples

Implementation via Cholesky Decomposition

For computational efficiency, SecAct uses Cholesky decomposition:

A = X'X + λI    (symmetric positive definite)
R = chol(A)     (A = R'R)
β = R⁻¹ (R')⁻¹ X' Y

This approach is: - Numerically stable: Exploits the SPD structure - Computationally efficient: O(n³/3) vs O(n³) for general inverse - Memory efficient: Only stores upper triangular R

Permutation Testing

SecAct assesses statistical significance through permutation testing:

1. Null hypothesis: No association between gene expression and secreted protein activity
2. Procedure:
   • Randomly permute sample labels in Y
   • Recompute β for each permutation
   • Build null distribution of coefficients

Permutation testing procedure

Z-score Calculation

The z-score quantifies how many standard deviations the observed coefficient deviates from the null:

\[z = \frac{\hat{\beta} - \mu_{null}}{\sigma_{null}}\]

Where: - \(\mu_{null}\): Mean of permutation coefficients - \(\sigma_{null}\): Standard deviation of permutation coefficients

P-value Computation

The empirical p-value is calculated as:

\[p = \frac{\sum_{i=1}^{n_{rand}} I(|\beta_i^{rand}| \geq |\hat{\beta}|) + 1}{n_{rand} + 1}\]

The “+1” correction ensures p-values are never exactly zero.

SecAct Signature Database

The SecAct signature matrix contains: - 1,170 secreted proteins - 7,919 downstream target genes - Curated from published literature and databases

library(SecAct)

# Download signature matrix (first time only)
SecAct.download.data()

# Load signature matrix
sig_path <- SecAct.get.sigmat.path()
sig_mat <- read.table(gzfile(sig_path), sep = "\t", header = TRUE, row.names = 1, nrows = 100)

cat("Signature matrix dimensions:\n")
cat("  Genes:", nrow(sig_mat), "(showing first 100)\n")
cat("  Secreted proteins:", ncol(sig_mat), "\n")

Signature Grouping

To reduce multicollinearity, SecAct groups highly correlated signatures:

# Demonstrate signature grouping concept
set.seed(123)
# Simulate correlation matrix for 20 signatures
n_sig <- 20
cor_mat <- matrix(runif(n_sig^2, 0.2, 0.9), n_sig, n_sig)
cor_mat <- (cor_mat + t(cor_mat)) / 2
diag(cor_mat) <- 1
rownames(cor_mat) <- colnames(cor_mat) <- paste0("SP", 1:n_sig)

# Hierarchical clustering
hc <- hclust(as.dist(1 - cor_mat), method = "complete")
plot(hc, main = "Signature Grouping by Correlation", xlab = "", sub = "")
abline(h = 0.1, col = "red", lty = 2)
text(15, 0.15, "Correlation threshold = 0.9", col = "red")

Hierarchical clustering of signature correlation

Lambda Selection

The regularization parameter λ controls the bias-variance tradeoff:

# Demonstrate lambda effect
lambdas <- c(1e3, 1e4, 1e5, 5e5, 1e6, 1e7)
Y_small <- Y[, 1:2]

results <- lapply(lambdas, function(l) {
  SecAct.inference.r(Y_small, lambda = l, nrand = 50)
})

# Plot coefficient ranges
coef_ranges <- sapply(results, function(r) range(r$beta))
par(mar = c(4, 4, 2, 1))
plot(log10(lambdas), coef_ranges[2,], type = "b", pch = 19, col = "blue",
     ylim = range(coef_ranges), xlab = "log10(lambda)", ylab = "Coefficient range",
     main = "Lambda Effect on Coefficients")
lines(log10(lambdas), coef_ranges[1,], type = "b", pch = 19, col = "red")
legend("topright", c("Max", "Min"), col = c("blue", "red"), pch = 19)
abline(v = log10(5e5), lty = 2, col = "gray")
text(log10(5e5), mean(coef_ranges), "Default", pos = 4)

Number of Permutations

More permutations yield more precise p-values:

Recommendation: Use 1000 permutations for publication-quality results.