Add reproducible GWAS fixtures and relational manifests

README.md

@@ -17,8 +17,27 @@ graph TD
     B --> D[CONCAT_CHUNKED_VCFS]
     B --> F[EXTRACT_SAMPLE_IDS]
     F --> G[GENERATE_PHENO_COV]
+    H[results/chunked_vcfs] --> I[FIXTURE_SELECT_SAMPLES]
+    H --> J[FIXTURE_THIN_CHROMOSOME]
+    I --> J
+    J --> K[FIXTURE_MERGE_CHROMOSOMES]
+    K --> L[FIXTURE_PCA]
+    K --> M[FIXTURE_CAUSAL_WEIGHTS]
+    M --> N[FIXTURE_SCORE_LIABILITY]
+    K --> N
+    L --> O[FIXTURE_PHENO_COVAR]
+    N --> O
+    I --> O
+    J --> P[FIXTURE_EXPORT_FORMATS]
+    K --> P
+    P --> Q[FIXTURE_COMPRESS_VCF]
 ```
 
+The two stages are independent. The upper stage downloads the 1000 Genomes release
+and chunks it; its output is already committed, so it is skipped with
+`--skip_source_generation`. The lower stage derives the GWAS pipeline fixtures from
+that committed output and never touches the network.
+
 ## Git clone the gwas pipeline test data
 
 If you want to get a local copy of the test data, you can either git clone the whole test data material, including all test data for all nf-core pipelnies, or if you want to save storage space you can clone the example data for one specific pipeline.
@@ -60,11 +79,30 @@ results/
 ├── pheno_cov/
 │   ├── example.pheno
 │   └── example.covar
+├── fixtures/
+│   ├── README.md
+│   ├── genotypes/
+│   │   ├── example_chr1.{pgen,pvar,psam,bed,bim,fam,vcf.gz,vcf.gz.tbi}
+│   │   ├── example_chr2.{...}
+│   │   ├── ...                     (one bundle per autosome, chr1..chr22)
+│   │   ├── example_chr22.{...}
+│   │   └── example_all.{...}
+│   └── pheno_cov/
+│       ├── example.pheno
+│       ├── example.covar
+│       ├── example.qcovar
+│       ├── example.catcovar
+│       ├── causal_weights.tsv
+│       ├── fixture_samples.tsv
+│       ├── fixture_sex.tsv
+│       └── fixture_summary.txt
 
 ```
 
 Each chromosome-specific VCF file (chr\*.vcf.gz) is accompanied by its corresponding tabix index (.vcf.gz.tbi), enabling efficient querying. A combined VCF and index are also included for downstream association tests or visualization.
 
+`results/fixtures/` holds the pipeline-level GWAS fixtures: 500 samples across all 22 autosomes, 4942 LD-thinned variants (`--indep-pairwise 100 10 0.2`), the same samples and variants in PLINK 2, PLINK 1 and VCF encodings, a headered phenotype file with a simulated signal-bearing quantitative trait and a simulated binary trait, and headered covariate files carrying sex, age and principal components. They are derived from `results/chunked_vcfs/` with no download step. See [`results/fixtures/README.md`](results/fixtures/README.md) for the dimensions, how the traits were simulated, and how to regenerate them.
+
 ## Support
 
 For further information or help, don't hesitate to get in touch on our [Slack organisation](https://nf-co.re/join/slack) (a tool for instant messaging).

main.nf

@@ -7,23 +7,87 @@ include { CONCAT_CHUNKED_VCFS } from './modules/concat_chunked_vcfs.nf'
 include { EXTRACT_SAMPLE_IDS } from './modules/extract_sample_ids.nf'
 include { GENERATE_PHENO_COV } from './modules/generate_pheno_cov.nf'
 include { INDEX_CHUNKED_VCFS } from './modules/index_chunked_vcfs.nf'
