PacBio physical endpoint model (2026-07-21)

Motivation

PacBio full-length libraries provide a physical expectation that an intact read should align close to both transcript ends. The prior empirical endpoint model learned a sample-specific 20x20 relative-gap histogram and regressed held-out PacBio accuracy, so PacBio auto disabled it. This experiment tests a constrained physical observation model instead.

Model

Each uniquely assigned alignment is generated by one of four latent states:

intact at both ends;
5'-truncated and 3'-intact;
5'-intact and 3'-truncated; or
broken/diffuse at both ends.

Intact gaps follow separate nucleotide-scale truncated exponential densities. Truncated gaps are diffuse over the physical transcript length. Mixture weights and the two intact tolerances are estimated by EM from unique reads, with a Dirichlet prior (80, 5, 5, 10) and 100 pseudo-observations anchoring each intact tolerance at 25 nt. Tolerances are constrained to 3--200 nt. Candidate likelihoods are normalized within each equivalence class and combined log-linearly with the existing logistic evidence.

The model is annotation-free beyond the transcript sequences and observed transcript-space alignments. It is opt-in as --coverage-ablation pacbio-physical-endpoint; the current default is not changed by this study.

Kinnex Figure 2

Model	D Pearson	D Spearman	D RMSE	E Pearson	E Spearman
no coverage	0.8273	0.8578	2.096	0.6396	0.6706
current auto	0.7530	0.8202	2.638	0.7015	0.7303
physical, weight 0.1	0.8301	0.8464	2.097	0.7060	0.7390
physical, weight 0.2	0.8385	0.8507	2.041	0.7067	0.7374
physical, weight 0.3	0.8396	0.8524	2.033	0.7068	0.7375

Weight 0.3 is a Pareto improvement over the current auto model and improves D Pearson/RMSE and E Pearson over no coverage. Across the six samples, the model learned approximately 88.5% intact, 9.6% 5'-truncated, 1.8% 3'-truncated, and 0.08% doubly broken mass. Both intact tolerances reached the conservative 3 nt lower bound, consistent with sharply defined PacBio transcript ends in these alignments.

LongBench PacBio 50k outer panel

At weight 0.3, changes relative to current PacBio auto across eight cell lines:

Metric	Mean change	Libraries improved
CCC	+0.00951	8/8
Spearman	+0.00029	6/8
RMSE	-73.48	8/8
MARD	-0.000207	8/8

Mean wall time rose from 2.05 to 3.47 seconds for a 50k library. Coverage-model time rose from 0.043 to 0.093 seconds; most of the total increase came from a longer EM (1.18 to 2.03 seconds). Peak RSS was unchanged within process noise (497.4 versus 493.7 MB).

PacBio 250k confirmations

Relative to the current no-endpoint PacBio auto kernel at weight 0.3:

Sample	CCC change	Spearman change	RMSE change	MARD change
H69	-0.00518	+0.00268	+17.07	-0.000699
H146	+0.00752	+0.00308	-84.83	-0.000687
H526	-0.00031	+0.00206	+0.99	-0.000569

All confirmations improve rank and relative deviation; calibration is mixed, with a meaningful H146 improvement, nearly neutral H526, and a modest H69 regression. Weights 0.1 and 0.2 did not remove the H69 regression and had worse mean CCC/RMSE, so weight 0.3 remains the best tested compromise.

Additional PacBio benchmarks

Two further truth-bearing PacBio datasets were rerun with the current release binary, SQUAREM, four threads, no alignment filters, and physical endpoint weight 0.3. This controls both model runs for code version and optimizer.

Independent human simulation (100k reads)

Model	Pearson	Spearman	CCC	RMSE	MARD
current auto	0.87278	0.70317	0.86562	0.26914	0.04233
physical endpoint	0.91557	0.72319	0.91220	0.21247	0.03123

The physical model improves every accuracy metric: RMSE falls by 21.1% and MARD by 26.2%. Wall time increases from 105.75 to 107.10 seconds. Coverage scoring increases from 0.110 to 0.210 seconds and EM from 4.94 to 5.82 seconds. Peak RSS does not increase (3068 versus 3053 MB; difference is process noise).

Public PacBio SIRV E0 FLNC

Model	RMSE	MARD	estimate CV
current auto	1.20682	0.42179	1.20682
physical endpoint	1.17425	0.39520	1.17425

E0 is equimolar, so Pearson, Spearman, and CCC are undefined and the evaluator reports them as zero. The physical model reduces RMSE by 2.7%, MARD by 6.3%, and unwanted estimate dispersion by 2.7%. Wall time is effectively unchanged (134.56 versus 134.81 seconds). Coverage scoring increases from 0.091 to 0.550 seconds, while EM is unchanged at 0.127 seconds and 52 evaluations. Peak RSS is unchanged within noise (147 versus 138 MB).

