Fix dichotomization method, add tiered precision, allele harmonization in first pass

- Add --use-beta-sign flag implementing the Southam et al. 2017 dichotomization
  method (sign of beta instead of p-value transform) for correlation estimation.
  The old method severely underestimates inter-study correlation under the null.

- Add tiered arbitrary precision: long double -> 200 -> 500 -> 1000 digits,
  with safe fallbacks for p-values too extreme for any tier.

- Add allele harmonization in first pass so that --use-beta-sign works correctly
  when effect alleles differ between input files.

- Fix missing return in initialise() (undefined behavior, caused crashes).

- Update README with new build instructions and --use-beta-sign option.

- Add test_data/ with GIANT Height analysis notebook, null simulation script,
  and generate.sh to reproduce all data files from scratch.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: arthur gilly

CHANGELOG.md

@@ -0,0 +1,68 @@
+# Changelog
+
+## v1.1.0
+
+### Bug fix: dichotomization method for correlation estimation
+
+The original implementation used `b_transform(p) = Phi^{-1}(1-p) <= 0` to dichotomize
+p-values before computing tetrachoric correlation between studies. This ignores the
+direction of effect and severely underestimates the inter-study correlation,
+leading to undercorrection for sample overlap.
+
+The paper (Southam et al. 2017) describes a different procedure:
+`z_k = Phi^{-1}(p_k/2) * sgn(beta_k)`, then dichotomize on the sign of the
+resulting z-score. This method correctly recovers the inter-study correlation
+under the null hypothesis.
+
+**Null simulation results** (100,000 variants under the null):
+
+| True rho | Old method rho-hat | New method rho-hat | Old lambda_GC | New lambda_GC |
+|----------|--------------------|--------------------|---------------|---------------|
+| 0.0      | -0.005             | 0.006              | 1.017         | 1.011         |
+| 0.1      | 0.008              | 0.097              | 1.090         | 1.042         |
+| 0.3      | 0.054              | 0.284              | 1.258         | 1.129         |
+| 0.5      | 0.159              | 0.475              | 1.375         | 1.181         |
+| 0.7      | 0.360              | 0.670              | 1.401         | 1.185         |
+
+The new method is enabled with `--use-beta-sign` and will become the default in a
+future version. The old method remains the default for backwards compatibility.
+
+Credit: Brady Ryan (University of Michigan) identified the discrepancy between the
+paper and the implementation.
+
+### Bug fix: allele harmonization in first pass
+
+The first pass (correlation matrix estimation) did not harmonize alleles between
+input files. When effect alleles were swapped between studies, the sign of beta
+was incorrect, producing wrong tetrachoric correlation estimates. The first pass
+now detects allele flips and negates beta accordingly, mirroring the allele
+harmonization already present in the second pass (meta-analysis).
+
+### Bug fix: missing return in `initialise()`
+
+The `initialise()` function (`int inline initialise(int, char**)`) had no `return`
+statement. This is undefined behavior in C++ and caused crashes (double-free) on
+some compilers with `-O3`.
+
+### Improvement: tiered arbitrary precision arithmetic
+
+The previous implementation had a single fallback from `long double` to
+`cpp_dec_float<200>`, which overflowed for p-values below ~1e-200.
+
+All four transform functions (`b_transform`, `z_transform` x2, `z_transform_fess`)
+now use a 4-tier precision cascade:
+
+```
+long double -> cpp_dec_float<200> -> cpp_dec_float<500> -> cpp_dec_float<1000> -> safe fallback
+```
+
+The final fallback caps z-scores at +/-38 (equivalent to p ~ 1e-315), which is
+reached only for the most extreme p-values in GWAS (e.g. GPC3 height locus at
+p = 2.8e-1135).
+
+### Build changes
+
+- Makefile now uses system Boost (`/usr/include`, `lib/`) instead of hardcoded
+  Sanger cluster paths. Original paths preserved as comments.
+- Removed `-static` flag (dynamic linking). For static builds, pass
+  `CFLAGS="-O3 -std=c++11 -static"` to make.

CLAUDE.md

