StrainR2 stalls during BBMap mapping on InSilicoSeq-generated paired-end metagenomes (no error, job never completes)

[mattgal11]

Hello, I am encountering a reproducible issue where StrainR2 stalls during the StrainR (BBMap) mapping step when using paired-end metagenomes generated with InSilicoSeq. The job does not error, but appears to go idle and never completes.

Summary

• StrainR2 stalls during the BBMap mapping stage
• Happens consistently for InSilicoSeq-generated paired-end metagenomes across multiple sizes (e.g., 20M - 100M paired reads)
• StrainR2 has previously worked for me on both smaller and larger real metagenomes

Environment

• StrainR2 v2.3.0
• BBMap version: 39.52 (via conda)
• Java: OpenJDK (conda)
• OS: Linux HPC (SLURM)
• Threads: typically 8 (tested 8-32)
• Memory: tested 8–64 GB

Input data

• Paired-end Illumina-style metagenome FASTQs generated with InSilicoSeq
• PE150 reads with file sizes varying from 1.2 Gb - 6.5 Gb in .fastq.gz format

What happens during StrainR run

• (separate) PreProcessR (completes successfully)
• fastp trimming (completes successfully)
• BBMap starts mapping and spawns threads
• Log output consistently ends at "Started 8 mapping threads"
• The output SAM file is created but stalls at exactly 64 KB each time and no further progress is made
• The job runs until the SLURM walltime limit is reached (typically set at 4/8 hours), despite previous runs completing in minutes
• An e.g., error log from a stalled 100M paired-end run is included in full below (though it is the same for smaller runs, e.g., 20M):

