diff --git a/scripts/cloud/tasks/p8.sh b/scripts/cloud/tasks/p8.sh index 2a28496..4be7354 100644 --- a/scripts/cloud/tasks/p8.sh +++ b/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 annotate --rename-chrs $W/chrmap.txt -Ou \ - | bcftools norm -f $FA -m-any --check-ref w -Ob -o $W/truth.{}.bcf 2> $W/truth.{}.log || echo "!! truth chr{} failed"' + 'bcftools norm -f $FA -m-any --check-ref w -Ob -o $W/truth.{}.bcf $W/nygc/raw.{}.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