diff --git a/CITATION.cff b/CITATION.cff index 3b8efca..9cb4ef4 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -2,7 +2,7 @@ cff-version: 1.2.0 message: "If you use VCF-RDFizer in your research, please cite it using the metadata below." title: "VCF-RDFizer" type: software -version: "3.0.2" +version: "3.0.3" authors: - name: "VCF-RDFizer maintainers" repository-code: "https://github.com/ecrum19/VCF-RDFizer" diff --git a/README.md b/README.md index b93c6ff..a43a104 100644 --- a/README.md +++ b/README.md @@ -347,9 +347,12 @@ preserves the per-sample vocabulary model: - every represented FORMAT slot is a `vcfc:FormatFieldValue`, linked to its declaration with `vcfc:declaredBy`; - `GT` is parsed into a `vcfc:Genotype` with an explicit `vcfc:phasingStatus` - and one `vcfc:GenotypeAlleleCall` per position, and `FT`, `PS`/`PSL`, `LAA`, - the copy-number and haplotype keys, and the `M`/`DPM`/`ADM` base-modification - families get their own resources; + — `vcfc:MixedPhasing` when its positions disagree — and one + `vcfc:GenotypeAlleleCall` per position carrying its own `vcfc:phaseIndicator`, + and `FT`, `PS`/`PSL`, `LAA` (with per-sample `vcfc:LocalAlleleMembership` + ordinals), the copy-number and haplotype keys, and the `M`/`DPM`/`ADM` + base-modification families in both their numeric and aliased spellings get + their own resources; - the VCF file declares `vcfc:representationProfile vcfc:ExpandedRepresentation`. With the default rules, these triples are appended directly from `records.tsv`; @@ -1080,7 +1083,7 @@ Safe termination: If you use VCF-RDFizer in a publication, please cite: -VCF-RDFizer maintainers. (2026). *VCF-RDFizer* (Version 3.0.2) [Computer software]. GitHub. https://github.com/ecrum19/VCF-RDFizer +VCF-RDFizer maintainers. (2026). *VCF-RDFizer* (Version 3.0.3) [Computer software]. GitHub. https://github.com/ecrum19/VCF-RDFizer BibTeX: @@ -1089,7 +1092,7 @@ BibTeX: author = {{VCF-RDFizer maintainers}}, title = {VCF-RDFizer}, year = {2026}, - version = {3.0.2}, + version = {3.0.3}, url = {https://github.com/ecrum19/VCF-RDFizer}, note = {Computer software} } diff --git a/conda-recipe/README.md b/conda-recipe/README.md index 51f7b4e..7dc88c7 100644 --- a/conda-recipe/README.md +++ b/conda-recipe/README.md @@ -7,11 +7,11 @@ do not submit this package to `staged-recipes`. ## Before submitting to conda-forge -1. Commit the version bump, then create and push a Git tag (for example `v3.0.2`). +1. Commit the version bump, then create and push a Git tag (for example `v3.0.3`). 2. Download the source tarball and compute sha256: ```bash curl -L -o vcf-rdfizer.tar.gz \ - https://github.com/ecrum19/VCF-RDFizer/archive/refs/tags/v3.0.2.tar.gz + https://github.com/ecrum19/VCF-RDFizer/archive/refs/tags/v3.0.3.tar.gz shasum -a 256 vcf-rdfizer.tar.gz ``` 3. Replace `version` and `sha256` in the feedstock's `recipe/meta.yaml`. diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index 7c41f55..de2ce94 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -1,5 +1,5 @@ {% set name = "vcf-rdfizer" %} -{% set version = "3.0.2" %} +{% set version = "3.0.3" %} package: name: {{ name|lower }} @@ -7,7 +7,7 @@ package: source: url: https://github.com/ecrum19/VCF-RDFizer/archive/refs/tags/v{{ version }}.tar.gz - sha256: 46c08df157269ec1e3b1c4666a20e37a7b1709b9aa8366265a2bae91d10ae637 + sha256: REPLACE_WITH_GITHUB_TARBALL_SHA256 build: noarch: python diff --git a/docs/architecture.md b/docs/architecture.md index 38df6e1..604940b 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -89,9 +89,9 @@ whose size is the product of its dimensions. | Emitter | Emits | Why not RML | | --- | --- | --- | -| `emit_sample_representation` | Expanded: `SampleCall`/`FormatFieldValue`/`Genotype`/`PhaseSet`/`LocalAlleleSet`/`BaseModification`. Condensed: `SampleSet`/`CohortCallMatrix`/`FormatValueVector`. `SampleSet` and `VCFSample` in both. | RML would first have to materialize a helper table of variants × samples (× FORMAT keys) — the largest intermediate the pipeline can produce | +| `emit_sample_representation` | Expanded: `SampleCall`/`FormatFieldValue`/`Genotype`/`PhaseSet`/`LocalAlleleSet`/`LocalAlleleMembership`/`BaseModification`. Condensed: `SampleSet`/`CohortCallMatrix`/`FormatValueVector`. `SampleSet` and `VCFSample` in both. | RML would first have to materialize a helper table of variants × samples (× FORMAT keys) — the largest intermediate the pipeline can produce | | `append_header_representation_rdf` | Header-line subclasses, the ordered `HeaderAttribute` resources every structured line needs, and the INFO/FORMAT/FILTER/ALT/contig/META/SAMPLE/PEDIGREE declaration properties | RML cannot choose a class per row, nor decompose one angle-bracketed cell into an ordered set of resources | -| `emit_record_detail` | `ID`, `ALT`, `QUAL`, `FILTER`, `infoRaw` (each may be the missing token); the allele layer (`ReferenceAllele`/`AltAllele`, breakends, symbolic types); structured `InfoFieldValue` with its `FieldValueItem` decomposition; and the SV carriers | RML cannot type the missing token per row, cannot switch a literal's datatype on a declared `Type`, and the allele and value-item layers would need variants × alleles and variants × INFO-keys helper tables | +| `emit_record_detail` | `ID`, `ALT`, `QUAL`, `FILTER`, `infoRaw` (each may be the missing token); their decompositions (`RecordIdentifier`, `filterStatus` with its `FilterCode` resources, the ordered `FormatKey` list); the allele layer (`ReferenceAllele`/`AltAllele`, breakends, symbolic types); structured `InfoFieldValue` with its `FieldValueItem` decomposition; and the SV carriers | RML cannot type the missing token per row, cannot switch a literal's datatype on a declared `Type`, and the allele, decomposition and value-item layers would need variants × alleles, variants × components and variants × INFO-keys helper tables | The exception count did not change with the move to the VCF Core vocabulary, but what each emitter covers grew considerably: the vocabulary models the VCF 4.5 diff --git a/docs/conversion.md b/docs/conversion.md index 2eda0bb..108b1b0 100644 --- a/docs/conversion.md +++ b/docs/conversion.md @@ -129,6 +129,23 @@ Float lexical space: | `.` or empty | `"."^^vcfc:Null` | | anything else | a plain string literal, deliberately kept rather than dropped, and reported by the SHACL layer | +Three of those cells are also **decomposed**, in both INFO representations, +because a single cell carries a list the vocabulary models as resources: + +- **ID.** Each semicolon-separated component becomes a + `vcfc:RecordIdentifier` at `…#record/{ROW_ID}/id/{N}`, linked with + `vcfc:hasIdentifier` and carrying `identifierValue` and `componentIndex`. + The raw `vcfc:recordId` stays beside it. +- **FILTER.** The call gets a `vcfc:filterStatus` of `vcfc:FiltersPassed` + (`PASS`), `vcfc:FiltersNotApplied` (`.`) or `vcfc:FiltersFailed`. In the + failed case each code also becomes a `vcfc:FilterCode` with + `filterCodeValue` and `componentIndex`, joined to its `##FILTER` declaration + with `vcfc:declaredByFilter` so a consumer can read the description without + re-reading the header. +- **FORMAT.** The key list becomes ordered `vcfc:FormatKey` resources carrying + `vcfc:fieldIndex` and `vcfc:declaredBy`, beside the raw `vcfc:formatRaw`. + The IRI shape follows the class's own `iriTemplate`. + `--info-representation structured` (the default) additionally emits: - **The allele layer.** `REF` and each `ALT` item become an ordered @@ -143,7 +160,8 @@ Float lexical space: `CHROM` names, with `vcfc:chromosome`. - **Structured INFO.** One `vcfc:InfoFieldValue` per record and key at `…#call/{ROW_ID}/info/{KEY}`, linked to the `##INFO` declaration through - `vcfc:declaredBy`. A single-valued field whose declared `Type` is `Integer` or + `vcfc:declaredBy` and carrying `vcfc:fieldIndex`, its position in the INFO + column. A single-valued field whose declared `Type` is `Integer` or `Float` also gets a typed `fieldValueInteger` / `fieldValueDecimal`; a `Flag` gets `vcfc:fieldValueBoolean true`. - **Value items.** A field whose declared `Number` gives its positions meaning — @@ -154,6 +172,11 @@ Float lexical space: …) record the width with `vcfc:tupleArity` so a consumer can regroup them. This closes the multi-valued-INFO gap that earlier versions recorded in `limitations.md`. + + The local-allele codes `LA`/`LR`/`LG` are a special case: they index a + *sample's* allele subset, and an INFO field has no sample. Their items are + therefore still materialized and indexed, but carry no `vcfc:forAllele` — + see the genotype layer below, where the subset exists. - **The SV carriers.** `SVLEN`, `SVCLAIM`, `IMPRECISE`, `NOVEL`, `END`, `EVENT` with `EVENTTYPE`, the confidence intervals (`CIPOS`/`CIEND` through FALDO `InRangePosition`, the rest as `vcfc:ConfidenceInterval`), and the @@ -224,11 +247,45 @@ instead of relying on repeated `sampleId` literals. **Expanded mode additionally parses the per-sample values** into the vocabulary's genotype layer: `GT` becomes a `vcfc:Genotype` with `genotypeString`, `ploidy`, `phasingStatus` and one `vcfc:GenotypeAlleleCall` per position (each either -pointing at the allele it called or flagged `isNoCall`); `FT` becomes -`sampleFilter`; `PS`/`PSL`/`PSO`/`PSQ` become a `vcfc:PhaseSet`; `LAA` becomes a -`vcfc:LocalAlleleSet`; `CN`/`CNQ`/`CNL`/`CNP` and `HAP`/`AHAP` become their -copy-number and haplotype properties; and the pattern-defined `M`/`DPM`/`ADM` -keys become `vcfc:BaseModification` resources pointing at a ChEBI residue. +pointing at the allele it called or flagged `isNoCall`, and each carrying its +own `vcfc:phaseIndicator`); `FT` becomes `sampleFilter`; `PS`/`PSL`/`PSO`/`PSQ` +become a `vcfc:PhaseSet`; `LAA` becomes a `vcfc:LocalAlleleSet`; +`CN`/`CNQ`/`CNL`/`CNP` and `HAP`/`AHAP` become their copy-number and haplotype +properties; and the pattern-defined `M`/`DPM`/`ADM` keys become +`vcfc:BaseModification` resources pointing at a ChEBI residue. + +Four details of that layer are worth stating explicitly: + +- **Local-allele fields are resolved through `LAA`, not positionally.** A + `Number=LA`/`LR` value list indexes the **sample's own local allele subset**, + not the record's ALT column. With `LAA=2,4` the local set is + {REF, ALT2, ALT4}, so the three values of an `LR` field belong to global + alleles 0, 2 and 4 — not 0, 1, 2. Each `vcfc:FieldValueItem` is joined with + `vcfc:forAllele` through that mapping. A sample that declares no `LAA` has no + local ordering to resolve against, so its items keep `vcfc:valueIndex` and get + no `vcfc:forAllele` at all: a wrong link is worse than an absent one, because + nothing in the graph would signal it. `LAA` is read before the per-key loop, + since a FORMAT string may legally list it after a field that depends on it. +- **Phasing is per position, not per genotype.** `vcfc:phasingStatus` is + `vcfc:Phased` or `vcfc:Unphased` only when every position agrees. A genotype + whose indicators disagree — `0|1/2` — gets `vcfc:MixedPhasing`, and the + per-call `vcfc:phaseIndicator` carries the precise reading. A single verdict + would misdescribe it. +- **Local-allele position belongs to the membership.** Alongside + `vcfc:hasLocalAllele`, each member of a `vcfc:LocalAlleleSet` gets a + `vcfc:LocalAlleleMembership` with `vcfc:localAllele` and `vcfc:localIndex`. + The ordinal cannot live on the allele resource, which is shared across + samples: two samples may give the same ALT different local positions. The + vocabulary deprecates `vcfc:localAlleleIndex` for that reason, and the + converter no longer emits it. +- **Base-modification aliases resolve.** VCF 4.5 reserves each modification key + in two spellings — `M[0-9]+[ACGTUN]` and a named alias such as `M5mC` for + `M27551C` — across all three families, 30 aliased keys in total. Both + spellings build the same carrier against the same ChEBI residue, while the + key keeps its written spelling in the resource IRI so the graph traces back + to the source column. A `Number=M` field is decomposed over the bases the + genotype's called alleles actually carry, in GT order, since that is the + order the specification defines it over. **Condensed mode deliberately stops at the vector.** Decomposing GT per sample is exactly the per-sample materialization the condensed profile exists to avoid. diff --git a/docs/limitations.md b/docs/limitations.md index f2d3f48..5402f53 100644 --- a/docs/limitations.md +++ b/docs/limitations.md @@ -80,7 +80,10 @@ declared assembly is recorded but not used for anything. target to the VCF Core vocabulary closed three gaps that used to be listed here: symbolic ALTs, breakends and `*` are now classified and parsed; `Number=A/R/G/P` values are decomposed into per-allele `vcfc:FieldValueItem` resources; and `GT` -is parsed into ordered allele calls with an explicit phasing status. See +is parsed into ordered allele calls with an explicit phasing status. The `ID`, +`FILTER` and `FORMAT` columns followed: each is now decomposed into ordered +resources beside its raw literal, and a genotype whose positions disagree about +phasing is `vcfc:MixedPhasing` rather than forced to a single verdict. See [`conversion.md`](conversion.md#4-stage-three--the-wrappers-own-emitters). What remains true: diff --git a/docs/rml-mappings.md b/docs/rml-mappings.md index 4c76062..66350fc 100644 --- a/docs/rml-mappings.md +++ b/docs/rml-mappings.md @@ -247,7 +247,7 @@ published vocabulary is worth knowing about: the ontology declares shape constrains them with `sh:datatype xsd:integer`, and a SHACL datatype constraint compares the datatype IRI exactly. This converter emits `xsd:integer`, which satisfies both. See -[`validation-migration-notes.md`](validation-migration-notes.md#3-ordinals-are-now-xsdinteger-not-xsdpositiveinteger). +[`validation-migration-notes.md`](validation-migration-notes.md#5-ordinals-are-xsdinteger). --- diff --git a/docs/roadmap.md b/docs/roadmap.md index 6b73f65..649efe3 100644 --- a/docs/roadmap.md +++ b/docs/roadmap.md @@ -33,7 +33,7 @@ a record with `ALT=.` satisfies both the shape and `vcfc:missingValuePolicy`. ### ~~1. The validation suite has not been migrated~~ — **done** -The suite passes: 415 tests, mutation score 96/113 (85%) across 60 mutations. +The suite passes: 707 tests, mutation score 96/113 (85%) across 60 mutations. The runner and the wrapper now share one vocabulary module instead of mirroring each other, and the census derives every new resource family from the VCF. Ten of the eighteen new mutations are recorded gaps, each with what would close it. diff --git a/docs/sample-representation-guide.md b/docs/sample-representation-guide.md index 60fc8dd..9473e94 100644 --- a/docs/sample-representation-guide.md +++ b/docs/sample-representation-guide.md @@ -231,8 +231,8 @@ a decoder or an application-level vector function. If you need to query genotypes directly, that is the argument for expanded mode. Everything *not* per-sample is identical in both profiles: the header layer, the -allele layer, structured INFO with its value items, and the SV carriers are all -emitted the same way. +allele layer, structured INFO with its value items, the SV carriers, and the +record-level ID/FILTER/FORMAT-key decompositions are all emitted the same way. ### Condensed growth diff --git a/docs/validation-methodology.md b/docs/validation-methodology.md index adc3cc5..d6b5e59 100644 --- a/docs/validation-methodology.md +++ b/docs/validation-methodology.md @@ -146,6 +146,19 @@ never declared it - making the oracle expect a `FilterDefinition` the graph could not contain. The oracle now reads the header block from the file itself, which is the same text the conversion reads. +The same normalisation reaches the data lines, and cost a second fix. The census +read each record from `str(variant)`, which restores a trailing FORMAT field the +file omitted: a sample written `0` under `GT:DP` comes back as `0:.`, so the +oracle counted a cell the file does not contain. That made the validator blind +to exactly the case VCF Core models on purpose - an emitter that *invented* the +dropped field would have looked correct. The census now reads the data lines +from the file's own text too, alongside the cyvcf2 reader, and asserts the two +stay aligned rather than silently miscounting if one skips a line. + +The rule generalises: **a parser's re-serialization is not the file.** Anything +the oracle counts should be read from the source text, not from what a library +renders back. + **When adding a query**, prefer datatype-*family* checks and lexical comparisons over anything that assumes a store's internal representation, and run the agreement script before trusting it. diff --git a/docs/validation-migration-notes.md b/docs/validation-migration-notes.md index 6f70518..21453f6 100644 --- a/docs/validation-migration-notes.md +++ b/docs/validation-migration-notes.md @@ -7,7 +7,7 @@ and follows each input's declared VCF version. The validation suite has been migrated with it. This file records what changed, what it found, and what is still open. -**Status: the suite passes.** 415 tests, no failures. The mutation score is +**Status: the suite passes.** 707 tests, no failures. The mutation score is 96/113 (85%) across 60 mutations, each listed with its expected detector. --- @@ -87,6 +87,33 @@ The fixture's `XSD_POSITIVE_INTEGER` is gone. VCF Core's datatype with the bound applied numerically, so the ontology range and the shapes no longer disagree. +### 6. The census counts the file's cells, not htslib's rendering of them + +`parse_vcf` read each data row from `str(variant)`. htslib normalises what it +re-serializes: a sample written `0` under `GT:DP` comes back as `0:.`, so the +census counted a FORMAT cell the file does not contain. That made the validator +blind to precisely the case VCF Core models deliberately — an emitter that +*invented* the dropped field would have looked correct. + +`read_vcf_data_columns` now reads the data lines from the file's own text, +alongside the cyvcf2 reader, the same way `read_vcf_header_text` already read +the header. The two iterators are asserted to stay aligned: a short read means +the reader skipped a line, and a silently wrong census is worse than an error. + +### 7. The record decompositions are derived, not mirrored + +`expected_census` gained the ID, FILTER and FORMAT-key decompositions, the +per-call `vcfc:phaseIndicator`, and `vcfc:fieldIndex` on both INFO entries and +FORMAT keys. `record_decomposition_counts` derives the ID and FILTER component +counts from the source columns — including the rule that `PASS` and the missing +token are statuses rather than failure codes, so they contribute no +`vcfc:FilterCode` — rather than asking the emitter what it produced. + +`typed_value_kind` now states the typed-value rule once for INFO and FORMAT +both. The previous INFO mirror used bare `int()`/`float()`, so it would have +counted a BCF-reserved integer and a non-finite float that the emitter +deliberately skips. + ## What the suite found Running it against the migrated emitters surfaced a real conversion defect that diff --git a/docs/validation.md b/docs/validation.md index 8787f7f..6eae925 100644 --- a/docs/validation.md +++ b/docs/validation.md @@ -97,6 +97,29 @@ vcf-rdfizer --mode full -i cohort.vcf.gz --validate \ --validation-query-timeout 14400 --validation-time-budget 86400 -o ./out ``` +### A startup probe that fails on a healthy engine + +Comunica's endpoint is probed twice at startup: `_await_bind` waits for the +port to answer at all, then `_await_warm` proves the endpoint can actually read +the source, so that readiness is charged to setup rather than billed to the +first query. + +The warm-up probe used to be a single request. Comunica answers the bind probe +and then hands the socket to its worker, so the next connection can be refused +for a moment while the process is perfectly healthy — and a refused connection +comes back from the kernel immediately, not after the warm-up budget. The result +was an engine that failed in milliseconds while reporting + +``` +comunica could not read /work/…/expanded.nt within 3600s +``` + +a timeout it had never waited for. The probe now retries until the budget is +genuinely spent, still failing at once if the server process exits or if the +endpoint *answers* with an error, and the duration in the message is measured +rather than read back from the flag. If you see that message now, the wait was +real: raise `--comunica-warmup-timeout`, or use `--engine qlever`. + ### A decode that reports success and writes garbage Every decode into N-Triples is now checked structurally before anything reads @@ -155,6 +178,15 @@ vcf-rdfizer --mode full \ For standalone validation, the input VCF and the RDF artifact must originate from the same conversion. +> **The validator ships inside the image.** `src/validation/` is copied into +> the Docker image at build time, so `--mode validation` runs the image's +> expectations, not the working tree's. A release that teaches the converter a +> new resource family teaches the runner to expect it in the same commit — but +> an image built before that commit will report the new predicates as +> unexpected. Pull the image matching your wrapper version (`docker pull +> ecrum19/vcf-rdfizer:latest`, or pin with `--image-version`) before reading a +> `MISMATCH` on Q9 or Q10 as a conversion defect. + ## Which artifact is validated `--rdf` accepts any artifact the pipeline produces: diff --git a/docs/vcf-coverage.md b/docs/vcf-coverage.md index 3feab54..b186869 100644 --- a/docs/vcf-coverage.md +++ b/docs/vcf-coverage.md @@ -40,13 +40,16 @@ rather than as silent confidence. Every gap is listed under | POS | `POS` | `vcfc:pos` (`xsd:integer`) | yes | `q01`, `q11`, `preflight_position_datatype` | `corrupt_pos`, `drop_pos`, `retype_pos_as_string` | | POS ↔ record binding | — | — | yes | `q11_record_digest` | `permute_pos` | | ID | `ID` | `vcfc:recordId` | yes | `q11_record_digest` | (covered by digest) | +| ID components | `ID` | `vcfc:hasIdentifier` → `vcfc:RecordIdentifier` (`identifierValue`, `componentIndex`) | yes | `q09_predicate_census`, `q10_class_census` | — | | REF | `REF` | `vcfc:ref` | yes | `q02`, `q03`, `q11` | `corrupt_alt` | | ALT | `ALT` | `vcfc:alt` | yes | `q02`, `q03`, `q11` | `corrupt_alt` | | REF/ALT ↔ record binding | — | — | yes | `q11_record_digest` | `permute_ref_alt` | | QUAL | `QUAL` | `vcfc:qual` (`xsd:decimal` / `vcfc:Null`) | yes | `q09`, `q11` | `drop_qual`, `drop_all_qual`, `corrupt_qual` | | FILTER | `FILTER` | `vcfc:filter` | yes | `q04` (exact lexical), `q11` | `drop_filter`, `corrupt_filter_lexical` | +| FILTER outcome | `FILTER` | `vcfc:filterStatus` → `vcfc:FiltersPassed` / `FiltersFailed` / `FiltersNotApplied` | yes | `q09_predicate_census`, `q10_class_census` | — | +| FILTER codes | `FILTER` | `vcfc:hasFilterCode` → `vcfc:FilterCode` (`filterCodeValue`, `componentIndex`, `declaredByFilter`) | yes | `q09_predicate_census`, `q10_class_census` | — | | INFO (raw) | `INFO` | `vcfc:infoRaw` | yes | `q11_record_digest` | `corrupt_info_raw` | -| INFO (structured) | `INFO` | `vcfc:hasInfoValue` → `vcfc:InfoFieldValue` → `vcfc:declaredBy` | yes | `q09`, `q12_info_value_digest` | `drop_info_value`, `corrupt_info_value`, `retype_info_value` | +| INFO (structured) | `INFO` | `vcfc:hasInfoValue` → `vcfc:InfoFieldValue` → `vcfc:declaredBy`, `vcfc:fieldIndex` | yes | `q09`, `q12_info_value_digest` | `drop_info_value`, `corrupt_info_value`, `retype_info_value` | INFO values carry `vcfc:fieldValue` plus a typed `fieldValueInteger` / `fieldValueDecimal` when the declaration says `Number=1` and the value parses; @@ -86,6 +89,7 @@ into ordered `vcfc:FieldValueItem` resources joined to their allele with | VCF element | RDF term | Represented | Validated by | Mutation | | --- | --- | --- | --- | --- | | FORMAT declaration | `vcfc:formatRaw` | yes | `q09_predicate_census` | — | +| FORMAT keys, in order | `vcfc:hasFormatKey` → `vcfc:FormatKey` (`fieldIndex`, `declaredBy`) | yes | `q09_predicate_census`, `q10_class_census` | — | | Sample identity | `vcfc:sampleId` / `vcfc:sampleName` | yes | `preflight_sample_gt_inventory`, `q09` | `drop_sample_call` | | GT values (lexical) | `vcfc:hasFormatValue` → `vcfc:fieldValue` | yes | `q05`, `q06`, `q13` | `flip_genotype` | | Non-GT FORMAT (DP, GQ, AD, PL…) | `vcfc:FormatFieldValue` / `vcfc:FormatValueVector` | yes | `q13_format_value_digest` | `drop_format_value_dp`, `corrupt_format_value_dp`, `corrupt_format_vector` | @@ -99,14 +103,16 @@ these values inside `vcfc:encodedValues`, which is the point of the profile. | VCF element | RDF term | Represented | Validated by | Mutation | | --- | --- | --- | --- | --- | | `GT` parsed | `vcfc:hasGenotype` → `vcfc:Genotype` (`genotypeString`, `ploidy`) | expanded only | — | — | -| Phasing (`\|` vs `/`) | `vcfc:phasingStatus` → `vcfc:Phased` / `vcfc:Unphased` | expanded only | — | — | +| Phasing (`\|` vs `/`) | `vcfc:phasingStatus` → `vcfc:Phased` / `vcfc:Unphased` / `vcfc:MixedPhasing` | expanded only | — | — | +| Per-position phase indicator | `vcfc:phaseIndicator` on each `vcfc:GenotypeAlleleCall` | expanded only | `q09_predicate_census` | — | | Each GT position | `vcfc:hasAlleleCall` → `vcfc:GenotypeAlleleCall` (`callIndex`, `isNoCall`, `calledAllele`) | expanded only | — | — | | `FT` | `vcfc:sampleFilter` | expanded only | — | — | | `PS`/`PSL`/`PSO`/`PSQ` | `vcfc:inPhaseSet` → `vcfc:PhaseSet` | expanded only | — | — | -| `LAA` | `vcfc:hasLocalAlleleSet` → `vcfc:LocalAlleleSet` | expanded only | — | — | +| `LAA` | `vcfc:hasLocalAlleleSet` → `vcfc:LocalAlleleSet` → `vcfc:hasLocalAllele` | expanded only | — | — | +| `LAA` membership order | `vcfc:hasLocalAlleleMembership` → `vcfc:LocalAlleleMembership` (`localAllele`, `localIndex`) | expanded only | — | — | | `CN`/`CNQ`/`CNL`/`CNP` | `vcfc:copyNumber` and friends | expanded only | — | — | | `HAP`/`AHAP` | `vcfc:haplotypeId`, `vcfc:ancestralHaplotypeId` | expanded only | — | — | -| `M*`/`DPM*`/`ADM*` | `vcfc:BaseModification` with a ChEBI `modifiedResidue` | expanded only | — | — | +| `M*`/`DPM*`/`ADM*`, numeric or aliased | `vcfc:BaseModification` with a ChEBI `modifiedResidue` | expanded only | — | — | ## Header section @@ -187,6 +193,13 @@ conformant graph. resources, the way `q11` already covers the record fields — or leaning on the vocabulary's own consistency SHACL profile, which checks several of these agreements directly. + + The ID, FILTER and FORMAT-key decompositions join this group rather than + forming a new one: `q09`/`q10` expect an exact count of each + `vcfc:RecordIdentifier`, `vcfc:FilterCode` and `vcfc:FormatKey`, derived from + the source columns, so dropping one or emitting a spurious one is caught — + but `identifierValue`, `filterCodeValue` and the component ordinals are not + compared against the file. No mutation targets them yet. 5. **Ordering is checked by shapes, not by queries.** `vcfc:lineIndex` and `vcfc:recordIndex` are counted, but `corrupt_record_index` is undetected by the query layer. The SPARQL SHACL profile enforces uniqueness, nondecreasing @@ -240,18 +253,26 @@ changes in the emitted graph. VCF Core 2.0.0 defines 338 terms. The converter emits from most of them, and its output is conformant, but the vocabulary can express more than the converter -currently derives. Ten families are entirely unemitted: +currently derives. + +Five of the ten families previously listed here as entirely unemitted are now +emitted. Two are complete; one term each remains in the other three: + +| Family | Emitted | Still unemitted | +| --- | --- | --- | +| Parsed FILTER | `filterStatus`, `FiltersPassed/Failed/NotApplied`, `FilterCode`, `hasFilterCode`, `filterCodeValue`, `declaredByFilter` | — | +| Local-allele membership | `LocalAlleleMembership`, `localAllele`, `localIndex` | — | +| Parsed FORMAT keys | `FormatKey`, `hasFormatKey`, `fieldIndex` | `keyPattern` | +| Parsed record IDs | `RecordIdentifier`, `hasIdentifier`, `identifierValue`, `componentIndex` | `aliasOf` — which identifier aliases which is not in the file | +| Mixed phasing | `MixedPhasing`, `phaseIndicator` | `allelePhaseSet` | + +Five families remain entirely unemitted: | Family | Terms | What it would add | | --- | --- | --- | -| Parsed FILTER | `filterStatus`, `FiltersPassed/Failed/NotApplied`, `FilterCode`, `hasFilterCode`, `declaredByFilter` | FILTER as resources joined to their `##FILTER` declarations, instead of only the raw string | | Padding semantics | `PaddingInterpretation`, `paddingRule`, `paddingSide`, `paddingAnchorPosition`, `paddingBaseCount` and their individuals | Which base of REF/ALT is the VCF padding base, and why | -| Parsed FORMAT keys | `FormatKey`, `hasFormatKey`, `fieldIndex`, `keyPattern` | The FORMAT key list as ordered resources, not only `formatRaw` | -| Parsed record IDs | `RecordIdentifier`, `hasIdentifier`, `identifierValue`, `aliasOf` | The semicolon-separated ID column split into individual identifiers | | Repeat units | `RepeatUnit`, `hasRepeatUnit`, `repeatUnitBases`, `rucConfidenceInterval`, `rbConfidenceInterval` | `RUB`-level detail below a repeat sequence | | Breakend mates | `mateBreakend`, `partnerBreakend`, `insertedSequence`, `isTelomereBreakend` | Resolving a breakend to the record describing its mate | -| Local-allele membership | `LocalAlleleMembership`, `localAllele`, `localIndex` | LAA membership with its own index, beside the current `LocalAlleleSet` | -| Mixed phasing | `MixedPhasing`, `phaseIndicator`, `allelePhaseSet` | Per-position phase indicators rather than one status per genotype | | Reference blocks | `hasReferenceBlock`, `blockAllele` | A record-level link to its gVCF block | | Raw/decoded values | `rawValue`, `decodedValue`, `percentEncodingPolicy`, `sampleFieldsRaw` | Percent-decoded companions beside the lossless source values | @@ -274,7 +295,7 @@ five version overlays — with **zero violations**, in both representation profiles, checked with pySHACL using the bundled ontology plus every `ontology/versions/*.ttl` and RDFS inference, which is the configuration the vocabulary's own `tests/validate_shacl.py` uses. See -[`validation-migration-notes.md`](validation-migration-notes.md#how-this-was-verified) +[`validation-migration-notes.md`](validation-migration-notes.md#how-conformance-was-verified) for exactly what was covered. The only results reported are `sh:Warning`s where an input VCF omits the diff --git a/pyproject.toml b/pyproject.toml index aa58357..df6a8d0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "vcf-rdfizer" -version = "3.0.2" +version = "3.0.3" description = "Docker-first VCF to RDF conversion targeting the VCF Core vocabulary, with compressed queryable representations (HDT, COTTAS), semantic validation, and data linking" readme = "README.md" requires-python = ">=3.10" diff --git a/src/validation/validation_runner.py b/src/validation/validation_runner.py index 5f9765c..3dfb94d 100644 --- a/src/validation/validation_runner.py +++ b/src/validation/validation_runner.py @@ -814,6 +814,35 @@ def expected_census( ): predicates[predicate] = 0 if parser[field] == METADATA_ABSENT else 1 + # ID and FILTER decompositions, and the record's declared FORMAT keys. + # All three are record-level, so both sample profiles carry them. + identifiers = parser.get("recordIdentifierCount", 0) + filter_codes = parser.get("filterCodeCount", 0) + format_keys = parser.get("formatKeyOccurrences", 0) + if identifiers: + classes[f"{VCFC}RecordIdentifier"] = identifiers + predicates[f"{VCFC}hasIdentifier"] = identifiers + predicates[f"{VCFC}identifierValue"] = identifiers + if filter_codes: + classes[f"{VCFC}FilterCode"] = filter_codes + predicates[f"{VCFC}hasFilterCode"] = filter_codes + predicates[f"{VCFC}filterCodeValue"] = filter_codes + predicates[f"{VCFC}declaredByFilter"] = parser.get("declaredFilterCodeCount", 0) + # componentIndex is shared by both decompositions. + if identifiers or filter_codes: + predicates[f"{VCFC}componentIndex"] = identifiers + filter_codes + # Every record gets a status, including PASS and the missing token. + predicates[f"{VCFC}filterStatus"] = records + if format_keys: + classes[f"{VCFC}FormatKey"] = format_keys + predicates[f"{VCFC}hasFormatKey"] = format_keys + predicates[f"{VCFC}fieldIndex"] = ( + predicates.get(f"{VCFC}fieldIndex", 0) + format_keys + ) + predicates[f"{VCFC}declaredBy"] = ( + predicates.get(f"{VCFC}declaredBy", 0) + format_keys + ) + # A sites-only VCF still declares a profile: vcfc:RepresentationProfileShape # requires exactly one on every file, and a file with no genotype columns # has no genotype data to condense. @@ -848,6 +877,11 @@ def expected_census( predicates[f"{VCFC}forSample"] = records * samples predicates[f"{VCFC}hasFormatValue"] = parser["formatValueSlots"] predicates[f"{VCFC}fieldValue"] = parser["nonEmptyFormatValues"] + # One preceding indicator per allele call, including the effective + # first one the source omits. + calls = parser.get("emittedGenotypePredicates", {}).get("callIndex", 0) + if calls: + predicates[f"{VCFC}phaseIndicator"] = calls # A single, non-missing Integer/Float FORMAT cell also gains the # typed companion property, exactly as an INFO value does. predicates[f"{VCFC}fieldValueInteger"] = ( @@ -912,6 +946,11 @@ def expected_census( predicates[f"{VCFC}fieldValue"] = ( predicates.get(f"{VCFC}fieldValue", 0) + values - parser["infoFlagCount"] ) + # Source order within the INFO column; the FORMAT keys add their own + # below, and both use the same predicate. + predicates[f"{VCFC}fieldIndex"] = ( + predicates.get(f"{VCFC}fieldIndex", 0) + values + ) # The allele layer, the value items, the SV carriers and the parsed # genotype layer travel with the structured INFO representation. @@ -1013,6 +1052,50 @@ def declared_field_types(raw_header: str, header_key: str = "INFO") -> dict[str, info_declared_types = declared_field_types +def record_decomposition_counts(columns, declared_filter_ids): + """Count the ID and FILTER components one record contributes. + + The emitter breaks both columns into ordered resources, so the census has to + count the same parts. PASS and the missing token are FILTER statuses rather + than failure codes and contribute no vcfc:FilterCode; a missing ID column + contributes no vcfc:RecordIdentifier. + + Returns (identifiers, filter codes, filter codes that cite a declaration). + """ + record_id = columns[2] if len(columns) > 2 else "." + identifiers = 0 + if record_id != ".": + identifiers = sum(1 for part in record_id.split(";") if part) + codes = declared = 0 + filter_column = columns[6] if len(columns) > 6 else "." + if filter_column not in ("PASS", "."): + for code in (part for part in filter_column.split(";") if part): + codes += 1 + if code in declared_filter_ids: + declared += 1 + return identifiers, codes, declared + + +def read_vcf_data_columns(vcf_path: Path): + """Yield each data line's columns, from the file's own text. + + Same reason as :func:`read_vcf_header_text`: htslib normalises what it + re-serializes. ``str(variant)`` restores a trailing FORMAT field the file + omitted -- a sample written ``0`` comes back as ``0:.`` -- so a census taken + from it counts cells the file does not contain. Worse, it makes the + validator blind to exactly the case VCF Core models deliberately: an emitter + that invented the dropped field would look correct. + """ + opener = gzip.open if vcf_path.name.endswith(".gz") else open + with opener(vcf_path, "rt", encoding="utf-8", errors="replace") as handle: + for line in handle: + if line.startswith("#"): + continue + stripped = line.rstrip("\r\n") + if stripped: + yield stripped.split("\t") + + def read_vcf_header_text(vcf_path: Path) -> str: """Read the header block straight from the VCF file. @@ -1195,12 +1278,15 @@ def parse_vcf( info_definitions: set[str] = set() info_values = info_flags = info_typed_integers = info_typed_decimals = 0 format_typed_integers = format_typed_decimals = 0 + declared_filter_ids = set(declared_field_types(_raw_header, "FILTER")) + record_identifiers = filter_codes = declared_filter_codes = 0 source_component = rml_uri_component(vcf_path.name) record_rows: list[list[str]] = [] records_with_format_column = records_with_format_keys = 0 format_key_occurrences = format_value_slots = non_empty_format_values = 0 distinct_format_keys: set[str] = set() _scan_started = time.monotonic() + source_columns = read_vcf_data_columns(vcf_path) try: for variant in reader: total_records += 1 @@ -1233,8 +1319,23 @@ def parse_vcf( if format_keys: records_with_format_column += 1 - columns = str(variant).rstrip("\r\n").split("\t") + # The file's own columns, not str(variant)'s: see + # read_vcf_data_columns. cyvcf2 yields records in file order, so the + # two iterators stay aligned; a short read means the reader skipped + # a line and the census would be silently wrong, so it is an error. + columns = next(source_columns, None) + if columns is None: + raise ValueError( + f"{vcf_path.name}: fewer data lines than parsed records at " + f"record {total_records}" + ) record_rows.append(columns) + _ids, _codes, _declared = record_decomposition_counts( + columns, declared_filter_ids + ) + record_identifiers += _ids + filter_codes += _codes + declared_filter_codes += _declared # The record digest is computed from the raw line so it matches the # lexical values the mapping puts in the graph, character for # character, with no round-trip through cyvcf2's typed accessors. @@ -1407,6 +1508,9 @@ def parse_vcf( "infoTypedDecimalCount": info_typed_decimals, "formatTypedIntegerCount": format_typed_integers, "formatTypedDecimalCount": format_typed_decimals, + "recordIdentifierCount": record_identifiers, + "filterCodeCount": filter_codes, + "declaredFilterCodeCount": declared_filter_codes, "q11_record_digest": [ {"bucket": bucket, "recordCount": int(count)} for bucket, count in sorted(digest_buckets.items()) @@ -2150,27 +2254,72 @@ def _await_bind(self, server_log: Path) -> None: f"{self.bind_timeout}s: {last_error}" ) + @staticmethod + def _is_connection_level_error(error: BaseException) -> bool: + """Is this the socket never opening, rather than the server answering? + + An HTTPError means the endpoint replied, so whatever it said is a real + result and must not be retried. Everything else at this layer -- a + refused, reset or aborted connection, or a socket that timed out -- is + about reaching the process, not about what it thinks of the query. + """ + import urllib.error + + if isinstance(error, urllib.error.HTTPError): + return False + if isinstance(error, urllib.error.URLError): + return isinstance(error.reason, (OSError, TimeoutError)) + return isinstance(error, (ConnectionError, TimeoutError)) + def _await_warm(self, server_log: Path) -> None: """Prove the endpoint can actually read the source, before the suite runs. Charged to setup rather than to the first query, so one query is not billed for readiness that every later query got for free. + + A refused connection is retried rather than fatal. ``_await_bind`` + returns as soon as the port answers one probe, but comunica hands the + socket to its worker after that reply, so the next connection can be + refused for a moment while the process is perfectly healthy. Treating + that as failure made the whole engine unusable on a timing accident -- + and reported it as a timeout that never happened, because a refused + connection comes back from the kernel at once rather than after the + warm-up budget. Only the server exiting, an answer from the endpoint, + or the budget genuinely running out ends this wait. """ assert self.endpoint is not None started = time.monotonic() - try: - self._post( - self.endpoint, "SELECT * WHERE { ?s ?p ?o } LIMIT 1", - timeout=self.warmup_timeout, - ) - except Exception as error: # noqa: BLE001 - surfaced as engine failure + deadline = started + self.warmup_timeout + last_error: BaseException | None = None + while True: self._assert_alive(server_log) - raise RuntimeError( - f"{self.name} could not read {self._endpoint_source_argument()} " - f"within {self.warmup_timeout}s (raise --comunica-warmup-timeout, " - f"or use --engine qlever): {error}" - ) from error - self.warmup_seconds = time.monotonic() - started + remaining = deadline - time.monotonic() + if remaining <= 0: + break + try: + self._post( + self.endpoint, "SELECT * WHERE { ?s ?p ?o } LIMIT 1", + timeout=remaining, + ) + except Exception as error: # noqa: BLE001 - surfaced as engine failure + self._assert_alive(server_log) + if not self._is_connection_level_error(error): + raise RuntimeError( + f"{self.name} could not read " + f"{self._endpoint_source_argument()}: {error}" + ) from error + last_error = error + time.sleep(min(0.5, max(0.0, deadline - time.monotonic()))) + else: + self.warmup_seconds = time.monotonic() - started + return + self._assert_alive(server_log) + raise RuntimeError( + f"{self.name} could not read {self._endpoint_source_argument()} " + f"within {time.monotonic() - started:.0f}s " + f"(raise --comunica-warmup-timeout, or use --engine qlever): " + f"{last_error or 'the endpoint never accepted a connection'}" + ) def _assert_alive(self, server_log: Path) -> None: if self.server is not None and self.server.poll() is not None: diff --git a/test/test_validation_engines_unit.py b/test/test_validation_engines_unit.py index 9a227da..d6c70d4 100644 --- a/test/test_validation_engines_unit.py +++ b/test/test_validation_engines_unit.py @@ -438,6 +438,141 @@ def fake_popen(command, stdout=None, **kwargs): engine.start() self.assertIn("index not found", str(caught.exception)) + def _comunica_engine_awaiting_warm(self, tmp_path, *, warmup_timeout): + """A started comunica engine, parked just before the warm-up probe.""" + source = tmp_path / "cohort.nt" + source.write_bytes(TRIPLES) + engine = V.build_engine( + "comunica", source, raw_dir=tmp_path, scratch=tmp_path, + options={"comunica_warmup_timeout": warmup_timeout}, + ) + engine.endpoint = f"http://127.0.0.1:{engine.port}/sparql" + engine.server = mock.MagicMock() + engine.server.poll.return_value = None + return engine + + def test_a_refused_connection_during_warmup_is_retried(self): + """The bind probe is not a promise that the next connection is accepted. + + comunica answers the bind probe and then hands the socket to its + worker, so a connection a moment later can be refused while the process + is perfectly healthy. That briefly-closed port used to fail the whole + engine - and the run it gated - on a timing accident. + """ + import urllib.error + + with tempfile.TemporaryDirectory() as td: + tmp_path = Path(td) + engine = self._comunica_engine_awaiting_warm(tmp_path, warmup_timeout=30) + attempts = [] + + def flaky_post(endpoint, query, *, timeout): + attempts.append(timeout) + if len(attempts) < 3: + raise urllib.error.URLError(ConnectionRefusedError(111, "refused")) + return b'{"results": {"bindings": []}}' + + with mock.patch.object(engine, "_post", flaky_post), \ + mock.patch.object(V.time, "sleep"): + engine._await_warm(tmp_path / "server.log") + + self.assertEqual(len(attempts), 3) + self.assertIsNotNone(engine.warmup_seconds) + # Each attempt is given what is left of the budget, never the whole + # of it again, so retrying cannot outlive --comunica-warmup-timeout. + self.assertLessEqual(attempts[-1], attempts[0]) + + def test_a_server_that_dies_during_warmup_fails_at_once(self): + """Retrying a refused port must not paper over a dead engine.""" + import urllib.error + + with tempfile.TemporaryDirectory() as td: + tmp_path = Path(td) + server_log = tmp_path / "server.log" + server_log.write_text("comunica: out of memory\n", encoding="utf-8") + engine = self._comunica_engine_awaiting_warm(tmp_path, warmup_timeout=3600) + engine.server.poll.return_value = 137 + engine.server.returncode = 137 + + def refused(endpoint, query, *, timeout): + raise urllib.error.URLError(ConnectionRefusedError(111, "refused")) + + started = time.monotonic() + with mock.patch.object(engine, "_post", refused): + with self.assertRaises(RuntimeError) as caught: + engine._await_warm(server_log) + + # Fails on the dead process, not after the 3600s budget. + self.assertLess(time.monotonic() - started, 30) + self.assertIn("137", str(caught.exception)) + self.assertIn("out of memory", str(caught.exception)) + + def test_an_answered_error_is_not_retried(self): + """An HTTP reply is the endpoint's verdict, not a connection problem. + + A source comunica cannot parse answers every probe the same way, so + retrying it would spin for the whole warm-up budget before reporting + what it already knew on the first attempt. + """ + import urllib.error + + with tempfile.TemporaryDirectory() as td: + tmp_path = Path(td) + engine = self._comunica_engine_awaiting_warm(tmp_path, warmup_timeout=3600) + attempts = [] + + def http_error(endpoint, query, *, timeout): + attempts.append(timeout) + raise urllib.error.HTTPError( + endpoint, 500, "Parse error", {}, None + ) + + with mock.patch.object(engine, "_post", http_error): + with self.assertRaises(RuntimeError) as caught: + engine._await_warm(tmp_path / "server.log") + + self.assertEqual(len(attempts), 1) + message = str(caught.exception) + self.assertIn("could not read", message) + # The old message claimed a timeout it had never waited for. + self.assertNotIn("3600s", message) + + def test_an_exhausted_warmup_budget_reports_the_time_it_waited(self): + """The duration in the message must be measured, not the flag's value. + + The original reported the configured budget whatever had happened, so a + connection refused in milliseconds was announced as an hour-long + timeout - pointing every reader at the wrong thing entirely. + """ + import itertools + import urllib.error + + with tempfile.TemporaryDirectory() as td: + tmp_path = Path(td) + engine = self._comunica_engine_awaiting_warm(tmp_path, warmup_timeout=30) + attempts = [] + + def refused(endpoint, query, *, timeout): + attempts.append(timeout) + raise urllib.error.URLError(ConnectionRefusedError(111, "refused")) + + # A clock that runs 10s per reading, so the budget expires after one + # attempt without the test waiting for any of it. + clock = itertools.count(0, 10) + with mock.patch.object(V.time, "monotonic", lambda: next(clock)), \ + mock.patch.object(V.time, "sleep"), \ + mock.patch.object(engine, "_post", refused): + with self.assertRaises(RuntimeError) as caught: + engine._await_warm(tmp_path / "server.log") + + self.assertEqual(len(attempts), 1) + message = str(caught.exception) + # 40s elapsed against a 30s budget: the figure is measured, and is + # not simply the flag read back out. + self.assertIn("within 40s", message) + self.assertIn("refused", message) + self.assertIn("--comunica-warmup-timeout", message) + def test_qlever_command_lines_are_overridable(self): """QLever's CLI varies by release, so both argv are user-overridable.""" with tempfile.TemporaryDirectory() as td: diff --git a/test/test_validation_oracle_unit.py b/test/test_validation_oracle_unit.py index 1b9cc70..374a927 100644 --- a/test/test_validation_oracle_unit.py +++ b/test/test_validation_oracle_unit.py @@ -192,6 +192,112 @@ def test_an_empty_allele_tuple_is_missing(self): self.assertEqual(V.classify_genotype((), has_gt=True), "MISSING") +class SourceDataColumnTests(VerboseTestCase): + """The census must see the file's cells, not htslib's rendering of them.""" + + def _write(self, body): + import tempfile, pathlib + path = pathlib.Path(tempfile.mkdtemp()) / "s.vcf" + path.write_text(body, encoding="utf-8") + return path + + def test_a_dropped_trailing_format_field_stays_dropped(self): + """VCF 4.5 lets a sample omit trailing FORMAT fields. + + cyvcf2 re-serializes that sample as '0:.', so a census taken from + str(variant) counts a cell the file does not contain -- and could never + catch an emitter that invented it. + """ + path = self._write( + "##fileformat=VCFv4.5\n" + '##FORMAT=\n' + '##FORMAT=\n' + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tS1\n" + "chr1\t1\t.\tA\tG\t.\t.\t.\tGT:DP\t0\n" + ) + columns = list(V.read_vcf_data_columns(path)) + self.assertEqual(len(columns), 1) + self.assertEqual(columns[0][9], "0") + + def test_header_lines_and_blank_lines_are_skipped(self): + """Only data lines reach the census.""" + path = self._write( + "##fileformat=VCFv4.5\n" + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n" + "chr1\t1\t.\tA\tG\t.\t.\t.\n" + "\n" + "chr1\t2\t.\tC\tT\t.\t.\t.\n" + ) + columns = list(V.read_vcf_data_columns(path)) + self.assertEqual([c[1] for c in columns], ["1", "2"]) + + +class TypedValueKindTests(VerboseTestCase): + """The rule deciding which cells gain a typed companion property. + + Stated independently of the emitter on purpose: an expectation computed by + the code under test would only prove it agrees with itself. + """ + + def test_a_single_integer_gains_an_integer_companion(self): + self.assertEqual(V.typed_value_kind("42", "Integer"), "Integer") + + def test_a_single_float_gains_a_decimal_companion(self): + self.assertEqual(V.typed_value_kind("0.5", "Float"), "Float") + + def test_a_string_field_gains_nothing(self): + self.assertIsNone(V.typed_value_kind("0/1", "String")) + + def test_a_comma_list_gains_nothing(self): + """A list is represented by value items, not one typed literal.""" + self.assertIsNone(V.typed_value_kind("1,2", "Integer")) + + def test_the_missing_token_gains_nothing(self): + self.assertIsNone(V.typed_value_kind(".", "Integer")) + + def test_a_BCF_reserved_integer_gains_nothing(self): + """VCF 4.5 reserves -2147483648..-2147483641; they keep only fieldValue.""" + for value in ("-2147483648", "-2147483641"): + with self.subTest(value=value): + self.assertIsNone(V.typed_value_kind(value, "Integer")) + self.assertEqual(V.typed_value_kind("-2147483640", "Integer"), "Integer") + + def test_a_non_finite_float_gains_nothing(self): + """INF and NaN have no xsd:decimal form.""" + for value in ("NaN", "inf", "-INF", "INFINITY"): + with self.subTest(value=value): + self.assertIsNone(V.typed_value_kind(value, "Float")) + + def test_a_value_that_does_not_parse_gains_nothing(self): + self.assertIsNone(V.typed_value_kind("twelve", "Integer")) + + +class DeclaredFieldTypeTests(VerboseTestCase): + """FORMAT declarations are read the same way INFO ones are.""" + + def test_format_declarations_are_read_when_asked_for(self): + header = "\n".join([ + '##INFO=', + '##FORMAT=', + '##FORMAT=', + ]) + self.assertEqual(V.declared_field_types(header, "FORMAT"), + {"AD": "Integer", "GT": "String"}) + + def test_asking_for_one_key_does_not_return_the_other(self): + header = "\n".join([ + '##INFO=', + '##FORMAT=', + ]) + self.assertEqual(V.declared_field_types(header, "INFO"), {"DP": "Integer"}) + self.assertEqual(V.declared_field_types(header, "FORMAT"), {"DP": "Float"}) + + def test_filter_declarations_expose_their_ids(self): + """A failure code links to its declaration, so the ids must be readable.""" + header = '##FILTER=' + self.assertEqual(set(V.declared_field_types(header, "FILTER")), {"q10"}) + + class InfoDeclaredTypeTests(VerboseTestCase): def test_each_declared_info_key_maps_to_its_type(self): """The INFO Type drives typed-value counting, so it is read per key.""" @@ -650,3 +756,35 @@ def test_no_progress_path_stays_none(self): if __name__ == "__main__": unittest.main() + + +class RecordDecompositionCountTests(VerboseTestCase): + """The census counts the same ID and FILTER parts the emitter builds.""" + + def _cols(self, record_id=".", filter_column="."): + return ["chr1", "1", record_id, "A", "G", ".", filter_column, ".", "GT", "0/1"] + + def test_semicolon_separated_identifiers_are_counted(self): + self.assertEqual(V.record_decomposition_counts(self._cols("idA;idB"), set()), + (2, 0, 0)) + + def test_the_missing_id_token_contributes_no_identifier(self): + self.assertEqual(V.record_decomposition_counts(self._cols("."), set()), (0, 0, 0)) + + def test_pass_and_missing_are_statuses_not_codes(self): + for value in ("PASS", "."): + with self.subTest(filter=value): + self.assertEqual( + V.record_decomposition_counts(self._cols(".", value), {"q10"}), + (0, 0, 0), + ) + + def test_failure_codes_are_counted_and_matched_against_declarations(self): + """Only a code with a ##FILTER line can cite one.""" + self.assertEqual( + V.record_decomposition_counts(self._cols(".", "q10;s50"), {"q10"}), + (0, 2, 1), + ) + + def test_a_truncated_record_is_treated_as_missing_rather_than_crashing(self): + self.assertEqual(V.record_decomposition_counts(["chr1", "1"], set()), (0, 0, 0)) diff --git a/test/test_vcf45_accessors_unit.py b/test/test_vcf45_accessors_unit.py new file mode 100644 index 0000000..b879832 --- /dev/null +++ b/test/test_vcf45_accessors_unit.py @@ -0,0 +1,363 @@ +"""Unit tests for the VCF 4.5 accessors added alongside the local-allele fix. + +Each test states the specification rule it pins, because several of these are +places where a plausible-looking wrong answer is the failure mode: a depth bound +to the wrong allele, or a modification value bound to the wrong base, comes back +well-formed and with full provenance attached. +""" +import unittest + +import vcf_rdfizer as R +import vcf_rdfizer_vocab as V +from test.helpers import VerboseTestCase + + +class LocalAlleleResolutionTests(VerboseTestCase): + """Number=LA/LR index the sample's own allele subset, not the ALT column.""" + + def test_LR_item_zero_is_the_reference(self): + """A Number=LR list runs REF first, then the local ALT alleles.""" + link = V.value_item_link("LR", 0, tuple_arity=None, local_alleles=[2, 4]) + self.assertEqual(link.allele_index, 0) + + def test_LR_items_resolve_through_LAA_not_by_position(self): + """With LAA=2,4 the list is REF, ALT2, ALT4 -- not REF, ALT1, ALT2.""" + got = [ + V.value_item_link("LR", i, tuple_arity=None, local_alleles=[2, 4]).allele_index + for i in range(3) + ] + self.assertEqual(got, [0, 2, 4]) + + def test_LA_items_skip_the_reference(self): + """Number=LA is one value per local ALT allele, so item 0 is LAA[0].""" + got = [ + V.value_item_link("LA", i, tuple_arity=None, local_alleles=[2, 4]).allele_index + for i in range(2) + ] + self.assertEqual(got, [2, 4]) + + def test_without_LAA_a_local_item_gets_no_allele(self): + """A wrong allele link is worse than none, so absence wins.""" + link = V.value_item_link("LR", 1, tuple_arity=None, local_alleles=None) + self.assertIsNone(link.allele_index) + self.assertTrue(link.is_empty) + + def test_an_item_past_the_end_of_LAA_gets_no_allele(self): + """A payload longer than the local set resolves no further.""" + link = V.value_item_link("LR", 5, tuple_arity=None, local_alleles=[2]) + self.assertIsNone(link.allele_index) + + def test_global_codes_are_untouched_by_a_local_allele_set(self): + """Number=R still indexes the ALT column even when LAA is present.""" + link = V.value_item_link("R", 1, tuple_arity=None, local_alleles=[2, 4]) + self.assertEqual(link.allele_index, 1) + + +class BaseModificationKeyTests(VerboseTestCase): + """VCF 4.5 reserves both a ChEBI-numeric spelling and a named alias.""" + + def test_the_numeric_form_parses(self): + parsed = V.parse_base_modification_key("M27551C") + self.assertEqual((parsed.family, parsed.chebi_id, parsed.base), ("M", "27551", "C")) + self.assertFalse(parsed.is_alias) + + def test_an_alias_resolves_to_the_same_residue(self): + """M5mC is documented as an alias for M27551C.""" + alias = V.parse_base_modification_key("M5mC") + numeric = V.parse_base_modification_key("M27551C") + self.assertEqual(alias.residue_uri, numeric.residue_uri) + self.assertTrue(alias.is_alias) + + def test_an_alias_keeps_its_written_spelling(self): + """The IRI must trace back to the column as the file wrote it.""" + parsed = V.parse_base_modification_key("M5mC") + self.assertEqual(parsed.key, "M5mC") + self.assertEqual(parsed.canonical_key, "M27551C") + + def test_all_three_families_share_one_modification_identity(self): + """M, DPM and ADM describe one modification, so they must name one.""" + ids = { + V.parse_base_modification_key(k).modification_id + for k in ("M5mC", "DPM5mC", "ADM5mC") + } + self.assertEqual(ids, {"27551C"}) + + def test_each_family_carries_its_own_value_property(self): + self.assertEqual( + [V.parse_base_modification_key(k).value_property + for k in ("M5mC", "DPM5mC", "ADM5mC")], + ["modificationFraction", "modificationDepth", "modificationAlleleDepth"], + ) + + def test_a_key_that_is_neither_spelling_is_rejected(self): + for key in ("NOTAKEY", "M5mCX", "GT", "M", ""): + with self.subTest(key=key): + self.assertIsNone(V.parse_base_modification_key(key)) + + +class BaseModificationPositionTests(VerboseTestCase): + """One value per base, either strand, that could carry the modification.""" + + def test_the_specification_worked_example(self): + """VCF 4.5: an allele of CGA has two M5mC values. + + The forward-strand C at the first base, and the reverse-strand C at the + second, whose forward base is G. + """ + self.assertEqual( + V.base_modification_positions([(0, "CGA")], "C"), + [(0, 0, "+"), (0, 1, "-")], + ) + + def test_bases_that_cannot_carry_it_yield_nothing(self): + self.assertEqual(V.base_modification_positions([(0, "AAA")], "C"), []) + + def test_N_yields_both_strands_negative_after_positive(self): + """The specification calls N out: every base could carry it.""" + self.assertEqual( + V.base_modification_positions([(0, "AC")], "N"), + [(0, 0, "+"), (0, 0, "-"), (0, 1, "+"), (0, 1, "-")], + ) + + def test_positions_run_across_alleles_in_genotype_order(self): + """Values are encoded over the concatenated genotype allele bases.""" + self.assertEqual( + V.base_modification_positions([(0, "C"), (2, "GC")], "C"), + [(0, 0, "+"), (2, 0, "-"), (2, 1, "+")], + ) + + def test_an_empty_allele_list_yields_nothing(self): + self.assertEqual(V.base_modification_positions([], "C"), []) + + +class GenotypeAlleleSequenceTests(VerboseTestCase): + """Number=M is encoded over the called alleles, in GT order.""" + + def test_alleles_come_back_in_genotype_order(self): + got = R._gt_allele_sequences(("2/0",), ("GT",), "G", "A,C,T") + self.assertEqual(got, [(2, "C"), (0, "G")]) + + def test_symbolic_alleles_contribute_nothing(self): + """The specification: symbolic alleles encode no modification values.""" + got = R._gt_allele_sequences(("0/3",), ("GT",), "G", "A,C,<*>") + self.assertEqual(got, [(0, "G")]) + + def test_no_call_positions_contribute_nothing(self): + got = R._gt_allele_sequences(("./1",), ("GT",), "G", "A") + self.assertEqual(got, [(1, "A")]) + + def test_a_record_without_GT_yields_nothing(self): + self.assertEqual(R._gt_allele_sequences(("12",), ("DP",), "G", "A"), []) + + +class PhaseIndicatorTests(VerboseTestCase): + """Each allele call carries the / or | that precedes it.""" + + def test_an_unphased_genotype_gives_every_call_a_slash(self): + parsed = V.parse_genotype("0/1") + self.assertEqual([c.phase_indicator for c in parsed.calls], ["/", "/"]) + + def test_a_phased_genotype_gives_every_call_a_bar(self): + parsed = V.parse_genotype("0|1") + self.assertEqual([c.phase_indicator for c in parsed.calls], ["|", "|"]) + + def test_mixed_indicators_are_kept_per_call(self): + """A single Phased/Unphased verdict would misdescribe this genotype. + + The first call's indicator is omitted in the source; it takes "/" here + because the genotype is not wholly phased. + """ + parsed = V.parse_genotype("0|1/2") + self.assertEqual([c.phase_indicator for c in parsed.calls], ["/", "|", "/"]) + + def test_a_bare_haploid_call_is_not_called_phased(self): + """There is no separator, so there is nothing to phase against. + + Defaulting to "|" here would contradict the genotype's own + phasingStatus, which is Unphased for exactly that reason. + """ + parsed = V.parse_genotype("0") + self.assertEqual([c.phase_indicator for c in parsed.calls], ["/"]) + self.assertFalse(parsed.is_phased) + + def test_an_explicit_leading_indicator_is_kept(self): + """VCF 4.4 added the leading indicator; it is not invented or dropped.""" + parsed = V.parse_genotype("|0/1") + self.assertEqual(parsed.calls[0].phase_indicator, "|") + + +if __name__ == "__main__": + unittest.main() + + +def _emit_to_list(): + """Collect emitted N-Triples lines, with the vocab prefix folded away.""" + out = [] + def emit(line): + out.append(line.strip().replace("https://w3id.org/vcf-core/vocab#", "vcfc:")) + return emit, out + + +def _objects(lines, predicate): + return [l.split("> ", 2)[-1].rstrip(" .").strip('"<>') + for l in lines if f"" in l] + + +class RecordIdentifierEmissionTests(VerboseTestCase): + """The ID column is a semicolon-separated list, kept in source order.""" + + def test_each_component_becomes_an_ordered_identifier(self): + emit, out = _emit_to_list() + stats = {"record_identifiers": 0} + R._emit_record_identifiers(emit, record_uri="F#record/1", value="idA;idB", stats=stats) + self.assertEqual(stats["record_identifiers"], 2) + self.assertEqual(_objects(out, "identifierValue"), ["idA", "idB"]) + self.assertEqual([o.split("^^")[0].strip('"') for o in _objects(out, "componentIndex")], + ["1", "2"]) + + def test_the_missing_token_is_a_status_not_an_identifier(self): + emit, out = _emit_to_list() + stats = {"record_identifiers": 0} + R._emit_record_identifiers(emit, record_uri="F#record/1", value=".", stats=stats) + self.assertEqual(out, []) + self.assertEqual(stats["record_identifiers"], 0) + + +class FilterCodeEmissionTests(VerboseTestCase): + """PASS and the missing token are FILTER statuses, not failure codes.""" + + def _run(self, value, definitions=None): + emit, out = _emit_to_list() + stats = {"filter_codes": 0} + R._emit_filter_codes(emit, subject_uri="F#call/1", value=value, + filter_definitions=definitions or {}, stats=stats) + return out, stats + + def test_pass_is_a_status_with_no_codes(self): + out, stats = self._run("PASS") + self.assertIn("vcfc:FiltersPassed", " ".join(out)) + self.assertEqual(stats["filter_codes"], 0) + + def test_the_missing_token_means_filters_were_not_applied(self): + out, stats = self._run(".") + self.assertIn("vcfc:FiltersNotApplied", " ".join(out)) + self.assertEqual(stats["filter_codes"], 0) + + def test_failure_codes_are_ordered_and_typed(self): + out, stats = self._run("q10;s50") + self.assertIn("vcfc:FiltersFailed", " ".join(out)) + self.assertEqual(stats["filter_codes"], 2) + self.assertEqual(_objects(out, "filterCodeValue"), ["q10", "s50"]) + + def test_a_declared_code_cites_its_header_line(self): + """So a consumer can read the description without re-reading the header.""" + out, _ = self._run("q10", {"q10": "F#header/line/5"}) + self.assertEqual(_objects(out, "declaredByFilter"), ["F#header/line/5"]) + + def test_an_undeclared_code_is_still_emitted_without_a_link(self): + out, stats = self._run("q10", {}) + self.assertEqual(stats["filter_codes"], 1) + self.assertEqual(_objects(out, "declaredByFilter"), []) + + +class BaseModificationEmissionTests(VerboseTestCase): + """One carrier per modification; the families contribute their values to it.""" + + def test_the_three_families_share_one_carrier(self): + emit, out = _emit_to_list() + stats = {"base_modifications": 0} + seen = set() + uris = set() + for key, value in (("M5mC", "0.9"), ("DPM5mC", "20"), ("ADM5mC", "18")): + uris.add(R._emit_base_modification( + emit, sample_uri="F#sample/1/S1", + modification=V.parse_base_modification_key(key), + value=value, stats=stats, emitted_modifications=seen)) + self.assertEqual(len(uris), 1, "M, DPM and ADM must name one modification") + self.assertEqual(stats["base_modifications"], 1) + joined = " ".join(out) + for prop in ("modificationFraction", "modificationDepth", "modificationAlleleDepth"): + self.assertIn(f"vcfc:{prop}", joined) + + def test_a_missing_value_emits_no_carrier(self): + emit, out = _emit_to_list() + stats = {"base_modifications": 0} + uri = R._emit_base_modification( + emit, sample_uri="F#sample/1/S1", + modification=V.parse_base_modification_key("M5mC"), + value=".", stats=stats, emitted_modifications=set()) + self.assertIsNone(uri) + self.assertEqual(out, []) + + +class ModificationItemTests(VerboseTestCase): + """A Number=M payload decomposes over the bases that can carry it.""" + + def _run(self, value, sequences, base="C"): + emit, out = _emit_to_list() + stats = {"value_items": 0, "modification_arity_mismatches": 0, + "modification_reverse_strand": 0} + R._emit_modification_items( + emit, parent_uri="F#fmt/M5mC", modification_uri="F#basemod/27551C", + modification=V.parse_base_modification_key("M5mC" if base == "C" else "MXaoN"), + value=value, record_uri="F#record/1", allele_sequences=sequences, stats=stats) + return out, stats + + def test_each_item_names_its_allele_offset_and_modification(self): + out, stats = self._run("0.9", [(0, "C")]) + self.assertEqual(stats["value_items"], 1) + self.assertEqual(_objects(out, "forAllele"), ["F#record/1/allele/0"]) + self.assertEqual(_objects(out, "forBaseModification"), ["F#basemod/27551C"]) + self.assertEqual([o.split("^^")[0].strip('"') for o in _objects(out, "modifiedBaseOffset")], + ["0"]) + + def test_a_reverse_strand_base_is_counted(self): + """CGA: the G at offset 1 carries the modification on the reverse strand.""" + out, stats = self._run("0.1,0.2", [(0, "CGA")]) + self.assertEqual(stats["value_items"], 2) + self.assertEqual(stats["modification_reverse_strand"], 1) + + def test_a_payload_that_disagrees_with_the_sequence_emits_nothing(self): + """A value on the wrong base is worse than a whole list on fieldValue.""" + out, stats = self._run("0.1,0.2,0.3", [(0, "C")]) + self.assertEqual(out, []) + self.assertEqual(stats["value_items"], 0) + self.assertEqual(stats["modification_arity_mismatches"], 1) + + def test_a_missing_payload_emits_nothing(self): + out, stats = self._run(".", [(0, "C")]) + self.assertEqual(out, []) + self.assertEqual(stats["modification_arity_mismatches"], 0) + + +class PhasingStatusEmissionTests(VerboseTestCase): + """A genotype whose indicators disagree is neither Phased nor Unphased.""" + + def _status(self, gt): + emit, out = _emit_to_list() + stats = {"genotypes": 0, "genotype_calls": 0} + R._emit_genotype(emit, sample_uri="F#sample/1/S1", record_uri="F#record/1", + value=gt, alt_count=3, stats=stats) + return [o for o in _objects(out, "phasingStatus")] + + def test_mixed_indicators_emit_MixedPhasing(self): + """vcfc:MixedPhasing exists for exactly this; the per-call indicators + carry the precise semantics.""" + self.assertEqual(self._status("0|1/2"), ["vcfc:MixedPhasing"]) + + def test_a_wholly_phased_genotype_is_Phased(self): + self.assertEqual(self._status("0|1"), ["vcfc:Phased"]) + + def test_a_wholly_unphased_genotype_is_Unphased(self): + self.assertEqual(self._status("0/1"), ["vcfc:Unphased"]) + + def test_a_malformed_genotype_produces_no_resource(self): + """The value survives on fieldValue; inventing a Genotype would fail + vcfc:GenotypeShape.""" + self.assertEqual(self._status("not-a-genotype"), []) + + +class MalformedGenotypeSequenceTests(VerboseTestCase): + def test_a_genotype_outside_the_lexical_space_yields_no_sequences(self): + """Number=M cannot be decomposed without a parsable GT.""" + self.assertEqual(R._gt_allele_sequences(("x/y",), ("GT",), "G", "A"), []) diff --git a/test/validation_fixtures.py b/test/validation_fixtures.py index 40b3430..7468f84 100644 --- a/test/validation_fixtures.py +++ b/test/validation_fixtures.py @@ -496,6 +496,32 @@ def _format_typed_shape() -> dict: } +def _record_decomposition_shape() -> dict: + """ID and FILTER components the emitter turns into ordered resources. + + Derived from the fixture's own columns, mirroring the rule rather than + calling the emitter: PASS and the missing token are FILTER statuses, so + neither contributes a failure code. + """ + declared_filters = set(_declared_ids("FILTER")) + identifiers = codes = declared_codes = 0 + for row in _record_rows(): + record_id = row[2] if len(row) > 2 else "." + if record_id != ".": + identifiers += sum(1 for part in record_id.split(";") if part) + filter_column = row[6] if len(row) > 6 else "." + if filter_column not in ("PASS", "."): + for code in (part for part in filter_column.split(";") if part): + codes += 1 + if code in declared_filters: + declared_codes += 1 + return { + "recordIdentifierCount": identifiers, + "filterCodeCount": codes, + "declaredFilterCodeCount": declared_codes, + } + + def _info_shape() -> dict: """INFO counters the census needs, mirroring the emitter's typing rules.""" declared = {} @@ -598,6 +624,7 @@ def parser_summary( "headerLineCount": len(HEADER_LINES), **_format_shape(), **_format_typed_shape(), + **_record_decomposition_shape(), **_info_shape(), "fileFormat": FILE_FORMAT, "referenceGenome": REFERENCE_GENOME, diff --git a/vcf_rdfizer.py b/vcf_rdfizer.py index d36502a..144cfe5 100644 --- a/vcf_rdfizer.py +++ b/vcf_rdfizer.py @@ -2432,9 +2432,17 @@ def _emit_genotype( f"<{genotype_uri}> <{_vocab('ploidy')}> " f"{_ordinal_literal(parsed.ploidy)} .\n" ) + indicators = {call.phase_indicator for call in parsed.calls} + if len(indicators) > 1: + # vcfc:MixedPhasing exists for exactly this: the per-call indicators + # carry the precise semantics, and a single Phased/Unphased verdict + # would misdescribe the genotype. + status = "MixedPhasing" + else: + status = "Phased" if parsed.is_phased else "Unphased" emit( f"<{genotype_uri}> <{_vocab('phasingStatus')}> " - f"<{_vocab('Phased' if parsed.is_phased else 'Unphased')}> .\n" + f"<{_vocab(status)}> .\n" ) for call in parsed.calls: call_uri = f"{genotype_uri}/call/{call.call_index}" @@ -2451,6 +2459,10 @@ def _emit_genotype( f"<{call_uri}> <{_vocab('isNoCall')}> " f'"{"true" if call.is_no_call else "false"}"^^<{XSD_BOOLEAN_URI}> .\n' ) + emit( + f"<{call_uri}> <{_vocab('phaseIndicator')}> " + f"{_ntriples_string_literal(call.phase_indicator)} .\n" + ) if call.allele_index is not None and call.allele_index <= alt_count: emit( f"<{call_uri}> <{_vocab('calledAllele')}> " @@ -2509,39 +2521,144 @@ def _emit_local_allele_set( for local_index, allele_index in enumerate(indices, start=1): allele_uri = f"{record_uri}/allele/{allele_index}" emit(f"<{set_uri}> <{_vocab('hasLocalAllele')}> <{allele_uri}> .\n") + # The local ordinal belongs to the membership, not to the allele: the + # allele resource is shared across samples, and two samples can give the + # same ALT different local positions. vcfc:localAlleleIndex is deprecated + # for exactly that reason. + membership_uri = f"{set_uri}/member/{local_index}" + emit( + f"<{set_uri}> <{_vocab('hasLocalAlleleMembership')}> " + f"<{membership_uri}> .\n" + ) + emit( + f"<{membership_uri}> <{RDF_TYPE_URI}> " + f"<{_vocab('LocalAlleleMembership')}> .\n" + ) + emit(f"<{membership_uri}> <{_vocab('localAllele')}> <{allele_uri}> .\n") emit( - f"<{allele_uri}> <{_vocab('localAlleleIndex')}> " + f"<{membership_uri}> <{_vocab('localIndex')}> " f"{_ordinal_literal(local_index)} .\n" ) stats["local_allele_sets"] += 1 -def _emit_base_modification( - emit, *, sample_uri: str, modification, value: str, stats: dict +def _gt_allele_sequences(sample_values, format_keys, ref: str, alt: str): + """This sample's called allele sequences, in GT order. + + VCF 4.5 encodes a Number=M field over the concatenated genotype allele + bases, so the order is the GT's, not the ALT column's. Missing and symbolic + alleles are dropped: the specification says they "encode no base + modification values". + """ + if "GT" not in format_keys: + return [] + parsed = vocab.parse_genotype(sample_values[format_keys.index("GT")]) + if parsed is None: + return [] + alleles = [ref] + ([] if vocab.is_missing(alt) else alt.split(",")) + out = [] + for call in parsed.calls: + index = call.allele_index + if index is None or index >= len(alleles): + continue + sequence = alleles[index] + if vocab.is_missing(sequence) or sequence.startswith("<"): + continue + out.append((index, sequence)) + return out + + +def _emit_modification_items( + emit, + *, + parent_uri: str, + modification_uri: str, + modification, + value: str, + record_uri: str, + allele_sequences, + stats: dict, ) -> None: - """Emit the vcfc:BaseModification carrier for one M/DPM/ADM FORMAT key.""" - if vocab.is_missing(value): + """Decompose a Number=M payload into one item per modifiable base. + + Each item names the allele the base belongs to, its offset within that + allele, and the modification it reports, so a query can reach one base's + value without splitting the comma list. When the payload and the positions + the sequences imply disagree, no item is emitted: a value bound to the wrong + base would be worse than leaving the list whole on vcfc:fieldValue. + """ + items = vocab.split_value_items(value) + if not items: return - modification_uri = f"{sample_uri}/basemod/{_rml_uri_component(modification.key)}" - emit( - f"<{modification_uri}> <{RDF_TYPE_URI}> " - f"<{_vocab('BaseModification')}> .\n" + positions = vocab.base_modification_positions( + allele_sequences, modification.base ) - emit( - f"<{modification_uri}> <{_vocab('modifiedResidue')}> " - f"<{modification.residue_uri}> .\n" - ) - # VCF 4.5 gives the base-modification FORMAT families a per-base offset; - # without an explicit one, position 0 is the modified base named by the key. - emit( - f"<{modification_uri}> <{_vocab('modifiedBaseOffset')}> " - f"{_ordinal_literal(0)} .\n" + if len(positions) != len(items): + stats["modification_arity_mismatches"] += 1 + return + for item_index, (item, (allele_index, offset, strand)) in enumerate( + zip(items, positions, strict=True) + ): + item_uri = f"{parent_uri}/value/{item_index}" + emit(f"<{parent_uri}> <{_vocab('hasValueItem')}> <{item_uri}> .\n") + emit(f"<{item_uri}> <{RDF_TYPE_URI}> <{_vocab('FieldValueItem')}> .\n") + emit( + f"<{item_uri}> <{_vocab('valueIndex')}> " + f"{_ordinal_literal(item_index)} .\n" + ) + emit(f"<{item_uri}> <{_vocab('itemValue')}> {_ntriples_literal(item)} .\n") + emit( + f"<{item_uri}> <{_vocab('forAllele')}> " + f"<{record_uri}/allele/{allele_index}> .\n" + ) + emit( + f"<{item_uri}> <{_vocab('forBaseModification')}> " + f"<{modification_uri}> .\n" + ) + emit( + f"<{item_uri}> <{_vocab('modifiedBaseOffset')}> " + f"{_ordinal_literal(offset)} .\n" + ) + stats["value_items"] += 1 + if strand == "-": + stats["modification_reverse_strand"] += 1 + + +def _emit_base_modification( + emit, + *, + sample_uri: str, + modification, + value: str, + stats: dict, + emitted_modifications: set, +) -> str | None: + """Emit the vcfc:BaseModification carrier for one M/DPM/ADM FORMAT key.""" + if vocab.is_missing(value): + return None + # Keyed by the modification itself, not by the FORMAT key: M, DPM and ADM + # describe one modification from three angles, and vcfc:BaseModification is + # the modification. Keying by column would split it across three subjects + # and no query could ask for a fraction and its depth together. + modification_uri = ( + f"{sample_uri}/basemod/{_rml_uri_component(modification.modification_id)}" ) + if modification_uri not in emitted_modifications: + emitted_modifications.add(modification_uri) + emit( + f"<{modification_uri}> <{RDF_TYPE_URI}> " + f"<{_vocab('BaseModification')}> .\n" + ) + emit( + f"<{modification_uri}> <{_vocab('modifiedResidue')}> " + f"<{modification.residue_uri}> .\n" + ) + stats["base_modifications"] += 1 emit( f"<{modification_uri}> <{_vocab(modification.value_property)}> " f"{_ntriples_literal(value)} .\n" ) - stats["base_modifications"] += 1 + return modification_uri def append_expanded_sample_rdf( @@ -2577,6 +2694,8 @@ def append_expanded_sample_rdf( "phase_sets": 0, "local_allele_sets": 0, "base_modifications": 0, + "modification_arity_mismatches": 0, + "modification_reverse_strand": 0, "value_items": 0, "triples": 0, "appended_bytes": 0, @@ -2662,6 +2781,18 @@ def produce(emit): sample_call_count += 1 phase_fields: dict[str, str] = {} + emitted_modifications: set[str] = set() + # LAA has to be read before any Number=LA/LR field, and the + # FORMAT string is free to list it afterwards. + local_alleles: list[int] = [] + if version.local_alleles and "LAA" in format_keys: + local_alleles = [ + index + for index in vocab.parse_local_allele_indices( + sample_values[format_keys.index("LAA")] + ) + if index <= alt_count + ] for format_index, format_key in enumerate(format_keys): format_value = sample_values[format_index] format_component = format_components.get(format_key) @@ -2720,6 +2851,7 @@ def produce(emit): }, version=version, stats=stats, + local_alleles=local_alleles, ) format_value_count += 1 @@ -2770,13 +2902,30 @@ def produce(emit): ) modification = base_modifications[format_key] if modification is not None: - _emit_base_modification( + modification_uri = _emit_base_modification( emit, sample_uri=sample_uri, modification=modification, value=format_value, stats=stats, + emitted_modifications=emitted_modifications, ) + if modification_uri is not None: + _emit_modification_items( + emit, + parent_uri=format_uri, + modification_uri=modification_uri, + modification=modification, + value=format_value, + record_uri=record_uri, + allele_sequences=_gt_allele_sequences( + sample_values, + format_keys, + record.ref, + record.alt, + ), + stats=stats, + ) if phase_fields: _emit_phase_set( @@ -3374,6 +3523,81 @@ def _emit_breakend(emit, allele_uri: str, breakend) -> None: ) +def _emit_record_identifiers(emit, *, record_uri: str, value: str, stats: dict) -> None: + """Decompose the semicolon-separated ID column into ordered components. + + The whole column stays on vcfc:recordId; this adds the parts, so a query can + ask for one identifier without splitting a string. The missing token is a + status, not an identifier, so '.' produces no component. + """ + if vocab.is_missing(value): + return + for position, component in enumerate( + (part for part in value.split(";") if part), start=1 + ): + identifier_uri = f"{record_uri}/id/{position}" + emit(f"<{record_uri}> <{_vocab('hasIdentifier')}> <{identifier_uri}> .\n") + emit( + f"<{identifier_uri}> <{RDF_TYPE_URI}> " + f"<{_vocab('RecordIdentifier')}> .\n" + ) + emit( + f"<{identifier_uri}> <{_vocab('identifierValue')}> " + f"{_ntriples_string_literal(component)} .\n" + ) + emit( + f"<{identifier_uri}> <{_vocab('componentIndex')}> " + f"{_ordinal_literal(position)} .\n" + ) + stats["record_identifiers"] += 1 + + +#: PASS and the missing token are FILTER statuses rather than failure codes. +_FILTER_STATUS = {"PASS": "FiltersPassed", ".": "FiltersNotApplied"} + + +def _emit_filter_codes( + emit, *, subject_uri: str, value: str, filter_definitions: dict, stats: dict +) -> None: + """Emit the FILTER/FT status and, when it failed, the ordered codes. + + A failure code links back to its '##FILTER' declaration where one exists, so + a consumer can read the description without re-reading the header. + """ + status = _FILTER_STATUS.get(value.strip()) + if status is not None: + emit( + f"<{subject_uri}> <{_vocab('filterStatus')}> " + f"<{_vocab(status)}> .\n" + ) + return + emit( + f"<{subject_uri}> <{_vocab('filterStatus')}> " + f"<{_vocab('FiltersFailed')}> .\n" + ) + for position, code in enumerate( + (part for part in value.split(";") if part), start=1 + ): + code_uri = f"{subject_uri}/filter/{position}" + emit(f"<{subject_uri}> <{_vocab('hasFilterCode')}> <{code_uri}> .\n") + emit(f"<{code_uri}> <{RDF_TYPE_URI}> <{_vocab('FilterCode')}> .\n") + emit( + f"<{code_uri}> <{_vocab('filterCodeValue')}> " + f"{_ntriples_string_literal(code)} .\n" + ) + emit( + f"<{code_uri}> <{_vocab('componentIndex')}> " + f"{_ordinal_literal(position)} .\n" + ) + declaration = filter_definitions.get(code) + if declaration is not None: + emit( + f"<{code_uri}> <{_vocab('declaredByFilter')}> " + f"<{declaration}> .\n" + ) + stats["filter_codes"] += 1 + + def _emit_value_items( emit, *, @@ -3384,6 +3608,7 @@ def _emit_value_items( allele_uris: dict[int, str], version: "vocab.VCFVersion", stats: dict, + local_alleles: list[int] | None = None, ) -> None: """Decompose a comma-separated payload into indexed vcfc:FieldValueItems. @@ -3421,7 +3646,9 @@ def _emit_value_items( f"<{item_uri}> <{_vocab('tupleArity')}> " f"{_ordinal_literal(tuple_arity)} .\n" ) - link = version.value_item_link(field_key, number, item_index) + link = version.value_item_link( + field_key, number, item_index, local_alleles=local_alleles + ) if link.allele_index is not None: allele_uri = allele_uris.get(link.allele_index) if allele_uri is not None: @@ -3751,12 +3978,17 @@ def append_record_detail_rdf( ) -> dict: """Append per-record detail the RML mapping cannot express. - Three groups of things land here, all of them per-record, so they share one + Four groups of things land here, all of them per-record, so they share one pass over ``records.tsv``: * The fixed fields whose datatype depends on the row -- ID, ALT, QUAL and FILTER. Each may be the VCF missing token, which the vocabulary requires as ``"."^^vcfc:Null``, and RML cannot switch datatype per row. + * The decompositions of those cells: ``vcfc:RecordIdentifier`` per + semicolon-separated ID, the ``vcfc:filterStatus`` and its + ``vcfc:FilterCode`` resources, and the ordered ``vcfc:FormatKey`` list. + A single source cell becomes several resources, which is the other thing + RML cannot do here. * The allele layer: ``vcfc:ReferenceAllele`` and ``vcfc:AltAllele`` resources parsed out of REF and ALT, with their index, kind, symbolic type and breakend components, plus the ``vcfc:chromosome`` link to the @@ -3771,6 +4003,8 @@ def append_record_detail_rdf( "qual_values": 0, "info_values": 0, "info_definitions": 0, + "record_identifiers": 0, + "filter_codes": 0, "alleles": 0, "breakends": 0, "value_items": 0, @@ -3800,6 +4034,18 @@ def append_record_detail_rdf( header_lines_tsv, "ALT" ).items() } + filter_definition_uris = { + identifier: definition.uri + for identifier, definition in _load_field_definitions( + header_lines_tsv, "FILTER" + ).items() + } + format_definition_uris = { + identifier: definition.uri + for identifier, definition in _load_field_definitions( + header_lines_tsv, "FORMAT" + ).items() + } with SampleRecordStream(records_tsv) as record_stream: if not record_stream.source_file: @@ -3834,6 +4080,40 @@ def produce(emit): f"<{call_uri}> <{_vocab('filter')}> " f"{_ntriples_literal(record.filter_value or '.')} .\n" ) + for key_position, format_key in enumerate(record.format_keys, start=1): + # IRI shape follows vcfc:FormatKey's own iriTemplate. + format_key_uri = f"{call_uri}/formatKey/{key_position}" + emit( + f"<{call_uri}> <{_vocab('hasFormatKey')}> " + f"<{format_key_uri}> .\n" + ) + emit( + f"<{format_key_uri}> <{RDF_TYPE_URI}> " + f"<{_vocab('FormatKey')}> .\n" + ) + emit( + f"<{format_key_uri}> <{_vocab('fieldIndex')}> " + f"{_ordinal_literal(key_position)} .\n" + ) + declaration = format_definition_uris.get(format_key) + if declaration is not None: + emit( + f"<{format_key_uri}> <{_vocab('declaredBy')}> " + f"<{declaration}> .\n" + ) + _emit_record_identifiers( + emit, + record_uri=record_uri, + value=record.record_id or ".", + stats=stats, + ) + _emit_filter_codes( + emit, + subject_uri=call_uri, + value=record.filter_value or ".", + filter_definitions=filter_definition_uris, + stats=stats, + ) emit( f"<{call_uri}> <{_vocab('infoRaw')}> " f"{_ntriples_literal(record.info or '.')} .\n" @@ -3900,7 +4180,7 @@ def produce(emit): ) info_entries = parse_info_entries(record.info) if emit_info else [] - for key, value in info_entries: + for entry_position, (key, value) in enumerate(info_entries, start=1): key_component = key_components.get(key) if key_component is None: key_component = _rml_uri_component(key) @@ -3947,6 +4227,10 @@ def produce(emit): f"<{info_uri}> <{_vocab('declaredBy')}> " f"<{definition.uri}> .\n" ) + emit( + f"<{info_uri}> <{_vocab('fieldIndex')}> " + f"{_ordinal_literal(entry_position)} .\n" + ) if value is None: emit( f"<{info_uri}> <{_vocab('fieldValueBoolean')}> " diff --git a/vcf_rdfizer_vocab.py b/vcf_rdfizer_vocab.py index 20c25b3..9b41106 100644 --- a/vcf_rdfizer_vocab.py +++ b/vcf_rdfizer_vocab.py @@ -254,13 +254,38 @@ def is_empty(self) -> bool: ) -def value_item_link(number: str, item_index: int, *, tuple_arity: int | None) -> ValueItemLink: +def _local_allele_index(position: int, local_alleles: list[int] | None) -> ValueItemLink: + """Resolve one position in a sample's local allele list to a global allele. + + ``local_alleles`` is that sample's LAA value: the one-based ALT indices it + declares, which are also the record-global allele indices, since global + index 0 is REF and ALT *i* is index *i*. Without an LAA the local ordering + is undefined, so the item gets no vcfc:forAllele rather than a positional + guess -- a wrong allele link is worse than none. + """ + if not local_alleles or position < 0 or position >= len(local_alleles): + return ValueItemLink() + return ValueItemLink(allele_index=local_alleles[position]) + + +def value_item_link( + number: str, + item_index: int, + *, + tuple_arity: int | None, + local_alleles: list[int] | None = None, +) -> ValueItemLink: """Resolve which allele or ordinal one comma-list item belongs to. ``item_index`` is zero-based within the comma list. When the field flattens fixed-width tuples, the allele stride is the tuple width, so items 0 and 1 of a per-ALT CIPOS list both belong to ALT allele 1. + The VCF 4.5 local-allele codes are indexed by the sample's own allele + subset, not by the record's ALT list: with ``LAA=2,4`` a ``Number=LR`` list + runs REF, ALT2, ALT4. Resolving them positionally would attach each value to + the wrong allele, so they go through ``local_alleles``. + This is the version-neutral form, driven only by the declared Number token. Prefer :meth:`VCFVersion.value_item_link`, which also knows which keys are positional in a given version despite declaring ``Number=.``. @@ -268,11 +293,20 @@ def value_item_link(number: str, item_index: int, *, tuple_arity: int | None) -> token = (number or "").strip() stride = tuple_arity or 1 group = item_index // stride + if token == "LA": + return _local_allele_index(group, local_alleles) + if token == "LR": + # Item 0 is the reference; the rest index the local ALT subset. + if group == 0: + return ValueItemLink(allele_index=0) + return _local_allele_index(group - 1, local_alleles) if token in _ARITY_PER_ALT: return ValueItemLink(allele_index=group + 1) if token in _ARITY_PER_ALLELE: return ValueItemLink(allele_index=group) if token in _ARITY_PER_GENOTYPE: + # Number=LG orders genotypes over the local allele set, so this ordinal + # is local to the sample; it is not a record-global genotype index. return ValueItemLink(genotype_index=group) if token == "P": return ValueItemLink(gt_allele_index=group) @@ -350,7 +384,13 @@ def is_positional(self, field_id: str, number: str) -> bool: return True return not value_item_link(number, 0, tuple_arity=None).is_empty - def value_item_link(self, field_id: str, number: str, item_index: int) -> ValueItemLink: + def value_item_link( + self, + field_id: str, + number: str, + item_index: int, + local_alleles: list[int] | None = None, + ) -> ValueItemLink: """Resolve one comma-list item's allele or ordinal, for this version.""" arity = self.tuple_arity(field_id) if arity is not None: @@ -359,7 +399,9 @@ def value_item_link(self, field_id: str, number: str, item_index: int) -> ValueI # has an index but nothing to point vcfc:forAllele at. return ValueItemLink() return ValueItemLink(allele_index=(item_index // arity) + 1) - return value_item_link(number, item_index, tuple_arity=None) + return value_item_link( + number, item_index, tuple_arity=None, local_alleles=local_alleles + ) def _version( @@ -709,6 +751,9 @@ class ParsedGenotypeCall: call_index: int #: The record-global allele index, or None when the position is a no-call. allele_index: int | None + #: The / or | preceding this position. VCF omits it before the first allele + #: except in 4.4+, so the effective one is carried here either way. + phase_indicator: str = "/" @property def is_no_call(self) -> bool: @@ -754,13 +799,23 @@ def parse_genotype(value: str) -> ParsedGenotype | None: allele_tokens = tokens[0::2] separators = tokens[1::2] + # The indicator before the first allele is usually absent. VCF phases a GT + # as a whole unless the indicators disagree, so an omitted one takes the + # value the rest of the string uses. + # Phased only when the string actually says so. A bare haploid call has no + # separator at all, and calling it phased would contradict the is_phased + # rule below, which holds that there is nothing to phase against. + default_indicator = "|" if "|" in value and "/" not in value else "/" + first_indicator = ("|" if leading_phased else "/") if value[0] in "|/" else default_indicator + calls = [] for call_index, token in enumerate(allele_tokens): + indicator = first_indicator if call_index == 0 else separators[call_index - 1] if token == "." or token == "": - calls.append(ParsedGenotypeCall(call_index, None)) + calls.append(ParsedGenotypeCall(call_index, None, indicator)) else: try: - calls.append(ParsedGenotypeCall(call_index, int(token))) + calls.append(ParsedGenotypeCall(call_index, int(token), indicator)) except ValueError: return None @@ -842,6 +897,24 @@ def genotype_index(allele_indices: list[int]) -> int: "ADM": "modificationAlleleDepth", } +#: VCF 4.5 also reserves a named alias for each common modification, so M5mC +#: and M27551C denote the same thing. The spec's reserved-FORMAT table lists +#: them as "Alias for ..."; the same ten modifications exist in all three +#: families. The vocabulary carries the correspondence as vcfc:aliasOf on each +#: reserved key, so a graph that resolves the alias agrees with it. +_BASE_MODIFICATION_ALIASES = { + "5mC": ("27551", "C"), + "5hmC": ("76792", "C"), + "5fC": ("76794", "C"), + "5caC": ("76793", "C"), + "5hmU": ("16964", "T"), + "5fU": ("80961", "T"), + "5caU": ("17477", "T"), + "6mA": ("28871", "A"), + "8oxoG": ("44605", "G"), + "XaoN": ("18107", "N"), +} + @dataclass(frozen=True) class ParsedBaseModification: @@ -852,28 +925,106 @@ class ParsedBaseModification: base: str #: The vcfc property that carries this family's value. value_property: str + #: The key as the file wrote it. An alias keeps its own spelling here so + #: the graph traces back to the source column, while the residue below + #: still resolves to the ChEBI the alias denotes. + source_key: str | None = None @property def residue_uri(self) -> str: return f"{CHEBI_NAMESPACE}{self.chebi_id}" @property - def key(self) -> str: + def canonical_key(self) -> str: return f"{self.family}{self.chebi_id}{self.base}" + @property + def modification_id(self) -> str: + """The modification itself, independent of which family reported it. + + M, DPM and ADM describe one modification from three angles -- fraction, + depth, allele depth -- so they must name one vcfc:BaseModification. + """ + return f"{self.chebi_id}{self.base}" + + @property + def key(self) -> str: + return self.source_key or self.canonical_key + + @property + def is_alias(self) -> bool: + return self.key != self.canonical_key + def parse_base_modification_key(key: str) -> ParsedBaseModification | None: - """Recognize a base-modification FORMAT key and resolve its ChEBI residue.""" + """Recognize a base-modification FORMAT key and resolve its ChEBI residue. + + Accepts both spellings VCF 4.5 reserves: the ChEBI-numeric form M27551C and + the named alias M5mC. Resolving the alias is what lets one query answer over + either spelling. + """ match = BASE_MODIFICATION_KEY_RE.match(key) - if not match: - return None - family = match.group("family") - return ParsedBaseModification( - family=family, - chebi_id=match.group("chebi"), - base=match.group("base"), - value_property=_BASE_MODIFICATION_PROPERTIES[family], - ) + if match: + family = match.group("family") + return ParsedBaseModification( + family=family, + chebi_id=match.group("chebi"), + base=match.group("base"), + value_property=_BASE_MODIFICATION_PROPERTIES[family], + source_key=key, + ) + for family in ("ADM", "DPM", "M"): # longest prefix first + if not key.startswith(family): + continue + alias = _BASE_MODIFICATION_ALIASES.get(key[len(family):]) + if alias is None: + continue + chebi_id, base = alias + return ParsedBaseModification( + family=family, + chebi_id=chebi_id, + base=base, + value_property=_BASE_MODIFICATION_PROPERTIES[family], + source_key=key, + ) + return None + + +#: Watson-Crick complements, for reading a modification off the reverse strand. +_COMPLEMENT = {"A": "T", "C": "G", "G": "C", "T": "A", "U": "A", "N": "N"} + + +def base_modification_positions( + allele_sequences: list[tuple[int, str]], modified_base: str +) -> list[tuple[int, int, str]]: + """Positions a Number=M field carries a value for, in source order. + + VCF 4.5: a Number=M field has one value per base, on either strand, of the + concatenated genotype allele sequences that could carry the modification, + ordered as the bases occur. The specification's own example is an allele of + CGA with two M5mC values -- the forward-strand C at the first base and the + reverse-strand C at the second, whose forward base is G. + + ``allele_sequences`` is (global allele index, sequence) in GT order, with + missing and symbolic alleles already dropped: the specification says they + "encode no base modification values". Returns (allele index, offset within + that allele, strand). + + N is the exception the specification calls out: every base could carry it, + so each position yields both strands, negative immediately after positive. + """ + base = (modified_base or "").upper() + out: list[tuple[int, int, str]] = [] + for allele_index, sequence in allele_sequences: + for offset, residue in enumerate(sequence.upper()): + if base == "N": + out.append((allele_index, offset, "+")) + out.append((allele_index, offset, "-")) + elif residue == base: + out.append((allele_index, offset, "+")) + elif _COMPLEMENT.get(residue) == base: + out.append((allele_index, offset, "-")) + return out def split_value_items(value: str) -> list[str]: