Auditing GFA traversals

Introduction

A pangenome graph usually claims to embed the assemblies it was built from: each named P path or W walk should spell out, base for base, the contig it is labelled with. That claim is easy to break and easy to miss. Graph construction and normalization rewrite segments, chunked builds are stitched back together, coordinates are lifted between references, and files are round-tripped through several tools before anyone looks at them. The graph stays well formed throughout, so structural validators stay quiet.

panpathAudit answers a narrower question: is the source sequence still recoverable from the traversal, and if not, exactly where does it stop matching? It reconstructs the sequence spelled by every named traversal, compares it against the source FASTA, and reports the first divergent base with its source, traversal and segment coordinates, along with digests and optional base-level alignment statistics.

Tools such as gfatools and odgi provide broader graph manipulation and analysis. panpathAudit instead provides a focused, reproducible audit of sequence preservation. It does not validate graph topology: links, overlaps beyond blunt joins, and connectivity are outside its scope. Results and diagnostics use Bioconductor DataFrame containers for integration with R workflows, while the bounded native implementation used by the command-line reference tool is retained. Source code and issues are available from the package repository.

Installation

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("panpathAudit")

You can install the development version of panpathAudit:

devtools::install_github("gkanogiannis/panpathAudit")

Building the package requires a C++20 compiler and zlib.

library(panpathAudit)

Input scope

Inputs may be plain or gzip-compressed FASTA and GFA files; compression is detected from file content, not from the file name. The package reads GFA P paths and W walks. Path overlaps must be blunt (* or 0M), and sequence symbols must use the IUPAC nucleotide alphabet. Graph topology is not validated.

Basic audit

fasta <- system.file("extdata", "example.fa", package = "panpathAudit")
gfa <- system.file("extdata", "example.gfa", package = "panpathAudit")
result <- auditGFA(fasta, gfa)
result
#> PanPathAuditResult
#>   identical: 1  divergent: 0 
#>   missing path: 0  missing source: 0 
#>   topology validated: no

Printing the object gives the headline counts. The per-traversal detail lives in the outcomes table.

auditOutcomes(result)[, c("identifier", "status", "record_type")]
#> DataFrame with 1 row and 3 columns
#>    identifier      status record_type
#>   <character> <character> <character>
#> 1      sample   IDENTICAL           P

The result separates traversal outcomes from summary, provenance and telemetry data. These accessors return SimpleList objects, so wrap them in as.list() when you want their contents printed.

as.list(auditSummary(result))
#> $counts
#>      identical      divergent   missing_path missing_source 
#>              1              0              0              0 
#> 
#> $percentages
#>      identical      divergent   missing_path missing_source 
#>            100              0              0              0 
#> 
#> $total
#> [1] 1
provenance <- auditProvenance(result)
as.list(provenance[c("gfa_version", "path_records", "walk_records")])
#> $gfa_version
#> [1] "1.1"
#> 
#> $path_records
#> [1] 1
#> 
#> $walk_records
#> [1] 0

Provenance also records the resolved input paths, the mapping and resource options in force, and a topology_validated flag that is always FALSE — a reminder of what the audit does not cover.

Audit statuses

Every source record and every named traversal is accounted for by exactly one of four statuses. All four are completed outcomes, not errors.

Status Meaning
IDENTICAL The traversal spells the source sequence exactly.
DIVERGENT Both sides exist but differ; see divergence columns.
MISSING_PATH A source record has no traversal of that name.
MISSING_SOURCE A traversal has no source record of that name.

A single audit can report all four at once. The following graph embeds one contig faithfully, one with a substitution, carries a path with no source, and leaves one source record unrepresented.

workdir <- tempfile("panpathAudit-demo-")
dir.create(workdir)

mixed_fasta <- file.path(workdir, "mixed.fa")
mixed_gfa <- file.path(workdir, "mixed.gfa")