+include { FIXTURE_SELECT_SAMPLES } from './modules/fixture_select_samples.nf'
+include { FIXTURE_THIN_CHROMOSOME } from './modules/fixture_thin_chromosome.nf'
+include { FIXTURE_MERGE_CHROMOSOMES } from './modules/fixture_merge_chromosomes.nf'
+include { FIXTURE_PCA } from './modules/fixture_pca.nf'
+include { FIXTURE_CAUSAL_WEIGHTS } from './modules/fixture_causal_weights.nf'
+include { FIXTURE_SCORE_LIABILITY } from './modules/fixture_score_liability.nf'
+include { FIXTURE_PHENO_COVAR } from './modules/fixture_pheno_covar.nf'
+include { FIXTURE_EXPORT_FORMATS } from './modules/fixture_export_formats.nf'
+include { FIXTURE_COMPRESS_VCF } from './modules/fixture_compress_vcf.nf'
 workflow {
-    // Run the download process
-    GENERATE_EXAMPLE_GENOTYPES_VCFS()
-
-    def vcfs_with_chr = GENERATE_EXAMPLE_GENOTYPES_VCFS.out.vcfs
-    	.flatten()
-    	.map { file ->
-        	def chr = file.name.toString().split("\\.")[1] // safer than `tokenize`
-        	tuple(chr, file)
-    	}
-
-    // Feed the tuples into the chunking process
-    CHUNK_VCFS(vcfs_with_chr)
-    chr1_chunked = CHUNK_VCFS.out.chunked_vcfs.filter { chr, file -> chr == 'chr1'}.map { chr, file -> file }
-    EXTRACT_SAMPLE_IDS(chr1_chunked)
-    GENERATE_PHENO_COV(EXTRACT_SAMPLE_IDS.out.sample_ids)
-    INDEX_CHUNKED_VCFS(CHUNK_VCFS.out.chunked_vcfs)
-    all_chunked_vcfs = CHUNK_VCFS.out.chunked_vcfs.map { chr, file -> file }.collect()
-    CONCAT_CHUNKED_VCFS(all_chunked_vcfs)
+    // Source-data stage: download the 1000 Genomes release and chunk it. Skipped
+    // with --skip_source_generation, because its output is already committed under
+    // results/chunked_vcfs and re-running it re-downloads tens of gigabytes.
+    if (!params.skip_source_generation) {
+        // Run the download process
+        GENERATE_EXAMPLE_GENOTYPES_VCFS()
+
+        def vcfs_with_chr = GENERATE_EXAMPLE_GENOTYPES_VCFS.out.vcfs
+            .flatten()
+            .map { file ->
+                def chr = file.name.toString().split("\\.")[1] // safer than `tokenize`
+                tuple(chr, file)
+            }
+
+        // Feed the tuples into the chunking process
+        CHUNK_VCFS(vcfs_with_chr)
+        chr1_chunked = CHUNK_VCFS.out.chunked_vcfs.filter { chr, file -> chr == 'chr1'}.map { chr, file -> file }
+        EXTRACT_SAMPLE_IDS(chr1_chunked)
+        GENERATE_PHENO_COV(EXTRACT_SAMPLE_IDS.out.sample_ids)
+        INDEX_CHUNKED_VCFS(CHUNK_VCFS.out.chunked_vcfs)
+        all_chunked_vcfs = CHUNK_VCFS.out.chunked_vcfs.map { chr, file -> file }.collect()
+        CONCAT_CHUNKED_VCFS(all_chunked_vcfs)
+    }
+
+    // GWAS fixture stage: derive the pipeline-level test fixtures from the
+    // already-committed chromosome-chunked VCFs. Nothing here reads the network.
+    def fixture_chrs = params.fixture_chromosomes.tokenize(',').collect { chr -> chr.trim() }
+    def fixture_sources = channel
+        .fromList(fixture_chrs)
+        .map { chr -> tuple(chr, file("${params.chunked_vcfs_dir}/${chr}_chunked.vcf.gz", checkIfExists: true)) }
+
+    // The sample subset is taken from the first fixture chromosome; every chunked
+    // VCF carries the same 2504 samples in the same order.
+    FIXTURE_SELECT_SAMPLES(
+        file("${params.chunked_vcfs_dir}/${fixture_chrs.first()}_chunked.vcf.gz", checkIfExists: true)
+    )
+
+    FIXTURE_THIN_CHROMOSOME(
+        fixture_sources,
+        FIXTURE_SELECT_SAMPLES.out.samples,
+        FIXTURE_SELECT_SAMPLES.out.sex
+    )
+
+    FIXTURE_MERGE_CHROMOSOMES(
+        fixture_chrs,
+        FIXTURE_THIN_CHROMOSOME.out.pgen
+            .map { prefix, pgen, pvar, psam -> [pgen, pvar, psam] }
+            .collect()
+    )
+
+    FIXTURE_PCA(FIXTURE_MERGE_CHROMOSOMES.out.pgen)
+
+    FIXTURE_CAUSAL_WEIGHTS(
+        FIXTURE_MERGE_CHROMOSOMES.out.pgen.map { prefix, pgen, pvar, psam -> pvar }
+    )
+
+    FIXTURE_SCORE_LIABILITY(
+        FIXTURE_MERGE_CHROMOSOMES.out.pgen,
+        FIXTURE_CAUSAL_WEIGHTS.out.weights
+    )
+
+    FIXTURE_PHENO_COVAR(
+        FIXTURE_SCORE_LIABILITY.out.sscore,
+        FIXTURE_PCA.out.eigenvec,
+        FIXTURE_SELECT_SAMPLES.out.sex
+    )
+
+    // Per-chromosome and combined bundles are both emitted: leave-one-chromosome-out
+    // methods need the per-chromosome split, single-cohort methods want the merge.
+    def fixture_bundles = FIXTURE_THIN_CHROMOSOME.out.pgen.mix(FIXTURE_MERGE_CHROMOSOMES.out.pgen)
+
+    FIXTURE_EXPORT_FORMATS(fixture_bundles)
+    FIXTURE_COMPRESS_VCF(FIXTURE_EXPORT_FORMATS.out.vcf)
 }

