P8 cloud task: stream NYGC VCFs through bcftools instead of downloading

The download-only phase (EBI ~15 MB/s) left the CPU idle long enough for the
idle alarm to terminate the first run. Streams are now cut to the manifest
positions (±100 bp) and the 120 individuals while genotyping runs.

scripts/cloud/tasks/p8.sh

@@ -21,13 +21,6 @@ CHROMS="$(seq 1 22) X"
 
 echo "=== downloads $(date -u +%T) ==="
 mkdir -p $W/nygc $W/gsa $W/giab
-( for c in $CHROMS; do
-    s=""; [[ $c == X ]] && s=".v2"
-    f=1kGP_high_coverage_Illumina.chr$c.filtered.SNV_INDEL_SV_phased_panel$s.vcf.gz
-    echo "$NYGC/$f $W/nygc/chr$c.vcf.gz"
-  done | xargs -P 12 -n 2 sh -c 'curl -sSfL --retry 5 -o "$1" "$0" || echo "!! failed $0"'
-  echo "NYGC downloaded $(du -sh $W/nygc | cut -f1) at $(date -u +%T)" ) > $W/nygc.log 2>&1 &
-NYGC_PID=$!
 ( for e in vcf.gz vcf.gz.tbi bed; do curl -sSfL --retry 5 -o $W/giab/hg001.$e $GIAB.$e; done ) &
 ( bash scripts/fetch_data.sh p6 ) > $W/fasta.log 2>&1 &
 FASTA_PID=$!
@@ -47,6 +40,42 @@ cd - > /dev/null
 BPM=$W/gsa/GSA-24v1-0_C2.bpm; CSV=$W/gsa/GSA-24v1-0_C2.csv; EGT=$W/gsa/GSA-24v1-0_C1_ClusterFile.egt
 ls -la $W/gsa
 
+# NYGC per-chromosome VCFs (~40 GB, EBI serves ~15 MB/s) are streamed and cut down on the fly to the manifest
+# positions (±100 bp, so left-normalised indels stay inside) and the GSE165845 individuals: no copy on disk, and
+# the CPU stays busy (a download-only phase would trip the idle alarm, which terminates the instance).
+.venv/bin/python - $BPM > $W/targets_chr.tsv <<'PY'
+import sys
+from simurg.formats.bpm import read_bpm
+m = read_bpm(sys.argv[1])
+iv = {}
+for c, p in zip(m.chroms, m.map_infos.tolist()):
+    if p > 0 and c in {str(i) for i in range(1, 23)} | {"X"}:
+        iv.setdefault(c, []).append((max(1, p - 100), p + 100))
+for c in [str(i) for i in range(1, 23)] + ["X"]:
+    cur = None
+    for s, e in sorted(iv.get(c, [])):
+        if cur and s <= cur[1]:
+            cur[1] = max(cur[1], e)
+            continue
+        if cur:
+            print(f"chr{c}\t{cur[0]}\t{cur[1]}")
+        cur = [s, e]
+    if cur:
+        print(f"chr{c}\t{cur[0]}\t{cur[1]}")
+PY
+wc -l < $W/targets_chr.tsv
+for c in $(seq 1 22) X Y M; do echo -e "chr$c\t$c"; done > $W/chrmap.txt
+cut -f2 scripts/p8/gse165845_samples.tsv | tail -n +2 | sort -u | paste -sd, > $W/inds.txt
+export W NYGC
+( for c in $CHROMS; do echo $c; done | xargs -P 23 -I{} sh -c \
+    's=""; [ {} = X ] && s=".v2"
+     curl -sSfL --retry 5 "$NYGC/1kGP_high_coverage_Illumina.chr{}.filtered.SNV_INDEL_SV_phased_panel$s.vcf.gz" \
+     | bcftools view -T $W/targets_chr.tsv -s $(cat $W/inds.txt) --force-samples -Ou \
+     | bcftools annotate --rename-chrs $W/chrmap.txt -Ob -o $W/nygc/raw.{}.bcf 2> $W/nygc/{}.log \
+     && echo "chr{} done $(date -u +%T)" || echo "!! chr{} failed"'
+  echo "NYGC streamed at $(date -u +%T)" ) > $W/nygc.log 2>&1 &
+NYGC_PID=$!
+
 /usr/bin/time -v .venv/bin/python scripts/p8/prepare.py --tsv scripts/p8/gse165845_samples.tsv --dl $W/geo \
   --out $W/idats --workers 24 2> $R/prepare_time.txt || { tail -30 $R/prepare_time.txt; exit 1; }
 cp $W/idats/samples_map.tsv $R/
@@ -70,17 +99,14 @@ ARR=$W/vcf/cohort.snv.bcf
 [[ -f $ARR ]] || { echo "!! no cohort BCF"; exit 1; }
 
 echo "=== truth subsets $(date -u +%T) ==="
-for c in $(seq 1 22) X Y M; do echo -e "chr$c\t$c"; done > $W/chrmap.txt
-cut -f3 $W/idats/samples_map.tsv | tail -n +2 | sort -u | paste -sd, > $W/inds.txt
 bcftools query -f '%CHROM\t%POS\n' $ARR | awk '{print "chr"$1"\t"$2}' | uniq > $W/sites_chr.tsv
 wc -l < $W/sites_chr.tsv
 bcftools norm -f $FA -m-any --check-ref w -Ob -o $W/array.norm.bcf $ARR 2> $R/norm_array.txt
 wait $NYGC_PID; cat $W/nygc.log; df -h /work | tail -1
-export W FA
+export FA
 echo $CHROMS | tr ' ' '\n' | xargs -P "$N" -I{} sh -c \
-  'bcftools view -T $W/sites_chr.tsv -s $(cat $W/inds.txt) --force-samples -Ou $W/nygc/chr{}.vcf.gz \
+  'bcftools norm -f $FA -m-any --check-ref w -Ob -o $W/truth.{}.bcf $W/nygc/raw.{}.bcf 2> $W/truth.{}.log \
-   | bcftools annotate --rename-chrs $W/chrmap.txt -Ou \
+   || echo "!! truth chr{} failed"'
-   | bcftools norm -f $FA -m-any --check-ref w -Ob -o $W/truth.{}.bcf 2> $W/truth.{}.log || echo "!! truth chr{} failed"'
 bcftools concat -Ob -o $W/truth.bcf $(for c in $CHROMS; do echo $W/truth.$c.bcf; done)
 bcftools view -T $W/sites_chr.tsv -Ou $W/giab/hg001.vcf.gz | bcftools annotate --rename-chrs $W/chrmap.txt -Ou \
   | bcftools norm -f $FA -m-any --check-ref w -Ob -o $W/giab.bcf