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 and 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.

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:

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.

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.

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:

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.

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:

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:

#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.

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.

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:

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 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:

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"):

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.