Calibration influence assessment

The three 250k LongBench comparisons were assessed with 1,000 paired gene-cluster bootstrap replicates. Clustering by gene preserves dependence among competing isoforms. The bootstrap resamples the paired baseline and physical estimates together and therefore measures uncertainty in their difference rather than uncertainty in either result alone.

Sample	Metric	Median delta	95% interval	P(improved)
H69	CCC	-0.00520	[-0.02500, 0.00669]	0.245
H69	RMSE	+16.91	[-60.04, 121.50]	0.366
H69	Spearman	+0.00267	[0.00161, 0.00367]	1.000
H69	MARD	-0.000699	[-0.000919, -0.000467]	1.000
H146	CCC	+0.00613	[-0.00972, 0.00897]	0.702
H146	RMSE	-84.18	[-197.67, 40.90]	0.726
H146	Spearman	+0.00301	[0.00185, 0.00432]	1.000
H146	MARD	-0.000681	[-0.000933, -0.000421]	1.000
H526	CCC	+0.00153	[-0.09599, 0.01459]	0.540
H526	RMSE	-4.70	[-218.31, 502.82]	0.537
H526	Spearman	+0.00204	[0.00102, 0.00307]	1.000
H526	MARD	-0.000572	[-0.000772, -0.000367]	1.000

Thus none of the CCC/RMSE changes is distinguishable from zero under the gene-cluster bootstrap, while the rank and relative-error improvements are consistent in every bootstrap replicate for all three samples.

Influence decomposition identified one common dominant outlier: ENST00000710941 (MALAT1-256), a 1,289-nt transcript with zero matched Illumina support. It receives little or no baseline abundance but substantial physical-model abundance. Its squared-error increase is 2.36 times the net H69 RMSE increase; improvements elsewhere cancel the remainder. In H526 it is 122 times the nearly-zero net RMSE increase. Removing only this transcript changes the model deltas as follows:

Sample	CCC delta	RMSE delta	Spearman delta	MARD delta
H69	+0.00171	-23.39	+0.00269	-0.000699
H146	+0.00837	-92.93	+0.00308	-0.000687
H526	+0.01201	-120.98	+0.00206	-0.000569

MALAT1 has many nested annotated isoforms. An alignment truncated relative to a long isoform can land close to both ends of a contained short isoform. The physical model then interprets geometrical endpoint agreement as evidence for the short isoform, even though the endpoint is not identifiable among the compatible transcript structures. This is a specific containment failure, not a reason to weaken endpoint evidence globally. The next refinement should gate endpoint evidence by within-equivalence-class endpoint identifiability and/or limit coverage-only creation of a previously unsupported nested isoform. This can be implemented from candidate transcript geometry without requiring a gene annotation.

The reproducible analysis is scripts/assess_pacbio_endpoint_calibration.py; raw JSON reports are under oarfish-evaluation-data/pacbio-physical-calibration-assessment-20260721.

Extreme nested-isoform creation guard

The influence result motivated a narrowly targeted safeguard. A transcript is eligible for protection only if it co-occurs in a read equivalence class with a compatible transcript at least 5 kb longer. The physical-model result is then compared with a full current PacBio auto/no-endpoint anchor. An increase is clamped only when all of the following hold:

the transcript has the extreme length-containment risk;
physical coverage creates more than 100 times its baseline abundance plus 10 counts per million; and
the resulting transcript would consume more than 0.25% of the library.

The anchor is used only for the post-fit guard and does not initialize the physical-model EM. Thus ordinary endpoint inference and its optimization path are unchanged. This uses only equivalence classes, transcript lengths, and sample-internal abundance estimates; no gene annotation is required.

PacBio 250k confirmation

Relative to current PacBio auto/no-endpoint, the guarded model improves every metric in all three samples:

Sample	CCC delta	Spearman delta	RMSE delta	MARD delta	Clamped transcripts
H69	+0.00066	+0.00269	-10.37	-0.000661	1
H146	+0.00752	+0.00308	-84.83	-0.000687	0
H526	+0.00907	+0.00206	-66.93	-0.000502	1

MALAT1-256 falls from 946 to 2.0 assigned reads in H69 and from 2,594 to 2.3 in H526. H146 is not clamped and retains the original improvement.

Preservation checks

