diff --git a/illumina_microarray/copy_egt_info_bcf_to_main.sh b/illumina_microarray/copy_egt_info_bcf_to_main.sh new file mode 100755 index 0000000..dcf132b --- /dev/null +++ b/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 diff --git a/illumina_microarray/extract_egt_info_bcf.sh b/illumina_microarray/extract_egt_info_bcf.sh new file mode 100755 index 0000000..eeabe70 --- /dev/null +++ b/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 <&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 diff --git a/references.py b/references.py index 928e335..684e041 100644 --- a/references.py +++ b/references.py @@ -878,6 +878,13 @@ 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 @@ -885,5 +892,5 @@ def is_folder(self) -> bool: GDA_8v1_0_D2_biallelic_snps_bed='GDA-8v1-0_D2-biallelic-snps.bed', ) - ) + ), ]