Docker_v1.2.0

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:

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,19 +826,29 @@ 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 $$
     fi
 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 <control><region_postfix> 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 <control><region_postfix>.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=()

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