IBS/IBD Neighbor Computation Example

This example demonstrates Stage 2 of the GRiD workflow: using computeIBSpbwt to generate the IBS neighbors file required by GRiD for haploid copy number estimation (Step 5).

The script automatically sources locus coordinates (chromosome and focal midpoint) derived during Stage 1 from data/pipeline_vars.sh.

02-run-ibs.sh — Stage 2: Compute IBS neighbors via PBWT

#!/bin/bash
#SBATCH --job-name=grid_ibs
#SBATCH --output=slurm/%x.out
#SBATCH --error=slurm/%x.err
#SBATCH --time=2-00:00:00
#SBATCH --cpus-per-task=4
#SBATCH --mem=100G

# Stage 2 of 3 — compute IBS neighbors (needed for haploid CN estimation).
# Submit after 01_prep_data.sh completes; sources data/pipeline_vars.sh for
# CHR/FOCAL_BP rather than re-deriving them, so this can't disagree with
# stage 1's locus coordinates.
#
# Usage: sbatch 02_run_ibs.sh   (from the same working directory as stage 1)

set -euo pipefail

module load samtools

NUM_NEIGHBORS=200
while [[ $# -gt 0 ]]; do
    case $1 in
        --num-neighbors) NUM_NEIGHBORS="$2"; shift 2 ;;
        -h|--help)
            echo "Usage: $0 [--num-neighbors N]"
            exit 0 ;;
        *)
            echo "Unknown argument: $1"
            exit 1 ;;
    esac
done

WORK_DIR="${WORK_DIR:-$(pwd)}"
LOG_DIR="$WORK_DIR/logs"
DATA_DIR="$WORK_DIR/data"
OUTPUT_DIR="$WORK_DIR/output"
QC_DIR="$WORK_DIR/qctool"
mkdir -p "$DATA_DIR" "$LOG_DIR" "$OUTPUT_DIR" "$WORK_DIR/slurm" "$QC_DIR"

THREADS="${SLURM_CPUS_PER_TASK:-$(nproc)}"

timestamp() { date '+%Y-%m-%d %H:%M:%S'; }
log() { echo "[$(timestamp)] $*" | tee -a "$LOG_DIR/IBS.log"; }

download_with_retry() {
    local url="$1" out="$2" attempt
    for attempt in 1 2 3 4 5; do
        if wget -q --tries=1 --timeout=60 -O "$out" "$url"; then
            return 0
        fi
        rm -f "$out"
        sleep $(( attempt * 3 + RANDOM % 3 ))
    done
    echo "ERROR: failed to download $url after 5 attempts" >&2
    return 1
}

# --- Load locus coordinates written by 01_prep_data.sh ---
VARS_FILE="$DATA_DIR/pipeline_vars.sh"
if [[ ! -f "$VARS_FILE" ]]; then
    echo "ERROR: $VARS_FILE not found — run 01_prep_data.sh first."
    exit 1
fi
source "$VARS_FILE"

if ! command -v "$WORK_DIR/qctool/qctool" &> /dev/null; then
    log "Installing qctool..."

    download_with_retry \
        "https://zenodo.org/records/22307601/files/qctool.gz?download=1" \
        "$QC_DIR/qctool.gz"

    gunzip -f "$QC_DIR/qctool.gz"

    chmod +x "$WORK_DIR/qctool/qctool"
fi

MISSING=()
for cmd in "$WORK_DIR/qctool/qctool" wget samtools; do
    command -v "$cmd" &>/dev/null || MISSING+=("$cmd")
