Documenting `sample` and `sample2` behaviour

[pathpipe]

Hello, thanks for developing seqkit! Opening this issue because I observed an unexpected number of sequences returned when using the -p option of sample and sample2 commands of seqkit v2.13.0:

E.g.

wget https://github.com/shenwei356/seqkit/raw/refs/tags/v2.13.0/tests/reads_1.fq.gz  # 2500 reads
seqkit sample2 -p 0.01 reads_1.fq.gz > test_output.fq  # Expect 25 reads
[INFO] 20 sequences outputted

I realise that you know about this behaviour because it's described here #518, but it's potentially misleading - I was certainly surprised when sampling the records from a FASTA/Q did not get the expected number of records back. After asking Claude AI how the functions work, I think it's just a case of carefully documenting the use cases.

I hope I'm right in assuming the sample command is now deprecated in favour of sample2. So I will focus on sample2 behaviour. This is what Claude said:

seqkit sample2

Designed specifically to provide exact count, uniform probability, low memory. It picks one of four branches:

1. By number, no -2: loads all records into memory, then does a partial Fisher-Yates shuffle — for i in [0, number), swap records[i] with records[j] where j = i + rand.Int63n(totalSeqs-i) — and outputs the first number records. Equivalent to an unbiased random permutation truncated to the sample size.
2. By proportion, no -2: same Bernoulli streaming trial as sample's proportion mode — approximate count, not exact.
3. By proportion + -2: first pass counts total records, converts proportion to an exact number = floor(seqNum * proportion), then falls into the reservoir path below.
4. By number + -2 (or derived from case 3): classic reservoir sampling via index selection rather than swapping records in a live reservoir array:

• Initialize savedIndices = [0, 1, ..., number-1].
• For each subsequent index i from number to seqNum-1: draw j = rand.Int63n(i+1); if j < number, overwrite savedIndices[j] = i.
• This is the standard Algorithm R reservoir sampling — each record ends up with equal probability number/seqNum of being selected.
• Sort the selected indices, then stream through the file a second time, emitting records whose index matches the sorted savedIndices list in order.

If this is accurate, my understanding is:

• If you don't care about an exact number of reads but want them to be proportional to input (only an approximate proprotion), use -p option
• If you care about an exact number of reads and proportionality (as well as memory efficiency), use -p -2 option
• If you want an exact number of reads regardless of input, use -n option.
• If you want an exact number of reads and memory efficiency (for large input), use -n -2 option.

As an example, tested the 2 pass for proportion and got the expected number of reads.

seqkit sample2 -2 -p 0.01 reads_1.fq.gz > test_output.fq  # Expect 25 reads
[INFO] first pass: counting seq number
[INFO] seq number: 2500
[INFO] sample 25/2500 by proportion (0.010000)
[INFO] second pass: reading and sampling
[INFO] 25 sequences outputted

For reference, seqtk seems to have some inexact behaviour too, but sample2 is an improvement on this if well documented (especially as the 2-pass mode doesn't work for proportions):

$ seqtk sample reads_1.fq.gz 0.01 | wc -l
96
$ seqtk sample reads_1.fq.gz 0.009 | wc -l
88
$ seqtk sample reads_1.fq.gz 0.011 | wc -l
100
$ seqtk sample -2 reads_1.fq.gz 0.01 | wc -l
[W::stk_sample] when sampling a fraction, option -2 is ignored.96

Many thanks!