modules/extract_sample_ids.nf

@@ -1,6 +1,6 @@
 process EXTRACT_SAMPLE_IDS {
 	container "community.wave.seqera.io/library/r-base:4.4.3--1e564c44feffeaa0"
-        publishDir params.outdir_pheno_cov, mode: 'symlink'
+        publishDir params.outdir_pheno_cov, mode: 'copy'
 
 	input:
 	path vcf_file

modules/fixture_causal_weights.nf

@@ -0,0 +1,55 @@
+// Draw the causal variants and their effect sizes for the simulated traits.
+//
+// Two disjoint causal sets are drawn from the merged fixture variants, one for the
+// quantitative trait and one for the binary trait's genetic liability, so the two
+// traits are not collinear. Every draw goes through a fixed seed with the RNG kind
+// pinned explicitly, so the weight file is byte-identical on every run.
+process FIXTURE_CAUSAL_WEIGHTS {
+	container "community.wave.seqera.io/library/r-base:4.4.3--1e564c44feffeaa0"
+	publishDir params.outdir_fixture_pheno_cov, mode: 'copy'
+
+	input:
+	path pvar
+
+	output:
+	path "causal_weights.tsv", emit: weights
+
+	script:
+	"""
+	#!/usr/bin/env Rscript
+	set.seed(${params.fixture_seed},
+			 kind = "Mersenne-Twister",
+			 normal.kind = "Inversion",
+			 sample.kind = "Rejection")
+
+	variants <- read.table("${pvar}", comment.char = "#", header = FALSE,
+						   stringsAsFactors = FALSE)
+	colnames(variants) <- c("CHROM", "POS", "ID", "REF", "ALT")
+	n_variants <- nrow(variants)
+
+	n_causal <- ${params.fixture_n_causal}
+	if (2 * n_causal > n_variants) {
+		stop("not enough variants for two disjoint causal sets")
+	}
+
+	drawn <- sample.int(n_variants, 2 * n_causal)
+	causal_qt <- sort(drawn[seq_len(n_causal)])
+	causal_bt <- sort(drawn[n_causal + seq_len(n_causal)])
+
+	weights <- data.frame(ID = variants[["ID"]],
+						  A1 = variants[["ALT"]],
+						  BETA_QT = 0,
+						  BETA_BT = 0,
+						  stringsAsFactors = FALSE)
+	weights[["BETA_QT"]][causal_qt] <- round(rnorm(n_causal), 6)
+	weights[["BETA_BT"]][causal_bt] <- round(rnorm(n_causal), 6)
+
+	write.table(weights, file = "causal_weights.tsv", sep = "\\t", quote = FALSE,
+				row.names = FALSE, col.names = TRUE)
+	"""
+
+	stub:
+	"""
+	printf 'ID\\tA1\\tBETA_QT\\tBETA_BT\\n' > causal_weights.tsv
+	"""
+}

modules/fixture_compress_vcf.nf

@@ -0,0 +1,27 @@
+// Block-compress and index the exported fixture VCFs.
+//
+// The ##fileDate header line plink2 writes carries the run date, which would make
+// the committed fixture differ on every regeneration; it is dropped so the VCF is
+// byte-reproducible.
+process FIXTURE_COMPRESS_VCF {
+	container "community.wave.seqera.io/library/bcftools_tabix_pip_tools:48085064a9189d8c"
+	publishDir params.outdir_fixture_genotypes, mode: 'copy'
+
+	input:
+	tuple val(prefix), path(vcf_file)
+
+	output:
+	tuple val(prefix), path("${prefix}.vcf.gz"), path("${prefix}.vcf.gz.tbi"), emit: vcf
+
+	script:
+	"""
+	grep -v '^##fileDate=' ${vcf_file} | bgzip -c > ${prefix}.vcf.gz
+	tabix -p vcf ${prefix}.vcf.gz
+	"""
+
+	stub:
+	"""
+	printf '##fileformat=VCFv4.3\\n#CHROM\\tPOS\\tID\\tREF\\tALT\\tQUAL\\tFILTER\\tINFO\\tFORMAT\\n' | bgzip -c > ${prefix}.vcf.gz
+	touch ${prefix}.vcf.gz.tbi
+	"""
+}

modules/fixture_export_formats.nf

@@ -0,0 +1,36 @@
+// Re-encode a PLINK 2 fixture bundle as PLINK 1 binary and as VCF, so all three
+// encodings describe exactly the same samples and variants.
+process FIXTURE_EXPORT_FORMATS {
+	container "community.wave.seqera.io/library/plink2:2.0.0a.6.9--e6710830a4b7f0c6"
+	publishDir params.outdir_fixture_genotypes, mode: 'copy', pattern: "*.{bed,bim,fam}"
+
+	input:
+	tuple val(prefix), path(pgen), path(pvar), path(psam)
+
+	output:
+	tuple val(prefix), path("${prefix}.bed"), path("${prefix}.bim"), path("${prefix}.fam"), emit: bed
+	tuple val(prefix), path("${prefix}.vcf"), emit: vcf
+
+	script:
+	"""
+	plink2 \\
+		--pfile ${prefix} \\
+		--threads 1 \\
+		--seed 1 \\
+		--make-bed \\
+		--out ${prefix}
+
+	plink2 \\
+		--pfile ${prefix} \\
+		--threads 1 \\
+		--seed 1 \\
+		--export vcf id-paste=iid \\
+		--out ${prefix}
+	"""
+
+	stub:
+	"""
+	touch ${prefix}.bed ${prefix}.bim ${prefix}.fam
+	printf '##fileformat=VCFv4.3\\n#CHROM\\tPOS\\tID\\tREF\\tALT\\tQUAL\\tFILTER\\tINFO\\tFORMAT\\n' > ${prefix}.vcf
+	"""
+}

modules/fixture_merge_chromosomes.nf

@@ -0,0 +1,36 @@
+// Merge the per-chromosome PLINK 2 fixtures into one multi-chromosome bundle.
+//
+// The merge list is built from the ordered chromosome list rather than from the
+// staged file order, so the merged variant order is reproducible.
+process FIXTURE_MERGE_CHROMOSOMES {
+	container "community.wave.seqera.io/library/plink2:2.0.0a.6.9--e6710830a4b7f0c6"
+	publishDir params.outdir_fixture_genotypes, mode: 'copy', pattern: "*.{pgen,pvar,psam}"
+
+	input:
+	val chrs
+	path psets
+
+	output:
+	tuple val("${params.fixture_prefix}_all"), path("${params.fixture_prefix}_all.pgen"), path("${params.fixture_prefix}_all.pvar"), path("${params.fixture_prefix}_all.psam"), emit: pgen
+
+	script:
+	def prefix = "${params.fixture_prefix}_all"
+	def merge_list = chrs.collect { chr -> "${params.fixture_prefix}_${chr}" }.join('\\n')
+	"""
+	printf '${merge_list}\\n' > merge_list.txt
+
+	plink2 \\
+		--pmerge-list merge_list.txt pfile \\
+		--threads 1 \\
+		--seed 1 \\
+		--out ${prefix}
+	"""
+
+	stub:
+	def prefix = "${params.fixture_prefix}_all"
+	"""
+	touch ${prefix}.pgen
+	printf '#CHROM\\tPOS\\tID\\tREF\\tALT\\n' > ${prefix}.pvar
+	printf '#FID\\tIID\\tSEX\\n' > ${prefix}.psam
+	"""
+}

modules/fixture_pca.nf

@@ -0,0 +1,29 @@
+// Principal components of the merged fixture genotypes, used as covariates.
+// Exact (non-approximate) PCA, so the result is a deterministic function of the
+// genotype matrix.
+process FIXTURE_PCA {
+	container "community.wave.seqera.io/library/plink2:2.0.0a.6.9--e6710830a4b7f0c6"
+
+	input:
+	tuple val(prefix), path(pgen), path(pvar), path(psam)
+
+	output:
+	path "${prefix}.eigenvec", emit: eigenvec
+	path "${prefix}.eigenval", emit: eigenval
+
+	script:
+	"""
+	plink2 \\
+		--pfile ${prefix} \\
+		--pca ${params.fixture_n_pcs_computed} \\
+		--threads 1 \\
+		--seed 1 \\
+		--out ${prefix}
+	"""
+
+	stub:
+	"""
+	printf '#FID\\tIID\\tPC1\\n' > ${prefix}.eigenvec
+	printf '0.0\\n' > ${prefix}.eigenval
+	"""
+}