Sambamba sort segfault + why sambamba?

[nikostr]

Hi! Thanks for this pipeline! I'm currently trying to get it to run, and mostly it seems to work, but I keep getting segfaults for sambamba sort as a part of F05_map_reads_for_blob. I have tried troubleshooting this, increasing memory and attempting to set the stack size using ulimit, though I'm unsure whether ulimit is propagating to the process on the cluster I'm using, but no luck. Any recommendations on handling this?

As a part of this I got to wondering why sambamba sort is being used. According to sambamba docs (https://github.com/biod/sambamba#introduction) samtools is faster for sorting, and sambamba seems to have seen little maintenance since august of 2023 (no closed issues, no commits to the repo since then), while samtools seems to be more actively maintained. Is there some reason in particular for using sambamba here? Would it be worth considering moving to samtools sort?

[diegomics]

Hey, thanks for pointing this out. Ok, let's swtich to samtools sort. I'm working on the next update already on the PR, but I will add this one too. Give me a couple of days ;)

[diegomics]

Please check if it's fixed with the new version (re-open if not!)

[nikostr]

I've now had the opportunity to run this, and the issue persists in that sambamba view crashes as well. Samtools is capable of replacing view, merge and flagstat pretty much out of the box if I understand things correctly. Maybe getting rid of the dependency entirely is the way to go here? When I fixed this myself relative to v0.4 I made the following changes:

From a91a9c222367723e6bbd56f61dab3c806d1a6adb Mon Sep 17 00:00:00 2001
From: nikostr <nikostr@users.noreply.github.com>
Date: Tue, 23 Jun 2026 13:15:56 +0200
Subject: [PATCH] Attempt at replacing sambamba sort with samtools sort

---
 workflow/rules/F_cont_analysis.smk | 33 ++++++++++++++----------------
 1 file changed, 15 insertions(+), 18 deletions(-)

diff --git a/workflow/rules/F_cont_analysis.smk b/workflow/rules/F_cont_analysis.smk
index db9df37..efeecc5 100644
--- a/workflow/rules/F_cont_analysis.smk
+++ b/workflow/rules/F_cont_analysis.smk
@@ -964,12 +964,11 @@ rule F05_map_reads_for_blob:
                 {input.asm} \
                 ${{ALL_READS[@]}} \
             | samtools view -bu -@ $VIEW_THREADS \
-            | sambamba sort \
-                -t $SORT_THREADS \
-                -m $(( {resources.mem_mb} * 60 / 100 ))M \
-                --tmpdir="$TEMP_DIR" \
-                -o {output.bam} \
-                /dev/stdin
+            | samtools sort \
+                -@ $SORT_THREADS \
+                -m $(( {resources.mem_mb} / SORT_THREADS * 60 / 100 ))M \
+                -T "$TEMP_DIR" \
+                -o {output.bam}
 
         elif [ "$READ_TYPE" = "illumina" ] || [ "$READ_TYPE" = "10x" ] || [ "$READ_TYPE" = "hic" ]; then
             # Short/paired reads
@@ -985,12 +984,11 @@ rule F05_map_reads_for_blob:
                     {input.asm} \
                     "${{ALL_READS[0]}}" "${{ALL_READS[1]}}" \
                 | samtools view -bu -@ $VIEW_THREADS \
-                | sambamba sort \
-                    -t $SORT_THREADS \
-                    -m $(( {resources.mem_mb} * 60 / 100 ))M \
-                    --tmpdir="$TEMP_DIR" \
-                    -o {output.bam} \
-                    /dev/stdin
+                | samtools sort \
+                    -@ $SORT_THREADS \
+                    -m $(( {resources.mem_mb} / SORT_THREADS * 60 / 100 ))M \
+                    -T "$TEMP_DIR" \
+                    -o {output.bam}
 
             elif [ "$((NUM_READS % 2))" -eq 0 ] && [ "$NUM_READS" -gt 2 ]; then
                 # Multiple pairs - map each, merge
@@ -1009,12 +1007,11 @@ rule F05_map_reads_for_blob:
                         {input.asm} \
                         "${{ALL_READS[$i]}}" "${{ALL_READS[$((i+1))]}}" \
                     | samtools view -bu -@ $VIEW_THREADS \
