Fine-mapping¶

c4 finemap fine-maps GWAS loci with SuSiE-RSS. For each locus it reports credible sets (small groups of variants that very likely contain a causal variant) and each variant's posterior inclusion probability (PIP).

It takes three steps:

Step	Command	What it does
1	`c4 finemap setup`	Builds an LD reference panel from a cohort's genotypes, or from 1000 Genomes
2	`c4 finemap loci`	Turns GWAS hits into regions to fine-map
3	`c4 finemap run`	Matches each region's GWAS variants to the panel and runs SuSiE
4	`c4 finemap visualize`	Plots each locus: p-values, PIPs and genes

You need plink2 for step 1, and R with susieR for step 3. On C4:

module load plink2 r

Quick start¶

# 1. A panel from the cohort the GWAS was run on (chr6 only)
c4 finemap setup --dna demo.imputed --panel-name demo --chr 6 --maf 0.001

# 2. ±500 kb around every p < 5e-8 variant
c4 finemap loci results.PHENO.glm.linear --build b38 -o loci.tsv

# 3. SuSiE on each locus
c4 finemap run results.PHENO.glm.linear --loci loci.tsv --panel demo --build b38 \
    --in-sample -o fm/pheno

Then look at fm/pheno.cs.tsv for the credible sets, and at fm/pheno.summary.tsv to check that every locus ran cleanly. To plot every locus that has a credible set:

c4 finemap visualize fm/pheno -o plots/

1. Build an LD panel¶

SuSiE needs the correlation (LD) between the variants in each locus. It gets it from a panel: a set of genotypes that setup converts to plink2 files.

c4 finemap avail

avail lists every genetic data set in the data catalog, named <cohort>.<entry>, plus any panels you have already built:

Genetic data in the catalog (use the first column as --dna):

          dna      kind format build chroms status source
 demo.imputed    cohort    vcf   b38   1-22  22/22 ✓
  demo.exome     cohort      -     -      -  not set
      kgp.b38 reference    vcf   b38   1-22  22/22 ✓ https://ftp.1000genomes.ebi.ac.uk/...
      kgp.b37 reference    vcf   b37   1-22  22/22 ✓ https://ftp.1000genomes.ebi.ac.uk/...

Then build a panel from one of them:

# The GWAS cohort's own genotypes (best, if the GWAS was run on them)
c4 finemap setup --dna demo.imputed --panel-name demo --chr 6

# 1000 Genomes, European samples only
c4 finemap setup --dna kgp.b38 --panel-name kgp_eur --pops EUR --chr 6

Use the GWAS cohort itself when you can

SuSiE-RSS is sensitive to differences between the panel's LD and the LD in the GWAS samples. If you have the genotypes the GWAS was run on, a panel built from them (in-sample LD) is much more reliable than 1000 Genomes.

Where it goes. Panels are written to $TMPDIR/<panel-name>/, one set of plink2 files per chromosome (chr6.pgen/.pvar/.psam), plus samples.keep and settings.json. On Slurm, $TMPDIR is deleted when the job ends, so build the panel in the same job that uses it (see Running it as a Slurm job). Use --dir to put it somewhere else.

Samples. The panel uses the cohort's samples for that data set from its ID map, minus patients marked drop=1. You can narrow that:

--keep file.txt keeps only those samples (plink format: IID, or FID IID).
--pops EUR (or --pops GBR,CEU) keeps populations or super-populations, for reference panels such as 1000 Genomes.

For 1000 Genomes, kgp.b38 (30x high coverage) and kgp.b37 (phase 3) both give the same 2,504 unrelated samples. The b38 VCFs' 698 extra samples, mostly relatives, are marked drop=1.

Variants. Variant IDs are CHROM:POS:REF:ALT. Imputed dosages (DS) are kept when the VCF has them, and LD is then computed from dosages rather than hard calls. Optional filters: --maf, --mac, --geno, --hwe, and --plink-args "..." for anything else plink2 takes. Records with more than two alleles are dropped (imputed VCFs are already split), as are repeated records of the same variant.

Running again. Running setup again with the same --panel-name adds chromosomes. If the panel was built with different settings, or a chromosome you asked for is already there, you're asked before anything is overwritten. --force skips the question, and batch jobs need it.

2. Define the loci¶

c4 finemap loci results.PHENO.glm.linear --build b38 -o loci.tsv

Every variant with p < --p (default 5e-8) gets a window of --window (default 500 kb) on each side. Overlapping windows are merged into one locus.

