gwseq-io 0.2.0

Rust library for processing bigWig, bigBed, BAM and HiC files
Documentation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
//! bbi headers: the common header, the zoom ladder, the total summary, autoSql.
//!
//! Field offsets are Supp. Tables
//! 5–7 of the bigWig paper; they are spelled out here rather than derived from
//! a struct layout, because the format's fields are not aligned and a
//! `#[repr(C)]` view of them would be undefined behaviour on half of them.

use crate::error::{Error, Result};
use crate::source::ByteSource;

/// Byte sizes of the fixed-width blocks (Supp. Tables 4–8, 14).
pub const BBI_HEADER_SIZE: u64 = 64;
pub const ZOOM_HEADER_SIZE: u64 = 24;
pub const TOTAL_SUMMARY_SIZE: u64 = 40;
pub const CHR_TREE_HEADER_SIZE: u64 = 32;
/// The R-tree index header the tree nodes follow.
pub const DATA_TREE_HEADER_SIZE: u64 = 48;

pub const BIGWIG_MAGIC_SWAPPED: u32 = 0x26FC_8F88;
pub const BIGBED_MAGIC_SWAPPED: u32 = 0xEBF2_8987;
pub const CHR_TREE_MAGIC: u32 = 0x78CA_8C91;
pub const CHR_TREE_MAGIC_SWAPPED: u32 = 0x91CA_8C78;
pub const DATA_TREE_MAGIC: u32 = 0x2468_ACE0;
pub const DATA_TREE_MAGIC_SWAPPED: u32 = 0xE0AC_6824;

pub const BBI_MIN_VERSION: u16 = 3;
pub const BBI_OUTPUT_VERSION: u16 = 4;

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum BbiKind {
    BigWig,
    BigBed,
}

impl BbiKind {
    pub fn as_str(self) -> &'static str {
        match self {
            BbiKind::BigWig => "bigwig",
            BbiKind::BigBed => "bigbed",
        }
    }

    pub fn is_bigbed(self) -> bool {
        matches!(self, BbiKind::BigBed)
    }
}

#[derive(Debug, Clone)]
pub struct BbiHeader {
    pub kind: BbiKind,
    pub version: u16,
    pub zoom_levels: u16,
    pub chr_tree_offset: u64,
    pub full_data_offset: u64,
    pub full_index_offset: u64,
    /// bigBed only.
    pub field_count: u16,
    pub defined_field_count: u16,
    pub auto_sql_offset: u64,
    pub total_summary_offset: u64,
    /// 0 when the file stores its blocks uncompressed.
    pub uncompress_buffer_size: u32,
}

#[derive(Debug, Clone, Copy)]
pub struct ZoomHeader {
    pub reduction_level: u32,
    pub data_offset: u64,
    pub index_offset: u64,
}

/// The summary block is optional: a file that carries none sets its offset to
/// 0, and the reader leaves this as it stands rather than parsing 40 bytes of
/// whatever happens to sit at the start of the file. Hence NaN extremes.
#[derive(Debug, Clone, Copy)]
pub struct TotalSummary {
    pub bases_covered: u64,
    pub min_value: f64,
    pub max_value: f64,
    pub sum_data: f64,
    pub sum_squared: f64,
}

impl Default for TotalSummary {
    fn default() -> Self {
        Self {
            bases_covered: 0,
            min_value: f64::NAN,
            max_value: f64::NAN,
            sum_data: 0.0,
            sum_squared: 0.0,
        }
    }
}

#[derive(Debug, Clone, Copy)]
pub struct ChrTreeHeader {
    pub block_size: u32,
    pub key_size: u32,
    pub val_size: u32,
    pub item_count: u64,
}

/// Read the common header.
///
/// A byte-swapped magic is recognised and **refused**: these readers do not
/// swap, and saying so is more useful than reporting the file as unrecognisable.
pub fn read_header(source: &dyn ByteSource) -> Result<BbiHeader> {
    let buf = source.read_exact_at(0, BBI_HEADER_SIZE as usize)?;
    let path = source.path();
    let mut c = crate::bytes::LeCursor::new(&buf, 0, path);

    let magic = c.read_u32()?;
    let kind = match magic {
        super::BIGWIG_MAGIC => BbiKind::BigWig,
        super::BIGBED_MAGIC => BbiKind::BigBed,
        BIGWIG_MAGIC_SWAPPED | BIGBED_MAGIC_SWAPPED => {
            return Err(Error::format(path, "incompatible endianness"))
        }
        _ => return Err(Error::format(path, "not a bigwig or bigbed file")),
    };

    let version = c.read_u16()?;
    if version < BBI_MIN_VERSION {
        return Err(Error::format(
            path,
            format!("bigwig or bigbed version {version} unsupported (>= {BBI_MIN_VERSION})"),
        ));
    }

    Ok(BbiHeader {
        kind,
        version,
        zoom_levels: c.read_u16()?,
        chr_tree_offset: c.read_u64()?,
        full_data_offset: c.read_u64()?,
        full_index_offset: c.read_u64()?,
        field_count: c.read_u16()?,
        defined_field_count: c.read_u16()?,
        auto_sql_offset: c.read_u64()?,
        total_summary_offset: c.read_u64()?,
        uncompress_buffer_size: c.read_u32()?,
        // Eight reserved bytes follow at offset 56. UCSC reads them as an
        // extension offset on version 4 files and follows the pointer when they
        // are not zero, so they cannot be repurposed.
    })
}