-                    | sambamba sort \
-                        -t $SORT_THREADS \
-                        -m $(( {resources.mem_mb} * 60 / 100 ))M \
-                        --tmpdir="$TEMP_DIR" \
-                        -o "$PARTIAL_BAM" \
-                        /dev/stdin
+                    | samtools sort \
+                        -@ $SORT_THREADS \
+                        -m $(( {resources.mem_mb} / SORT_THREADS * 60 / 100 ))M \
+                        -T "$TEMP_DIR" \
+                        -o "$PARTIAL_BAM"
 
                     PARTIAL_BAMS+=("$PARTIAL_BAM")
                 done
-- 
2.39.5

and

From 3dd470fe383f3a3c3ec8941ac94075a6570674d3 Mon Sep 17 00:00:00 2001
From: nikostr <nikostr@users.noreply.github.com>
Date: Wed, 24 Jun 2026 10:57:41 +0200
Subject: [PATCH] Replace sambamba with samtools entirely

---
 workflow/rules/F_cont_analysis.smk | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/workflow/rules/F_cont_analysis.smk b/workflow/rules/F_cont_analysis.smk
index efeecc5..a8f8480 100644
--- a/workflow/rules/F_cont_analysis.smk
+++ b/workflow/rules/F_cont_analysis.smk
@@ -943,7 +943,7 @@ rule F05_map_reads_for_blob:
         [ "$MAP_THREADS" -lt 2 ] && MAP_THREADS=2
         [ "$SORT_THREADS" -lt 2 ] && SORT_THREADS=2
 
-        echo "[GEP2] Thread split: minimap2=$MAP_THREADS, samtools=$VIEW_THREADS, sambamba=$SORT_THREADS"
+        echo "[GEP2] Thread split: minimap2=$MAP_THREADS, samtools view=$VIEW_THREADS, samtools sort=$SORT_THREADS"
 
         # Determine read inputs based on read type
         READ_TYPE="{params.read_type}"
@@ -1017,8 +1017,8 @@ rule F05_map_reads_for_blob:
                 done
 
                 echo "[GEP2] Merging $PAIR_IDX partial BAMs..."
-                sambamba merge \
-                    -t {threads} \
+                samtools merge \
+                    -@ {threads} \
                     {output.bam} \
                     "${{PARTIAL_BAMS[@]}}"
 
@@ -1041,13 +1041,13 @@ rule F05_map_reads_for_blob:
         echo "[GEP2] BAM size: $((BAM_SIZE / 1024 / 1024)) MB"
 
         # Quick validation: check BAM header
-        if ! sambamba view -H {output.bam} > /dev/null 2>&1; then
+        if ! samtools view -H {output.bam} > /dev/null 2>&1; then
             echo "[GEP2] ERROR: Output BAM file is not valid" >&2
             exit 1
         fi
 
         # Count mapped reads
-        MAPPED=$(sambamba flagstat -t {threads} {output.bam} 2>/dev/null | head -1 | awk '{{print $1}}')
+        MAPPED=$(samtools flagstat -@ {threads} {output.bam} 2>/dev/null | head -1 | awk '{{print $1}}')
         echo "[GEP2] Total alignments: $MAPPED"
 
         echo "[GEP2] Read mapping complete: {output.bam}"
-- 
2.39.5

Note that the commit message of the last patch is slightly misleading as it actually does not completely get rid of sambamba.

[diegomics]

Hi, I re-opened this. Removing sambamba completely is at the very top of the to-do list. The 0.7 update should land in a couple of weeks including the full swtich to samtools.

[diegomics]

Btw, I couldn't replicate the sambamba segfault with my data here, nevertheless moving to samtools is the right move. @mxmmns we can wrap the new version with the super report and this alone to speed up things

[nikostr]

Yeah, I'm not sure how to replicate it either which makes this whole thing a pain. It's very consistently crashing on my cluster, and sambamba seems to be known to segfault. As said, it might be a ulimit thing, but checking this is tricky for me as it's running through SLURM and snakemake and using containers - I'm even unsure of how to check that a set ulimit propagates under these conditions, haha. Very thankful that using samtools resolves the issue!