Algorithm and Methodology

Algorithm Overview

The CNV detection follows these steps:

1. Gene Annotation: Map genes to chromosomal positions using hg20 annotations
2. Expression Binning: Aggregate expression into genomic bins (220 bins per chromosome)
3. Baseline Normalization: Use normal cells or synthetic baseline as reference
4. Smoothing: Apply Kalman filtering to reduce noise
5. Clustering: Hierarchical clustering to identify aneuploid populations

Mathematical Framework

For a gene \(g\) at genomic position \(p\), the relative copy number is estimated as:

\[\text{CNV}_{g} = \log_2\left(\frac{E_g^{tumor}}{E_g^{baseline}} + 1\right)\]

where \(E_g\) represents the normalized expression level.

CNV profile across chromosomes

Why MCA?

Unlike PCA which assumes continuous data, MCA is designed for categorical/count data and provides:

• Joint embedding of cells and genes
• Natural handling of sparse scRNA-seq data
• Interpretable gene-cell associations

Mathematical Formulation

Given a normalized expression matrix \(\mathbf{X}\) with \(n\) cells and \(p\) genes:

1. Construct indicator matrix \(\mathbf{Z}\) by discretizing expression levels
2. Compute correspondence matrix: \[\mathbf{S} = \mathbf{D}_r^{-1/2}(\mathbf{Z} - \mathbf{rc}^T)\mathbf{D}_c^{-1/2}\]
3. Perform SVD: \(\mathbf{S} = \mathbf{U\Sigma V}^T\)
4. Extract coordinates:
   • Cell coordinates: \(\mathbf{F} = \mathbf{D}_r^{-1/2}\mathbf{U\Sigma}\)
   • Gene coordinates: \(\mathbf{G} = \mathbf{D}_c^{-1/2}\mathbf{V\Sigma}\)

MCA embedding showing cell-gene correspondence

C++ Implementation

scPharm implements MCA using RcppArmadillo for computational efficiency:

// Simplified MCA computation (actual implementation in src/mca.cpp)
arma::mat SparseMCAStep1(arma::sp_mat& X) {
    // Column sums and total
    arma::vec col_sum = arma::vec(arma::sum(X, 0).t());
    double total = arma::accu(X);
    
    // Compute standardized residuals
    arma::mat Z = compute_residuals(X, col_sum, total);
    
    return Z;
}

Distance Metric

The association between cell \(i\) and gene \(j\) is computed as:

\[d_{ij} = \sqrt{\sum_{k=1}^{K} (F_{ik} - G_{jk})^2}\]

where \(K\) is the number of MCA components.

Signature Extraction

For each cell, the top \(n\) genes (by distance) form its identity signature:

\[\text{Signature}_i = \{g_1, g_2, ..., g_n\} \text{ where } d_{i,g_1} < d_{i,g_2} < ... < d_{i,g_n}\]

GSEA Algorithm

1. Rank genes by their association with the cell
2. Calculate running enrichment score (ES): \[ES = \max_{j}\left(\sum_{g_i \in S, i \leq j} \sqrt{\frac{N-|S|}{|S|}} - \sum_{g_i \notin S, i \leq j} \sqrt{\frac{|S|}{N-|S|}}\right)\]
3. Normalize to obtain NES (Normalized Enrichment Score)

GSEA running enrichment score

Drug Sensitivity Gene Sets

Drug sensitivity gene sets are derived from GDSC2 database:

• For each drug, genes are ranked by correlation with IC50
• Top correlated genes → Sensitivity signature
• Bottom correlated genes → Resistance signature

Gaussian Mixture Model

Cell NES values are modeled as a mixture of Gaussian distributions:

\[P(\text{NES}) = \pi_S \cdot \mathcal{N}(\mu_S, \sigma_S^2) + \pi_R \cdot \mathcal{N}(\mu_R, \sigma_R^2) + \pi_O \cdot \mathcal{N}(\mu_O, \sigma_O^2)\]

Gaussian Mixture Model for cell classification

Classification Criteria

Cells are classified based on NES thresholds:

Drug Prioritization Score (Dr)

\[Dr = 1 - \left(\alpha \cdot \frac{n_S}{n_{tumor}} + \beta \cdot \overline{\text{NES}}_S\right)\]

where: - \(n_S\) = number of sensitive cells - \(n_{tumor}\) = total tumor cells - \(\overline{\text{NES}}_S\) = mean NES of sensitive cells - \(\alpha, \beta\) = weighting parameters

Interpretation: Lower Dr indicates better drug candidate.

Drug Side Effect Score (Dse)

\[Dse = \frac{n_{S,adj}}{n_{adj}} \cdot \overline{\text{NES}}_{S,adj}\]

where subscript \(adj\) denotes adjacent (non-tumor) cells.

Interpretation: Higher Dse indicates more potential side effects.