One thing that is noticable on fasta data is the compressability of DNA sequence is just poor with deflate.

On human chromosome 1 libdeflate maximum compression gets 1.94 bits per base. It's barely better than a naive 2-bit encoding. Zlib's best is only 2.2bpb.

I plotted both zlib (top) and libdeflate (bottom) together. Each point is a compression level from 1 to 9. We see level 4 on libdeflate is a sharp corner and beyond that it's really not worth while trying harder. With zlib (which ideally no one would use) it's somewhere between 4, 5 or 6. I think therefore the default compression level should perhaps be different here.

Also, I noticed saying --output-format-option level=0 is equivalent to level=6 rather than uncompressed. Maybe this is generic and happens elsewhere too? It doesn't with samtools view but that could well be writing the hfile arguments to work around the problem. (Edit: thinking about it it's kind of irrelevant anyway as there's no benefit to writing uncompressed bgzf - you'd just do what we already do and not add the .gz suffix. So it's a quirk, but not a problem provided to change the man page to say level 1 to 9 instead of 0 to 9.)

Curiously, if you drive zlib correctly then this also turns things on its head. zlib can be better than libdeflate which can be better than zstd! I expect zstd -19 would win on size, and there are lots of options in zstd to tweak too. Note also this is whole file compression rather than block based, but it's just an easy test to implement.

The key is just accepting that LZ matches are hard to find in fasta files so there's little benefit (unless we're printing up many entries and they're all related). So RLE mode is usually best bang-for-buck. Zlib lets us set this directly, and ironically it offers far better compression ratios than the default strategy, even at maximum compression level.

      		       	    size       user time
zlib Z_DEFAULT_STRATEGY	-1  79259417    4.27
zlib Z_DEFAULT_STRATEGY	-6  68323107   29.91
zlib Z_DEFAULT_STRATEGY	-9  65413881  170.36
zlib Z_HUFFMAN_ONLY	    67423830    3.37
zlib Z_RLE		    64770965	3.94 <<<
libdeflate gzip -1	    73115740	1.92
libdeflate gzip -3	    65588478	5.79
libdeflate gzip -4	    62213062	8.32
libdeflate gzip -6	    61686376   16.35
libdeflate gzip -9	    61470656   70.74
libdeflate gzip -11	    60058576   67.89
zstd -1     		    71152620    1.67
zstd -3	    		    70563644	2.15
zstd -5	    		    69284690	4.95
zstd -9	    		    67189974   15.10
zstd -12		    65468432   39.34
zstd -15		    61001117  130.88
zstd -19		    56960008  241.22
xz -1			    69552440   23.44
xz -3			    67042544   46.98
xz -5			    56234960  251.54
xz -9			    54810092  368.98
ransx16pr -o4		    60466010    1.01
ransx16pr -o5		    58583693	1.30 <<< CRAM entropy encoder

I'm not sure what we should do with this information though. Probably just ignore it. Our recommended approach is libdeflate anyway and that basically means just sticking with a modest compression level as the default to hit the sweet spot.

Edit: added a few more entries, including a demonstration of what pure entropy encoding gives instead of LZ encoding. We can see hunting for LZ matches does help in the end (zstd -19 or xz -5 onwards), but it's extremely costly. I still think light weight compression (given this will be block based anyway so long-range matches wouldn't ever exist) is the correct setting.

It looks good, but I still think the default compression level ought to be 4 rather than 6. The default is just a really bad trade-off in CPU vs size for fasta data.

@ seq4[samtools.../samtools]130; for l in `seq 1 9`;do echo -e -n "$l\t";/usr/bin/time -f '%U' samtools faidx --output-fmt-option level=$l -@4 $HREF38 chr1 -o /tmp/_.gz 2>&1 | tr '\012' '\011';wc -c < /tmp/_.gz;done
1	2.23	73528585
2	3.27	71830032
3	5.16	65682361
4	8.91	62372818
5	16.09	61960402
6	35.22	62021583
7	51.57	62072212
8	53.30	61391026
9	88.69	61043614

@ seq4[samtools.../samtools]; /usr/bin/time -f '%U' samtools faidx  -@4 $HREF38 chr1 -o /tmp/_.gz 2>&1 | tr '\012' '\011';wc -c < /tmp/_.gz
35.08	62021583

Suggestion:

diff --git a/faidx.c b/faidx.c
index 3fb70ca4..dc2afd9a 100644
--- a/faidx.c
+++ b/faidx.c
@@ -456,7 +456,7 @@ int faidx_core(int argc, char *argv[], enum fai_format_options format)
         } else {
             hts_opt *opts = (hts_opt *)(out.gopt->out.specific);
             char mode[13] = "w";
-            int level = -1;
+            int level = 4;
             while (opts) {
                 if (opts->opt == HTS_OPT_COMPRESSION_LEVEL) {   //compression level
                     level = opts->val.i;

Set the default level as 4 instead of -1/zlib default.
Man page updated for the same.

There's some trailing white space in the man page that I was about to just force push up prior to merge (we don't need a discussion on it) and one minor slight rewording about default compression levels to avoid an implication that the default output style is compressed.

However I also spotted something with the kstring so I'll comment on that elsewhere. The man page change was this:

diff --git a/doc/samtools-faidx.1 b/doc/samtools-faidx.1
index fffb0add..99497710 100644
--- a/doc/samtools-faidx.1
+++ b/doc/samtools-faidx.1
@@ -58,7 +58,7 @@ retrieved and printed to stdout in the FASTA format.
 
 The input and output can be files compressed in the
 .B BGZF
-format. The default output compression level is 4.
+format. When output is compressed, the default compression level is 4.
 
 The sequences in the input file should all have different names.
 If they do not, indexing will emit a warning about duplicate sequences and
@@ -74,7 +74,7 @@ any extracted subsequence will be in FASTA format.
 .TP 8
 .BI "-o, --output " FILE
 Write FASTA to file rather than to stdout. With .gz, .bgz, .bgzf file extensions
-the output will be 
+the output will be
 .B BGZF
 compressed.
 .TP

(You can't see the "the output will be " space, but it's there.) Trivial and we don't strictly check man pages, only code, but I spotted it in diff so thought I'd just tidy it during merging.

updated makefile and made changes on kstring* to kstring.