Skip to content

conversion.md

Latest commit

 

History

History
491 lines (405 loc) · 26.7 KB

conversion.md

File metadata and controls

491 lines (405 loc) · 26.7 KB
Raw

Conversion: VCF to RDF

What the conversion actually does to a VCF, in order, and what it decides on your behalf. The companion documents are architecture.md (how the pieces fit), rml-mappings.md (how to change the mapping), and vcf-coverage.md (which VCF element ends up where).

1. Input acceptance

--input takes a single file or a directory.

  • A file must be named *.vcf or *.vcf.gz. Nothing else is accepted — not .bcf, not .vcf.bgz, not a bare .gz. Extension, not content, decides.
  • A directory is enumerated one level deep, sorted, at the moment the run starts. Files that appear later are not picked up. Subdirectories are ignored.
  • Each input is processed independently and end to end before the next one begins. A failure is isolated to its input.

The output basename is the filename with .vcf / .vcf.gz removed. Two inputs that reduce to the same basename would collide, and the pre-flight collision check rejects the run before Docker starts.

2. Stage one — VCF to TSV

src/vcf_as_tsv.sh reads the VCF (through gzip -dc when compressed) in one awk pass and writes three tab-separated tables.

Output One row per Columns
<sample>.file_metadata.tsv file (exactly one row) SOURCE_FILE, FILE_FORMAT, FILE_DATE, SOURCE_SOFTWARE, REFERENCE_GENOME, HEADER_COUNT, RECORD_COUNT
<sample>.header_lines.tsv ## meta-information line SOURCE_FILE, HEADER_INDEX, HEADER_KEY, HEADER_VALUE, RAW_LINE
<sample>.records.tsv VCF data line SOURCE_FILE, ROW_ID, CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO, FORMAT, (sample payload)

Two more paths exist for compatibility, sample_calls.tsv and sample_format_values.tsv. Under the shipped mapping they are written header-only — see §5.

What this stage does, precisely:

  • ## lines are split on the first =. Everything after it is HEADER_VALUE; RAW_LINE keeps the original text minus the leading ##. No structural parsing happens here.
  • Four keys are lifted into file_metadata.tsv by case-insensitive match: fileformat, filedate, source, reference.
  • HEADER_INDEX counts ## lines and the #CHROM line, 1-based.
  • ROW_ID is a 1-based counter over data lines, stable within one file. It is the join key for every record-level IRI.
  • Trailing carriage returns are stripped from every field.
  • All sample columns are joined into one trailing field separated by single spaces, and runs of whitespace are collapsed. The column's name is the whitespace-joined sample IDs from #CHROM, or the literal SAMPLES when the VCF declares none — which is why a mapping cannot reference it by a fixed name.

Limitations of this stage, stated plainly. There is no htslib here. The parser is regex-and-field-index awk, which means:

  • The #CHROM line is recognised only when it is tab-delimited. A space-delimited header line is matched by no rule, so sample column names are lost and the records header falls back to SAMPLES — silently.
  • Nothing validates VCF spec conformance. A malformed file produces a malformed graph rather than an error, and the first thing that notices is the validation suite.
  • A data line with fewer than 8 columns yields empty strings for the missing fields rather than an error.
  • Because sample fields are whitespace-normalized, a value containing a literal space would be corrupted. The VCF specification forbids spaces in data fields, so this is only reachable with an already-invalid file — but it is not detected.

3. Stage two — TSV to RDF with RMLStreamer

src/run_conversion.sh runs RMLStreamer 2.5.0 over the TSVs using the mapping at --rules (default rules/default_rules.ttl), then normalizes the Flink part files and merges them into one aggregate.

The mapping's five csvw:url values are literal container paths (/data/tsv/records.tsv and friends) that the wrapper rewrites per input to /data/tsv/<sample>.records.tsv. That rewrite is why the contract in rml-mappings.md insists those strings stay verbatim.

--rdf-storage-mode decides how the aggregate is assembled:

Mode Behaviour
plain Parts are merged into one uncompressed <sample>.nt
space-optimized Each part is streamed through gzip into one <sample>.nt.gz and the source part is deleted immediately, so a full uncompressed copy never exists

The output is N-Triples. No named graphs are produced, and no blank nodes: every class in the vocabulary declares an IRI template, and preflight_blank_nodes treats a blank node as a validation failure precisely because it means a term map produced no IRI.

--spark-partitions is a parallelism hint passed through to RMLStreamer. It does not change the output, only how many parts are produced before the merge.

4. Stage three — the wrapper's own emitters

Three classes of triple cannot come from RML, and are appended to the aggregate by the host process before compression. See architecture.md for why. The rule is: RML carries every field whose datatype is the same for every row; the wrapper carries everything else — anything that may be the missing token, and anything that has to be decomposed out of a single source cell.

