mixqtl branch: run on current NumPy/PyTorch and match mixqtl.R nominal output - #210
Open
fbyukgo wants to merge 1 commit into
Open
mixqtl branch: run on current NumPy/PyTorch and match mixqtl.R nominal output#210fbyukgo wants to merge 1 commit into
fbyukgo wants to merge 1 commit into
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
We've been using the mixqtl branch on our data (~250 donors) and made a small set of fixes so it runs on current NumPy/PyTorch and reproduces mixqtl.R's nominal output. Sharing them back in case they're useful.
What's here:
Dependency updates (P1–P2): NumPy 2.0 (np.NaN→np.nan, np.float→float) and PyTorch 2.0 (torch.solve/qr/symeig → their torch.linalg equivalents). No behaviour change — just needed to import/run.
mixqtl.R parity in map_nominal/trc_calc (P4–P7, these change results):
NaN-initialise the meta/asc arrays — ASE-absent genes were emitting uninitialised memory as beta_meta (a bug).
Meta rule follows rlib_meta.R: trc-fallback when ASE is absent, and the n_cutoff gate (branch had neither).
trc_calc fits the covariate offset on all count-observed individuals before the count filter (as R does), instead of re-fitting on the survivor subset — the old order is unstable/biased on low-expression genes. (Also fixes an SE bug where the phenotype shape gave dof = -1.)
p-value convention matched to R (z above n_cutoff, else t with df = n).
Two of these (NaN-init, and the covariate order) are genuine correctness fixes, not just cosmetic alignment.
Validation. With trc_cutoff matched to R, beta_trc/beta_asc/beta_meta agree to ~4 dp (r ≈ 0.9997, ~0.02% sign flips, scale ratios ~1.0, p-values match after P7), and it runs ~100× faster than the R split job. The concordance script we used is included as a regression check.
One note: count_threshold is a call-site literal in map_nominal; might be worth exposing as a parameter (we left the existing default untouched).
Happy to split this differently, and — the branch is nominal-only, so if there's interest we'd be glad to contribute our permutation-based eGene layer for the ASE-integrated statistic as a separate PR.
Full patch descriptions:
P1 (NumPy 2): np.NaN→np.nan (genotypeio.py, mixqtl.py, eigenmt.py); np.float(→float( (eigenmt.py).
P2 (PyTorch 2): torch.solve(B,A)→torch.linalg.solve(A,B); torch.qr→torch.linalg.qr; torch.symeig(…,eigenvectors=False)→torch.linalg.eigvalsh (mixqtl.py, core.py, eigenmt.py). Equivalent numerics.
P4 (bug): meta/asc arrays were np.empty and only written when ASE existed → ASE-absent genes emitted uninitialised memory as beta_meta. Now np.full(…, nan) (and -1 for int columns).
P5 (meta rule = rlib_meta.R::my_meta_fast_): start from the larger-n arm; per-SNP fall back to the other arm where the primary is NaN; inverse-variance combine only when both arms have ≥ n_cutoff (15) samples; when ASE is absent, meta = trc. Branch previously combined only when ASE existed, with no n-gate and no fallback.
P6 (bug / sounder estimator): fit the covariate offset on all count-observed individuals before the count filter, subtract it, then drop low-count individuals, then regress on genotype with no covariate re-fit (as mixqtl.R). Old order re-fit covariates on the survivor subset → rank-deficient when survivors < #covariates and biased (survivors selected on expression). Also: pass phenotype 1-D so dof = n-2 (was giving -1 → NaN SE).
P7 (p-value = get_pval_fast_): z when the relevant n > n_cutoff, else Student-t with df = n; meta's effective n tracked per row. Closes the last -log10 p gap on small-n genes.