The test/test_expr tool is both a battery of unit tests (when run without arguments), but also a general purpose expression calculator.

Examples:

$ ./test/test_expr '2*3 + 4*5'
26
$ ./test/test_expr '3&2 > 5&1';  # This differs to C precedence.  It is equiv to (3&2) > (5&1)
1
$ ./test/test_expr '"foo" =~ "[aeiou]{2}"'
1
$ ./test/test_expr '"a\\b\"c"'
a\b"c
$ ./test/test_expr '17/8'
2.125

Things that may need work.

• Windows. Specifically mingw doesn't seem to come wth a regular expression library. See https://stackoverflow.com/questions/45242741/how-use-regex-library-using-c-with-mingw-gcc-on-linux

• Naming of fields. This match the SAM source, but maybe we should explicitly match the SAM specification. For example "mqual" to "mapq" and "isize" to "tlen"? Thoughts? Alternatively we could go with much longer more descriptive names (as is done in sambamba).

I think we should use the names in the specification, as the ones in the source are probably less well known. Or possibly we could allow aliases.

The lack of a regex library in Windows is a bit annoying. We could paper over it by using PCRE, but I suspect it might result in regexps that work differently between Windows and everything else. Other alternatives available in MSYS2 might be TRE or Henry Spencer's library.

For now I'll just do some configure checks and get it to build without regex support. It's useful anyway, and people can kludge in pcre themselves via CPPFLAGS and LDFLAGS if they wish in which case configure will reenable the support.

Actually, bcftools uses regcomp, and has a Windows solution already. So it doesn't look to be a big problem after all.

Updates.

• Check for regcomp function (-lregex) for Mingw.

• Rename isize to tlen and mqual to mapq (as per SAM)

• Add rnext and pnext keywords (but keeping mrname and mpos as aliases), also as per SAM.

• Renamed tid to refid and mtid to mrefid. Loosely based on BAM spec, although that has next_refid. I don't like that as we're mixing both ordering and with/without underscores: pnext vs next_refid. There's no obvious ?next phrase to go along side pnext and rnext, hence keeping the "m" prefix terms.

While this is specifically for SAM filtering should we be designing it with an eye for adding VCF filters later, even if we don't want to add them now? For example, we could have hts_filter_t instead of sam_filter_t, just filter for the filter-setting option and HTS_OPT_FILTER for the enumerated type. Attempting to set a filter on VCF files could be met with a polite refusal until we get around to implementing it.

I actually started with HTS_OPT_FILTER but renamed it because the filter is explicitly for SAM with sam keywords. Indeed it was thinking about other format filters that made me rename it to HTS_OPT_SAM_FILTER.

Good point about data types though. I'll rename those.

The design was always to separate language from data-specific keywords. Hence expr.c as its own file.

Note however the bcftools filtering is far more powerful, with functions, arrays, etc.

Wow — I thought I saw this gestating as a samtools-specific thing…

Is the reason for moving this to htslib to make the functionality available to other programs if they added appropriate function calls, or to make it available to other programs (via user-settable format options on the command line) even without their knowledge or changes being made to their code?

The latter might not work too well as not too many programs make input file format options accessible to the user, so they would need source code changes to enable this anyway.

The former (i.e., making this available via an additional API call like int sam_passes_filter(const bam1_t *, filter_expr_something_t *) that programs would add explicitly into their read loops) has the advantage that it is more explicit and more flexible.

(Example of the lack of flexibility of the latter within-sam_read1() approach: samtools view would not be able to use this as a filtering option whose rejected reads are diverted to the -U file.)

At present, the PR implements the latter via hts_set_opt() and adding a rather hacky loop inside sam_read1(). Fortunately just about everything, including most iterators, does its reading via sam_read1(). However filters would not be applied for e.g. a CRAM file read via a multi-region iterator (i.e., when cram_readrec() is used).

Re naming, and potentially implementing this for non-SAM htsFiles in future: the hts_set_filter_expression() function added takes an expression as a string and does not appear to be SAM-specific. So it would appear to be (in principle) usable for SAM/VCF/whatever and users might expect it to internally check the fp's file format and parse the given expression appropriately (or return “not implemented for this format”). Similarly the hts_opt_apply() option just takes a string, independent of the file type.

So if pursuing the implicit filtering-as-format-option approach, I think either there should be HTS_OPT_SAM_FILTER + sam_set_filter_expression() and equivalents in due course for VCF; or there should be HTS_OPT_FILTER + hts_set_filter_expression() that check applicability in the innermost implementation (i.e., inside hts_set_filter_expression()). Probably the latter, as the argument here is just a string, not a sam_filter_expr_t or something. At the moment, there's a mixture of the two.

Thanks for the comments. I agree with the general thrust of it and I think it's relatively easy modifications to make it a little more flexible so it can encompass both the way it works now as well as a more programmatic API instead.