emit_record_detail — the row-dependent fixed fields, alleles, INFO

Four fixed fields and the raw INFO string are always emitted here, because each may be the VCF missing token and the vocabulary requires that as "."^^vcfc:Null:

Field Emitted as
ID the lexical value, or "."^^vcfc:Null
ALT the lexical value, or "."^^vcfc:Null
FILTER the lexical value, or "."^^vcfc:Null
INFO (vcfc:infoRaw) the lexical value, or "."^^vcfc:Null
QUAL see below

QUAL has its own table because its shape is a disjunction over the whole VCF Float lexical space:

QUAL value Emitted as
a finite number "<lexical form>"^^xsd:decimal — the source lexical form, so no precision is gained or lost
INF / INFINITY / NAN, any case "<lexical form>"^^vcfc:VCFFloat — these are in the VCF Float lexical space but in no XSD numeric one
. 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 vcfc:ReferenceAllele / vcfc:AltAllele with alleleIndex (0 for REF, then 1..n in source order), alleleValue and alleleKind. The kind is the VCF syntactic category: a base sequence, the overlapping-deletion *, the missing ., an angle-bracketed symbolic ID, the unspecified <*>/<NON_REF>, or a breakend. A symbolic allele is joined to its ##ALT declaration with declaredByAlt and, when its ID is one of the reserved codes, given a vcfc:svType. A breakend is parsed into its orientation, replacement string and mate position. The record is also linked to the contig declaration its 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 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 — A, R, LA, LR (one per allele), G, LG (one per genotype), P (one per GT allele) — is decomposed into ordered vcfc:FieldValueItem resources, each joined to the allele it describes with vcfc:forAllele, or carrying its genotype/GT ordinal. Keys that flatten fixed-width tuples (CIPOS, MEINFO, …) 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 RN/RUS/RUL/RUC/RB tandem-repeat structure.

Those carriers are emitted once per record, after the INFO values, rather than key by key — several of them need more than one key to be well formed. A vcfc:VariantEvent needs both EVENT and EVENTTYPE; a reference block needs both END and POS; a tandem repeat needs RN before its repeat values have a grouping. Where a record does not supply the whole set, the carrier is not emitted and the values stay available as ordinary INFO values, rather than producing a resource that would fail its SHACL shape.

--info-representation raw emits only the opaque vcfc:infoRaw string, and drops the value items with it.

The allele layer is not tied to the INFO representation. Two layers join to <record>/allele/<index>: the structured INFO value items, and the expanded sample layer's per-call vcfc:calledAllele. Either one alone requires the alleles, so the layer is emitted when --info-representation structured or --sample-representation expanded is in force, and left out only for raw + condensed, where nothing references it.

append_header_representation_rdf — structured headers

--header-representation structured (the default) types each ## line with its vocabulary subclass (INFOHeaderLine, ContigHeaderLine, FILTERHeaderLine, ALTHeaderLine, MetaHeaderLine, SampleHeaderLine, PedigreeHeaderLine, AssemblyHeaderLine, PedigreeDBHeaderLine, …) and lifts its attributes into their own properties. The attribute parser respects quoting and backslash escapes, so a Description containing a comma is handled correctly.

Every structured line also gets one ordered vcfc:HeaderAttribute per source attribute, carrying attributeKey, attributeValue and a one-based attributeIndex. This is not optional decoration: vcfc:StructuredHeaderLineShape requires at least one vcfc:hasAttribute, so a line typed as one of the structured subclasses without them is non-conformant. It also keeps implementation-defined attributes queryable, which the dedicated properties alone cannot do.

An unrecognised ## key is still typed as vcfc:StructuredHeaderLine or vcfc:UnstructuredHeaderLine according to its value's form, which keeps it queryable without inventing a subclass the vocabulary does not define. An unstructured line is only typed when it has a value, because vcfc:UnstructuredHeaderLineShape requires one.

##fileDate is emitted as xsd:date when its form is recognisable, and verbatim as a plain literal otherwise — again preserved rather than dropped, and reported by SHACL.

--header-representation basic keeps only the base header-line triples the RML mapping emits. It is the fastest and smallest header form, but it does not produce a SHACL-conformant graph, because the structured lines then carry no attributes and no subclass.

emit_sample_representation — genotypes

Exactly one sample emitter runs per conversion, chosen by --sample-representation. Full treatment in sample-representation-guide.md; the short version:

Mode Shape Growth
expanded (default) one vcfc:SampleCall per record × sample, one vcfc:FormatFieldValue per FORMAT key, and the parsed genotype layer below ≈ variants × samples × FORMAT fields
condensed one reusable vcfc:SampleSet, one vcfc:CohortCallMatrix per call, one vcfc:FormatValueVector per FORMAT key holding a tab-separated vcfc:VCFTextVector ≈ samples + variants × FORMAT fields

