API reference: c4.finemap¶
The functions behind c4 finemap; see Fine-mapping for how
they fit together.
Fine-mapping: per-session LD reference panels built from the data catalog, the loci to fine-map, and SuSiE-RSS on each locus.
c4 finemap setup converts a cohort's genetic data (a dna: entry in
catalog.yaml, named <cohort>.<entry> e.g. demo.imputed) into one plink2
pgen fileset per chromosome:
<dir>/<panel name>/
chr6.pgen / .pvar / .psam variant IDs are CHROM:POS:REF:ALT
chr6.log plink2 log
samples.keep the samples kept (#IID)
settings.json source, build, filters, and the chromosomes built
Samples are Cohort.samples() for that data set (the id_map's samples,
minus patients marked drop=1), optionally narrowed to populations (--pops, for
a reference panel such as kgp) and to a --keep file. A reference panel's VCFs
can hold more samples than that (kgp b38: 3,202 VCF samples, 2,504 kept).
available_dna(catalog=None)
¶
One row per dna: entry in the catalog, with the name to pass as --dna.
setup_panel(dna, name, out_dir=None, chrom=None, keep=None, pops=None, maf=None, mac=None, geno=None, hwe=None, plink_args=None, dosage=None, overwrite=False, plink2='plink2', threads=None, memory=None, catalog=None)
¶
Build (or add chromosomes to) an LD reference panel at <out_dir>/<name>/.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
dna
|
str
|
Genetic data set as |
required |
name
|
str
|
Panel name; the directory created under |
required |
out_dir
|
Optional[Path]
|
Where to put the panel. Default: $TMPDIR. |
None
|
chrom
|
Optional[str]
|
Chromosome(s) to build, e.g. "6" or "6,7". Default: 1-22. |
None
|
keep
|
Optional[Path]
|
plink-style sample file narrowing the id_map samples. |
None
|
pops
|
Optional[list]
|
Populations or super-populations to keep, from the cohort's pheno pop/super_pop columns (e.g. ["EUR"] for kgp). |
None
|
maf
|
Optional[float]
|
plink2 --maf (unset = not applied; likewise mac, geno, hwe). |
None
|
mac
|
Optional[int]
|
plink2 --mac. |
None
|
geno
|
Optional[float]
|
plink2 --geno. |
None
|
hwe
|
Optional[float]
|
plink2 --hwe. |
None
|
plink_args
|
Optional[list]
|
Extra plink2 arguments, as a list. |
None
|
dosage
|
Optional[bool]
|
Import FORMAT/DS dosages. Default: when the VCF declares DS. |
None
|
overwrite
|
bool
|
Replace a panel built with different settings, or rebuild chromosomes that are already there. Otherwise PanelExistsError. |
False
|
Returns:
| Type | Description |
|---|---|
dict
|
The panel's settings (also written to settings.json), with a per- |
dict
|
chromosome record of variant and sample counts. |
list_panels(root)
¶
Settings of every panel under root (directories with a settings.json).
find_loci(gwas_files, out, p_threshold=5e-08, window=500000, build=None, max_width=3000000, skip_mhc=False, traits=None, trait_col=None, chrom_col=None, pos_col=None, p_col=None)
¶
Read GWAS files, define loci with define_loci(), and write them to out (TSV).
Returns (loci, counts) with counts as from read_hits().
define_loci(hits, window=500000, build=None, max_width=3000000)
¶
Regions to fine-map: ±window around every hit, overlapping windows
merged (pooled across traits, so a region doesn't depend on which trait
found it).
With build, every window touching the extended MHC is merged with the
whole MHC, so chr6 gets one MHC locus rather than overlapping slices of it.
Loci wider than max_width are flagged, not split.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
hits
|
DataFrame
|
CHROM, POS, ID, P, TRAIT, as from read_hits(). |
required |
Returns:
| Type | Description |
|---|---|
DataFrame
|
One row per locus, IDs from coordinates ( |
DataFrame
|
with the lead variant across traits, the number of distinct |
DataFrame
|
significant variants, and each trait's minimum p. |
read_hits(gwas_files, p_threshold=5e-08, traits=None, trait_col=None, chrom_col=None, pos_col=None, p_col=None)
¶
Variants with p < p_threshold in one or more GWAS files.
Each file is one trait, named by traits or from the file name (the
phenotype of a plink2 .glm file); with trait_col, the trait is
read from that column instead, so one file can hold several.
Returns (hits, counts): hits has CHROM, POS, ID, P, TRAIT; counts is {file: (variants read, hits)}.
run_finemap(gwas_file, loci_file, panel, build, out, panel_root=None, locus_ids=None, n=None, L=10, coverage=0.95, min_abs_corr=0.5, estimate_residual_variance=None, drop_ambiguous=False, max_variants=20000, in_sample=False, catalog=None, chrom_col=None, pos_col=None, ref_col=None, alt_col=None, p_col=None, z_col=None, beta_col=None, se_col=None, n_col=None, ea_col=None, rscript='Rscript')
¶
Fine-map every locus in loci_file with SuSiE-RSS, using GWAS z-scores
and LD from a panel made by setup_panel().
Every GWAS variant in a locus is used, whatever its p-value. Variants are matched to the panel on CHROM:POS and alleles (REF/ALT swaps allowed), z is oriented to the panel ALT allele, and LD comes from the panel's dosages (hard calls if the panel was built without DS).
in_sample says the GWAS was run on the panel's own genotypes (refused
for a reference panel such as kgp). SuSiE then estimates the residual
variance (unless estimate_residual_variance is False),
drop_ambiguous is refused (strand can't differ), and loci whose N
differs from the panel's sample count, or with kriging flags, get a note.
estimate_residual_variance defaults to True in-sample, else False.
Writes
{out}.cs.tsv one row per credible set: purity, lead by PIP, lead and min p by GWAS p {out}.cs_members.tsv every CS member with its full GWAS row, PIP and alpha {out}.pip.tsv every fine-mapped variant: z, p, PIP, CS, kriging diagnostics {out}.summary.tsv one row per locus: status, variant counts, drop reasons, N, s, n_cs {out}.dropped.tsv GWAS variants in a locus that weren't fine-mapped, with DROP_REASON {out}.run.json inputs, panel settings, parameters and versions
Returns:
| Type | Description |
|---|---|
dict
|
Run summary dict. |
fine_map_locus(locus, gwas, gwas_chrom, gwas_pos, cols, eff, panel_dir, pvar, n_samples, n=None, dosage=False, rscript='Rscript', L=10, coverage=0.95, min_abs_corr=0.5, estimate_residual_variance=False, drop_ambiguous=False, max_variants=20000, in_sample=False)
¶
Fine-map one locus. Returns {"summary": dict, "pip": DataFrame, "cs": [rows], "members": [DataFrames], "dropped": DataFrame}. Problems that only affect this locus are reported in summary["status"] rather than raised; things worth a look go in summary["notes"].
With in_sample (the GWAS was run on the panel's own genotypes), a
sample size that differs from the panel's and any kriging flags are noted,
since with in-sample LD neither should happen.
prepare_locus(gwas, cols, eff, panel, drop_ambiguous=False)
¶
Match one locus's GWAS rows to panel variants and orient z to the panel ALT allele.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
gwas
|
DataFrame
|
GWAS rows in the locus. |
required |
panel
|
DataFrame
|
panel variants in the locus (PANEL_ID, CHROM, POS, REF, ALT). |
required |
Returns:
| Type | Description |
|---|---|
tuple
|
(matched, dropped). matched has the GWAS columns plus PANEL_ID, MATCH, |
tuple
|
Z (oriented to ALT), P, N and AMBIGUOUS_STRAND; dropped has the GWAS |
tuple
|
columns plus DROP_REASON. |
effect_columns(columns, cols, z_col=None, beta_col=None, se_col=None, n_col=None, ea_col=None)
¶
Which GWAS columns give z, N and the effect allele.
z is read directly (Z_STAT, Z, T_STAT), else BETA/SE, else log(OR)/SE. The effect allele is A1 / EFFECT_ALLELE / ALLELE1, else the ALT column.
Returns {"z", "beta", "or", "se", "n", "ea": column or None, "z_from": text}.
Locus plots for c4 finemap visualize: GWAS p-values, PIPs and genes.
One figure per locus, three panels sharing the x-axis (Mb):
-log10(p) coloured by r² with a credible-set lead (LocusZoom bins);
credible-set members get one marker shape per set, leads are
purple diamonds, kriging-flagged variants have a black ring
PIP same encoding
genes from the build's GTF in the catalog (resources with
gene_annotation: gtf)
Everything comes from the outputs of c4 finemap run ({prefix}.pip.tsv,
.cs.tsv, .summary.tsv, .run.json), so the LD panel isn't needed.
visualize(prefix, out_dir, locus_ids=None, include_all=False, full_locus=False, fmt='png', color_by='cs', genes=True, gtf=None, all_biotypes=False, catalog=None)
¶
One plot per locus of a c4 finemap run, written to
<out_dir>/<chrom>_<start>-<end>.<fmt>.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
prefix
|
Path
|
The run's output prefix (as given to |
required |
locus_ids
|
Optional[list]
|
Only these loci. Default: every locus with a credible set. |
None
|
include_all
|
bool
|
Also plot loci that ran but have no credible set. |
False
|
full_locus
|
bool
|
Show the whole locus window, not just the span of its fine-mapped variants (plus a margin). |
False
|
color_by
|
str
|
"cs" (r² with the nearest credible-set lead) or "top" (r² with the top-PIP variant). |
'cs'
|
genes
|
bool
|
Draw the gene track, from |
True
|
Returns:
| Type | Description |
|---|---|
list
|
The paths written. |
plot_locus(pip, cs, locus_id, genes=None, summary=None, build=None, color_by='cs')
¶
Plot one locus from its rows of pip.tsv and cs.tsv. color_by is "cs"
(r² with the nearest credible-set lead) or "top" (r² with the top-PIP
variant). genes is None for no gene track. Returns the matplotlib figure.
read_gtf_genes(path, chroms=None, biotypes=('protein_coding',))
¶
Gene records from a GTF (plain or gzipped): chrom, start, end, strand,
name, biotype. Only chroms (if given) and biotypes (None = all).