pub fn read_zoom_headers(source: &dyn ByteSource, count: u16) -> Result<Vec<ZoomHeader>> {
    if count == 0 {
        return Ok(Vec::new());
    }
    let len = count as u64 * ZOOM_HEADER_SIZE;
    let buf = source.read_exact_at(BBI_HEADER_SIZE, len as usize)?;
    let mut c = crate::bytes::LeCursor::new(&buf, BBI_HEADER_SIZE, source.path());
    let mut headers = Vec::with_capacity(count as usize);
    for _ in 0..count {
        let reduction_level = c.read_u32()?;
        c.skip(4)?; // reserved
        headers.push(ZoomHeader {
            reduction_level,
            data_offset: c.read_u64()?,
            index_offset: c.read_u64()?,
        });
    }
    Ok(headers)
}

/// Read the total summary, or the default when the file declares none.
pub fn read_total_summary(source: &dyn ByteSource, offset: u64) -> Result<TotalSummary> {
    if offset == 0 {
        return Ok(TotalSummary::default());
    }
    let buf = source.read_exact_at(offset, TOTAL_SUMMARY_SIZE as usize)?;
    let mut c = crate::bytes::LeCursor::new(&buf, offset, source.path());
    Ok(TotalSummary {
        bases_covered: c.read_u64()?,
        min_value: c.read_f64()?,
        max_value: c.read_f64()?,
        sum_data: c.read_f64()?,
        sum_squared: c.read_f64()?,
    })
}

/// The R-tree index header that precedes the tree nodes.
///
/// Only the magic is checked here: everything else the walk needs, it reads off
/// the nodes themselves.
pub fn check_data_tree_magic(source: &dyn ByteSource, offset: u64) -> Result<()> {
    let buf = source.read_exact_at(offset, 4)?;
    let magic = u32::from_le_bytes([buf[0], buf[1], buf[2], buf[3]]);
    match magic {
        DATA_TREE_MAGIC => Ok(()),
        DATA_TREE_MAGIC_SWAPPED => Err(Error::format(
            source.path(),
            "incompatible endianness (data tree)",
        )),
        _ => Err(Error::format(
            source.path(),
            "invalid data tree magic number",
        )),
    }
}

/// The standard predefined BED field names, in order (Supp. Table 1), used to
/// name the columns of a file carrying no autoSql of its own.
pub const BED_FIELD_NAMES: &[&str] = &[
    "chrom",
    "chromStart",
    "chromEnd",
    "name",
    "score",
    "strand",
    "thickStart",
    "thickEnd",
    "itemRgb",
    "blockCount",
    "blockSizes",
    "blockStarts",
];

/// The bigBed column declaration. Field order matters — the first three are the
/// coordinates — so this is an `IndexMap`, not a `HashMap`.
///
/// An offset of 0 means the file carries no autoSql, and the standard BED names
/// stand in for as many columns as its header declares.
///
/// The text is one field per line, `type name1, name2;`. The parser takes
/// everything up to the first semicolon on a line as the field list, which is
/// why the writer must keep semicolons out of comments: one there would be read
/// as another field and the count would stop matching the header.
pub fn read_auto_sql(
    source: &dyn ByteSource,
    offset: u64,
    field_count: u16,
) -> Result<indexmap::IndexMap<String, String>> {
    if offset == 0 {
        return Ok(default_bed_fields(field_count as usize));
    }
    let text = read_nul_terminated(source, offset)?;

    let mut fields = indexmap::IndexMap::new();
    for line in text.lines() {
        // `<type> <names>;` — the first whitespace-delimited token is the type,
        // everything up to the first `;` is the comma-separated name list.
        let Some((head, rest)) = split_once_whitespace(line.trim()) else {
            continue;
        };
        let Some(list) = rest.split(';').next().filter(|l| !l.trim().is_empty()) else {
            continue;
        };
        if !rest.contains(';') {
            continue;
        }
        for name in list.split(',') {
            // The declaration may carry an array size — `int[blockCount] sizes`
            // — and a trailing comment; the name is the first token either way.
            let name = name.split_whitespace().next().unwrap_or("");
            if !name.is_empty() {
                fields.insert(name.to_string(), head.to_string());
            }
        }
    }

    validate_bed_field_names(&fields, source.path())?;
    if fields.len() != field_count as usize {
        return Err(Error::format(
            source.path(),
            format!(
                "field count {field_count} does not match autosql field count {}",
                fields.len()
            ),
        ));
    }
    Ok(fields)
}

