run_filtering.sh

#!/bin/bash
# ──────────────────────────────────────────────────────────────────────────
#  End-to-end clinical filtering with Pangolin splice scoring.
#
#  Pass 1: filtering_r.pl emits <proband>.pangolin_input.csv (rare, panel,
#          MANE, whitelisted "structural-pass" variants) for any proband that
#          does not yet have a Pangolin score map.
#  Pangolin: compute de-novo splice scores on those variants (GPU).
#  Parse:    build <proband>.pangolin.tsv (chr-pos-ref-alt -> max|delta|).
#  Pass 2:   filtering_r.pl produces <proband>.<panel>.candidatos using the
#            scores (pangolin_score column + splice rescue: >= 0.2 whitelisted
#            splice consequences, >= 0.5 probes/others; PP3_Supporting >= 0.2).
#  Cleanup:  all Pangolin scratch is deleted; only <proband>.<panel>.candidatos
#            and the annotated VCFs (*.germline.vep.vcf.gz + .tbi + _summary.html)
#            survive. Pangolin is therefore recomputed every run (cheap — only
#            the few hundred structural-pass candidates).
# ──────────────────────────────────────────────────────────────────────────
set -euo pipefail
# Optional arg $1 = genes-of-interest file (one symbol per line) forwarded to
# filtering_r.pl; if omitted the default panel (g4e-2026.txt) is used.
# Working directory holds the *.germline.vep.vcf.gz data + reference files
# (defaults to this script's own directory; override with $WORKDIR).
GENES="${1:-}"
_CF_DIR="$(dirname "$(readlink -f "$0")")"            # resolve before cd
_CF_SITE="$_CF_DIR/site.sh"
# Canonicalize the panel path BEFORE cd — otherwise a relative panel resolves
# against $WORKDIR (or silently falls back to a same-named file in the repo).
[[ -n "$GENES" ]] && GENES="$(readlink -f "$GENES")"
cd "${WORKDIR:-$_CF_DIR}"

# Optional: force specific sample(s) as proband, overriding filename-based
# auto-discovery.  e.g.  PROBAND="FAM002-M" bash run_filtering.sh
#                        PROBAND="FAM002-P FAM002-M" bash run_filtering.sh genes.txt
PROBAND_ARGS=()
for _p in ${PROBAND:-}; do PROBAND_ARGS+=(--proband "$_p"); done
GENES_ARGS=()
[[ -n "$GENES" ]] && GENES_ARGS=(--list "$GENES")   # filtering_r.pl takes the panel via -l/--list only
FWD=("${PROBAND_ARGS[@]+"${PROBAND_ARGS[@]}"}" "${GENES_ARGS[@]+"${GENES_ARGS[@]}"}")

# Configurable environment — all paths come from site.sh (override in an untracked
# site.env; see README section 0). Sourced via the SCRIPT's real directory because we
# have already cd'd into $WORKDIR above.
source "$_CF_SITE"
source "$CONDA_BASE/etc/profile.d/conda.sh"
conda activate "$PANGOLIN_ENV"

FA="$PANGOLIN_FASTA"
DB="$PANGOLIN_DB"

echo "===== Pass 1: emit Pangolin candidate inputs ====="
perl "$_CF_DIR/filtering_r.pl" "${FWD[@]}"

echo "===== Pangolin scoring ====="
shopt -s nullglob
for csv in *.pangolin_input.csv; do
    proband="${csv%.pangolin_input.csv}"
    tsv="$proband.pangolin.tsv"
    # Zero structural-pass variants: write the (legitimately) empty score map
    # directly instead of invoking Pangolin on an empty set.
    if [ "$(wc -l < "$csv")" -le 1 ]; then
        echo "[pangolin] $proband: no variants to score — writing empty score map"
        : > "$tsv"
        continue
    fi
    echo "[pangolin] scoring $proband ($(($(wc -l < "$csv") - 1)) variants) ..."
    pangolin "$csv" "$FA" "$DB" "$proband.pangolin" -c CHROM,POS,REF,ALT
    perl "$_CF_DIR/parse_pangolin.pl" "$proband.pangolin.csv" > "$tsv"
    scored="$(wc -l < "$tsv")"
    echo "[pangolin] -> $tsv ($scored variants scored)"
    # An empty score map must NOT survive into pass 2. filtering_r.pl switches to its
    # final pass on the mere EXISTENCE of this file, so a Pangolin run that exits 0 but
    # scores nothing (missing GPU, bad refs, empty parse) would produce a complete,
    # successful-looking candidatos table with every pangolin_score blank — the splice
    # rescue arm silently dead and BP7 never firing — and cleanup would then delete the
    # evidence. Remove the stub and fail loudly instead.
    if [ "$scored" -eq 0 ] && [ "$(wc -l < "$csv")" -gt 1 ]; then
        rm -f -- "$tsv"
        echo "[pangolin] FATAL: 0 of $(( $(wc -l < "$csv") - 1 )) variants scored for $proband." >&2
        echo "[pangolin]        Removed the empty $tsv so the next run re-scores instead of" >&2
        echo "[pangolin]        silently reporting without splice evidence. Check the GPU/env/refs." >&2
        exit 1
    fi
done

echo "===== Pass 2: final filtering with splice scores ====="
perl "$_CF_DIR/filtering_r.pl" "${FWD[@]}"

# Keep only the final tables + annotated VCFs; drop all regenerable Pangolin
# scratch. Runs only on success (set -e aborts earlier), so a failed run leaves
# intermediates in place for debugging.
echo "===== Cleanup: removing Pangolin intermediates ====="
rm -f -- *.pangolin_input.csv *.pangolin.csv *.pangolin.tsv *.pangolin.md5

echo "===== done ====="