diff --git a/README.md b/README.md index d69a334..ed5ea60 100644 --- a/README.md +++ b/README.md @@ -59,6 +59,12 @@ To browse releases instead, https://zenodo.org/records/13694210/latest opens the RUFUS requires the following data to run: 1) A subject sample in FASTQ/BAM/CRAM/generator format. BAM/CRAM may be unaligned when using whole genome mode; FASTQ is whole-genome only and cannot be combined with `-R/--region`. If the sample is split across several files, pass `-s` once per file — they are treated as one sample, not as separate subjects. 2) At least one of: one or more control samples (`-c`, same formats as the subject), or an exclude hash (`-e`). Multiple distinct controls are supported, e.g. mother and father for a trio. Supplying only `-e` — typically pre-built 1000G/control hashes — is single-sample mode. + +BAM, CRAM and generator inputs may be freely mixed across `-s` and `-c` — they all feed the same read stream internally. FASTQ must stand alone: if any input is a FASTQ, all of them must be. RUFUS filters reads from the FASTQ mate files directly in that case, so reads from any BAM/CRAM/generator alongside them would be counted but never filtered, quietly costing calls. + +A generator file is a shell script that writes SAM to stdout (e.g. a single line `samtools view -h -F 3328 /path/sample.bam`); RUFUS runs it to obtain reads. Generators are whole-genome only — `-R/--region` is not applied to them, so the SLURM launcher rejects them in windowed (`-w`) mode. + +The launcher scans generator files for the paths they reference and bind-mounts those directories automatically, reporting what it added. The scan is best-effort: it can only use paths that appear literally in the file and exist on the host, so a path built at runtime (`$DATA/sample.bam`) or supplied through the environment (samtools' `REF_PATH`/`REF_CACHE` for CRAM decode) will be missed — bind those yourself with `-d`. To keep such gaps from surfacing hours into a queued job, the launcher then runs each generator inside the container and refuses to submit unless it produces SAM, printing the generator's own error output. That check is skipped with a warning if no container runtime is on the submit host's PATH. 3) A reference fasta file (this must be indexed by BWA) - for use in reporting the called variants. *It's recommended to provide the BWA indexes in the same directory as the reference to save time creating them during the RUFUS run.*\ \ To create the BWA indexes, run the following commands: diff --git a/runRufus.sh b/runRufus.sh index 55b80f7..5a9810d 100755 --- a/runRufus.sh +++ b/runRufus.sh @@ -491,6 +491,34 @@ make_jelly_hash () fi } +# Normalize a jellyfish -s size token (e.g. 64G, 500M, 1000000) to the power-of-two slot +# count jellyfish actually allocates, so two sizes compare the way "jellyfish merge" compares +# them. Echoes the normalized integer, or nothing if the token is unparseable. +normalize_hash_size () +{ + local tok="$1" num unit bytes p + [[ "$tok" =~ ^([0-9]+)([GgMmKk]?)$ ]] || { echo ""; return; } + num="${BASH_REMATCH[1]}"; unit="${BASH_REMATCH[2]}" + case "$unit" in + G|g) bytes=$(( num * 1024 * 1024 * 1024 ));; + M|m) bytes=$(( num * 1024 * 1024 ));; + K|k) bytes=$(( num * 1024 ));; + *) bytes=$num;; + esac + p=1 + while [ "$p" -lt "$bytes" ]; do p=$(( p * 2 )); done + echo "$p" +} + +# Read the -s (hash size) a pre-built Jhash was created with, straight from its header via +# "jellyfish info", normalized to the allocated power-of-two. Empty if it can't be determined. +read_built_hash_size () +{ + local hash="$1" tok + tok=$($modifiedJelly info "$hash" 2>/dev/null | sed -n 's/.* -s \([0-9]\+[GMKgmk]\?\) .*/\1/p' | head -1) + [ -n "$tok" ] && normalize_hash_size "$tok" +} + check_empty_hashes () { local region_arg="$1" @@ -758,7 +786,13 @@ fi ProbandExtension="${ProbandFileName##*.}" ProbandGenerator="${ProbandFileName}${region_postfix}.generator" -# Build concatenated generator from all subject files +# Build concatenated generator from all subject files. +# Exactly one command in the generator may emit a SAM header: the body is run as a single +# stream (`bash "$ProbandGenerator" | samtools ...`) and samtools aborts on a second @HD +# mid-stream -- collate discards the whole stream, so Filter sees zero reads. _subj_hdr +# carries the header flag for the first emitting command and is cleared thereafter; the +# FASTQ block below continues the same flag so a bam+fastq mix stays single-headered. +_subj_hdr="-h " > "$ProbandGenerator" for subject in "${_arg_subjects[@]}" do @@ -777,7 +811,8 @@ do _region_exit_reason="missing_subject_bam_index" exit 1 fi - echo "samtools view -h -@ 8 -F 3328 $subject $_arg_region" >> "$ProbandGenerator" + echo "samtools view ${_subj_hdr}-@ 8 -F 3328 $subject $_arg_region" >> "$ProbandGenerator" + _subj_hdr="" elif [[ "$subjectExtension" == "cram" ]] then if [[ ! -e "$subject".crai ]] @@ -791,11 +826,20 @@ do echo "ERROR cram reference not provided for cram input" kill -9 $$ fi - echo "samtools view -h -@ 8 -F 3328 -T $_arg_cramref $subject $_arg_region" >> "$ProbandGenerator" + echo "samtools view ${_subj_hdr}-@ 8 -F 3328 -T $_arg_cramref $subject $_arg_region" >> "$ProbandGenerator" + _subj_hdr="" _arg_ref="$_arg_cramref" elif [[ "$subjectExtension" == "generator" ]] then - cat "$subject" >> "$ProbandGenerator" + # A pre-built generator carries its own header-emitting command. Keep it only if it + # lands first; otherwise strip the header flags as it is appended. + if [ -n "$_subj_hdr" ] + then + cat "$subject" >> "$ProbandGenerator" + else + sed -e 's/^\(samtools view\) -h /\1 /' -e 's/ header$//' "$subject" >> "$ProbandGenerator" + fi + _subj_hdr="" else echo "unknown error during generator generation, killing run with non-zero exit status" kill -9 $$ @@ -803,7 +847,8 @@ do done # FASTQ subject(s): whole-genome only -- unaligned reads cannot be region-scoped. The loop above -# skipped them; build the generator here as one @HD header (on the first file) + unmapped SAM records. +# skipped them; build the generator here as unmapped SAM records, headed by a single @HD if no +# bam/cram subject above has already claimed it (_subj_hdr). if [ ${#_arg_subject_fastqs[@]} -gt 0 ]; then if [ -n "$_arg_region" ]; then echo "ERROR: FASTQ input is whole-genome only and cannot be region-scoped; remove -R/--region (or supply an aligned bam/cram)." @@ -814,10 +859,9 @@ if [ ${#_arg_subject_fastqs[@]} -gt 0 ]; then echo "ERROR: FASTQ subject input requires a reference via -r/--ref." exit 1 fi - _fq_first=1 for fq in "${_arg_subject_fastqs[@]}"; do [ -e "$fq" ] || { echo "FASTQ subject file $fq does not exist; killing run"; kill -9 $$; } - _hdr=""; [ "$_fq_first" -eq 1 ] && _hdr=" header"; _fq_first=0 + _hdr=""; [ -n "$_subj_hdr" ] && _hdr=" header"; _subj_hdr="" if [[ "$fq" == *.gz ]]; then echo "perl $RDIR/scripts/FastqToSam.pl <(zcat $fq)$_hdr" >> "$ProbandGenerator" else @@ -887,8 +931,31 @@ do _arg_ref="$_arg_cramref" elif [[ "$parentExtension" = "generator" ]] then + # The historical name is as a FULL path, and it is load-bearing: + # RunJellyForRUFUS.sh early-returns when $GEN.Jhash exists, so a pre-built control hash + # placed next to the input as .Jhash (e.g. a DSA hash symlinked to + # DSA_SMHT004.1.generator.wg.Jhash) is picked up and jellyfish is skipped entirely. The + # generator body is never executed in that case, which is the point -- the hash already + # exists and the reads it came from may not even be on this filesystem. + # + # A caller may equally supply a pre-scoped generator at that path. Only when neither is + # present does the generator actually have to run, and only then is a runnable copy + # materialised -- in the working directory, not next to the user's input. + # + # The generator is used verbatim: -R/--region is NOT applied to it, exactly as for generator + # subjects above. Scoping a generator to a region is the caller's responsibility. parentGenerator="${parent}${region_postfix}" - ParentGenerators+=("$parentGenerator") + if [ ! -e "$parentGenerator" ] && [ ! -e "${parentGenerator}.Jhash" ] + then + if [[ ! -e "$parent" ]] + then + echo "The control generator file $parent does not exist; killing run with non-zero exit status" + kill -9 $$ + fi + parentGenerator="${parentFileName}${region_postfix}.generator" + cat "$parent" > "$parentGenerator" + fi + ParentGenerators+=("$parentGenerator") fi done ################################################################# @@ -1068,6 +1135,47 @@ done ################################################## +############__PREFLIGHT: HASH SIZE MATCH__################ +# Jellyfish can only merge/diff hashes built at the same -s. The subject hash is built fresh here +# (hours for a whole-genome sample), but pre-built control/DSA/exclude hashes carry a fixed size +# from when they were made. If those disagree with the subject size, "$modifiedJelly merge" aborts +# with "Can't merge hash with different size", leaving an empty HashList that only surfaces much +# later as the misleading "No mutant hashes pulled from fastqs". Catch it here, in seconds, before +# paying for the subject build. +# +# Only pre-built hashes that already exist on disk can mismatch; control generators counted in this +# run are built at the subject size and match by construction, so those (not yet on disk) are skipped. + +# Intended subject hash size -- mirror make_jelly_hash: -hs override, else 16G whole-genome / 1G region. +if [ -n "$_arg_hash_size" ]; then + _subject_hash_size="$_arg_hash_size" +elif [ -z "$_arg_region" ]; then + _subject_hash_size="16G" +else + _subject_hash_size="1G" +fi +_subject_slots=$(normalize_hash_size "$_subject_hash_size") + +_hash_size_mismatch=0 +for _control_hash in $(echo $parentsString) $(echo $parentsExcludeString); do + [ -f "$_control_hash" ] || continue # not-yet-built generator control -> will match by construction + _control_slots=$(read_built_hash_size "$_control_hash") + [ -n "$_control_slots" ] || continue # size unreadable -> don't block the run + if [ "$_control_slots" != "$_subject_slots" ]; then + echo "ERROR: hash size mismatch. The subject hash will be built at -s $_subject_hash_size, but a pre-built control/exclude hash was built at a different -s:" >&2 + echo " $_control_hash" >&2 + _hash_size_mismatch=1 + fi +done +if [ "$_hash_size_mismatch" -ne 0 ]; then + echo "Jellyfish cannot merge hashes of different sizes, so the k-mer subtraction would silently yield zero mutant k-mers." >&2 + echo "Fix: re-run with -hs/--hash_size set to the control's size (e.g. -hs 64G), or rebuild the control(s) at -s $_subject_hash_size." >&2 + _region_exit_reason="hash_size_mismatch" + exit 100 +fi +######################################################## + + ####################__GENERATE_JHASH_FILES_FROM_JELLYFISH__##################### _region_exit_reason="jellyfish_stage" CONTROL_EXIT_FILES=() diff --git a/scripts/RunJellyForRUFUS.sh b/scripts/RunJellyForRUFUS.sh index 208ebf9..de47c8f 100755 --- a/scripts/RunJellyForRUFUS.sh +++ b/scripts/RunJellyForRUFUS.sh @@ -1,5 +1,8 @@ #!/bin/bash set -e +# pipefail so the feeder pipeline below reports a failure in `bash "$GEN"` and not just in the +# samtools stage that terminates it. +set -o pipefail GEN=$1 K=$2 T=$3 @@ -25,7 +28,9 @@ else mkfifo "$FIFO_FQ" # bash "$GEN" | "$RDIR/bin/PassThroughSamCheck" "$GEN.Jelly.chr" > "$FIFO_FQ" & - # samtools fastq validated bit-identical to PassThroughSamCheck for counting (job 16682513); generator now emits -h so the header is present. + # samtools fastq validated bit-identical to PassThroughSamCheck for counting (job 16682513); the + # generator's first command emits -h so the header is present -- and only its first, since + # samtools aborts on a second @HD mid-stream. bash "$GEN" | samtools fastq -@ "$T" - > "$FIFO_FQ" & FEEDER=$! @@ -59,9 +64,22 @@ else exit 2 fi - wait + # The feeder's status is not jellyfish's. If the feeder dies partway -- a truncated stream, + # an unreadable input, a malformed generator -- jellyfish sees a clean EOF on the FIFO and + # reports success over however many reads happened to arrive, so an under-counted hash + # looks identical to a complete one. Reap it explicitly and treat a failure as a tool + # failure (2), discarding the partial hash so a rerun cannot pick it up via the + # skip-if-exists check at the top. + wait "$FEEDER" + feeder_rc=$? set -e + if [ "$feeder_rc" -ne 0 ]; then + rm -f "$GEN.Jhash" + echo "ERROR: read feeder failed (exit $feeder_rc) for $GEN; k-mer counts would be incomplete" >&2 + exit 2 + fi + # A zero exit with no output file means jellyfish died without reporting it. if [ ! -s "$GEN.Jhash" ]; then echo "ERROR: jellyfish count reported success but produced no $GEN.Jhash" >&2 diff --git a/singularity/launch_utilities/arg_parser.sh b/singularity/launch_utilities/arg_parser.sh index 8cf3b4a..0225999 100644 --- a/singularity/launch_utilities/arg_parser.sh +++ b/singularity/launch_utilities/arg_parser.sh @@ -10,8 +10,9 @@ DEFAULT_WG_MEM_PER_JOB="150G" usage() { echo "Usage: $0 [-s subject1,subject2,...] [-c control1,control2,control3...] [-b genome_build] [-a slurm_account] [-p slurm_partition] ...options" echo "Required Arguments:" - echo "-s subject(s) A single subject or comma-delimited array of multiple subject BAM/CRAM files (full paths)" - echo "-c control(s) A single control or comma-delimited array of multiple controls (full paths)" + echo "-s subject(s) A single subject or comma-delimited array of multiple subject BAM/CRAM/FASTQ/generator files (full paths)" + echo "-c control(s) A single control or comma-delimited array of multiple controls (full paths, same formats as -s)" + echo " BAM, CRAM and generator inputs may be mixed; FASTQ inputs must be used alone." echo "-b genome_build The desired genome build; currently only supports GRCh38" echo "-r reference Full path to the reference file matching the genome build" echo "-a slurm_account The account for the slurm job" @@ -35,6 +36,7 @@ usage() { echo "-C cpus_per_call How many cpus to allot to each rufus calling stage job; default 40 for entire genome; 12 for 1MB windows" echo "-d dev_binds Comma-delimited list of host:container bind mounts for dev testing (e.g., /local/runRufus.sh:/opt/RUFUS/runRufus.sh)" echo "-P par_low_cov_threshold Control k-mer count ceiling below which a variant is flagged as low-coverage-parent/inherited (default 7; set to 0 to disable, e.g. when using an assembly as the control)" + echo "-H hash_size jellyfish hash size (-s) for the k-mer count step, e.g. 64G. MUST match the -s any pre-built control/DSA/exclude hash was built with, or the merge fails. Maps to runRufus -hs; RUFUS defaults to 16G whole-genome / 1G windowed if unset." echo "-h help Print usage" echo "" echo "Output files are written to the current working directory." @@ -67,6 +69,7 @@ MEM_PER_JOB="" CPUS_PER_JOB="" PAR_LOW_COV_THRESHOLD_RUFUS_ARG="7" DEV_BIND_MOUNTS_ARG=() +HASH_SIZE_RUFUS_ARG="" # Parse command line options using getopts # @@ -78,7 +81,7 @@ DEV_BIND_MOUNTS_ARG=() # STILL MISWIRED: -h carries a colon ("h:") so it demands an argument. Bare `-h` never reaches # the h) case; it falls to the missing-argument branch, printing "Option -h requires an argument" # before the usage text. Usage still prints, so this is cosmetic. Fix is to drop the colon. -while getopts ":s:c:b:a:p:r:m:w:e:l:q:t:f:x:y:z:M:C:K:G:D:V:d:P:h" opt; do +while getopts ":s:c:b:a:p:r:m:w:e:l:q:t:f:x:y:z:M:C:K:G:D:V:d:P:H:h" opt; do case ${opt} in s) IFS=',' read -r -a SUBJECTS_RUFUS_ARG <<< "$OPTARG" @@ -156,6 +159,13 @@ while getopts ":s:c:b:a:p:r:m:w:e:l:q:t:f:x:y:z:M:C:K:G:D:V:d:P:h" opt; do fi PAR_LOW_COV_THRESHOLD_RUFUS_ARG=$OPTARG ;; + H) + if ! [[ "$OPTARG" =~ ^[0-9]+[GMKgmk]?$ ]]; then + echo "ERROR: -H hash_size must be a jellyfish hash size, e.g. 64G, 500M, or a plain integer." >&2 + exit 1 + fi + HASH_SIZE_RUFUS_ARG=$OPTARG + ;; h) usage ;; @@ -199,37 +209,57 @@ for control in "${CONTROLS_RUFUS_ARG[@]}"; do fi done -# Validate that subject and control files are all the same type (bam, cram, or fastq) +# Validate input file types. +# +# This used to require every subject and control to be the exact same type. That was broader than +# anything downstream actually needs: bam, cram and generator inputs all funnel into the same +# concatenated generator in runRufus.sh and the Filter stage streams that generator, so those three +# are indistinguishable by the time reads are pulled. Only FASTQ has to stand alone -- see below. get_input_type() { case "$1" in *.cram) echo "cram" ;; *.bam) echo "bam" ;; *.fastq.gz|*.fq.gz) echo "fastq" ;; *.fastq|*.fq) echo "fastq" ;; + *.generator) echo "generator" ;; *) echo "unknown" ;; esac } -# Validate type of first subject, then ensure all subjects + controls match -SUBJECT_TYPE=$(get_input_type "${SUBJECTS_RUFUS_ARG[0]}") -if [ "$SUBJECT_TYPE" == "unknown" ]; then - echo "ERROR: subject file ${SUBJECTS_RUFUS_ARG[0]} has an unrecognized file type. Supported types: .bam, .cram, .fastq, .fq, .fastq.gz, .fq.gz" >&2 +HAS_FASTQ_INPUT="false" +HAS_NON_FASTQ_INPUT="false" +HAS_GENERATOR_INPUT="false" +for input in "${SUBJECTS_RUFUS_ARG[@]}" "${CONTROLS_RUFUS_ARG[@]}"; do + case "$(get_input_type "$input")" in + fastq) + HAS_FASTQ_INPUT="true" + ;; + generator) + HAS_GENERATOR_INPUT="true" + HAS_NON_FASTQ_INPUT="true" + ;; + bam|cram) + HAS_NON_FASTQ_INPUT="true" + ;; + *) + echo "ERROR: input file $input has an unrecognized file type. Supported types: .bam, .cram, .generator, .fastq, .fq, .fastq.gz, .fq.gz" >&2 + exit 1 + ;; + esac +done + +# FASTQ must be exclusive. A FASTQ subject populates _arg_fastqA/_arg_fastqB in runRufus.sh, and the +# Filter stage then reads ONLY those two mate files -- the generator holding the bam/cram/generator +# reads is never filtered, even though k-mer counting did span it. That yields a HashList whose +# k-mers have no reads to assemble from: fewer calls, no error message. Output naming desyncs as +# well, since post_process is handed subject[0] while runRufus names off the first non-FASTQ subject. +if [ "$HAS_FASTQ_INPUT" == "true" ] && [ "$HAS_NON_FASTQ_INPUT" == "true" ]; then + echo "ERROR: FASTQ inputs cannot be combined with BAM/CRAM/generator inputs. RUFUS filters reads" >&2 + echo " from the FASTQ mate files alone in that case, so reads from the other inputs would be" >&2 + echo " counted but never filtered, silently costing calls. Pass all inputs as FASTQ, or" >&2 + echo " convert the FASTQ to BAM/generator first. BAM, CRAM and generator inputs may be mixed." >&2 exit 1 fi -for subject in "${SUBJECTS_RUFUS_ARG[@]:1}"; do - SUBJ_TYPE=$(get_input_type "$subject") - if [ "$SUBJ_TYPE" != "$SUBJECT_TYPE" ]; then - echo "ERROR: all subject files must be the same type, but first subject is ${SUBJECT_TYPE} and $subject is ${SUBJ_TYPE}. Please ensure all inputs are either all BAMs, all CRAMs, or all FASTQs." >&2 - exit 1 - fi -done -for control in "${CONTROLS_RUFUS_ARG[@]}"; do - CTRL_TYPE=$(get_input_type "$control") - if [ "$CTRL_TYPE" != "$SUBJECT_TYPE" ]; then - echo "ERROR: all subject and control files must be the same type, but subject is ${SUBJECT_TYPE} and control $control is ${CTRL_TYPE}. Please ensure all inputs are either all BAMs, all CRAMs, or all FASTQs." >&2 - exit 1 - fi -done # Check that reference file exists if [ ! -f "$REFERENCE_RUFUS_ARG" ]; then @@ -290,6 +320,75 @@ else THREAD_LIMIT_RUFUS_ARG=${THREAD_LIMIT_RUFUS_ARG:-10} fi +# Guard the per-job memory against the jellyfish hash floor. jellyfish pre-faults the ENTIRE hash +# array at startup, so the -s size sets a hard RAM floor (not a peak that ramps): if -M is below it +# the count is OOM-killed in the first minute -- a fast, confusing failure that otherwise only shows +# up after submission. Catch it here instead. +# +# _hash_size_to_pow2_gb: parse a jellyfish -s token to the GiB of the power-of-two array jellyfish +# actually rounds up to (so -H 48G is judged as the 64G it allocates, not 48). +_hash_size_to_pow2_gb() { + local tok="$1" num unit bytes p + [[ "$tok" =~ ^([0-9]+)([GgMmKk]?)$ ]] || { echo 0; return; } + num="${BASH_REMATCH[1]}"; unit="${BASH_REMATCH[2]}" + case "$unit" in + G|g) bytes=$(( num * 1024 * 1024 * 1024 ));; + M|m) bytes=$(( num * 1024 * 1024 ));; + K|k) bytes=$(( num * 1024 ));; + *) bytes=$num;; + esac + p=1 + while [ "$p" -lt "$bytes" ]; do p=$(( p * 2 )); done + echo $(( p / (1024 * 1024 * 1024) )) +} +# Approx RAM (GiB) jellyfish -m 25 pre-faults for a hash of the given -s, from measured `jellyfish +# mem` points (16G->64, 32G->123, 64G->237, 128G->~460); 4 GiB per GiB-of-hash elsewhere, which is +# conservative (never lands under the true floor) for the large whole-genome sizes this guards. +_hash_mem_floor_gb() { + local g; g=$(_hash_size_to_pow2_gb "$1") + case "$g" in + 16) echo 64;; 32) echo 123;; 64) echo 237;; 128) echo 460;; + *) echo $(( g * 4 ));; + esac +} +# Parse a SLURM memory string (e.g. 300G, 300000M, 1T) to GiB; 0 if it has no unit we recognize. +_slurm_mem_to_gb() { + local tok="$1" num unit + [[ "$tok" =~ ^([0-9]+)([GgMmTt])$ ]] || { echo 0; return; } + num="${BASH_REMATCH[1]}"; unit="${BASH_REMATCH[2]}" + case "$unit" in T|t) echo $(( num * 1024 ));; G|g) echo "$num";; M|m) echo $(( num / 1024 ));; esac +} + +# Effective count-step hash size: -H override, else RUFUS's own default (16G whole-genome, 1G window). +if [ -n "$HASH_SIZE_RUFUS_ARG" ]; then + _effective_hash_size="$HASH_SIZE_RUFUS_ARG" +elif [ "$WINDOW_SIZE_RUFUS_ARG" -eq 0 ]; then + _effective_hash_size="16G" +else + _effective_hash_size="1G" +fi +_hash_floor_gb=$(_hash_mem_floor_gb "$_effective_hash_size") +_mem_gb=$(_slurm_mem_to_gb "$MEM_PER_JOB") +if [ "$_mem_gb" -gt 0 ] && [ "$_hash_floor_gb" -gt "$_mem_gb" ]; then + echo "ERROR: per-job memory (-M ${MEM_PER_JOB}) is below the jellyfish hash-size floor for -s ${_effective_hash_size}." >&2 + echo " jellyfish pre-faults the whole ${_effective_hash_size} array (~${_hash_floor_gb} GiB) at startup and would be" >&2 + echo " OOM-killed within the first minute of the count step." >&2 + echo " Raise -M to at least $(( _hash_floor_gb + 40 ))G (headroom for reads + spill), or lower -H." >&2 + echo " Whole-genome floors: -s 16G ~64, 32G ~123, 64G ~237 GiB." >&2 + exit 1 +fi + +# Generator inputs are whole-genome only, for the same reason FASTQ is: runRufus.sh cannot scope one +# to a region. A generator is an arbitrary shell command producing SAM, so `-R` is simply not applied +# to it -- the contents are used verbatim. In windowed mode every array task would therefore count +# and call the entire genome, burning thousands of node-hours to produce identical per-window VCFs. +if [ "$HAS_GENERATOR_INPUT" == "true" ] && [ "$WINDOW_SIZE_RUFUS_ARG" -ne 0 ]; then + echo "ERROR: generator inputs cannot be used with windowed mode (-w); they are whole-genome only." >&2 + echo " A generator cannot be region-scoped, so every window would re-run the whole genome." >&2 + echo " Drop -w to run whole-genome, or supply the sample as an indexed BAM/CRAM." >&2 + exit 1 +fi + if [ "$THREAD_LIMIT_RUFUS_ARG" -ge "$CPUS_PER_JOB" ]; then echo "ERROR: thread limit ($THREAD_LIMIT_RUFUS_ARG) must be less than cpus per job ($CPUS_PER_JOB)." >&2 exit 1 @@ -340,6 +439,54 @@ if [ ${#DEV_BIND_MOUNTS_ARG[@]} -gt 0 ]; then echo "DEV MODE: additional bind mounts:${DEV_BIND_ARGS}" fi +# Directories that must never be bind-mounted from the host, including anything beneath them. +# A bind shadows whatever the container has at that path, so binding /opt would hide the entire +# /opt/RUFUS install and binding /usr would swap the container's toolchain for the host's. +# Generators legitimately reference paths under these (e.g. /opt/RUFUS/scripts/FastqToSam.pl), +# so they have to be filtered out silently rather than treated as an error. +BIND_DENYLIST=(/ /bin /boot /dev /etc /lib /lib64 /opt /proc /root /run /sbin /srv /sys /usr /var) + +is_denied_bind_dir() { + local dir="$1" denied + for denied in "${BIND_DENYLIST[@]}"; do + if [ "$denied" == "/" ]; then + [ "$dir" == "/" ] && return 0 + else + [ "$dir" == "$denied" ] && return 0 + [[ "$dir" == "$denied"/* ]] && return 0 + fi + done + return 1 +} + +# Echo the directories referenced *inside* a generator file, one per line. +# +# A generator is a shell script RUFUS executes (`bash `) to produce SAM, so its data +# dependencies live in the file body rather than on the command line -- the launcher would +# otherwise bind the generator itself and none of the data it reads. +# +# Deliberately best-effort and strictly additive: a token is used only if it resolves to something +# that exists on the host. Paths assembled at runtime ("$DATA/sample.bam" yields the non-existent +# "/sample.bam"), globs, and process substitutions are skipped silently, leaving the bind set +# exactly as it was before. Nothing here can turn a working bind set into a broken one, and the +# container preflight in setup_slurm.sh is what catches whatever this misses. +generator_referenced_dirs() { + local gen="$1" + local token resolved dir + while IFS= read -r token; do + [ -n "$token" ] || continue + [ -e "$token" ] || continue + resolved="$(realpath "$token" 2>/dev/null)" || continue + if [ -d "$resolved" ]; then + dir="$resolved" + else + dir="$(dirname "$resolved")" + fi + is_denied_bind_dir "$dir" && continue + echo "$dir" + done < <(grep -o "/[^[:space:]'\";|&<>()\`]*" "$gen" 2>/dev/null | sort -u) +} + # Collect unique parent directories for all input files to use as bind mounts. # Singularity --bind preserves host paths inside the container (no remapping needed). collect_bind_dirs() { @@ -388,12 +535,36 @@ collect_bind_dirs() { fi done + # Generator inputs carry their data references inside the file, so bind those dirs too. + local gen_added=() + for f in "${SUBJECTS_RUFUS_ARG[@]}" "${CONTROLS_RUFUS_ARG[@]}"; do + [ "$(get_input_type "$f")" == "generator" ] || continue + local gd + while IFS= read -r gd; do + [ -n "$gd" ] || continue + if [ -z "${seen_dirs[$gd]+x}" ]; then + seen_dirs["$gd"]=1 + dirs+=("$gd") + gen_added+=("$gd") + fi + done < <(generator_referenced_dirs "$f") + done + # collect_bind_dirs runs inside a command substitution, so this cannot set a flag the caller + # would see; the caller silences the repeat by setting _REPORTED_GEN_BINDS after the first call. + if [ ${#gen_added[@]} -gt 0 ] && [ -z "${_REPORTED_GEN_BINDS:-}" ]; then + echo "Binding directories referenced inside generator input(s): ${gen_added[*]}" >&2 + echo " (best-effort scan; add any it missed with -d host:container)" >&2 + fi + # Join with commas local IFS=',' echo "${dirs[*]}" } BIND_MOUNTS="$(collect_bind_dirs)" +# setup_slurm.sh recomputes BIND_MOUNTS once the S3 hash dirs are known; the generator-bind notice +# above has already been shown, so suppress it on that second pass. +_REPORTED_GEN_BINDS=1 # Export variables for use in the main script export BIND_MOUNTS @@ -421,4 +592,5 @@ export CONTROL_HASH_VERSION export MEM_PER_JOB export CPUS_PER_JOB export PAR_LOW_COV_THRESHOLD_RUFUS_ARG -export DEV_BIND_ARGS \ No newline at end of file +export DEV_BIND_ARGS +export HASH_SIZE_RUFUS_ARG \ No newline at end of file diff --git a/singularity/setup_slurm.sh b/singularity/setup_slurm.sh index f7b504e..d51098a 100644 --- a/singularity/setup_slurm.sh +++ b/singularity/setup_slurm.sh @@ -69,6 +69,111 @@ fi # Re-compute bind mounts now that hash dirs may have been set by S3 downloads BIND_MOUNTS="$(collect_bind_dirs)" +# Preflight every generator input by running it inside the container. +# +# A generator is an arbitrary shell script RUFUS executes to obtain reads, so not all of its +# dependencies are discoverable by reading it: paths assembled from environment variables, the +# REF_PATH/REF_CACHE lookup samtools uses to decode CRAM without an explicit -T, and any binary it +# shells out to. collect_bind_dirs() binds what it can find statically; this catches the rest. +# +# Without it those failures surface inside a queued job -- in windowed mode, across every task in +# the array -- hours after submission. Running the generator here costs seconds and surfaces the +# generator's own error text. +# +# Uses a detected runtime rather than the hard-coded `singularity` of the generated scripts (see +# the TODO at the top of this file): this is new setup-time code, so it can do the right thing +# without touching the generated-script machinery that deferral is about. If no runtime is on PATH +# the check is skipped with a warning -- setup has never required one, and refusing to generate +# scripts on a submit host without a container runtime would be a regression. +preflight_generators() { + local generators=() + local f + for f in "${SUBJECTS_RUFUS_ARG[@]}" "${CONTROLS_RUFUS_ARG[@]}"; do + [ "$(get_input_type "$f")" == "generator" ] && generators+=("$f") + done + [ ${#generators[@]} -eq 0 ] && return 0 + + # setup_slurm.sh is normally invoked from *inside* the container: + # apptainer exec rufus.sif bash /opt/RUFUS/singularity/setup_slurm.sh ... + # as every README example shows. No container runtime exists inside the image, and none is + # needed -- this is already the environment the generator will run in, samtools included -- so + # run the generator directly. Only enter the container when setup is running on the host. + # + # Note the binds differ between the two cases: in-container we inherit whatever the caller's + # exec bound (at CHPC the singularity/apptainer module sets *_BINDPATH=/scratch,/uufs), whereas + # the generated job scripts use BIND_MOUNTS. A generator that reads from a path bound only in + # the job environment can therefore fail here; the error text below says so. + local -a run_prefix=() + if [ -n "${APPTAINER_CONTAINER:-}${SINGULARITY_CONTAINER:-}" ] || [ -d /.singularity.d ]; then + : # already inside the container -- run_prefix stays empty + else + local runtime sif + runtime="$(command -v apptainer || command -v singularity)" || runtime="" + if [ -z "$runtime" ]; then + echo "WARNING: setup is not running inside a container and neither apptainer nor" >&2 + echo " singularity is on PATH; skipping the generator preflight. Generator" >&2 + echo " errors will not surface until the jobs run." >&2 + return 0 + fi + sif="${CONTAINER_PATH_RUFUS_ARG:-rufus.sif}" + if [ ! -f "$sif" ]; then + echo "WARNING: container $sif not found; skipping the generator preflight." >&2 + return 0 + fi + run_prefix=("$runtime" exec --bind "${BIND_MOUNTS}${DEV_BIND_ARGS}" "$sif") + fi + + local gen out rc checked=0 + for gen in "${generators[@]}"; do + # A generator with a pre-built hash beside it is never executed: runRufus.sh keeps the + # naming and RunJellyForRUFUS.sh skips jellyfish when + # .Jhash exists. That is how pre-built DSA/control hashes are + # supplied, and such a generator is legitimately empty, so there is nothing to preflight. + # Generators are rejected in windowed mode, so .wg is the only postfix reachable here. + if [ -e "${gen}.wg.Jhash" ]; then + echo "Skipping preflight for $(basename "$gen"): pre-built hash ${gen}.wg.Jhash will be used instead." + continue + fi + checked=$((checked + 1)) + + # head closes the pipe once it has enough to judge, so the generator is not run to + # completion; its exit status is therefore not meaningful and the output is what we check. + out="$(timeout 120 "${run_prefix[@]}" bash -c "bash '$gen' 2>&1 | head -20" 2>&1)" + rc=$? + + if [ $rc -eq 124 ]; then + echo "ERROR: generator $gen produced no output within 120s inside the container." >&2 + echo " A generator that blocks this long is usually waiting on a reference or" >&2 + echo " index that is not bound into the container." >&2 + exit 1 + fi + + if echo "$out" | grep -qE '^@(HD|SQ|RG|PG|CO)[[:space:]]'; then + continue + fi + if echo "$out" | awk -F'\t' 'NF>=11 { found=1; exit } END { exit !found }'; then + continue + fi + + echo "ERROR: generator $gen did not produce SAM inside the container." >&2 + echo " RUFUS runs generators with 'bash ' and expects SAM on stdout." >&2 + echo " Output was:" >&2 + if [ -n "$out" ]; then + echo "$out" | sed 's/^/ /' >&2 + else + echo " (no output)" >&2 + fi + echo " An empty generator, or one whose data paths are not visible here, fails this" >&2 + echo " way. If the paths are only bound in the job environment, bind them for setup" >&2 + echo " too (SINGULARITY_BINDPATH/--bind) or add them with -d." >&2 + exit 1 + done + if [ "$checked" -gt 0 ]; then + echo "Generator preflight passed ($checked of ${#generators[@]} generator(s) produced SAM in the container)." + fi +} +preflight_generators + WORKING_DIR=$(pwd) echo -en "##RUFUS_callCommand=" > rufus.cmd @@ -115,13 +220,32 @@ function write_out_rest_of_rufus_args() { echo -en "\$HASH_ARGS " >> rufus.cmd fi - # Use -cr for CRAM inputs, -r otherwise (check first subject) - local ref_flag="-r" - if [[ "${SUBJECTS_RUFUS_ARG[0]}" == *.cram ]]; then - ref_flag="-cr" + # -cr is needed if ANY input is a CRAM, subject or control: runRufus.sh kills the run the moment + # it decodes a .cram with _arg_cramref unset. Keying this off the first subject alone broke every + # mixed set (e.g. -s x.generator -c y.cram, or a BAM subject with a CRAM control). + # + # Both flags are emitted rather than swapping one for the other. runRufus.sh only assigns + # _arg_ref from _arg_cramref inside its per-file CRAM branches, which run after it computes + # _arg_ref_cat="${_arg_ref%.*}". With -cr alone and a non-CRAM subject, _arg_ref_cat would be + # empty at that point and the BWA prefix would fall back to the reference path instead of the + # extension-stripped prefix. Passing -r as well sets _arg_ref up front; both take the same path. + local ref_flags="-r $REFERENCE_RUFUS_ARG" + local _input + for _input in "${SUBJECTS_RUFUS_ARG[@]}" "${CONTROLS_RUFUS_ARG[@]}"; do + if [[ "$_input" == *.cram ]]; then + ref_flags="-r $REFERENCE_RUFUS_ARG -cr $REFERENCE_RUFUS_ARG" + break + fi + done + echo -en "$ref_flags -m $KMER_DEPTH_CUTOFF_RUFUS_ARG -k 25 -t $THREAD_LIMIT_RUFUS_ARG -L -vs " >> $RUFUS_SLURM_SCRIPT + echo -en "$ref_flags -m $KMER_DEPTH_CUTOFF_RUFUS_ARG -k 25 -t $THREAD_LIMIT_RUFUS_ARG -L -vs " >> rufus.cmd + + # Hash size (-hs) for the k-mer count step. Must match the -s of any pre-built control/DSA/exclude + # hash or runRufus.sh's merge aborts; left unset, RUFUS applies its own default (16G wg / 1G window). + if [ -n "$HASH_SIZE_RUFUS_ARG" ]; then + echo -en "-hs $HASH_SIZE_RUFUS_ARG " >> $RUFUS_SLURM_SCRIPT + echo -en "-hs $HASH_SIZE_RUFUS_ARG " >> rufus.cmd fi - echo -en "$ref_flag $REFERENCE_RUFUS_ARG -m $KMER_DEPTH_CUTOFF_RUFUS_ARG -k 25 -t $THREAD_LIMIT_RUFUS_ARG -L -vs " >> $RUFUS_SLURM_SCRIPT - echo -en "$ref_flag $REFERENCE_RUFUS_ARG -m $KMER_DEPTH_CUTOFF_RUFUS_ARG -k 25 -t $THREAD_LIMIT_RUFUS_ARG -L -vs " >> rufus.cmd if [ "${PAR_LOW_COV_THRESHOLD_RUFUS_ARG}" != "7" ]; then echo -en "-plct $PAR_LOW_COV_THRESHOLD_RUFUS_ARG " >> $RUFUS_SLURM_SCRIPT diff --git a/tests/functional/cases/f7_multi_subject.sh b/tests/functional/cases/f7_multi_subject.sh new file mode 100644 index 0000000..41c9995 --- /dev/null +++ b/tests/functional/cases/f7_multi_subject.sh @@ -0,0 +1,223 @@ +#!/bin/bash +# F7 — MULTI-SUBJECT INPUT (-s repeated). +# +# A sample sequenced across several files (different centers, flowcells, or split lanes) is passed +# as repeated -s arguments; setup_slurm.sh emits exactly that form from its comma-separated -s list. +# runRufus.sh concatenates one read-emitting command per subject into a single generator, and the +# generator body is run as ONE stream (`bash generator | samtools ...`). So every subject file must +# contribute its reads while exactly ONE command emits a SAM header -- samtools aborts on a second +# @HD mid-stream, and because `collate` buffers before emitting, it discards the WHOLE stream rather +# than truncating it. The visible symptom is a clean-looking run that calls nothing at all, in every +# region ("no_reads_passed_filter"), with the real error buried in stderr. +# +# The test splits the somatic tumor fixture into two files by read name -- a true partition, so the +# union is exactly the original single file -- and asserts a multi-subject run reproduces the +# single-file result: +# PARTITION the two parts are non-empty and their read counts sum to the original +# HEADER the multi-subject generator has 2 read commands but emits exactly 1 @HD +# RECALL every planted somatic is recovered from the split input +# CONCORDANCE split call set == single-file call set (bam), and the cram arm matches too +# +# The HEADER assertion is what pinpoints a regression: without it, a reintroduced duplicate header +# shows up only as "no calls", which is indistinguishable from a dozen unrelated failures. F3 already +# establishes bam/cram call-set equality, so the cram arm here isolates the -T/-cr generator line +# rather than re-testing format concordance. +# +# Run directly: bash f7_multi_subject.sh +# Or submit: sbatch f7_multi_subject.sh +# EXTRA_BIND=host:container[,...] overlays uncommitted code (dev loop) before the image is rebuilt. +#SBATCH --account=marth-rw +#SBATCH --partition=marth-rw +#SBATCH --job-name=rufus_f7_multi_subject +#SBATCH --nodes=1 +#SBATCH --cpus-per-task=8 +#SBATCH --mem=32G +#SBATCH --time=02:00:00 +#SBATCH --output=%x_%j.out +#SBATCH --error=%x_%j.err +set -euo pipefail + +# sbatch copies this script to the spool dir, so BASH_SOURCE-derived paths break. See F4. +HERE="${FUNCTIONAL_CASES_DIR:-${SLURM_SUBMIT_DIR:-$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)}}" +FIX="$(cd "$HERE/../fixtures" 2>/dev/null && pwd)" || { echo "ERROR: cannot locate fixtures from HERE=$HERE"; exit 1; } +DATA=/uufs/chpc.utah.edu/common/HIPAA/u0746015/marth_software/RUFUS/resources/reg_test_files +SIF=${SIF:-$DATA/rufus_dev.sif} +OUTROOT=${OUTROOT:-$DATA/runs/f7_multi_subject} +THREADS=${SLURM_CPUS_PER_TASK:-8} +WINDOW=10 +EXTRA_BIND=${EXTRA_BIND:-} + +REF=$FIX/ref/tiny.fa +SOMATIC=$FIX/designed/somatic.vcf # expected somatic loci +SPLIT=$OUTROOT/_split # derived inputs, rebuilt every run + +module load apptainer 2>/dev/null || true +command -v apptainer >/dev/null || { echo "ERROR: apptainer unavailable"; exit 1; } +[ -f "$SIF" ] || { echo "ERROR: SIF not found: $SIF"; exit 1; } +for f in "$REF" "$FIX/somatic/tumor.bam" "$FIX/somatic/normal.bam" "$FIX/somatic/normal.cram" "$SOMATIC"; do + [ -f "$f" ] || { echo "ERROR: fixture missing: $f (run make_fixtures.sh)"; exit 1; } +done + +rm -rf "$OUTROOT"; mkdir -p "$SPLIT" +BINDS="$FIX,$OUTROOT,$DATA"; [ -n "$EXTRA_BIND" ] && BINDS="$BINDS,$EXTRA_BIND" + +fail() { echo "RESULT: FAIL — $*"; exit 1; } +insif() { apptainer exec --bind "$BINDS" "$SIF" "$@"; } + +echo "=== F7 multi-subject input (-s repeated) | $(date) ===" +[ -n "$EXTRA_BIND" ] && echo " OVERLAY : $EXTRA_BIND (testing an uncommitted local fix)" +echo + +# --------------------------------------------------------------------------------------------- +# Build the split. Reads are assigned by QNAME so both mates land in the same part: the parts are +# a true partition of the original, which is what makes "split == single file" a fair comparison. +# Coordinate order is preserved within each part, but the concatenated stream is no longer globally +# sorted -- RUFUS name-sorts the filtered mates before assembly, so that is expected to wash out. +# --------------------------------------------------------------------------------------------- +echo "-- building split fixtures from tumor.bam --" +# All of the setup runs in ONE container invocation. Container startup dominates this step -- it has +# been seen at ~5 min per exec on a busy node -- so issuing eleven of them made a 9-minute test take +# well over an hour. The work itself is seconds. Written as a file rather than an inline `bash -c` +# so the awk program needs no second level of shell quoting. +cat > "$SPLIT/build_split.sh" <<'SPLITEOF' +#!/bin/bash +# Partition SRC into two BAMs by QNAME (both mates land in the same part), add CRAM copies of each, +# and record read counts so the caller can verify the parts really are a partition. +set -euo pipefail +SRC=$1; REF=$2; OUT=$3 +for p in 1 2; do + samtools view -h "$SRC" \ + | awk -v part="$p" 'BEGIN{OFS="\t"} /^@/{print; next} { if (!($1 in m)) m[$1]=(++n % 2); if (m[$1]==part%2) print }' \ + | samtools view -b -o "$OUT/tumor.part$p.bam" - + samtools index "$OUT/tumor.part$p.bam" + samtools view -C -T "$REF" -o "$OUT/tumor.part$p.cram" "$OUT/tumor.part$p.bam" + samtools index "$OUT/tumor.part$p.cram" +done +{ samtools view -c "$SRC" + samtools view -c "$OUT/tumor.part1.bam" + samtools view -c "$OUT/tumor.part2.bam"; } > "$OUT/counts.txt" +SPLITEOF +insif bash "$SPLIT/build_split.sh" "$FIX/somatic/tumor.bam" "$REF" "$SPLIT" + +# PARTITION: non-empty parts that sum to the original. A split that silently produced one empty +# file would make the multi-subject run trivially equal to a single-subject run and prove nothing. +N_ALL=$(awk 'NR==1{print $1}' "$SPLIT/counts.txt") +N_P1=$(awk 'NR==2{print $1}' "$SPLIT/counts.txt") +N_P2=$(awk 'NR==3{print $1}' "$SPLIT/counts.txt") +[ -n "$N_ALL" ] && [ -n "$N_P1" ] && [ -n "$N_P2" ] || fail "split build did not produce read counts — see $SPLIT" +echo " reads: original=$N_ALL part1=$N_P1 part2=$N_P2 (sum=$((N_P1 + N_P2)))" +[ "$N_P1" -gt 0 ] && [ "$N_P2" -gt 0 ] || fail "split produced an empty part (part1=$N_P1 part2=$N_P2)" +[ "$((N_P1 + N_P2))" -eq "$N_ALL" ] || fail "split is not a partition: $N_P1 + $N_P2 != $N_ALL" + +# --------------------------------------------------------------------------------------------- +# run_case