/// Split on the whitespace that separates an autoSql field's type from its
/// names.
///
/// Not simply the first run of it: a real bigBed declares
/// `enum("+", "-", ".") strand;`, and `set(a, b) flags;`, whose types carry
/// spaces inside their parentheses. Whitespace inside a bracket or a quoted
/// string is part of the type, and the split is the first run outside both — so
/// the line above is the type `enum("+", "-", ".")` and the name `strand`,
/// rather than the type `enum("+",` and a field count one short.
fn split_once_whitespace(line: &str) -> Option<(&str, &str)> {
    let (mut depth, mut quote) = (0i32, None::<char>);
    let mut head_end = None;
    for (at, c) in line.char_indices() {
        match quote {
            Some(q) => {
                if c == q {
                    quote = None;
                }
            }
            None => match c {
                '"' | '\'' => quote = Some(c),
                '(' | '[' | '{' => depth += 1,
                ')' | ']' | '}' => depth -= 1,
                _ if depth <= 0 && c.is_whitespace() => {
                    head_end = Some(at);
                    break;
                }
                _ => {}
            },
        }
    }
    let (head, rest) = line.split_at(head_end?);
    Some((head, rest.trim_start()))
}

/// The autoSql text is NUL-terminated and its length is nowhere in the header,
/// so it is read in chunks until the terminator turns up.
fn read_nul_terminated(source: &dyn ByteSource, offset: u64) -> Result<String> {
    const CHUNK: usize = 4096;
    // An autoSql block describing a few hundred columns is a few kilobytes; a
    // file claiming megabytes of it is corrupt, not verbose.
    const MAX: usize = 1 << 20;
    let mut buf: Vec<u8> = Vec::new();
    loop {
        let chunk = source.read_at(offset + buf.len() as u64, CHUNK)?;
        if chunk.is_empty() {
            return Err(Error::corrupt(
                source.path(),
                offset,
                "autosql block is not NUL-terminated before the end of the file",
            ));
        }
        if let Some(at) = memchr::memchr(0, &chunk) {
            buf.extend_from_slice(&chunk[..at]);
            break;
        }
        buf.extend_from_slice(&chunk);
        if buf.len() > MAX {
            return Err(Error::corrupt(
                source.path(),
                offset,
                format!("autosql block is longer than {MAX} bytes"),
            ));
        }
    }
    String::from_utf8(buf)
        .map_err(|_| Error::corrupt(source.path(), offset, "autosql is not valid UTF-8"))
}

/// Check that the first three fields are the chromosome, the start and the end,
/// in that order.
///
/// Shared by the reader, which validates the autoSql it parses, and by the
/// writer, which validates the fields it was handed before turning them into
/// autoSql: a description this accepts on the way out must be one it accepts on
/// the way back in.
///
/// `chrom` may be spelled `chr`, since a BED written from this library names
/// its first column that, and the start and end may or may not carry the
/// `chrom` prefix. Hand-written matchers rather than regexes — three fixed
/// shapes do not need a regex engine.
pub fn validate_bed_field_names(
    fields: &indexmap::IndexMap<String, String>,
    path: &str,
) -> Result<()> {
    let bad = || Error::format(path, "missing or misplaced chr, start or end in autosql");
    if fields.len() < 3 {
        return Err(bad());
    }
    let name = |i: usize| fields.get_index(i).map(|(k, _)| k.as_str()).unwrap_or("");
    if !is_chr_name(name(0)) || !is_coord_name(name(1), "start") || !is_coord_name(name(2), "end") {
        return Err(bad());
    }
    Ok(())
}

/// `chr`, `chrom`, and either with an `id` or `name` suffix: `chr_id`,
/// `chromName`, `chrom_id`, …
fn is_chr_name(name: &str) -> bool {
    let lower = name.to_ascii_lowercase();
    let rest = match strip_chr_prefix(&lower) {
        Some(rest) => rest,
        None => return false,
    };
    let rest = rest.strip_prefix('_').unwrap_or(rest);
    matches!(rest, "" | "id" | "name")
}

/// `start`/`end`, with an optional `chr`/`chrom` prefix: `chromStart`,
/// `chrom_end`, `start`, …
fn is_coord_name(name: &str, suffix: &str) -> bool {
    let lower = name.to_ascii_lowercase();
    if lower == suffix {
        return true;
    }
    let Some(rest) = strip_chr_prefix(&lower) else {
        return false;
    };
    rest.strip_prefix('_').unwrap_or(rest) == suffix
}

