updated to make consensus work with sam.gz along with threads and ite…#2269
Conversation
| thread; filters requiring header also needs this irrespective of file | ||
| type; it is easy to copy from tdata[0] rather than reading and parsing; | ||
| will be released along with file close */ | ||
| tdata[i].fp->bam_header = sam_hdr_dup(tdata[0].fp->bam_header); |
There was a problem hiding this comment.
htsFile ought to be an opaque type, but is not for implementation reasons (mostly that sam_itr_next()/etc are macros). It seems a shame to introduce a direct reference to the bam_header field, which has not been (ab)used from outside htslib previously.
There was a problem hiding this comment.
Yes, it would be good to have an interface around this in HTSlib. It could probably also take a reference to the existing header instead of duplicating it, which would be more memory efficient.
There was a problem hiding this comment.
This code could call ignore = sam_hdr_read(tdata[i].fp) (and would probably need to destroy the header returned too) to set up bam_header without new interfaces.
But that sort of API function would be a useful addition too.
There was a problem hiding this comment.
Not loading the header for every file descriptor was deliberate as it's costly in memory. However I obviously didn't think about some filters that need the header.
I note the sam_hdr_incr_ref htslib interface exists but appears to never be used in htslib nor samtools/bcftools, which is a pity.
It originated from io_lib and with a scram_set_header function which is used in the scramble tool to share a header between input and output file descriptors.
Htslib lacks the set-header function, but this is the sort of scenario it would be beneficial for.
There was a problem hiding this comment.
Yes, sam_hdr_read will work but it adds the overhead of parsing header which I was trying to avoid.
But accessing bam_header directly is worse than that!
I will add get/set header and use.
There was a problem hiding this comment.
Parsing headers 4 or 16 times (or whatever small integer number of threads is specified) doesn't seem like a major overhead. 🤷
This field is not always set, so it should perhaps be a more abstracted API than simple getter/setters.
|
Added get/set methods to htslib and used it. Depends on htslib PR 1959. |
Consensus crashes when sam.gz is used as input along with threads and iterators.
This happens due to the data validation attempt in sam_parse1 where fp->bam_header, which is null, is used for header.
To avoid this, each fp in thread data is set with a copy of already read header.
This is made for bam and cram as well, for filters to work fine where an expression needs to access the header (like rname, mrname, library etc.)