Merging multiple files if no reads are found in one of the files fails

[rkostadi]

When using samtools merge with multiple bam files (-b bamfilelist.txt) and
with selecting a specific region with e.g. -R chr1:1-1000,
I noticed that if one of the files (fileX.bam) is missing any reads within the region,
samtools merge will fail with the message:
[bam_merge_core2] failed to read first record from fileX.bam

I think I corrected this behavior by
replacing lines 1338-1339 in bam_sort.c :

print_error(cmd, "failed to read first record from \"%s\"", fn[i]);
goto fail;

with:

print_error(cmd, "failed to read first record from \"%s\", skipping file during merge", fn[i]);
h->pos = HEAP_EMPTY;
bam_destroy1(h->b);
h->b = NULL;

This appears to fix the problem, so that the new behavior is that
fileX.bam is skipped, and the rest of the merge proceeds with the rest of the files.

Do you think that this fix is correct and should be the default behavior of samtools merge?

[daviesrob]

In theory the 'no reads in region' case should be covered by line 1333:

else if (res == -1 && (!iter[i] || iter[i]->finished)) {

while 1338-1339 is supposed to handle errors (e.g. corrupt data, I/O problems) when reading the file.

We need to make a test case to see why line 1333 isn't doing what it's supposed to here.

[jenniferliddle]

I've tried to reproduce this problem, without success. Could you:

• produce a set of (preferably small) files that show the problem?
• tell me exactly what command line options you are using
• tell me which version of samtools you're using (samtools --version)
  @rkostadi

[rkostadi]

Hi,

This is reproducible with samtools version:

samtools 1.3.1-50-gd02c177
Using htslib 1.3.2-221-g9e2fbac

and version:

samtools 1.3.1
Using htslib 1.3.2

test1.bam contains 2 reads that occur at
chr1 725432
chr1 725464
test2.bam contains 2 reads that occur at
chr1 10024
chr1 10065

This fails to produce the merged file "test4.bam":

[rumen@vishera bin]$ ./samtools --version
samtools 1.3.1-50-gd02c177
Using htslib 1.3.2-221-g9e2fbac
Copyright (C) 2017 Genome Research Ltd.
[rumen@vishera bin]$ ./samtools merge -R chr1:725000-726000 test4.bam test1.bam test2.bam
samtools merge: failed to read first record from "test2.bam"

This works fine, since the region captures all of the reads from test1.bam and test2.bam:

[rumen@vishera bin]$ ./samtools merge -R chr1:1-726000 test4.bam test1.bam test2.bam
[rumen@vishera bin]$

Attached are the test files.
test.zip

[jenniferliddle]

Thanks for this. I'll take a look....

[jkbonfield]

Thanks for the test data. It was very useful in reproducing the problem (now fixed).

More debateable is why this error doesn't occur with other iterator uses. It appears that merge is the only function that validates iter->finished as well as return -1. I am unsure if this is an incorrect check or whether everything else has poor error checking - both are equally likely and there is no documentation on the iterator API sadly.

[rkostadi]

Thanks @jkbonfield, this works much better than my hack.