Skip to content

stats: wire verbs

The facts on this page are the same facts mzLibRust and mzLibR will publish, because all three render them from one spec per wire verb. The parameter, result, error and caveat tables are generated; edit the spec, not this page. How the reference facts are generated explains the mechanism. For the Python signatures and result classes, see the API reference. For worked tasks, see the Differential abundance guide.

Python spells stdin as a list argument. Each verb reads its list input on stdin, and pyMzLib writes it for you: stats fit's stdin is coefficients=, stats adjust's is p_values= (with None for an untested entry) and stats meta's is studies=. The stdin rows below describe the wire.

Performance: each call starts one bridge process, about 120 ms before mzLib does any work. stats fit then fits every feature in one pass: a 10,000-feature by 24-sample table takes about 0.6 s on one thread and about 0.4 s with threads=-1. Test every coefficient you need in one call, since the model is fitted once and moderated per coefficient. Every number is mzLib's; the bridge parses the tables and shapes the result, and pyMzLib does no arithmetic.

fit

Fit one linear model per feature of a feature-by-sample table and test each named coefficient with limma's empirical-Bayes moderated t, Benjamini-Hochberg adjusted, in one call.

Wire verb stats fit · Python pymzlib.stats.fit

Wraps

Parameters

Parameter Type Default Unit Range Meaning
responses path required — — A TSV with a header row. The first column is the feature id (unique, any header); every other column is one sample, named by its header. Cells are used as given, so log-transform intensities first. A blank cell, NA or NaN (any case) is missing and is omitted from that feature's fit, never imputed. A 0 is an OBSERVATION. Blank lines are skipped.
design path required — — A TSV with a header row and one row per sample. The first column is the sample name and must match the responses header exactly (rows are matched by name, so their order does not matter); every other column is one coefficient, named by its header, and every cell must be a finite number. The design is used as written: include an intercept column if you want one.
stdin string[] required — >= 1 line; each a design column name, none repeated The coefficients to test, one design column name per line. Each gets its own moderated test and its own Benjamini-Hochberg family. The model always fits every design column.
trend flag false — — Let the prior variance trend with each feature's average response (limma's trend = TRUE), as a natural cubic spline.
spline_basis (wire --spline-basis) int 4 basis functions >= 1; only with --trend Basis functions of the trend spline, intercept included. mzLib caps it at the distinct average responses and the fitted features minus one; prior.spline_basis_count reports what was used.
threads int 1 threads >= 1, or -1 for every core Features fitted at once. The bridge default is 1 (BULK.md); mzLib's own default is -1. Results are identical at any value.

Returns

Top-level fields:

Field Type Unit Null? Meaning
responses_file string — never Absolute path of the responses table.
design_file string — never Absolute path of the design table.
feature_count int features never Feature rows read.
sample_count int samples never Sample columns read.
sample_names string[] — never In responses header order, which is the order the design rows were put in.
coefficient_names string[] — never Every design column, in order.
tested_coefficients string[] — never The stdin coefficients, in order; the table holds one block of feature_count rows per entry.
trend bool — never Whether --trend was given.
threads int threads never The thread count used.
missing_count int cells never Response cells that were blank, NA or NaN.
zero_count int cells never Response cells that are exactly 0, fitted as observations.
status_counts map<string,int> features never Features per status: fitted, too_few_observations, rank_deficient.
residual_df_differ bool — never True when the fitted features do not all have the same residual df (missing values omitted). Then default limma (>= 3.61, eBayes legacy = FALSE) would use a different prior estimator and report different moderated statistics; these follow legacy = TRUE. A caveat says so.
prior object — never The fitted variance prior, shared by every tested coefficient; fields under result.tables.prior.
row_count int rows never feature_count x the number of tested coefficients.
column_names string[] — never The table's columns, in order.
columns table — never One row per (tested coefficient, feature): coefficient blocks in stdin order, features in file order.
caveats string[] — never Three static caveats, plus one each when residual_df_differ, when any feature was not fitted (with counts), when any cell is 0, and when the prior df is infinite.

Per-row fields (one value per record, in columns or each list entry):