On the eight-library LongBench 50k PacBio panel, mean CCC improves by 0.00163, Spearman by 0.00175, RMSE by 7.10, and MARD by 0.000394 relative to current auto. Spearman and MARD improve in 8/8 libraries, RMSE in 7/8, and CCC in 5/8. Zero or one transcript is clamped per library.
The independent 100k human simulation is unchanged to displayed precision: CCC 0.91220, Spearman 0.72318, RMSE 0.21247, and MARD 0.03123.
Public PacBio SIRV E0 is exactly unchanged: RMSE 1.17425 and MARD 0.39520. No transcript is protected or clamped.
Kinnex Figure 2 is unchanged: Panel-D Pearson 0.83963, Spearman 0.85258, and RMSE 2.03261; Panel-E Pearson 0.70683 and Spearman 0.73735. No SIRV is protected or clamped.

The cost is the extra baseline EM used to calibrate the guard. At 250k it adds 7.6--8.8 seconds; mean LongBench 50k wall time increases by 2.35 seconds over current auto. Peak memory remains effectively unchanged. Mapping-dominated raw-read runs change less: the independent simulation takes 110.4 seconds and public SIRV 134.7 seconds.

This is the first tested form that preserves the independent and Kinnex gains, eliminates the influential MALAT containment failure, and makes all four 250k metrics improve in all three confirmation samples. Final outputs are under oarfish-evaluation-data/pacbio-physical-extreme-creation-final-20260721*, pacbio-physical-extreme-creation-final-other-20260721, and kinnex-wtc11/pacbio-physical-extreme-creation-final-20260721.

EM optimizer comparison

The guarded physical model was compared with ordinary EM, SQUAREM, and DAAREM on H69 at 50k and 250k reads. The internal no-endpoint calibration anchor uses SQUAREM in all arms; this comparison varies the user-selected final optimizer.

Depth	Optimizer	Final EM seconds	Evaluations	Converged	Wall seconds
50k	ordinary	1.09	1000	no	3.66
50k	SQUAREM	1.64	1000	no	4.16
50k	DAAREM	6.61	1000	no	9.15
250k	ordinary	4.79	1000	no	13.56
250k	SQUAREM	5.73	1000	no	13.67
250k	DAAREM	11.61	1000	no	19.84

Ordinary EM and DAAREM produce effectively identical H69 estimates. SQUAREM ends at a slightly different point because all three hit the evaluation limit. Acceleration cannot help when it does not reduce evaluations; its more complex evaluations then add overhead. On Kinnex Day 0 replicate 1, all three converge at the mandatory minimum of 51 evaluations and are indistinguishable: final EM is 0.419--0.420 seconds and wall time is 3.94--3.97 seconds.

Therefore SQUAREM is not a runtime win for the physical endpoint model in the tested cases. Ordinary EM is the fastest or tied choice. DAAREM should not be recommended here. Raw outputs are under oarfish-evaluation-data/pacbio-physical-em-comparison-20260721.

Ordinary-EM calibration anchor

The internal PacBio no-endpoint calibration anchor was subsequently changed from SQUAREM to ordinary EM. The user-facing final optimizer already defaults to ordinary EM. In a controlled H69 comparison where the final optimizer was ordinary in both arms:

Depth	Anchor	Coverage/anchor phase	Final EM	Wall time
50k	SQUAREM	1.746 s	1.092 s	3.66 s
50k	ordinary	1.221 s	1.093 s	3.14 s
250k	SQUAREM	7.307 s	4.793 s	13.56 s
250k	ordinary	5.978 s	5.568 s	13.02 s

The 250k final-EM timing difference is run-to-run noise and masks part of the 1.33-second anchor saving. Quantification is unchanged to practical precision: H69 CCC changes by less than 1.3e-7, Spearman by less than 8.1e-7, RMSE by less than 0.002, and MARD by less than 4e-9 at both depths.

With ordinary EM in both stages, the full 250k panel still improves all four metrics in 3/3 samples (mean CCC +0.00578, Spearman +0.00288, RMSE -54.62, MARD -0.000473). The LongBench 50k panel retains mean CCC +0.00164, Spearman +0.00205, RMSE -7.07, and MARD -0.000328; Spearman and MARD improve in 8/8 libraries and RMSE in 7/8. Outputs are under oarfish-evaluation-data/pacbio-physical-ordinary-anchor-final-*.

Parallel physical EM

The original parallel M-step used a shared atomic floating-point accumulator for every transcript. Enabling it at four threads was counterproductive: atomic contention approximately doubled EM time. Following the current Salmon EM layout, the M-step now partitions reads into one contiguous shard per worker, accumulates into private transcript vectors, and reduces those vectors in parallel over transcripts. The private vectors are allocated once and reused across fixed-point evaluations. Invariant alignment-score, coverage, and KDE terms are also packed once before inference instead of being recomputed twice per alignment per iteration. Singleton reads bypass the denominator because their only candidate necessarily receives one count. See Salmon's packed.rs and lib.rs.

