Public update docs

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)

docs/.gitignore

@@ -0,0 +1,7 @@
+tests/.nextflow/
+tests/.nextflow.log*
+tests/work/
+tests/*.log
+tests/test_output/
+tests/nf_log.txt
+

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.

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 <<EOF
+vcf_template: '\${launchDir}/test_output/g1k/data/1000g.chr{CHR}'
+keep: '\${launchDir}/test_output/g1k/data/EUR.txt'
+outdir: 'test_output/g1k/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
+EOF
+
+# Export ldzip path if needed
+export LDZIP=$(pwd)/../../cpp/bin/ldzip
+
+# Run pipeline
+nextflow run ../../pipelines/wholeGenomeLD/main.nf -params-file test_output/g1k/1000g.yaml -resume
+
+# Verify output
+Rscript -e "
+library(LDZipMatrix)
+ld <- LDZipMatrix('test_output/g1k/output/whole_genome/EUR')
+result <- fetchLD(ld, '20:64331475:C:T', '20:64333832:C:A')
+cat('LD between 20:64331475:C:T and 20:64333832:C:A:', result, '\n')
+"
+
+echo "1000 Genomes test completed successfully!"

docs/tests/test_ukbb.sh

@@ -0,0 +1,41 @@
+#!/bin/bash
+# Test script for UKBB tutorial quick run
+
+set -e
+
+# Download single chunk for testing
+mkdir -p test_output/ukbb/data
+
+S3_BUCKET="s3://broad-alkesgroup-ukbb-ld/UKBB_LD"
+chunk="chr20_1_3000001"
+
+aws s3 cp --no-sign-request "${S3_BUCKET}/${chunk}.npz" test_output/ukbb/data/
+aws s3 cp --no-sign-request "${S3_BUCKET}/${chunk}.gz" test_output/ukbb/data/
+
+# Create YAML config
+cat > test_output/ukbb/ukbb.yaml <<EOF
+npz_template: '\${launchDir}/test_output/ukbb/data/chr{CHR}_{CHUNK}.npz'
+outdir: 'test_output/ukbb/output'
+prefix: 'european_ukbb'
+chroms: '20'
+npz_ld_type: 'UNPHASED_R'
+concat_pairwise: true
+min: 0.1
+bits: 8
+EOF
+
+# Export ldzip path if needed
+export LDZIP=$(pwd)/../../cpp/bin/ldzip
+
+# Run pipeline
+nextflow run ../../pipelines/wholeGenomeLD/main.nf -params-file test_output/ukbb/ukbb.yaml -resume
+
+# Verify output
+Rscript -e "
+library(LDZipMatrix)
+ld <- LDZipMatrix('test_output/ukbb/output/whole_genome/european_ukbb')
+result <- fetchLD(ld, 'rs995008', 'rs4813467')
+cat('LD between rs995008 and rs4813467:', result, '\n')
+"
+
+echo "UKBB test completed successfully!"

docs/tests/test_ukbb_susie.sh

@@ -0,0 +1,87 @@
+#!/bin/bash
+# Test script for UKBB SuSiE analysis tutorial
+
+set -e
+
+# Create directory
+mkdir -p test_output/ukbb_susie
+
+# Download GWAS summary statistics
+wget -P test_output/ukbb_susie 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
+wget -P test_output/ukbb_susie https://broad-ukb-sumstats-us-east-1.s3.amazonaws.com/round2/annotations/variants.tsv.bgz
+
+# Define region (GATA3 locus)
+CHR=10
+START=7750000
+END=8350000
+REGION="chr${CHR}_${START}_${END}"
+
+mkdir -p test_output/ukbb_susie/regions
+
+# Extract GWAS summary statistics
+bgzip -dc test_output/ukbb_susie/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
+      }' > 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!"