fn strip_chr_prefix(lower: &str) -> Option<&str> {
    let rest = lower.strip_prefix("chr")?;
    Some(rest.strip_prefix("om").unwrap_or(rest))
}

/// The standard BED column names for a given width — `field13` and up past the
/// twelve the format names.
pub fn default_bed_fields(col_count: usize) -> indexmap::IndexMap<String, String> {
    (0..col_count)
        .map(|i| {
            let name = BED_FIELD_NAMES
                .get(i)
                .map(|n| (*n).to_string())
                .unwrap_or_else(|| format!("field{}", i + 1));
            (name, "string".to_string())
        })
        .collect()
}

/// Narrow a count or an offset to the 32 bits a bbi file stores it on.
///
/// Every coordinate, count and index the writer puts in a
/// file goes through this: the format is 32-bit throughout, and a value that
/// does not fit has to be refused where it is written rather than truncated
/// into a file that reads back as something else.
pub fn to_bbi_u32(value: i64, what: &str) -> Result<u32> {
    if !(0..=0xFFFF_FFFF).contains(&value) {
        return Err(Error::invalid(format!(
            "{what} {value} does not fit the 32 bits a bigwig or bigbed file stores it on"
        )));
    }
    Ok(value as u32)
}

/// Narrow an `f64` accumulator to the `f32` a zoom record stores it on,
/// clamping to the range `f32` can represent.
///
/// An `as` cast saturates to infinity, and a reader would carry that value
/// into every mean it computes. Clamping to the extremes keeps a summary
/// finite.
pub fn to_bbi_f32(value: f64) -> f32 {
    const LIMIT: f64 = f32::MAX as f64;
    if value > LIMIT {
        return f32::MAX;
    }
    if value < -LIMIT {
        return f32::MIN;
    }
    value as f32
}

/// The 64 bytes of the common header (Supp. Table 5), **with the magic left as
/// zero**.
///
/// The writer stamps the magic separately and last, so a write that dies partway
/// leaves a file no reader accepts rather than a header pointing at data that
/// was never written.
pub fn write_header(header: &BbiHeader) -> Result<Vec<u8>> {
    let mut out = Vec::with_capacity(BBI_HEADER_SIZE as usize);
    out.extend_from_slice(&0u32.to_le_bytes()); // magic, stamped separately
    out.extend_from_slice(&header.version.to_le_bytes());
    out.extend_from_slice(&header.zoom_levels.to_le_bytes());
    out.extend_from_slice(&header.chr_tree_offset.to_le_bytes());
    out.extend_from_slice(&header.full_data_offset.to_le_bytes());
    out.extend_from_slice(&header.full_index_offset.to_le_bytes());
    out.extend_from_slice(&header.field_count.to_le_bytes());
    out.extend_from_slice(&header.defined_field_count.to_le_bytes());
    out.extend_from_slice(&header.auto_sql_offset.to_le_bytes());
    out.extend_from_slice(&header.total_summary_offset.to_le_bytes());
    out.extend_from_slice(&header.uncompress_buffer_size.to_le_bytes());
    out.extend_from_slice(&[0u8; 8]); // reserved
    debug_assert_eq!(out.len(), BBI_HEADER_SIZE as usize);
    Ok(out)
}

/// The 24 bytes of one zoom header (Supp. Table 6).
pub fn write_zoom_header(header: &ZoomHeader) -> Vec<u8> {
    let mut out = Vec::with_capacity(ZOOM_HEADER_SIZE as usize);
    out.extend_from_slice(&header.reduction_level.to_le_bytes());
    out.extend_from_slice(&0u32.to_le_bytes()); // reserved
    out.extend_from_slice(&header.data_offset.to_le_bytes());
    out.extend_from_slice(&header.index_offset.to_le_bytes());
    out
}

/// The 40 bytes of the total summary (Supp. Table 7).
///
/// A file carrying no data writes zero extremes rather than the NaN the struct
/// defaults to: the block is present, so a reader parses it, and NaN extremes
/// would propagate into everything computed from them.
pub fn write_total_summary(summary: &TotalSummary) -> Vec<u8> {
    let empty = summary.bases_covered == 0;
    let mut out = Vec::with_capacity(TOTAL_SUMMARY_SIZE as usize);
    out.extend_from_slice(&summary.bases_covered.to_le_bytes());
    out.extend_from_slice(&if empty { 0.0 } else { summary.min_value }.to_le_bytes());
    out.extend_from_slice(&if empty { 0.0 } else { summary.max_value }.to_le_bytes());
    out.extend_from_slice(&summary.sum_data.to_le_bytes());
    out.extend_from_slice(&summary.sum_squared.to_le_bytes());
    out
}

/// The field types an autoSql declaration may name.
pub const BED_FIELD_TYPES: &[&str] = &["string", "int", "uint", "float"];