Both modes emit the file-level vcfc:SampleSet and its ordered vcfc:VCFSample members, and link them from the #CHROM vcfc:ColumnHeaderLine with vcfc:hasGenotypeColumns. That is one resource per sample column for the whole file, not per record, so it costs nothing at cohort scale — and it gives the expanded profile's SampleCall a reusable identity through vcfc:forSample 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, 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. Every value stays recoverable by decoding a vector against its FORMAT definition and the matrix SampleSet.

A GT value that is not in the vcfc:GenotypeString lexical space produces no Genotype resource at all — the value remains present as the FORMAT field's fieldValue, and inventing a genotype for it would create a resource that fails vcfc:GenotypeShape.

The file declares which shape it carries via vcfc:representationProfile, so a consumer can branch on it before querying. The SPARQL SHACL profile enforces that a graph never mixes the two.

4a. VCF versions

VCF Core ships one SHACL overlay per VCF version (4.1 through 4.5), each scoped by the file's vcfc:fileFormat, plus an optional vcfc:VCF4xFile subclass that activates that version's gate. The conversion follows the version of each input.

Detection is automatic

##fileformat is required to be the first line of a conforming VCF, and src/vcf_as_tsv.sh already lifts it into file_metadata.tsv, so detection costs one small read and needs nothing from you. It happens per input, not per run: a directory may hold files of different versions and each gets its own mapping and emitter behaviour.

--vcf-version overrides it, for a file whose declaration is missing or wrong:

vcf-rdfizer --mode full -i ./cohort.vcf.gz --vcf-version 4.3 -o ./results

The run log states which version was used and how it was chosen — (detected), (forced), or a warning that the declaration was not recognized.

VCF 4.0 is not recognized. VCF Core claims no 4.0 conformance overlay, so a 4.0 file converts with the newest supported rules — nothing representable is dropped — but no vcfc:VCF4xFile class is emitted, so the graph never claims a version gate that cannot be checked. The same applies to a missing or malformed ##fileformat line. All three cases are reported, never assumed silently.

What the version changes

4.1 4.2 4.3 4.4 4.5
Number=R — ✓ ✓ ✓ ✓
Number=P (FORMAT) — — — ✓ ✓
Number=LA/LR/LG/M (FORMAT) — — — — ✓
CIPOS/CIEND one pair per record one pair per record one pair per record one pair per ALT one pair per ALT
CILEN/CICN — — — ✓ ✓
EVENT one per record one per record one per record one per ALT one per ALT
EVENTTYPE — — — ✓ ✓
Local alleles (LAA) — — — — ✓
Base modifications (M/DPM/ADM) — — — — ✓
Leading GT phase indicator — — — ✓ ✓

Three consequences in the emitted graph:

  • vcfc:FieldValueItem linkage. Before 4.4 a CIPOS list is one pair describing the record, so its items get a vcfc:valueIndex and a vcfc:tupleArity but no vcfc:forAllele — there is no single allele to point at, and each ALT allele of a multi-allelic record shares the interval. From 4.4 the list repeats per ALT and each item joins its own allele.
  • Tuple items are always materialized. The version overlays count a tuple key's value items against the ALT count, so CIPOS and friends are decomposed even though 4.4 and 4.5 declare them Number=., which on its own would say "no positional meaning".
  • Families a version does not define are not invented. A LAA column in a 4.4 file, or an M27551C column in a 4.3 file, keeps its raw vcfc:fieldValue and gets no vcfc:LocalAlleleSet or vcfc:BaseModification resource.

What it does not change

The version does not affect parsing leniency. A GT with a leading phase indicator in a 4.1 file is still parsed into a vcfc:Genotype; the version overlay reports it as non-conformant, which is the right division of labour — the converter transcribes, the validator judges. The same holds for a Number=R declaration in a 4.1 file, or a reserved key declared with the wrong arity: the graph records what the file says, and SHACL says whether the file was right.

5. The two helper tables

sample_calls.tsv and sample_format_values.tsv exist so that a mapping can express genotypes in RML if it wants to. Under the shipped mapping they stay header-only and the wrapper streams genotypes directly, because materializing them means one row per variant × sample (× FORMAT key) — the largest intermediate the pipeline can produce.

A custom mapping that consumes them in any other way forces full materialization, and is rejected in --sample-representation condensed, which would otherwise emit both genotype representations into one graph.

6. IRI templates

Every IRI is derived from the source filename and the row counter, so a conversion is deterministic and re-running it produces byte-comparable subjects.