--p only decides where the loci are. Every variant inside a locus is fine-mapped, whatever its p-value. SuSiE needs the non-significant variants to model the LD.
Windows are merged by distance, not by LD. Two nearby signals end up in one locus, and SuSiE separates them. This is different from clumping, which would split them.
Locus IDs are coordinates, such as 6:31500000-32500000, so they stay the same between runs and when you add traits.
The MHC is one locus. With --build, any window that touches the extended MHC (b37 chr6:28,477,797-33,448,354; b38 chr6:28,510,120-33,480,577) is merged into a single MHC locus, marked is_mhc. SuSiE results there are hard to interpret because of the MHC's extreme LD; --skip-mhc leaves it out.
Wide loci are flagged, not split. Loci wider than --max-width (3 Mb) are marked exceeds_max_width.

Several traits. Give several GWAS files and the loci are pooled across them. A locus doesn't depend on which trait found it, which is what colocalization needs later.

c4 finemap loci trait1.tsv trait2.tsv --trait T1 --trait T2 --build b38 -o loci.tsv

Trait names come from the plink2 phenotype in the file name, or from --trait. For a single file holding several traits, use --trait-col.

loci.tsv has one row per locus:

Column	Meaning
`locus_id`, `chrom`, `start`, `end`, `width`	The region
`lead_id`, `lead_p`, `lead_trait`	Most significant variant in the region, across traits
`n_sig`	Number of distinct variants with p < `--p`
`traits`, `trait_min_p`	Each trait with a hit, and its minimum p
`is_mhc`, `exceeds_max_width`	Flags

You can edit this file, or write your own with chrom, start and end columns, to fine-map other regions.

3. Run SuSiE¶

c4 finemap run results.PHENO.glm.linear --loci loci.tsv --panel demo --build b38 \
    --in-sample -o fm/pheno

For each locus, run:

Takes every GWAS variant in the region.
Matches it to a panel variant on position and alleles. An exact REF/ALT match is preferred; a swapped match is accepted when it's the only one. Preferring exact matches matters for indels: an insertion and a deletion at the same position, such as C/CA and CA/C, have the same two alleles.
Orients z to the panel's ALT allele, so z and LD use the same allele.
Computes LD from the panel's dosages.
Runs susie_rss in R, along with two checks that the z-scores fit the LD (see Checking the results).

If one locus fails or is skipped, the others still run; summary.tsv says what happened to each one.

Use unfiltered summary statistics. SuSiE needs every variant in a locus, including the null ones, to model the LD. Sumstats cut down to significant or suggestive hits make PIPs and credible sets overconfident. run warns when the GWAS's largest p-value is below 0.1, and notes any locus whose variants all have p < 1e-3.

--build must match the panel. The run stops if the GWAS build and the panel's build (from its settings.json) differ.

Where z, N and the effect allele come from. z is read from Z_STAT, Z or T_STAT, else computed as BETA/SE, else log(OR)/SE. The effect allele is A1, EFFECT_ALLELE or ALLELE1, else ALT. The sample size is the median OBS_CT or N of the locus's variants, unless you give --n. The --*-col options override any of these. plink2 .glm files work as they are.

--in-sample. Use it when the GWAS was run on the panel's own genotypes. It's refused for reference panels such as 1000 Genomes. With it:

SuSiE estimates the residual variance, as susieR recommends for in-sample LD. --fixed-residual-variance turns that off.
--drop-ambiguous is refused, because the strand can't differ between the GWAS and the panel.
A locus whose GWAS N differs from the panel's sample count by more than 5%, or that has kriging flags, gets a note in summary.tsv. With true in-sample LD neither should happen.

Other options. --locus ID runs only that locus, and can be repeated. --L is the maximum number of signals per locus (default 10). --coverage (default 0.95) and --min-abs-corr (default 0.5, the purity cut-off) control credible sets. --max-variants (default 20,000) skips very large loci.

Output files¶

All files start with the -o prefix.

File	One row per
`.cs.tsv`	Credible set: size, coverage, purity, log10 Bayes factor, lead variant by PIP, lead variant and minimum p by GWAS p
`.cs_members.tsv`	Credible set member, with its full GWAS row, PIP and alpha
`.pip.tsv`	Fine-mapped variant: z, p, PIP, credible set, r² with the top variant and with the nearest credible-set lead, kriging diagnostics, strand-ambiguity flag
`.summary.tsv`	Locus: status, variant counts, drop reasons, N used, LD source, `s`, number of credible sets, notes
`.dropped.tsv`	GWAS variant in a locus that wasn't fine-mapped, with `DROP_REASON`
`.run.json`	The run: inputs, panel settings, parameters, susieR version

Credible set IDs are <locus_id>:L<k>, for example 6:31500000-32500000:L1.

Reading cs.tsv.

Purity (min_abs_corr, mean_abs_corr) is the LD between the set's members. Sets below --min-abs-corr are already removed. A pure set means its variants are hard to tell apart, not that the signal is weak.
Two leads. lead_pip_variant is SuSiE's best guess at the causal variant. lead_p_variant and min_p are the most significant GWAS variant in the set. They often differ, and both are useful.
Full GWAS rows. cs_members.tsv has the full GWAS row for every member, so you don't have to look members up again.