/// An autoSql declaration and the two counts the header states alongside it.
#[derive(Debug, Clone)]
pub struct BedAutoSql {
    pub text: String,
    pub field_count: u16,
    /// How many leading columns are the standard BED ones, which is what tells
    /// a browser a bigBed is a bed4, a bed6 or a bed12.
    pub defined_field_count: u16,
}

/// What the standard columns are for, as UCSC describes them (Supp. Table 2).
///
/// A browser shows these, so the ones it already knows are named the way it
/// expects and the rest fall back to the name the caller gave.
const BED_FIELD_COMMENTS: &[&str] = &[
    "Reference sequence chromosome or scaffold",
    "Start position in chromosome",
    "End position in chromosome",
    "Name of item",
    "Score from 0-1000",
    "+ or - for strand",
    "Start of where display should be thick",
    "End of where display should be thick",
    "Used as itemRgb as of 2004-11-22",
    "Number of blocks",
    "Comma separated list of block sizes",
    "Start positions relative to chromStart",
];

/// Turn a writer's declared columns into the autoSql its header points at.
///
/// The first three are the coordinates whatever the caller called them, and go
/// out under the standard names and types: a browser reading the file expects
/// `chrom`, `chromStart`, `chromEnd`, and the writer's own names are only ever
/// used to address the values it is handed.
pub fn build_auto_sql(fields: &indexmap::IndexMap<String, String>) -> Result<BedAutoSql> {
    validate_bed_field_names(fields, "fields")?;
    for (name, kind) in fields {
        if !BED_FIELD_TYPES.contains(&kind.as_str()) {
            return Err(Error::invalid(format!(
                "field type {kind} of {name} invalid (string, int, uint, float)"
            )));
        }
        // autoSql separates a type from its names with whitespace, its names
        // from each other with commas, and a declaration from the next with a
        // semicolon. A name carrying any of those writes a block this library's
        // own reader parses into a different number of fields than the header
        // declares — a file written here and then refused here.
        if let Some(bad) = name
            .chars()
            .find(|c| c.is_whitespace() || matches!(c, ';' | ',' | '(' | ')' | '"'))
        {
            return Err(Error::invalid(format!(
                "field name {name:?} contains {bad:?}, which autosql uses as punctuation                  (names may not hold whitespace, ';', ',', parentheses or quotes)"
            )));
        }
        if name.is_empty() {
            return Err(Error::invalid("a field name may not be empty"));
        }
    }

    let named = |index: usize, given: &str| -> String {
        if index < 3 {
            BED_FIELD_NAMES[index].to_string()
        } else {
            given.to_string()
        }
    };

    let mut text = String::from("table gwseqBed\n\"gwseq_io bed entries\"\n(\n");
    for (index, (given, kind)) in fields.iter().enumerate() {
        let name = named(index, given);
        let standard =
            index < 3 || (index < BED_FIELD_NAMES.len() && given == BED_FIELD_NAMES[index]);
        let kind = if index < 3 {
            if index == 0 {
                "string"
            } else {
                "uint"
            }
        } else {
            kind.as_str()
        };
        let comment = if standard {
            BED_FIELD_COMMENTS[index]
        } else {
            &name
        };
        text.push_str(&format!("{kind} {name};\t\"{comment}\"\n"));
    }
    text.push_str(")\n");

    let mut defined_field_count = 0u16;
    for (index, (given, _)) in fields.iter().enumerate() {
        let name = named(index, given);
        if index >= BED_FIELD_NAMES.len() || name != BED_FIELD_NAMES[index] {
            break;
        }
        defined_field_count += 1;
    }

    Ok(BedAutoSql {
        text,
        field_count: to_bbi_u32(fields.len() as i64, "fieldCount")? as u16,
        defined_field_count,
    })
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::source::testing::MemorySource;

    /// A minimal but valid 64-byte header, plus whatever follows.
    fn header_bytes(magic: u32, version: u16, zoom_levels: u16) -> Vec<u8> {
        let mut b = Vec::new();
        b.extend_from_slice(&magic.to_le_bytes());
        b.extend_from_slice(&version.to_le_bytes());
        b.extend_from_slice(&zoom_levels.to_le_bytes());
        b.extend_from_slice(&1000u64.to_le_bytes()); // chromTreeOffset
        b.extend_from_slice(&2000u64.to_le_bytes()); // fullDataOffset
        b.extend_from_slice(&3000u64.to_le_bytes()); // fullIndexOffset
        b.extend_from_slice(&9u16.to_le_bytes()); // fieldCount
        b.extend_from_slice(&6u16.to_le_bytes()); // definedFieldCount
        b.extend_from_slice(&400u64.to_le_bytes()); // autoSqlOffset
        b.extend_from_slice(&500u64.to_le_bytes()); // totalSummaryOffset
        b.extend_from_slice(&32768u32.to_le_bytes()); // uncompressBufSize
        b.extend_from_slice(&0u64.to_le_bytes()); // reserved
        assert_eq!(b.len(), 64);
        b
    }

    #[test]
    fn reads_a_bigwig_header() {
        let source = MemorySource::new(header_bytes(super::super::BIGWIG_MAGIC, 4, 0));
        let h = read_header(&source).unwrap();
        assert_eq!(h.kind, BbiKind::BigWig);
        assert_eq!(h.version, 4);
        assert_eq!(h.chr_tree_offset, 1000);
        assert_eq!(h.full_index_offset, 3000);
        assert_eq!(h.field_count, 9);
        assert_eq!(h.uncompress_buffer_size, 32768);
    }

    #[test]
    fn a_bigbed_magic_gives_the_other_kind() {
        let source = MemorySource::new(header_bytes(super::super::BIGBED_MAGIC, 4, 0));
        assert_eq!(read_header(&source).unwrap().kind, BbiKind::BigBed);
    }

    #[test]
    fn a_byte_swapped_file_is_refused_by_name() {
        for magic in [BIGWIG_MAGIC_SWAPPED, BIGBED_MAGIC_SWAPPED] {
            let source = MemorySource::new(header_bytes(magic, 4, 0));
            let err = read_header(&source).unwrap_err().to_string();
            assert!(err.contains("incompatible endianness"), "{err}");
        }
    }

    #[test]
    fn something_else_entirely_is_not_a_bbi_file() {
        let source = MemorySource::new(header_bytes(0xDEAD_BEEF, 4, 0));
        let err = read_header(&source).unwrap_err().to_string();
        assert!(err.contains("not a bigwig or bigbed file"), "{err}");
    }

    #[test]
    fn an_old_version_names_the_floor() {
        let source = MemorySource::new(header_bytes(super::super::BIGWIG_MAGIC, 2, 0));
        let err = read_header(&source).unwrap_err().to_string();
        assert!(err.contains("version 2 unsupported (>= 3)"), "{err}");
    }

    #[test]
    fn a_truncated_header_is_corrupt_not_a_panic() {
        let source = MemorySource::new(vec![0u8; 40]);
        assert!(matches!(
            read_header(&source),
            Err(crate::error::Error::Corrupt { .. })
        ));
    }

    #[test]
    fn zoom_headers_follow_the_common_one() {
        let mut b = header_bytes(super::super::BIGWIG_MAGIC, 4, 2);
        for (reduction, data, index) in [(10u32, 100u64, 200u64), (40, 300, 400)] {
            b.extend_from_slice(&reduction.to_le_bytes());
            b.extend_from_slice(&0u32.to_le_bytes());
            b.extend_from_slice(&data.to_le_bytes());
            b.extend_from_slice(&index.to_le_bytes());
        }
        let source = MemorySource::new(b);
        let zooms = read_zoom_headers(&source, 2).unwrap();
        assert_eq!(zooms.len(), 2);
        assert_eq!(zooms[0].reduction_level, 10);
        assert_eq!(zooms[1].index_offset, 400);
        assert!(read_zoom_headers(&source, 0).unwrap().is_empty());
    }

    #[test]
    fn a_summary_offset_of_zero_means_the_file_carries_none() {
        let source = MemorySource::new(vec![0u8; 8]);
        let s = read_total_summary(&source, 0).unwrap();
        assert_eq!(s.bases_covered, 0);
        assert!(s.min_value.is_nan() && s.max_value.is_nan());
    }

    #[test]
    fn reads_a_total_summary() {
        let mut b = vec![0u8; 64];
        b.extend_from_slice(&1234u64.to_le_bytes());
        b.extend_from_slice(&(-1.5f64).to_le_bytes());
        b.extend_from_slice(&9.25f64.to_le_bytes());
        b.extend_from_slice(&100.0f64.to_le_bytes());
        b.extend_from_slice(&500.0f64.to_le_bytes());
        let source = MemorySource::new(b);
        let s = read_total_summary(&source, 64).unwrap();
        assert_eq!(s.bases_covered, 1234);
        assert_eq!((s.min_value, s.max_value), (-1.5, 9.25));
        assert_eq!((s.sum_data, s.sum_squared), (100.0, 500.0));
    }

    #[test]
    fn the_data_tree_magic_is_checked_both_ways() {
        let source = MemorySource::new(DATA_TREE_MAGIC.to_le_bytes().to_vec());
        assert!(check_data_tree_magic(&source, 0).is_ok());
        let source = MemorySource::new(DATA_TREE_MAGIC_SWAPPED.to_le_bytes().to_vec());
        let err = check_data_tree_magic(&source, 0).unwrap_err().to_string();
        assert!(err.contains("incompatible endianness"), "{err}");
        let source = MemorySource::new(vec![1, 2, 3, 4]);
        assert!(check_data_tree_magic(&source, 0).is_err());
    }
}

