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. Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv
This commit is contained in:
1 parent
60a7ae5ada
commit
9f36f6443b
1 file changed
+39
-13
+39
-13
@@ -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
|
||||
|
||||
Reference in new issue
Block a user