kemal
3291609577
scripts/p8/prepare.py lays out the 360 GEO IDATs under their header barcode/position; scripts/p8/validate_p8.py computes concordance, NRD, no-call rate, per-class and per-MAF sensitivity/specificity/precision/F1 vs NYGC 30x and GIAB HG001, replicate concordance, systematic/allele-swap site errors, per-chip batch table and duplicate detection. Cloud userdata can skip the Omni2.5 archive (OMNI=0).
97 lines
5.4 KiB
Bash
97 lines
5.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 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=$!
|
||
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
|
||
|
||
/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) ==="
|
||
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
|
||
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 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) ==="
|