Data Preparation Example

This example covers Stage 1 of the GRiD workflow: fetching reference inputs, parsing locus coordinates, streaming regional CRAM slices from 1000 Genomes, and building sample manifests.

It sets up the necessary data directory hierarchy and exports coordinate variables to data/pipeline_vars.sh so downstream steps (IBS computation and GRiD configuration) operate on identical locus bounds.

01-prep-data.sh — Stage 1: Reference, CRAM streaming, and manifest prep
  1#!/bin/bash
  2#SBATCH --job-name=grid_prep
  3#SBATCH --output=slurm/%x.out
  4#SBATCH --error=slurm/%x.err
  5#SBATCH --time=4-00:00:00
  6#SBATCH --cpus-per-task=8
  7#SBATCH --mem=200G
  8
  9# Stage 1 of 3 — download reference/panel/regions/repeat-mask, stream LPA
 10# CRAMs for all 1000G samples, and build the GRiD-facing samples file.
 11#
 12# Submit this first. Writes data/pipeline_vars.sh, which stages 2 and 3
 13# source to pick up the LPA locus coordinates without re-deriving them.
 14#
 15# Usage: sbatch 01_prep_data.sh   (run from the shared working directory —
 16# stages 2 and 3 must be submitted from this same directory)
 17
 18set -euo pipefail
 19
 20module load samtools
 21
 22WORK_DIR="$(pwd)"
 23CRAM_DIR="$WORK_DIR/crams"
 24LOG_DIR="$WORK_DIR/logs"
 25DATA_DIR="$WORK_DIR/data"
 26
 27mkdir -p "$CRAM_DIR" "$LOG_DIR" "$DATA_DIR" "$WORK_DIR/slurm"
 28
 29THREADS="${SLURM_CPUS_PER_TASK:-$(nproc)}"
 30# Cap concurrent EBI connections — higher concurrency has caused intermittent
 31# "could not find directory listing" failures against this mirror before.
 32STREAM_JOBS=$(( THREADS > 4 ? 4 : THREADS ))
 33
 34BASE_URL="https://ftp.1000genomes.ebi.ac.uk/vol1/ftp"
 35PANEL_URL="${BASE_URL}/release/20130502/integrated_call_samples_v3.20130502.ALL.panel"
 36REF_URL="https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz"
 37REGIONS_FILE_URL="https://raw.githubusercontent.com/caterer-z-t/GRiD/main/files/734_possible_coding_vntr_regions.IBD2R_gt_0.25.uniq.txt"
 38REPEAT_MASK_URL="https://raw.githubusercontent.com/alexliyihao/vntrwrap/main/normalize_mosdepth/external_source/repeat_mask_list.hg38.ucsc_bed"
 39
 40timestamp() { date '+%Y-%m-%d %H:%M:%S'; }
 41log() { echo "[$(timestamp)] $*" | tee -a "$LOG_DIR/prep.log"; }
 42
 43download_with_retry() {
 44    local url="$1" out="$2" attempt
 45    for attempt in 1 2 3 4 5; do
 46        if wget -q --tries=1 --timeout=60 -O "$out" "$url"; then
 47            return 0
 48        fi
 49        rm -f "$out"
 50        sleep $(( attempt * 3 + RANDOM % 3 ))
 51    done
 52    echo "ERROR: failed to download $url after 5 attempts" >&2
 53    return 1
 54}
 55
 56log "=== Phase 0: static inputs and locus coordinates ==="
 57
 58if [[ ! -f "$DATA_DIR/regions.txt" ]]; then
 59    log "Downloading VNTR regions file..."
 60    download_with_retry "$REGIONS_FILE_URL" "$DATA_DIR/regions.txt"
 61fi
 62
 63read -r CHR START END < <(awk '$7=="LPA" {print $1, $2, $3; exit}' "$DATA_DIR/regions.txt") || true
 64if [[ -z "${CHR:-}" || -z "${START:-}" || -z "${END:-}" ]]; then
 65    echo "ERROR: Could not parse LPA coordinates from regions.txt"
 66    exit 1
 67fi
 68REGION="chr${CHR}:${START}-${END}"
 69FOCAL_BP=$(( (START + END) / 2 ))
 70
 71# Persist derived coordinates so stages 2 and 3 don't have to re-derive
 72# (and can't silently disagree with) these values.
 73cat > "$DATA_DIR/pipeline_vars.sh" << VARS
 74export CHR="$CHR"
 75export START="$START"
 76export END="$END"
 77export REGION="$REGION"
 78export FOCAL_BP="$FOCAL_BP"
 79VARS
 80
 81log "Locus: LPA KIV-2 ($REGION, hg38, focal_bp=$FOCAL_BP)"
 82log "Wrote shared coordinates to $DATA_DIR/pipeline_vars.sh"
 83
 84REPEAT_MASK="$DATA_DIR/repeat_mask.bed"
 85if [[ ! -f "$REPEAT_MASK" ]]; then
 86    log "Downloading repeat mask file..."
 87    download_with_retry "$REPEAT_MASK_URL" "$REPEAT_MASK"
 88fi
 89
 90log "=== Phase 1: reference genome ==="
 91REF_FA="$DATA_DIR/GRCh38_no_alt.fa"
 92if [[ ! -f "$REF_FA" ]]; then
 93    log "Downloading GRCh38 reference genome..."
 94    download_with_retry "$REF_URL" "$DATA_DIR/GRCh38_no_alt.fa.gz"
 95    gunzip "$DATA_DIR/GRCh38_no_alt.fa.gz"
 96    samtools faidx "$REF_FA"
 97    log "Reference genome ready: $REF_FA"
 98else
 99    log "Reference genome already present — skipping download."