@@ -0,0 +1,50 @@
+# CLAUDE.md
+
+This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
+
+## Project Overview
+
+METACARPA (META-analysis in C++ Accounting for Relatedness, using arbitrary Precision Arithmetic) is a C++11 bioinformatics tool for meta-analysing genetic association studies with overlapping or related samples when the degree of overlap is unknown. It combines the p-value correction method of Province and Borecki (2013) with the effect-size method of Lin and Sullivan (2009).
+
+## Build
+
+The project is a single C++ source file compiled via Make. Boost 1.60.0+ is required (program_options, serialization, multiprecision, uBLAS, math).
+
+```bash
+# Build (produces static Linux x64 binary at src/metacarpa)
+make -C src/
+
+# The Makefile IDIR and LDIR variables must point to your Boost headers and libs
+```
+
+Compiler flags: `-O3 -std=c++11 -static`. Links against: `boost_program_options`, `boost_serialization`, `pthread`.
+
+## Architecture
+
+The entire implementation lives in a single file: `src/metacarpa.cpp` (~1850 lines).
+
+**Two-phase algorithm:**
+1. **Matrix calculation** — Reads all input files, identifies overlapping variant positions across studies, and computes a variance-covariance correlation matrix using tetrachoric correlation of binarized p-values. The matrix is serialized to disk (`.matrix.txt`) for reuse across runs.
+2. **Single-point meta-analysis** — Re-reads input files and performs per-variant meta-analysis using the precomputed correlation matrix, producing corrected effect sizes, Z-scores, and p-values.
+
+**Key components in `metacarpa.cpp`:**
+- `studies_correlation` class — Stores and manages correlation matrices indexed by study "mask" (bitmask of which studies contribute to a variant). Handles Boost serialization (XML/text).
+- `count_data` struct (nested in `studies_correlation`) — 2x2 contingency tables for tetrachoric correlation.
+- `output_record` struct — Holds per-variant meta-analysis results.
+- `meta_analyse()` — Core computation: inverts correlation submatrix, computes corrected Z-scores and effect sizes.
+- `z_transform()` / `z_transform_fess()` — P-value to Z-score conversion with automatic fallback to `cpp_dec_float<200>` arbitrary precision for extremely small p-values (< ~1e-29).
+- `InvertMatrix()` — LU-decomposition matrix inversion via Boost uBLAS.
+- `tetrachoric()` — Estimates correlation between studies from binary-transformed p-value vectors.
+
+**Data flow:** Parse CLI args → open sorted input files → first pass builds correlation matrix → serialize matrix → second pass performs per-variant meta-analysis → write tab-delimited output.
+
+## Key Constraints
+
+- Input files must be sorted by chromosome then position.
+- Alleles are case-sensitive (expects capitals); allele flipping is supported but strand flipping is not.
+- Sample sizes can be per-variant (via `--size-col`) or constant (appended to filename: `-I file,N`).
+- The correlation matrix can be precomputed with `-x` (stop after matrix) and reused with `-m` for different traits or chunked data.
+
+## No Test Suite
+
+There is no automated test framework. Validation is done externally via the [metacarpa-simulation](https://bitbucket.org/agilly/metacarpa-simulation) project.

README.md

@@ -15,9 +15,24 @@ METACARPA is compiled as a static Linux x64 executable. You can download it in t
 
 ### Compilation
 
-The src directory contains a single source file and Makefile. The [Makefile](src/Makefile) is written for a compiler that supports C++ 11. Please note that you will also need a copy of the [Boost](http://www.boost.org) mathematical libraries on your machine (the binaries were compiled using Boost 1.60.0).
+The src directory contains a single source file and Makefile. A C++11 compiler and [Boost](http://www.boost.org) 1.60.0+ are required (headers + `program_options` and `serialization` libraries).
 
-After compiling the Boost libraries, you will be provided with the `-I` and `-L` arguments to change in the Makefile.
+On Debian/Ubuntu:
+```bash
+sudo apt-get install libboost-program-options-dev libboost-serialization-dev
+```
+
+Then build:
+```bash
+cd src && make
+```
+
+The Makefile expects Boost headers in `/usr/include` and libraries in `lib/` (relative to `src/`). Adjust `IDIR` and `LDIR` in the Makefile if your Boost installation is elsewhere.
+
+For a static binary (no runtime dependencies):
+```bash
+cd src && c++ -O3 -std=c++11 -static -I/usr/include -L/usr/lib/x86_64-linux-gnu -o metacarpa metacarpa.cpp -lboost_program_options -lboost_serialization -lpthread
+```
 
 ## Program Options
 
@@ -54,6 +69,9 @@ Options description :
   -x [ --stop ]         Stop METACARPA after generating the matrix.
   -d [ --debug ]        Toggles an extremely verbose output, for debugging
                         purposes only.
+  --use-beta-sign       Use sign of beta for dichotomization in correlation
+                        matrix calculation (Southam et al. 2017), instead of
+                        the default p-value transform. Recommended.
 ```
 
 ## Data preparation

src/Makefile

@@ -1,9 +1,11 @@
 #IDIR=/nfs/users/nfs_a/ag15/boost_1_55_0
-IDIR=/nfs/team144/software/boost_1.6/boost_1_60_0
+#IDIR=/nfs/team144/software/boost_1.6/boost_1_60_0
+IDIR=/usr/include
 CC=c++
 #LDIR=/nfs/users/nfs_a/ag15/boost_1_55_0/stage/lib
-LDIR=/nfs/team144/software/boost_1.6/boost_1_60_0/stage/lib
-CFLAGS=-O3 -std=c++11 -static -I$(IDIR) -L$(LDIR)
+#LDIR=/nfs/team144/software/boost_1.6/boost_1_60_0/stage/lib
+LDIR=lib
+CFLAGS=-O3 -std=c++11 -I$(IDIR) -L$(LDIR)
 LIBS=-lboost_program_options -lboost_serialization -lpthread
 
 metacapa: metacarpa.cpp