done
if [[ ${#MISSING[@]} -gt 0 ]]; then
    log "ERROR: missing required tools: ${MISSING[*]}"
    exit 1
fi

COMPUTE_IBS="$WORK_DIR/computeIBSpbwt"
if [[ ! -f "$COMPUTE_IBS" ]]; then
    log "Downloading computeIBSpbwt binary..."
    download_with_retry \
        "https://raw.githubusercontent.com/mhujoel/STRs/main/cpp_files/computeIBSpbwt" \
        "$COMPUTE_IBS"
    chmod +x "$COMPUTE_IBS"
fi

PHASED_VCF_URL="https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV/1kGP_high_coverage_Illumina.chr${CHR}.filtered.SNV_INDEL_SV_phased_panel.vcf.gz"
GENETIC_MAP_URL="https://alkesgroup.broadinstitute.org/Eagle/downloads/tables/genetic_map_hg38_withX.txt.gz"

LPA_VCF="$DATA_DIR/1kGP_chr${CHR}_LPA_phased.vcf.gz"
BGEN_FILE="$DATA_DIR/1kGP_chr${CHR}_LPA_phased.bgen"
SAMPLE_FILE="$DATA_DIR/1kGP_chr${CHR}_LPA_phased.sample"
GENETIC_MAP="$DATA_DIR/genetic_map_hg38_chr${CHR}.txt"
OUTPUT_FILE="$DATA_DIR/ibs_neighbors_chr${CHR}.tsv.gz"

# --- Stream LPA region from the 1000G phased VCF ---
if [[ ! -f "$LPA_VCF" ]]; then
    log "Streaming LPA region from 1000G phased VCF..."
    tabix -h "$PHASED_VCF_URL" "$REGION" | bgzip > "$LPA_VCF"
    tabix -p vcf "$LPA_VCF"
    find "$WORK_DIR" -maxdepth 1 -name "*.tbi" ! -path "${LPA_VCF}.tbi" -delete
    log "LPA phased VCF ready: $LPA_VCF"
else
    log "LPA phased VCF already present — skipping."
fi

# --- Convert phased VCF -> BGEN v1.2 ---
if [[ ! -f "$BGEN_FILE" ]]; then
    log "Converting phased VCF to BGEN v1.2..."
    "$WORK_DIR/qctool/qctool" \
        -g "$LPA_VCF" \
        -og "$BGEN_FILE" \
        -os "$SAMPLE_FILE" \
        -ofiletype bgen_v1.2 \
        -bgen-bits 16 \
        -bgen-permitted-input-rounding-error 0 \
        >> "$LOG_DIR/qctool.log" 2>&1
    log "BGEN file ready: $BGEN_FILE"
else
    log "BGEN file already present — skipping."
fi

# --- Download + decompress + filter Eagle genetic map ---
# The map's chrom column is bare (e.g. "6"), not "chr6" — filtering on
# "chr${CHR}" here would silently match nothing and produce a header-only
# file, so this must stay as the bare chromosome number.
if [[ ! -f "$GENETIC_MAP" ]]; then
    log "Downloading Eagle hg38 genetic map..."
    wget -q -O - "$GENETIC_MAP_URL" \
        | gunzip \
        | awk -v c="${CHR}" 'NR==1 || $1==c' \
        > "$GENETIC_MAP"
    log "Genetic map ready: $GENETIC_MAP"
else
    log "Genetic map already present — skipping."
fi

# --- Run computeIBSpbwt ---
log "Computing IBS neighbors..."
log "  CHR=$CHR  FOCAL_BP=$FOCAL_BP  NEIGHBORS=$NUM_NEIGHBORS  THREADS=$THREADS"

"$COMPUTE_IBS" \
    "$CHR"          \
    "$FOCAL_BP"     \
    "$BGEN_FILE"    \
    "$SAMPLE_FILE"  \
    "$GENETIC_MAP"  \
    "$NUM_NEIGHBORS"\
    "$THREADS"      \
    "$OUTPUT_FILE"

log "IBS neighbors written to: $OUTPUT_FILE"
log "Next: submit 03_generate_config_and_run.sh"

Key Parameters & Environment Variables

The coordinates and variables used by 02-run-ibs.sh are sourced dynamically from Stage 1 (data/pipeline_vars.sh) or can be passed via command line options:

Parameter / Variable	Description
--num-neighbors	Number of IBS neighbors per haplotype (default: 200).
CHR	Target chromosome (e.g., 6, exported by Stage 1).
FOCAL_BP	Base pair position at center of VNTR (e.g., LPA KIV-2 midpoint in hg38, exported by Stage 1).
COMPUTE_IBS	Path to the auto-downloaded or local computeIBSpbwt binary.
PHASED_VCF_URL	1000 Genomes high-coverage phased VCF URL streamed for the target region via tabix.
GENETIC_MAP	Eagle hg38 genetic map auto-downloaded and filtered for chromosome CHR.
OUTPUT_FILE	Compressed output path (data/ibs_neighbors_chr${CHR}.tsv.gz).

Usage

Submit to SLURM after 01-prep-data.sh completes:

sbatch examples/02-run-ibs.sh

Override the default number of neighbors:

sbatch examples/02-run-ibs.sh --num-neighbors 100

Passing Output to GRiD

Stage 3 (03-run-grid.sh) automatically detects data/ibs_neighbors_chr${CHR}.tsv.gz upon completion and embeds the path directly into config.yaml:

compute_haploid_genotypes:
  run: True
  output_file_prefix: "haploid_genotypes"
  ibs_output: "data/ibs_neighbors_chr6.tsv.gz"
  min_neighbors: 1
  max_neighbors: 10
  n_iters: 100
  method: "ibs"