Skip to content

Preprocessing

Optional. If your genotypes are already filtered to the SNPs you want, skip this page: acquire custom and run are enough.

If you are starting from raw biobank PLINK files, two commands sit between acquire and run. Each reads a cohort directory and writes a new one, so they compose and can be re-run:

acquire  →  preprocess  →  subsample  →  run
            filter SNPs    choose the fit samples

preprocess filters SNPs and never removes a sample; subsample chooses samples and never touches a SNP. Every published flow runs them in that order, so the MAF filter sees the whole cohort before any subset is drawn.

The cohort directory

What acquire writes and what preprocess, subsample and run read:

<dir>/
  config.yaml            names every file below, relative to <dir>
  colormap.json          whole_cohort   (or colormap_fit.json + colormap_project.json for projection)
  data/
    fit_subset.{bed,bim,fam}
    project_subset.{bed,bim,fam}
    labels.csv           whole_cohort   (or fit_labels.csv + project_labels.csv for projection)

A geographic_coords file, when the input names one, is carried through both commands -- filtered to the output's project .fam -- since neither removes a project-side sample; a carried embedding section has its input_mode and landmarking keys (knn, t, n_landmark, random_landmarking) dropped, with a warning per key, whenever the output's preset differs from the input's, so they cannot silently override what the new preset sets.

config.yaml is the entry point, which is why preprocess and subsample take configs rather than PLINK prefixes: the config says where the labels and colormaps are. A label file must cover at least half of the .fam it is paired with — the rule acquire applies — and samples without a row are warned about and drawn grey; extra rows are fine. Both commands check this before doing any work, so a label file that describes the wrong cohort fails in seconds, not after the shell has run. Each command rewrites the labels in its output directory from its input, so they always describe the genotypes beside them.

preprocess — filter SNPs

manifold-genetics preprocess cohort/config.yaml --out filtered/
manifold-genetics preprocess ref/config.yaml biobank/config.yaml --out proj/ --preset harmonise

One config filters a cohort: its own fit and project sets are the two sides. Two configs intersect two cohorts: the first supplies the fit set and its labels, the second the project set, and the output config has preset: projection. Either way the result is a new cohort directory, ready for run, subsample, or another preprocess.

Under the hood this runs preprocess_cross_projection.sh, the shell script that produced the published figures, shipped inside the package. It needs bash, plink2 and plink v1.9 — both binaries are fetched by manifold-genetics setup, or on first use where there is internet.

Presets

preset what runs reproduces
--preset intersect-only indel removal, missingness (--geno), position intersection the examples/generic flows
--preset harmonise WRayner/TOPMed allele checks, GIAB difficult regions, HLA, deduplication, MAF on the reference side only (--skip-project-maf), LD pruning, intersection; --cleanup the All of Us flow
(none) the shell's defaults: every step on, MAF on both sides; add --skip-* flags the UK Biobank flow is --skip-wrayner --skip-project-maf — harmonise minus WRayner, without --cleanup

Every step has a --skip-* flag (--skip-wrayner, --skip-giab, --skip-hla, --skip-ld-prune, --skip-dedup, --skip-maf, --skip-geno, --skip-project-maf), and the thresholds are --maf, --geno, --ld-window, --ld-step, --ld-r2. A skip flag added to a preset is honoured on top of it.

--skip-geno exists for a specific reason. The public HGDP+1KGP panel that acquire hgdp fetches is already filtered and LD-pruned, so preprocessing it should change nothing — but intersect-only still recomputes missingness on the fit and project subsets, and 747 of its 172,152 SNPs (0.43%) newly exceed 5% on those smaller sets. --preset intersect-only --skip-geno on that panel is exactly lossless: same 172,152 SNPs, same samples on both sides. The same is true of the intersected UK Biobank cohort.

Note that the public panel is not a drop-in reference for a harmonised cross-projection: the UK Biobank flag set on it removes 41% of its SNPs (GIAB and HLA exclusion, then LD re-pruning of what is left). The published UK Biobank and All of Us runs each used a different HGDP+1KGP file — the a private gnomAD-derived panel and the workbench's own archive respectively — and a reproduction has to start from the same one.

