PairwiseMappingFormat.jl

This package is for reading files of the PAF (Pairwise mApping Format) format, 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 type wraps an IO or a BufferIO.AbstractBufReader, and is an iterator of PAFRecord objects:

reader = PAFReader(open(path_to_paf))
records = collect(reader)
@assert isempty(reader)
close(reader)

@assert typeof(records) == Vector{PAFRecord}

Similar to the common open(path) do io-syntax in Julia, PAFReader takes an optional f::Function as a first argument, in which case it will apply f to the returned PAFReader object, and close the underlying io when f returns (or errors):

julia> PAFReader(open(path_to_paf)) do reader
           for record in reader
               println(record.qlen)
           end
       end
301156
299273
288659

The PAFReader constructor takes an optional keyword copy, which defaults to true. If it is false, the PAFRecords returned by iterating the reader, or by try_next!, 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:

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 PAFRecord object represents a single line in a PAF file. The individual columns of the PAF line is obtained by accessing properties of the records:

julia> record = PAFReader(first, open(path_to_paf));

julia> println(record.qname)
query1

julia> println(record.qlen)
301156

julia> println(record.mapq)
0

Note that Base.getproperty is overloaded for PAFRecord, so the properties are public and stable across major versions, but may not reflect the actual underlying fields as they are stored in the PAFRecord.

Fields of `PAFRecord`

qname::StringView{A}. A is an implementation detail. 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::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 Int64. Rightmost (i.e. highest) position of the alignment in the query and target, respectively. This does not take strandedness into account. 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::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.
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::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, query_coverage and target_coverage.

Besides being iterated from PAFReader, records can also be created by parsing from a bytes-like object such as a String:

data = "contig_11\t288659\t27\t288653\t+\tCP004047.1\t6701780\t5027225\t5315851\t288618\t288626\t0\ttp:A:P\tcm:i:28871\ts1:i:288618\ts2:i:288618\tdv:f:0.0000\trl:i:57"

record = parse(PAFRecord, data)
@assert record isa PAFRecord

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. You can check if a record is mapped or unmapped with the function is_mapped:

julia> unmapped = parse(PAFRecord, "myquery\t5032251\t0\t0\t*\t*\t0\t0\t0\t0\t0\t0");

julia> is_mapped(unmapped)
false

julia> println(unmapped.qname)
myquery

julia> println(unmapped.qlen)
5032251

Accessing any other property of an unmapped record throws a PairwiseMappingFormat.UnmappedError, as does calling target_coverage, query_coverage or aln_identity on it:

julia> unmapped.tname
ERROR: UnmappedError: Cannot get `tname` of an unmapped PAFRecord. Check `is_mapped` first.
[...]

julia> query_coverage(unmapped)
ERROR: UnmappedError: Cannot get `query_coverage` of an unmapped PAFRecord. Check `is_mapped` first.
[...]

Low-level interface

Iterating PAFReaders, and the parse function will throw a PairwiseMappingFormat.ParserException if the data is invalid:

julia> parse(PAFRecord, "not a PAF line")
ERROR: Error when parsing PAF record on line 1, near byte number 14 in line: Not enough tab-separated fields. Each line must have at least 12 fields
[...]

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.

The public, unexported try_parse function can be used instead of Base.parse. It either returns a valid PAFRecord, or else returns the ParserException, but does not throw the exception:

julia> const PAF = PairwiseMappingFormat;

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 try_next! function. This function will either:

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.

reader = PAFReader(open(path_to_paf))

# Emits `PAFRecord` for valid records
records = [try_next!(reader) for i in 1:3]
@assert records isa Vector{PAFRecord}

# Indefinitely emits `nothing` once the reader is empty
nothings = [try_next!(reader) for i in 1:10]
@assert nothings isa Vector{Nothing}

close(reader)

reader = PAFReader(IOBuffer("not a PAF record\n" * open(readline, path_to_paf)))
err = try_next!(reader)
@assert err isa PAF.ParserException

# The invalid line was skipped, so the next record can be read
@assert try_next!(reader) isa PAFRecord
@assert try_next!(reader) === nothing

The ParserException type contains the error type as an Enum, and the line where the exception occurred. These can be obtained with the .kind and .line properties.

julia> err.line
1