#[cfg(test)]
mod write_tests {
    use super::*;
    use crate::source::testing::MemorySource;

    fn fields(pairs: &[(&str, &str)]) -> indexmap::IndexMap<String, String> {
        pairs
            .iter()
            .map(|(a, b)| ((*a).to_string(), (*b).to_string()))
            .collect()
    }

    #[test]
    fn a_written_header_reads_back_as_itself() {
        let header = BbiHeader {
            kind: BbiKind::BigBed,
            version: BBI_OUTPUT_VERSION,
            zoom_levels: 7,
            chr_tree_offset: 4096,
            full_data_offset: 512,
            full_index_offset: 2048,
            field_count: 6,
            defined_field_count: 4,
            auto_sql_offset: 368,
            total_summary_offset: 400,
            uncompress_buffer_size: 32768,
        };
        let mut bytes = write_header(&header).unwrap();
        // The magic is stamped separately, so the bytes as written are refused.
        let err = read_header(&MemorySource::new(bytes.clone())).unwrap_err();
        assert!(err.to_string().contains("not a bigwig or bigbed"), "{err}");

        bytes[..4].copy_from_slice(&super::super::BIGBED_MAGIC.to_le_bytes());
        let read = read_header(&MemorySource::new(bytes)).unwrap();
        assert_eq!(read.kind, BbiKind::BigBed);
        assert_eq!(read.version, header.version);
        assert_eq!(read.zoom_levels, header.zoom_levels);
        assert_eq!(read.chr_tree_offset, header.chr_tree_offset);
        assert_eq!(read.full_data_offset, header.full_data_offset);
        assert_eq!(read.full_index_offset, header.full_index_offset);
        assert_eq!(read.field_count, header.field_count);
        assert_eq!(read.defined_field_count, header.defined_field_count);
        assert_eq!(read.auto_sql_offset, header.auto_sql_offset);
        assert_eq!(read.total_summary_offset, header.total_summary_offset);
        assert_eq!(read.uncompress_buffer_size, header.uncompress_buffer_size);
    }

