Algorithm and Methodology

The Double-Dipping Problem

In standard scRNA-seq analysis pipelines:

1. Clustering: Cells are grouped based on gene expression patterns
2. Differential Expression Testing: Genes are tested for differential expression between clusters

This creates a statistical problem known as double-dipping: the same data is used both to define groups and to test for differences between them. This leads to:

• Inflated test statistics
• Overly optimistic p-values
• Increased false discovery rates
• Over-clustering: finding more clusters than biologically meaningful

The Knockoff Filter

The knockoff filter (Barber & Candès, 2015) is a powerful method for variable selection with FDR control. Given original features \(X_1, \ldots, X_p\), we construct knockoff features \(\tilde{X}_1, \ldots, \tilde{X}_p\) that:

1. Preserve the correlation structure of the original features
2. Are exchangeable with the originals conditional on the response
3. Are known a priori to be null (not associated with the response)

The W Statistic

For each gene \(j\), we compute the knockoff W statistic:

\[W_j = -\log_{10}(p_j^{\text{original}}) - (-\log_{10}(p_j^{\text{knockoff}}))\]

Where: - \(p_j^{\text{original}}\) is the p-value for gene \(j\) in differential expression testing - \(p_j^{\text{knockoff}}\) is the p-value for the corresponding knockoff gene

FDR-Controlled Selection

The knockoff filter selects genes where \(W_j \geq T\), where the threshold \(T\) is chosen to control FDR at level \(q\):

\[T = \min\left\{t > 0: \frac{1 + \#\{j: W_j \leq -t\}}{\#\{j: W_j \geq t\}} \leq q\right\}\]

This provides rigorous FDR control without distributional assumptions.

Stage 1: Synthetic Null Variable Generation

Given the expression matrix \(X \in \mathbb{R}^{n \times p}\) (cells × genes), we generate synthetic “knockoff” genes \(\tilde{X} \in \mathbb{R}^{n \times p}\) by:

1. Estimating marginal distributions for each gene
2. Sampling from the estimated distribution to create knockoffs

Stage 2: Joint Analysis

The augmented matrix \([X, \tilde{X}] \in \mathbb{R}^{n \times 2p}\) undergoes:

1. Normalization: Log-normalization
2. Feature Selection: Selecting top variable features from both original and knockoff genes
3. Dimensionality Reduction: PCA on the augmented data
4. Clustering: Louvain or Leiden algorithm

Stage 3: Iterative Calibration

Algorithm: FindClustersRecall
─────────────────────────────────────────
Input: Seurat object, resolution_start, q (target FDR)
Output: Calibrated cluster assignments

1. Generate knockoff features
2. resolution ← resolution_start
3. REPEAT:
   a. Cluster augmented data at current resolution
   b. FOR each cluster pair (i, j):
      - Compute differential expression (original vs knockoff)
      - Apply knockoff filter at FDR = q
      - IF no genes selected:
           Mark pair as "not significantly different"
           BREAK
   c. IF all pairs significantly different:
      RETURN cluster assignments
   d. resolution ← resolution × (1 - reduction_percentage)
4. RETURN cluster assignments

Maximum Likelihood Estimation

For count data \(x_1, \ldots, x_n\) following ZIP(\(\lambda\), \(\pi\)):

\[P(X = k) = \begin{cases} \pi + (1-\pi)e^{-\lambda} & k = 0 \\ (1-\pi)\frac{\lambda^k e^{-\lambda}}{k!} & k > 0 \end{cases}\]

The MLE has a closed-form solution using the Lambert W function:

\[\hat{\lambda} = W_0(-\gamma e^{-\gamma}) + \gamma\]

where \(\gamma = \bar{x}/(1-r_0)\) and \(r_0\) is the observed zero proportion.

Implementation

library(recall)

# Simulate ZIP data
set.seed(42)
true_lambda <- 2.5
true_pi <- 0.3

data <- rzipoisson(1000, lambda = true_lambda, prop.zero = true_pi)

# Estimate parameters
params <- estimate_zi_poisson(data)
print(params)
# $lambda.hat: ~2.5
# $pi.hat: ~0.3

Gaussian Copula

1. Fit marginal distributions to each gene
2. Transform to uniform marginals via probability integral transform
3. Fit Gaussian copula to capture dependencies
4. Simulate by sampling from copula and transforming back

# Use copula-based knockoffs for better correlation modeling
pbmc <- FindClustersRecall(
  pbmc,
  null_method = "NB-copula",  # Negative Binomial with Gaussian copula
  cores = 4
)

FDR Control

Under the knockoff framework, recall provides finite-sample FDR control:

\[\text{FDR} = \mathbb{E}\left[\frac{V}{R \vee 1}\right] \leq q\]

where \(V\) is the number of false discoveries and \(R\) is the total number of discoveries.

Comparison with Count Splitting