faidx/fqidx indexes along with output creation#2125
Conversation
30b3ab5 to
8233445
Compare
There was a problem hiding this comment.
I still have more testing to do, but it looks like it works.
It was more complex than I envisaged though and I wonder really what the benefit of building the index on the fly is given for .gz we then have to run bgzf_index_dump anyway to sequentially reread the data and generate a .gzi. I benchmarked it and it is still marginally quicker to build the index this way, so given we've now got the code I don't think it's worth removing it, but it would probably have been simpler to just add a call to fai_build at the end and let it do all the work.
Edit: ignore the bgzf_index_dump comment - it doesn't need to reread the data as it's already doing it on the fly due to the earlier bgzf_index_build_init call. So yes, this is a better implementation than the naive one, albeit more complex.
| .TP 8 | ||
| .BI "-o, --output " FILE | ||
| Write FASTQ to file rather than to stdout. | ||
| Write FASTQ to file rather than to stdout. With .gz, .bgz, .bgzf file extensions |
There was a problem hiding this comment.
I thought about adding an "and" in the extension list, but maybe flipping it around is more readable. Eg:
"If FILE ends with .gz, .bgz or .bgzf then it will be BGZF compressed."
There was a problem hiding this comment.
applied to both faidx and fqidx.
| # define ABS(x) ((x)>=0?(x):-(x)) | ||
| #endif | ||
|
|
||
| /// holds the indexing info for each read name and offsets |
There was a problem hiding this comment.
Having both an idx and an idx_t type name is a little confusing. Maybe rename the first to idx_row or idx_entry to indicate single idx entry vs whole index?
| return out->isbgzip ? bgzf_write(out->bgzf_fp, buffer, length) : | ||
| fwrite(buffer, 1, length, out->fp); | ||
| if ((in) && ((in->n + 1) >= in->m)) { //+1 to ensure space to save seq offset for next one | ||
| int newlen = in->m < 1 ? 16 : in->m << 1; //doubles when reallocation needed |
There was a problem hiding this comment.
I was debating whether int is the correct type here. Should it be unsigned (along with n and m)? Or size_t?
The answer is "probably", but it's all a moot point given BAM and CRAM headers would fall over well before getting remotely close to that point due to the size of the SQ lines. The pragmatist wins out, unless Rob strongly disagrees.
| returns updated input, if reallocation is required | ||
| */ | ||
| static inline size_t wrappedwrite(output *out, const char *buffer, size_t length) | ||
| static inline idx* allocidx(idx* in) |
There was a problem hiding this comment.
This function is a little peculiar. The comment above it implies that the idx pointer may be changing and when calling we need something like:
if ((in = allocidx(in)) == NULL) {err...}
However the input pointer never changes. The contents do, but the pointer does not. The natural return value therefore would be 0 vs <0 for error as returning a pointer is misleading as to what's going on.
However the code calling this is complex with lots of things like:
out->idxdata->indx[out->idxdata->n].line_length = ...
I feel therefore a more natural type for this function would be to return idx_t * (or idx_entry *, see above) with a return &in->indx[in->n++]. Then we can assign it when calling and use that pointer when filling it out. Ie:
idx_t *e; // index entry
if (out->gopt->write_index && !(e = allocidx(out->idxdata))) {
// ... err
}
// ...
e->name = tmp;
e->seq_offset += len; //last offset + length of name line
e->seq_length = seq_len;
This makes the code considerably easier to read than out->idxdata->indx[out->idxdata->n].name etc.
There was a problem hiding this comment.
I'm also thinking that allocating for n+1 is confusing. It took a lot of understanding what's going on here.
However it appears while filling out entry n it stores the cumulated offset of where we've got so far in entry n so we can then append to that when we set it via ...->seq_offset += len. It feels like it would be much more explicit to simply have an accumulated offset in out->idxdata itself. Eg ->seq_offset = out->idxdata->offset + len and then assign to that instead of allocating +1 in the index table and storing the data in there.
| enum fai_format_options format) { | ||
| hts_pos_t seq_len, wrap_len = length, len = 0; | ||
| char *seq = NULL, *qual = NULL; | ||
| char *seq = NULL, *qual = NULL, *tmp = NULL; |
There was a problem hiding this comment.
There's no need to initialise tmp to NULL here, and infact it could simply be moved down to the one block where it's used.
| } else { //fasta, next sequence's temp offset | ||
| out->idxdata->indx[out->idxdata->n + 1].seq_offset = out->idxdata->indx[out->idxdata->n].seq_offset + | ||
| seq_len + seq_len / wrap_len + ((seq_len % wrap_len) ? 1 : 0); | ||
| ++out->idxdata->n; |
There was a problem hiding this comment.
Why are there two places incrementing out->idxdata->n? It's confusing as fasta does it in one place and fastq in another, leading to code duplication. Ultimately though we only have one index entry regardless of format.
Again, if we had a single function to allocate and return the next index entry (allocidx()) then this could increment n itself.
|
After our discussions yesterday on simplifying the code structure, I made my own commit on top of this to refactor it. You can see the changes here: https://github.com/jkbonfield/samtools/tree/writeidx There were many changes, some stylistic and some refactoring, but no actual bug fixes. The code you have here does work, albeit perhaps overly complex. Hence feel free to push back on the changes if you wish or we can discuss style things in our weekly meeting with the wider group. I also double checked Windows. You have added the assumption that the number of characters per line is always strictly 1 more than the number of bases per line when you're building the index. I was initially concerned this may break on windows, but thanks to here and here the files we create are always in binary mode so the windows \r\n line endings cannot happen on files we create ourselves. Hence dumping the index is valid as-is, but I added a comment to that effect given it was a concern that needed investigation. |
|
Thanks James, I have used the refactored one. Also updated are the doc and size_t changes. |
Lots of small changes. - Rename idx_t to idx_entry as it's confusing having both idx and idx_t structs. - Remove the +1 allocation to idx_t as all we're using the next entry for is to track the cumulative file offset. Added this to idx instead. - Allocidx now returns an idx_t instead of the original (unchanged) pointer. This makes the calling code far simpler to read by removing many layers of redirection. - Allocidx also increments in->n to avoid needing two different places to increment this in the calling code. - Fixed some of the long lines so they're more readable on an 80 column wide terminal, which is the norm for most of samtools/htslib coding. - Restructured the write_output function so the on-the-fly index code is, more or less, all in the one block rather than in many separate blocks. - Write_output now uses literal strings instead of reusing the kstring for "+\n". - Removed EXIT_SUCCESS and EXIT_FAILURE bloat from write_line, as it hugely simplifies the calling code (shorter lines and easier error checking). They were never intended for internal functions anyway, but for returning from main().
|
I just rebased it and amended the commit message (I could have squashed, but it's probably easier to keep the history given it's to and fro between two authors). I will merge once the checks pass. Thanks. |
faidx/fqidx index can be created along with sequence retrieval, with --write-index option
uncompressed output was created through FILE handle and updated to create using BGZF uncompressed writing
fqidx man page update for compressed output creation details.
Resolves #2118