Wow — I thought I saw this gestating as a samtools-specific thing…

Is the reason for moving this to htslib to make the functionality available to other programs if they added appropriate function calls, or to make it available to other programs (via user-settable format options on the command line) even without their knowledge or changes being made to their code?

It started off a week ago in samtools as a draft. However there are considerable advantages to putting it into htslib.

• We don't need to make adjustments to every single samtools sub-command to permit filtering options. Everything that currently has an input-fmt-option magically gains this ability, even including test_view.
• We don't need to use samtools view | samtools sub-command2 as a way of filtering. That can have advantages (and it's still an option), but it also has disadvantages.
• It's available to other users as htslib.

Note this doesn't mean it won't also be available by easier and friendlier command line options. I'll probably still keep something like samtools view -z "filter-expression".

The latter might not work too well as not too many programs make input file format options accessible to the user, so they would need source code changes to enable this anyway.

All of our programs do. Any third party tool that doesn't make these options available isn't going to break. It simply won't be able to support filtering. That'd be true though if we added it to samtools view via its own dedicated command line option too, so I don't see any difference there.

The former (i.e., making this available via an additional API call like int sam_passes_filter(const bam1_t *, filter_expr_something_t *) that programs would add explicitly into their read loops) has the advantage that it is more explicit and more flexible.

What you're describing is already how it works, but it's just a matter of deciding whether we wish to expose that interface. Right now the option is calling hts_set_filter_expression(fp, expr). We could make that public as a friendlier programmatic interface.

That in turn uses sam_filter_t *sam_filter_init(const char *str) and int sam_filter_eval(...) in expr.c that creates a filter and checks whether a bam record passes that filter. Again we could make these public if desired. I think they need renaming though from sam to hts and create a second sam_filter function that explicitly passes in the bam_sym_lookup function.

(Example of the lack of flexibility of the latter within-sam_read1() approach: samtools view would not be able to use this as a filtering option whose rejected reads are diverted to the -U file.)

Why does it have to be either/or?

Yes you have more flexibility by calling sam_read1 yourself and then optionally calling e.g. sam_passes_filter on one or more sets of filters. However you lose the automatic filtering that can be added to any tool without needing to change code. I don't see it as one or the other. It would make sense to be able to specify and automatic filter, as well as provide the ability to do more custom on-the-fly ones precisely where you wish in your code.

At present, the PR implements the latter via hts_set_opt() and adding a rather hacky loop inside sam_read1(). Fortunately just about everything, including most iterators, does its reading via sam_read1(). However filters would not be applied for e.g. a CRAM file read via a multi-region iterator (i.e., when cram_readrec() is used).

Gah, I chose sam_read1 as I thought it was the one central place that everything went via. I wasn't aware the cram iterator had it's own reading code that avoided that function. :-(

Obviously I need to fix that.

I know it's a little hacky by design. The choice was either a few lines added to the top and end, a mammoth reindentation of the entire function making PR checking hard, or restructuring into multiple functions (that's probably long overdue tbh!). I went with the simplest to review, but I'll consider a wider restructure then.

Re naming, and potentially implementing this for non-SAM htsFiles in future: the hts_set_filter_expression() function added takes an expression as a string and does not appear to be SAM-specific. So it would appear to be (in principle) usable for SAM/VCF/whatever and users might expect it to internally check the fp's file format and parse the given expression appropriately (or return “not implemented for this format”). Similarly the hts_opt_apply() option just takes a string, independent of the file type.

Agreed, and Rob's comments above were similar. I deliberately designed it to have a non-SAM specific component with a view to using in other formats in the future, but somewhere along the line forgot this in my function naming. At present it's too simple to replace the bcftools filtering, but it doesn't have to arrive fully complete on day 1 as long as the design doesn't forbid future work.

I'll edit expr.c to use hts prefixes instead of sam prefixes. Good point on doing introspection on the file format type to decide which format specific symbol-lookup function to pass in. That may mean going back to the generalised option is better again. I'll make the other changes and then see how it's looking.

Hopefully this now addresses all the comments. Thanks for raising them.

It's now more public than it was before, to facilitate being able to call in and do your own filtering where and when you want instead of being required to use the hts_opt string interface. This exposes more, so I've done a lot of data type and function renaming.

The general structure is:

* denotes new public API

hts_expr.c    	  General filter expression language
       *hts_filter_init	       Compile
       *hts_filter_free	       Free
       *hts_filter_eval	       Validate

htslib/hts.h
	"hts_filter_t *filter;" a new member of htsFile

htslib/hts_expr.h  Header for things exposed in hts_expr.c.  Automatically included by htslib/hts.h
	hts_filter_* funcs above
	struct hts_filter_t hts_filter_t	Opaque type
	typedef struct {...} hts_expr_val_t

