MPAQT Quantification: Implemented Methods

Scope

MPAQT integrates short-read equivalence-class counts with optional long-read transcript counts, positional-bias weights, UMI normalization, and transcript-level priors. This article describes the algorithm implemented in package version 2.4.0.

The internal function is named run_em_algorithm(), but the implementation is not a classical EM routine that materializes an E-step assignment matrix. Within each iteration it performs sequential transcript-wise coordinate Newton updates and immediately updates the expected equivalence-class counts.

Notation and Short-Read Model

Let:

• $n_i$ be the observed count for short-read equivalence class (EC) $i$;
• $\beta_j \ge 0$ be the raw abundance for transcript $j$;
• $P_{ij}$ be the stored EC-by-transcript design/exposure value; and
• $y_i = \sum_j P_{ij}\beta_j$ be the expected EC count.

For a default bulk index, the nonzero x values in index$p_matrices[[j]] are counts produced by simulated reads, not columns that are automatically normalized to sum to one. Each per-transcript table stores the EC indices i, the values x, and the associated 5’ and 3’ distances.

Ignoring constants, the short-read Poisson negative log-likelihood is

$$\operatorname{NLL}_{SR}(\boldsymbol\beta) = \sum_i \left[y_i - n_i\log(y_i)\right].$$

Optional UMI normalization

do_umi_correction = TRUE normalizes each transcript’s P column to sum to one. Without positional weights, or with umi_correction_timing = "pre", this is done before the coordinate iterations. With positional weights and umi_correction_timing = "post", MPAQT multiplies by the positional weights first and normalizes the resulting weighted column during each iteration.

This P-matrix normalization is distinct from UMI molecule deduplication in mpaqt_prepare_short_reads_sc(). Read preparation defaults to do_umi_dedup = FALSE; users who want molecule-level EC counts must enable deduplication explicitly.

Sequential Coordinate Updates

MPAQT initializes all $\beta_j$ and $y_i$ to zero. It then visits transcripts in index order during every iteration.

First iteration

Because the Poisson expressions are unstable at an all-zero initialization, the first iteration uses a Gaussian residual approximation. For the current transcript column $p_j$,

$$\Delta_j = \frac{\sum_i (n_i-y_i)P_{ij}} {\sum_i P_{ij}^2}.$$

Later iterations

For subsequent iterations, with a small numerical constant $\epsilon$, the implemented short-read gradient and positive curvature terms are

$$g_j = \sum_i P_{ij} \left(1-\frac{n_i+\epsilon}{y_i+\epsilon}\right)$$

and

$$h_j = \sum_i (n_i+\epsilon) \left(\frac{P_{ij}}{y_i+\epsilon}\right)^2.$$

The unregularized coordinate step is $\Delta_j=-g_j/h_j$. Non-finite steps are replaced with zero, and the lower bound $\Delta_j \ge -\beta_j+10^{-10}$ preserves positive abundances. After each step, the implementation immediately applies

$$\beta_j \leftarrow \beta_j+\Delta_j, \qquad y_i \leftarrow y_i+P_{ij}\Delta_j.$$

Consequently, later transcripts in an iteration see the updates already made to earlier transcripts. MPAQT does not construct or store latent assignment probabilities $\gamma_{ij}$.

Log-Scale Regularization

Starting after the first iteration, positive transcript abundances receive a quadratic penalty in log space:

$$R_j = \frac{(\log\beta_j-\mu_j)^2}{2\sigma^2}.$$

The penalty contributions added to the coordinate gradient and curvature are

$$g_{R,j}=\frac{\log\beta_j-\mu_j}{\beta_j\sigma^2}$$

and

$$h_{R,j}=\frac{-\log\beta_j+\mu_j+1}{\beta_j^2\sigma^2}.$$

The regularized Newton step is bounded relative to the unregularized step and the step toward the prior mode, $\exp(\mu_j)-\beta_j$. When it lies outside that interval, MPAQT minimizes the one-dimensional penalized Poisson objective with stats::optimize() over the bounded interval.

Without a user prior, $\mu$ and $\sigma^2$ are updated from the mean and mean squared deviation of positive log abundances after each iteration. The variance is bounded below by $10^{-6}$ in post-quantification.

Long-Read Integration

When long-read counts are present, MPAQT fits a Poisson GLM after every coordinate sweep:

$$c_j \sim \operatorname{Poisson}(\beta_j q_j), \qquad q_j=\exp(C_j\boldsymbol\theta).$$

Here $c_j$ is the long-read transcript count and $C_j$ is the row of index covariates. The current index covariates are an intercept, log GC ratio, log transcript length, and a protein-coding indicator. Sequencing-depth scaling is absorbed by the fitted intercept. Although the result field is named lr_coverage_probs, $q_j$ is a positive rate/exposure multiplier and is not constrained to the interval $[0,1]$.

During an abundance coordinate update, a long-read observation is appended only when $c_j>0$. A transcript with zero long-read counts therefore does not receive a direct zero-count long-read coordinate term. Zero counts still participate in the GLM fit and in the diagnostic likelihood described below.

Prior Specifications

prior_model accepts the following forms:

Value	Implemented prior mean behavior
NULL	Uses the global mean of positive log abundances
"shrinkage"	Uses transcript prior predictions initialized to zero from prior_start
"long_read"	Uses log(long_read_count + 1) as transcript prior predictions
data frame/data table	Requires transcript_id; remaining columns are GPBoost covariates
list with mean	Uses a full-length or transcript-named custom mean vector