    #[test]
    fn a_written_zoom_header_and_summary_read_back_as_themselves() {
        let zooms = [
            ZoomHeader {
                reduction_level: 10,
                data_offset: 100,
                index_offset: 200,
            },
            ZoomHeader {
                reduction_level: 40,
                data_offset: 300,
                index_offset: 400,
            },
        ];
        let mut bytes = vec![0u8; BBI_HEADER_SIZE as usize];
        for z in &zooms {
            bytes.extend_from_slice(&write_zoom_header(z));
        }
        let read = read_zoom_headers(&MemorySource::new(bytes), 2).unwrap();
        assert_eq!(read.len(), 2);
        assert_eq!(read[1].reduction_level, 40);
        assert_eq!(read[1].data_offset, 300);
        assert_eq!(read[1].index_offset, 400);

        let summary = TotalSummary {
            bases_covered: 1234,
            min_value: -1.5,
            max_value: 9.25,
            sum_data: 1000.0,
            sum_squared: 5000.0,
        };
        let read = read_total_summary(&MemorySource::new(write_total_summary(&summary)), 0);
        // Offset 0 means "no summary" to a reader, so it is read at 1 instead.
        assert_eq!(read.unwrap().bases_covered, 0);
        let mut padded = vec![0u8];
        padded.extend_from_slice(&write_total_summary(&summary));
        let read = read_total_summary(&MemorySource::new(padded), 1).unwrap();
        assert_eq!(read.bases_covered, 1234);
        assert_eq!(read.min_value, -1.5);
        assert_eq!(read.sum_squared, 5000.0);
    }

    #[test]
    fn an_empty_summary_writes_zero_extremes_rather_than_nan() {
        let bytes = write_total_summary(&TotalSummary::default());
        let mut padded = vec![0u8];
        padded.extend_from_slice(&bytes);
        let read = read_total_summary(&MemorySource::new(padded), 1).unwrap();
        assert_eq!(read.min_value, 0.0);
        assert_eq!(read.max_value, 0.0);
    }

    #[test]
    fn built_autosql_reads_back_as_the_fields_it_describes() {
        // The round trip is the point: a writer that emitted a declaration this
        // library's own reader refuses would produce an unreadable file.
        let declared = fields(&[
            ("chr", "string"),
            ("start", "uint"),
            ("end", "uint"),
            ("name", "string"),
            ("score", "uint"),
        ]);
        let sql = build_auto_sql(&declared).unwrap();
        assert_eq!(sql.field_count, 5);
        // chr/chromStart/chromEnd/name/score are the first five standard ones.
        assert_eq!(sql.defined_field_count, 5);

        let mut bytes = vec![0u8];
        bytes.extend_from_slice(sql.text.as_bytes());
        bytes.push(0);
        let read = read_auto_sql(&MemorySource::new(bytes), 1, 5).unwrap();
        assert_eq!(
            read.keys().map(String::as_str).collect::<Vec<_>>(),
            ["chrom", "chromStart", "chromEnd", "name", "score"]
        );
        assert_eq!(read["chrom"], "string");
        assert_eq!(read["score"], "uint");
    }

