Close CRAM reference

[rwdavies]

Hey, have been using seqLib successfully for a while, but I just ran into an issue. I have a function, and within this function, I create a BamReader object, get reads out, etc. This works well, except I've recently noticed that when I use CRAM files and set the reference using BamReader::SetCramReference, the reference is not closed when calling BamReader::Close(), or when the BamReader object goes out of scope when the function ends. So when I call my function on several thousand CRAM files, the code crashes because I don't have enough file handles. Can you please implement a method to close the cram reference?

Error message wise, when running I get this

connect: Network is unreachable

one for each file (this is run on a cluster without internet access), then once I run out of file handles, I get

connect: Network is unreachable
test-data/mouse_data/mm10_2016_10_02.fa: Too many open files
Failed to populate reference for id 10
Unable to fetch reference #10 10001343..10999757
Failure to decode slice

[walaj]

Thanks for finding this bug. Indeed, I never call cram_close to close the reference.

I added this now to BamReader.h. It seems to work on both CRAM or BAM files without any issues. Can you let me know if it solves the problem for your case (with several thousand CRAM files), which I don't have for testing?

[rwdavies]

Hmm, weird, when I compile with the change, it stops working entirely for BAM and CRAM. It passes your tests?

This is usage like

    SeqLib::BamReader reader;
    reader.Open(file_name);
//stuff
reader.Close();

and it now crashes on the Close()

[rwdavies]

So I didn't get any further on the above. Same problem, when I actively call Close(), everything breaks, for BAMs and CRAMs, but when I don't call Close(), the file limit problem persists and things blow up when I use lots of CRAMs

I'm having trouble making a minimal C++ only example using the example code you've listed on the README though? So "make seqtools" works for me normally, but when I try to make a new Frankenstein "seqtools.cpp" with the code from one of the examples, I get the following
log.txt
seqtools.txt

[rwdavies]

OK, now I've got what should be reproducible. With the commit I get a segfault when trying to close the reader

bash test.txt &> log.txt

log.txt
test.txt
seqtools.txt

[walaj]

Before you spend more time on this, I realized that I had made an error with the commit to fix this issue (and erased the commit).

It turns out that the original code does close the cram reference. I put a print statement in the htslib function that closes the cram reference, and confirmed it will indeed close with either calling the Close command or if the object goes out of scope.

So that leaves me puzzled about your error. Could it be that your readers are not going out of scope or being closed soon enough?

[rwdavies]

Ah, you deleted the commit. I thought I was crazy for awhile there.

Anyway, I'm still getting the bug. Here's what should be reproducible, using one small main and one function that loads some reads out of a cram. This occurs even though I am using Close. The example runs fine with internet access, and only occurs without internet access when the file limit runs out. I'm not 100% sure if it's actually something different between the two machines, the head nodes (with internet) and the compute nodes (without internet) are otherwise fairly similar, although some shared libraries might not be available.

bash test.setup.txt &> log.setup.txt
bash test.run.txt &> log.withinternet.txt

and then on a cluster node without internet

bash test.run.txt &> log.withoutinternet.txt

log.withoutinternet.txt
log.withinternet.txt
log.setup.txt
seqtools.txt
test.run.txt
test.setup.txt

[walaj]

Thanks for sending along an example. Right now, it may take me some time to get to this but I am curious to get this resolved. Do you have any issues with similar pipelines that query a purely HTSlib based tool, like SAMtools? I am wondering if the issue is with SeqLib handling or with HTSlib itself.