diff --git a/.github/workflows/unittest.yml b/.github/workflows/unittest.yml index 031c8cb..63d767a 100644 --- a/.github/workflows/unittest.yml +++ b/.github/workflows/unittest.yml @@ -16,6 +16,9 @@ jobs: julia-version: ['1', '1.11'] os: [ubuntu-latest, macOS-latest, windows-latest] experimental: [false] + # Runner's native architecture. Must be a matrix key, so that the x86 entry below + # overrides it and becomes a new job, instead of being merged into an existing job. + arch: [default] include: # Include nightly, but experimental, so it's allowed to fail without # failing CI. @@ -23,6 +26,11 @@ jobs: os: ubuntu-latest experimental: true fail_ci_if_error: false + # 32-bit Linux (i686), which setup-julia calls x86 + - julia-version: '1' + os: ubuntu-latest + arch: x86 + experimental: false steps: - name: Checkout Repository @@ -31,6 +39,7 @@ jobs: uses: julia-actions/setup-julia@latest with: version: ${{ matrix.julia-version }} + arch: ${{ matrix.arch || 'default' }} - name: Run Tests uses: julia-actions/julia-runtest@latest - name: Create CodeCov diff --git a/.gitignore b/.gitignore index db690aa..d9d8f9b 100644 --- a/.gitignore +++ b/.gitignore @@ -28,3 +28,5 @@ Manifest.toml TODO.md *.code-workspace + +/REVIEW.md diff --git a/AGENTS.md b/AGENTS.md new file mode 100644 index 0000000..fc2bbc7 --- /dev/null +++ b/AGENTS.md @@ -0,0 +1,14 @@ +# AGENTS.md +This file contains instructions for AI agents working with this repository. + +Purpose: PairwiseMappingFormat.jl provides parsing for the bioinformatics PAF format. + +Design goals: +* Be correct, even in edge cases. Unrealistic edge cases like 2 GiB lines or 4 GiB references may not be optimized for, but should be handled. +* Be performant + +Scope: +* Writing or mutating records is intentionally not supported yet. + +Guidelines for developing +* After making changes, format with `runic -i .`. If Runic is not installed, alert the user and do not attempt to install it. diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..5083b5c --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,75 @@ +# Changelog +All notable changes to this project will be documented in this file. + +The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), +and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). + +## v0.2.0 +### Breaking changes +- Allow XAMAuxData 0.2. Since `aux_data` returns a `XAMAuxData.SAM.Auxiliary`, + the breaking changes of XAMAuxData 0.2 apply to auxiliary data of PAF records. + For example, `H` values with lowercase hex digits are now loaded as an `Error`, + and `length` and `haskey` throw on malformed auxiliary data. + See the migration guide in the documentation of XAMAuxData.jl. +- `PAFReader` now reads from a `BufferIO.AbstractBufReader`, and its type + parameter is the type of the buffered reader. An `IO` passed to `PAFReader` is + wrapped in a `BufferIO.BufReader`, and so becomes `PAFReader{BufReader{T}}`. + Reading from an `AbstractBufReader` avoids copying the data into the reader's own buffer. +- A `PAFRecord` may now be a view into the memory it was parsed from, instead + of owning a copy of the data: + - With `copy=false`, the record iterated by `PAFReader` and returned by `try_next!` + is a view into the reader's buffer, and is invalidated when the reader is advanced. + Using an invalidated record is not detected. + - `try_parse` and `parse` return a view into the input, so mutating the input + also mutates the record. + - `copy(record)` returns a record that owns its data. +- `try_next!` now respects the `copy` keyword of `PAFReader`. With the default + `copy=true`, it returns a new record which owns its data. +- `PAFRecord` is now immutable. With `copy=false`, iterating a `PAFReader` + no longer returns the same object on each iteration, but a new record each time. + On Julia versions before 1.14, returning a `PAFRecord` from a non-inlined function, + such as `try_next!` or `iterate`, allocates. +- The `PAFRecord(buffer_size::Int)` constructor has been removed. +- `aux_data` now returns an immutable `Auxiliary{ImmutableMemoryView{UInt8}}`, + so the auxiliary data of records can no longer be modified. +- Reading from a closed `PAFReader` now throws, instead of + behaving as if the reader is empty. +- Accessing any property other than `qname` and `qlen` of an unmapped + record now throws the new exception `PairwiseMappingFormat.UnmappedError`, as does calling `query_coverage`, `target_coverage` + or `aln_identity` on it. Previously, `tname` and `is_rc` returned `nothing`, and the + remaining properties returned arbitrary values. As a result, `tname` is now always + a `StringView`, and `is_rc` always a `Bool`. +- Fields 6-12 of unmapped records are now validated like those of mapped + records, except that the integer fields may be zero. Previously, their content was ignored. +- `parse` and `try_parse` now strip one optional trailing `\n` or `\r\n`, + and reject any other `\n` with the new error kind `EmbeddedNewline`. + Previously, a trailing newline was part of the last field, and an embedded newline + was accepted in names and auxiliary data. +- In `PAFReader`, the first line consisting only of whitespace now marks the + end of the records, after which all input must be whitespace. + Previously, only an empty line did, and a line of e.g. a single space was a + `TooFewFields` error. +- When `try_next!` returns a `ParserException`, or iterating a `PAFReader` + throws one, the reader now advances past the invalid line, so the next call continues + from the following line. Previously, the reader returned the same error forever. + This allows skipping invalid lines. For an `EmptyNonTailingLine` error, the reader + skips all consecutive whitespace-only lines, and the error is reported at the first + of them, instead of at the first non-whitespace byte after them. + +### Changed +- Allow MemoryViews 0.5. +- Sequence lengths and positions up to `typemax(UInt32)` (about 4.3 Gbp) are now supported, + up from `typemax(Int32)`. Larger integers are rejected with `IntegerOverflow`. +- `parse(PAFRecord, x)` now accepts the same inputs as `try_parse`, i.e. any `x` that + implements `MemoryView`, such as `Vector{UInt8}`. Previously, it only accepted + `String` and `SubString{String}`. +- Empty query and target intervals, i.e. where the start equals the end, are now + accepted anywhere within the sequence. Previously, they were rejected with + `InvalidZero` at the start of the sequence, where the end is zero, but accepted elsewhere. + +### Fixed +- `isempty(reader)` no longer consumes a record from a `PAFReader`. It may invalidate + the most recently returned record if that record is a view into the reader's buffer, + i.e. if the reader was constructed with `copy=false`. + `PAFReader` now implements `Base.isdone`, which never invalidates records, and returns + `missing` when it can't determine whether the reader is done without doing so. diff --git a/Project.toml b/Project.toml index d64435c..f6c3cb4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,9 +1,10 @@ name = "PairwiseMappingFormat" uuid = "30ee2946-b067-41c0-bed2-0900b78dbe59" -version = "0.1.2" -authors = ["Jakob Nybo Nissen "] +version = "0.2.0" +authors = ["Jakob Nybo Andersen "] [deps] +BufferIO = "ff366017-094f-4144-bd70-a828a0d67152" MemoryViews = "a791c907-b98b-4e44-8f4d-e4c2362c6b2f" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" StringViews = "354b36f9-a18e-4713-926e-db85100087ba" @@ -11,16 +12,18 @@ XAMAuxData = "e99d641e-1821-45d7-9150-ecb7bf333fe1" [compat] Aqua = "0.8.7" +BufferIO = "0.3" FormatSpecimens = "1.3" -MemoryViews = "0.2, 0.3, 0.4" +MemoryViews = "0.2, 0.3, 0.4, 0.5" PrecompileTools = "1.2.1" -StringViews = "1.3.3" +StringViews = "2.0.0" Test = "1.11" -XAMAuxData = "0.1" +XAMAuxData = "0.2" julia = "1.11" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +BufferIO = "ff366017-094f-4144-bd70-a828a0d67152" FormatSpecimens = "3372ea36-2a1a-11e9-3eb7-996970b6ffbd" MemoryViews = "a791c907-b98b-4e44-8f4d-e4c2362c6b2f" StringViews = "354b36f9-a18e-4713-926e-db85100087ba" @@ -28,4 +31,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" XAMAuxData = "e99d641e-1821-45d7-9150-ecb7bf333fe1" [targets] -test = ["Test", "FormatSpecimens", "MemoryViews", "StringViews", "XAMAuxData", "Aqua"] +test = ["Test", "FormatSpecimens", "MemoryViews", "StringViews", "XAMAuxData", "Aqua", "BufferIO"] diff --git a/docs/make.jl b/docs/make.jl index ef6f1de..dd4989a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -14,7 +14,7 @@ makedocs(; sitename = "PairwiseMappingFormat.jl", modules = [PairwiseMappingFormat], pages = ["Home" => "index.md"], - authors = "Jakob Nybo Nissen", + authors = "Jakob Nybo Andersen", checkdocs = :public, remotes = nothing, ) diff --git a/docs/src/index.md b/docs/src/index.md index 695211d..2f93939 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -16,14 +16,15 @@ which is a simple tab-delimited format used by e.g. minimap2 and strobealign to store records representing pairwise alignment between biological sequences. ## Reader -The [`PAFReader`](@ref) type wraps an `IO`, and is an iterator of [`PAFRecord`](@ref) objects: +The [`PAFReader`](@ref) type wraps an `IO` or a `BufferIO.AbstractBufReader`, +and is an iterator of [`PAFRecord`](@ref) objects: ```jldoctest; output = false reader = PAFReader(open(path_to_paf)) records = collect(reader) +@assert isempty(reader) close(reader) -@assert isempty(reader) @assert typeof(records) == Vector{PAFRecord} # output @@ -46,23 +47,25 @@ julia> PAFReader(open(path_to_paf)) do reader ``` The [`PAFReader`](@ref) constructor takes an optional keyword `copy`, which defaults to `true`. -If it is `false`, the record will iterate the _same_ [`PAFRecord`](@ref) object, overwriting -it in-place. -This reduces allocations and gives a slight speed increase, but may result in bugs -if records, or references of records of previous iterations are stored and -unexpectedly overwritten: +If it is `false`, the [`PAFRecord`](@ref)s returned by iterating the reader, or by +[`try_next!`](@ref), are views into the reader's buffer, instead of owning a copy of their data. +Such a record is only valid until the reader is advanced, or otherwise mutated, which +may happen by other non-advancing calls such as `peek` or `isempty` on the reader. +This reduces allocations and copying and gives a slight speed increase, but may result +in bugs if records from previous iterations are stored. +Use `copy` on a record to obtain a record that owns its data: ```jldoctest -julia> records = PAFReader(collect, open(path_to_paf); copy=false); - -julia> println(map(i -> i.qlen, records)) # NB: All the same record! -[288659, 288659, 288659] +julia> PAFReader(open(path_to_paf); copy=false) do reader + [copy(record) for record in reader if record.qlen > 290000] + end |> length +2 ``` At the moments, readers do not support seeking. ## Record -The mutable [`PAFRecord`](@ref) object represents a single line in a PAF file. +The [`PAFRecord`](@ref) object represents a single line in a PAF file. The individual columns of the PAF line is obtained by accessing properties of the records: ```jldoctest record @@ -87,27 +90,31 @@ underlying fields as they are stored in the `PAFRecord`. Names of the query sequences When the sequences are from a FASTA file, this is typically the _identifier_ of the FASTA records, i.e. the part of the header before the first space. -* `tname::Union{Nothing, StringView{A}}`. Same as `qname`, but for the target (i.e. subject) - sequence. May be `nothing` if the record is unmapped. -* `qlen` and `tlen`. Of type `Int`. Length of the full query and target sequence. -* `qstart` and `tstart`. Of type `Int`. Leftmost (i.e. lowest) position of the alignment +* `tname::StringView{A}`. Same as `qname`, but for the target (i.e. subject) + sequence. +* `qlen` and `tlen`. Of type `Int64`. Length of the full query and target sequence. + Lengths and positions up to `typemax(UInt32)`, about 4.3 Gbp, are supported. +* `qstart` and `tstart`. Of type `Int64`. Leftmost (i.e. lowest) position of the alignment in the query and target, respectively. This does not take strandedness into account. -* `qend` and `tend`. Of type `Int`. Rightmost (i.e. highest) position of the alignment +* `qend` and `tend`. Of type `Int64`. Rightmost (i.e. highest) position of the alignment in the query and target, respectively. This does not take strandedness into account. - As such, we always have `record.qend - record.qtart) >= 0`, and likewise for the target. -* `alnlen::Int`. Length of alignment. That includes gaps in query and target, and may therefore + The interval may be empty, in which case `record.qstart:record.qend` is an empty range. + As such, we always have `record.qend - record.qstart >= -1`, and likewise for the target. +* `alnlen::Int64`. Length of alignment. That includes gaps in query and target, and may therefore not match either `qend - qstart + 1`, or that of the target. -* `matches::Int`. Number of residue matches (not mismatches) in the alignment. - `matches / alnlen` is the alignment identity and is between 0 and 1. +* `matches::Int64`. Number of residue matches (not mismatches) in the alignment. + `matches / alnlen` is the alignment identity, which is between 0 and 1 for well-formed + records. The parser does not check that `matches <= alnlen`, or that `alnlen > 0`. + See [`aln_identity`](@ref). * `mapq::Union{Int, Nothing}`. Mapping quality, from 0:254, where 254 is the better quality. When calibrated, the probability the mapping is correct should be 10^(-mapq / 10). A value of `nothing` indicates the mapping quality is unavailable. -* `is_rc::Union{Bool, Nothing}` indicates whether the alignment matches the target on the forward (`false`) or reverse-complement `true` strand. A missing alignment is `nothing`. +* `is_rc::Bool` indicates whether the alignment matches the target on the forward (`false`) or reverse-complement `true` strand. The PAF format does not contain fields for alignment identity, query coverage or target coverage. However, these can be approximated with the functions [`aln_identity`](@ref), -[`query_coverage`](@ref) and [`query_coverage`](@ref). +[`query_coverage`](@ref) and [`target_coverage`](@ref). Besides being iterated from [`PAFReader`](@ref), records can also be created by parsing from a bytes-like object such as a `String`: @@ -122,41 +129,34 @@ record = parse(PAFRecord, data) ``` -Rarely, PAF records are unmapped. An unmapped record contain the query id, the query length, and any [auxiliary fields](@ref aux). -The target name and the strand of an unmapped return is `nothing`. -You can check if a record is mapped or unmapped with the function `is_mapped`: +Some PAF records are unmapped, for example those written by minimap2 with the +`--paf-no-hit` flag. An unmapped record only contains the query name, the query length, +and any [auxiliary fields](@ref aux). +You can check if a record is mapped or unmapped with the function [`is_mapped`](@ref): ```jldoctest unmapped -julia> unmapped = parse(PAFRecord, "myquery\t5032251\t9\t11\t*\t*\t5\t2\t1\t7\t4\t0"); +julia> unmapped = parse(PAFRecord, "myquery\t5032251\t0\t0\t*\t*\t0\t0\t0\t0\t0\t0"); -julia> @assert !is_mapped(unmapped) +julia> is_mapped(unmapped) +false julia> println(unmapped.qname) myquery julia> println(unmapped.qlen) 5032251 - -julia> @assert isnothing(unmapped.tname) - -julia> # The strand is given by is_rc, and is nothing for unmapped records - @assert isnothing(unmapped.is_rc) ``` -All other properties of unmapped records may contain arbitrary data, and should not be relied on. -For this reason, the functions [`target_coverage`](@ref), [`query_coverage`](@ref) and [`aln_identity`](@ref) -may give nonsense answers for unmapped records: +Accessing any other property of an unmapped record throws a [`PairwiseMappingFormat.UnmappedError`](@ref), +as does calling [`target_coverage`](@ref), [`query_coverage`](@ref) or [`aln_identity`](@ref) on it: ```jldoctest unmapped -julia> # Note! These results are arbitrary and not guaranteed to be stable - # across minor releases of this package - println(target_coverage(unmapped)) -Inf - -julia> println(query_coverage(unmapped)) -1.9871822768776836e-7 +julia> unmapped.tname +ERROR: UnmappedError: Cannot get `tname` of an unmapped PAFRecord. Check `is_mapped` first. +[...] -julia> println(aln_identity(unmapped)) -NaN +julia> query_coverage(unmapped) +ERROR: UnmappedError: Cannot get `query_coverage` of an unmapped PAFRecord. Check `is_mapped` first. +[...] ``` ## Low-level interface @@ -167,7 +167,7 @@ ERROR: Error when parsing PAF record on line 1, near byte number 14 in line: Not [...] ``` -In some circuomstances, throwing exceptions may not be acceptable, and so PairwiseMappingFormat.jl +In some circumstances, throwing exceptions may not be acceptable, and so PairwiseMappingFormat.jl provide functionality for returning errors as values. These values can then be checked to programmatically handle error states. @@ -182,16 +182,14 @@ julia> println(PAF.try_parse("not a PAF line")) PairwiseMappingFormat.ParserException(PairwiseMappingFormat.Errors.TooFewFields, 14, 1) ``` -Similarly, the next record of a `PAFReader` may be obtained with the unexported [`try_next!`](@ref) +Similarly, the next record of a `PAFReader` may be obtained with the [`try_next!`](@ref) function. This function will either: -* Return `nothing` if the reader has no more records -* Return a `ParserException` if the next record is invalid -* Return a `PAFRecord` if there is another valid record. - Only in this case it will advance the position of the reader. - -In other words, if a call to [`try_next!`](@ref) returns `nothing` or a `ParserException`, -any future calls will return the same value forever: +* Return `nothing` if the reader has no more records. Any future calls will + also return `nothing`. +* Return a `ParserException` if the next record is invalid, and advance the reader + past the invalid line. Hence, invalid lines can be skipped by ignoring the error. +* Return a `PAFRecord` if there is another valid record, and advance the reader. ```jldoctest public; output = false reader = PAFReader(open(path_to_paf)) @@ -206,12 +204,13 @@ nothings = [try_next!(reader) for i in 1:10] close(reader) -reader = PAFReader(IOBuffer("not a PAF record")) +reader = PAFReader(IOBuffer("not a PAF record\n" * open(readline, path_to_paf))) err = try_next!(reader) - @assert err isa PAF.ParserException -err2 = try_next!(reader) -@assert err == err2 + +# The invalid line was skipped, so the next record can be read +@assert try_next!(reader) isa PAFRecord +@assert try_next!(reader) === nothing # output @@ -228,9 +227,16 @@ julia> err.kind TooFewFields::Err = 0 ``` -The precise value of this `.kind` object for any given error condition is subject -to change across minor versions and new values may be introduced. -This is because the same error may be detected in multiple different ways. +The `.kind` is intended for diagnostics, such as logging or counting errors, and should +not be relied on to handle errors programmatically. +An invalid line is often invalid in several ways, and which of them is reported depends +on implementation details. Therefore, the kind reported for a given invalid input is +deterministic within a version of this package, but may change in any minor release, +and new kinds may be added. +In contrast, whether input is rejected, and the line the error is reported on, are stable. + +The error message also includes the approximate position of the error within the line. +This position is only meant for humans to read, and is not part of the API. ## [Auxiliary fields](@id aux) [`PAFRecord`](@ref)s may contain extra auxiliary fields at the end of the records, @@ -266,19 +272,15 @@ julia> aux["km"] julia> haskey(aux, "cm") true - -julia> aux["AB"] = 5.5; - -julia> haskey(aux, "AB") -true ``` ## Misc info -* The `PAFReader` cannot handle trailing whitespace, except trailing newlines at - the end of the file. This is because trailing whitespace may be significant - in some records, e.g. if it ends with an auxiliary field of the element type - `Z`. - A PAF record with trailing whitespace is considered invalid. +* Lines end with `\n` or `\r\n`. Other whitespace at the end of a record is not + stripped, because it may be significant, e.g. if the record ends with an auxiliary + field of the element type `Z`. Instead, it is part of the last field, and may make it + invalid. +* The first line in a file consisting only of whitespace marks the end of the records. + All remaining input must then be whitespace. ## Reference @@ -291,8 +293,12 @@ aln_identity query_coverage target_coverage try_next! +Base.isempty(::PAFReader) +Base.peek(::PAFReader) is_mapped PairwiseMappingFormat.Errors PairwiseMappingFormat.try_parse +Base.parse(::Type{PAFRecord}, ::Any) PairwiseMappingFormat.ParserException +PairwiseMappingFormat.UnmappedError ``` diff --git a/src/PairwiseMappingFormat.jl b/src/PairwiseMappingFormat.jl index b2980d5..feeaf78 100644 --- a/src/PairwiseMappingFormat.jl +++ b/src/PairwiseMappingFormat.jl @@ -1,7 +1,7 @@ """ module PairwiseMappingFormat -Parse files of the PairwisemAppingFormat (PAF) format. +Parse files in the Pairwise mApping Format (PAF). """ module PairwiseMappingFormat @@ -9,7 +9,7 @@ export PAFReader, PAFRecord, aln_identity, aux_data, try_next!, query_coverage, target_coverage, is_mapped # Not exported, because this very PAF-specific, and operates on foreign types. -public try_parse, ParserException, Errors, Err +public try_parse, ParserException, UnmappedError, Errors, Err using PrecompileTools: @compile_workload, @setup_workload, @recompile_invalidations @@ -22,16 +22,21 @@ using PrecompileTools: @compile_workload, @setup_workload, @recompile_invalidati using StringViews: StringView end -using MemoryViews: MemoryView, ImmutableMemoryView +using MemoryViews: MemoryViews, ImmutableMemoryView using XAMAuxData.SAM: Auxiliary +using BufferIO: AbstractBufReader, BufReader, ExpandReader, get_buffer, fill_buffer, consume """ - PAFRecord(buffer_size::Int) + PAFRecord -Mutable type representing a single PAF line. The content of the record +Immutable type representing a single PAF line. The content of the record is accessed through its properties. +A `PAFRecord` does not necessarily own its data, but may be a view into the +memory it was parsed from, e.g. the buffer of a [`PAFReader`](@ref) with `copy=false`. +Use `copy` to obtain a record that owns its data. + # Examples ```jldoctest julia> record = PAFReader(first, open(path_to_paf)); @@ -41,65 +46,93 @@ julia> record.qname julia> record.qlen 301156 + +julia> startswith(String(record), "query1\t301156\t") +true ``` See also: [`PAFReader`](@ref) # Extended help -The following properties may be used: -* `qname::StringView`. The query name, May be empty, and contain any bytes +The following properties may be used. For unmapped records, only `qname` and `qlen` +are available, and accessing any other property throws an [`UnmappedError`](@ref). +Use [`is_mapped`](@ref) to check if a record is mapped. +* `qname::StringView`. The query name. May be empty, and contain any bytes except `\\t` and `\\n`. -* `tname::Union{StringView, Nothing}`. Target name. Like `qname`, but is `nothing` - if and only if the record is unmapped. -* `qlen` and `tlen::Int`, and gives the length of the query and target - sequences, respectively. This must be > 0. -* `qstart`, `qend`, `tstart` and `tend::Int`, and give the starting and +* `tname::StringView`. Target name. Like `qname`. +* `qlen` and `tlen::Int64` give the length of the query and target + sequences, respectively. This must be > 0, and at most `typemax(UInt32)` (about 4.3 Gbp). +* `qstart`, `qend`, `tstart` and `tend::Int64` give the starting and ending positions of the alignments on the query and target strand. - These uses one-based, closed interval indexing as usual in Julia, but unlike - the underlying PAF format. - The ending positions are always >= the starting ones, and the ending positions - are always <= the query/target lengths. -* `matches::Int` gives the number of nucleotides / residues which match + These use one-based, closed interval indexing as usual in Julia, but unlike + the underlying PAF format. The intervals `qstart:qend` or `tstart:tend` may be empty. +* `matches::Int64` gives the number of nucleotides / residues which match (i.e. are equal) in the alignment. -* `alnlen::Int` gives the length of the alignment -* `mapq::Union{Int, Nothing}` gives the mapping quality, or `nothing` if +* `alnlen::Int64` gives the length of the alignment +* `mapq::Union{Int, Nothing}` gives the mapping quality, or `nothing` if this information is unavailable. -* `is_rc::Union{Bool,Nothing}` tells if the query and target strands are reverse-complement - relative to each other. Is `nothing` if the record is unmapped. -""" -mutable struct PAFRecord - # Data contains query name, subject name, AUX data, - # without delimiters and with no trailing whitespace. - const data::Vector{UInt8} - qname_len::Int32 - tname_len::Int32 - qlen::Int32 # always >= 1 - # Note: qstart/qend (and tstart/tend) are zero-indexed semi-open interval - # in PAF format, but stored in this struct one-based. - qstart::Int32 # always >= 1, except when unmapped - qend::Int32 # always >= qstart - tlen::Int32 # always >= 1, except when unmapped - tstart::Int32 # always >= 1, except when unmapped - tend::Int32 # always >= tstart - matches::Int32 - alnlen::Int32 # always >= 1, except when unmapped +* `is_rc::Bool` tells if the query and target strands are reverse-complement + relative to each other. + +The line the record was parsed from, without the trailing newline, can be obtained +with `MemoryView(record)`, which returns an `ImmutableMemoryView{UInt8}` view of the +record's data, or `String(record)`, which copies it into a new `String`. +Unlike `MemoryView`, `String` and `string` do not alias the record's memory, so they +remain valid if the underlying memory is mutated. +For a non-copying string, use `StringView(MemoryView(record))`. +""" +struct PAFRecord + # The line the record was parsed from, without the trailing newline. + # May be a view into memory not owned by the record, e.g. a reader's buffer. + data::ImmutableMemoryView{UInt8} + qname_len::Int32 # query name begins at index 1 + tname_start::Int32 + tname_len::Int32 # always zero when unmapped + aux_start::Int32 # lastindex(data) + 1 if there is no aux data + # Positions and lengths are stored as UInt32, so sequences of up to + # typemax(UInt32) bytes are supported. + qlen::UInt32 # always >= 1 + # Note: qstart/qend (and tstart/tend) are zero-indexed semi-open intervals + # in the PAF format, and are stored as such. The properties are one-based. + # The start is stored zero-based, such that the conversion can't overflow. + # The fields below, and tname, are arbitrary when unmapped, and not accessible + # through the properties of the record. + qstart::UInt32 + qend::UInt32 # always >= qstart + tlen::UInt32 # always >= 1 + tstart::UInt32 + tend::UInt32 # always >= tstart + matches::UInt32 + alnlen::UInt32 # may be zero, and may be smaller than matches mapq::UInt8 # stored as 0xff for missing - strand::UInt8 # 0x00: *, 0x01: -, 0x02: + + strand::UInt8 # 0x00: *, 0x01: +, 0x02: - # 16 bits of padding, which we can use for something else later end # Since the interface of the records are defined in terms of its -# properties, we need to be able to abstract over them -function Base.getproperty(record::PAFRecord, sym::Symbol) +# properties, we need to be able to abstract over them. +# Positions and lengths are Int64, not Int, since they may not fit in a 32-bit Int. +# Note: This is inlined so `sym` is constant folded, leaving only the branch for +# that property, and the cheap check for mapped-ness. +@inline function Base.getproperty(record::PAFRecord, sym::Symbol) return if sym === :qname qname(record) - elseif sym == :tname + elseif sym === :qlen + Int64(getfield(record, :qlen)) + elseif sym === :tname + check_mapped(record, sym) tname(record) - elseif sym == :mapq + elseif sym === :mapq + check_mapped(record, sym) mapq(record) - elseif sym in (:qlen, :qstart, :qend, :tlen, :tstart, :tend, :matches, :alnlen) - Int(getfield(record, sym)) - elseif sym == :is_rc + elseif sym in (:qstart, :tstart) + check_mapped(record, sym) + Int64(getfield(record, sym)) + 1 + elseif sym in (:qend, :tlen, :tend, :matches, :alnlen) + check_mapped(record, sym) + Int64(getfield(record, sym)) + elseif sym === :is_rc + check_mapped(record, sym) is_rc(record) else error(lazy"Type PAFRecord has no property $sym") @@ -107,30 +140,70 @@ function Base.getproperty(record::PAFRecord, sym::Symbol) end """ - is_mapped(x::PAFRecord) -> Bool + is_mapped(x::PAFRecord)::Bool -Compute whether the `PAFRecord` is mapped. An unmapped record will have -the properties `is_rc` and `tname` unavailable. -The properties `qname` and `qlen`, and the auxiliary data of an unmapped record -can be relied on, but the remaining properties contain arbitrary data. +Compute whether the `PAFRecord` is mapped. +Only the properties `qname` and `qlen`, and the auxiliary data of an unmapped record +are available. Accessing any other property of an unmapped record, or calling +[`query_coverage`](@ref), [`target_coverage`](@ref) or [`aln_identity`](@ref) on it, +throws an [`UnmappedError`](@ref). Note that some PAF parsers do not handle unmapped records correctly, so be wary when writing unmapped records. # Examples -```jldoctes +```jldoctest julia> is_mapped(record) true + +julia> unmapped = parse(PAFRecord, "read1\\t1523\\t0\\t0\\t*\\t*\\t0\\t0\\t0\\t0\\t0\\t0"); + +julia> is_mapped(unmapped) +false ``` See also: [`PAFRecord`](@ref) """ is_mapped(x::PAFRecord) = !iszero(getfield(x, :strand)) -function is_rc(record::PAFRecord)::Union{Bool, Nothing} - x = getfield(record, :strand) - return iszero(x) ? nothing : x == 0x02 +""" + UnmappedError + +The exception thrown when accessing a property of an unmapped [`PAFRecord`](@ref) which +only exists for mapped records, such as `tname`, or when calling a function which +requires a mapped record, such as [`query_coverage`](@ref). +The name of the property or function may be accessed with the `.name` field. + +# Examples +```jldoctest +julia> unmapped = parse(PAFRecord, "read1\\t1523\\t0\\t0\\t*\\t*\\t0\\t0\\t0\\t0\\t0\\t0"); + +julia> err = try unmapped.tname catch e e end; + +julia> err isa PAF.UnmappedError +true + +julia> err.name +:tname +``` + +See also: [`is_mapped`](@ref) +""" +struct UnmappedError <: Exception + name::Symbol +end + +function Base.showerror(io::IO, err::UnmappedError) + return print(io, "UnmappedError: Cannot get `", err.name, "` of an unmapped PAFRecord. Check `is_mapped` first.") end +# Keep the error path out of line, so the check is cheap to inline +@noinline throw_unmapped(name::Symbol) = throw(UnmappedError(name)) + +check_mapped(record::PAFRecord, name::Symbol) = (is_mapped(record) || throw_unmapped(name); nothing) + +# Note: The functions below do not check if the record is mapped. +is_rc(record::PAFRecord) = getfield(record, :strand) == 0x02 + # For tab completion function Base.propertynames(::PAFRecord) return ( @@ -149,11 +222,6 @@ function Base.propertynames(::PAFRecord) ) end -function PAFRecord(size::Int = 0) - data = Vector{UInt8}(undef, max(size, 0)) - return PAFRecord(data, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0xff, 0x00) -end - function Base.show(io::IO, ::MIME"text/plain", record::PAFRecord) buf = IOBuffer() is_mapped(record) || print(buf, "Unmapped ") @@ -187,20 +255,29 @@ end # Return StringView to not allocate function qname(record::PAFRecord) - return StringView( - ImmutableMemoryView(getfield(record, :data))[1:(getfield(record, :qname_len))], - ) + return StringView(getfield(record, :data)[1:(getfield(record, :qname_len))]) end # Return StringView to not allocate -function tname(record::PAFRecord)::Union{StringView, Nothing} - return if is_mapped(record) - ql = getfield(record, :qname_len) - span = (ql + 1):(ql + getfield(record, :tname_len)) - StringView(ImmutableMemoryView(getfield(record, :data))[span]) - else - nothing - end +function tname(record::PAFRecord) + start = getfield(record, :tname_start) + span = start:(start + getfield(record, :tname_len) - 1) + return StringView(getfield(record, :data)[span]) +end + +# The original line, without the trailing newline +MemoryViews.MemoryView(record::PAFRecord) = getfield(record, :data) +Base.String(record::PAFRecord) = String(getfield(record, :data)) +Base.print(io::IO, record::PAFRecord) = (write(io, getfield(record, :data)); nothing) + +# Records are equal if their lines are byte-equal, regardless of whether they share memory. +Base.:(==)(a::PAFRecord, b::PAFRecord) = getfield(a, :data) == getfield(b, :data) + +# Salt the hash, such that a record does not hash like its line as a string. +# Hash through StringView, because it is faster than the generic array hash. +const RECORD_HASH_SALT = 0x8a2e5d1f6b3c9047 % UInt +function Base.hash(record::PAFRecord, h::UInt) + return hash(StringView(getfield(record, :data)), h ⊻ RECORD_HASH_SALT) end # Note: We return nothing instead of missing because fuck missing @@ -210,11 +287,14 @@ function mapq(record::PAFRecord)::Union{Int, Nothing} end """ - aln_identity(rec::Record) -> Float64 + aln_identity(record::PAFRecord)::Float64 Return the nucleotide / residue identity, defined as the number of matches divided by the alignment length. -This is a number between 0 and 1. +For well-formed records, this is a number between 0 and 1. However, the parser does not +check that the number of matches is at most the alignment length, so this may be larger +than 1. If the alignment length is zero, this returns `NaN` or `Inf`. +Throws an [`UnmappedError`](@ref) if the record is unmapped. # Examples ```jldoctest @@ -225,13 +305,15 @@ julia> aln_identity(record) See also: [`query_coverage`](@ref), [`target_coverage`](@ref) """ function aln_identity(record::PAFRecord)::Float64 + check_mapped(record, :aln_identity) return record.matches / record.alnlen end """ - query_coverage(rec::Record) -> Float64 + query_coverage(record::PAFRecord)::Float64 Compute the fraction of the query covered by the alignment. +Throws an [`UnmappedError`](@ref) if the record is unmapped. # Examples ```jldoctest @@ -242,13 +324,15 @@ julia> query_coverage(record) See also: [`target_coverage`](@ref), [`aln_identity`](@ref) """ function query_coverage(record::PAFRecord)::Float64 - return (record.qend - record.qstart + Int32(1)) / record.qlen + check_mapped(record, :query_coverage) + return (record.qend - record.qstart + 1) / record.qlen end """ - target_coverage(rec::Record) -> Float64 + target_coverage(record::PAFRecord)::Float64 Compute the fraction of the target covered by the alignment. +Throws an [`UnmappedError`](@ref) if the record is unmapped. # Examples ```jldoctest @@ -259,11 +343,12 @@ julia> target_coverage(record) See also: [`query_coverage`](@ref), [`aln_identity`](@ref) """ function target_coverage(record::PAFRecord)::Float64 - return (record.tend - record.tstart + Int32(1)) / record.tlen + check_mapped(record, :target_coverage) + return (record.tend - record.tstart + 1) / record.tlen end """ - aux_data(rec::Record) -> XAMAuxData.SAM.Auxiliary + aux_data(record::PAFRecord)::XAMAuxData.SAM.Auxiliary Return a lazily evaluated and validated `SAM.Auxiliary`, which is an `AbstractDict` of the auxiliary fields at the end of the record. @@ -280,17 +365,12 @@ julia> length(aux) julia> aux["cm"] 29990 -julia> aux["k1"] = "add a new aux string like this"; - julia> haskey(aux, "k1") -true +false ``` """ function aux_data(record::PAFRecord) - return Auxiliary( - getfield(record, :data), - getfield(record, :qname_len) + getfield(record, :tname_len) + 1, - ) + return Auxiliary(getfield(record, :data)[getfield(record, :aux_start):end]) end """ @@ -299,6 +379,10 @@ end Wrapper module used to contain the instances of error types, and not pollute the namespace of the package. +The error kinds are intended for diagnostics. Which kind is reported for a given invalid +input may change in any minor release, and new kinds may be added. +See [`PairwiseMappingFormat.ParserException`](@ref). + # Examples ```jldoctest julia> print(PAF.Errors.InvalidZero) @@ -320,6 +404,8 @@ module Errors EmptyInteger InvalidZero PositionOutOfBounds + LineTooLong + EmbeddedNewline end end # module Errors @@ -331,8 +417,20 @@ using .Errors: Err The exception type thrown by iterating [`PAFReader`](@ref), or a failed `parse` of a [`PAFRecord`](@ref). The functions `try_parse` and `try_next!` may return (not throw) values of this type. -The line the error occurred may be accessed with the `.line` field. -The error kind may be accessed with the `.kind` field. + +The line the error occurred on may be accessed with the `.line` property. +For `try_parse` and `parse`, this is always 1. + +The `.kind` property is a [`PairwiseMappingFormat.Errors.Err`](@ref PairwiseMappingFormat.Errors) +describing why the line was rejected. It is intended for diagnostics, such as logging or +counting errors: An invalid line is often invalid in several ways, and which of them is +reported depends on implementation details. Therefore, the kind reported for a given +invalid input is deterministic within a version of this package, but may change in any +minor release, and new kinds may be added. +In contrast, whether input is rejected, and the line the error is reported on, are stable. + +The error message printed by `showerror` includes the approximate position of the error +within the line. This position is not part of the API. # Examples ```jldoctest @@ -348,7 +446,8 @@ TooFewFields struct ParserException <: Exception kind::Errors.Err index_in_line::Int32 - line::Int + # Int64, not Int, such that it does not overflow for huge files on 32-bit platforms + line::Int64 end ParserException(index::Integer, kind::Errors.Err) = ParserException(kind, index, 1) @@ -378,7 +477,7 @@ function Base.showerror(io::IO, err::ParserException) elseif kind == Errors.BadInteger "Integer field contains characters out of 0-9 range" elseif kind == Errors.InvalidStrand - "The strand column must be '+' or '-'" + "The strand column must be '+' or '-', or '*' for unmapped records" elseif kind == Errors.EmptyNonTailingLine "Whitespace-only line(s) followed by non-whitespace" elseif kind == Errors.BackwardsIndices @@ -389,6 +488,10 @@ function Base.showerror(io::IO, err::ParserException) "Parsed an integer to zero which cannot meaningfully be zero" elseif kind == Errors.PositionOutOfBounds "Query or target start or end position is longer than query/target length" + elseif kind == Errors.LineTooLong + "Line is too long. PairwiseMappingFormat.jl can't handle lines of 2147483647 bytes or longer" + elseif kind == Errors.EmbeddedNewline + "Line contains a newline. Only a single trailing \"\\n\" or \"\\r\\n\" is allowed" end print(buf, s) return print(io, String(take!(buf))) @@ -404,11 +507,21 @@ macro var"?"(expr) end """ - try_parse(x) -> Union{PAFRecord, ParserException} + try_parse(x)::Union{PAFRecord, ParserException} Try parsing `x` into a `PAFRecord`. `x` may be any type that implements `MemoryView`, such as `String` or `Vector{UInt8}`. +`x` must contain a single line, optionally ending with `"\\n"` or `"\\r\\n"`. +Any other `\\n` in `x` gives an error of kind `EmbeddedNewline`. +Other trailing whitespace, including a lone `\\r`, is not stripped: +It is part of the last field, and may make it invalid. + +!!! warning + The resulting record is a view into the memory of `x`, and does not copy it. + If `x` is mutated, the record is mutated too, and may become invalid. + Use `copy` on the record to obtain a record that owns its data. + # Examples ```jldoctest julia> valid_line = open(readline, path_to_paf); @@ -420,41 +533,83 @@ julia> typeof(PAF.try_parse("invalid string")) PairwiseMappingFormat.ParserException ``` -See also: [`PAFRecord`](@ref), [`try_next!`](@ref) +See also: [`PAFRecord`](@ref), [`try_next!`](@ref), [`parse`](@ref Base.parse(::Type{PAFRecord}, ::Any)) """ try_parse(x) = try_parse(ImmutableMemoryView(x)) -# At least 19 bytes of input are stored in the fixed fields, so -# we need at most 19 fewer bytes in the data field of the vector. -try_parse(x::ImmutableMemoryView{UInt8}) = parse_line!(PAFRecord(length(x) - 19), x) -function Base.parse(::Type{PAFRecord}, s::Union{String, SubString{String}}) - y = try_parse(s) +# Strip one optional trailing \n or \r\n +function strip_newline(x::ImmutableMemoryView{UInt8})::ImmutableMemoryView{UInt8} + (isempty(x) || @inbounds(x[end]) != UInt8('\n')) && return x + has_cr = length(x) > 1 && @inbounds(x[end - 1]) == UInt8('\r') + return @inbounds x[1:(end - 1 - has_cr)] +end + +function try_parse(x::ImmutableMemoryView{UInt8}) + x = strip_newline(x) + # Check the length before searching for newlines, since error positions are Int32 + @? check_line_length(x) + newline = findfirst(==(UInt8('\n')), x) + newline === nothing || return ParserException(newline, Errors.EmbeddedNewline) + return parse_line(x) +end + +""" + parse(::Type{PAFRecord}, x)::PAFRecord + +Parse `x` into a `PAFRecord`, or throw a [`ParserException`](@ref +PairwiseMappingFormat.ParserException) if `x` is not a valid PAF line. +`x` may be any type that implements `MemoryView`, such as `String` or `Vector{UInt8}`. +See [`try_parse`](@ref PairwiseMappingFormat.try_parse) for what counts as a line. + +!!! warning + The resulting record is a view into the memory of `x`, and does not copy it. + If `x` is mutated, the record is mutated too, and may become invalid. + Use `copy` on the record to obtain a record that owns its data. + +# Examples +```jldoctest +julia> record = parse(PAFRecord, codeunits("q\\t10\\t0\\t5\\t+\\tt\\t20\\t5\\t10\\t5\\t5\\t60\\n")); + +julia> (record.qname, record.tname, record.mapq) +("q", "t", 60) +``` +""" +function Base.parse(::Type{PAFRecord}, x) + y = try_parse(x) return y isa PAFRecord ? y : throw(y) end -function parse_line!( - record::PAFRecord, - mem::ImmutableMemoryView{UInt8}, - )::Union{PAFRecord, ParserException} - if lastindex(mem) > typemax(Int32) - error("PairwiseMappingFormat.jl can't handle lines longer than 2147483647 bytes") +# Indices into the line, and lastindex(line) + 1, where the aux data of a record may begin, +# are stored as Int32. Every index can be converted to Int32 if this returns `nothing`. +function check_line_length(mem::ImmutableMemoryView{UInt8})::Union{Nothing, ParserException} + return if lastindex(mem) >= typemax(Int32) + ParserException(typemax(Int32), Errors.LineTooLong) + else + nothing end +end + +function parse_line(mem::ImmutableMemoryView{UInt8})::Union{PAFRecord, ParserException} + # All the indices converted to Int32 by `% Int32` below can't be truncated, because + # of this check. + @? check_line_length(mem) # Read fields incrementally. `i` is the index to read from, beginning at 1. # The parse_x functions will return the next index to begin from (after the tab) (qname, i) = @? parse_str(mem, 1) (qlen, i) = @? parse_int(mem, i, false) + # Intervals may be empty, i.e. start == end, so the end may be zero (qstart, i) = @? parse_int(mem, i, true) - (qend, i) = @? parse_int(mem, i, false) + (qend, i) = @? parse_int(mem, i, true) qend ≥ qstart || return ParserException(i % Int32, Errors.BackwardsIndices) qend > qlen && return ParserException(i % Int32, Errors.PositionOutOfBounds) - # Load the strand field. We need at least 14 more bytes to encode 7 more mandatory fields - # plus 7 more tabs. - i + 14 > lastindex(mem) && + # Load the strand field. The rest of the line is at least 14 bytes: The strand, the + # tab after it, an empty tname, and six single-digit integers separated by six tabs. + lastindex(mem) - i < 13 && return ParserException(lastindex(mem) % Int32, Errors.TooFewFields) b = @inbounds mem[i] strand = if b == UInt8('*') - return finish_unmapped!(record, mem, qname, qlen, i + 2) + 0x00 elseif b == UInt8('-') 0x02 elseif b == UInt8('+') @@ -462,14 +617,17 @@ function parse_line!( else return ParserException(i % Int32, Errors.InvalidStrand) end + # The byte after the strand is in bounds by the length check above. If it is not a tab, + # the strand field is longer than one byte, like "+-". @inbounds mem[i + 1] == UInt8('\t') || - return ParserException((i + 1) % Int32, Errors.TooFewFields) + return ParserException((i + 1) % Int32, Errors.InvalidStrand) + iszero(strand) && return finish_unmapped(mem, qname, qlen, i + 2) # Load rest of the fields (tname, i) = @? parse_str(mem, i + 2) # i + 2 because strand plus tab took two bytes (tlen, i) = @? parse_int(mem, i, false) (tstart, i) = @? parse_int(mem, i, true) - (tend, i) = @? parse_int(mem, i, false) + (tend, i) = @? parse_int(mem, i, true) tend ≥ tstart || return ParserException(i % Int32, Errors.BackwardsIndices) tend > tlen && return ParserException(i % Int32, Errors.PositionOutOfBounds) (matches, i) = @? parse_int(mem, i, true) @@ -483,85 +641,49 @@ function parse_line!( mapq % UInt8 end - # Copy data to the .data field - query name, target name, AUX data - auxlen = lastindex(mem) - i + 1 - filled = length(qname) + length(tname) + auxlen - data = getfield(record, :data) - length(data) == filled || resize!(data, filled) - doff = 1 - dataview = MemoryView(data) - unsafe_copyto!(dataview, mem[qname]) - doff += length(qname) - unsafe_copyto!(dataview[doff:end], mem[tname]) - doff += length(tname) - iszero(auxlen) || unsafe_copyto!(dataview[doff:end], mem[i:end]) - - # Fill in fields of the struct - # Note: PAF uses zero-based semi-open indexing like e.g. Python, - # so compensate for that in qstart and qend. - record.qname_len = length(qname) % Int32 - record.tname_len = length(tname) % Int32 - record.qlen = qlen - record.qstart = qstart + one(qstart) - record.qend = qend - record.tlen = tlen - record.tstart = tstart + one(tstart) - record.tend = tend - record.matches = matches - record.alnlen = alnlen - record.mapq = mapq - record.strand = strand - - return record -end - -# If the record is unmapped we don't even bother reading -# the rest of the fields, so take this fast path + # The aux data begins at `i`. + return PAFRecord( + mem, + length(qname) % Int32, + first(tname) % Int32, + length(tname) % Int32, + i % Int32, + qlen, + qstart, + qend, + tlen, + tstart, + tend, + matches, + alnlen, + mapq, + strand, + ) +end + +# Parse the remaining fields of an unmapped record. +# These are validated like for mapped records, except that the integers may be zero, +# as written by e.g. minimap2's `--paf-no-hit`. Their values are not stored. # i is the byte index of the strand plus two -function finish_unmapped!( - record::PAFRecord, +function finish_unmapped( mem::ImmutableMemoryView{UInt8}, qname::UnitRange{Int}, - qlen::Int32, + qlen::UInt32, i::Int, - ) - # Determine the position of aux data - for _ in 1:6 - i = findnext(==(UInt8('\t')), mem, i) - i = if isnothing(i) - return ParserException(lastindex(mem) % Int32, Errors.TooFewFields) - else - i + 1 - end - end - i = findnext(==(UInt8('\t')), mem, i) - aux_end = lastindex(mem) - aux_start = isnothing(i) ? aux_end + 1 : i + 1 - aux = aux_start:aux_end - - # Copy data - filled = length(qname) + length(aux) - data = getfield(record, :data) - dataview = MemoryView(data) - length(data) == filled || resize!(data, filled) - unsafe_copyto!(dataview, mem[qname]) - unsafe_copyto!(dataview[(length(qname) + 1):end], mem[aux]) - - # Fill in fields - record.qname_len = length(qname) % Int32 - record.tname_len = 0 - record.qlen = qlen - record.qstart = 0 - record.qend = 0 - record.tlen = 0 - record.tstart = 0 - record.tend = 0 - record.matches = 0 - record.alnlen = 0 - record.mapq = 0xff - record.strand = 0x00 - - return record + )::Union{PAFRecord, ParserException} + (_, i) = @? parse_str(mem, i) + (tlen, i) = @? parse_int(mem, i, true) + (tstart, i) = @? parse_int(mem, i, true) + (tend, i) = @? parse_int(mem, i, true) + tend ≥ tstart || return ParserException(i % Int32, Errors.BackwardsIndices) + tend > tlen && return ParserException(i % Int32, Errors.PositionOutOfBounds) + (_, i) = @? parse_int(mem, i, true) # matches + (_, i) = @? parse_int(mem, i, true) # alnlen + (mapq, i) = @? parse_int(mem, i, true, true) + mapq > 255 && return ParserException(i % Int32, Errors.IntegerOverflow) + + # The aux data begins at `i`. + return PAFRecord(mem, length(qname) % Int32, 1, 0, i % Int32, qlen, 0, 0, 0, 0, 0, 0, 0, 0xff, 0x00) end # Just get the index of the next \t @@ -577,27 +699,74 @@ function parse_str( end end -# TODO: 40% of runtime is spent in this function. Microoptimize it +# Return the number of leading bytes of `x` (in memory order) which are ASCII digits, +# from 0 to 8, where `x` is 8 bytes loaded as little-endian. +function n_leading_digits(x::UInt64)::Int + # A byte b is a digit iff the high nibbles of both b and b + 6 are 3. + # Adding 6 may carry into the next byte if b >= 0xfa, but that only corrupts bytes + # after a non-digit, which are not counted. + t = (x & 0xf0f0f0f0f0f0f0f0) | (((x + 0x0606060606060606) & 0xf0f0f0f0f0f0f0f0) >> 4) + return trailing_zeros(t ⊻ 0x3333333333333333) >>> 3 +end + +# Compute the value of the first `nd` digits of `x`, with `nd` in 0:8, and `x` +# as in `n_leading_digits`. +function swar_digits(x::UInt64, nd::Int)::UInt64 + # Subtracting may borrow from bytes after the digits, which are shifted out. + # Shifting moves the digits to the most significant bytes, with zeros in front. + x = (x - 0x3030303030303030) << ((8 - nd) << 3) + x = ((x * 10) + (x >> 8)) & 0x00ff00ff00ff00ff + x = ((x * 100) + (x >> 16)) & 0x0000ffff0000ffff + return ((x * 10000) + (x >> 32)) & 0x00000000ffffffff +end + +# Load 8 bytes beginning at index i as little-endian. Index i + 7 must be in bounds. +function load_u64(v::ImmutableMemoryView{UInt8}, i::Int)::UInt64 + return GC.@preserve v ltoh(unsafe_load(Ptr{UInt64}(pointer(v, i)))) +end + function parse_int( v::ImmutableMemoryView{UInt8}, from::Int, allow_zero::Bool, at_end::Bool = false, - )::Union{Tuple{Int32, Int}, ParserException} - n = Int32(0) - i = from - for outer i in from:lastindex(v) + )::Union{Tuple{UInt32, Int}, ParserException} + # Fast path: Parse integers of up to 7 digits from 8 bytes at once, avoiding a + # data-dependent loop. Longer integers, i.e. >= 10 Mbp, continue in the loop below. + # A similar second fast path for 8-10 digits was measured to be slower. + # Accumulate in an UInt64, not Int, which can overflow on 32-bit platforms + n = UInt64(0) + start = from + # Note: Written so it can't overflow, unlike `from + 7 <= lastindex(v)`. + if lastindex(v) - from >= 7 + x = load_u64(v, from) + nd = n_leading_digits(x) + n = swar_digits(x, nd) + if nd < 8 + # The byte after the digits is in bounds, since there are 8 bytes + i = from + nd + @inbounds(v[i]) == UInt8('\t') || return ParserException(i % Int32, Errors.BadInteger) + iszero(nd) && return ParserException(i % Int32, Errors.EmptyInteger) + (!allow_zero & iszero(n)) && return ParserException(i % Int32, Errors.InvalidZero) + # At most 7 digits, so this fits in an UInt32 + return (n % UInt32, i + 1) + end + start = from + 8 + end + # Slow path. `n` can't overflow, because we check against typemax(UInt32) after + # every digit. The result therefore fits in an UInt32. + for i in start:lastindex(v) b = @inbounds v[i] # If we hit a tab, this field is over. if b == UInt8('\t') i == from && return ParserException(i % Int32, Errors.EmptyInteger) (!allow_zero & iszero(n)) && return ParserException(i % Int32, Errors.InvalidZero) - return (n, i + 1) + return (n % UInt32, i + 1) end b ∈ 0x30:0x39 || return ParserException(i % Int32, Errors.BadInteger) - n < 0 && return ParserException(i % Int32, Errors.IntegerOverflow) - n = Int32(10) * n + (b - Int32(48)) + n = 10 * n + (b - 0x30) + n > typemax(UInt32) && return ParserException(i % Int32, Errors.IntegerOverflow) end # Note: It is not possible to get an InvalidZero error here, because # the last field (mapq) does allow zero, and for any other fields, @@ -607,28 +776,42 @@ function parse_int( return if !at_end ParserException(lastindex(v) % Int32, Errors.TooFewFields) elseif from > lastindex(v) - ParserException(i % Int32, Errors.EmptyInteger) + ParserException(from % Int32, Errors.EmptyInteger) else - (n, i + 1) + (n % UInt32, lastindex(v) + 1) end end +# The line of the record returned by `PAFRecord()`: An unmapped record with an empty +# query name, as written by minimap2's `--paf-no-hit`. +# All empty records alias this memory, which is fine, since it is never mutated. +const EMPTY_PAF_CONTENT = ImmutableMemoryView(Memory{UInt8}(codeunits("\t1\t0\t0\t*\t*\t0\t0\t0\t0\t0\t0"))) + +# Parsed once, such that its fields are always consistent with its line. +# This must be defined after the parsing functions, since they are called here. +const EMPTY_RECORD = parse(PAFRecord, EMPTY_PAF_CONTENT) + +PAFRecord() = EMPTY_RECORD + """ - PAFReader(io::IO; copy::Bool=true) + PAFReader(io::Union{IO, AbstractBufReader}; copy::Bool=true) Construct a `PAFReader`, an iterator of [`PAFRecord`](@ref) read from `io`. Iterating this object returns a `PAFRecord` for each valid line of input, and throws a [`ParserException`](@ref) when an invalid record is found. -For efficiency, the auxiliary fields are not validated until they are accessed. -If `copy` is false, the *same* record will be returned on each iteration, -with its content being overwritten. This removes a few allocations per iteration, -but may cause problems if records of old iterations are stored. +If `copy` is false, the records returned by iterating the reader, or by [`try_next!`](@ref), +are views into the reader's buffer, and are only valid until the reader is mutated, +which happens when iterating it, or calling `try_next!`, `peek` or `isempty` on it. +This removes a few allocations and copies per iteration, but records must not be used +after they are invalidated. + +See also: [`PAFRecord`](@ref), [`try_next!`](@ref) # Examples ```jldoctest julia> reader = PAFReader(open(path_to_paf)); typeof(reader) -PAFReader{IOStream} +PAFReader{BufferIO.BufReader{IOStream}} julia> typeof(first(reader)) PAFRecord @@ -643,41 +826,52 @@ julia> PAFReader(open(path_to_paf)) do reader 288659 julia> PAFReader(open(path_to_paf); copy=false) do reader - fst = first(reader) - all(reader) do record - record === fst - end + sum(record -> record.qlen, reader) end -true +889088 ``` -See also: [`PAFRecord`](@ref), [`try_next!`](@ref) +# Extended help +After throwing a `ParserException`, the reader advances over the malformed record, +so iteration may be consumed after catching the exception. + +For efficiency, the auxiliary fields are not validated until they are accessed. + +Lines end with `"\\n"` or `"\\r\\n"`. Other whitespace at the end of a line is not stripped: +It is part of the last field, and may make it invalid. + """ -mutable struct PAFReader{I <: IO} - const io::I - const record::PAFRecord - mem::Memory{UInt8} - # Index of first unused byte in buffer. Equal to filled + 1 - # if the buffer has no usable data. - index::Int32 - # mem[1:filled] contain usable data - filled::Int32 - line::Int - # Whether to iterate copies of the record stored in this struct +mutable struct PAFReader{I <: AbstractBufReader} + # TODO: If Julia ever gets computed field types, we can have this field be typed + # `is_expandable(I) ? I : ExpandReader{I}`. + const io::ExpandReader{I} + # Number of bytes of the last parsed line, including its newline, which have not + # yet been consumed from `io`. Consuming is deferred to the next call to `try_next!`, + # such that the record, which is a view into the buffer of `io`, stays valid. + to_consume::Int + # Int64, not Int, such that it does not overflow for huge files on 32-bit platforms + line::Int64 + # Whether to iterate copies of the parsed records copy::Bool + # Whether a record which is a view into the buffer of `io` has been returned to the + # user, and is still valid. While this is set, `isdone` must not fill the buffer. + # Only ever set when `copy` is false. + view_live::Bool +end + +function PAFReader(io::AbstractBufReader; copy::Bool = true) + return PAFReader{typeof(io)}(ExpandReader(io), 0, 1, copy, false) end function PAFReader(io::IO; buf_size::Int = 2^16, copy::Bool = true) - mem = Memory{UInt8}(undef, max(buf_size, 16)) - rec = PAFRecord(512) - return PAFReader{typeof(io)}(io, rec, mem, 1, 0, 1, copy) + return PAFReader(BufReader(io, max(buf_size, 16)); copy) end Base.IteratorSize(::Type{<:PAFReader}) = Base.SizeUnknown() Base.eltype(::Type{<:PAFReader}) = PAFRecord -Base.close(reader::PAFReader) = close(reader.io) # TODO: Docstring +Base.close(reader::PAFReader) = close(reader.io) -function PAFReader(f::Function, io::IO; kwargs...) +function PAFReader(f::Function, io::Union{IO, AbstractBufReader}; kwargs...) reader = PAFReader(io; kwargs...) return try f(reader) @@ -690,7 +884,9 @@ function Base.copy(record::PAFRecord) return PAFRecord( copy(getfield(record, :data)), getfield(record, :qname_len), + getfield(record, :tname_start), getfield(record, :tname_len), + getfield(record, :aux_start), getfield(record, :qlen), getfield(record, :qstart), getfield(record, :qend), @@ -711,96 +907,24 @@ function Base.iterate(reader::PAFReader, (::Nothing) = nothing) elseif res === nothing nothing else - rec = reader.copy ? copy(res) : res - (rec, nothing) + (res, nothing) end end -# This is never called during normal operations, and only if a single line -# does not fit in the default buffer, which is highly unlikely. -function reallocate!(reader::PAFReader) - mem = reader.mem - newlen = length(mem) * 2 - if newlen >= typemax(Int32) - error("PAFReader encountered a line with more than 2147483646 bytes.") - end - new = Memory{UInt8}(undef, newlen) - unsafe_copyto!(new, 1, mem, 1, reader.filled) - reader.mem = new - return reader -end - -# Shift the data inside the buffer such that the usable -# data begins at index 1, thus making room for more -# data at the end of the buffer -function shift_buffer!(reader::PAFReader) - isone(reader.index) && return reader - mem = reader.mem - len = reader.filled - reader.index + 1 - copyto!(mem, 1, mem, reader.index, len) - reader.index = 1 - reader.filled = len - return reader -end - -# Read data into the buffer. Returns zero if and only if underlying -# stream is EOF (or buggy). -function read_buffer!(reader::PAFReader)::Int - mem = reader.mem - remaining = length(mem) - reader.filled - # If no space, try to shift data in the buffer to make room - if iszero(remaining) - shift_buffer!(reader) - remaining = length(mem) - reader.filled - # If still no space, the buffer is truly full, and we need to - # resize it - if iszero(remaining) - reallocate!(reader) - mem = reader.mem # it was change in reallocate! - remaining = length(mem) - reader.filled - end - end - - @assert !iszero(remaining) - # This part annoys me a lot. Julia has no efficient API to read bytes - # into a byte buffer. How on Earth can that be true? - # * read! will throw EOF if the memory is too short, and there is no - # way to query how many bytes are available - # * readbytes! is not optimised for Memory - # * unsafe_read has almost no promised semantics so practically can't be used. - # The approach here uses unsafe_read in a manner similar to TranscodingStreams.jl. - # It uses bytesavailable to check how many bytes to read and avoid buffer overflows etc, - # and when there are zero bytes available, it reads a single byte, hoping that - # more bytes will be available. - n_available = bytesavailable(reader.io) - if iszero(n_available) && !eof(reader.io) - # Read one byte and hope more bytes become available - byte = read(reader.io, UInt8) - @inbounds mem[reader.filled + 1] = byte - reader.filled += 1 - remaining -= 1 - n_available = bytesavailable(reader.io) - end - # Read with unsafe_read - n_bytes = min(n_available, remaining) - if !iszero(n_bytes) - GC.@preserve mem unsafe_read(reader.io, pointer(mem, reader.filled + 1), n_bytes) - reader.filled += n_bytes - end - return n_bytes -end - """ - try_next!(reader::PAFReader) -> Union{Nothing, PAFRecord, ParserException} + try_next!(reader::PAFReader)::Union{Nothing, PAFRecord, ParserException} Try to read a record from the [`PAFReader`](@ref) and advance the reader. This may yield one of three options: * If the reader is empty, return `nothing` and do not advance the reader * If the next record is invalid, return (do not throw) a `ParserException` and - do not advance the reader + advance the reader past the invalid line. The next call continues from the next line, + so invalid lines may be skipped by ignoring the error. * If the next record is valid, return it as a [`PAFRecord`](@ref) and advance the reader -Even if the reader itself is not advanced, it may still fill its internal buffer, advancing -the `IO` object it wraps. +If the reader was constructed with `copy=true` (the default), the returned record owns +its data. If `copy=false`, the record is a view into the reader's buffer, which is +invalidated when the reader is advanced, and may be invalidated by calling +[`isempty`](@ref Base.isempty(::PAFReader)) on the reader. See [`PAFReader`](@ref). # Examples ```jldoctest @@ -815,79 +939,201 @@ true julia> typeof([try_next!(reader) for i in 1:100]) Vector{Nothing} (alias for Array{Nothing, 1}) -julia> close(reader); reader = PAFReader(IOBuffer("bad data")); +julia> close(reader); reader = PAFReader(IOBuffer("bad data\\n" * open(readline, path_to_paf))); julia> try_next!(reader) isa PAF.ParserException true -julia> all(i -> try_next!(reader) isa PAF.ParserException, 1:100) +julia> try_next!(reader) isa PAFRecord +true + +julia> try_next!(reader) === nothing true ``` """ function try_next!(reader::PAFReader)::Union{Nothing, PAFRecord, ParserException} - # n_scanned: The number of bytes in buffer we have already memchrd searching for a newline - n_scanned = 0 - newline = findfirst( - ==(0x0a), - ImmutableMemoryView(reader.mem)[(reader.index + n_scanned):(reader.filled % Int)], - ) - n_bytes_read = 0 + return expose(reader, next_view!(reader)) +end + +""" + peek(reader::PAFReader)::Union{Nothing, PAFRecord, ParserException} + +Return what [`try_next!`](@ref) would return, without advancing the reader. + +While this does not advance the reader, if the reader was constructed with `copy=false`, +it invalidates any record previously returned by the reader. Similarly, records +returned by `peek` on a reader with `copy=false` are invalidated like records +produced by iterating the reader. + +# Examples +```jldoctest +julia> reader = PAFReader(open(path_to_paf)); + +julia> peek(reader).qname +"query1" + +julia> try_next!(reader).qname +"query1" + +julia> peek(reader).qname +"another query" + +julia> close(reader) +``` +""" +function Base.peek(reader::PAFReader)::Union{Nothing, PAFRecord, ParserException} + # Invalidates the live view, if any, such that the buffer may be filled + consume_previous!(reader) + line = reader.line + res = next_view!(reader) + # Undo the advancement. Since the previous line was consumed above, the + # buffer begins at the peeked line, and nothing needs to be consumed. + reader.to_consume = 0 + reader.line = line + return expose(reader, res) +end + +# Return a parsed view to the user, either as a copy or as a live view +function expose(reader::PAFReader, res::Union{Nothing, PAFRecord, ParserException}) + res isa PAFRecord || return res + # A copy owns its data, so the view never reaches the user + reader.copy && return copy(res) + reader.view_live = true + return res +end + +# Like `try_next!`, but always returns a view, which is not yet exposed to the user +function next_view!(reader::PAFReader)::Union{Nothing, PAFRecord, ParserException} + consume_previous!(reader) + (buffer, last_index, line_len) = buffer_line(reader.io) + line = @inbounds buffer[1:last_index] + # Inlining avoids returning the record through a union in memory, and copying it. + # This was measured to make the reader about 8% faster. + res = @inline parse_line(line) + # A whitespace-only line, which can never be a valid record, is the end of the + # stream if the rest of it is whitespace. Checking only on errors keeps the + # check off the hot path. + if res isa ParserException && all(is_whitespace, line) + tail = check_tail(reader, buffer) + tail === nothing && return nothing + # Skip all the whitespace-only lines, and report the error at the first of them + (n_bytes, n_lines) = tail + err = ParserException(Errors.EmptyNonTailingLine, Int32(1), reader.line) + reader.to_consume = n_bytes + reader.line += n_lines + return err + end + res = if res isa ParserException + ParserException(res.kind, getfield(res, :index_in_line), reader.line) + else + res + end + # Advance past the line, even if it is invalid, such that the next call + # continues from the next line + reader.to_consume = line_len + reader.line += 1 + return res +end + +# Consume the previous line. This invalidates the record if it is a view into the buffer. +function consume_previous!(reader::PAFReader) + if !iszero(reader.to_consume) + consume(reader.io, reader.to_consume) + reader.to_consume = 0 + end + reader.view_live = false + return nothing +end + +# Fill the buffer until it contains a full line, or EOF is reached. +# Return the buffer, the index of the last byte of the line before the \n (or \r\n), +# and the number of bytes in the line including the newline. +function buffer_line(io::ExpandReader)::Tuple{ImmutableMemoryView{UInt8}, Int, Int} + buffer = get_buffer(io) + newline = findfirst(==(0x0a), buffer) # Keep reading until we find a newline or reach EOF while newline === nothing - n_scanned += reader.filled - reader.index - n_scanned + 1 - n_bytes_read = read_buffer!(reader) + # The number of bytes in buffer we have already memchrd searching for a newline + n_scanned = length(buffer) + n_filled = fill_buffer(io) + # Filling invalidates the old view, even if no bytes were added + buffer = get_buffer(io) # This signals EOF - we might still have data in the buffer, if the underlying stream # ends with a record that doesn't end in a newline - iszero(n_bytes_read) && break - newline = findfirst( - ==(0x0a), - ImmutableMemoryView(reader.mem)[(reader.index + n_scanned):(reader.filled % Int)], - ) + iszero(n_filled) && return (buffer, length(buffer), length(buffer)) + newline = findnext(==(0x0a), buffer, n_scanned + 1) end - # If a newline was found, the current value of newline was relative to where - # the last search began. So, update it to account for the bytes searched in - # previous iterations of search for tab - newline === nothing || (newline += reader.index + n_scanned - 1) - mem = reader.mem - has_windows_newline = false # Last index before \n, and before \r if the newline is \r\n as on Windows - last_index = if newline === nothing - reader.filled - else - if isone(newline) - 0 - else - has_windows_newline = @inbounds(mem[newline - 1]) == UInt8('\r') - newline - (1 + has_windows_newline) - end - end - parse_mem = ImmutableMemoryView(mem)[reader.index:last_index] - # If there is no memory to parse (and we have already verified that - # there are no more bytes in the IO), we end - if isempty(parse_mem) - line = reader.line - col = 0 - @inbounds for b in ImmutableMemoryView(mem)[reader.index:reader.filled] - col += 1 - line += b == UInt8('\n') - col *= b != UInt8('\n') - if b ∉ (UInt8('\r'), UInt8('\n'), UInt8('\t'), UInt8(' ')) - return ParserException(Errors.EmptyNonTailingLine, col % Int32, line) + has_windows_newline = !isone(newline) && @inbounds(buffer[newline - 1]) == UInt8('\r') + return (buffer, newline - (1 + has_windows_newline), newline) +end + +is_whitespace(b::UInt8) = b ∈ (UInt8('\r'), UInt8('\n'), UInt8('\t'), UInt8(' ')) + +# Called when the next line is whitespace-only, with `buffer` beginning at that line. +# Return `nothing` if the rest of the stream is whitespace. Else, return the number of +# bytes and lines of the whitespace-only lines before the first line with non-whitespace. +# The buffer may end after the line, so keep filling it until EOF. +# This does not mutate `reader`, except by filling its buffer. +function check_tail(reader::PAFReader, buffer::ImmutableMemoryView{UInt8})::Union{Nothing, Tuple{Int, Int}} + io = reader.io + n_lines = 0 + # Index of the first byte after the last newline seen + line_start = 1 + n_scanned = 0 + while true + for i in (n_scanned + 1):length(buffer) + b = @inbounds buffer[i] + if b == UInt8('\n') + n_lines += 1 + line_start = i + 1 + elseif !is_whitespace(b) + return (line_start - 1, n_lines) end end - return nothing + n_scanned = length(buffer) + iszero(fill_buffer(io)) && return nothing + buffer = get_buffer(io) end - # If this results in an error, we return the error. - res = parse_line!(reader.record, parse_mem) - res isa ParserException && - return ParserException(res.kind, getfield(res, :index_in_line), reader.line) - # Else, we update the reader. This means `next!` will continuously return - # an error if an error is found. - # If the newline was found, we compensate for the newline by skipping that byte in - # the next iteration, hence the extra `newline !== nothing`. - reader.index = last_index + 1 + has_windows_newline + (newline !== nothing) - reader.line += 1 - return res + return +end + +# This does not advance the reader, and never invalidates a returned record. +# If that is not possible, because a returned record is a view into the buffer and the +# buffer must be filled to determine the answer, return `missing`. +function Base.isdone(reader::PAFReader, _...)::Union{Bool, Missing} + rest = get_buffer(reader.io)[(reader.to_consume + 1):end] + # The reader is done exactly when the rest of the stream is whitespace: + # If the next line contains non-whitespace, `try_next!` returns a record or an error. + # Else, it is the end of the stream only if the rest of it is whitespace. + # This usually returns at the first byte. + all(is_whitespace, rest) || return false + # If no view is live, we can freely consume the previous line and fill the buffer. + # Else, we can't tell if the stream continues with non-whitespace. + return reader.view_live ? missing : isdone_filling!(reader) +end + +function isdone_filling!(reader::PAFReader)::Bool + consume_previous!(reader) + return check_tail(reader, get_buffer(reader.io)) === nothing +end + +""" + isempty(reader::PAFReader)::Bool + +Check whether `reader` has no more records, i.e. whether [`try_next!`](@ref) would return +`nothing`. This does not advance the reader, but if the reader was constructed with +`copy=false`, it may invalidate the previous record by mutating the reader's inner state. + +In contrast, `Base.isdone(reader)`, never invalidates a record, +and returns `missing` if it can't determine whether the reader is done +without invalidating previous records. +""" +function Base.isempty(reader::PAFReader) + done = Base.isdone(reader) + done === missing || return done + # Invalidate the live view, such that `isdone` may fill the buffer + return isdone_filling!(reader) end @setup_workload begin diff --git a/test/runtests.jl b/test/runtests.jl index 6d1b5bb..2a312bc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,7 +6,10 @@ using Test using FormatSpecimens using MemoryViews using Aqua +using BufferIO: CursorReader +using BufferIO.Testing: BoundedReader, ChunkedReader +using XAMAuxData: XAMAuxData using XAMAuxData.SAM: AuxTag function cmp_aux(record::PAFRecord, other) @@ -15,7 +18,7 @@ function cmp_aux(record::PAFRecord, other) @test length(aux) == length(other) for ((kaux, vaux), (ko, vo)) in zip(aux, other) @test kaux == AuxTag(ko) - if vaux isa Union{Real, AbstractString, Char} + if vaux isa Union{Real, AbstractString, Char, XAMAuxData.Error} @test vaux == vo elseif vaux isa AbstractVector @test eltype(vaux) == eltype(vo) @@ -44,6 +47,10 @@ end @test rec1.alnlen == 301142 @test rec1.mapq == 0 + # tname and is_rc are type stable, since they throw for unmapped records + @test (@inferred (r -> r.tname)(rec1)) isa AbstractString + @test (@inferred (r -> r.is_rc)(rec1)) isa Bool + @test_throws Exception rec1.foobar @test_throws Exception rec1.strand @test_throws Exception rec1.data @@ -58,6 +65,54 @@ end end end +@testset "Raw line" begin + line = "q\t10\t0\t5\t+\tt\t20\t5\t10\t5\t5\t60\tab:Z:x" + for input in (line * "\n", line * "\r\n", line) + rec = parse(PAFRecord, input) + @test MemoryView(rec) isa ImmutableMemoryView{UInt8} + @test MemoryView(rec) == codeunits(line) + @test String(rec) isa String + @test String(rec) == line + @test string(rec) isa String + @test string(rec) == line + @test sprint(print, rec) == line + end + unmapped = "read1\t1523\t0\t0\t*\t*\t0\t0\t0\t0\t0\t0" + @test String(parse(PAFRecord, unmapped)) == unmapped + + # String and string copy the data, MemoryView does not + data = collect(codeunits(line)) + rec = parse(PAFRecord, data) + (s1, s2, mem) = (String(rec), string(rec), MemoryView(rec)) + data[1] = UInt8('Q') + @test s1 == s2 == line + @test first(mem) == UInt8('Q') +end + +@testset "Equality and hashing" begin + line = "q\t10\t0\t5\t+\tt\t20\t5\t10\t5\t5\t60\tab:Z:x" + a = parse(PAFRecord, line) + for b in ( + parse(PAFRecord, line * "\n"), + parse(PAFRecord, line * "\r\n"), + parse(PAFRecord, collect(codeunits(line))), + copy(a), + ) + @test a == b + @test isequal(a, b) + @test hash(a) == hash(b) + end + c = parse(PAFRecord, line * "\tcd:i:1") + @test a != c + @test hash(a) != hash(c) + @test hash(a) != hash(line) + @test length(Set([a, copy(a), c])) == 2 + + unmapped = "read1\t1523\t0\t0\t*\t*\t0\t0\t0\t0\t0\t0" + @test parse(PAFRecord, unmapped) == parse(PAFRecord, unmapped * "\n") + @test parse(PAFRecord, unmapped) != a +end + @testset "Coverage and identity" begin function get_rec(qlen::Int, tlen::Int, matches::Int, alnlen::Int) fields = split( @@ -83,6 +138,17 @@ end @test isapprox(target_coverage(get_rec(17, 16, 15, 16)), 6 / 16) end +function with_replaced_str(s::String, field::Integer, new::String) + fields = split(s, '\t') + fields[field] = new + return join(fields, "\t") +end + +function with_replaced(s::String, field::Integer, new::String) + y = try_parse(with_replaced_str(s, field, new)) + return y isa PairwiseMappingFormat.ParserException ? y.kind : y +end + @testset "Edge case records" begin (rec1, rec2, rec3) = PAFReader(open(joinpath(path_of_format("PAF"), "good2.paf"))) do reader y = first(reader, 3) @@ -91,15 +157,21 @@ end end @test rec1.qname == "中国科学院大学" @test rec1.qlen == typemax(Int32) - cmp_aux(rec1, ["tp" => 'S', "cm" => 29990, "s1" => 301142, "dv" => -13.2442f-18, "rl" => 0, "xX" => [0xea, 0x19, 0xfe, 0x1b], "K1" => Int8[15, 9, 22, 127, -12, -128, 0]]) + # XAMAuxData 0.2 rejects lowercase hex digits + xX = pkgversion(XAMAuxData) < v"0.2" ? [0xea, 0x19, 0xfe, 0x1b] : XAMAuxData.Errors.InvalidHex + cmp_aux(rec1, ["tp" => 'S', "cm" => 29990, "s1" => 301142, "dv" => -13.2442f-18, "rl" => 0, "xX" => xX, "K1" => Int8[15, 9, 22, 127, -12, -128, 0]]) cmp_aux(rec2, ["dv" => 1.243f-4, "rl" => 38, "kk" => "!=[]! And then~~ "]) cmp_aux(rec3, ["dv" => +17.0f21, "rl" => 38]) # Empty record rec = PAFReader(first, open(joinpath(path_of_format("PAF"), "good3.paf"))) - @test isnothing(rec.tname) - @test isnothing(rec.is_rc) @test !is_mapped(rec) + for sym in (:tname, :is_rc, :mapq, :qstart, :qend, :tlen, :tstart, :tend, :matches, :alnlen) + @test_throws PairwiseMappingFormat.UnmappedError(sym) getproperty(rec, sym) + end + for f in (query_coverage, target_coverage, aln_identity) + @test_throws PairwiseMappingFormat.UnmappedError(Symbol(f)) f(rec) + end @test rec.qname == "query" @test rec.qlen == 56 @test Dict(aux_data(rec)) == Dict{AuxTag, Any}( @@ -110,15 +182,62 @@ end "rl" => 0, ) + # Unmapped record as written by minimap2 --paf-no-hit, with zeros + rec = parse(PAFRecord, "read1\t1523\t0\t0\t*\t*\t0\t0\t0\t0\t0\t0") + @test !is_mapped(rec) + @test rec.qname == "read1" + @test rec.qlen == 1523 + # Unmapped record with too few fieds @test try_parse("query\t123\t1\t100\t*\t*\t5\t19").kind == Errors.TooFewFields -end -function with_replaced(s::String, field::Integer, new::String) - fields = split(s, '\t') - fields[field] = new - y = try_parse(join(fields, "\t")) - return y isa PairwiseMappingFormat.ParserException ? y.kind : y + # The fields of unmapped records are validated like for mapped records, + # except that integers may be zero + unmapped = "read1\t1523\t0\t0\t*\t*\t0\t0\t0\t0\t0\t0" + for field in 6:12 + @test with_replaced(unmapped, field, "0") isa PAFRecord + end + # Use a longer target name, so emptying a field does not make the line too short + long_unmapped = with_replaced_str(unmapped, 6, "some_target_name") + for field in 7:12 + @test with_replaced(long_unmapped, field, "") == Errors.EmptyInteger + @test with_replaced(unmapped, field, "-1") == Errors.BadInteger + @test with_replaced(unmapped, field, "abc") == Errors.BadInteger + @test with_replaced(unmapped, field, "4294967296") == Errors.IntegerOverflow + end + @test with_replaced(unmapped, 12, "255") isa PAFRecord + @test with_replaced(unmapped, 12, "256") == Errors.IntegerOverflow + @test with_replaced(unmapped, 9, "1") == Errors.PositionOutOfBounds + @test try_parse("read1\t1523\t0\t0\t*\t*\t10\t5\t4\t0\t0\t0").kind == Errors.BackwardsIndices + @test try_parse("read1\t1523\t0\t0\t*\t*\t10\t0\t5\t0\t0\t0") isa PAFRecord + @test try_parse("read1\t1523\t0\t0\t*x*\t10\t0\t5\t0\t0\t0").kind == Errors.InvalidStrand + # Aux data is found after the mapq field + rec = try_parse("read1\t1523\t0\t0\t*\t*\t0\t0\t0\t0\t0\t0\tcm:i:4") + @test Dict(aux_data(rec)) == Dict{AuxTag, Any}("cm" => 4) + + # Shortest possible records: Empty names and single-digit integers + rec = parse(PAFRecord, "\t1\t0\t1\t+\t\t1\t0\t1\t1\t1\t0") + @test rec.qname == rec.tname == "" + @test (rec.qlen, rec.tlen, rec.qstart, rec.qend, rec.mapq) == (1, 1, 1, 1, 0) + @test isempty(aux_data(rec)) + rec = parse(PAFRecord, "\t1\t0\t0\t*\t\t0\t0\t0\t0\t0\t0") + @test !is_mapped(rec) + @test rec.qname == "" + # Removing any byte from them makes them invalid + for line in ("\t1\t0\t1\t+\t\t1\t0\t1\t1\t1\t0", "\t1\t0\t0\t*\t\t0\t0\t0\t0\t0\t0") + for i in 1:ncodeunits(line) + @test try_parse(line[1:(i - 1)] * line[(i + 1):end]) isa PairwiseMappingFormat.ParserException + end + end + + # The empty record is consistent with its line + rec = PAFRecord() + @test !is_mapped(rec) + @test rec.qname == "" + @test rec.qlen == 1 + @test isempty(aux_data(rec)) + @test parse(PAFRecord, String(rec)) == rec + @test String(rec) == "\t1\t0\t0\t*\t*\t0\t0\t0\t0\t0\t0" end @testset "Errors" begin @@ -132,10 +251,14 @@ end @test with_replaced(good, allow_zero, "0") isa PAFRecord end - for not_allow_zero in [2, 4, 7, 9] + for not_allow_zero in [2, 7] @test with_replaced(good, not_allow_zero, "0") == Errors.InvalidZero end + # Ends may be zero, but only if the start is too, giving an empty interval + @test with_replaced(good, 4, "0") == Errors.BackwardsIndices + @test with_replaced(good, 9, "0") == Errors.BackwardsIndices + # qend and tend must be >= qstart tstart @test with_replaced(good, 3, "1000000") == Errors.BackwardsIndices @test with_replaced(good, 4, "2") == Errors.BackwardsIndices @@ -151,6 +274,10 @@ end # Char not valid @test with_replaced(good, 5, "/") == Errors.InvalidStrand + # The strand must be exactly one byte + for strand in ("", "+-", "++", "--", "-+", "**", "+ ", "+\r") + @test with_replaced(good, 5, strand) == Errors.InvalidStrand + end # Too few columns for i in 0:11 @@ -166,6 +293,86 @@ end # Mapq too high @test with_replaced(good, 12, "256") == Errors.IntegerOverflow + + # Integer overflow, including on the last digit + @test with_replaced(good, 10, "4294967295").matches == typemax(UInt32) + @test with_replaced(good, 10, "4294967296") == Errors.IntegerOverflow + @test with_replaced(good, 10, "9999999999") == Errors.IntegerOverflow + @test with_replaced(good, 10, "99999999999999999999") == Errors.IntegerOverflow + @test with_replaced(good, 12, "4294967296") == Errors.IntegerOverflow + @test with_replaced(good, 2, "4294967296") == Errors.IntegerOverflow + + # Sequences longer than typemax(Int32), up to typemax(UInt32) + m = typemax(UInt32) + rec = try_parse("q\t$m\t$(m - 1)\t$m\t+\tt\t$m\t0\t$m\t1\t1\t0") + @test (rec.qlen, rec.qstart, rec.qend) == (m, m, m) + # Int64, not Int, so the values fit on 32-bit platforms + for sym in (:qlen, :qstart, :qend, :tlen, :tstart, :tend, :matches, :alnlen) + @test getproperty(rec, sym) isa Int64 + end + @test (rec.tlen, rec.tstart, rec.tend) == (m, 1, m) + @test query_coverage(rec) == 1 / m + @test target_coverage(rec) == 1.0 + + # Empty interval at the very end of the sequence doesn't overflow + rec = try_parse("q\t$m\t$m\t$m\t+\tt\t5\t0\t5\t1\t1\t0") + @test (rec.qstart, rec.qend) == (Int64(m) + 1, m) +end + +@testset "Integer parsing" begin + # Integers are parsed 8 bytes at a time, unless they have more than 7 digits, + # or are close to the end of the line. Test both cases for all digit counts. + good = "my_qname\t301156\t3\t301145\t+\tmy_tname\t6701780\t2764860\t3066002\t301142\t301142\t0" + m = Int64(typemax(UInt32)) + alnlen(x) = x isa PAFRecord ? x.alnlen : x + for n_digits in 1:11, n_zeros in 0:9 + s = "0"^n_zeros * first("9182736450918273645", n_digits) + n = parse(Int64, s) + expected = n > m ? Errors.IntegerOverflow : n + # alnlen is close to the end of the line, except when there is aux data + @test alnlen(with_replaced(good, 11, s)) == expected + @test alnlen(with_replaced(good * "\tcm:i:1", 11, s)) == expected + # Same with a bad byte after the integer + @test with_replaced(good, 11, s * "/") == (n > m ? Errors.IntegerOverflow : Errors.BadInteger) + end + + # A bad byte at every position of an integer, including bytes which may carry + # into the next byte when adding to them. + for n_digits in 1:12, pos in 1:n_digits, bad in (0x2f, 0x3a, 0x20, 0x00, 0xfa, 0xff, UInt8('a')) + bytes = collect(codeunits(first("123456789123", n_digits))) + bytes[pos] = bad + s = String(bytes) + for field in (10, 11, 12) + # 1234567891 fits in an UInt32, so only 11 digits overflow before the bad byte + expected = pos > 11 ? Errors.IntegerOverflow : Errors.BadInteger + @test with_replaced(good, field, s) == expected + end + end +end + +@testset "Empty intervals" begin + # Empty intervals are accepted anywhere within the sequence, including at its + # start and end. They become empty one-based ranges. + for (start, stop) in [(0, 0), (5, 5), (10, 10)] + rec = parse(PAFRecord, "q\t10\t$start\t$stop\t+\tt\t10\t$start\t$stop\t0\t0\t0") + @test (rec.qstart, rec.qend) == (start + 1, stop) + @test (rec.tstart, rec.tend) == (start + 1, stop) + @test isempty(rec.qstart:rec.qend) + @test isempty("0123456789"[rec.tstart:rec.tend]) + @test query_coverage(rec) == 0.0 + @test target_coverage(rec) == 0.0 + @test isnan(aln_identity(rec)) + end + + # But not out of bounds + @test try_parse("q\t10\t11\t11\t+\tt\t10\t0\t0\t0\t0\t0").kind == Errors.PositionOutOfBounds + @test try_parse("q\t10\t0\t0\t+\tt\t10\t11\t11\t0\t0\t0").kind == Errors.PositionOutOfBounds + + # matches and alnlen are not validated against each other + rec = parse(PAFRecord, "q\t10\t0\t5\t+\tt\t10\t0\t5\t7\t3\t0") + @test aln_identity(rec) == 7 / 3 + rec = parse(PAFRecord, "q\t10\t0\t5\t+\tt\t10\t0\t5\t2\t0\t0") + @test aln_identity(rec) == Inf end @testset "Bad reader" begin @@ -193,6 +400,53 @@ end end end +@testset "AbstractBufReader" begin + good = "my_qname\t301156\t3\t301145\t+\tmy_tname\t6701780\t2764860\t3066002\t301142\t301142\t0" + other = replace(good, "my_qname" => "other") + data = join([good, other, good], '\n') + + # Records iterated with copy=true do not alias the reader's buffer + records = PAFReader(collect, CursorReader(data)) + @test [r.qname for r in records] == ["my_qname", "other", "my_qname"] + + # With copy=false, each record is a distinct view into the reader's buffer. + # CursorReader never moves its data, so the views stay valid here. + PAFReader(CursorReader(data); copy = false) do reader + rec = try_next!(reader) + @test rec.qname == "my_qname" + rec2 = try_next!(reader) + @test rec2 !== rec + @test rec2.qname == "other" + end + records = PAFReader(collect, CursorReader(data); copy = false) + @test [r.qname for r in records] == ["my_qname", "other", "my_qname"] + + # try_next! respects copy: With copy=true, records stay valid after advancing, + # even when the reader's buffer is overwritten + reader = PAFReader(ChunkedReader([good * "\n", other * "\n"])) + rec = try_next!(reader) + @test try_next!(reader).qname == "other" + @test rec.qname == "my_qname" + + # The last line may lack a newline, also when filling the buffer at EOF + # invalidates views of it, and when the buffer can't hold the whole line + data = good * "\n" * other + for max_buffer_size in (4, 200), copy in (true, false) + reader = PAFReader(BoundedReader(data, max_buffer_size; max_fill = 3); copy) + @test try_next!(reader).qname == "my_qname" + @test try_next!(reader).qname == "other" + @test try_next!(reader) === nothing + end + + # Parsing a vector aliases it, but copying the record does not + v = collect(codeunits(good)) + rec = try_parse(v) + copied = copy(rec) + v[1:2] .= UInt8('X') + @test rec.qname == "XX_qname" + @test copied.qname == "my_qname" +end + # We assume that aux data is correctly tested in its own package, here we just test # that the aux data has the right span and byte content. function test_aux_content(record_data::String, target_aux::String) @@ -209,11 +463,48 @@ end test_aux_content(good * stringaux_with_whitespace, goodaux * stringaux_with_whitespace) end -# Note: We currently DONT trim whitespace off ends of records, for the sole reason that -# whitespace is allowed in string AuxData. -# This is very annoying, because it makes parsing less robust to erroneous trailing whitespace. -# However, I can't think of a good way to deal with this. -# TODO: Find a better method +@testset "Parsing lines" begin + good = "my_qname\t301156\t3\t301145\t+\tmy_tname\t6701780\t2764860\t3066002\t301142\t301142\t0" + expected = parse(PAFRecord, good) + same(rec) = rec isa PAFRecord && getfield(rec, :data) == getfield(expected, :data) + + # parse accepts the same inputs as try_parse + for x in (good, SubString(good, 1), codeunits(good), collect(codeunits(good)), MemoryView(good)) + @test same(parse(PAFRecord, x)) + @test same(try_parse(x)) + end + @test_throws PairwiseMappingFormat.ParserException parse(PAFRecord, codeunits("abc")) + + # One optional trailing newline is stripped, also before aux data + @test same(try_parse(good * "\n")) + @test same(try_parse(good * "\r\n")) + @test aux_data(parse(PAFRecord, good * "\tab:Z:x\r\n"))["ab"] == "x" + + # Any other newline is an error + @test try_parse(good * "\n\n").kind == Errors.EmbeddedNewline + @test try_parse(good * "\r\n\n").kind == Errors.EmbeddedNewline + @test try_parse(good * "\n" * good).kind == Errors.EmbeddedNewline + @test try_parse(replace(good, "my_qname" => "my\nqname")).kind == Errors.EmbeddedNewline + @test try_parse(good * "\tab:Z:x\ny").kind == Errors.EmbeddedNewline + @test try_parse("\n").kind == Errors.TooFewFields + + # An embedded \r is data + @test parse(PAFRecord, replace(good, "my_qname" => "my\rqname")).qname == "my\rqname" + + # Other trailing whitespace is not stripped, but part of the last field + @test try_parse(good * " ").kind == Errors.BadInteger + @test try_parse(good * "\r").kind == Errors.BadInteger + # A trailing tab begins an empty aux field, which is only invalid when accessed + @test try_parse(good * "\t") isa PAFRecord + @test try_parse(good * " \n").kind == Errors.BadInteger + @test aux_data(parse(PAFRecord, good * "\tab:Z:x \r\n"))["ab"] == "x " + # A trailing \r is kept in aux data, where it is invalid for Z values + @test collect(MemoryView(aux_data(parse(PAFRecord, good * "\tab:Z:x\r")))) == codeunits("ab:Z:x\r") +end + +# Note: We don't trim whitespace off ends of records, because whitespace is allowed in +# string AuxData. Trailing whitespace is part of the last field, and may make it invalid. +# After the records, a file may end with whitespace-only lines. @testset "Trailing whitespace" begin data = "my_qname 301156 3 301145 + my_tname 6701780 2764860 3066002 301142 301142 0 an:Z: dljdl \r\n" data *= "some 5 3 4 + other 9 5 8 2 2 0\r\n" @@ -231,6 +522,243 @@ end @test first(reader, 2) isa Vector{PAFRecord} @test try_next!(reader).kind == Errors.EmptyNonTailingLine end + + # Blank line in the middle, where the buffer ends right after the blank line + good = "my_qname\t301156\t3\t301145\t+\tmy_tname\t6701780\t2764860\t3066002\t301142\t301142\t0" + reader = PAFReader(ChunkedReader([good * "\n\n", good * "\n" * good * "\n"])) + @test try_next!(reader) isa PAFRecord + @test try_next!(reader).kind == Errors.EmptyNonTailingLine + + # Trailing whitespace delivered in multiple chunks is still fine + reader = PAFReader(ChunkedReader([good * "\n\n", "\r\n", "\n"])) + @test length(collect(reader)) == 1 + + # The first whitespace-only line ends the records, whether or not it is empty + for tail in ("\n \n", "\n \n\n", "\n\t\r\n \n", "\n ", "\n\r", "\n\r\r\n", "\n\n \n") + @test length(PAFReader(collect, IOBuffer(good * tail))) == 1 + reader = PAFReader(IOBuffer(good * tail)) + @test try_next!(reader) isa PAFRecord + @test isempty(reader) + @test try_next!(reader) === nothing + end + @test isempty(PAFReader(IOBuffer(" \n\t\n"))) + + # ... after which only whitespace may follow + for tail in ("\n \n" * good, "\n\t\nx", "\n \r\n\n" * good * "\n") + reader = PAFReader(IOBuffer(good * tail)) + @test try_next!(reader) isa PAFRecord + @test !isempty(reader) + @test try_next!(reader).kind == Errors.EmptyNonTailingLine + end + + # Trailing whitespace on a record is not stripped + reader = PAFReader(IOBuffer(good * " \n")) + @test try_next!(reader).kind == Errors.BadInteger + reader = PAFReader(IOBuffer(good * "\tab:Z:x \n" * good * "\t\n")) + @test aux_data(try_next!(reader))["ab"] == "x " + @test try_next!(reader) isa PAFRecord # empty aux field is only invalid when accessed + + # A line with non-whitespace after leading whitespace is an ordinary error + reader = PAFReader(IOBuffer(good * "\n \tx\n")) + @test try_next!(reader) isa PAFRecord + @test try_next!(reader).kind == Errors.BadInteger +end + +@testset "Skipping invalid lines" begin + good = "my_qname\t301156\t3\t301145\t+\tmy_tname\t6701780\t2764860\t3066002\t301142\t301142\t0" + kinds(results) = [r isa PAFRecord ? :record : (r.kind, r.line) for r in results] + function read_all(reader) + results = [] + while (res = try_next!(reader)) !== nothing + push!(results, res) + end + return results + end + + # try_next! advances past invalid lines + data = join([good, "bad", good, good * "\tx", good * " ", good], '\n') + @test kinds(read_all(PAFReader(IOBuffer(data)))) == + [:record, (Errors.TooFewFields, 2), :record, :record, (Errors.BadInteger, 5), :record] + + # Including a truncated last line, as from a killed job + reader = PAFReader(IOBuffer(good * "\n" * good[1:30])) + @test try_next!(reader) isa PAFRecord + @test try_next!(reader).line == 2 + @test try_next!(reader) === nothing + @test isempty(reader) + + # A run of whitespace-only lines followed by non-whitespace is one error, + # reported at the first whitespace-only line. Reading continues after the run. + data = good * "\n\n \r\n\t\n" * good * "\n\nbad\n\n" + @test kinds(read_all(PAFReader(IOBuffer(data)))) == [ + :record, (Errors.EmptyNonTailingLine, 2), :record, + (Errors.EmptyNonTailingLine, 6), (Errors.TooFewFields, 7), + ] + + # Iteration throws, but also advances, so it can be resumed after catching + reader = PAFReader(IOBuffer(join([good, "bad", good], '\n'))) + @test iterate(reader) isa Tuple + @test_throws PairwiseMappingFormat.ParserException iterate(reader) + @test iterate(reader) isa Tuple + @test iterate(reader) === nothing + + # Mixing iteration and try_next! gives the same results + reader = PAFReader(IOBuffer(join(["bad", "bad", good], '\n'))) + @test_throws PairwiseMappingFormat.ParserException iterate(reader) + @test try_next!(reader).line == 2 + @test only(collect(reader)).qname == "my_qname" +end + +@testset "peek" begin + good = open(readline, joinpath(@__DIR__, "example.paf")) + other = replace(good, "query1" => "other") + data = good * "\nbad\n\n \n" * other * "\r\n\n \n" + results(r) = r isa PAFRecord ? String(r) : r isa PairwiseMappingFormat.ParserException ? (r.kind, r.line) : r + expected = PAFReader(IOBuffer(data)) do reader + out = [] + while true + res = try_next!(reader) + push!(out, results(res)) + res === nothing && return out + end + end + @test length(expected) == 5 # good, bad, empty non-tailing line, other, nothing + + for copy in (true, false), chunk in (1, 3, 1000), n_peeks in (0, 1, 2) + chunks = [data[i:min(end, i + chunk - 1)] for i in 1:chunk:lastindex(data)] + reader = PAFReader(ChunkedReader(chunks); copy) + for exp in expected + for _ in 1:n_peeks + @test results(peek(reader)) == exp + end + @test results(try_next!(reader)) == exp + end + @test peek(reader) === nothing + @test isempty(reader) + end + + # A peeked record is valid until the reader is advanced, also with copy=false + PAFReader(IOBuffer(good * "\n" * other); copy = false) do reader + rec = peek(reader) + @test !isempty(reader) + @test Base.isdone(reader) === false + @test String(rec) == good + @test String(try_next!(reader)) == good + @test String(peek(reader)) == other + @test String(try_next!(reader)) == other + @test peek(reader) === nothing + end + + # With copy=true, peeked records own their data + PAFReader(IOBuffer(good * "\n" * other)) do reader + rec = peek(reader) + try_next!(reader) + try_next!(reader) + @test String(rec) == good + end +end + +@testset "isempty" begin + good = "my_qname\t301156\t3\t301145\t+\tmy_tname\t6701780\t2764860\t3066002\t301142\t301142\t0" + + # isempty does not consume records + reader = PAFReader(IOBuffer(join(fill(good, 3), '\n'))) + @test !isempty(reader) + @test !isempty(reader) + @test length(collect(reader)) == 3 + @test isempty(reader) + + # isempty does not invalidate records iterated with copy=true, even when it needs to + # fill the buffer + expected = parse(PAFRecord, good) + for chunks in ([good * "\n", good * "\n"], [good * "\n", "\n", "\r\n"], [good]) + reader = PAFReader(ChunkedReader(chunks)) + (rec, _) = iterate(reader) + while !isempty(reader) + @test rec.qname == expected.qname + @test rec.tname == expected.tname + @test aux_data(rec) == aux_data(expected) + (rec, _) = iterate(reader) + end + @test rec.qname == expected.qname + @test rec.tname == expected.tname + @test iterate(reader) === nothing + end + + # isdone does not fill the buffer when a view returned by try_next! is live, + # but isempty does + reader = PAFReader(ChunkedReader([good * "\n", good * "\n"]); copy = false) + rec = try_next!(reader) + @test Base.isdone(reader) === missing + @test rec.qname == "my_qname" + @test !isempty(reader) + @test try_next!(reader) isa PAFRecord + @test isempty(reader) + + inputs = [ + "", "\n", "\r\n", "\n\n \t\r\n", good, good * "\n", good * "\r\n\n", + good * "\n\n" * good, good * "\nbad", good * "\n \n", good * "\n\r", + good * "\n\r\r\n", " ", "\nV", good * "\n\n\n\nx", + good * "\n \n" * good, good * "\n\t \r\n ", " \tx", good * "\n \n\tx", + "bad\n" * good, good * "\nbad\n" * good * "\n", "bad\r\nbad\n\n \n" * good, + ] + # Empty chunks are not allowed, since fill_buffer returning zero signals EOF + function chunkings(data) + splits = [filter!(!isempty, [data[1:i], data[(i + 1):end]]) for i in 0:length(data)] + return push!(splits, [string(c) for c in data]) + end + # Summarize the results of reading all of `reader`, continuing past errors. + # Each call must advance the reader, so the number of results is bounded. + function read_all(reader, max_results) + results = [] + while (res = try_next!(reader)) !== nothing + push!(results, res isa PAFRecord ? copy(getfield(res, :data)) : res) + length(results) > max_results && error("Reader did not advance") + end + return results + end + for data in inputs + expected_results = read_all(PAFReader(IOBuffer(data)), ncodeunits(data)) + for chunks in chunkings(data) + @test read_all(PAFReader(ChunkedReader(chunks)), ncodeunits(data)) == expected_results + end + end + + for data in inputs, chunks in chunkings(data), copy in (true, false) + # isdone and isempty agree with try_next!, and isdone never invalidates + # the record returned by try_next!. With copy=true, no view is live, so + # isdone is never missing, and isempty does not invalidate the record either. + # Reading continues past errors. + reader = PAFReader(ChunkedReader(chunks); copy) + (rec, expected) = (nothing, nothing) + while true + maybe_done = Base.isdone(reader) + copy && @test maybe_done isa Bool + rec === nothing || @test getfield(rec, :data) == getfield(expected, :data) + done = isempty(reader) + copy && rec !== nothing && @test getfield(rec, :data) == getfield(expected, :data) + maybe_done === missing || @test maybe_done == done + res = try_next!(reader) + @test done == (res === nothing) + res === nothing && break + (rec, expected) = res isa PAFRecord ? (res, Base.copy(res)) : (nothing, nothing) + end + + # When iterating with copy=true, no view is live, so isdone is never missing + reader = PAFReader(ChunkedReader(chunks)) + while true + done = Base.isdone(reader) + @test done isa Bool + it = try + iterate(reader) + catch e + e isa PairwiseMappingFormat.ParserException || rethrow() + e + end + @test done == (it === nothing) + it === nothing && break + end + end end @testset "Aqua" begin