Resource Template
VCFFile file://{SOURCE_FILE}
VCFHeader file://{SOURCE_FILE}#header
HeaderLine file://{SOURCE_FILE}#header/line/{HEADER_INDEX}
VCFRecord file://{SOURCE_FILE}#record/{ROW_ID}
VariantCall file://{SOURCE_FILE}#call/{ROW_ID}
ColumnHeaderLine file://{SOURCE_FILE}#header/columns
HeaderAttribute file://{SOURCE_FILE}#header/line/{HEADER_INDEX}/attribute/{N}
SampleDeclaration (##SAMPLE/##PEDIGREE) file://{SOURCE_FILE}#header/sample/{ID}
InfoFieldValue file://{SOURCE_FILE}#call/{ROW_ID}/info/{KEY}
Allele (REF at 0, ALT from 1) file://{SOURCE_FILE}#record/{ROW_ID}/allele/{N}
FieldValueItem <parent value IRI>/value/{N}
VariantEvent file://{SOURCE_FILE}#record/{ROW_ID}/event/{EVENT}
SampleSet (both profiles) file://{SOURCE_FILE}#samples
VCFSample (both profiles) file://{SOURCE_FILE}#samples/{SAMPLE}
SampleCall (expanded) file://{SOURCE_FILE}#sample/{ROW_ID}/{SAMPLE}
FormatFieldValue (expanded) file://{SOURCE_FILE}#sample/{ROW_ID}/{SAMPLE}/fmt/{KEY}
Genotype (expanded) file://{SOURCE_FILE}#sample/{ROW_ID}/{SAMPLE}/genotype
GenotypeAlleleCall (expanded) …/genotype/call/{N}
PhaseSet (expanded) file://{SOURCE_FILE}#sample/{ROW_ID}/{SAMPLE}/phaseset
CohortCallMatrix (condensed) file://{SOURCE_FILE}#call/{ROW_ID}/matrix
FormatValueVector (condensed) file://{SOURCE_FILE}#call/{ROW_ID}/matrix/fmt/{KEY}

These follow the vcfc:iriTemplate patterns the vocabulary recommends for each class. Where the vocabulary declares no template — the allele-layer, genotype and SV resources — the IRI is minted beneath the resource it belongs to, so the whole graph stays derivable from the filename and the row counter alone.

{SOURCE_FILE} is the basename, not a path, so a graph does not encode where the VCF happened to live. The consequence is that two different VCFs with the same filename mint the same IRIs; if you convert chr1/data.vcf and chr2/data.vcf, their graphs will collide when merged. Rename before converting, or keep the graphs separate.

7. Missing values

The vocabulary's vcfc:missingValuePolicy says a missing token should be "."^^vcfc:Null, and the conversion follows it. preflight_missing_token_conformance reports a plain "." literal as an anomaly, and --strict-conformance promotes that from a report to a failure.

Not every "." is a missing value, and the check excludes the three places the vocabulary requires a bare dot: vcfc:fieldNumber (Number=. is VCF's variable-cardinality arity), vcfc:genotypeString (a fully missing call is literally . or ./.), and the vcfc:attributeValue of a header attribute whose vcfc:attributeKey is Number — the structured-header layer carries each declaration's attributes verbatim, so Number=. appears there too, and the shapes require that value to be xsd:string. A bare dot anywhere else, including on any other header attribute, is still reported.

This used to conflict with the published SHACL shapes, which constrained vcfc:alt to sh:datatype xsd:string and so rejected every REF-only and gVCF-style record. The vocabulary fixed that: vcfc:VCFRecordShape now accepts xsd:string or vcfc:Null for alt, and vcfc:VariantCallShape accepts the full VCF Float lexical space for qual. Emitting the missing token as vcfc:Null is now both the policy and the conformant choice.

CHROM, POS and REF are the three fixed fields VCF 4.5 forbids from being the missing token, which is exactly why they can stay in the RML mapping while ID, ALT, QUAL and FILTER cannot.

8. What conversion does not do

  • No normalization. No left-alignment, no trimming, no multi-allelic splitting. ALT=A,T stays one record with one alt literal. This is deliberate — the graph is a faithful transcription of the file — but it means the graph is not directly joinable with normalized external resources. This is the central problem the data-linking design has to solve.
  • No reference checking. REF is not verified against any genome.
  • No cross-file merging. Each VCF produces its own graph.

Two entries that used to be in this list no longer belong here, because the move to the VCF Core vocabulary added both:

  • Structural variants are modelled. Symbolic ALTs, breakends and * are classified by vcfc:alleleKind, given a vcfc:svType where the identifier is a reserved code, and parsed into their components; END, SVLEN, SVCLAIM, EVENT/EVENTTYPE, the CI* intervals and the tandem-repeat keys all have meaning. The raw literals are still there — the structure is added alongside, not instead.
  • Genotypes are interpreted in the expanded profile: GT is parsed into ordered allele calls with an explicit vcfc:phasingStatus, so | versus / is modelled rather than only preserved. The condensed profile keeps them lexical by design.

See also