julia> err.kind
TooFewFields::Err = 0

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

PAFRecords may contain extra auxiliary fields at the end of the records, similar to SAM records. Any auxiliary data is stored in the PAFRecord, and lazily parsed and validated as they are accessed. This means that a PAFRecord may contain invalid auxiliary data.

The aux_data function constructs a lazy SAM.Auxiliary object from the package XAMAuxData.jl. The precise semantics of this Auxiliary type is a little complicated, and can be found in the documentation of XAMAuxData.jl. Here is a small example:

julia> line = first(eachline(path_to_paf));

julia> aux_string = join(split(line, '\t')[13:end], '\t')
"tp:A:P\tcm:i:29990\tdv:f:0.0000\tkm:Z:some text here"

julia> record = parse(PAFRecord, line);

julia> aux = aux_data(record);

julia> isvalid(aux)
true

julia> aux["tp"]
'P': ASCII/Unicode U+0050 (category Lu: Letter, uppercase)

julia> aux["km"]
"some text here"

julia> haskey(aux, "cm")
true

Misc info

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

PairwiseMappingFormat.PairwiseMappingFormat — Module

module PairwiseMappingFormat

Parse files in the Pairwise mApping Format (PAF).

PairwiseMappingFormat.PAFReader — Type

PAFReader(io::Union{IO, AbstractBufReader}; copy::Bool=true)

Construct a PAFReader, an iterator of PAFRecord read from io. Iterating this object returns a PAFRecord for each valid line of input, and throws a ParserException when an invalid record is found.

If copy is false, the records returned by iterating the reader, or by try_next!, 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: PAFReader

Extended help

The following properties may be used. For unmapped records, only qname and qlen are available, and accessing any other property throws an UnmappedError. Use is_mapped to check if a record is mapped.

qname::StringView. The query name. May be empty, and contain any bytes except \t and \n.
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 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::Int64 gives the length of the alignment
mapq::Union{Int, Nothing} gives the mapping quality, or nothing if this information is unavailable.
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)).

PairwiseMappingFormat.aux_data — Function

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. For more details about this object, read the documentation of the package XAMAuxData.jl.

Examples

julia> aux = aux_data(record);

julia> length(aux)
4

julia> aux["cm"]
29990

julia> haskey(aux, "k1")
false

PairwiseMappingFormat.aln_identity — Function

aln_identity(record::PAFRecord)::Float64

Return the nucleotide / residue identity, defined as the number of matches divided by the alignment length. 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 if the record is unmapped.

Examples

julia> aln_identity(record)
1.0

PairwiseMappingFormat.query_coverage — Function

query_coverage(record::PAFRecord)::Float64

Compute the fraction of the query covered by the alignment. Throws an UnmappedError if the record is unmapped.

Examples

julia> query_coverage(record)
0.9999535124653004

PairwiseMappingFormat.target_coverage — Function

target_coverage(record::PAFRecord)::Float64

Compute the fraction of the target covered by the alignment. Throws an UnmappedError if the record is unmapped.

Examples

julia> target_coverage(record)
0.044934629307437725

PairwiseMappingFormat.try_next! — Function

try_next!(reader::PAFReader)::Union{Nothing, PAFRecord, ParserException}

Try to read a record from the PAFReader 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 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 and advance the reader

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 on the reader. See PAFReader.

Examples

julia> reader = PAFReader(open(path_to_paf));

julia> try_next!(reader) isa PAFRecord
true

julia> [try_next!(reader) for i in 1:2] isa Vector{PAFRecord}
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\n" * open(readline, path_to_paf)));

julia> try_next!(reader) isa PAF.ParserException
true

julia> try_next!(reader) isa PAFRecord
true

julia> try_next!(reader) === nothing
true

Base.isempty — Method

isempty(reader::PAFReader)::Bool

Check whether reader has no more records, i.e. whether try_next! 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.

Base.peek — Method

peek(reader::PAFReader)::Union{Nothing, PAFRecord, ParserException}

Return what try_next! 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

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)

PairwiseMappingFormat.is_mapped — Function

is_mapped(x::PAFRecord)::Bool

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, target_coverage or aln_identity on it, throws an UnmappedError. Note that some PAF parsers do not handle unmapped records correctly, so be wary when writing unmapped records.

Examples

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