Field Type Unit Null? Meaning
feature string — never The feature id from the responses file's first column.
coefficient string — never The coefficient this row tests.
status string — never fitted | too_few_observations (fewer observed samples than coefficients plus one) | rank_deficient (the design restricted to this feature's observed samples is not of full rank). Only fitted rows carry statistics or enter the BH family.
estimate float response units per coefficient unit (e.g. log2 fold change for a 0/1 group column on log2 intensities) yes: status is not fitted The least-squares coefficient; moderation does not change it.
standard_error float response units per coefficient unit yes: not fitted Moderated standard error: sqrt(posterior_variance) x the coefficient's unscaled standard deviation.
t float standard errors yes: not fitted Moderated t-statistic.
df_total float degrees of freedom yes: not fitted Degrees of freedom of the moderated t: df_residual + prior df, capped at the pooled residual df.
p_value float fraction (0 to 1) yes: not fitted Two-sided p-value.
bh_adjusted float fraction (0 to 1) yes: not fitted (not in the family) Benjamini-Hochberg adjusted p-value over this coefficient's fitted features. Not a target-decoy q-value.
posterior_variance float squared response units yes: not fitted The moderated residual variance (d0 s0^2 + d s^2) / (d0 + d).
prior_variance float squared response units yes: not fitted (the feature did not enter the prior) s0^2 for this feature: one value for every feature without --trend, varying with average_response with it.
sigma float response units yes: not fitted Residual standard deviation, sqrt(RSS / df_residual), before moderation.
df_residual int — yes: not fitted Observed samples minus coefficients.
observed int samples never Samples with a finite response for this feature.
average_response float response units yes: the feature has no finite response at all Mean of the feature's observed responses: the covariate --trend fits on.

Fields of prior:

Field Type Unit Null? Meaning
df float degrees of freedom yes: infinite: df_infinite is true Prior degrees of freedom d0.
df_infinite bool — never True when the residual variances are no more dispersed than sampling alone explains, so every feature takes the prior variance.
trended bool — never Whether s0^2 varies with average_response.
spline_basis_count int basis functions never Basis functions the trend used, intercept included; 1 without a trend.
scale float squared response units yes: trended: s0^2 is per feature, in the prior_variance column s0^2, the prior variance.

Errors

Kind pyMzLib raises When
usage UsageError responses or design missing, blank or not found; a table empty, header-only or with too few columns; a row whose field count differs from its header; an empty or repeated feature id, sample or coefficient name; a response cell that is not a number (other than blank, NA, NaN); a design cell that is not a finite number
usage UsageError the design's samples are not exactly the responses' sample columns (both differences are named); no coefficient on stdin, one named twice, or one not in the design; --spline-basis without --trend, or below 1; --threads 0 or below -1
usage UsageError mzLib refuses the input: the design is not of full column rank, or fewer than 2 features could be fitted, so no variance prior exists (mzLib's ArgumentException, reclassified)
correctness BridgeError never expected
service_unavailable ServiceUnavailableError never raised by this verb

Caveats

  • This is limma's LEGACY empirical-Bayes estimator (eBayes(legacy = TRUE)). Since limma 3.61, eBayes defaults to a different estimator whenever residual df differ between features, which omitting missing values causes; residual_df_differ says when that applies.
  • Not implemented: robust = TRUE, contrasts, the B-statistic, observation weights. Write a contrast as its own design column instead.
  • A 0 is an observation. Log-transform intensities first, and blank any cell that means 'not measured'.
  • bh_adjusted adjusts for the false discovery rate among the features tested. It is not a target-decoy q-value for identifications.

Example

fit = pymzlib.stats.fit('malat_dilution_log2_vs_fluc.tsv', 'malat_dilution_design.tsv', ['malat_250ng', 'malat_125ng'])

Recorded from the live bridge as stats_fit_malat.json. The doctests in the API reference replay these recordings in CI.

fit = pymzlib.stats.fit('limma_reference_responses.tsv', 'limma_reference_design.tsv', ['age_decades'])

Recorded from the live bridge as stats_fit_limma.json. The doctests in the API reference replay these recordings in CI.

Same verb in other bindings

Python (pyMzLib) Rust (mzLibRust) R (mzLibR)
pymzlib.stats.fit mzlib::stats::fit_with with FitOptions stats_fit

References

Since

Wire protocol pyMzLib mzLibRust mzLibR
1 0.3.0 not yet shipped not yet shipped

Not yet verified

  • Rust and R spellings are proposals until those bindings ship the verb.
  • stats_fit_malat.json is the largest recorded fixture (about 258 KB): 492 features x 2 coefficients. It is real data (mzLib's RNA oligo table, #1388), derived by a rule held in StatisticsTests.MalatDilutionFixture_IsLog2OfMzLibsOligoTable.

adjust

Adjust a list of p-values for the false discovery rate with Benjamini-Hochberg, keeping every input position and leaving untested entries out of the family.

Wire verb stats adjust · Python pymzlib.stats.adjust

Wraps

Parameters

Parameter Type Default Unit Range Meaning
stdin string[] required — >= 1 line; each in [0, 1], blank, NA or NaN One p-value per line. BLANK LINES ARE KEPT (as in sdrf parse-age) so output row i is line i; one trailing newline does not add a line. A blank, NA or NaN line is a feature that was not tested: it stays null and is not counted in m. A UTF-8 BOM is dropped.

Returns

Top-level fields:

Field Type Unit Null? Meaning
line_count int lines never Lines read from stdin; the length of every column.
tested_count int p-values never m: the finite p-values, the family the adjustment is over.
row_count int rows never Equal to line_count.
column_names string[] — never The table's columns.
columns table — never One row per input line, in order.
caveats string[] — never Two static caveats.

Per-row fields (one value per record, in columns or each list entry):

Field Type Unit Null? Meaning
p_value float fraction (0 to 1) yes: the line was blank, NA or NaN: not tested The p-value as read.
bh_adjusted float fraction (0 to 1) yes: not tested min over j >= i of p_(j) m / j, capped at 1, for the i-th smallest p-value.

Errors

Kind pyMzLib raises When
usage UsageError stdin empty; a line that is not a number (other than blank, NA, NaN); a p-value outside [0, 1] (mzLib's refusal, reclassified)
correctness BridgeError never expected
service_unavailable ServiceUnavailableError never raised by this verb

Caveats

  • Benjamini-Hochberg controls the false discovery rate among the features tested. It is not a target-decoy q-value for identifications.
  • m is the number of finite p-values. Leave a line blank, or write NA, for a feature that was not tested rather than writing 1, which would enlarge the family.
  • stats fit already reports bh_adjusted per coefficient; use this verb for p-values from anywhere else.

Example

adjusted = pymzlib.stats.adjust([0.0002, 0.004, 0.019, None, 0.031, 0.2, float('nan'), 0.74])

Recorded from the live bridge as stats_adjust.json. The doctests in the API reference replay these recordings in CI.

Same verb in other bindings

Python (pyMzLib) Rust (mzLibRust) R (mzLibR)
pymzlib.stats.adjust mzlib::stats::adjust stats_adjust

References

Since

Wire protocol pyMzLib mzLibRust mzLibR
1 0.3.0 not yet shipped not yet shipped

Not yet verified

  • Rust and R spellings are proposals until those bindings ship the verb.

meta

Pool one effect size per study into a random-effects estimate per feature (DerSimonian-Laird), with heterogeneity, direction agreement and the largest leave-one-out change.

Wire verb stats meta · Python pymzlib.stats.meta

Wraps

Parameters

Parameter Type Default Unit Range Meaning
stdin string[] required — >= 1 line One study per line, tab-separated: feature, estimate, standard_error. No header. A feature's lines need not be adjacent; features are reported in order of first appearance. The estimate must be finite and the standard error finite and positive. Blank lines are skipped.
confidence float 0.95 fraction (0 to 1) (0, 1) Coverage of the reported interval, as a fraction.

Returns

Top-level fields:

Field Type Unit Null? Meaning
study_count int studies never Study lines read, over every feature.
feature_count int features never Distinct features; the table's row count.
confidence float fraction (0 to 1) never The coverage used, as a fraction.
row_count int rows never Equal to feature_count.
column_names string[] — never The table's columns.
columns table — never One row per feature.
caveats string[] — never Three static caveats.

Per-row fields (one value per record, in columns or each list entry):

Field Type Unit Null? Meaning
feature string — never As given on stdin.
studies int studies never Studies pooled for this feature.
estimate float units of the study estimates never Random-effects pooled estimate.
standard_error float units of the study estimates never Standard error of the pooled estimate.
confidence_low float units of the study estimates never Lower bound of the interval: estimate - z standard_error.
confidence_high float units of the study estimates never Upper bound of the interval.
p_value float fraction (0 to 1) never Two-sided p-value of the pooled estimate against 0, normal reference. Not adjusted across features.
tau2 float squared units of the study estimates never Between-study variance (DerSimonian-Laird); 0 for one study, and exactly 0 whenever q <= studies - 1.
q float chi-square statistic (studies - 1 degrees of freedom) yes: one study: no heterogeneity to measure Cochran's Q.
i_squared float fraction (0 to 1) yes: one study Share of the variation attributed to heterogeneity, as a fraction in [0, 1].
direction_agree int studies never Studies whose estimate has the pooled estimate's sign (a 0 estimate counts for neither).
direction_disagree int studies never Studies with the opposite sign.
leave_one_out_max_delta float units of the study estimates yes: one study Largest absolute change in the pooled estimate from dropping any one study.

Errors

Kind pyMzLib raises When
usage UsageError stdin empty; a line without exactly 3 tab-separated fields; an empty feature name; an estimate that is not a finite number; a standard error that is not finite and positive (mzLib's refusal, reclassified, naming the feature); --confidence not strictly between 0 and 1, or given without a value
correctness BridgeError never expected
service_unavailable ServiceUnavailableError never raised by this verb

Caveats

  • DerSimonian-Laird with a normal reference: the interval and p-value do not widen for few studies (no Hartung-Knapp adjustment), so treat a pooled estimate from two or three studies with care.
  • A feature with one study is reported with tau2 = 0 and q, i_squared and leave_one_out_max_delta null.
  • Each feature is pooled on its own and p_value is not adjusted across features; send the p_value column to stats adjust.

Example

pooled = pymzlib.stats.meta(studies)

Recorded from the live bridge as stats_meta_metafor.json. The doctests in the API reference replay these recordings in CI.

Same verb in other bindings

Python (pyMzLib) Rust (mzLibRust) R (mzLibR)
pymzlib.stats.meta mzlib::stats::meta_with with MetaOptions stats_meta

References

Since

Wire protocol pyMzLib mzLibRust mzLibR
1 0.3.0 not yet shipped not yet shipped

Not yet verified

  • Rust and R spellings are proposals until those bindings ship the verb.
  • The example's studies are mzLib's metafor reference cases (Test/StatisticalModels/ReferenceData/dl_inputs.tsv), copied verbatim to pkg/python/tests/fixtures/stats/metafor_dl_inputs.tsv. They are constructed cases, not proteomics data.