What needs internet

intersect-only needs nothing beyond plink2. harmonise needs three references — the GIAB difficult-regions bed, the WRayner HRC-1000G-check-bim.pl checker, and the TOPMed reference panel (about 4 GB, almost all of it TOPMed). The checker ships with the package; the shell downloads the other two into the tool cache when they are missing. On a cluster whose compute nodes have no internet, prefetch them on a login node:

manifold-genetics setup --preprocessing

The WRayner checker (v4.3.0, MIT-licensed, by William Rayner) ships with the package, because its upstream URL now returns 404: preprocess installs the bundled copy and needs no network for it. If WRayner runs but writes no corrections — a truncated or wrong TOPMed panel, say — preprocess stops rather than carry on unharmonised. --skip-wrayner skips the step, which is what the UK Biobank flow does.

Resources

Set --threads (default: SLURM_CPUS_PER_TASK, else 4) and --memory in MB (default 100000). Reproducing the UK Biobank projection — 486,748 samples against a 120,849-SNP panel — peaked at 85 GB, so run that on a node with at least 128 GB and leave --memory at its default. --temp-dir moves the scratch space (default OUT/data/temp); --cleanup deletes intermediates as they are consumed, which matters for a biobank-sized .bed.

Troubleshooting

preprocess resumes: a step whose output already exists and is complete is skipped. After an out-of-memory kill, a truncated .bed under OUT/data/temp/ is not mistaken for a finished output — an incomplete intermediate is detected and redone on the next run. Every run also checks OUT/data/temp/preprocess-inputs.json, a sentinel recording the genotype prefixes and every flag that changes what gets filtered (--preset, --maf, --geno, the LD settings, the --skip-* flags, --fit-has-chr-prefix); if it disagrees with this run, preprocess refuses rather than mix intermediates made under different settings — --force does not override this, since it only rewrites the config and labels. Use a new --out, or delete OUT/data/temp, and rerun.

If the two cohorts share fewer than 50,000 SNPs the run aborts rather than producing a projection nobody should trust; --min-common-snps changes the threshold, and usually the right fix is upstream (different builds, or a chr prefix on one side — see --fit-has-chr-prefix).

subsample — choose the fit samples

manifold-genetics subsample proj/config.yaml --out 10k/ \
    --group "self_described_ancestry=^British$:10000" \
    --group "self_described_ancestry=^Irish$:5000"

Reads a cohort directory and writes one whose fit set is a chosen subset of its project set. The project set is linked, not copied, the labels are filtered to match, and the output config uses the subsample preset: fit on the subset, embed the subset, with random landmarking. Given a projection config, the input's fit side (the reference panel) is dropped — the output fits on a subset of the biobank.

Choose the samples by exactly one of:

  • --group COLUMN=PATTERN:COUNT, repeatable. PATTERN is a case-insensitive regex searched for in COLUMN of the label file — a substring match, so British also takes Black or Black British; anchor it (^British$) to mean the whole value. COUNT samples are drawn from the matches with --seed (default 42), or all of them if there are fewer. A sample is taken at most once across groups. --include-rest appends every sample no group matched. For All of Us:

    manifold-genetics subsample proj/config.yaml --out 10k_wbh/ \
    --group "race_ethnicity=White|European:10000" \
    --group "race_ethnicity=Black or African American:10000" \
    --group "race_ethnicity=Hispanic or Latino:10000" \
    --include-rest --seed 42

  • --fit-samples FILE: a FID IID list chosen elsewhere.

    manifold-genetics subsample proj/config.yaml --fit-samples fit_samples.txt --out 10k/

  • --geosketch N: geometric sketching (Hie et al. 2019) — a fit set that covers the cohort's genetic variation evenly instead of following its majority, chosen without any label. It replaces a hand-picked --group balance. Geosketch cannot run on genotypes (hundreds of thousands of columns), so it runs on PCs, in three steps:

    1. PCA is fitted on --sketch-pool samples drawn at random (default 100,000, or everyone if the cohort is smaller) and every sample is projected in — --n-pcs PCs, default 20. Being random, the pool is as imbalanced as the cohort, and so are these PCs.
    2. N samples are sketched in those PCs. They are saved as OUT/sketch_pca.csv; a rerun reuses them rather than fitting again.
    3. run fits a fresh PCA on the N chosen samples, as for any subsample, and embeds them. The imbalanced PCs from step 1 are used for nothing but the choice.
    manifold-genetics subsample proj/config.yaml --geosketch 60000 --out sketch_60k/
