Problem Formulation

Given:

• Bulk expression matrix \(\mathbf{X} \in \mathbb{R}^{n \times p}\) (n samples × p genes)
• Phenotype vector \(\mathbf{y}\) (continuous, binary, or survival)
• Single-cell expression matrix \(\mathbf{S} \in \mathbb{R}^{m \times p}\) (m cells × p genes)

The goal is to find gene weights \(\boldsymbol{\beta}\) that associate gene expression with phenotype, then apply these weights to single-cell data to compute per-cell risk scores.

Network-Regularized Sparse Regression

scPAS uses the APML0 (Augmented and Penalized Minimization L0) algorithm with network regularization:

\[ \hat{\boldsymbol{\beta}} = \arg\min_{\boldsymbol{\beta}} \left\{ L(\boldsymbol{\beta}; \mathbf{X}, \mathbf{y}) + \lambda_1 \|\boldsymbol{\beta}\|_1 + \lambda_2 \boldsymbol{\beta}^T \mathbf{L} \boldsymbol{\beta} \right\} \]

Where:

• \(L(\boldsymbol{\beta})\) is the loss function (depends on phenotype type)
• \(\lambda_1 \|\boldsymbol{\beta}\|_1\) is the L1 penalty (LASSO) for sparsity
• \(\lambda_2 \boldsymbol{\beta}^T \mathbf{L} \boldsymbol{\beta}\) is the Laplacian penalty for network regularization
• \(\mathbf{L}\) is the Laplacian matrix of the gene-gene network

Effect of Network Regularization

Loss Functions by Phenotype Type

Shared Nearest Neighbor (SNN) Network

scPAS constructs a gene-gene similarity network from single-cell data using the SNN algorithm:

SNN Network Construction

The network construction process:

1. Calculate gene correlations from single-cell expression
2. Find k-nearest neighbors for each gene
3. Compute SNN similarity based on shared neighbors
4. Threshold to create binary adjacency matrix

Per-Cell Risk Score

Once the model is trained, the risk score for each cell is computed as:

\[ RS_j = \sum_{g=1}^{p} \hat{\beta}_g \cdot \tilde{S}_{jg} \]

where:

• \(RS_j\) is the risk score for cell \(j\)
• \(\hat{\beta}_g\) is the learned coefficient for gene \(g\)
• \(\tilde{S}_{jg}\) is the standardized expression of gene \(g\) in cell \(j\)

Risk Score Distribution

Normalized Risk Score

The raw risk score is converted to a Z-statistic:

\[ NRS_j = \frac{RS_j - \mu_{bg}}{\sigma_{bg}} \]

where \(\mu_{bg}\) and \(\sigma_{bg}\) are the mean and standard deviation of the background distribution estimated from permutation.

Permutation Test

To assess significance, scPAS performs a permutation test:

Permutation Test Principle

Algorithm:

1. For each permutation \(b = 1, \ldots, B\):
   • Randomly shuffle the gene coefficients \(\boldsymbol{\beta}\)
   • Calculate permuted risk scores for all cells
2. Compute two-tailed P-value:

\[ p_j = \frac{1}{B} \sum_{b=1}^{B} \mathbb{I}(|RS_j^{(b)}| \geq |RS_j|) \]

FDR Correction

Multiple testing correction using Benjamini-Hochberg procedure:

\[ FDR_j = \min\left(1, \frac{p_j \cdot m}{\text{rank}(p_j)}\right) \]

where \(m\) is the total number of cells.

Sparse Matrix Operations

scPAS uses efficient sparse matrix operations for large-scale single-cell data:

# Efficient correlation calculation for sparse matrices
sparse.cor <- function(x) {
  # Uses optimized algorithm that avoids dense conversion
  # Handles numerical precision issues
  # Returns proper correlation matrix
}

# Efficient row scaling
sparse_row_scale <- function(x, center = TRUE, scale = TRUE) {
  # Row-wise standardization
  # Preserves sparsity when only scaling (not centering)
}

Parallel Computing

For large permutation counts, scPAS supports parallel processing:

result <- scPAS(
  bulk_dataset = bulk_data,
  sc_dataset = sc_obj,
  phenotype = phenotype,
  permutation_times = 5000,
  n_cores = 4  # Use 4 CPU cores
)