100fi
101
102log "=== Phase 2: sample list ==="
103PANEL="$DATA_DIR/1000G_panel.txt"
104SAMPLES_FILE="$DATA_DIR/streaming_manifest.txt"
105if [[ ! -f "$PANEL" ]]; then
106    log "Downloading 1000G panel file..."
107    download_with_retry "$PANEL_URL" "$PANEL"
108fi
109awk 'NR>1 {print $1, $2}' "$PANEL" > "$SAMPLES_FILE"
110TOTAL=$(wc -l < "$SAMPLES_FILE")
111log "Streaming manifest built: $TOTAL samples"
112
113log "=== Phase 3: stream LPA-region CRAMs (high-coverage, 30x) ==="
114FAILED_LOG="$LOG_DIR/failed_samples.txt"
115: > "$FAILED_LOG"
116
117stream_cram() {
118    local sample="$1" pop="$2" ref="$3" region="$4" out_dir="$5"
119    local out="${out_dir}/${sample}.cram"
120
121    if [[ -f "$out" && -f "${out}.crai" ]]; then
122        return 0
123    fi
124
125    local dir_url="https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000_genomes_project/data/${pop}/${sample}/alignment/"
126
127    local listing="" attempt
128    for attempt in 1 2 3 4 5; do
129        listing=$(wget -qO- --tries=1 --timeout=30 "$dir_url" 2>/dev/null || true)
130        [[ -n "$listing" ]] && break
131        sleep $(( attempt * 3 + RANDOM % 3 ))
132    done
133
134    if [[ -z "$listing" ]]; then
135        echo "ERROR: failed to fetch directory listing for $sample after 5 attempts: $dir_url" >&2
136        echo "$sample" >> "$FAILED_LOG"
137        return 1
138    fi
139
140    local full_filename
141    full_filename=$(echo "$listing" | grep -oE "${sample}\.alt_bwamem_GRCh38DH\.[0-9]+\.${pop}\.low_coverage\.cram" | head -n 1)
142
143    if [[ -z "$full_filename" ]]; then
144        echo "ERROR: sample $sample has a listing but no matching high_coverage CRAM at $dir_url" >&2
145        echo "$sample" >> "$FAILED_LOG"
146        return 1
147    fi
148
149    local url="${dir_url}${full_filename}"
150
151    # Each worker gets its own scratch directory so:
152    #  (a) htslib's remote-index caching (which can drop a .crai named after
153    #      the remote URL into the current working directory) can't collide
154    #      with other concurrent workers or pollute $WORK_DIR, and
155    #  (b) a stalled/hung network connection is bounded by an explicit
156    #      timeout instead of silently occupying a worker slot forever —
157    #      the previous version had NO timeout on this step (only on the
158    #      earlier directory-listing wget), which is the likely reason most
159    #      samples in this run never completed AND never got logged as
160    #      failed: a stuck pipe neither succeeds nor errors.
161    local scratch
162    scratch=$(mktemp -d "${out_dir}/.scratch_${sample}.XXXXXX")
163
164    (
165        cd "$scratch" &&
166        timeout 600 bash -c "samtools view -T '$ref' -b '$url' '$region' | samtools sort -@ 2 -o '${sample}.cram'" &&
167        samtools index "${sample}.cram"
168    )
169    local status=$?
170
171    if [[ $status -eq 0 && -f "$scratch/${sample}.cram" && -f "$scratch/${sample}.cram.crai" ]]; then
172        mv "$scratch/${sample}.cram" "$out"
173        mv "$scratch/${sample}.cram.crai" "${out}.crai"
174        rm -rf "$scratch"
175        return 0
176    fi
177
178    rm -rf "$scratch"
179    if [[ $status -eq 124 ]]; then
180        echo "ERROR: streaming timed out (600s) for $sample" >&2
181    else
182        echo "ERROR: streaming/sort/index failed for $sample (exit $status)" >&2
183    fi
184    echo "$sample" >> "$FAILED_LOG"
185    return 1
186}
187export -f stream_cram
188
189xargs -L 1 -P "$STREAM_JOBS" bash -c '
190    ref="$1"; region="$2"; out_dir="$3"; sample="$4"; pop="$5"
191    stream_cram "$sample" "$pop" "$ref" "$region" "$out_dir"
192' _ "$REF_FA" "$REGION" "$CRAM_DIR" < "$SAMPLES_FILE" >> "$LOG_DIR/stream_crams.log" 2>&1 || true
193
194N_ACTUAL=$(find "$CRAM_DIR" -maxdepth 1 -name "*.cram" | wc -l | tr -d ' ')
195N_FAILED=$(wc -l < "$FAILED_LOG" | tr -d ' ')
196N_UNACCOUNTED=$(( TOTAL - N_ACTUAL - N_FAILED ))
197log "Streaming complete: $N_ACTUAL CRAM files in $CRAM_DIR, $N_FAILED logged failures, $N_UNACCOUNTED unaccounted for"
198if [[ "$N_UNACCOUNTED" -ne 0 ]]; then
199    log "WARNING: unaccounted-for samples means something neither completed nor logged a failure — check $LOG_DIR/stream_crams.log"
200fi
201
202# Clean up any stray remote-index cache files htslib may have dropped in the
203# working directory during remote CRAM streaming.
204find "$WORK_DIR" -maxdepth 1 \( -name "*.crai" -o -name "*.tbi" \) -delete
205
206log "=== Phase 4: build GRiD samples file (one sample ID per line) ==="
207GRID_SAMPLES_FILE="$DATA_DIR/grid_samples.txt"
208: > "$GRID_SAMPLES_FILE"
209for cram_path in "$CRAM_DIR"/*.cram; do
210    [[ -f "$cram_path" ]] || continue
211    filename=$(basename "$cram_path")
212    echo "${filename%.cram}" >> "$GRID_SAMPLES_FILE"
213done
214
215N_GRID_SAMPLES=$(wc -l < "$GRID_SAMPLES_FILE" | tr -d ' ')
216if [[ "$N_GRID_SAMPLES" -eq 0 ]]; then
217    echo "ERROR: no CRAM files found in $CRAM_DIR — nothing to pass to GRiD."
218    exit 1
219fi
220log "GRiD samples file ready: $GRID_SAMPLES_FILE ($N_GRID_SAMPLES samples)"
221
222log "=== Prep complete ==="
223log "Next: submit 02_run_ibs.sh, then 03_generate_config_and_run.sh"

What This Script Does

  1. Locus Coordinates: Downloads the coding VNTR regions file, extracts LPA coordinates, and writes data/pipeline_vars.sh.

  2. Reference Genome: Downloads and indexes the GRCh38 no-alt reference genome via samtools faidx.

  3. Sample Manifest: Fetches the 1000 Genomes panel file and extracts sample IDs.

  4. CRAM Streaming: Parallel-streams LPA regional slices directly from EBI using samtools view and indexes them locally into crams/.

  5. GRiD Manifest: Writes data/grid_samples.txt containing valid, non-empty sample IDs for processing.

Usage

Submit to SLURM:

sbatch examples/01-prep-data.sh

Or execute locally:

bash examples/01-prep-data.sh