`Read1 before filtering:
total reads: 100000000
total bases: 12600000000
Q20 bases: 12085270862(95.9148%)
Q30 bases: 11588108014(91.9691%)
Q40 bases: 0(0%)

Read2 before filtering:
total reads: 100000000
total bases: 12600000000
Q20 bases: 11881205624(94.2953%)
Q30 bases: 11268220785(89.4303%)
Q40 bases: 0(0%)

Read1 after filtering:
total reads: 99251208
total bases: 12505502489
Q20 bases: 11998033847(95.942%)
Q30 bases: 11505770740(92.0057%)
Q40 bases: 0(0%)

Read2 after filtering:
total reads: 99251208
total bases: 12505501372
Q20 bases: 11828835065(94.5891%)
Q30 bases: 11233653434(89.8297%)
Q40 bases: 0(0%)

Filtering result:
reads passed filter: 198502416
reads failed due to low quality: 1493410
reads failed due to too many N: 314
reads failed due to too short: 3860
reads with adapter trimmed: 11610
bases trimmed due to adapters: 630196

Duplication rate: 1.29851%

Insert size peak (evaluated by paired-end reads): 221

JSON report: strainR2_results/HiSeq_100M_K12_90_Sakai_10_rep3
HTML report: strainR2_results/HiSeq_100M_K12_90_Sakai_10_rep3

fastp -i final_metagenomes/HiSeq_100M_K12_90_Sakai_10_rep3_R1.fq.gz -o strainR2_results/HiSeq_100M_K12_90_Sakai_10_rep3/tmp/forward.fastq.gz -I final_metagenomes/HiSeq_100M_K12_90_Sakai_10_rep3_R2.fq.gz -O strainR2_results/HiSeq_100M_K12_90_Sakai_10_rep3/tmp/reverse.fastq.gz --trim_poly_g --json strainR2_results/HiSeq_100M_K12_90_Sakai_10_rep3 --html strainR2_results/HiSeq_100M_K12_90_Sakai_10_rep3 --length_required 50 --n_base_limit 0 --thread 8
fastp v1.0.1, time used: 342 seconds
java -ea --add-modules jdk.incubator.vector -Xmx8g -Xms8g -cp /well/bag/users/wfr051/.conda/envs/strainr2/opt/bbmap-39.52-0/current/ align2.BBMap build=1 overwrite=true fastareadlen=500 in=strainR2_results/HiSeq_100M_K12_90_Sakai_10_rep3/tmp/forward.fastq.gz in2=strainR2_results/HiSeq_100M_K12_90_Sakai_10_rep3/tmp/reverse.fastq.gz ref=StrainR2_ISS_metagenome_db/BBindex/BBIndex.fasta outm=strainR2_results/HiSeq_100M_K12_90_Sakai_10_rep3/sample.sam rpkm=strainR2_results/HiSeq_100M_K12_90_Sakai_10_rep3/sample.rpkm threads=8 deterministic=t averagepairdist=200 -Xmx8g perfectmode=t local=f ambiguous=toss pairedonly=t nodisk=t
WARNING: Using incubator modules: jdk.incubator.vector
Executing align2.BBMap [build=1, overwrite=true, fastareadlen=500, in=strainR2_results/HiSeq_100M_K12_90_Sakai_10_rep3/tmp/forward.fastq.gz, in2=strainR2_results/HiSeq_100M_K12_90_Sakai_10_rep3/tmp/reverse.fastq.gz, ref=StrainR2_ISS_metagenome_db/BBindex/BBIndex.fasta, outm=strainR2_results/HiSeq_100M_K12_90_Sakai_10_rep3/sample.sam, rpkm=strainR2_results/HiSeq_100M_K12_90_Sakai_10_rep3/sample.rpkm, threads=8, deterministic=t, averagepairdist=200, -Xmx8g, perfectmode=t, local=f, ambiguous=toss, pairedonly=t, nodisk=t]
Version 39.52

Ambiguously mapped reads will be considered unmapped.
Executing dna.FastaToChromArrays2 [StrainR2_ISS_metagenome_db/BBindex/BBIndex.fasta, 1, writeinthread=false, genscaffoldinfo=true, retain, waitforwriting=false, gz=true, maxlen=536670912, writechroms=false, minscaf=1, midpad=300, startpad=8000, stoppad=8000, nodisk=true]

Set genScaffoldInfo=true
Set genome to 1

Loaded Reference: 0.004 seconds.
Loading index for chunk 1-1, build 1
Indexing threads started for block 0-1
Indexing threads finished for block 0-1
Generated Index: 9.066 seconds.
Analyzed Index: 3.030 seconds.
Started output stream: 0.479 seconds.
Cleared Memory: 0.151 seconds.
Processing reads in paired-ended mode.
Started read stream.
Started 8 mapping threads.`

Additional context

I have also tried:

• Increasing and decreasing memory allocations
• Changing thread counts
• Rebuilding PreProcessR reference database with different parameters (different numbers of strains, simplified FASTA headers)
• Verifying simulated FASTQ read lengths and pairing (though these have been used successfully elsewhere before)
• Using compressed and uncompressed FASTQs and renaming file extensions (e.g., .fq.gz vs .fastq.gz)

None of these resolved the issue. As mentioned StrainR2 has worked for me previously on varying file sizes (~33M human-depleted reads - ~1.8Gb paired-end .fastq.gz) in a matter of minutes, so I am not sure what is causing the stalling that I am encountering. Thanks for your time/help!

[BisanzLab]

Thanks for this detailed summary! Would you be willing to share one of your smallest in silico samples that creates this issue so we can investigate further?

[mattgal11]

HiSeq_20M_K12_50_Sakai_50_rep1_R1_sub250k.fastq.gz

HiSeq_20M_K12_50_Sakai_50_rep1_R2_sub250k.fastq.gz

Of course - I have attached a paired-end InSilicoSeq subsample (250k read pairs) that still reproduces the stalling behaviour during the BBMap step on my end. Let me know if one of the full, larger files would be more helpful. For these runs the reference database was built successfully from 24 E. coli genomes using PreProcessR defaults.