v2 CN: truth validation script (NYGC 30x SV freeze V3) and cloud task cnv

2026-09-29 01:39:49 +01:00
1 parent 0bf5576049
commit b6e8177e5a
@@ -0,0 +1,60 @@
+#!/bin/bash
+# v2 CN/LOH validation on GSE165845 (360 GSA arrays, 120 1000 Genomes individuals x 3): simurg cn call vs the NYGC
+# 30x structural-variant callset (1KGP_3202 GATK-SV/svtools freeze V3) and between replicates.
+# Args: <unused data dir>  <results dir>. Run with OMNI=0.
+set -uo pipefail
+OUT=$2
+W=/work/cnv; R=$OUT/cnv; mkdir -p $W $R
+N=$(nproc)
+export DEBIAN_FRONTEND=noninteractive
+apt-get install -y -q bcftools tabix gcc > /dev/null
+curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y --profile minimal > /dev/null 2>&1
+export PATH="$HOME/.cargo/bin:$PATH"
+VIRTUAL_ENV=$PWD/.venv uv pip install -q ./rust/simurg-core || echo "!! kernels not built (pure Python)"
+
+SV=https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20210124.SV_Illumina_Integration/1KGP_3202.gatksv_svtools_novelins.freeze_V3.wAF.vcf.gz
+( curl -sSfL --retry 5 -o $W/sv.vcf.gz $SV && curl -sSfL --retry 5 -o $W/sv.vcf.gz.tbi $SV.tbi \
+  && echo "SV callset $(du -h $W/sv.vcf.gz | cut -f1) at $(date -u +%T)" ) > $W/sv.log 2>&1 &
+SV_PID=$!
+( bash scripts/fetch_data.sh p6 ) > $W/fasta.log 2>&1 &
+FASTA_PID=$!
+
+echo "=== data $(date -u +%T) ==="
+GSA=https://webdata.illumina.com/downloads/productfiles/global-screening-array/v1-0
+mkdir -p $W/gsa && cd $W/gsa
+for z in infinium-global-screening-array-24-v1-0-c2-manifest-file-bpm-build38.zip \
+         infinium-global-screening-array-v1-0-c1-cluster-file.zip; do
+  curl -sSfL --retry 5 -O "$GSA/$z" && unzip -oq $z
+done
+cd - > /dev/null
+BPM=$W/gsa/GSA-24v1-0_C2.bpm; EGT=$W/gsa/GSA-24v1-0_C1_ClusterFile.egt
+.venv/bin/python scripts/p8/prepare.py --tsv scripts/p8/gse165845_samples.tsv --dl $W/geo --out $W/idats --workers 24 \
+  | tail -2 || exit 1
+
+echo "=== genotype $(date -u +%T) ==="
+.venv/bin/simurg genotype --bpm $BPM --egt $EGT --idats $W/idats --out $W/gt --threads $N --quiet
+cp $W/idats/samples_map.tsv $W/gt/sample_summary.tsv $R/
+
+echo "=== model $(date -u +%T) ==="
+wait $FASTA_PID; tail -1 $W/fasta.log
+/usr/bin/time -v .venv/bin/simurg cn build-model --bpm $BPM --fasta data/p6/GRCh38.primary.fa --egt $EGT \
+  --reference-gtcs $W/gt --max-reference 120 --out $W/gsa1.simurg-cyto 2> $R/model_time.txt | grep -v "GC windows"
+grep -E "Elapsed|Maximum resident" $R/model_time.txt
+
+echo "=== cn call $(date -u +%T) ==="
+/usr/bin/time -v .venv/bin/simurg cn call --bpm $BPM --model $W/gsa1.simurg-cyto --gtcs $W/gt \
+  --fasta data/p6/GRCh38.primary.fa --out $W/cn --threads $N --quiet 2> $R/call_time.txt
+grep -E "Elapsed|Maximum resident" $R/call_time.txt
+.venv/bin/simurg cn annotate --vcfs $W/cn --quiet || echo "annotate failed"
+cp $W/cn/cn_sample_summary.tsv $W/cn/cn_call_manifest.json $R/
+
+echo "=== truth $(date -u +%T) ==="
+wait $SV_PID; cat $W/sv.log
+cut -f3 $W/idats/samples_map.tsv | tail -n +2 | sort -u > $W/inds.txt
+bcftools view -S $W/inds.txt --force-samples -i 'INFO/SVTYPE="DEL" || INFO/SVTYPE="DUP"' -Ou $W/sv.vcf.gz \
+  | bcftools query -i 'GT="alt"' -f '[%CHROM\t%POS\t%INFO/END\t%INFO/SVTYPE\t%SAMPLE\n]' > $W/truth.tsv
+wc -l < $W/truth.tsv
+.venv/bin/python scripts/v2/cn_validate_truth.py --calls $W/cn --map $W/idats/samples_map.tsv --truth $W/truth.tsv \
+  --bpm $BPM --out $R/validation
+gzip -c $W/truth.tsv > $R/truth.tsv.gz
+echo "=== task done $(date -u +%T) ==="