writeLines(c(
    ">same", "ACGT",
    ">changed", "AAAA",
    ">absent", "NN"
), mixed_fasta)

writeLines(c(
    "H\tVN:Z:1.1",
    "S\t1\tACGT",
    "S\t2\tAATA",
    "S\t3\tGG",
    "P\tsame\t1+\t*",
    "P\tchanged\t2+\t*",
    "P\textra\t3+\t*"
), mixed_gfa)
mixed <- auditGFA(mixed_fasta, mixed_gfa)
auditOutcomes(mixed)[, c(
    "identifier", "status", "record_type", "source_length", "graph_length"
)]
#> DataFrame with 4 rows and 5 columns
#>    identifier         status record_type source_length graph_length
#>   <character>    <character> <character>     <numeric>    <numeric>
#> 1      absent   MISSING_PATH          NA             2           NA
#> 2     changed      DIVERGENT           P             4            4
#> 3       extra MISSING_SOURCE           P            NA            2
#> 4        same      IDENTICAL           P             4            4
auditSummary(mixed)$counts
#>      identical      divergent   missing_path missing_source 
#>              1              1              1              1

source_length is NA for MISSING_SOURCE and graph_length is NA for MISSING_PATH; columns that do not apply to a row are always NA. See ?PanPathAuditResult for the full column catalogue.

Reading a divergence

The point of the audit is the DIVERGENT row. Consider a source contig ACGA against a two-segment path that spells AC + GT:

divergent_fasta <- file.path(workdir, "divergent.fa")
divergent_gfa <- file.path(workdir, "divergent.gfa")

writeLines(c(">sample", "ACGA"), divergent_fasta)
writeLines(c(
    "H\tVN:Z:1.1",
    "S\t1\tAC",
    "S\t2\tGT",
    "P\tsample\t1+,2+\t*"
), divergent_gfa)
divergent <- auditGFA(divergent_fasta, divergent_gfa)
auditOutcomes(divergent)[, c(
    "status", "divergence_kind", "source_position_1based", "segment",
    "orientation", "traversal_position_1based", "segment_position_1based"
)]
#> DataFrame with 1 row and 7 columns
#>        status divergence_kind source_position_1based     segment orientation
#>   <character>     <character>              <numeric> <character> <character>
#> 1   DIVERGENT            BASE                      4           2           +
#>   traversal_position_1based segment_position_1based
#>                   <numeric>               <numeric>
#> 1                         2                       2

The source reads ACGA and the path spells ACGT, so the first three bases agree and the fourth does not. The report locates that base three ways, all of them one-based, as the _1based suffixes indicate: it is base 4 of the source contig, it falls in segment 2 (GT), and within that segment it is base 2 — the T the graph supplies where the source has A. traversal_position_1based locates that same divergence along the segment as the path reads it. This step is +, so the path reads segment 2 forwards and the two coordinates coincide. On a reverse step they do not, which is why both are reported; Reverse traversals below shows that case.

divergence_kind distinguishes three failure shapes:

  • BASE — both sides supply a base at that position and they differ.
  • END — one side ran out first.
  • RANGE_LENGTH — a walk reconstructed to a different length than its declared range.

When the graph is exhausted before the source is, no segment can be named:

truncated_fasta <- file.path(workdir, "truncated.fa")
truncated_gfa <- file.path(workdir, "truncated.gfa")

writeLines(c(">sample", "ACGT"), truncated_fasta)
writeLines(c("H\tVN:Z:1.1", "S\t1\tACG", "P\tsample\t1+\t*"), truncated_gfa)

auditOutcomes(auditGFA(truncated_fasta, truncated_gfa))[, c(
    "divergence_kind", "source_position_1based", "source_length",
    "graph_length", "segment"
)]
#> DataFrame with 1 row and 5 columns
#>   divergence_kind source_position_1based source_length graph_length     segment
#>       <character>              <numeric>     <numeric>    <numeric> <character>
#> 1             END                      4             4            3          NA