hts.c
	hts_opt_add / hts_opt_apply: "filter" -> HTS_OPT_FILTER
	hts_set_opt: hts_set_filter_expression call
	hts_close: hts_filter_free call
      * hts_set_filter_expression (easier programmatic interface than hts_opt_add)

sam.c
	bam_sym_lookup		Handles SAM columns; callback for hts_filter
       *sam_passes_filter	Wraps up hts_filter_eval + bam_sym_lookup
	sam_read1		Calls hts_filter_eval
	cram_readrec		Calls hts_filter_eval

I need some help to make it work on Windows though.

I don't understand the nuances of our symbol export code. I've added HTSLIB_EXPORT macros, but it's still not sufficient for Windows. I included hts_defs.h too which defines it, but that's not enough either. Guidance please!

The error is:

Checking shared library exports
hts_filter_eval
hts_filter_free
hts_filter_init
Error: Found unexported symbols (listed above)

Thoughts on some of the terms we filter on.

• library (this doesn't exist yet). Samtools can filter on this, but we can only do RG.

• rlen is from bam_cigar2rlen so it's simply the number of reference bases consumed by CIGAR.

• My qlen is quite simply length(seq). If we had functions we could do that explicitly.

  This is not the same as bam_cigar2qlen however as unmapped data returns 0 there. Otherwise obviously it would be.

  Note samtools view -m documents qlen as the length of query sequence consumed by CIGAR, but inexplicably explicitly codes it to also include hard-clips. There's no comment or obvious reason why. Maybe it's just a long standing bug?

• seq and qual (don't exist yet). There's not a huge amount of point in this, but potentially people may wish to (inefficiently) select on them.

Should I add seq_len as a token, or bite the bullet and do basic functions, although that's more work. I think it'd make more sense though for qlen and rlen to be defined consistently.

There are other potentially useful functions. Eg we could have min(qual), max(qual), avg(qual), etc. I was delaying that for another PR though given the bigger they are the longer (much much longer) they ever take to be reviewed and merged and it doesn't all have to be done from day 1.

In absence of any comment either way, I decided to see how complex it is to add basic functions and it's easy enough. Not an entirely clean design, but good enough for now. (It accepts length("foo") but not length ("foo").)

So I implemented the list of thoughts above. They're not brilliantly efficient, but if you really want to do filtering on sequence by regexp then you now can! By quality may be a tad heroic though...

I think your appveyor failure here is caused by a lack of #define HTS_BUILDING_LIBRARY in the new hts_expr.c file. It's needed on Windows to set up the macros that export the symbols, as opposed to importing them.

See for example https://github.com/samtools/htslib/blob/238600/hts.c#L26

Also, please could you make the list of header files in the Makefile dependency list be in the same order as they are in the .c file? It makes comparing them easier.

Why does it have to be either/or?

Yes you have more flexibility by calling sam_read1 yourself and then optionally calling e.g. sam_passes_filter on one or more sets of filters. However you lose the automatic filtering that can be added to any tool without needing to change code.

It doesn't have to be either/or, but providing automatic filtering inside sam_read1() does involve making rather hacky changes and diluting its core purpose of reading one SAM record (the clue is, as they say, in the name). The HTSlib API is a very explicit API (consider e.g. the addition of an explicit header argument to sam_read1() compared to the old samtools API function), and this is very implicit filtering.

However with the explicit version added and enabling more advanced usage, FWIW I'm happy enough either way.

The splitting out of sam_read1() into subfunctions looks clearer; some more comments:

• Perhaps the new helpers could be static inline.
• The cram_readrec() changes would benefit from looking more like those in sam_read1(): s/filtered/pass_filter/ and needs to call sam_passes_filter() only if there's an active filter.
• sam_read1() case sam can lose the extra braces now.

The problem is we need a transparent way to retrofit this to other tools. I don't think requiring code changes to them is justified purely by the desire to be tidy.

Also why is the clue in the name? sam_read1 is named thus not because it reads 1 record, but because it's version 1 of the API (as in returning bam1_t). Or at least that's what I thought it meant, although I haven't ever seen it spelt out anywhere.

I don't think the core purpose has changed either. It does fetch a single record, one at a time. It's simply now one record matching your criteria if you specified them. At it's heart nothing has changed, but optionally you have more control over it. I don't see that as a bad thing.

Thanks for comments on tidying various bits up though. I'll do these today.

Thanks John, and great catch on the lack of check in cram_readrec. That lead to a crash. (Note to self - we need to do some coverage analysis and work out which paths are not exercised by the test harness at some point.)

Hopefully I've addressed all the issues now.

Pushed some new changes. Thanks for the comments.

Rebased, squashed and merged as 999d181 to c7c433f