manifold-genetics run sketch_60k/config.yaml

    On UK Biobank (487k samples, 120k SNPs) the step-1 fit on 100,000 took about 1.5 h and the projection 15 min; --memory-gb raises the budget above the default 8 GB.

  • --geosketch N --pca CSV: sketch in PCA coordinates you already have, restricted to the project .fam first. --n-pcs limits how many columns of the CSV are used (default: all).

    manifold-genetics subsample proj/config.yaml --geosketch 50000 \
    --pca 10k/outputs/pca/project_pca_20.csv --out sketch/

Worked examples

UK Biobank projected onto HGDP+1KGP

UK Biobank is acquire custom with your own paths — there is nothing for the package to fetch and no label schema to derive; the label CSV needs a sample_id column plus one column per grouping (see Formats). Then intersect the two, skipping WRayner and MAF filtering on the biobank side, as the published UK Biobank preprocessing did:

manifold-genetics acquire custom --fit-plink /path/to/hgdp_tgp --labels hgdp_metadata.csv --out ref/
manifold-genetics acquire custom --fit-plink /path/to/ukb_array --labels ukbb_labels.csv --out ukbb/
manifold-genetics preprocess ref/config.yaml ukbb/config.yaml --skip-wrayner --skip-project-maf --out proj/
manifold-genetics run proj/config.yaml

This reproduces the published cohort exactly from the raw inputs: 3,340 × 120,849 on the fit side and 486,748 × 120,849 on the project side. Peak memory was 85 GB; use a node with at least 128 GB.

All of Us, in the Researcher Workbench

The aou_prepare workflow does this step on a workflow machine and is the tested route. By hand, it is four commands, each resumable, in a terminal or a notebook cell prefixed with !. The workbench has internet, so harmonise can fetch GIAB and TOPMed as it goes; the WRayner checker ships with the package. acquire hgdp writes the workbench's HGDP+1KGP panel with chr-prefixed chromosome names, as the published flow did, hence --fit-has-chr-prefix.

manifold-genetics acquire hgdp --archive gs://fc-secure-47ccf5a8-b9ba-460a-aa03-dea8d260953b/Data/1KGPHGDP.tar.gz --out ref/
manifold-genetics acquire aou --out aou/
manifold-genetics preprocess ref/config.yaml aou/config.yaml --preset harmonise --fit-has-chr-prefix --out proj/
manifold-genetics run proj/config.yaml

acquire aou needs the aou extra and checks the workbench environment first (GOOGLE_PROJECT, WORKSPACE_CDR, gsutil, bq), naming everything missing at once. It has run inside the workbench on the v9 release, as part of aou_prepare. Because that archive carries no relatedness metadata, acquire hgdp --archive fits on all of its samples, where the public panel fits on the 3,400 unrelated.

A biobank subsample

Subsetting comes after preprocessing, so every filter sees the whole cohort:

manifold-genetics preprocess ref/config.yaml ukbb/config.yaml --skip-wrayner --skip-project-maf --out proj/
manifold-genetics subsample proj/config.yaml --out 10k/ \
    --group "self_described_ancestry=^British$:10000" \
    --group "self_described_ancestry=^Irish$:5000"
manifold-genetics run 10k/config.yaml

Or choose the fit set without labels, by geometric sketch (see --geosketch):

manifold-genetics subsample proj/config.yaml --geosketch 60000 --out sketch_60k/
manifold-genetics run sketch_60k/config.yaml