Metadata-Version: 2.4
Name: gwas-sumstats-parser
Version: 0.5.0
Summary: Vocabulary, parsing and standardisation for GWAS summary-statistics files
Author-email: Chris Finan <c.finan@ucl.ac.uk>
License-Expression: GPL-3.0-or-later
Project-URL: Homepage, https://cfinan.gitlab.io/gwas-sumstats-parser
Project-URL: Repository, https://gitlab.com/cfinan/gwas-sumstats-parser
Classifier: Programming Language :: Python
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: Implementation :: CPython
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Requires-Python: <3.14,>=3.11
Description-Content-Type: text/markdown
License-File: LICENSE.txt
Requires-Dist: numpy
Requires-Dist: scipy
Requires-Dist: mpmath
Provides-Extra: io
Requires-Dist: rapidgzip; extra == "io"
Requires-Dist: pysam>=0.17; extra == "io"
Provides-Extra: frame
Requires-Dist: pandas; extra == "frame"
Requires-Dist: tqdm; extra == "frame"
Provides-Extra: dev
Requires-Dist: pytest; extra == "dev"
Requires-Dist: pytest-dependency; extra == "dev"
Requires-Dist: flake8; extra == "dev"
Requires-Dist: bump2version; extra == "dev"
Requires-Dist: build; extra == "dev"
Requires-Dist: twine; extra == "dev"
Requires-Dist: ruff; extra == "dev"
Requires-Dist: mypy; extra == "dev"
Requires-Dist: pandas; extra == "dev"
Requires-Dist: tqdm; extra == "dev"
Dynamic: license-file

# gwas-sumstats-parser

__version__: `0.5.0`