    #[test]
    fn a_non_standard_column_stops_the_defined_field_count() {
        let declared = fields(&[
            ("chr", "string"),
            ("start", "uint"),
            ("end", "uint"),
            ("pvalue", "float"),
            ("score", "uint"),
        ]);
        let sql = build_auto_sql(&declared).unwrap();
        // Three coordinates, then a column BED does not know.
        assert_eq!(sql.defined_field_count, 3);
        assert_eq!(sql.field_count, 5);
        // And it keeps its own name and type, with itself as the comment.
        assert!(
            sql.text.contains("float pvalue;\t\"pvalue\"\n"),
            "{}",
            sql.text
        );
    }

    #[test]
    fn a_bad_field_name_or_type_is_refused_before_anything_is_written() {
        let err = build_auto_sql(&fields(&[
            ("name", "string"),
            ("start", "uint"),
            ("end", "uint"),
        ]))
        .unwrap_err()
        .to_string();
        assert!(err.contains("missing or misplaced"), "{err}");
        let err = build_auto_sql(&fields(&[
            ("chr", "string"),
            ("start", "uint"),
            ("end", "uint"),
            ("x", "double"),
        ]))
        .unwrap_err()
        .to_string();
        assert!(err.contains("field type double"), "{err}");
    }

    #[test]
    fn a_coordinate_past_32_bits_is_refused_rather_than_truncated() {
        assert_eq!(to_bbi_u32(4_294_967_295, "chromEnd").unwrap(), u32::MAX);
        let err = to_bbi_u32(4_294_967_296, "chromEnd")
            .unwrap_err()
            .to_string();
        assert!(err.contains("chromEnd 4294967296"), "{err}");
        assert!(to_bbi_u32(-1, "chromStart").is_err());
    }

    #[test]
    fn narrowing_a_summary_clamps_rather_than_reaching_infinity() {
        assert_eq!(to_bbi_f32(1.0e300), f32::MAX);
        assert_eq!(to_bbi_f32(-1.0e300), f32::MIN);
        assert_eq!(to_bbi_f32(1.5), 1.5f32);
    }
}

#[cfg(test)]
mod autosql_tests {
    use super::*;

    fn fields(pairs: &[(&str, &str)]) -> indexmap::IndexMap<String, String> {
        pairs
            .iter()
            .map(|(k, v)| ((*k).to_string(), (*v).to_string()))
            .collect()
    }

    /// A real bigBed declares `enum("+", "-", ".") strand;`, whose type carries
    /// spaces inside its parentheses. Split on the first whitespace, that line
    /// gave the type `enum("+",` and one field too few — and the whole file was
    /// then refused for a field-count mismatch.
    #[test]
    fn a_type_with_spaces_inside_its_brackets_is_one_token() {
        for (line, kind, rest) in [
            (
                "enum(\"+\", \"-\", \".\") strand;",
                "enum(\"+\", \"-\", \".\")",
                "strand;",
            ),
            ("set(a, b) flags;", "set(a, b)", "flags;"),
            (
                "int[blockCount] blockSizes;",
                "int[blockCount]",
                "blockSizes;",
            ),
            ("uint  score;", "uint", "score;"),
            ("string name; comment", "string", "name; comment"),
        ] {
            assert_eq!(split_once_whitespace(line), Some((kind, rest)), "{line}");
        }
        // Nothing to split on at all.
        assert_eq!(split_once_whitespace("uint"), None);
    }

    /// A name carrying autoSql's own punctuation writes a block this reader
    /// parses into a different field count than the header declares — a file
    /// written here and refused here.
    #[test]
    fn a_field_name_may_not_hold_autosql_punctuation() {
        for bad in [
            "two words",
            "semi;colon",
            "com,ma",
            "paren(s)",
            "quo\"te",
            "",
        ] {
            let f = fields(&[
                ("chr", "string"),
                ("start", "uint"),
                ("end", "uint"),
                (bad, "string"),
            ]);
            assert!(
                build_auto_sql(&f).is_err(),
                "{bad:?} was accepted as a field name"
            );
        }
        let ok = fields(&[
            ("chr", "string"),
            ("start", "uint"),
            ("end", "uint"),
            ("itemRgb", "string"),
        ]);
        assert!(build_auto_sql(&ok).is_ok());
    }
}