Reverse traversals

When a path traverses a segment on the reverse strand, the two position columns come apart, which is why both are reported. Here the divergence is at source base 5, inside the reverse-complemented segment 2 (TTGA, four bases long):

reverse_fasta <- file.path(workdir, "reverse.fa")
reverse_gfa <- file.path(workdir, "reverse.gfa")

writeLines(c(">sampleB#1#chr1", "ACGTGCAACC"), reverse_fasta)
writeLines(c(
    "H\tVN:Z:1.0",
    "S\t1\tACGT", "S\t2\tTTGA", "S\t3\tCC",
    "P\tsampleB#1#chr1\t1+,2-,3+\t*"
), reverse_gfa)

reverse <- auditGFA(reverse_fasta, reverse_gfa)
auditOutcomes(reverse)[, c(
    "source_position_1based", "segment", "orientation",
    "traversal_position_1based", "segment_position_1based"
)]
#> DataFrame with 1 row and 5 columns
#>   source_position_1based     segment orientation traversal_position_1based
#>                <numeric> <character> <character>                 <numeric>
#> 1                      5           2           -                         1
#>   segment_position_1based
#>                 <numeric>
#> 1                       4

traversal_position_1based is 1 — the divergence is at the first base the path reads out of that segment — while segment_position_1based is 4, because reading TTGA backwards starts at its last base. Use the traversal coordinate to reason about the path, and the segment coordinate to find the base in the S record.

Anatomy of a reverse-strand divergence. The highlighted column and both position rulers are taken from the audit result above; the segment blocks are the graph defined in this section.

Anatomy of a reverse-strand divergence. The highlighted column and both position rulers are taken from the audit result above; the segment blocks are the graph defined in this section.

Digests

Every row also carries source_sha256, graph_sha256 and their BLAKE3 counterparts (plus addressed_* variants for walks). They let you confirm an IDENTICAL verdict independently, or compare an audit against another tool’s output, without re-reading either file.

as.data.frame(auditOutcomes(divergent)[, c("source_sha256", "graph_sha256")])
#>                                                      source_sha256
#> 1 66b270c0d57529c39a5f52293606adb9a193b2c4e71137c32612b1df8fb695b1
#>                                                       graph_sha256
#> 1 1dff3e84fe7877e0673b69bbddcf40124e396e3f9943dd890c91b6a09adb9af0

Walks and mappings

W records identify a sequence as sample#haplotype#sequence and address a half-open, zero-based range of it. Assembly FASTA headers are rarely written in that form, so without help the two sides do not meet: the same file audited against a walk graph reports the source as MISSING_PATH and the walk as MISSING_SOURCE.

walk_gfa <- system.file("extdata", "walk.gfa", package = "panpathAudit")

plain_fasta <- file.path(workdir, "chr1.fa")
writeLines(c(">chr1", "ACGTGT"), plain_fasta)

auditOutcomes(auditGFA(plain_fasta, walk_gfa))[, c(
    "identifier", "status", "record_type"
)]
#> DataFrame with 2 rows and 3 columns
#>      identifier         status record_type
#>     <character>    <character> <character>
#> 1          chr1   MISSING_PATH          NA
#> 2 sample#1#chr1 MISSING_SOURCE           W

A mapping supplies the missing naming convention, with one row per FASTA file.

Mode Required fields Identifier
exact none FASTA record name
pansn sample, numeric haplotype sample#haplotype#name
prefix prefix prefix plus record name
mapping <- data.frame(
    mode = "pansn", sample = "sample", haplotype = "1", prefix = ""
)
mapped <- auditGFA(plain_fasta, walk_gfa, mapping = mapping)
auditOutcomes(mapped)[, c("identifier", "status", "record_type")]
#> DataFrame with 1 row and 3 columns
#>      identifier      status record_type
#>     <character> <character> <character>
#> 1 sample#1#chr1   IDENTICAL           W

