diff --git a/README.md b/README.md index d6f3e7e..a3f12c8 100644 --- a/README.md +++ b/README.md @@ -78,8 +78,14 @@ The Nextflow pipeline automates creation of a **whole-genome** compressed LD arc - **I have a PLINK LD matrix and want to create a `.ldzip` file. What should I do?** Build the **C++ `ldzip` binary** and run the `compress` command. Go to: [C++ Binary](#c-binary) -- **I have PLINK pgen files and want to build whole-genome `.ldzip` outputs in a pipeline. What should I do?** - Use the **Nextflow** workflow. Go to: [Nextflow](#nextflow) +- **I have PLINK pgen/VCF files and want to build whole-genome `.ldzip` outputs in a pipeline. What should I do?** + Use the **Nextflow** workflow. Go to: [Nextflow](#nextflow) or see the [1000 Genomes tutorial](docs/tutorials/g1k-tutorial.md) for an example using VCF files. + +- **I have NPZ files (e.g., UK Biobank LD matrices) and want to build whole-genome `.ldzip` outputs. What should I do?** + Use the **Nextflow** workflow with NPZ input mode. See the [UK Biobank tutorial](docs/tutorials/ukbb-tutorial.md) for an example. + +- **How do I run SuSiE fine-mapping using LDZip-compressed LD matrices?** + See the [UKBB SuSiE analysis tutorial](docs/tutorials/ukbb-susie-analysis.md) for a complete example using UK Biobank GWAS data and LDZip-compressed LD matrices. - **I already have a `.ldzip` file and want to convert it back to my own format. What should I do?** Build the **C++ `ldzip` binary** and run the `decompress` command. Go to: [C++ Binary](#c-binary) diff --git a/docs/.gitignore b/docs/.gitignore new file mode 100644 index 0000000..79ceae7 --- /dev/null +++ b/docs/.gitignore @@ -0,0 +1,7 @@ +tests/.nextflow/ +tests/.nextflow.log* +tests/work/ +tests/*.log +tests/test_output/ +tests/nf_log.txt + diff --git a/docs/README.md b/docs/README.md new file mode 100644 index 0000000..5867e5b --- /dev/null +++ b/docs/README.md @@ -0,0 +1,15 @@ +# LDZip Documentation + +This directory contains tutorials and examples for working with LDZip-compressed LD matrices. + +## Tutorials + +- **[g1k-tutorial.md](tutorials/g1k-tutorial.md)** - Generate whole-genome LDZip matrices from 1000 Genomes VCF files. Demonstrates creation of LDZip matrices from raw genotype data. + +- **[ukbb-tutorial.md](tutorials/ukbb-tutorial.md)** - Generate whole-genome LDZip matrices from UK Biobank pre-computed LD matrices. Demonstrates creation of LDZip matrices from pre-computed LD matrices stored as NPZ files. + +- **[ukbb-susie-analysis.md](tutorials/ukbb-susie-analysis.md)** - Perform SuSiE fine-mapping using LDZip-compressed LD matrices. + +## Tests + +The [tests/](tests/) directory contains executable test scripts for each tutorial that validate the quick-run workflows. diff --git a/docs/tests/test_g1k.sh b/docs/tests/test_g1k.sh new file mode 100755 index 0000000..d8dd61b --- /dev/null +++ b/docs/tests/test_g1k.sh @@ -0,0 +1,52 @@ +#!/bin/bash +# Test script for 1000 Genomes tutorial quick run + +set -e + +# Download VCF for chr20 +mkdir -p test_output/g1k/data + +BASE_URL="ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20201028_3202_phased" +FILE_PREFIX="CCDG_14151_B01_GRM_WGS_2020-08-05" +chr=20 + +wget -O test_output/g1k/data/1000g.chr${chr}.vcf.gz ${BASE_URL}/${FILE_PREFIX}_chr${chr}.filtered.shapeit2-duohmm-phased.vcf.gz +wget -O test_output/g1k/data/1000g.chr${chr}.vcf.gz.tbi ${BASE_URL}/${FILE_PREFIX}_chr${chr}.filtered.shapeit2-duohmm-phased.vcf.gz.tbi + +# Download sample panel +wget -P test_output/g1k/data https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call_samples_v3.20130502.ALL.panel + +# Retrieve list of EUR samples +awk '$3=="EUR" {print $1}' test_output/g1k/data/integrated_call_samples_v3.20130502.ALL.panel > test_output/g1k/data/EUR.txt + +# Create YAML config +cat > test_output/g1k/1000g.yaml < test_output/ukbb/ukbb.yaml <= start && a[2] <= end) print + }' > test_output/ukbb_susie/regions/${REGION}_gwas.tsv + +# Extract variant annotations +tabix -s 2 -b 3 -e 3 -S 1 test_output/ukbb_susie/variants.tsv.bgz +bgzip -dc test_output/ukbb_susie/variants.tsv.bgz | head -1 > test_output/ukbb_susie/regions/${REGION}_vars.tsv +tabix test_output/ukbb_susie/variants.tsv.bgz ${CHR}:${START}-${END} >> test_output/ukbb_susie/regions/${REGION}_vars.tsv + +# Run analysis +# Note: This requires ukbb_ld/european_std to exist. +# Generate it first by running test_ukbb.sh or following the UKBB tutorial. +# See prerequisites in docs/tutorials/ukbb-susie-analysis.md +Rscript -e " +library(data.table) +library(susieR) +library(LDZipMatrix) + +chr <- 10; start <- 7750000; end <- 8350000 +region_name <- sprintf('chr%d_%d_%d', chr, start, end) + +gwas_file <- sprintf('test_output/ukbb_susie/regions/%s_gwas.tsv', region_name) +gwas <- fread(gwas_file) +print(sprintf('Loaded %d GWAS variants', nrow(gwas))) + +vars_file <- sprintf('test_output/ukbb_susie/regions/%s_vars.tsv', region_name) +vars <- fread(vars_file) +print(sprintf('Loaded %d variant annotations', nrow(vars))) + +gwas <- merge(gwas, vars[, .(variant, rsid)], by = 'variant') + +region_query <- sprintf('%d:%d-%d', chr, start, end) +ld_pointer <- LDZipMatrix('ukbb_ld/european_std') +ld_mat <- fetchLD(ld_pointer, region_query, region_query) +print(sprintf('Loaded %d x %d LD matrix', nrow(ld_mat), ncol(ld_mat))) + +common_rsids <- intersect(gwas\$rsid, rownames(ld_mat)) +print(sprintf('%d variants in common', length(common_rsids))) + +gwas <- gwas[match(common_rsids, gwas\$rsid)] +ld_mat <- ld_mat[common_rsids, common_rsids] +print(sprintf('Final: %d variants', nrow(gwas))) + +z <- gwas\$beta / gwas\$se +res <- susie_rss(z, ld_mat, n = median(gwas\$n_complete_samples), L = 10) + +gwas[, pip := res\$pip] +cs <- susie_get_cs(res, coverage = 0.95) + +print(sprintf('Found %d credible sets', length(cs\$cs))) +if (length(cs\$cs) > 0) { + for (i in seq_along(cs\$cs)) { + top_idx <- cs\$cs[[i]][which.max(res\$pip[cs\$cs[[i]]])] + print(sprintf('Credible Set %d: %d variants, top variant %s (PIP=%.3f, p=%.2e)', + i, length(cs\$cs[[i]]), gwas\$rsid[top_idx], res\$pip[top_idx], gwas\$pval[top_idx])) + } +} +" + +echo "UKBB SuSiE analysis test completed successfully!" diff --git a/docs/tutorials/g1k-tutorial.md b/docs/tutorials/g1k-tutorial.md index ebc7939..b53e200 100644 --- a/docs/tutorials/g1k-tutorial.md +++ b/docs/tutorials/g1k-tutorial.md @@ -1,41 +1,92 @@ # Generating Whole-Genome LDZip Matrix from 1000 Genomes -This tutorial walks through generating whole-genome linkage disequilibrium (LD) matrices from 1000 Genomes Project (1000G) data using the Nextflow pipeline provided in this repository. The same workflow can be applied to other large-scale datasets such as UK Biobank or TOPMed. - -**See also:** [UK Biobank tutorial](ukbb-tutorial.md) for generating LD matrices from pre-computed NPZ files. - ## Overview +This tutorial demonstrates generating whole-genome linkage disequilibrium (LD) matrices from 1000 Genomes Project (1000G) phased VCF data. The same workflow can be applied to other large-scale datasets such as UK Biobank or TOPMed. + The workflow consists of two main steps: 1. Download phased VCFs for each chromosome -2. Run a Nextflow pipeline to compute a whole genome LDZip matrix +2. Run a Nextflow pipeline to compute a whole genome LDZip matrix ---- +A quick run option is provided to test a subset of chromosomes and ensure the workflow works correctly before scaling up to the full genome. + +**See also:** [UK Biobank tutorial](ukbb-tutorial.md) for generating LD matrices from pre-computed NPZ files. ## Prerequisites - [`nextflow`](https://www.nextflow.io/docs/latest/install.html) - [`plink2`](https://www.cog-genomics.org/plink/2.0/) -- Sufficient compute (multi-core recommended) -- Adequate storage for final/intermediate files +- HPC cluster (recommended for whole-genome processing) + +--- + +## Setup: Clone and Build + +Before running the pipeline, clone the repository and build the required binaries: + +```bash +# Clone LDZip repository +git clone git@github.com:23andMe/LDZip.git +cd LDZip + +# Build C++ binary +cd cpp +make +cd .. + +# Install R package +cd R +make install +cd .. + +# Return to working directory +cd .. +``` + +After setup, the `ldzip` binary will be at `LDZip/cpp/bin/ldzip` and the R package `LDZipMatrix` will be installed. --- ## Step 1: Download 1000G Phased VCFs +### Quick run + +Download VCF for chromosome 20 for testing, then proceed to Step 2. + +```bash +mkdir -p data + +# Download VCF for chr20 +BASE_URL="ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20201028_3202_phased" +FILE_PREFIX="CCDG_14151_B01_GRM_WGS_2020-08-05" +chr=20 + +wget -O data/1000g.chr${chr}.vcf.gz ${BASE_URL}/${FILE_PREFIX}_chr${chr}.filtered.shapeit2-duohmm-phased.vcf.gz +wget -O data/1000g.chr${chr}.vcf.gz.tbi ${BASE_URL}/${FILE_PREFIX}_chr${chr}.filtered.shapeit2-duohmm-phased.vcf.gz.tbi + +# Download sample panel (2504 samples, phase3) +wget -P data https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call_samples_v3.20130502.ALL.panel +``` + +### Full run + +Download VCFs for all chromosomes. + ```bash mkdir -p data -# download VCFs +# Download VCFs for all chromosomes +BASE_URL="ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20201028_3202_phased" +FILE_PREFIX="CCDG_14151_B01_GRM_WGS_2020-08-05" + for chr in {1..22}; do - wget -O data/1000g.chr${chr}.vcf.gz ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20201028_3202_phased/CCDG_14151_B01_GRM_WGS_2020-08-05_chr${chr}.filtered.shapeit2-duohmm-phased.vcf.gz - wget -O data/1000g.chr${chr}.vcf.gz.tbi ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20201028_3202_phased/CCDG_14151_B01_GRM_WGS_2020-08-05_chr${chr}.filtered.shapeit2-duohmm-phased.vcf.gz.tbi + wget -O data/1000g.chr${chr}.vcf.gz ${BASE_URL}/${FILE_PREFIX}_chr${chr}.filtered.shapeit2-duohmm-phased.vcf.gz + wget -O data/1000g.chr${chr}.vcf.gz.tbi ${BASE_URL}/${FILE_PREFIX}_chr${chr}.filtered.shapeit2-duohmm-phased.vcf.gz.tbi done -# download sample panel (2504 samples, phase3) -wget https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call_samples_v3.20130502.ALL.panel - +# Download sample panel (2504 samples, phase3) +wget -P data https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call_samples_v3.20130502.ALL.panel ``` @@ -43,98 +94,102 @@ wget https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call ## Step 2: Run Nextflow LD Pipeline for European Samples +Retrieve list of EUR samples: -- Clone the LDZip repository: +```bash +awk '$3=="EUR" {print $1}' data/integrated_call_samples_v3.20130502.ALL.panel > data/EUR.txt +``` - ```bash - git clone git@github.com:23andMe/LDZip.git - ``` +### Quick run + +**File: 1000g.yaml** +```yaml +vcf_template: '${launchDir}/data/1000g.chr{CHR}' +keep: '${launchDir}/data/EUR.txt' +outdir: 'output' +ld_command: '--r-unphased ref-based cols=id,ref,alt' +prefix: 'EUR' +chroms: '20' +ld_window_kb: 1000 +ld_window_r2: 0.01 +min_col: 'UNPHASED_R' +ld_threads: 1 +chunk_size_kb: 20000 +overlap_size_kb: 1000 +``` -- Retrieve list of EUR samples +If `plink2` or `ldzip` are NOT available in your `$PATH`, export their paths: - ```bash - awk '$3=="EUR" {print $1}' integrated_call_samples_v3.20130502.ALL.panel > data/EUR.txt - ``` +```bash +export PLINK2=/path/to/plink2 +export LDZIP=$(pwd)/LDZip/cpp/bin/ldzip +``` -- Create YAML file of input parameters +Run as follows: - **File: 1000g.yaml** - ```yaml - vcf_template: '${launchDir}/data/1000g.chr{CHR}' - keep: '${launchDir}/data/EUR.txt' - outdir: 'output' - ld_command: '--r-unphased ref-based cols=id,ref,alt' - prefix: 'EUR' - chroms: '20,21,22' - ld_window_kb: 1 - ld_window_r2: 0.2 - min_col: 'UNPHASED_R' - ld_threads: 1 - chunk_size_kb: 20000 - overlap_size_kb: 1000 - ``` +```bash +nextflow run LDZip/pipelines/wholeGenomeLD/main.nf -params-file 1000g.yaml -resume +``` -- If `plink2` or `ldzip` are NOT available in your `$PATH`, export their paths before running the pipeline (might need to use `set -x` or `setenv` depending on the shell): +Go to Step 3 to test whether it worked correctly. - ```bash - export PLINK2=/path/to/plink2 - export LDZIP=$(pwd)/LDZip/cpp/bin/ldzip - ``` +### Full run -- Run the pipeline: +To run on all chromosomes, update the YAML: + +```yaml +chroms: '1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22' +ld_threads: 8 +``` - ```bash - nextflow run LDZip/pipelines/wholeGenomeLD/main.nf -params-file 1000g.yaml -resume - ``` +For whole-genome processing, you might need an HPC cluster. For example, if using SLURM: - **Tip:** For detailed parameter descriptions and validation rules, run: - ```bash - nextflow run LDZip/pipelines/wholeGenomeLD/main.nf --help - ``` +**File: slurm.config** +```groovy +params.partition = "example_partition" - --- +process.executor = "slurm" +process.cpus = 1 +process.errorStrategy = 'retry' +process.maxRetries = 5 +process.queue = params.partition -- **Note** +process.withName: ldPlink { + cpus = params.ld_threads +} +executor.perCpuMemAllocation = true +``` - The above configuration runs locally for chromosomes 20-22 only (`chroms`), with a small LD window (`ld_window_kb: 1`) and a higher r² threshold (`ld_window_r2: 0.2`) so the pipeline runs quickly (typically ~5 minutes). This serves as a sanity check to ensure the inputs are correct and the workflow completes end-to-end with reasonable outputs. +Run as follows: - To run LDZip on all chromosomes, update these parameters to more reasonable values. Note that `overlap_size_kb` must be greater than or equal to `ld_window_kb`; otherwise, the overlapping regions will not fully cover the specified LD window. +```bash +nextflow run LDZip/pipelines/wholeGenomeLD/main.nf -params-file 1000g.yaml -C slurm.config -resume +``` +For other HPC environments, refer to the Nextflow executor [guidelines](https://www.nextflow.io/docs/latest/executor.html). - ```yaml - chroms: '1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22' - ld_window_kb: 1000 - ld_window_r2: 0.01 - ld_threads: 8 - chunk_size_kb: 20000 - overlap_size_kb: 1000 - ``` +--- -- **Running on HPC** +## Step 3: Verify Output - If you are running with a higher resolution as above, you would want to use a HPC. E.g. if you were using a SLURM cluster you could pass the following config file to your command line +After successful completion, verify the LD matrix by querying two variants: - **File: slurm.config** - ```groovy - params.partition = "example_partition" +```r +library(LDZipMatrix) - process.executor = "slurm" - process.cpus = 1 - process.errorStrategy = 'retry' - process.maxRetries = 5 - process.queue = params.partition +# Load the LD matrix +ld <- LDZipMatrix("output/whole_genome/EUR") - process.withName: ldPlink { - cpus = params.ld_threads - } - executor.perCpuMemAllocation = true - ``` +# Query LD between two variants +fetchLD(ld, "20:64331475:C:T", "20:64333832:C:A") +# [1] 0.976378 - ```bash - nextflow run LDZip/pipelines/wholeGenomeLD/main.nf -params-file 1000g.yaml -C slurm.config -resume - ``` +# Benchmark query time +system.time(fetchLD(ld, "20:64331475:C:T", "20:64333832:C:A")) +# user system elapsed +# 0.003 0.000 0.010 +``` - For other HPC environments, refer to the Nextflow executor [guidelines](https://www.nextflow.io/docs/latest/executor.html). --- ## Parameters diff --git a/docs/tutorials/ukbb-susie-analysis.md b/docs/tutorials/ukbb-susie-analysis.md new file mode 100644 index 0000000..b3ba1c2 --- /dev/null +++ b/docs/tutorials/ukbb-susie-analysis.md @@ -0,0 +1,136 @@ +# UK Biobank SuSiE Analysis with LDZip-compressed LD Matrices + +This guide demonstrates how to perform fine-mapping on UK Biobank GWAS data (from [Neale Lab](http://www.nealelab.is/uk-biobank)) using [SuSiE](https://github.com/stephenslab/susieR) with LDZip-compressed LD matrices (from [Alkes Price Lab](https://labs.icahn.mssm.edu/minervalab/resources/data-ark/ukbb_ld/)). + +**Prerequisites**: +- R packages (`data.table`, `susieR`, `LDZipMatrix`) +- Whole-genome LD file for UKBB at `ukbb_ld/european_std` ([see UKBB LD tutorial](ukbb-tutorial.md)) +- Command-line tools (`bgzip`, `tabix`) + +## Workflow Overview + +This example fine-maps **Asthma (phenotype 20002_1111)** at the ***GATA3* locus (chr10:7.75-8.35 Mb, GRCh38)**. + +1. Download GWAS summary statistics and variant annotations +2. Extract region from whole-genome files +3. Run SuSiE analysis in R + +## Step 1: Download UK Biobank Data for given Phenotype + +```bash +# Create directory +mkdir -p ukbb + +# Download GWAS summary statistics (example: Asthma, phenotype 20002_1111) +wget -P ukbb https://broad-ukb-sumstats-us-east-1.s3.amazonaws.com/round2/additive-tsvs/20002_1111.gwas.imputed_v3.both_sexes.tsv.bgz + +# Download variant annotations with rsIDs +wget -P ukbb https://broad-ukb-sumstats-us-east-1.s3.amazonaws.com/round2/annotations/variants.tsv.bgz +``` + +## Step 2: Extract Region + +```bash +# Define region (GATA3 locus: chr10:7.75-8.35 Mb, GRCh38) +CHR=10 +START=7750000 +END=8350000 +REGION="chr${CHR}_${START}_${END}" + +# Create output directory +mkdir -p ukbb/regions + +# Extract GWAS summary statistics using awk +bgzip -dc ukbb/20002_1111.gwas.imputed_v3.both_sexes.tsv.bgz | \ + awk -F'\t' -v chr="$CHR" -v start="$START" -v end="$END" ' + NR==1 {print; next} + { + split($1, a, ":") + if (a[1] == chr && a[2] >= start && a[2] <= end) print + }' > ukbb/regions/${REGION}_gwas.tsv + +# Extract variant annotations using tabix (will create index if needed) +tabix -s 2 -b 3 -e 3 -S 1 ukbb/variants.tsv.bgz +bgzip -dc ukbb/variants.tsv.bgz | head -1 > ukbb/regions/${REGION}_vars.tsv +tabix ukbb/variants.tsv.bgz ${CHR}:${START}-${END} >> ukbb/regions/${REGION}_vars.tsv +``` + +## Step 3: Run Analysis in R + +The R script loads the pre-extracted files from disk (fast) and performs the analysis: + +```r +#!/usr/bin/env Rscript +# UK Biobank SuSiE Fine-mapping with LDZip +# Note: Run region extraction first to subset data + +library(data.table) +library(susieR) +library(LDZipMatrix) + +# Define region (GATA3 locus: chr10:7.75-8.35 Mb, GRCh38) +chr <- 10; start <- 7750000; end <- 8350000 +region_name <- sprintf("chr%d_%d_%d", chr, start, end) + +# Load pre-extracted GWAS data +gwas_file <- sprintf("ukbb/regions/%s_gwas.tsv", region_name) +gwas <- fread(gwas_file) +print(sprintf("Loaded %d GWAS variants", nrow(gwas))) + +# Load pre-extracted variant annotations +vars_file <- sprintf("ukbb/regions/%s_vars.tsv", region_name) +vars <- fread(vars_file) +print(sprintf("Loaded %d variant annotations", nrow(vars))) + +# Merge to decorate with rsIDs +gwas <- merge(gwas, vars[, .(variant, rsid)], by = "variant") + +# Load LD matrix from LDZip +region_query <- sprintf("%d:%d-%d", chr, start, end) +ld_pointer <- LDZipMatrix("ukbb_ld/european_std") +ld_mat <- fetchLD(ld_pointer, region_query, region_query) +print(sprintf("Loaded %d x %d LD matrix", nrow(ld_mat), ncol(ld_mat))) + +# Intersect and subset data +common_rsids <- intersect(gwas$rsid, rownames(ld_mat)) +print(sprintf("%d variants in common", length(common_rsids))) + +gwas <- gwas[match(common_rsids, gwas$rsid)] +ld_mat <- ld_mat[common_rsids, common_rsids] +print(sprintf("Final: %d variants", nrow(gwas))) + +# Run SuSiE +z <- gwas$beta / gwas$se +res <- susie_rss(z, ld_mat, n = median(gwas$n_complete_samples), L = 10) + +# Extract credible sets +gwas[, pip := res$pip] +cs <- susie_get_cs(res, coverage = 0.95) + +# Summary of results +print(sprintf("Found %d credible sets", length(cs$cs))) +if (length(cs$cs) > 0) { + for (i in seq_along(cs$cs)) { + top_idx <- cs$cs[[i]][which.max(res$pip[cs$cs[[i]]])] + print(sprintf("Credible Set %d: %d variants, top variant %s (PIP=%.3f, p=%.2e)", + i, length(cs$cs[[i]]), gwas$rsid[top_idx], res$pip[top_idx], gwas$pval[top_idx])) + } +} + +# Expected output: +# [1] "Loaded 3704 GWAS variants" +# [1] "Loaded 3704 variant annotations" +# [1] "Loaded 4957 x 4957 LD matrix" +# [1] "3633 variants in common" +# [1] "Final: 3633 variants" +# [1] "Found 3 credible sets" +# [1] "Credible Set 1: 1 variants, top variant rs11567923 (PIP=0.961, p=1.39e-10)" +# [1] "Credible Set 2: 51 variants, top variant rs117158080 (PIP=0.080, p=1.21e-07)" +# [1] "Credible Set 3: 9 variants, top variant rs263424 (PIP=0.406, p=2.49e-05)" +``` + +## References + +- [Neale Lab UK Biobank GWAS](http://www.nealelab.is/uk-biobank) +- [Alkes Price Lab UK Biobank LD](https://labs.icahn.mssm.edu/minervalab/resources/data-ark/ukbb_ld/) +- [SuSiE (Sum of Single Effects)](https://github.com/stephenslab/susieR) diff --git a/docs/tutorials/ukbb-tutorial.md b/docs/tutorials/ukbb-tutorial.md index 700724a..63660a6 100644 --- a/docs/tutorials/ukbb-tutorial.md +++ b/docs/tutorials/ukbb-tutorial.md @@ -1,8 +1,30 @@ # Generating Whole-Genome LDZip Matrix from UK Biobank LD Matrices -This tutorial walks through generating whole-genome linkage disequilibrium (LD) matrices from pre-computed UK Biobank LD matrices (UKBB-LD) using the Nextflow pipeline provided in this repository. +## Overview -**See also:** [1000 Genomes tutorial](g1k-tutorial.md) for generating LD matrices from raw VCF files. +UKBB-LD provides summary linkage disequilibrium (LD) matrices computed from UK Biobank based on N=337K British-ancestry individuals. The LD matrices were computed by the [Alkes Price group at Harvard](https://labs.icahn.mssm.edu/minervalab/resources/data-ark/ukbb_ld/) and are publicly available in AWS S3. The LD information is stored as 2,763 3Mb-long regions spanning the entire genome in NPZ format. + +If you want to use the LDZip-compressed LD matrix generated from this tutorial for SuSiE fine-mapping, see the [UKBB SuSiE analysis tutorial](ukbb-susie-analysis.md). See also the [1000 Genomes tutorial](g1k-tutorial.md) for generating LD matrices from raw VCF files. + +The workflow consists of two main steps: + +1. Download pre-computed NPZ LD matrices from AWS S3 +2. Run a Nextflow pipeline to compress and concatenate into a whole genome LDZip matrix + +A quick run option is provided to test a single chromosome chunk and ensure the workflow works correctly before scaling up to the full genome. + +**Note:** The pipeline automatically handles several data-specific processing steps that are specific to the Alkes Price lab UKBB-LD data format. If you use this on other data, use it with care and verify the output: +- **Boundary variant removal**: The first variant in each 3Mb region (at position ending in 00001) is automatically removed to prevent duplication when concatenating overlapping regions. +- **Multi-allelic variant handling**: RSIDs that appear multiple times (different alleles at same position) are made unique by appending `:REF:ALT` to the variant ID. +- **Exact duplicate removal**: Variants with identical `CHR:POS:REF:ALT` are deduplicated to ensure each variant appears only once in the final matrix. +- **Diagonal correction**: Matrix diagonal values are set to 1.0 (some NPZ files have 0.5 on the diagonal, which is corrected during processing). +- **Matrix asymmetry**: Overlapping regions in the input NPZ files can have slightly different raw LD values for the same variant pair. For example, depending on the raw input file, the pair (`rs2275806`, `rs76602649`) might have a value of `0.4848628044` (in `chr10_6000001_9000001.npz` and `chr10_7000001_10000001.npz`) or `0.4848628640` (in `chr10_8000001_11000001.npz`). When these values fall on opposite sides of a quantization boundary, the final quantized values become different, leading to small asymmetries where `LD[i,j] ≠ LD[j,i]`. This can trigger `XtX is not symmetric` warnings in downstream tools like SuSiE but should not affect the final results. + +## Prerequisites + +- [`nextflow`](https://www.nextflow.io/docs/latest/install.html) +- AWS CLI (for downloading data) +- HPC cluster (recommended for whole-genome processing) --- @@ -20,7 +42,7 @@ cd cpp make cd .. -# Install R package (installs dependencies, builds, and installs) +# Install R package cd R make install cd .. @@ -33,195 +55,106 @@ After setup, the `ldzip` binary will be at `LDZip/cpp/bin/ldzip` and the R packa --- -## Quick Start: 10-Minute Test Run +## Step 1: Download UKBB-LD NPZ Files + +### Quick run -Before downloading all 2,763 files, test with the first 4 chromosome 1 files: +Download single chunk for testing, then proceed to Step 2. ```bash mkdir -p data/ukbb -# Download first 4 files for quick test (chr1 chunks: 0-12Mb) +# Download one chr20 file (chr20: 0-3Mb) S3_BUCKET="s3://broad-alkesgroup-ukbb-ld/UKBB_LD" -chunks=(chr1_1_3000001 chr1_3000001_6000001 chr1_6000001_9000001 chr1_9000001_12000001) - -for chunk in "${chunks[@]}"; do - aws s3 cp --no-sign-request "${S3_BUCKET}/${chunk}.npz" data/ukbb/ & - aws s3 cp --no-sign-request "${S3_BUCKET}/${chunk}.gz" data/ukbb/ & -done - -wait +chunk="chr20_1_3000001" -# Create test config (ukbb_test.yaml) -cat > ukbb_test.yaml << 'EOF' -npz_template: '${launchDir}/data/ukbb/chr{CHR}_{CHUNK}.npz' -outdir: 'output_test' -prefix: 'european_ukbb_test' -chroms: '1' -npz_ld_type: 'UNPHASED_R' -concat_pairwise: true -min: 0.1 -bits: 8 -EOF - -# Run test pipeline -nextflow run LDZip/pipelines/wholeGenomeLD/main.nf -params-file ukbb_test.yaml -resume +aws s3 cp --no-sign-request "${S3_BUCKET}/${chunk}.npz" data/ukbb/ +aws s3 cp --no-sign-request "${S3_BUCKET}/${chunk}.gz" data/ukbb/ ``` -This completes in ~10 minutes and produces a compressed LD matrix for the first 12Mb of chr1 in `output_test/whole_genome/`. - ---- - -## Full Pipeline Overview - -UKBB-LD provides summary linkage disequilibrium (LD) matrices computed from UK Biobank based on N=337K British-ancestry individuals. The LD matrices were computed by the [Alkes Price group at Harvard](https://labs.icahn.mssm.edu/minervalab/resources/data-ark/ukbb_ld/) and are publicly available in AWS S3. The LD information is stored as 2,763 3Mb-long regions spanning the entire genome in NPZ format. - -The workflow consists of two main steps: - -1. Download pre-computed NPZ LD matrices from AWS S3 -2. Run a Nextflow pipeline to compress and concatenate into a whole genome LDZip matrix - ---- - -## Prerequisites +### Full run -- [`nextflow`](https://www.nextflow.io/docs/latest/install.html) -- [`ldzip`](../../README.md#installation) (C++ binary from this repository) -- AWS CLI (for downloading data) -- Sufficient compute (multi-core recommended) -- Adequate storage for final/intermediate files (~15GB for final output) - ---- - -## Step 1: Download UKBB-LD NPZ Files +Download all chunks for all chromosomes. ```bash mkdir -p data/ukbb -# Download all NPZ files from AWS S3 (no AWS credentials required) +# Download all 2,763 NPZ files from AWS S3 (no AWS credentials required) aws s3 sync --no-sign-request \ s3://broad-alkesgroup-ukbb-ld/UKBB_LD/ \ data/ukbb/ ``` -This will download 2,763 NPZ files (one per 3Mb region) with accompanying `.gz` variant metadata files. - --- ## Step 2: Run Nextflow LD Pipeline -- Clone the LDZip repository: - - ```bash - git clone git@github.com:23andMe/LDZip.git - ``` - -- Create YAML file of input parameters - - **File: ukbb.yaml** - ```yaml - npz_template: '${launchDir}/data/ukbb/chr{CHR}_{CHUNK}.npz' - outdir: 'output' - prefix: 'european_ukbb' - chroms: '20,21,22' - npz_ld_type: 'UNPHASED_R' - concat_pairwise: true - min: 0.1 - bits: 8 - ``` - - **Note:** The `npz_template` should point to where you downloaded the files in Step 1. Adjust the path as needed. - -- If `ldzip` is NOT available in your `$PATH`, export its path before running the pipeline: - - ```bash - export LDZIP=$(pwd)/LDZip/cpp/bin/ldzip - ``` - -- Run the pipeline: - - ```bash - nextflow run LDZip/pipelines/wholeGenomeLD/main.nf -params-file ukbb.yaml -resume - ``` - - **Tip:** For detailed parameter descriptions and validation rules, run: - ```bash - nextflow run LDZip/pipelines/wholeGenomeLD/main.nf --help - ``` - ---- +### Quick run -- **Note** +**File: ukbb.yaml** +```yaml +npz_template: '${launchDir}/data/ukbb/chr{CHR}_{CHUNK}.npz' +outdir: 'output' +prefix: 'european_ukbb' +chroms: '20' +npz_ld_type: 'UNPHASED_R' +concat_pairwise: true # IMPORTANT: Required for this dataset since there are regions that overlap with 3 consecutive chunks +min: 0.1 +bits: 8 +``` - The above configuration runs locally for chromosomes 20-22 only (`chroms`) so the pipeline completes quickly (typically ~10 minutes). This serves as a sanity check to ensure the inputs are correct and the workflow completes end-to-end with reasonable outputs. +If `ldzip` is NOT available in your `$PATH`, export its path: - To run LDZip on all chromosomes, update to: +```bash +export LDZIP=$(pwd)/LDZip/cpp/bin/ldzip +``` - ```yaml - chroms: '1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22' - ``` +Run as follows: - **Important:** Set `concat_pairwise: true` for this dataset because the 3Mb regions have overlapping boundaries that cannot be handled by direct concatenation. +```bash +nextflow run LDZip/pipelines/wholeGenomeLD/main.nf -params-file ukbb.yaml -resume +``` -- **Running on HPC** +This runs locally and requires at least 16GB RAM. - For whole-genome processing, you may want to use an HPC cluster. For example, if using SLURM: +Go to Step 3 to test whether it worked correctly. - **File: slurm.config** - ```groovy - params.partition = "example_partition" +### Full run - process.executor = "slurm" - process.cpus = 1 - process.errorStrategy = 'retry' - process.maxRetries = 5 - process.queue = params.partition +To run on all chromosomes, update the YAML: - executor.perCpuMemAllocation = true - ``` +```yaml +chroms: '1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22' +``` - ```bash - nextflow run LDZip/pipelines/wholeGenomeLD/main.nf -params-file ukbb.yaml -C slurm.config -resume - ``` +For whole-genome processing, you might need an HPC cluster. For example, if using SLURM: - For other HPC environments, refer to the Nextflow executor [guidelines](https://www.nextflow.io/docs/latest/executor.html). +**File: slurm.config** +```groovy +params.partition = "example_partition" ---- +process.executor = "slurm" +process.cpus = 1 +process.errorStrategy = 'retry' +process.maxRetries = 5 +process.queue = params.partition -## Expected Output +executor.perCpuMemAllocation = true +``` -After successful completion, the final output will contain: +Run as follows: -``` -output/whole_genome/ -├── european_ukbb.i.bin # Row indices (compressed) -├── european_ukbb.i.bin.index # Index for row indices -├── european_ukbb.meta.json # Metadata -├── european_ukbb.p.bin # Variant positions -├── european_ukbb.sqlite # Variant index database -├── european_ukbb.vars.txt # Variant information -├── european_ukbb.x.UNPHASED_R.bin # LD values (compressed) -└── european_ukbb.x.UNPHASED_R.bin.index # Index for LD values +```bash +nextflow run LDZip/pipelines/wholeGenomeLD/main.nf -params-file ukbb.yaml -C slurm.config -resume ``` -**File sizes for whole genome (all chromosomes):** -``` -5.6G european_ukbb.i.bin -912K european_ukbb.i.bin.index -4.0K european_ukbb.meta.json -148M european_ukbb.p.bin -1.2G european_ukbb.sqlite -510M european_ukbb.vars.txt -5.6G european_ukbb.x.UNPHASED_R.bin -908K european_ukbb.x.UNPHASED_R.bin.index ---- -13G total -``` +For other HPC environments, refer to the Nextflow executor [guidelines](https://www.nextflow.io/docs/latest/executor.html). --- -## Using the LDZip Matrix in R +## Step 3: Verify Output -Once the pipeline completes, you can query the compressed LD matrix using the R package: +After successful completion, verify the LD matrix by querying two variants: ```r library(LDZipMatrix) @@ -230,11 +163,11 @@ library(LDZipMatrix) ld <- LDZipMatrix("output/whole_genome/european_ukbb") # Query LD between two variants -fetchLD(ld, "rs133036", "rs6001980") -# [1] 0.3700787 +fetchLD(ld, "rs995008", "rs4813467") +# [1] 0.5234156 # Benchmark query time -system.time(fetchLD(ld, "rs133036", "rs6001980")) +system.time(fetchLD(ld, "rs995008", "rs4813467")) # user system elapsed # 0.003 0.000 0.010 ```