61 lines
3.2 KiB
Bash
61 lines
3.2 KiB
Bash
#!/bin/bash
|
|
# v2 CN/LOH validation on GSE165845 (360 GSA arrays, 120 1000 Genomes individuals x 3): simurg cn call vs the NYGC
|
|
# 30x structural-variant callset (1KGP_3202 GATK-SV/svtools freeze V3) and between replicates.
|
|
# Args: <unused data dir> <results dir>. Run with OMNI=0.
|
|
set -uo pipefail
|
|
OUT=$2
|
|
W=/work/cnv; R=$OUT/cnv; mkdir -p $W $R
|
|
N=$(nproc)
|
|
export DEBIAN_FRONTEND=noninteractive
|
|
apt-get install -y -q bcftools tabix gcc > /dev/null
|
|
curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y --profile minimal > /dev/null 2>&1
|
|
export PATH="$HOME/.cargo/bin:$PATH"
|
|
VIRTUAL_ENV=$PWD/.venv uv pip install -q ./rust/simurg-core || echo "!! kernels not built (pure Python)"
|
|
|
|
SV=https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20210124.SV_Illumina_Integration/1KGP_3202.gatksv_svtools_novelins.freeze_V3.wAF.vcf.gz
|
|
( curl -sSfL --retry 5 -o $W/sv.vcf.gz $SV && curl -sSfL --retry 5 -o $W/sv.vcf.gz.tbi $SV.tbi \
|
|
&& echo "SV callset $(du -h $W/sv.vcf.gz | cut -f1) at $(date -u +%T)" ) > $W/sv.log 2>&1 &
|
|
SV_PID=$!
|
|
( bash scripts/fetch_data.sh p6 ) > $W/fasta.log 2>&1 &
|
|
FASTA_PID=$!
|
|
|
|
echo "=== data $(date -u +%T) ==="
|
|
GSA=https://webdata.illumina.com/downloads/productfiles/global-screening-array/v1-0
|
|
mkdir -p $W/gsa && cd $W/gsa
|
|
for z in infinium-global-screening-array-24-v1-0-c2-manifest-file-bpm-build38.zip \
|
|
infinium-global-screening-array-v1-0-c1-cluster-file.zip; do
|
|
curl -sSfL --retry 5 -O "$GSA/$z" && unzip -oq $z
|
|
done
|
|
cd - > /dev/null
|
|
BPM=$W/gsa/GSA-24v1-0_C2.bpm; EGT=$W/gsa/GSA-24v1-0_C1_ClusterFile.egt
|
|
.venv/bin/python scripts/p8/prepare.py --tsv scripts/p8/gse165845_samples.tsv --dl $W/geo --out $W/idats --workers 24 \
|
|
| tail -2 || exit 1
|
|
|
|
echo "=== genotype $(date -u +%T) ==="
|
|
.venv/bin/simurg genotype --bpm $BPM --egt $EGT --idats $W/idats --out $W/gt --threads $N --quiet
|
|
cp $W/idats/samples_map.tsv $W/gt/sample_summary.tsv $R/
|
|
|
|
echo "=== model $(date -u +%T) ==="
|
|
wait $FASTA_PID; tail -1 $W/fasta.log
|
|
/usr/bin/time -v .venv/bin/simurg cn build-model --bpm $BPM --fasta data/p6/GRCh38.primary.fa --egt $EGT \
|
|
--reference-gtcs $W/gt --max-reference 120 --out $W/gsa1.simurg-cyto 2> $R/model_time.txt | grep -v "GC windows"
|
|
grep -E "Elapsed|Maximum resident" $R/model_time.txt
|
|
|
|
echo "=== cn call $(date -u +%T) ==="
|
|
/usr/bin/time -v .venv/bin/simurg cn call --bpm $BPM --model $W/gsa1.simurg-cyto --gtcs $W/gt \
|
|
--fasta data/p6/GRCh38.primary.fa --out $W/cn --threads $N --quiet 2> $R/call_time.txt
|
|
grep -E "Elapsed|Maximum resident" $R/call_time.txt
|
|
.venv/bin/simurg cn annotate --vcfs $W/cn --quiet || echo "annotate failed"
|
|
cp $W/cn/cn_sample_summary.tsv $W/cn/cn_call_manifest.json $R/
|
|
|
|
echo "=== truth $(date -u +%T) ==="
|
|
wait $SV_PID; cat $W/sv.log
|
|
cut -f3 $W/idats/samples_map.tsv | tail -n +2 | sort -u > $W/inds.txt
|
|
bcftools view -S $W/inds.txt --force-samples -i 'INFO/SVTYPE="DEL" || INFO/SVTYPE="DUP"' -Ou $W/sv.vcf.gz \
|
|
| bcftools query -i 'GT="alt"' -f '[%CHROM\t%POS\t%INFO/END\t%INFO/SVTYPE\t%SAMPLE\n]' > $W/truth.tsv
|
|
wc -l < $W/truth.tsv
|
|
.venv/bin/python scripts/v2/cn_validate_truth.py --calls $W/cn --map $W/idats/samples_map.tsv --truth $W/truth.tsv \
|
|
--bpm $BPM --out $R/validation
|
|
gzip -c $W/truth.tsv > $R/truth.tsv.gz
|
|
echo "=== task done $(date -u +%T) ==="
|