kemal
522429b731
Found by the P10 fresh-VM run: `simurg inspect … | head` raised BrokenPipeError, reported as an internal error, and stopped examples/gsa_demo.sh (set -o pipefail). run_guarded now flushes stdout inside the guard and maps BrokenPipeError to ExitCode.SIGPIPE (141) without a message. p10.sh: SKIP_DOCKER/SKIP_CONDA, clearer digest report; p10t.sh runs build+pip+tutorial only.
71 lines
3.6 KiB
Bash
71 lines
3.6 KiB
Bash
#!/usr/bin/env bash
|
|
# End-to-end SimurgArray run on the public GSA v4 PGx demo set (36 samples):
|
|
# IDAT → GTC + sample QC → HTML QC report → VCF (per sample and cohort BCF) → PLINK / bedGraph / final report.
|
|
#
|
|
# Prerequisites (from the repository root):
|
|
# uv venv && uv pip install -e .
|
|
# bash scripts/fetch_data.sh p0 # demo IDATs + manifest + cluster file (~1.4 GB)
|
|
# bash scripts/fetch_data.sh p6 # GRCh38 FASTA for VCF output (~0.9 GB download)
|
|
#
|
|
# Usage: bash examples/gsa_demo.sh [OUT_DIR] [N_SAMPLES] (defaults: results/gsa_demo, all 36)
|
|
# Takes ~10 min for 36 samples on a Raspberry Pi 5 (2 worker processes); much less on a desktop.
|
|
set -euo pipefail
|
|
cd "$(dirname "$0")/.."
|
|
OUT=${1:-results/gsa_demo}
|
|
N=${2:-0}
|
|
THREADS=${THREADS:-2}
|
|
SIMURG=${SIMURG:-.venv/bin/simurg}
|
|
|
|
D=data/p0/gsa_v4_pgx
|
|
BPM=$D/GSA-PGx-48v4-0_20079540_E2_bpm/GSA-PGx-48v4-0_20079540_E2.bpm
|
|
CSV=$D/GSA-PGx-48v4-0_20079540_E2_csv/GSA-PGx-48v4-0_20079540_E2.csv
|
|
EGT=$D/GSA-PGx-48v4-0_20079540_E2_ClusterFile/GSA-PGx-48v4-0_20079540_E2_ClusterFile.egt
|
|
IDATS=$D/idats_Demo_36/idats_Demo_36
|
|
FASTA=data/p6/GRCh38.primary.fa
|
|
for f in "$BPM" "$CSV" "$EGT" "$IDATS" "$FASTA"; do
|
|
[[ -e $f ]] || { echo "missing $f — run scripts/fetch_data.sh p0 and p6 first"; exit 1; }
|
|
done
|
|
mkdir -p "$OUT"
|
|
|
|
# Optional subset: symlink the first N IDAT pairs into $OUT/idats (same <barcode>/<barcode>_<pos>_<ch>.idat layout).
|
|
if [[ $N -gt 0 ]]; then
|
|
rm -rf "$OUT/idats"
|
|
find "$IDATS" -name '*_Grn.idat' | sort | sed -n "1,${N}p" | while read -r g; do
|
|
b=$(basename "$(dirname "$g")"); mkdir -p "$OUT/idats/$b"
|
|
ln -s "$(realpath "$g")" "$OUT/idats/$b/"; ln -s "$(realpath "${g%_Grn.idat}_Red.idat")" "$OUT/idats/$b/"
|
|
done
|
|
IDATS=$OUT/idats
|
|
fi
|
|
|
|
echo "== 1. inspect the inputs"
|
|
$SIMURG inspect "$BPM" "$EGT" | sed -n 1,20p # sed reads to the end: no SIGPIPE under pipefail
|
|
|
|
echo "== 2. genotype: GTC v5 per sample, sample_summary.tsv, control-probe QC, contamination, relatedness"
|
|
$SIMURG genotype --bpm "$BPM" --egt "$EGT" --csv "$CSV" --idats "$IDATS" --out "$OUT/gtc" \
|
|
--qc-profile genotyping --threads "$THREADS" --force
|
|
column -t -s$'\t' "$OUT/gtc/sample_summary.tsv" | cut -c1-150 | sed -n 1,8p
|
|
|
|
echo "== 3. HTML QC report + QC table"
|
|
$SIMURG qc report --data "$OUT/gtc" --out "$OUT/qc_report" --output-format both
|
|
|
|
echo "== 4. VCF: one bgzipped, tabix-indexed VCF per sample, then one cohort BCF"
|
|
$SIMURG gtc-to-vcf --bpm "$BPM" --csv "$CSV" --fasta "$FASTA" --gtcs "$OUT/gtc" --out "$OUT/vcf" \
|
|
--threads "$THREADS" --force
|
|
$SIMURG gtc-to-vcf --bpm "$BPM" --csv "$CSV" --fasta "$FASTA" --gtcs "$OUT/gtc" --out "$OUT/bcf" --merge \
|
|
--threads "$THREADS" --force
|
|
|
|
echo "== 5. exports"
|
|
$SIMURG export plink --vcfs "$OUT/vcf" --gtcs "$OUT/gtc" --out "$OUT/plink" --format both --force
|
|
$SIMURG export bedgraph --bpm "$BPM" --gtcs "$OUT/gtc" --out "$OUT/bedgraph" --chrom-style ucsc --force
|
|
$SIMURG export finalreport --bpm "$BPM" --csv "$CSV" --egt "$EGT" --gtcs "$OUT/gtc" --out "$OUT/finalreport" --force
|
|
|
|
echo "== 6. differential test: re-call with the older normalization and compare field by field"
|
|
$SIMURG genotype --bpm "$BPM" --egt "$EGT" --idats "$IDATS" --out "$OUT/gtc_norm112" --norm-version 1.1.2 \
|
|
--no-controls --no-cohort-qc --threads "$THREADS" --quiet --force
|
|
$SIMURG compare "$OUT/gtc" "$OUT/gtc_norm112" --bpm "$BPM" --out "$OUT/compare" --title "norm 1.2.0 vs 1.1.2" \
|
|
--threads "$THREADS" --force || true
|
|
|
|
echo
|
|
echo "Done. Open $OUT/qc_report/*.html and $OUT/compare/dashboard.html in a browser."
|
|
find "$OUT" -maxdepth 2 -type f \( -name '*.html' -o -name '*.bcf' -o -name '*.pgen' -o -name '*.bed' \) | sort
|