Add GDA EGT cluster INFO BCF source

illumina_microarray/copy_egt_info_bcf_to_main.sh

@@ -0,0 +1,14 @@
+#!/usr/bin/env bash
+
+# Promote EGT-derived sample-less BCF + CSI from cpg-common-test to cpg-common-main.
+# Run via analysis-runner at --access-level full --dataset common.
+
+set -ex
+
+gcloud storage cp \
+    gs://cpg-common-test/references/illumina_microarray/GDA-8v1-0_D1_ClusterFile_info.bcf \
+    gs://cpg-common-main/references/illumina_microarray/GDA-8v1-0_D1_ClusterFile_info.bcf
+
+gcloud storage cp \
+    gs://cpg-common-test/references/illumina_microarray/GDA-8v1-0_D1_ClusterFile_info.bcf.csi \
+    gs://cpg-common-main/references/illumina_microarray/GDA-8v1-0_D1_ClusterFile_info.bcf.csi

illumina_microarray/extract_egt_info_bcf.sh

@@ -0,0 +1,105 @@
+#!/usr/bin/env bash
+#
+# Re-transcribe the per-variant INFO payload of an Illumina EGT cluster file
+# into a sample-less BCF using `bcftools +gtc2vcf` with no `--gtcs`. The output
+# is the "expected" cluster geometry (GenTrain_Score, Cluster_Sep, N_AA/AB/BB,
+# meanTHETA_*, devTHETA_*, meanR_*, devR_*) keyed on (CHROM, POS, REF, ALT) —
+# the reference half of any observed-vs-EGT comparison.
+#
+# Run locally (or inside the bcftools:1.23-1 image) against the canonical
+# BPM / EGT / fasta references staged in cpg-common-main. Staging the output
+# to cpg-common-test and promoting to cpg-common-main are handled out-of-band
+# (see illumina_microarray/copy_egt_info_bcf_to_main.sh for the test->main
+# promote step).
+#
+# Requirements: bcftools (with the gtc2vcf plugin) on PATH.
+#
+# Example:
+#   ./extract_egt_info_bcf.sh \
+#     --bpm   /path/to/GDA-8v1-0_D2.bpm \
+#     --egt   /path/to/GDA-8v1-0_D1_ClusterFile.egt \
+#     --fasta /path/to/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna \
+#     --out-dir ./out
+
+set -euo pipefail
+
+OUT_NAME='GDA-8v1-0_D1_ClusterFile_info.bcf'
+
+BPM=''
+EGT=''
+FASTA=''
+OUT_DIR='./out'
+
+usage() {
+    cat <<EOF
+Usage:
+  $0 --bpm BPM --egt EGT --fasta FASTA [--out-dir DIR]
+
+  --bpm      Path to the Illumina BPM manifest (.bpm)
+  --egt      Path to the Illumina EGT cluster file (.egt)
+  --fasta    Path to the reference fasta (.fna); .fai must sit alongside
+  --out-dir  Local output directory (default: ./out)
+EOF
+}
+
+while [[ $# -gt 0 ]]; do
+    case "$1" in
+        --bpm)     BPM="$2"; shift 2 ;;
+        --egt)     EGT="$2"; shift 2 ;;
+        --fasta)   FASTA="$2"; shift 2 ;;
+        --out-dir) OUT_DIR="$2"; shift 2 ;;
+        -h|--help) usage; exit 0 ;;
+        *) echo "Unknown argument: $1" >&2; usage >&2; exit 2 ;;
+    esac
+done
+
+mkdir -p "$OUT_DIR"
+OUT_BCF="${OUT_DIR}/${OUT_NAME}"
+OUT_CSI="${OUT_BCF}.csi"
+METADATA_TSV="${OUT_DIR}/${OUT_NAME%.bcf}.metadata.tsv"
+
+if [[ -z "$BPM" || -z "$EGT" || -z "$FASTA" ]]; then
+    echo "ERROR: --bpm, --egt, and --fasta are required." >&2
+    usage >&2
+    exit 2
+fi
+
+for f in "$BPM" "$EGT" "$FASTA" "${FASTA}.fai"; do
+    if [[ ! -f "$f" ]]; then
+        echo "ERROR: missing input: $f" >&2
+        exit 1
+    fi
+done
+
+if ! command -v bcftools >/dev/null 2>&1; then
+    echo "ERROR: bcftools not on PATH." >&2
+    exit 1
+fi
+
+if ! bcftools plugin -l 2>/dev/null | grep -q gtc2vcf; then
+    echo "ERROR: bcftools is missing the gtc2vcf plugin. Use the bcftools:1.23-1 image." >&2
+    exit 1
+fi
+
+# Mirrors src/popgen_genotyping/jobs/gtc_to_bcfs_job.py — same flags + norm
+# + sort — but with no `--gtcs`, so the resulting record set has only INFO
+# (and an empty FORMAT column, which `view -G` strips defensively).
+echo "Running bcftools +gtc2vcf (no GTCs) -> norm -> sort -> view -G | bcf ..."
+bcftools +gtc2vcf \
+    --no-version \
+    --do-not-check-bpm \
+    --bpm "$BPM" \
+    --egt "$EGT" \
+    --fasta-ref "$FASTA" \
+    --extra "$METADATA_TSV" \
+  | bcftools norm -m -both --no-version -c x -f "$FASTA" \
+  | bcftools sort -T "$(mktemp -d)" \
+  | bcftools view -G -O b -o "$OUT_BCF" --write-index=csi
+
+echo "Wrote ${OUT_BCF}"
+echo "Wrote ${OUT_CSI}"
+echo "Wrote ${METADATA_TSV}"
+
+if ! bcftools view -h "$OUT_BCF" | grep -q 'ID=GenTrain_Score'; then
+    echo "WARNING: INFO/GenTrain_Score not declared in header — verify gtc2vcf output." >&2
+fi

references.py

@@ -878,12 +878,19 @@ def is_folder(self) -> bool:
             GCA_000001405_15_GRCh38_no_alt_analysis_set_fna='GCA_000001405.15_GRCh38_no_alt_analysis_set.fna',
             GCA_000001405_15_GRCh38_no_alt_analysis_set_fna_fai='GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.fai',
             GDA_8v1_0_D1_ClusterFile_egt='GDA-8v1-0_D1_ClusterFile.egt',
+            # Per-variant INFO from the GDA-8v1-0_D1 EGT cluster file,
+            # This is the "expected" cluster geometry
+            # (GenTrain_Score, Cluster_Sep, N_AA/AB/BB, meanTHETA_*, devTHETA_*,
+            # meanR_*, devR_*) used by the popgen-genotyping SNP QC report as the
+            # reference half of every observed-vs-EGT comparison.
+            GDA_8v1_0_D1_ClusterFile_egt_info_bcf='GDA-8v1-0_D1_ClusterFile_info.bcf',
+            GDA_8v1_0_D1_ClusterFile_egt_info_bcf_index='GDA-8v1-0_D1_ClusterFile_info.bcf.csi',
             GDA_8v1_0_D2_bpm='GDA-8v1-0_D2.bpm',
             # Biallelic SNV BED derived from the BPM, polarised against the
             # GRCh38 fasta. Produced by
             # illumina_microarray/run_generate_bed_from_illumina_bpm.py.
             GDA_8v1_0_D2_biallelic_snps_bed='GDA-8v1-0_D2-biallelic-snps.bed',
         )
 
-    )
+    ),
 ]