[xyg123]

Summary

Adds a biosample_enrichment step to pts. For each GWAS study and each tissue or cell type in a specificity panel, it tests whether the genes L2G prioritises are more specific to that biosample than gene length and gene density explain.

There is one linear regression per study × biosample, fitted over protein-coding genes:

• Outcome: the gene's highest L2G score in the study (0 if unscored).
• Predictor: the gene's specificity_score in the biosample.
• Covariates: log10 gene length and the number of TSSs within ±1 Mb.
• Standard errors: genes in the same 1 Mb region are counted as one.
• Test: one-sided (enrichment only), with Bonferroni and BH corrections across each study's biosamples.

Optionally, studies with an identical diseaseIds set can be pooled (pool_by_disease, off by default).

The code is split in two:

• pts/src/pts/transformers/biosample_enrichment.py handles the data (loading, joining, shaping, writing).
• pts/src/pts/transformers/enrichment/regression.py does only the numbers (numpy arrays in, numbers out).

Main:

1. biosample_enrichment() is called by the transform task from pts/config.yaml. It reads the settings: panels, pooling, kind of standard error and window sizes.

2. Loading and preparing inputs (biosample_enrichment.py). All inputs are read lazily with scan_dataset.

Function	Input	Produces
gene_covariates()	target	Per gene: biotype, chromosome, TSS, log10 length, TSS density within ±1 Mb
gwas_studies()	study	GWAS studies only, each assigned to a unit: the study itself, or its disease-set group when pooling
study_gene_scores()	l2g_prediction + credible_set	The highest L2G per study and gene, linked through studyLocusId
(inline)	—	The highest L2G per unit and gene, and which studies belong to each unit (for nStudiesPooled)

3. Per panel.

• specificity_matrix() pivots the panel's baseline_expression rows into a gene × biosample table.
• score_panel() then:
  • restricts to protein-coding genes with specificity in the panel;
  • standardises the covariates and assigns each gene to a 1 Mb region;
  • builds a regression.Design once for the panel;
  • fits units in batches of 512, which keeps memory to about 80 MB per batch.

4. The regression (regression.py).

• Design.__init__ runs once per panel. It removes from each biosample's specificity what length and density explain. By Frisch-Waugh-Lovell, this gives the same slope as the full model and lets every study reuse the adjustment. It also precomputes per-region totals.
• Design.fit runs once per batch. The slope for the whole batch is one matrix product. The standard error comes from _noise_by_region (default), _noise_by_gene (HC0) or the classical formula.
• upper_tail_p, bonferroni and benjamini_hochberg turn z-scores into p-values and corrected p-values.

5. Writing. Each unit's results are expanded back to its member studies. The panels are concatenated, cast to OUTPUT_SCHEMA, sorted, and written with write_dataset to output/biosample_enrichment: one row per study × biosample.

Other changes

• pts/config.yaml: the new step, with a GTEx tissue panel.
• orchestration/.../unified_pipeline.yaml: a pts_biosample_enrichment node, also added to pts_release_metrics.depends_on. It is published but no other step reads it.
• croissant: the distribution entry and the biosample_enrichment record set.
• pts/pyproject.toml: statsmodels added to the test group only, as the reference for the tests.

Validation

• Tests: 25 new tests.
• Matches the prototype: on the same CELLEX specificity, the port matches the gentropy notebook prototype.
• Release specificity vs CELLEX: the release baseline_expression specificity gives a per-study Spearman of 0.985 (median) against CELLEX.
• Full release: run on platform-2610-2 outputs. 58,326 studies × 50 GTEx tissues gives 2.9M rows in about 2 min with a peak RSS of 1.9 GB on one machine. For IBD (GCST90292538), the top 5 tissues are the same as the prototype's.

In one picture

config.yaml ──> biosample_enrichment()
                 ├─ gene_covariates(target)                 → length, density, positions
                 ├─ gwas_studies(study)                     → study → unit (pooling)
                 ├─ study_gene_scores(l2g, credible_set)    → y: max L2G per study/unit × gene
                 └─ for each panel:
                      specificity_matrix(baseline_expression) → x: gene × biosample
                      score_panel()
                        ├─ Design(covariates, x, regions)   → adjust x once     [regression.py]
                        └─ per batch of 512 units:
                             Design.fit(y)                  → slope + SE        [regression.py]
                               └─ _noise_by_region / _noise_by_gene
                             z, upper_tail_p, BH, Bonferroni                    [regression.py]
                 └─ write_dataset() → output/biosample_enrichment

[project-defiant]

The L2G already requires cs, study & target upstream, no need to add these, will make the dag render more complex then needed.