For a data-frame prior, transcript versions are trimmed when matching the covariate table to the index. Beginning at prior_start, GPBoost is fitted after calculation of the diagnostic likelihood and after the convergence check. Its predictions and variance therefore affect the next iteration.

To ensure a data-frame prior is used, choose convergence_start later than prior_start. The single-cell GTEx example uses:

results <- mpaqt_quant_sc(
  index = index,
  sr_counts_list = sr_list,
  prior_model = prior_data,
  prior_start = 25L,
  convergence_start = 50L
)

Positional-Bias Prequantification

Positional-bias correction is optional. When positional_bias is "3p" or "5p", mpaqt_quant() first calls mpaqt_prequant().

The selected transcript-EC distances $d_{ij}$ are divided into global quantile bins. All bin weights start at one. After the default 20-iteration warm-up, MPAQT builds an EC-by-bin expected-count matrix $S$ and minimizes

$$\operatorname{NLL}_{bias}(\mathbf w) =\sum_i\left[\lambda_i-n_i\log(\lambda_i)\right], \qquad \boldsymbol\lambda=S\exp(\log\mathbf w),$$

using L-BFGS-B. At each update, log weights are bounded to the current log weights plus or minus one. For identifiability, fitted weights are divided by the first-bin weight and abundances are multiplied by that same value.

The "3p" and "5p" settings select the distance origin; they do not impose a monotonic enrichment direction. mpaqt_prequant() returns positional weights, quantile boundaries, bias metadata, convergence information, and initial abundances. It does not accept long-read counts or fit a transcript prior.

Diagnostic Likelihood and Convergence

The scalar stored as log_likelihood is a convergence diagnostic. It is not the complete penalized objective used by the coordinate steps.

The diagnostic begins with

$$D_{SR}=\sum_i[n_i\log(y_i+\epsilon)-y_i]$$

and adds $-T\log(\sigma^2)/2$. It does not include the transcript-specific quadratic terms $(\log\beta_j-\mu_j)^2/(2\sigma^2)$.

When long reads are present, the implementation also adds

$$\sum_j[c_j\log(\beta_jq_j+\epsilon)-\beta_jq_j]$$

and the fitted GLM’s logLik() value. These contributions should therefore be interpreted as the current diagnostic calculation, not as a single separately defined likelihood optimized exactly by every coordinate update.

After convergence_start, the algorithm declares convergence only when the diagnostic change is positive and smaller than tolerance:

$$0 < D^{(t)}-D^{(t-1)} < \texttt{tolerance}.$$

There is no separate relative-abundance convergence criterion.

Two-Phase API

The main mpaqt_quant() wrapper runs prequantification only when positional bias is requested. The same work can be run explicitly:

prequant <- mpaqt_prequant(
  index = index,
  sr_counts = sr_counts,
  positional_bias = "3p",
  n_bins = 50L
)

mpaqt_save_prequant(prequant, "results/mpaqt.prequant.rds")

result <- mpaqt_postquant(
  index = index,
  sr_counts = sr_counts,
  lr_counts = lr_counts,
  prequant = prequant,
  prior_model = prior_data,
  prior_start = 25L,
  convergence_start = 50L,
  compute_uncertainty = TRUE
)

Long-read observations and priors belong to mpaqt_postquant(); they are not inputs or outputs of mpaqt_prequant().

Uncertainty

With compute_uncertainty = TRUE, MPAQT returns a diagonal, short-read-only approximation. For each positive transcript abundance,

$$I_j=\sum_i\frac{P_{ij}^2n_i}{\max(y_i,10^{-10})^2}, \qquad SE_j=I_j^{-1/2}.$$

This calculation does not include off-diagonal transcript covariance, long-read information, prior curvature, or explicit positional-weight terms. It returns standard errors only; the package does not currently return confidence intervals.

Returned Abundances and TPM

An mpaqt_quant_result stores raw coordinate estimates in abundances and a separately scaled vector in tpm:

$$\operatorname{TPM}_j=10^6\frac{\beta_j}{\sum_k\beta_k}.$$

Use abundances(result) and tpm(result) to keep these values distinct. The tpm field is calculated for every normalize setting. normalize = "depth" also stores the million-scaled values in normalized_counts; "tpm" and "none" retain raw abundances in that field.

Defaults by Entry Point

Parameter	mpaqt_quant()	mpaqt_prequant()	mpaqt_postquant()	mpaqt_quant_sc()
max_iter	100	100	100	100
tolerance	1e-4	1e-1	1e-4	1e-4
n_bins	50	50	from prequant	50
weight_update_start	internal prequant default	20	not applicable	internal prequant default
prior_start	25	not applicable	50	50
convergence_start	25	25	25	25
do_umi_correction	FALSE	not applicable	FALSE	TRUE
umi_correction_timing	"post"	not applicable	"post"	"post"
compute_uncertainty	FALSE	not applicable	FALSE	not exposed (FALSE)

mpaqt_quant() passes its public post-quantification defaults to mpaqt_postquant(), which is why its prior_start default differs from a direct call to mpaqt_postquant().

References

1. Bray, N. L., et al. (2016). Near-optimal probabilistic RNA-seq quantification. Nature Biotechnology.
2. Chen, Y., et al. (2021). Context-aware transcript quantification from long read RNA-seq data with bambu. Nature Methods.
3. Sigrist, F. (2022). Gaussian process boosting. Journal of Machine Learning Research.