Checking the results¶

Check summary.tsv before trusting a locus's credible sets.

Column	What to look for
`status`	`ok`, or why the locus was skipped or failed
`n_gwas`, `n_matched`, `n_variants`	GWAS variants in the region, matched to the panel, and fine-mapped. A large gap means much of the locus is missing from the panel
`dropped`	Counts by reason (see below)
`s`	From `estimate_s_rss`: how inconsistent the z-scores and LD are. Near 0 is good; above about 0.1–0.2, read the locus with caution
`n_kriging_flagged`	Variants whose z doesn't fit the LD (see below)
`notes`	Warnings: every variant has p < 1e-3 (the GWAS looks filtered by p-value); with `--in-sample`, an N mismatch or kriging flags

Kriging flags. kriging_rss predicts each variant's z-score from its neighbours' z-scores and the panel LD. A variant is flagged (kriging_flag in pip.tsv) when its observed z fits much better with its alleles flipped (logLR > 2) and it has a real signal (|z| > 2). Common causes:

One or two flagged variants with a strong z: an allele flip in the sumstats, a strand-ambiguous SNP (check ambiguous_strand), or a badly imputed variant.
Many flagged variants across the locus: the panel doesn't match the GWAS samples, for example a different ancestry. Use in-sample LD if you can.

Why variants are dropped (DROP_REASON in dropped.tsv):

Reason	Meaning
`not_in_panel`	No panel variant at that position
`allele_mismatch`	A panel variant at that position, but with different alleles
`ambiguous`	Matches more than one panel variant, and the alleles can't decide which. With REF/ALT given this is rare; with only `A1` it happens at multi-allelic sites
`duplicate_in_gwas`	Another GWAS row matched the same panel variant; the more significant one is kept
`missing_p`, `missing_pos`, `missing_z`	Missing or unreadable values
`effect_allele_mismatch`	The effect allele is neither panel allele
`monomorphic_in_panel`	No variation among the panel's samples, so it has no LD
`strand_ambiguous`	A/T or C/G SNP, dropped because of `--drop-ambiguous`

Plotting a locus¶

c4 finemap visualize fm/pheno -o plots/                       # every locus with a credible set
c4 finemap visualize fm/pheno -o plots/ --locus 6:31500000-32500000 --format pdf

Each plot has three panels sharing the x-axis:

−log₁₀(p), a local Manhattan plot, with the genome-wide significance line.
PIP, SuSiE's posterior inclusion probability.
Genes, protein-coding genes from the catalog's Ensembl GTF for the run's build (--all-biotypes for every gene, --gtf for another file, --no-genes to leave the track out).

In both data panels:

Colour is r² with the nearest credible-set lead, in the usual LocusZoom bins. Each signal gets its own red cloud around its own lead. --color-by top uses r² with the single top-PIP variant instead.
Shape: credible-set members get one marker per set (triangle for L1, square for L2, …). Each set's lead (top PIP) is a purple diamond, labelled with its ID.
Black ring: the variant is kriging-flagged (see Checking the results).

The title gives the number of variants, N, the number of credible sets, s, any kriging flags, and whether the LD was in-sample. The MHC is shaded on chr6.

The x-axis spans the locus's fine-mapped variants plus a small margin; --full-locus shows the whole locus window instead.

Plots are drawn from the run's output files, so the LD panel doesn't need to exist any more: you can make them after the Slurm job has finished, or on another machine. --all also plots loci without a credible set.

Running it as a Slurm job¶

Because the panel lives in $TMPDIR, build it and use it in the same job:

#!/bin/bash
#SBATCH --job-name=finemap
#SBATCH --cpus-per-task=4
#SBATCH --mem=32G
#SBATCH --time=12:00:00

module load plink2 r

c4 finemap setup --dna demo.imputed --panel-name demo --force
c4 finemap loci results.PHENO.glm.linear --build b38 -o loci.tsv
c4 finemap run results.PHENO.glm.linear --loci loci.tsv --panel demo --build b38 \
    --in-sample -o fm/pheno
c4 finemap visualize fm/pheno -o plots/

To save time, build only the chromosomes your loci are on, e.g. --chr 6,19.

From Python¶

The same steps are available as functions; see the API reference.

from c4.finemap import setup_panel, find_loci, run_finemap

setup_panel("demo.imputed", "demo", chrom="6")
loci, _ = find_loci(["results.PHENO.glm.linear"], "loci.tsv", build="b38")
result = run_finemap("results.PHENO.glm.linear", "loci.tsv", "demo", "b38", "fm/pheno",
                     in_sample=True)
result["cs"]          # the credible sets, as a DataFrame