kemal
aa063ccfae
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.
123 lines
6.4 KiB
Bash
123 lines
6.4 KiB
Bash
#!/bin/bash
|
||
# Phase P8 scientific validation on GEO GSE165845 (GSA-MD v1, 120 1000 Genomes individuals × 3 replicates):
|
||
# 1. fetch IDATs (GEO), GSA v1 C2 build38 BPM/CSV + C1 EGT, GRCh38 FASTA, NYGC 30x per-chromosome VCFs, GIAB HG001;
|
||
# 2. `simurg genotype` → `simurg gtc-to-vcf --merge`;
|
||
# 3. truth subsets at the array positions (bcftools, split + left-normalised like the array records);
|
||
# 4. scripts/p8/validate_p8.py (plan §16.5 metrics).
|
||
# Args: <unused data dir> <results dir>. Run with OMNI=0 (the Omni2.5 archive is not needed).
|
||
set -uo pipefail
|
||
OUT=$2
|
||
W=/work/p8; mkdir -p $W $OUT/p8
|
||
R=$OUT/p8
|
||
N=$(nproc)
|
||
export DEBIAN_FRONTEND=noninteractive
|
||
apt-get install -y -q bcftools tabix > /dev/null
|
||
bcftools --version | head -1
|
||
|
||
GSA=https://webdata.illumina.com/downloads/productfiles/global-screening-array/v1-0
|
||
NYGC=https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV
|
||
GIAB=https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/NA12878_HG001/NISTv4.2.1/GRCh38/HG001_GRCh38_1_22_v4.2.1_benchmark
|
||
CHROMS="$(seq 1 22) X"
|
||
|
||
echo "=== downloads $(date -u +%T) ==="
|
||
mkdir -p $W/nygc $W/gsa $W/giab
|
||
( 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=$!
|
||
cd $W/gsa
|
||
for z in infinium-global-screening-array-24-v1-0-c2-manifest-file-bpm-build38.zip \
|
||
infinium-global-screening-array-24-v1-0-c2-manifest-file-csv-build38.zip \
|
||
infinium-global-screening-array-v1-0-c1-cluster-file.zip; do
|
||
curl -sSfL --retry 5 -O "$GSA/$z"
|
||
done
|
||
sha256sum -c - <<'EOF' || { echo "!! GSA product file checksum mismatch"; exit 1; }
|
||
10ba8e692f6f14f9c50b3bed00d92a96af7ef8aa873df115abc52b5ae7c209b8 infinium-global-screening-array-24-v1-0-c2-manifest-file-bpm-build38.zip
|
||
fed9af6de6a8f000b42f44ff1aeecc9fde2a4ae29c834c542bedb2f28879a815 infinium-global-screening-array-24-v1-0-c2-manifest-file-csv-build38.zip
|
||
8bd33c286feded0e77b93cacc8533c27dd879a06e97ef9a575dff2037348cb67 infinium-global-screening-array-v1-0-c1-cluster-file.zip
|
||
EOF
|
||
for z in *.zip; do unzip -oq $z; done
|
||
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/
|
||
|
||
echo "=== genotype $(date -u +%T) ==="
|
||
GT=$W/gt
|
||
/usr/bin/time -v .venv/bin/simurg genotype --bpm $BPM --egt $EGT --idats $W/idats --out $GT \
|
||
--threads "$N" --quiet 2> $R/genotype_time.txt || echo "genotype exit $?"
|
||
tail -22 $R/genotype_time.txt
|
||
cp $GT/sample_summary.tsv $GT/relatedness.tsv $GT/controls_qc.tsv $GT/run_manifest.json $R/ 2>/dev/null
|
||
cp -r $GT/logs $R/genotype_logs 2>/dev/null
|
||
|
||
wait $FASTA_PID; tail -2 $W/fasta.log
|
||
FA=data/p6/GRCh38.primary.fa
|
||
echo "=== gtc-to-vcf $(date -u +%T) ==="
|
||
/usr/bin/time -v .venv/bin/simurg gtc-to-vcf --bpm $BPM --csv $CSV --fasta $FA --gtcs $GT --out $W/vcf --merge \
|
||
--threads "$N" --quiet 2> $R/vcf_time.txt || echo "gtc-to-vcf exit $?"
|
||
tail -22 $R/vcf_time.txt
|
||
cp $W/vcf/run_manifest.json $R/vcf_run_manifest.json; cp -r $W/vcf/logs $R/vcf_logs 2>/dev/null
|
||
ARR=$W/vcf/cohort.snv.bcf
|
||
[[ -f $ARR ]] || { echo "!! no cohort BCF"; exit 1; }
|
||
|
||
echo "=== truth subsets $(date -u +%T) ==="
|
||
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 FA
|
||
echo $CHROMS | tr ' ' '\n' | xargs -P "$N" -I{} sh -c \
|
||
'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
|
||
echo "array $(bcftools view -H $W/array.norm.bcf | wc -l) records, truth $(bcftools view -H $W/truth.bcf | wc -l)," \
|
||
"giab $(bcftools view -H $W/giab.bcf | wc -l)"
|
||
|
||
echo "=== validate $(date -u +%T) ==="
|
||
/usr/bin/time -v .venv/bin/python scripts/p8/validate_p8.py --array $W/array.norm.bcf --truth $W/truth.bcf \
|
||
--giab $W/giab.bcf --giab-bed $W/giab/hg001.bed --map $W/idats/samples_map.tsv --gt-dir $GT --out $R \
|
||
2> $R/validate_time.txt || echo "validation exit $?"
|
||
tail -22 $R/validate_time.txt
|
||
gzip -f $R/replicate_pairs.tsv $R/systematic_sites.tsv 2>/dev/null
|
||
echo "=== task done $(date -u +%T) ==="
|