Because a walk addresses only part of its source, a walk outcome reports both the range it covers and how much of the contig that range actually reaches.

walk_fasta <- system.file("extdata", "walk.fa", package = "panpathAudit")
walk_result <- auditGFA(walk_fasta, walk_gfa)
auditOutcomes(walk_result)[, c(
    "identifier", "status", "source_start_0based",
    "source_end_0based_exclusive", "addressed_source_bases", "source_length"
)]
#> DataFrame with 1 row and 6 columns
#>      identifier      status source_start_0based source_end_0based_exclusive
#>     <character> <character>           <numeric>                   <numeric>
#> 1 sample#1#chr1   IDENTICAL                   0                           6
#>   addressed_source_bases source_length
#>                <numeric>     <numeric>
#> 1                      4             6

Mind the two conventions in that one row. source_start_0based and source_end_0based_exclusive retain GFA’s zero-based, half-open ranges, so 0 and 6 span the whole six-base contig. The graph reaches only four of those bases: this example holds two walks, covering [0, 2) and [4, 6), leaving the two bases in between unaddressed. Every other displayed position, including all the divergence coordinates above, is one-based.

Walk coverage and the two coordinate conventions. Edge numbers are zero-based and half-open; base numbers are one-based. Ranges and totals are taken from the outcome row above.

Walk coverage and the two coordinate conventions. Edge numbers are zero-based and half-open; base numbers are one-based. Ranges and totals are taken from the outcome row above.

Comprehensive statistics

statistics = "comprehensive" adds a base-level alignment of every divergent pair and aggregates the result globally, by source and by record type.

detailed <- auditGFA(mixed_fasta, mixed_gfa, statistics = "comprehensive")
names(auditStatistics(detailed))
#> [1] "global"         "by_source"      "by_record_type"

global <- auditStatistics(detailed)$global
global$bases
#> $matched
#> [1] 7
#> 
#> $substituted
#> [1] 1
#> 
#> $source_only
#> [1] 0
#> 
#> $graph_only
#> [1] 0
#> 
#> $not_embedded
#> [1] 0
#> 
#> $missing_path
#> [1] 2
#> 
#> $missing_source
#> [1] 2
#> 
#> $unaligned_divergent_source
#> [1] 0
#> 
#> $unaligned_divergent_graph
#> [1] 0

The buckets account for every base the audit saw, on both sides. matched and substituted count aligned columns, each pairing one source base with one graph base; the remaining buckets count bases on one side only.

  • matched, substituted — aligned columns that agree or disagree.
  • source_only, graph_only — insertions and deletions in the alignment.
  • not_embedded — source bases no walk range addresses.
  • missing_path, missing_source — bases in records with no counterpart.
  • unaligned_divergent_* — divergent bases left unaligned by a resource limit.

Reading the example: of the 10 source bases, 7 align and match, 1 aligns and is substituted, and 2 sit in absent, which has no path. Of the 10 graph bases, the same 8 align and 2 sit in extra, which has no source. Per-row alignment detail sits alongside in the outcomes table.

How the comprehensive base buckets partition each side. Both bars are drawn from the values printed above and sum to the reported source and graph base totals.

How the comprehensive base buckets partition each side. Both bars are drawn from the values printed above and sum to the reported source and graph base totals.

auditOutcomes(detailed)[, c(
    "identifier", "status", "matched", "substituted", "edit_distance",
    "alignment_status"
)]
#> DataFrame with 4 rows and 6 columns
#>    identifier         status   matched substituted edit_distance
#>   <character>    <character> <numeric>   <numeric>     <numeric>
#> 1      absent   MISSING_PATH         0           0            NA
#> 2     changed      DIVERGENT         3           1             1
#> 3       extra MISSING_SOURCE         0           0            NA
#> 4        same      IDENTICAL         4           0             0
#>   alignment_status
#>        <character>
#> 1   not_applicable
#> 2         complete
#> 3   not_applicable
#> 4     not_required