To avoid parallel setup overhead, the physical anchor and final EM use this path only with more than one requested thread and at least 100,000 aligned reads; smaller libraries remain serial.

At four threads on the 250k confirmation panel:

Sample	Serial wall	Parallel wall	Saving	Serial anchor	Parallel anchor	Serial final EM	Parallel final EM
H69	13.02	7.66	41%	5.98	3.30	5.57	2.79
H146	11.66	7.15	39%	5.01	3.02	5.26	2.70
H526	13.93	7.92	43%	6.12	3.34	6.19	2.97

Counts and all reported accuracy metrics are identical to displayed precision, and peak RSS is unchanged within measurement noise. The 50k threshold smoke test remains serial (3.06 versus 3.14 seconds in the prior run). Raw outputs are under oarfish-evaluation-data/pacbio-physical-adaptive-parallel-final-*.

The Salmon-style reusable shards and packed invariant weights improve the four-thread result further:

Sample	First parallel wall	Packed wall	Additional saving	Packed anchor	Packed final EM
H69	7.66	5.25	31%	2.01	1.60
H146	7.15	4.93	31%	1.96	1.54
H526	7.92	5.35	32%	2.11	1.66

Relative to the original serial implementation, total wall time is down 58--62%. Peak RSS remains approximately 494 MB: the packed weight vector is offset by replacing three fat slice pointers per read with one alignment slice and a weight offset. Accuracy is unchanged apart from floating-point reduction order. The 50k serial smoke test remains 3.02 seconds with identical metrics.

H69 thread scaling after packing shows diminishing returns:

Threads	Wall	Anchor	Final EM	Peak RSS
1 (pre-parallel serial)	13.02	5.98	5.57	~494 MB
2	7.11	2.67	2.21	490 MB
4	5.25	2.01	1.60	494 MB
8	4.99	2.07	1.44	492 MB
16	4.93	1.84	1.47	505 MB

Four threads retain the best efficiency; eight threads provide a small wall time improvement, while sixteen provide essentially no further gain. Final raw outputs are under oarfish-evaluation-data/pacbio-physical-salmon-parallel-final-* and the pacbio-physical-salmon-scaling-* directories.

The same packed invariant representation is used by ordinary serial inference, not only by the parallel M-step. On the H69 50k threshold case this reduces the physical run from 3.02 to 2.26 seconds; anchor EM falls from 1.14 to 0.78 seconds and final EM from 1.06 to 0.65 seconds. Accuracy is unchanged. A regular coverage-free serial run takes 1.51 seconds total, including 0.66 seconds of EM. Bootstrap replicates retain their sampling-specific iterator and are not yet packed; they are already parallelized across replicates. Raw serial results are under oarfish-evaluation-data/serial-packed-h69-pb-50k-20260721.

Promotion to the PacBio auto default

After the nested-isoform guard and packed-EM validation, PacBio and PacBio HiFi now select the guarded physical endpoint kernel when users request --coverage-model auto with the default full configuration. The three-sample 250k confirmation reproduces the explicit physical-model estimates and guard decisions: H69 and H526 each clamp the same one extreme nested transcript, H146 clamps none, and accuracy metrics match the validated opt-in runs. Four-thread wall times are 4.98, 4.91, and 5.29 seconds respectively. ONT cDNA and dRNA smoke tests retain their prior adaptive and competing-risks kernels. Outputs are under oarfish-evaluation-data/pacbio-physical-auto-default-250k-20260721 and coverage-auto-default-nonpacbio-smoke-20260721.

Decision

The guarded physical model is now the PacBio kernel selected by --coverage-model auto. It is the first endpoint candidate to improve Kinnex D and E simultaneously, improves all four mean LongBench 50k metrics, and improves every usable metric in two additional PacBio datasets. The evidence generalizes across human simulation, real human LongBench libraries, Kinnex SIRV mixtures, and an independent public SIRV FLNC library. The original default blocker was rare coverage-only creation of endpoint-unidentifiable nested isoforms, dominated by MALAT1-256; the abundance-anchor guard removes that failure while preserving the gains. Explicit ablation selections remain available for controlled comparisons.

Raw outputs are under oarfish-evaluation-data/kinnex-wtc11/pacbio-physical-endpoint-20260721, pacbio-physical-longbench-20260721, and pacbio-physical-confirm-250k-20260721*. The additional controlled runs are in oarfish-evaluation-data/pacbio-physical-other-20260721.