gwas-sumstats-parser is a small Python library that reads GWAS
summary-statistics flat files into a standard column vocabulary, parses and
validates the values, and yields mutable scalar list rows. It is designed as
the shared ingest layer for
[merit](https://gitlab.com/cfinan/merit) and
[gwas-norm](https://gitlab.com/cfinan/gwas-norm): both consume it, neither
depends on the other.

It provides:

1. **A column vocabulary.** The canonical names for the quantities a
   summary-statistics file can carry (chromosome, position, alleles, effect
   size, standard error, p-value, allele frequency, sample sizes and so on).
   Callers explicitly map these names to source headers; aliases are not guessed.
2. **Value parsers.** Chromosome naming, allele strings, missing-value tokens
   and numeric fields, each
   with a stated rule for what is accepted, what is repaired and what is
   rejected.
3. **Scalar list-row reading.** The current reader handles directly supplied
   linear and ratio effects and retains probability and its negative base-ten
   logarithm. Sample counts can be derived from case and control counts or
   supplied globally. Allele frequencies can be derived from diploid allele
   counts or minor-allele information. Probability conversion and missing-value
   repair are supported, along with confidence intervals, t statistics,
   quantitative scores, approximate log odds and calibrated direction effects.
   No array bindings, pandas dependency or compiled parser are used.
4. **Explicit row failures.** The base reader raises a shared error with the
   original value and source row. The dropping reader skips bad rows and
   exposes failure records and counts grouped by failing method.

## Opening input files

The convenience opener detects compression from file contents, so a plain
file ending in `.gz` remains plain and gzip data ending in `.txt` is decoded.
It supports plain text and gzip, including block gzip (bgzip), and rejects
zstd and lz4 with named errors.

```python
from gwas_sumstats_parser import Reader
from gwas_sumstats_parser.opening import open_source

mapping = {
    "chr_name": "CHR", "start_pos": "POS",
    "effect_allele": "A1", "other_allele": "A2",
    "effect_size": "BETA", "standard_error": "SE", "pvalue": "P",
}
with open_source("summary-statistics.tsv.gz") as source:
    reader = Reader(source, mapping, effect_type="beta")
    for row in reader:
        print(row)
```

Automatic decoding (`codec="auto"`) uses rapidgzip when installed and falls
back to standard-library gzip otherwise. Select `codec="gzip"` to use the
standard library or `codec="rapidgzip"` to require the faster decoder. The
latter raises an installation error if unavailable and the file is gzip.
Set `encoding` when the file uses an encoding other than UTF-8.

Set `threads` to give rapidgzip a thread count; leaving it unset lets
rapidgzip choose. One rule covers decoder settings: a setting is ignored
whenever its decoder does not read the file. For plain text and indexed
input, `threads` and `codec="rapidgzip"` are ignored without any message,
because nothing was lost. When standard-library gzip reads a gzip file and
`threads` was given, the opener issues an `IgnoredSettingWarning`, because
the expected speed-up did not happen. Check `source.codec` to see which
decoder was used, and silence the warning with:

```python
import warnings
from gwas_sumstats_parser.opening import IgnoredSettingWarning

warnings.simplefilter("ignore", IgnoredSettingWarning)
```

A neighbouring `.tbi` or `.csi` genomic index selects pysam, taking precedence
over the gzip codec choice. Indexed input must use bgzip; ordinary gzip next
to an index is refused. Missing pysam or an unusable index raises an error;
it never silently starts an unindexed scan. The source exposes `indexed`,
`codec`, and the indexed chromosome names in `contigs`. Its low-level
`fetch(reference, start, end)` returns raw records with tabix's zero-based,
half-open boundaries and exact chromosome names. Full iteration includes
the file's header. Caller-owned streams have no implicit index access.

### Reading genomic regions

The reader's `fetch` method yields the rows inside zero or more requested
regions. A region is a chromosome, a start and an end. **Start and end are
one-based and both are included**, on every kind of source. With an indexed
source the reader itself subtracts one from the start for tabix, so callers
never adapt coordinates.

```python
with open_source("summary-statistics.tsv.gz") as source:
    reader = Reader(source, mapping, effect_type="beta")
    for row in reader.fetch(("chr1", 100, 200), ("2", 500, 600)):
        print(row)
```

- The requested chromosome is normalised exactly as row chromosomes are,
  including any `chromosome_synonyms`, so `chr1`, `CHR1`, `01` and `1` request
  the same rows, whatever spelling the file or its index uses.
- `fetch()` with no regions is plain iteration: same rows, same repair
  counts, same drop counts.
- The order and overlap of the requests never change the result. The reader
  sorts the requests into file order and merges requests on one chromosome
  that overlap or touch, so each matching row comes back once, in file
  order. An indexed source is asked in the order its index lists the
  chromosomes. Both kinds of source therefore return the same rows in the
  same order, provided the file is sorted the way an index requires.
- An indexed source serves regions from its index and can be queried as
  often as needed. An unindexed source is scanned. The scan looks up each
  line's chromosome first, reads the position only for a requested
  chromosome, and skips a line placed outside every region without parsing
  or reporting it. A line whose chromosome or position cannot be read is
  parsed in full, so its failure is reported.
- A data line with fewer cells than the header, including a blank line, is
  a row failure like any other. The one exception follows the rule above:
  during a region query, a short line that still has its chromosome and
  position cells, and is placed outside every region by them, is skipped.
  The reader then issues a `SkippedShortLineWarning`, importable from
  `gwas_sumstats_parser.reader`, naming the line and its cell count. A line
  with more cells than the header is read normally.
- A scan can be repeated. When a scan has finished, by reaching the end of
  the file or by stopping on an error, the next scan first returns the
  stream to the position it held when the reader received it and skips the
  header again, so row numbers restart at one and any comment lines you
  skipped stay skipped. A loop you leave early with `break` is not finished:
  the next loop carries on from the following line, as a Python file does.
  A stream that cannot go back, such as standard input, serves one finished
  scan; asking for another raises `ParseError`. To let the reader go back,
  skip leading comments with `readline()`, not with `next()`, because Python
  text files stop reporting their position after `next()`.
- A malformed request (positions that are not whole numbers, a start below
  one, an end before its start, an unusable chromosome) raises `ParseError`
  at the call, before any row is read.
- Walking a genome in chunks is the caller's job: it is a loop over `fetch`.

The source row index column, `source_row_idx`, is the position of the record
among the records the source returned to the reader, counted from one. For
an unindexed source that is the data line number after the header, also
during a region scan. For an indexed region read it counts the records the
index returned, in order, across the regions of one `fetch` call, so it is
not a file line number. If the mapping names a file column for
`source_row_idx`, that column's whole-number value is used directly.

### Filtering rows

`fetch` also takes two optional row filters on parsed values, by keyword.
A third filter, on the raw text of a column the package does not parse, is
described under "Files that hold many studies" below.

```python
wanted = [("chr1", 100, "A", "G"), ("2", 600, "C", "T")]
for row in reader.fetch(("chr1", 1, 5000), max_pvalue="5e-8", variants=wanted):
    print(row)
```

- **Maximum probability, `max_pvalue`.** A row is kept when its probability
  is at most the threshold.
  - A row exactly at the threshold is kept. A threshold given as a string
    is converted by the same rule as a row's p-value, so the two logarithms
    are equal. A float threshold agrees with the string for every value an
    ordinary float holds at full precision, which is anything at or above
    about 2.2e-308: the float `0.3` and the string `"0.3"` both keep a row
    written `3e-1`.
  - The comparison is made on the negative base-ten logarithm (the
    `mlog10_pvalue` column), not on the probability itself.
  - The threshold may be a float or a decimal string. Give very small
    thresholds as a string: the float `1e-400` is zero in Python and is
    refused, while the string `"1e-400"` is read with decimal arithmetic
    and works. It then keeps a row at `1e-450` and rejects a row at
    `1e-350`, although both probabilities are stored as `0.0`. Decimal
    arithmetic is used for every p-value and threshold too small for a
    float to hold at full precision; larger ones use an ordinary float.
  - The threshold is always a probability, also when the reader was told
    the file holds logged probabilities.
- **Variant membership, `variants`.** A row is kept when it is one of the
  listed variants. Each variant is a tuple of chromosome, position, effect
  allele and other allele.
  - The request goes through the same chromosome, position and allele
    parsers as the rows, including any `chromosome_synonyms` and the
    reader's allele alphabet, so `("01", "150", "a", "g")` names the same
    variant as `("chr1", 150, "A", "G")`.
  - The two alleles are compared as an unordered pair. A row whose effect
    and other allele are swapped relative to the request is a member. The
    row is returned as the file holds it; nothing is flipped.
  - An empty list keeps nothing. `None`, the default, applies no filter.
  - A reader built with `require_other_allele=False` and no other-allele
    column cannot compare variants and raises `ParseError`.
- **Both together.** A row must pass both filters (logical AND).
- **Filters run after parsing.** A filtered read repairs, reports and
  numbers rows exactly as an unfiltered read does. A malformed row is
  raised, or recorded by the dropping reader, even if it would have been
  filtered out, and repair and drop counts match an unfiltered read.
  Filtering never changes a row's `source_row_idx`.
- **Filters never choose what is read.** Membership does not consult a
  genomic index. Building region queries from a variant list, and walking a
  genome in chunks, remain the caller's job; pass regions as well to avoid
  scanning the whole file.
- A malformed threshold or variant raises `ParseError` at the call, before
  any row is read.

### Files that hold many studies

Some files hold many studies side by side. Expression summary statistics
are the usual case: for every gene, all variants within about a megabase,
so one genomic position carries a row for every gene whose window covers
it. A region query alone returns all of those genes mixed together.

Name the study column in the mapping, under any output name that is not a
vocabulary name. The reader copies that cell to the output as raw text,
and `fetch` can keep only the rows you want:

```python
mapping = {
    "chr_name": "chrom",
    "start_pos": "pos",
    # ... the other vocabulary columns ...
    "gene": "gene_id",  # not a vocabulary name: carried through as text
}
reader = DroppingReader(source, mapping, effect_type="beta")
wanted = {"gene": {"ENSG00000139618", "ENSG00000141510"}}
for row in reader.fetch(("chr13", 32_000_000, 33_000_000), column_values=wanted):
    print(row)
```

- **Carried-through columns.** A mapping key that is not a vocabulary name
  names a carried-through column.
  - The cell is copied exactly as the file holds it. Nothing is parsed, and
    missing-value text such as `NA` or a blank stays as it is.
  - Carried-through columns come after every other output column, in the
    order the mapping lists them, in every retention mode.
  - The key may be the file's own header for that column
    (`"gene_id": "gene_id"`) or a new name (`"gene": "gene_id"`). It may
    not be the header of a *different* column in the file, because the
    output would then mislabel the data.
  - **A key that nearly spells a vocabulary name is refused**, so that a
    typing mistake such as `imputaton_info` is caught and not silently
    carried through as text. "Nearly" means within one character (one
    added, dropped or changed, or two neighbours swapped) once case and
    punctuation are ignored. The error names the vocabulary column. One
    realistic name is caught by this: `qvalue` is one letter from
    `pvalue`. Use another output name, for example `qval`.
- **The column-value filter, `column_values`.** A dictionary from
  carried-through column name to the wanted texts. The name is the mapping
  key (`gene` above), not the file's header (`gene_id`).
  - A row is kept when the cell equals one of the wanted texts once
    leading and trailing whitespace is removed from both, so `ENSG1` and
    `ENSG1 ` match each other. Whitespace means anything Python's
    `str.strip()` removes, which includes tabs and non-breaking spaces. Nothing else is changed: case and inner
    spaces must match as written. The cell is still copied to the output
    exactly as the file holds it, padding included. With two columns named,
    both must match.
  - It works on its own, scanning the whole file, and together with
    regions, `max_pvalue` and `variants`. A row must pass everything given.
  - It may not name a vocabulary column. An unknown name, an empty
    collection of wanted texts, a wanted value that is not text, or a
    wanted value that is empty once trimmed, such as `""` or `" "`, raises
    `ParseError` at the call, before any row is read.
- **This filter runs before parsing, unlike the other two.** A row whose
  cell is not wanted is skipped without being parsed.
  - That is what makes it cheap on a file with thousands of studies.
  - The consequence: **a malformed row of an unwanted study is not
    reported.** It is neither raised nor recorded by the dropping reader.
    Only malformed rows of the wanted studies are. A row outside every
    requested region is treated the same way.
  - A line too short to hold the filtered cell cannot be judged, so it is
    refused for its width and reported.
  - Skipped rows still count towards `source_row_idx`.

The opener owns its returned source. Exiting its context manager
closes both the text stream and the index handle;
calling `source.close()` has the same effect. Always close the source, or
use the context manager: rapidgzip aborts the whole Python process if one of
its files is still open when the interpreter shuts down. The reader does not
close caller-owned streams. The caller still positions the source at the intended
header, including skipping any leading comments when needed.

The detection functions `detect_compression`, `detect_index`, and
`detect_opener` in `gwas_sumstats_parser.opening` inspect files without loading
optional dependencies. They report the compression family, neighbouring index
path, and opening family (`plain`, `gzip`, or `pysam`), respectively. Decoder
selection occurs when the source is opened.

## Checking the genome assembly

`check_assembly` reports which of two genome builds an indexed file's
coordinates follow. It is standalone: reading a file never runs it, and it
changes nothing in the file or the reader. It needs a block-compressed file
with a neighbouring tabix index and refuses any other input with an error
saying an index is required; it never scans an unindexed file.

```python
from gwas_sumstats_parser import check_assembly

result = check_assembly("study.tsv.gz", mapping, effect_type="beta")
print(result.decision)      # "GRCh37", "GRCh38", "ambiguous" or "unknown"
print(result.examined)      # probes found at their position in either build
print(result.hits)          # {"GRCh37": 17, "GRCh38": 3}
print(result.scores)        # each build's hits divided by examined
print(result.failed_rows)   # fetched rows skipped because they did not parse
print(result.species, result.builds)
```

The check looks up a table of probe variants, each with a position in build
A and in build B, through the reader's region query. A probe is *examined*
when a row sits at its position in either build, and is a *hit* for a build
when a row sits at that build's position and carries the probe's two alleles
in either order. Chromosome spellings are normalised, so `chr1` in the file
matches `1` in the table.

The mapping and effect type are the reader's, and the reader's required
quantities apply. Further reader settings that describe the file, such as
`logged_pvalue=True`, `chrpos_spec` or `chromosome_synonyms`, are passed as
extra named arguments. Two settings are fixed and refused if supplied: the
other allele is always required, and rows are always read with the widest
allele alphabet. A row that still fails to parse is skipped, counts as no
evidence, and is counted in `failed_rows`.

**Undecided outcomes.** The decision is a build name only when the evidence
supports it:

- `"unknown"` when fewer than `min_examined` probes (default 5) were
  examined, or the better score is below `match_threshold` (default 0.8);
- `"ambiguous"` when both conditions pass but the two scores differ by less
  than `confidence_delta_threshold` (default 0.1);
- otherwise the better-scoring build. A value equal to its threshold passes.

Treat both undecided outcomes as "not confirmed", never as agreement with a
declared build. The evidence fields say why the check was undecided.

**Choosing the probe table.** Without a `probes` argument the check uses the
human GRCh37/GRCh38 table shipped inside the package. For another species,
or different probes, pass a path; nothing else selects a table:

```python
result = check_assembly(
    "mouse_study.tsv.gz", mapping, effect_type="beta", probes="mouse_probes.tsv"
)
```

A probe table is tab-separated. Its first line names the species, its second
names exactly two builds in A then B order, and its header holds at least
these columns:

```text
#species	mouse
#builds	GRCm38	GRCm39
probe_id	chr_a	pos_a	chr_b	pos_b	ref	alt	coverage_score
p1	1	1000	1	1500	A	G	1.0
```

Probe alleles may use the widest allele alphabet, the one file rows are read
with. A supplied table without a species or without two builds is an error; it
never falls back to the human labels. `load_probe_table` reads and validates
a table on its own.

**The bundled human table.** It holds 110 common single-base variants, five
on each of chromosomes 1 to 22. Every one sits at a different position in
GRCh37 and GRCh38, none has a palindromic allele pair, and neither build has
any known variant at a probe's position in the other build, so a file on one
build matches that build only. The variants are on three widely used
genotyping arrays, which makes them likely to be present in a real file: in
a survey of nine real files the worst file lacked 1 of the 110. How the table
was derived is recorded beside it in
`gwas_sumstats_parser/probe_data/assembly_probes.PROVENANCE.md`, and the
survey is `resources/agent-docs/surveys/001-assembly-probe-presence.md`.
Chromosomes X and Y carry no probes.

## Reading a text stream

The caller opens and decodes the stream and skips any leading comments. The
reader consumes its header and accepts standard-library CSV settings, using
tabs by default. It does not open paths.

```python
from io import StringIO
from gwas_sumstats_parser import Reader

stream = StringIO("CHR\tPOS\tEA\tOA\tBETA\tSE\tP\n1\t42\tA\tG\t2\t0.5\t1e-400\n")
mapping = {
    "chr_name": "CHR", "start_pos": "POS",
    "effect_allele": "EA", "other_allele": "OA",
    "effect_size": "BETA", "standard_error": "SE", "pvalue": "P",
}
reader = Reader(stream, mapping, effect_type="beta")
for row in reader:
    probability = row[reader.column_indices["pvalue"]]  # 0.0
    log_probability = row[reader.column_indices["mlog10_pvalue"]]  # 400.0
```

The published `column_names` list gives the output order; `column_indices`
maps each name to its position. Add records to the ordered `COLUMNS` tuple
using `Column(name, type, missing)` and pass it as `columns` to reserve
consumer columns at their defaults. Use `logged_pvalue=True` when the mapped
probability input already contains its negative base-ten logarithm.

Choose returned columns with the `retention` argument:

- `retention="mapped"` is the default: keep mapped columns, including inputs
  consumed by a calculation.
- `retention="lean"` removes mapped inputs consumed to produce another
  quantity. For example, a combined coordinate or confidence bounds disappear
  after they supply required output values. Unused bounds and mapped variant
  identifiers survive.
- `retention="everything"` also appends unmapped source columns, in source
  order, as raw strings under their original headers. Missing-looking raw text
  is preserved. Headers colliding with vocabulary or caller-added column names
  are rejected before data rows.

All modes preserve required values, both probability outputs, source row
indices and caller-added defaults. Vocabulary columns retain template order;
retention changes output columns without disabling calculations.

Use `DroppingReader` with the same configuration to retain usable rows after
a failure. Its `failures` attribute contains source indices, raw rows,
messages, method names and values; `drop_counts` counts each method's failures.
The library installs no logging handlers.
Consumers can catch `ParseError`, read its positional message through the
read-only `error_msg` property, and identify the failing callable with
`error_func.__name__`. The original exception, offending value, source row
index and raw cells remain available. Consumers own bad-row files, failure
thresholds and output annotations. When adopting the reader, account for its
per-row repair policy and do not replace zero probability with a floor: the
preserved logarithm remains the informative value for an underflowed float.

An absent probability column can be derived from effect and standard error,
or from a t-statistic using the large-sample normal approximation. A mapped
probability takes priority. After executing the selected routes, the reader
attempts to repair missing probability, effect and standard error, in that
order. Each missing quantity receives one attempt. Successful row repairs
are counted by quantity in `reader.repair_counts`, without extra output flags.
Standard error can be derived from effect and probability. Effect repair
requires an explicit sign from the effect column under `direction_beta` or
`direction_log_or`; accepted signs are `+`, `-`, `1`, `+1` and `-1`.

All derived probabilities have both an ordinary float and a negative
base-ten logarithm. Probability-driven calculations use the logarithm even
when the ordinary float underflows to zero. Scalar conversion functions live
in `gwas_sumstats_parser.stats`. Probability from signed z or t uses twice the
normal survival function of its absolute value. The t conversion is a
large-sample approximation; it does not apply a degrees-of-freedom correction.
Underflow triggers advancing arbitrary-precision attempts at 1,000, 2,000,
4,000, 6,000, 8,000 and 10,000 decimal places, then a shared parsing error if
the result remains unusable.

Alleles default to the strict A/C/G/T alphabet, with a lone deletion hyphen
also accepted. Set `allele_level="indel_tokens"` to additionally accept lone
D and I tokens, `allele_level="iupac"` for ambiguity letters, or
`allele_level="iupac_indel"` for ambiguity letters and a lone I token.
All levels strip surrounding whitespace, uppercase letters, and reject
embedded hyphens.

Chromosome names are normalized without a species table. Pass
`chromosome_synonyms={"chr01": "one"}` to restrict input to normalized
dictionary keys and return the corresponding names. Missing keys raise the
shared parsing error.

Map a combined-coordinate column using `chrpos`, and describe its layout
with `chrpos_spec`, a format string in the style of a date format. A
percent sign and a letter capture one part: `%c` the chromosome, `%s` the
start position, `%e` the end position, `%f` the effect allele, `%o` the
other allele, `%i` a variant identifier, and `%%` a literal percent sign.
Every other character is literal text that must appear exactly, and the
whole cell must match from start to end. The common four-part cell
`chr1:42:A:G`, where the first of the two alleles is the effect allele, is
read by
`chrpos_spec="%c:%s:%f:%o"`. GTEx's `chr1_13550_G_A_b38`, which lists the
reference allele before the alternative and reports its slope against the
alternative, is read by `chrpos_spec="%c_%s_%o_%f_b38"`: the build tag is
literal text, so a cell without it is refused.

There is no default: when `chrpos` is mapped the format string is required,
so the caller always states which part is the effect allele and the package
never guesses. The chromosome part accepts any run of non-whitespace characters,
so a prefix glued to it, such as `build37-chr1`, is part of the chromosome;
put a prefix you want dropped into the literal text instead. A malformed
format string is refused when the reader is created, naming the fault and
its position. Separate mapped coordinate, allele and variant-identifier
columns take priority when both forms supply the same quantity.

Sample totals prefer a mapped `number_of_samples` column, then the sum of
mapped `number_of_cases` and `number_of_controls`. If neither route is
available, pass `number_of_samples=701` to supply one total for every row.
Counts must be nonnegative integers; zero is valid. Missing or invalid operands
in a selected route raise the shared error. A row whose mapped total, case
count or control count is missing takes the supplied global value instead, and
later routes, such as frequency from an allele count, use it. Such rows are
counted under `number_of_samples` in `reader.repair_counts`. Without a global
value the total stays at 0, the missing value for counts; a zero written in the
file is kept as a real count. The global value must be a whole number of one
or more; anything else, including 0, is refused when the reader is built.

Effect-allele frequency prefers its mapped value, then effect-allele count,
then minor-allele frequency, then minor-allele count. Count routes divide by
twice the sample total, which must be positive and can itself be derived.
Minor-allele routes require a minor allele identity: matching the effect allele
keeps the frequency; a different identity complements it. All frequencies must
be finite and between zero and one. For example, 160 effect-allele copies in
400 diploid people gives 0.2. Derived frequencies are included when their
prerequisites are available; frequency is not a new mandatory input.

The other allele is required by default. Set `require_other_allele=False` to
read a stream without it. Construction resolves required quantities before
reading any data row. Inspect `reader.plan.steps` for the ordered targets,
selected routes and preserved method names; `reader.plan.consumed_inputs`
records mapped columns used to derive other quantities.

For odds, risk and hazard ratios, choose `or`, `rr` or `hr`: positive ratios
become their natural logarithms. Already logged ratios use `log_or`, `log_rr`
or `log_hr`; these and linear effects (`beta`) retain their supplied values.
The standard-error input must already be on the normalized effect scale and
must be positive. A signed t statistic can supply an absent effect by
multiplication with standard error, or an absent standard error by division
of effect by t. Derived effects are already normalized and are not logged
again. Direct values take priority; contradictory effect/t signs fail when
used to derive standard error.

For signed association scores, map the score as `effect_size`. With
`effect_type="z_score_cc"`, the retained conversion estimates a
standardised quantitative per-allele effect. The inherited `cc` label means
correlation coefficient, but the frequency-adjusted output is not claimed
to be a Pearson correlation. This quantitative conversion uses
`SE = 1 / sqrt(2*f*(1-f)*(N+z*z))` and `effect = z*SE`, where `f` is
frequency and `N` is total sample count. The total must be positive and
frequency strictly between zero and one. Both can come from supported
count/frequency derivations. A supplied standard error retains priority.

With `effect_type="z_score_log_or"`, the output effect is a natural-log odds
ratio. A supplied log-odds standard error takes priority and the effect is
its product with the signed score. Otherwise, positive case and control
counts (`C`, `U`) and effect-allele frequency (`f`) give the approximation
`SE = 1 / sqrt(2*f*(1-f)*C*U/(C+U))`, then `log OR = z*SE`.
Frequency must be strictly between zero and one; a total count and frequency
alone do not suffice. Supported frequency derivations may feed this route.
For example, score 2, frequency .5, and 48 cases plus 48 controls yield
SE approximately .288675 and log OR approximately .577350.

This logistic approximation assumes additive genotypes, Hardy-Weinberg
equilibrium and small effects. It does not recover the original fitted
logistic coefficient; strong covariates can change the fitted standard error.
It never substitutes the quantitative-trait calculation. The original signed
score supplies absent probability for both score types. Missing mapped cells
still use the fixed repair pass described above. Scalar score conversions are
available in `gwas_sumstats_parser.case_control`.

For direction-only input, select `direction_beta` for linear effects or
`direction_log_or` for log odds. For quantitative direction input, prefer
supplied standard error, then the eligible sample-count approximation, then
calibrated Burgess. The inherited `set_effect_size_p` and
`set_standard_error_p` methods reconstruct signed z from logarithmic
probability and direction, then use the quantitative-score formula above.
They require a positive total count and frequency strictly between zero and
one. They are eligible for `direction_beta` with an omitted or continuous
trait, never as a log-odds or general binary-trait fallback. A newly available
count can therefore change the chosen method; counts unused by the chosen
Burgess method do not change that calculation.

If an ordinary standard-error source is
unavailable, supply the reference calibration as `average_variation=A`
and map effect-allele frequency. The Burgess approximation is
`SE = A / sqrt(f*(1-f))` and `effect = signed_z * SE`, using the preserved
probability logarithm even when the ordinary probability underflows.
Calibration must be finite and positive, and frequency strictly between
zero and one. No sample count is required or used by this calculation.
An ordinary route, including supplied standard error, retains priority;
an unused Burgess alternative does not require calibration.

The helper `gwas_sumstats_parser.burgess.average_variation` accepts paired
sequences of reference standard errors and allele frequencies. It averages
the finite `SE * sqrt(f*(1-f))` terms. The reader never calls this helper
automatically; it performs no file access and needs no pandas. Calibration
must come from the same study and effect scale as the target variants.
It incorporates study precision and is not the phenotype standard deviation.
The approximation assumes comparable precision across reference and target
variants, including sufficiently similar sample sizes. The caller is
responsible for study provenance and scale compatibility.

Optionally declare the trait with `trait_type`; omission preserves selection
from the input effect interpretation without guessing a trait from columns.
The allowed combinations are:

| Trait | Input interpretations |
| --- | --- |
| `continuous` | `beta`, `direction_beta`, `z_score_cc` |
| `binary` | `beta`, `direction_beta`, `or`, `rr`, `log_or`, `log_rr`, `z_score_log_or`, `direction_log_or` |
| `ordinal` or `count` | `beta`, `direction_beta` |
| `time_to_event` | `hr`, `log_hr` |

An incompatible combination fails before data rows, even if operands for
its formula are available. Declaring a binary trait does not convert a
linear effect to log odds. Risk and hazard ratios keep their respective
logged risk and logged hazard meanings; the trait declaration never changes
the output scale. The parser performs no automatic unit conversion.

Confidence intervals can supply an absent effect or standard error. Map separate
lower and upper bounds as `ci_lower` and `ci_upper`, or a combined interval such
as `(1,3)` as `ci_combined`. Ratio bounds are logged before conversion. The
midpoint gives the effect; the half-width divided by the normal quantile gives
the standard error. An upper bound and an effect also suffice for standard
error. Set `ci_coverage` to the interval's central coverage, default 0.95;
values must be strictly between zero and one. Supplied standard error takes
priority, and unused mapped interval inputs are not marked as consumed.

## Reading into a pandas data frame

The optional module `gwas_sumstats_parser.frame` holds one function,
`read_gwas`, which reads a file and returns two pandas data frames: the
parsed rows, and one line for each row the parser refused.

```python
from gwas_sumstats_parser.frame import read_gwas

mapping = {
    "chr_name": "chromosome",
    "start_pos": "position",
    "effect_allele": "effect_allele",
    "other_allele": "other_allele",
    "effect_size": "beta",
    "standard_error": "se",
    "pvalue": "p",
}
rows, failures = read_gwas(
    "study.tsv.gz",
    mapping,
    effect_type="beta",
    max_pvalue="5e-8",
    verbose=True,
)
```

- It is a thin wrapper. It opens the path with the package's opener, reads
  with the reader that skips bad rows, and builds the two data frames. It
  does not remove duplicates and sets no index column.
- `regions`, `max_pvalue`, `variants` and `column_values` mean exactly what
  they mean under "Reading genomic regions", "Filtering rows" and "Files
  that hold many studies" above. `regions` is a list
  of `(chromosome, start, end)`. Any other named argument is passed to the
  reader unchanged, for example `allele_level="iupac_indel"`.
- `failures` has the columns `source_row_idx`, `message`, `method_name`,
  `value` and `raw_row`. A refused row inside a requested region is
  reported; one outside every requested region is skipped unread. A row
  whose chromosome or position cannot be read cannot be placed, so a scan
  of an unindexed file reports it even during a region query. The p-value
  and variant filters never narrow the failures.
- `verbose=True` shows a counter of the rows kept so far, and logs the
  number of rows kept and refused at the end. The counter has no percentage
  and no time remaining, because the number of rows in a file is unknown
  until it has been read. The default is `False`, which prints nothing.
- The rows are held in memory. On a real file of 33.5 million rows expect
  tens of gigabytes; pass `regions` or `max_pvalue` to read less.

pandas and tqdm are not installed with the package. Install them with the
`frame` extra:

```bash
pip install "gwas-sumstats-parser[frame]"
```

Importing `gwas_sumstats_parser` on its own never imports pandas, so a
package that depends on this one does not inherit it.

## Worked examples

Five Jupyter notebooks under `resources/examples/` show the package on small
made-up files, and are built into the online documentation under "Example
code". Read them in this order:

1. `reading_a_file.ipynb` - map column headers, open a file and read rows.
2. `handling_bad_rows.ipynb` - stop at a bad row, or skip bad rows and keep
   a record of each.
3. `regions_and_filters.ipynb` - read genomic regions from plain and indexed
   files, and keep rows by p-value or by a list of variants.
4. `deriving_missing_quantities.ipynb` - standard errors, p-values,
   frequencies and effects that the file does not carry.
5. `checking_the_genome_build.ipynb` - find out whether a file follows
   GRCh37 or GRCh38.

## What it is not

- It is **not** a normalisation pipeline. Liftover, rsID mapping, variant
  validation, bad-row files and output formatting stay in gwas-norm.
- It is **not** an analysis package. Instrument selection, Mendelian
  randomisation and colocalisation stay in merit.
- It is **not** related to the EBI GWAS Catalog's `gwas-sumstats-tools`,
  `gwas-sumstats-validator` or `gwas-sumstats-harmoniser`, despite the
  shared prefix. Those tools target the GWAS Catalog's own submission
  format. This package targets the heterogeneous files that arrive from
  consortia and individual studies, and standardises them for downstream
  Python analysis.

There is [online](https://cfinan.gitlab.io/gwas-sumstats-parser/index.html)
documentation for gwas-sumstats-parser.

## Installation instructions

Installation can be via pip or conda:

Pip:

```
pip install gwas-sumstats-parser
```

Conda:

```
conda install -c conda-forge -c bioconda -c cfin gwas-sumstats-parser
```

Runtime dependencies are `numpy`, `scipy` and `mpmath`. The optional `io`
extra installs `rapidgzip` and `pysam>=0.17` for faster gzip decoding and
indexed input:

```bash
pip install "gwas-sumstats-parser[io]"
```

The optional `frame` extra installs `pandas` and `tqdm` for the data frame
function `read_gwas` (see "Reading into a pandas data frame"):

```bash
pip install "gwas-sumstats-parser[frame]"
```

Reading caller-owned streams does not require those extras or either
consumer package.

Supported CPython versions are 3.11 through 3.13; development is on 3.13.

### Developer install

Clone the repository and install it in editable mode with the development
extras:

```
git clone git@gitlab.com:cfinan/gwas-sumstats-parser.git
cd gwas-sumstats-parser
python -m pip install -e ".[dev]"
```

Then run the tests from the repository root:

```
pytest tests
```

## Command endpoints

None yet. The package is a library first; command-line endpoints will be
listed here when they are added.

## Project documentation

Design, plans, specifications and known issues live under
`resources/agent-docs/` and are tracked in git. Start with
`resources/agent-docs/PROJECT.md`.