alignment_status records why a row does or does not carry alignment numbers: not_required for an identical pair, complete for a finished alignment, not_applicable where one side is absent, and limit_exceeded or memory_limit_exceeded where a bound stopped the work. Alignment is bounded by alignmentMaxCells and memoryMiB; an unavailable alignment remains a completed outcome and is reported in the unaligned_divergent_* buckets rather than failing the audit. auditTelemetry() reports what the run actually used.

as.list(auditTelemetry(detailed))
#> $limit_bytes
#> [1] 1073741824
#> 
#> $peak_tracked_bytes
#> [1] 10
#> 
#> $oversized_contigs
#> [1] 0
#> 
#> $threads
#> [1] 1

Set threads to audit records in parallel; results are deterministic and do not depend on the thread count.

Diagnostics

Argument, input, and operational errors raise panpathAudit_argument_error, panpathAudit_input_error, and panpathAudit_operational_error respectively. All three inherit from panpathAudit_error, so one handler catches any of them, and each carries a diagnostics DataFrame rather than only a message.

diagnostics <- tryCatch(
    auditGFA(file.path(tempdir(), "missing.fa"), gfa),
    panpathAudit_error = function(error) error$diagnostics
)
diagnostics[, c("category", "code", "message")]
#> DataFrame with 1 row and 3 columns
#>      category        code           message
#>   <character> <character>       <character>
#> 1      source  FASTA_READ cannot open input

Alongside category, code and message, the table carries path, line, identifier, position and symbol wherever the failure can be pinned to a location; unrelated fields are NA. Diagnostic positions are one-based. A malformed input reports every problem it found, not just the first, so one call is enough to see the full picture.

Citation

citation("panpathAudit")
#> To cite panpathAudit:
#> 
#>   Gkanogiannis A (2026). "Panpath-Audit: Sequence-Level Verification of
#>   Embedded Pangenome Paths Against Source Assemblies." _Preprints_,
#>   *202608.1872*. doi:10.20944/preprints202608.1872.v1
#>   <https://doi.org/10.20944/preprints202608.1872.v1>. Version 1,
#>   <https://doi.org/10.20944/preprints202608.1872.v1>.
#> 
#> A BibTeX entry for LaTeX users is
#> 
#>   @Article{,
#>     title = {Panpath-Audit: Sequence-Level Verification of Embedded Pangenome Paths Against Source Assemblies},
#>     author = {Anestis Gkanogiannis},
#>     year = {2026},
#>     journal = {Preprints},
#>     volume = {202608.1872},
#>     note = {Version 1},
#>     doi = {10.20944/preprints202608.1872.v1},
#>     url = {https://doi.org/10.20944/preprints202608.1872.v1},
#>   }

References

Session information

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] panpathAudit_0.99.2 BiocStyle_2.41.0   
#> 
#> loaded via a namespace (and not attached):
#>  [1] cli_3.6.6            knitr_1.52           rlang_1.3.0         
#>  [4] xfun_0.61            otel_0.2.0           generics_0.1.4      
#>  [7] jsonlite_2.0.0       S4Vectors_0.51.10    buildtools_1.0.0    
#> [10] htmltools_0.5.9      maketools_1.3.2      sys_3.4.3           
#> [13] sass_0.4.10          stats4_4.6.1         rmarkdown_2.32      
#> [16] evaluate_1.0.5       jquerylib_0.1.4      fastmap_1.2.0       
#> [19] yaml_2.3.12          lifecycle_1.0.5      BiocManager_1.30.27 
#> [22] compiler_4.6.1       Rcpp_1.1.2           digest_0.6.39       
#> [25] R6_2.6.1             bslib_0.12.0         tools_4.6.1         
#> [28] BiocGenerics_0.59.12 cachem_1.1.0