Variants and files API

For examples, see file loading, sample queries, and transforms.

Reference genomes

varcode.Genome(ensembl_release: Any = None, *, fasta: Optional[Any] = None, verify: bool = True)

Pyensembl Genome + optional chromosome FASTA.

Two source layers under one object:

• self.fasta — optional chromosome FASTA. When attached, :meth:sequence reads from it directly; :meth:reference_base and :meth:reference_range prefer it before falling back to transcript cDNA.
• Wrapped pyensembl Genome — transcript annotations + cDNA. Always present; provides the fallback for exonic positions when no FASTA is attached.

The two methods have intentionally different fall-through semantics:

• :meth:sequence — chromosome FASTA only. Returns "" when no FASTA is attached. Mirrors the proposed pyensembl.Genome.sequence() shape from openvax/pyensembl#337 so the eventual upstream migration is mechanical.
• :meth:reference_base / :meth:reference_range — tiered. FASTA first when attached; otherwise transcript cDNA. Use these when you want "whatever varcode can tell you" about a position.

PARAMETER	DESCRIPTION
ensembl_release	Release identifier — passes through :func:varcode.reference.infer_genome, so the same shapes accepted by load_vcf(genome=...) work here (int, reference-name string, pyensembl.Genome). Passing another :class:Genome rewraps idempotently (inheriting fasta unless overridden). TYPE: Any DEFAULT: None
fasta	Optional chromosome FASTA: Path string — opened with pyfaidx (only required on this path). pyfaidx.Fasta — used as-is. Any object supporting fa[contig][start:end] and returning a string or an object with .seq. None — no chromosome-level access; features fall back to transcript cDNA. varcode does not take ownership of the FASTA object — when the caller passes a pre-opened pyfaidx.Fasta, the caller is responsible for its lifetime. Closing the FASTA after passing it to Genome will cause subsequent lookups to fail. TYPE: Optional[Any] DEFAULT: None
verify	When True (default) and fasta is provided, spot-check a few exonic positions against pyensembl's transcript cDNA to catch mislabeled FASTAs (e.g. GRCh37 attached to a GRCh38 release). Iterates self.transcripts() which on a fresh process may trigger lazy DB construction; pass verify=False to defer that cost. TYPE: bool DEFAULT: True

Examples:

>>> import varcode
>>> g = varcode.Genome(81, fasta="/path/to/GRCh38.fa")
>>> vc = varcode.load_vcf("tumor.vcf", genome=g)
>>> g.sequence("7", 117_480_000, 117_480_050)
'AC...'

Notes

Returned sequence is always uppercase. Soft-masked (lowercase) repeat annotations from the FASTA are dropped — varcode treats all bases uniformly. Callers that need the soft-masking signal should read the FASTA directly.

Equality and hashing fall back to object identity — two Genome instances wrapping the same pyensembl release compare unequal. This is intentional: a wrapper with an attached FASTA and one without are different objects in any sense that matters for varcode features, even if they share the same reference_name.

Source code in varcode/genome.py

def __init__(
        self,
        ensembl_release: Any = None,
        *,
        fasta: Optional[Any] = None,
        verify: bool = True):
    fasta_was_provided = fasta is not None
    if isinstance(ensembl_release, Genome):
        # Idempotent rewrap. Inherits ._ensembl always; inherits
        # .fasta unless the caller is providing a fresh one.
        self._ensembl = ensembl_release._ensembl
        if fasta_was_provided:
            self.fasta = _resolve_fasta(fasta, self._ensembl)
        else:
            self.fasta = ensembl_release.fasta
    else:
        if ensembl_release is None:
            raise ValueError(
                "varcode.Genome requires an ensembl_release argument "
                "(int release number, reference-name string, or a "
                "pyensembl.Genome / varcode.Genome instance).")
        self._ensembl, _ = infer_genome(ensembl_release)
        self.fasta = (_resolve_fasta(fasta, self._ensembl)
                      if fasta_was_provided else None)

    # Verify only when the caller freshly provided a FASTA.
    # Inheriting a verified FASTA on rewrap doesn't need a re-check.
    if (verify and fasta_was_provided
            and self.fasta is not None
            and self._ensembl is not None):
        _verify_fasta_against_transcripts(self._ensembl, self.fasta)

    # Per-instance warn-once for missing-reference signals from
    # consumers (cryptic_exons, etc.). Bound to this Genome's
    # lifetime — no module-level cache to invalidate.
    self._missing_reference_warned = False

__getattr__(name)

Delegate everything else to the wrapped pyensembl Genome.

Called only when normal attribute lookup fails, so explicit attributes (_ensembl, fasta, _missing_reference_warned) still win. Underscore-prefixed names are not delegated — this prevents infinite recursion when Python's copy / pickle machinery probes for __deepcopy__ / __reduce__ / __getstate__ etc., and keeps the wrapper's private state from accidentally pulling private state from the wrapped object.

Source code in varcode/genome.py

def __getattr__(self, name):
    """Delegate everything else to the wrapped pyensembl Genome.

    Called only when normal attribute lookup fails, so explicit
    attributes (``_ensembl``, ``fasta``, ``_missing_reference_warned``)
    still win. Underscore-prefixed names are not delegated — this
    prevents infinite recursion when Python's copy / pickle
    machinery probes for ``__deepcopy__`` / ``__reduce__`` /
    ``__getstate__`` etc., and keeps the wrapper's private state
    from accidentally pulling private state from the wrapped
    object.
    """
    if name.startswith("_"):
        raise AttributeError(name)
    ensembl = self.__dict__.get("_ensembl")
    if ensembl is None:
        raise AttributeError(
            "varcode.Genome has no wrapped pyensembl Genome; "
            "cannot resolve attribute %r" % name)
    return getattr(ensembl, name)

__dir__()

Include attributes of the wrapped pyensembl Genome so dir(genome) and IDE auto-completion surface the full delegated API, not just the wrapper's own attributes.

Source code in varcode/genome.py

def __dir__(self):
    """Include attributes of the wrapped pyensembl Genome so
    ``dir(genome)`` and IDE auto-completion surface the full
    delegated API, not just the wrapper's own attributes."""
    own = set(super().__dir__())
    ensembl = self.__dict__.get("_ensembl")
    if ensembl is not None:
        own.update(dir(ensembl))
    return sorted(own)

sequence(contig: str, start: int, end: int) -> str

Chromosome FASTA sequence on the + strand.

1-based inclusive coordinates. Returns "" when no FASTA is attached or the contig / range isn't covered.

FASTA-only — does not fall back to transcript cDNA. Use :meth:reference_base / :meth:reference_range for the tiered lookup. The split is deliberate: this method mirrors the proposed pyensembl.Genome.sequence() API so callers who want raw chromosome bases (and the upstream migration path) can use it unambiguously.

Source code in varcode/genome.py

def sequence(self, contig: str, start: int, end: int) -> str:
    """Chromosome FASTA sequence on the ``+`` strand.

    1-based inclusive coordinates. Returns ``""`` when no FASTA is
    attached or the contig / range isn't covered.

    FASTA-only — does **not** fall back to transcript cDNA. Use
    :meth:`reference_base` / :meth:`reference_range` for the
    tiered lookup. The split is deliberate: this method mirrors
    the proposed ``pyensembl.Genome.sequence()`` API so callers
    who want raw chromosome bases (and the upstream migration
    path) can use it unambiguously.
    """
    if self.fasta is None:
        return ""
    return _fasta_range(self.fasta, contig, start, end)

reference_base(contig: str, position: int) -> str

Tiered + strand base lookup (FASTA → transcript cDNA → "").

Convenience method delegating to :func:varcode.genome_sequence.reference_base.

Source code in varcode/genome.py

def reference_base(self, contig: str, position: int) -> str:
    """Tiered ``+`` strand base lookup (FASTA → transcript cDNA → ``""``).

    Convenience method delegating to
    :func:`varcode.genome_sequence.reference_base`.
    """
    return _reference_base_lookup(self, contig, position)

reference_range(contig: str, start: int, end: int) -> str

Tiered + strand range lookup. All-or-nothing — returns "" if any position in the range is uncovered by the chosen source.

Source code in varcode/genome.py

def reference_range(self, contig: str, start: int, end: int) -> str:
    """Tiered ``+`` strand range lookup. All-or-nothing — returns
    ``""`` if any position in the range is uncovered by the chosen
    source."""
    return _reference_range_lookup(self, contig, start, end)

Variants

varcode.Variant(contig, start, ref, alt, genome=None, ensembl=None, allow_extended_nucleotides=False, normalize_contig_names=True, convert_ucsc_contig_names=None)

Bases: Serializable

Construct a Variant object.

PARAMETER	DESCRIPTION
contig	Chromosome that this variant is on TYPE: str
start	1-based position on the chromosome of first reference nucleotide TYPE: int
ref	Reference nucleotide(s) TYPE: str
alt	Alternate nucleotide(s) TYPE: str
genome	Name of reference genome, Ensembl release number, or object derived from pyensembl.Genome. Default to latest available release of GRCh38 TYPE: Genome, EnsemblRelease, or str, or int DEFAULT: None
ensembl	Previous name used instead of 'genome', the two arguments should be mutually exclusive. TYPE: Genome, EnsemblRelease, or str, or int (DEPRECATED) DEFAULT: None
allow_extended_nucleotides	Extended nucleotides include 'Y' for pyrimidies or 'N' for any base TYPE: bool DEFAULT: False
normalize_contig_names	By default the contig name will be normalized by converting integers to strings (e.g. 1 -> "1"), and converting any letters after "chr" to uppercase (e.g. "chrx" -> "chrX"). If you don't want this behavior then pass normalize_contig_name=False. TYPE: bool DEFAULT: True
convert_ucsc_contig_names	Setting this argument to True causes UCSC chromosome names to be coverted, such as "chr1" to "1". If the default value (None) is used then it defaults to whether or not a UCSC genome was pass in for the 'genome' argument. TYPE: bool DEFAULT: None

Source code in varcode/variant.py

def __init__(
        self,
        contig,
        start,
        ref,
        alt,
        genome=None,
        ensembl=None,
        allow_extended_nucleotides=False,
        normalize_contig_names=True,
        convert_ucsc_contig_names=None):
    """
    Construct a Variant object.

    Parameters
    ----------
    contig : str
        Chromosome that this variant is on

    start : int
        1-based position on the chromosome of first reference nucleotide

    ref : str
        Reference nucleotide(s)

    alt : str
        Alternate nucleotide(s)

    genome : Genome, EnsemblRelease, or str, or int
        Name of reference genome, Ensembl release number, or object
        derived from pyensembl.Genome. Default to latest available release
        of GRCh38

    ensembl : Genome, EnsemblRelease, or str, or int (DEPRECATED)
        Previous name used instead of 'genome', the two arguments should
        be mutually exclusive.

    allow_extended_nucleotides : bool
        Extended nucleotides include 'Y' for pyrimidies or 'N' for any base

    normalize_contig_names : bool
        By default the contig name will be normalized by converting integers
        to strings (e.g. 1 -> "1"), and converting any letters after "chr"
        to uppercase (e.g. "chrx" -> "chrX"). If you don't want
        this behavior then pass normalize_contig_name=False.

    convert_ucsc_contig_names : bool, optional
        Setting this argument to True causes UCSC chromosome names to be
        coverted, such as "chr1" to "1". If the default value (None) is used
        then it defaults to whether or not a UCSC genome was pass in for
        the 'genome' argument.
    """

    # first initialize the fields we use to cache lists of overlapping
    # pyensembl Gene and Transcript objects, or their properties such
    # as names/IDs
    self._genes = None
    self._transcripts = None
    self._gene_ids = None
    self._gene_names = None

    # Provenance: when this variant was produced by a
    # ``varcode.transforms`` operation (e.g. ``pair_breakends``,
    # ``combine_cis_snvs``), the source rows are recorded here.
    # Empty tuple for variants loaded directly from a VCF/MAF.
    # Not part of hash/equality — identity is still
    # (contig, start, ref, alt, reference_name).
    self.source_variants: tuple = ()

    # store the options which affect how properties of this variant
    # may be changed/transformed
    self.normalize_contig_names = normalize_contig_names
    self.allow_extended_nucleotides = allow_extended_nucleotides

    # if genome not specified, try the old name 'ensembl'
    # if ensembl is also None, then default to "GRCh38"
    if genome is None and ensembl is None:
        genome = "GRCh38"
    elif genome is None:
        genome = ensembl


    # user might supply Ensembl release as an integer, reference name,
    # or pyensembl.Genome object
    self.original_genome = genome
    self.genome, self.original_genome_was_ucsc = infer_genome(genome)

    self.reference_name = self.genome.reference_name
    if self.original_genome_was_ucsc:
        self.original_reference_name = ensembl_to_ucsc_reference_names[
            self.reference_name]
    else:
        self.original_reference_name = self.reference_name

    self.original_contig = contig

    if convert_ucsc_contig_names is None:
        self.convert_ucsc_contig_names = self.original_genome_was_ucsc
    else:
        self.convert_ucsc_contig_names = convert_ucsc_contig_names

    self.contig = self._normalize_contig_name(contig)

    if ref != alt and ref in STANDARD_NUCLEOTIDES and alt in STANDARD_NUCLEOTIDES:
        # Optimization for common case.
        self.original_ref = self.ref = ref
        self.original_alt = self.alt = alt
        self.original_start = self.start = self.end = int(start)
        return

    # we want to preserve the ref/alt/pos both as they appeared in the
    # original VCF or MAF file but also normalize variants to get rid
    # of shared prefixes/suffixes between the ref and alt nucleotide
    # strings e.g. g.10 CTT>T can be normalized into g.10delCT
    #
    # The normalized variant properties go into fields
    #    Variant.{original_ref, original_alt, original_pos}
    # whereas the trimmed fields are:
    #    Variant.{ref, alt, start, end}

    # the original entries must preserve the number of nucleotides in
    # ref and alt but we still want to normalize e.g. '-' and '.' into ''
    self.original_ref = normalize_nucleotide_string(
        ref,
        allow_extended_nucleotides=allow_extended_nucleotides)
    self.original_alt = normalize_nucleotide_string(
        alt,
        allow_extended_nucleotides=allow_extended_nucleotides)
    self.original_start = int(start)

    # normalize the variant by trimming any shared prefix or suffix
    # between ref and alt nucleotide sequences and then
    # offset the variant position in a strand-dependent manner
    (trimmed_ref, trimmed_alt, prefix, _) = \
        trim_shared_flanking_strings(self.original_ref, self.original_alt)

    self.ref = trimmed_ref
    self.alt = trimmed_alt

    if len(trimmed_ref) == 0:
        # insertions must be treated differently since the meaning of a
        # position for an insertion is:
        #   "insert the alt nucleotides after this position"
        #
        # Aside: what if both trimmed ref and alt strings are empty?
        # This means we had a "null" variant, probably from a VCF
        # generated by force-calling mutations which weren't actually
        # found in the sample.
        # Null variants are interepted as inserting zero nucleotides
        # after the whole reference sequence.
        #
        # Start and end both are base-1 nucleotide position before
        # insertion.
        self.start = self.original_start + max(0, len(prefix) - 1)
        self.end = self.start
    else:
        # for substitutions and deletions the [start:end] interval is
        # an inclusive selection of reference nucleotides
        self.start = self.original_start + len(prefix)
        self.end = self.start + len(trimmed_ref) - 1

ensembl property

Deprecated alias for Variant.genome

RETURNS	DESCRIPTION
Genome

trimmed_ref property

Eventually the field Variant.ref will store the reference nucleotides as given in a VCF or MAF and trimming of any shared prefix/suffix between ref and alt will be done via the properties trimmed_ref and trimmed_alt.

trimmed_alt property

Eventually the field Variant.ref will store the reference nucleotides as given in a VCF or MAF and trimming of any shared prefix/suffix between ref and alt will be done via the properties trimmed_ref and trimmed_alt.

trimmed_base1_start property

Currently the field Variant.start carries the base-1 starting position adjusted by trimming any shared prefix between Variant.ref and Variant.alt. Eventually this trimming should be done more explicitly via trimmed_* properties.

trimmed_base1_end property

Currently the field Variant.end carries the base-1 "last" position of this variant, adjusted by trimming any shared suffix between Variant.ref and Variant.alt. Eventually this trimming should be done more explicitly via trimmed_* properties.

short_description property

HGVS nomenclature for genomic variants More info: http://www.hgvs.org/mutnomen/

coding_transcripts property

Protein coding transcripts

genes property

Return Gene object for all genes which overlap this variant.

gene_ids property

Return IDs of all genes which overlap this variant. Calling this method is significantly cheaper than calling Variant.genes(), which has to issue many more queries to construct each Gene object.

gene_names property

Return names of all genes which overlap this variant. Calling this method is significantly cheaper than calling Variant.genes(), which has to issue many more queries to construct each Gene object.

coding_genes property

Protein coding transcripts

is_insertion property

Does this variant represent the insertion of nucleotides into the reference genome?

is_deletion property

Does this variant represent the deletion of nucleotides from the reference genome?

is_indel property

Is this variant either an insertion or deletion?

is_snv property

Is the variant a single nucleotide variant

is_transition property

Is this variant and pyrimidine to pyrimidine change or purine to purine change

is_transversion property

Is this variant a pyrimidine to purine change or vice versa

__lt__(other)

Variants are ordered by locus.

Source code in varcode/variant.py

def __lt__(self, other):
    """
    Variants are ordered by locus.
    """
    require_instance(other, Variant, name="variant")
    if self.contig == other.contig:
        return self.start < other.start
    return self.contig < other.contig

to_dict()

We want the original values (un-normalized) field values while serializing since normalization will happen in init.

Source code in varcode/variant.py

def to_dict(self):
    """
    We want the original values (un-normalized) field values while
    serializing since normalization will happen in __init__.
    """
    return dict(
        contig=self.original_contig,
        start=self.original_start,
        ref=self.original_ref,
        alt=self.original_alt,
        genome=self.original_genome,
        allow_extended_nucleotides=self.allow_extended_nucleotides,
        normalize_contig_names=self.normalize_contig_names,
        convert_ucsc_contig_names=self.convert_ucsc_contig_names)

effects(raise_on_error=True, annotator=None, phase_resolver=None, rna_resolver=None, germline=None)

Predict the variant's effects on overlapping transcripts.

Splice-disrupting effects are always wrapped in a :class:varcode.splice_outcomes.SpliceOutcomeSet carrying the candidate mechanisms (normal splicing, exon skipping, intron retention, cryptic splice). Mechanism candidates are computed lazily, so the wrap is cheap. See openvax/varcode#262, #391.

PARAMETER	DESCRIPTION
raise_on_error	If True, raise on annotation errors; if False, capture per-transcript errors as Failure effects. Failed initial gene/transcript lookups return an empty EffectCollection and log the error. TYPE: bool DEFAULT: True
annotator	Per-call annotator override. None uses the currently configured default ("fast" today; swappable via :func:varcode.set_default_annotator or :func:varcode.use_annotator). String names are resolved against the registry. See openvax/varcode#271. TYPE: str, EffectAnnotator, or None DEFAULT: None
phase_resolver	Optional phase-evidence source (typically a :class:~varcode.phasing.MolecularPhaseResolver wrapping an upstream RNA-phasing tool, or :class:~varcode.phasing.VCFPhaseResolver). When provided and the resolver has an observed mutant transcript for (self, transcript), the returned effect's mutant_transcript is populated with that observed transcript — the protein is the protein actually observed in RNA rather than one inferred from the reference. See openvax/varcode#269. TYPE: resolver object or None DEFAULT: None
rna_resolver	Optional RNA-observed-outcome source. When provided, any :class:~varcode.MultiOutcomeEffect in the result is refined with observed candidates from the resolver. Splice mechanism sets are replaced by reconciled copies that track added RNA candidates and excluded DNA predictions; other multi-outcome effects keep the additive candidate behavior. See openvax/varcode#259. TYPE: RNAEvidenceResolver or None DEFAULT: None
germline	Optional patient-germline context. When non-empty, every per-transcript effect is computed against the patient's germline-applied transcript instead of the reference. Codons / splice signals where germline overlaps the somatic and phase is unknown produce a :class:PhaseCandidateSet with one outcome per haplotype hypothesis; LOH at germline het positions sets effect.is_loh = True. None or :meth:GermlineContext.empty falls through to today's reference-relative behaviour byte-identically. See openvax/varcode#268. TYPE: GermlineContext or None DEFAULT: None

Source code in varcode/variant.py

def effects(
        self, raise_on_error=True,
        annotator=None, phase_resolver=None, rna_resolver=None,
        germline=None):
    """Predict the variant's effects on overlapping transcripts.

    Splice-disrupting effects are always wrapped in a
    :class:`varcode.splice_outcomes.SpliceOutcomeSet` carrying the
    candidate mechanisms (normal splicing, exon skipping, intron
    retention, cryptic splice). Mechanism candidates are computed
    lazily, so the wrap is cheap. See openvax/varcode#262, #391.

    Parameters
    ----------
    raise_on_error : bool
        If True, raise on annotation errors; if False, capture
        per-transcript errors as Failure effects. Failed initial
        gene/transcript lookups return an empty EffectCollection
        and log the error.

    annotator : str, EffectAnnotator, or None
        Per-call annotator override. ``None`` uses the currently
        configured default (``"fast"`` today; swappable via
        :func:`varcode.set_default_annotator` or
        :func:`varcode.use_annotator`). String names are resolved
        against the registry. See openvax/varcode#271.

    phase_resolver : resolver object or None
        Optional phase-evidence source (typically a
        :class:`~varcode.phasing.MolecularPhaseResolver` wrapping an
        upstream RNA-phasing tool, or
        :class:`~varcode.phasing.VCFPhaseResolver`). When provided
        and the resolver has an observed mutant transcript for
        ``(self, transcript)``, the returned effect's
        ``mutant_transcript`` is populated with that observed
        transcript — the protein is the protein actually observed
        in RNA rather than one inferred from the reference. See
        openvax/varcode#269.

    rna_resolver : RNAEvidenceResolver or None
        Optional RNA-observed-outcome source. When provided, any
        :class:`~varcode.MultiOutcomeEffect` in the result is
        refined with observed candidates from the resolver. Splice
        mechanism sets are replaced by reconciled copies that track
        added RNA candidates and excluded DNA predictions; other
        multi-outcome effects keep the additive candidate behavior.
        See openvax/varcode#259.

    germline : GermlineContext or None
        Optional patient-germline context. When non-empty, every
        per-transcript effect is computed against the patient's
        germline-applied transcript instead of the reference.
        Codons / splice signals where germline overlaps the
        somatic and phase is unknown produce a
        :class:`PhaseCandidateSet` with one outcome per
        haplotype hypothesis; LOH at germline het positions sets
        ``effect.is_loh = True``. ``None`` or
        :meth:`GermlineContext.empty` falls through to today's
        reference-relative behaviour byte-identically. See
        openvax/varcode#268.
    """
    effects = predict_variant_effects(
        variant=self,
        raise_on_error=raise_on_error,
        annotator=annotator,
        germline=germline,
        phase_resolver=phase_resolver,
    )
    if phase_resolver is not None:
        from .phasing import apply_phase_resolver_to_effects
        apply_phase_resolver_to_effects(effects, phase_resolver)
    if rna_resolver is not None:
        from .rna_evidence import apply_rna_evidence_to_effects
        apply_rna_evidence_to_effects(effects, rna_resolver)
    return effects

effect_on_transcript(transcript, annotator=None, germline=None, phase_resolver=None)

Annotate one transcript using the same selection as :meth:effects.

Source code in varcode/variant.py

def effect_on_transcript(
        self, transcript, annotator=None, germline=None, phase_resolver=None):
    """Annotate one transcript using the same selection as :meth:`effects`."""
    return predict_variant_effect_on_transcript(
        self, transcript, annotator=annotator, germline=germline,
        phase_resolver=phase_resolver)

clone_without_ucsc_data()

Clone this variant but discarding the original format of its genome and contig: useful if we want to mix hg19 and GRCh37 variants.

RETURNS	DESCRIPTION
Variant

Source code in varcode/variant.py

def clone_without_ucsc_data(self):
    """
    Clone this variant but discarding the original format of its genome
    and contig: useful if we want to mix hg19 and GRCh37 variants.

    Returns
    -------
    Variant
    """
    return Variant(
        contig=self.contig,
        start=self.original_start,
        ref=self.original_ref,
        alt=self.original_alt,
        genome=self.genome,
        allow_extended_nucleotides=self.allow_extended_nucleotides,
        normalize_contig_names=self.normalize_contig_names,
        convert_ucsc_contig_names=False)

varcode.VariantCollection(variants, distinct=True, sort_key=variant_ascending_position_sort_key, sources=None, source_to_metadata_dict={})

Bases: Collection

Construct a VariantCollection from a list of Variant records.

PARAMETER	DESCRIPTION
variants	Variant objects contained in this VariantCollection TYPE: iterable
distinct	Don't keep repeated variants TYPE: bool DEFAULT: True
sort_key	TYPE: callable DEFAULT: variant_ascending_position_sort_key
sources	Optional set of source names, may be larger than those for which we have metadata dictionaries. TYPE: set DEFAULT: None
source_to_metadata_dict	Dictionary mapping each source name (e.g. VCF path) to a dictionary from metadata attributes to values. TYPE: dict DEFAULT: {}

Source code in varcode/variant_collection.py

def __init__(
        self,
        variants,
        distinct=True,
        sort_key=variant_ascending_position_sort_key,
        sources=None,
        source_to_metadata_dict={}):
    """
    Construct a VariantCollection from a list of Variant records.

    Parameters
    ----------
    variants : iterable
        Variant objects contained in this VariantCollection

    distinct : bool
        Don't keep repeated variants

    sort_key : callable

    sources : set
        Optional set of source names, may be larger than those for
        which we have metadata dictionaries.

    source_to_metadata_dict : dict
        Dictionary mapping each source name (e.g. VCF path) to a dictionary
        from metadata attributes to values.
    """
    self.source_to_metadata_dict = source_to_metadata_dict
    if sources is None:
        sources = set(source_to_metadata_dict.keys())
    if any(source not in sources for source in source_to_metadata_dict.keys()):
        raise ValueError(
            "Mismatch between sources=%s and keys of source_to_metadata_dict=%s" % (
                sources,
                set(source_to_metadata_dict.keys())))
    Collection.__init__(
        self,
        elements=variants,
        distinct=distinct,
        sort_key=sort_key,
        sources=sources)
    # Keep self.variants in sync with the Collection's post-sort,
    # post-dedup elements so that iterating `vc` and reading
    # `vc.variants` produce the same order.  See openvax/varcode#220.
    self.variants = self.elements

metadata property

The most common usage of a VariantCollection is loading a single VCF, in which case it's annoying to have to always specify that path when accessing metadata fields. This property is meant to both maintain backward compatibility with old versions of Varcode and make the common case easier.

samples property

Sorted list of sample names present in the collection's sample_info metadata (empty if no VCFs with sample columns were loaded).

to_dict()

Since Collection.to_dict() returns a state dictionary with an 'elements' field we have to rename it to 'variants'.

Source code in varcode/variant_collection.py

def to_dict(self):
    """
    Since Collection.to_dict() returns a state dictionary with an
    'elements' field we have to rename it to 'variants'.
    """
    return dict(
        variants=self.variants,
        distinct=self.distinct,
        sort_key=self.sort_key,
        sources=self.sources,
        source_to_metadata_dict=self.source_to_metadata_dict)

clone_with_new_elements(new_elements)

Create another VariantCollection of the same class and with same state (including metadata) but possibly different entries.

Warning: metadata is a dictionary keyed by variants. This method leaves that dictionary as-is, which may result in extraneous entries or missing entries.

Source code in varcode/variant_collection.py

def clone_with_new_elements(self, new_elements):
    """
    Create another VariantCollection of the same class and with
    same state (including metadata) but possibly different entries.

    Warning: metadata is a dictionary keyed by variants. This method
    leaves that dictionary as-is, which may result in extraneous entries
    or missing entries.
    """
    kwargs = self.to_dict()
    kwargs["variants"] = new_elements
    return self.from_dict(kwargs)

effects(raise_on_error=True, annotator=None, phase_resolver=None, rna_resolver=None, germline=None, validate_reference=True)

Splice-disrupting effects are always wrapped in a :class:varcode.splice_outcomes.SpliceOutcomeSet carrying the candidate mechanisms (always-on as of varcode 6.0). Mechanism candidates are computed lazily — only the cheap NormalSplicing candidate is built eagerly. See openvax/varcode#262, #391.

PARAMETER	DESCRIPTION
raise_on_error	If exception is raised while determining effect of variant on a transcript, should it be raised? This default is True, meaning errors result in raised exceptions, otherwise they are only logged. TYPE: bool DEFAULT: True
annotator	Per-call annotator override applied to every variant in the collection. See :meth:Variant.effects and openvax/varcode#271. TYPE: str, EffectAnnotator, or None DEFAULT: None
phase_resolver	Optional phase-evidence source (e.g. a :class:~varcode.phasing.MolecularPhaseResolver wrapping an upstream RNA-phasing tool, or :class:~varcode.phasing.VCFPhaseResolver). When provided, any effect whose (variant, transcript) is covered by an observed mutant transcript has its mutant_transcript populated with that observed transcript. See openvax/varcode#269. TYPE: resolver object or None DEFAULT: None
rna_resolver	Optional RNA-observed-outcome source. When provided, any :class:~varcode.MultiOutcomeEffect in the result is refined with observed candidates from the resolver. Splice mechanism sets are reconciled into replacement sets; other multi-outcome effects keep the additive candidate behavior. See openvax/varcode#259. TYPE: RNAEvidenceResolver or None DEFAULT: None
germline	Optional patient-germline context. When non-empty, every per-transcript effect is computed against the patient's germline-applied transcript instead of the reference. See :meth:Variant.effects and openvax/varcode#268. TYPE: GermlineContext or None DEFAULT: None
validate_reference	Cross-check that the germline context's reference build matches this collection's reference build before running annotation. Hard error on mismatch. Set to False if you've already lifted over and know the builds agree. Ignored when germline is None. TYPE: bool DEFAULT: True

Source code in varcode/variant_collection.py

def effects(
        self, raise_on_error=True,
        annotator=None, phase_resolver=None, rna_resolver=None,
        germline=None, validate_reference=True):
    """
    Splice-disrupting effects are always wrapped in a
    :class:`varcode.splice_outcomes.SpliceOutcomeSet` carrying the
    candidate mechanisms (always-on as of varcode 6.0). Mechanism
    candidates are computed lazily — only the cheap
    ``NormalSplicing`` candidate is built eagerly. See
    openvax/varcode#262, #391.

    Parameters
    ----------
    raise_on_error : bool, optional
        If exception is raised while determining effect of variant on a
        transcript, should it be raised? This default is True, meaning
        errors result in raised exceptions, otherwise they are only logged.

    annotator : str, EffectAnnotator, or None
        Per-call annotator override applied to every variant in
        the collection. See :meth:`Variant.effects` and
        openvax/varcode#271.

    phase_resolver : resolver object or None
        Optional phase-evidence source (e.g. a
        :class:`~varcode.phasing.MolecularPhaseResolver` wrapping an
        upstream RNA-phasing tool, or
        :class:`~varcode.phasing.VCFPhaseResolver`). When provided,
        any effect whose ``(variant, transcript)`` is covered by an
        observed mutant transcript has its ``mutant_transcript``
        populated with that observed transcript. See
        openvax/varcode#269.

    rna_resolver : RNAEvidenceResolver or None
        Optional RNA-observed-outcome source. When provided, any
        :class:`~varcode.MultiOutcomeEffect` in the result is
        refined with observed candidates from the resolver. Splice
        mechanism sets are reconciled into replacement sets; other
        multi-outcome effects keep the additive candidate behavior.
        See openvax/varcode#259.

    germline : GermlineContext or None
        Optional patient-germline context. When non-empty, every
        per-transcript effect is computed against the patient's
        germline-applied transcript instead of the reference.
        See :meth:`Variant.effects` and openvax/varcode#268.

    validate_reference : bool, default True
        Cross-check that the germline context's reference build
        matches this collection's reference build before running
        annotation. Hard error on mismatch. Set to False if
        you've already lifted over and know the builds agree.
        Ignored when ``germline`` is ``None``.
    """
    from datetime import datetime, timezone

    from .annotators.registry import resolve_annotator
    from .phasing import (
        apply_phase_resolver_to_effects,
        build_haplotype_effects,
    )

    # Pre-flight: cross-VCF build mismatch surfaces here as one
    # readable error rather than per-variant ReferenceMismatchError
    # noise. The empty/None context is a no-op.
    if germline:
        germline.validate_against(
            self, validate_reference=validate_reference)

    annotator_instance = resolve_annotator(annotator)
    per_variant = [
        effect
        for variant in self
        for effect in variant.effects(
            raise_on_error=raise_on_error,
            annotator=annotator,
            germline=germline,
            phase_resolver=phase_resolver,
        )
    ]
    if phase_resolver is not None:
        apply_phase_resolver_to_effects(per_variant, phase_resolver)
        # Joint cis-variant effects sit alongside the per-variant
        # ones (#269). Consumers pick whichever granularity they
        # need; top-priority sort still reflects the highest
        # individual effect severity.
        per_variant.extend(build_haplotype_effects(
            self, per_variant, phase_resolver))
    if rna_resolver is not None:
        from .rna_evidence import apply_rna_evidence_to_effects
        apply_rna_evidence_to_effects(per_variant, rna_resolver)
    return EffectCollection(per_variant,
        annotator=getattr(annotator_instance, "name", None),
        annotator_version=getattr(annotator_instance, "version", None),
        annotated_at=datetime.now(timezone.utc).isoformat(
            timespec="seconds"),
    )

reference_names()

All distinct reference names used by Variants in this collection.

RETURNS	DESCRIPTION
set of str

Source code in varcode/variant_collection.py

@memoize
def reference_names(self):
    """
    All distinct reference names used by Variants in this
    collection.

    Returns
    -------
    set of str
    """
    return {variant.reference_name for variant in self}

original_reference_names()

Similar to reference_names but preserves UCSC references, so that a variant collection derived from an hg19 VCF would return {"hg19"} instead of {"GRCh37"}.

RETURNS	DESCRIPTION
set of str

Source code in varcode/variant_collection.py

@memoize
def original_reference_names(self):
    """
    Similar to `reference_names` but preserves UCSC references,
    so that a variant collection derived from an hg19 VCF would
    return {"hg19"} instead of {"GRCh37"}.

    Returns
    -------
    set of str
    """
    return {variant.original_reference_name for variant in self}

groupby_gene_name()

Group variants by the gene names they overlap, which may put each variant in multiple groups.

Source code in varcode/variant_collection.py

def groupby_gene_name(self):
    """
    Group variants by the gene names they overlap, which may put each
    variant in multiple groups.
    """
    return self.multi_groupby(lambda x: x.gene_names)

gene_counts()

Returns number of elements overlapping each gene name. Expects the derived class (VariantCollection or EffectCollection) to have an implementation of groupby_gene_name.

Source code in varcode/variant_collection.py

def gene_counts(self):
    """
    Returns number of elements overlapping each gene name. Expects the
    derived class (VariantCollection or EffectCollection) to have
    an implementation of groupby_gene_name.
    """
    return {
        gene_name: len(group)
        for (gene_name, group)
        in self.groupby_gene_name().items()
    }

filter_by_transcript_expression(transcript_expression_dict, min_expression_value=0.0)

Filters variants down to those which have overlap a transcript whose expression value in the transcript_expression_dict argument is greater than min_expression_value.

PARAMETER	DESCRIPTION
transcript_expression_dict	Dictionary mapping Ensembl transcript IDs to expression estimates (either FPKM or TPM) TYPE: dict
min_expression_value	Threshold above which we'll keep an effect in the result collection TYPE: float DEFAULT: 0.0

Source code in varcode/variant_collection.py

def filter_by_transcript_expression(
        self,
        transcript_expression_dict,
        min_expression_value=0.0):
    """
    Filters variants down to those which have overlap a transcript whose
    expression value in the transcript_expression_dict argument is greater
    than min_expression_value.

    Parameters
    ----------
    transcript_expression_dict : dict
        Dictionary mapping Ensembl transcript IDs to expression estimates
        (either FPKM or TPM)

    min_expression_value : float
        Threshold above which we'll keep an effect in the result collection
    """
    return self.filter_any_above_threshold(
        multi_key_fn=lambda variant: variant.transcript_ids,
        value_dict=transcript_expression_dict,
        threshold=min_expression_value)

filter_by_gene_expression(gene_expression_dict, min_expression_value=0.0)

Filters variants down to those which have overlap a gene whose expression value in the transcript_expression_dict argument is greater than min_expression_value.

PARAMETER	DESCRIPTION
gene_expression_dict	Dictionary mapping Ensembl gene IDs to expression estimates (either FPKM or TPM) TYPE: dict
min_expression_value	Threshold above which we'll keep an effect in the result collection TYPE: float DEFAULT: 0.0

Source code in varcode/variant_collection.py

def filter_by_gene_expression(
        self,
        gene_expression_dict,
        min_expression_value=0.0):
    """
    Filters variants down to those which have overlap a gene whose
    expression value in the transcript_expression_dict argument is greater
    than min_expression_value.

    Parameters
    ----------
    gene_expression_dict : dict
        Dictionary mapping Ensembl gene IDs to expression estimates
        (either FPKM or TPM)

    min_expression_value : float
        Threshold above which we'll keep an effect in the result collection
    """
    return self.filter_any_above_threshold(
        multi_key_fn=lambda effect: effect.gene_ids,
        value_dict=gene_expression_dict,
        threshold=min_expression_value)

exactly_equal(other)

Comparison between VariantCollection instances that takes into account the info field of Variant instances.

RETURNS	DESCRIPTION
True if the variants in this collection equal the variants in the other
collection. The Variant.info fields are included in the comparison.

Source code in varcode/variant_collection.py

def exactly_equal(self, other):
    '''
    Comparison between VariantCollection instances that takes into account
    the info field of Variant instances.

    Returns
    ----------
    True if the variants in this collection equal the variants in the other
    collection. The Variant.info fields are included in the comparison.
    '''
    return (
        self.__class__ == other.__class__ and
        len(self) == len(other) and
        all(x.exactly_equal(y) for (x, y) in zip(self, other)))

union(*others, **kwargs)

Returns the union of variants in a several VariantCollection objects.

Source code in varcode/variant_collection.py

def union(self, *others, **kwargs):
    """
    Returns the union of variants in a several VariantCollection objects.
    """
    return self._combine_variant_collections(
        combine_fn=set.union,
        variant_collections=(self,) + others,
        kwargs=kwargs)

intersection(*others, **kwargs)

Returns the intersection of variants in several VariantCollection objects.

Source code in varcode/variant_collection.py

def intersection(self, *others, **kwargs):
    """
    Returns the intersection of variants in several VariantCollection objects.
    """
    return self._combine_variant_collections(
        combine_fn=set.intersection,
        variant_collections=(self,) + others,
        kwargs=kwargs)

difference(*others, **kwargs)

Returns variants present in this collection but not in any of the others.

Source code in varcode/variant_collection.py

def difference(self, *others, **kwargs):
    """
    Returns variants present in this collection but not in any of the others.
    """
    return self._combine_variant_collections(
        combine_fn=set.difference,
        variant_collections=(self,) + others,
        kwargs=kwargs)

to_dataframe()

Build a DataFrame from this variant collection.

Source code in varcode/variant_collection.py

def to_dataframe(self):
    """Build a DataFrame from this variant collection."""
    structural_columns = STRUCTURAL_VARIANT_COLUMNS if any(
        getattr(v, "is_structural", False) for v in self) else ()
    def row_from_variant(variant):
        row = OrderedDict([
            ("chr", variant.contig),
            ("start", variant.original_start),
            ("ref", variant.original_ref),
            ("alt", variant.original_alt),
            ("gene_name", ";".join(variant.gene_names)),
            ("gene_id", ";".join(variant.gene_ids))
        ])
        row.update((name, getattr(variant, name, None)) for name in structural_columns)
        return row
    rows = [row_from_variant(v) for v in self]
    # Always return a DataFrame with the expected columns, even
    # when empty, so downstream code (CSV round-trip, joins) doesn't
    # have to special-case len == 0.
    return pd.DataFrame.from_records(rows, columns=self._DATAFRAME_COLUMNS + structural_columns)

to_csv(path, include_header=True)

Write this collection to CSV.

PARAMETER	DESCRIPTION
path	Output path. TYPE: str
include_header	If True (default), prepend # key=value metadata lines with varcode version and reference genome so the file can be read back via from_csv without supplying a genome explicitly. Pass include_header=False for fast consumers that don't tolerate comment lines. TYPE: bool DEFAULT: True

Source code in varcode/variant_collection.py

def to_csv(self, path, include_header=True):
    """Write this collection to CSV.

    Parameters
    ----------
    path : str
        Output path.

    include_header : bool
        If True (default), prepend ``# key=value`` metadata lines
        with varcode version and reference genome so the file can
        be read back via ``from_csv`` without supplying a genome
        explicitly. Pass ``include_header=False`` for fast
        consumers that don't tolerate comment lines.
    """
    df = self.to_dataframe()
    if not include_header:
        df.to_csv(path, index=False)
        return

    metadata = OrderedDict()
    metadata["varcode_version"] = _varcode_version
    metadata["reference_name"] = self._serialized_reference_name()
    with open(path, "w") as f:
        write_metadata_header(f, metadata)
        df.to_csv(f, index=False)

from_csv(path, genome=None, distinct=True, sort_key=variant_ascending_position_sort_key) classmethod

Rebuild a VariantCollection from a CSV previously written by VariantCollection.to_csv().

The CSV round-trip is human-readable and easy to inspect. For byte-for-byte round-trip or for faster loading of large collections (≳10k variants), prefer from_json — CSV parsing plus per-row Variant construction is significantly slower.

PARAMETER	DESCRIPTION
path	Path to the CSV file. Lines starting with '#' are treated as comments and parsed as key=value metadata. TYPE: str
genome	Reference genome to associate with the loaded variants. If None, the reference is read from the CSV's metadata header (# reference_name=...). If neither is available, raises ValueError. TYPE: pyensembl.Genome, str, int, or None DEFAULT: None
distinct	Drop duplicate variants (same as the constructor). TYPE: bool DEFAULT: True
sort_key	Sort key for the resulting collection. TYPE: callable DEFAULT: variant_ascending_position_sort_key

RETURNS	DESCRIPTION
VariantCollection

Source code in varcode/variant_collection.py

@classmethod
def from_csv(
        cls,
        path,
        genome=None,
        distinct=True,
        sort_key=variant_ascending_position_sort_key):
    """Rebuild a VariantCollection from a CSV previously written by
    ``VariantCollection.to_csv()``.

    The CSV round-trip is human-readable and easy to inspect. For
    byte-for-byte round-trip or for faster loading of large
    collections (≳10k variants), prefer ``from_json`` — CSV parsing
    plus per-row Variant construction is significantly slower.

    Parameters
    ----------
    path : str
        Path to the CSV file. Lines starting with '#' are treated as
        comments and parsed as ``key=value`` metadata.

    genome : pyensembl.Genome, str, int, or None
        Reference genome to associate with the loaded variants. If
        ``None``, the reference is read from the CSV's metadata
        header (``# reference_name=...``). If neither is available,
        raises ``ValueError``.

    distinct : bool
        Drop duplicate variants (same as the constructor).

    sort_key : callable
        Sort key for the resulting collection.

    Returns
    -------
    VariantCollection
    """
    header = read_metadata_header(path)
    warn_on_version_drift(header, _varcode_version, path)
    if genome is None:
        genome = header.get("reference_name")
    if genome is None:
        raise ValueError(
            "from_csv needs a reference genome: pass the `genome` "
            "argument explicitly, or write the CSV with "
            "`to_csv(include_header=True)` so `# reference_name=...` "
            "is recorded in the header. Neither was found at %s." % path)

    # Accept either "chr" or "contig" as the contig column so
    # CSVs are interchangeable between VariantCollection and
    # EffectCollection (openvax/varcode#274). Declaring both in
    # dtype is a no-op for whichever column is absent.
    df = pd.read_csv(
        path, comment="#", dtype={"chr": str, "contig": str})
    reject_structural_csv(df)
    contig_col = resolve_contig_column(df.columns)
    if contig_col is None:
        raise ValueError(
            "CSV at %s is missing a contig column: expected one of %s."
            % (path, list(CONTIG_COLUMN_ALIASES)))
    required = {"start", "ref", "alt"}
    missing = required - set(df.columns)
    if missing:
        raise ValueError(
            "CSV at %s is missing required columns: %s" % (
                path, sorted(missing)))

    # Extract the columns we need by name and coerce types up front,
    # then iterate with zip. This is robust against column reordering
    # and extra columns (unlike itertuples, which depends on
    # attribute access to valid-identifier column names in fixed
    # positions) and avoids the per-row overhead of iterrows.
    contigs = df[contig_col].astype(str)
    starts = df["start"].astype(int)
    refs = df["ref"].fillna("")
    alts = df["alt"].fillna("")

    variants = [
        Variant(
            contig=contig,
            start=start,
            ref=ref,
            alt=alt,
            genome=genome,
        )
        for contig, start, ref, alt in zip(contigs, starts, refs, alts)
    ]
    return cls(
        variants=variants,
        distinct=distinct,
        sort_key=sort_key,
    )

has_sample_data()

True if the collection has any per-sample genotype info.

Source code in varcode/variant_collection.py

def has_sample_data(self):
    """True if the collection has any per-sample genotype info."""
    return len(self.samples) > 0

genotype(variant, sample)

Return the Genotype for sample at variant.

PARAMETER	DESCRIPTION
variant	TYPE: Variant
sample	TYPE: str

RETURNS	DESCRIPTION
Genotype or None	None if the variant has no sample_info metadata at all (e.g. it was constructed directly rather than loaded from a multi-sample VCF).

RAISES	DESCRIPTION
SampleNotFoundError	If the variant's metadata exists but doesn't include the requested sample. Subclass of KeyError.

Source code in varcode/variant_collection.py

def genotype(self, variant, sample):
    """Return the ``Genotype`` for ``sample`` at ``variant``.

    Parameters
    ----------
    variant : Variant
    sample : str

    Returns
    -------
    Genotype or None
        ``None`` if the variant has no sample_info metadata at all
        (e.g. it was constructed directly rather than loaded from
        a multi-sample VCF).

    Raises
    ------
    SampleNotFoundError
        If the variant's metadata exists but doesn't include the
        requested sample. Subclass of ``KeyError``.
    """
    meta = self._metadata_for(variant)
    if meta is None:
        return None
    sample_info = meta.get("sample_info")
    if sample_info is None:
        return None
    if sample not in sample_info:
        raise SampleNotFoundError(
            "Sample %r not found in %s. Available samples: %s" % (
                sample, variant, sorted(sample_info.keys())))
    return Genotype.from_sample_info(sample_info[sample])

zygosity(variant, sample)

Zygosity of the given sample at the given variant.

Multi-allelic aware: at a site split into multiple Variants, each asks "does this sample carry this alt?".

Source code in varcode/variant_collection.py

def zygosity(self, variant, sample):
    """Zygosity of the given sample at the given variant.

    Multi-allelic aware: at a site split into multiple Variants,
    each asks "does this sample carry *this* alt?".
    """
    gt = self.genotype(variant, sample)
    if gt is None:
        return Zygosity.MISSING
    return gt.zygosity_for_alt(self._alt_index_for(variant))

for_sample(sample)

Return a VariantCollection restricted to variants where sample carries the alt (heterozygous or homozygous). Useful for multi-sample VCFs where not every row is called in every sample.

Source code in varcode/variant_collection.py

def for_sample(self, sample):
    """Return a VariantCollection restricted to variants where
    ``sample`` carries the alt (heterozygous or homozygous). Useful
    for multi-sample VCFs where not every row is called in every
    sample.
    """
    return self._filter_by_zygosity(
        sample,
        keep=lambda z: z in (Zygosity.HETEROZYGOUS, Zygosity.HOMOZYGOUS),
    )

heterozygous_in(sample)

Variants where sample is heterozygous for this variant's alt.

Source code in varcode/variant_collection.py

def heterozygous_in(self, sample):
    """Variants where ``sample`` is heterozygous for this variant's alt."""
    return self._filter_by_zygosity(
        sample,
        keep=lambda z: z is Zygosity.HETEROZYGOUS,
    )

homozygous_alt_in(sample)

Variants where sample is homozygous for this variant's alt.

Source code in varcode/variant_collection.py

def homozygous_alt_in(self, sample):
    """Variants where ``sample`` is homozygous for this variant's alt."""
    return self._filter_by_zygosity(
        sample,
        keep=lambda z: z is Zygosity.HOMOZYGOUS,
    )

varcode.StructuralVariant(contig: str, start: int, sv_type: str, end: Optional[int] = None, alt: Optional[str] = None, ref: str = 'N', mate_contig: Optional[str] = None, mate_start: Optional[int] = None, mate_orientation: Optional[str] = None, ci_start: Optional[Tuple[int, int]] = None, ci_end: Optional[Tuple[int, int]] = None, alt_assembly: Optional[str] = None, info: Optional[Mapping[str, Any]] = None, genome=None, ensembl=None, normalize_contig_names: bool = True, convert_ucsc_contig_names=None, affected_start: Optional[int] = None, affected_end: Optional[int] = None)

Bases: Variant

A structural variant — deletion, duplication, inversion, insertion, CNV, or breakend — too large or too complex to represent as a simple ref/alt nucleotide pair.

Subclasses :class:Variant so isinstance(v, Variant) still works; downstream code that handles variant kinds generically (effect collections, serialization) sees a :class:Variant and the shared contract still applies. The SV-specific fields (:attr:sv_type, :attr:end, breakend mate fields) live here and are consulted by SV-aware code.

The SV position model:

• :attr:start — 1-based event/record position (matches VCF POS). For parsed symbolic spans this is the retained padding base.
• :attr:end — 1-based inclusive end. For a DEL/DUP/INV/CNV this is the SV endpoint on the same contig. For an INS it equals start (insertions are zero-width in reference coords). For a BND, end == start and the other breakpoint lives in :attr:mate_contig / :attr:mate_start.
• :attr:affected_start / :attr:affected_end — inclusive bases changed by a span event. These normally equal start / end; parsed symbolic spans and paired breakend events use them to exclude retained VCF padding bases from exon and mutant-sequence annotation.

PARAMETER	DESCRIPTION
contig	Chromosome of the (first) breakpoint. TYPE: str
start	1-based start position. TYPE: int
sv_type	One of :data:SV_TYPES. TYPE: str
end	1-based inclusive end position. Defaults to start for zero-width SVs (INS, BND). TYPE: int DEFAULT: None
alt	Original ALT field from the VCF — <DEL>, <INS:ME:ALU>, G]17:198982], etc. Kept so round-trip to VCF is possible. Defaults to "<{sv_type}>". TYPE: str DEFAULT: None
ref	Original REF base (usually one nucleotide, the anchor). Defaults to "N". TYPE: str DEFAULT: 'N'
mate_contig	For BND: the mate breakpoint's chromosome. Normalized the same way as contig (e.g. "chr4" -> "4" when converting UCSC names). TYPE: str DEFAULT: None
mate_start	For BND: the mate breakpoint's position. TYPE: int DEFAULT: None
mate_orientation	For BND: one of "[[", "[]", "][, "]]" encoding the VCF 4.1 breakend strand + direction shorthand (first bracket = preceding; second = following). See VCF §5.4 for the full grammar. When alt isn't a breakend, "[[" / "]]" still tell the annotator which side of the mate is kept. TYPE: str DEFAULT: None
ci_start	Confidence interval around start (VCF CIPOS). TYPE: (int, int) DEFAULT: None
ci_end	Confidence interval around end (VCF CIEND). TYPE: (int, int) DEFAULT: None
alt_assembly	Caller-supplied assembled sequence of the rearranged allele. Hook for long-read / targeted-assembly pipelines. The SV annotator can prefer this over inferring from breakpoints. TYPE: str DEFAULT: None
info	Open-ended bag for extra VCF INFO fields the core class doesn't model (HOMLEN, SVMETHOD, MATEID, etc.). Kept as a Mapping so callers can pass whatever shape their caller produces. TYPE: Mapping[str, Any] DEFAULT: None
genome	Same meaning as :class:Variant. DEFAULT: None
ensembl	Same meaning as :class:Variant. DEFAULT: None
normalize_contig_names	Same meaning as :class:Variant. DEFAULT: None
convert_ucsc_contig_names	Same meaning as :class:Variant. DEFAULT: None
affected_start	Inclusive affected-region coordinates. Defaults to start and end respectively. Primarily used for typed breakend pairs whose event coordinates include a retained padding base. TYPE: int DEFAULT: None
affected_end	Inclusive affected-region coordinates. Defaults to start and end respectively. Primarily used for typed breakend pairs whose event coordinates include a retained padding base. TYPE: int DEFAULT: None

Source code in varcode/structural_variant.py

def __init__(
        self,
        contig: str,
        start: int,
        sv_type: str,
        end: Optional[int] = None,
        alt: Optional[str] = None,
        ref: str = "N",
        mate_contig: Optional[str] = None,
        mate_start: Optional[int] = None,
        mate_orientation: Optional[str] = None,
        ci_start: Optional[Tuple[int, int]] = None,
        ci_end: Optional[Tuple[int, int]] = None,
        alt_assembly: Optional[str] = None,
        info: Optional[Mapping[str, Any]] = None,
        genome=None,
        ensembl=None,
        normalize_contig_names: bool = True,
        convert_ucsc_contig_names=None,
        affected_start: Optional[int] = None,
        affected_end: Optional[int] = None):
    if sv_type not in SV_TYPES:
        raise ValueError(
            "Unknown sv_type %r (expected one of %s)"
            % (sv_type, sorted(SV_TYPES)))
    if end is None:
        end = start

    # Initialize the base Variant with a placeholder ref/alt so the
    # nucleotide-normalization path doesn't reject <DEL>-style
    # symbolic alleles. Restore the actual record fields below; placeholder
    # normalization must not leak into public alleles or exports (#417).
    Variant.__init__(
        self,
        contig=contig,
        start=start,
        ref=ref if ref else "N",
        alt="A",
        genome=genome,
        ensembl=ensembl,
        allow_extended_nucleotides=True,
        normalize_contig_names=normalize_contig_names,
        convert_ucsc_contig_names=convert_ucsc_contig_names)

    # Override end position — base Variant ignores end for SNVs
    # but we need it as an explicit SV endpoint.
    self.end = int(end)
    self.affected_start = int(
        start if affected_start is None else affected_start)
    self.affected_end = int(
        end if affected_end is None else affected_end)

    # SV-specific fields.
    self.sv_type = sv_type
    self.mate_contig = (
        self._normalize_contig_name(mate_contig)
        if mate_contig is not None else None)
    self.mate_start = int(mate_start) if mate_start is not None else None
    self.mate_orientation = mate_orientation
    self.ci_start = tuple(ci_start) if ci_start is not None else None
    self.ci_end = tuple(ci_end) if ci_end is not None else None
    self.alt_assembly = alt_assembly
    self.info = dict(info) if info is not None else {}

    # Preserve the original symbolic ALT string so round-tripping
    # and downstream consumers see what the VCF said.
    self._sv_alt = alt if alt is not None else "<%s>" % sv_type
    self.ref = self.original_ref
    self.alt = self.original_alt = self._sv_alt
    self.start = self.original_start = int(start)

    # Caches, like Variant's overlapping-gene / transcript caches:
    # the junction list derived from the fields above, and results
    # the SV annotator derives from the variant alone (mate-locus
    # transcripts, cryptic-exon candidates), which would otherwise be
    # recomputed for every transcript the variant is annotated on.
    self._junctions = None
    self._annotation_cache = {}

symbolic_alt: str property

The original symbolic / breakend ALT string from the VCF.

junctions: Tuple[Tuple[Breakend, Breakend], ...] property

The novel adjacencies this variant creates, each a pair of :class:Breakend ends.

One junction for a breakend record (its own position joined to its mate), and one for a deletion or tandem duplication. A symbolic inversion has two, since it joins both of its ends; an inversion built from one breakend pair has only the junction that pair observed. Empty for insertions, CNVs and single breakends, which have no second end to join.

Positions follow the VCF symbolic-allele convention that start is the base before the event, so a deletion joins start (keeping the left side) to end + 1 (keeping the right side), and a tandem duplication joins end to start + 1.

breakpoints: Tuple[Tuple[str, int], ...] property

Distinct (contig, position) breakpoints of this variant's junctions, for consumers that read sequence around them (cryptic-exon and splice-window scans). Falls back to the variant's own start when it has no junction.

length: Optional[int] property

Length of the SV span in reference coordinates, when well-defined. None for breakends (the span depends on the mate, which may be on another contig).

to_dict()

The constructor arguments, so serialization (JSON, pickle) round-trips every SV field rather than only the base locus.

Source code in varcode/structural_variant.py

def to_dict(self):
    """The constructor arguments, so serialization (JSON, pickle)
    round-trips every SV field rather than only the base locus."""
    return dict(
        contig=self.original_contig,
        start=self.original_start,
        sv_type=self.sv_type,
        end=self.end,
        affected_start=self.affected_start,
        affected_end=self.affected_end,
        alt=self._sv_alt,
        ref=self.original_ref,
        mate_contig=self.mate_contig,
        mate_start=self.mate_start,
        mate_orientation=self.mate_orientation,
        ci_start=self.ci_start,
        ci_end=self.ci_end,
        alt_assembly=self.alt_assembly,
        info=self.info,
        genome=self.original_genome,
        normalize_contig_names=self.normalize_contig_names,
        convert_ucsc_contig_names=self.convert_ucsc_contig_names)

varcode.parse_symbolic_alt(contig: str, start: int, ref: str, alt: str, info=None, genome=None, normalize_contig_names: bool = True, convert_ucsc_contig_names: Optional[bool] = None) -> Optional[StructuralVariant]

Parse a single symbolic or breakend ALT into a :class:StructuralVariant. Returns None if the ALT is not symbolic (the caller keeps handling it as a simple variant).

info is an optional mapping (e.g. a pyvcf INFO dict) that may carry END, SVTYPE, CIPOS, CIEND, MATEID, etc. The parser reads those when present but doesn't require them — the ALT shape alone is enough to distinguish symbolic from breakend from inline.

normalize_contig_names and convert_ucsc_contig_names mean the same as for :class:~varcode.Variant and apply to both contig and a breakend's mate contig. The VCF loader forwards its own settings.

Source code in varcode/sv_allele_parser.py

def parse_symbolic_alt(
        contig: str,
        start: int,
        ref: str,
        alt: str,
        info=None,
        genome=None,
        normalize_contig_names: bool = True,
        convert_ucsc_contig_names: Optional[bool] = None,
) -> Optional[StructuralVariant]:
    """Parse a single symbolic or breakend ALT into a
    :class:`StructuralVariant`. Returns ``None`` if the ALT is not
    symbolic (the caller keeps handling it as a simple variant).

    ``info`` is an optional mapping (e.g. a ``pyvcf`` INFO dict) that
    may carry ``END``, ``SVTYPE``, ``CIPOS``, ``CIEND``, ``MATEID``,
    etc. The parser reads those when present but doesn't require
    them — the ALT shape alone is enough to distinguish symbolic
    from breakend from inline.

    ``normalize_contig_names`` and ``convert_ucsc_contig_names`` mean
    the same as for :class:`~varcode.Variant` and apply to both
    ``contig`` and a breakend's mate contig. The VCF loader forwards
    its own settings.
    """
    if not alt:
        return None

    # Symbolic allele: <DEL>, <DUP>, <INS:ME:ALU>, <CN0>, etc.
    m = _SYMBOLIC_RE.match(alt)
    if m:
        token = m.group(1).upper()
        # Top-level type is the first colon-delimited token; e.g.
        # ``INS:ME:ALU`` → ``INS``. The rest stays in ``info`` under
        # the custom ``symbolic_subtype`` key.
        top, _, subtype = token.partition(":")

        # Copy-number alleles: <CN0>, <CN1>, <CN2>, ..., <CNV>. We collapse
        # to sv_type="CNV" but preserve the integer count separately so
        # downstream code can distinguish a deletion (CN0 / CN1 in a
        # diploid) from a duplication (CN3+) without re-parsing the ALT.
        copy_number = None
        cn_match = _CN_TOKEN_RE.match(top)
        if cn_match:
            cn_digits = cn_match.group("n")
            if cn_digits:
                copy_number = int(cn_digits)
            top = "CNV"

        if top not in SV_TYPES:
            # Unknown symbolic — keep the raw token as subtype, fall
            # back to BND so the caller still gets a usable object.
            subtype = token
            top = "BND"

        end = _extract_info(info, "END") or start
        svtype_from_info = _extract_info(info, "SVTYPE")
        ci_start = _extract_info(info, "CIPOS")
        ci_end = _extract_info(info, "CIEND")

        # INFO/SVTYPE overrides when present and recognized.
        if svtype_from_info and svtype_from_info.upper() in SV_TYPES:
            top = svtype_from_info.upper()

        # Inserted-sequence hint for INS rows. Long-read callers ship
        # the assembled inserted sequence as INSSEQ (Manta) or
        # SVINSSEQ (some Sniffles/PBSV variants). Both are accepted;
        # the resolved sequence flows through StructuralVariant's
        # ``alt_assembly`` slot, where the SV annotator already prefers
        # it over inferring from breakpoint coordinates alone.
        alt_assembly = None
        if top == "INS":
            alt_assembly = (
                _extract_info_scalar(info, "INSSEQ")
                or _extract_info_scalar(info, "SVINSSEQ"))

        sv_info = {}
        if subtype:
            sv_info["symbolic_subtype"] = subtype
        if copy_number is not None:
            sv_info["copy_number"] = copy_number

        # VCF symbolic POS is the retained padding base before the event.
        # Keep record/junction coordinates unchanged; span consumers already
        # use affected_start/end for paired breakends (#404).
        affected_start = None
        if top in {"DEL", "DUP", "INV", "CNV"}:
            if int(end) <= int(start):
                raise ValueError("A symbolic %s requires END greater than POS" % top)
            affected_start = int(start) + 1

        return StructuralVariant(
            contig=contig,
            start=int(start),
            end=int(end),
            affected_start=affected_start,
            sv_type=top,
            alt=alt,
            ref=ref or "N",
            alt_assembly=alt_assembly,
            ci_start=tuple(ci_start) if ci_start else None,
            ci_end=tuple(ci_end) if ci_end else None,
            info=sv_info,
            genome=genome,
            normalize_contig_names=normalize_contig_names,
            convert_ucsc_contig_names=convert_ucsc_contig_names,
        )

    # Breakend: t[CHROM:POS[ etc. Four orientation shapes per §5.4.
    m = _BREAKEND_RE.match(alt)
    if m:
        prefix = m.group("prefix")
        suffix = m.group("suffix")
        open_br = m.group("open")
        close_br = m.group("close")
        mate_contig = m.group("mate_contig")
        mate_pos = int(m.group("mate_pos"))
        # Two-letter orientation token captures which side the REF
        # base is on (prefix = joined-after, suffix = joined-before)
        # and the bracket direction (mate strand).
        orientation = open_br + close_br
        # One record is one breakend, even when the caller labels it
        # SVTYPE=DEL / DUP / INV: that label describes the event both
        # halves make together, which
        # :func:`~varcode.transforms.pair_breakends` builds from the
        # pair. The label rides along in ``info`` for it to read.
        return StructuralVariant(
            contig=contig,
            start=int(start),
            end=int(start),
            sv_type="BND",
            alt=alt,
            ref=ref or "N",
            mate_contig=mate_contig,
            mate_start=mate_pos,
            mate_orientation=orientation,
            info={"mateid": (_extract_info(info, "MATEID")
                             or _extract_info(info, "PARID")),
                  "bnd_anchor": prefix or suffix,
                  "svtype": _extract_info_scalar(info, "SVTYPE")},
            genome=genome,
            normalize_contig_names=normalize_contig_names,
            convert_ucsc_contig_names=convert_ucsc_contig_names,
        )

    # Single breakend: ``.ACGT`` / ``ACGT.``. There's no mate to point
    # at, so it's a BND with ``mate_contig=None``, which the annotator
    # treats like a breakend into intergenic space.
    m = _SINGLE_BREAKEND_RE.match(alt)
    if m:
        return StructuralVariant(
            contig=contig,
            start=int(start),
            end=int(start),
            sv_type="BND",
            alt=alt,
            ref=ref or "N",
            info={"single_breakend": True,
                  "bnd_anchor": (m.group("joined_before")
                                 or m.group("joined_after")),
                  "svtype": _extract_info_scalar(info, "SVTYPE")},
            genome=genome,
            normalize_contig_names=normalize_contig_names,
            convert_ucsc_contig_names=convert_ucsc_contig_names,
        )

    # Spanning-deletion placeholder. VCF 4.2+ uses ``*`` to mean "this
    # position is covered by a deletion recorded on another row".
    # We drop these like before — there's no SV object to construct.
    return None

varcode.SV_TYPES = frozenset({'DEL', 'DUP', 'INV', 'INS', 'CNV', 'BND'}) module-attribute

Genotypes

varcode.Genotype(raw_gt: str, alleles: Tuple[Optional[int], ...], phased: bool = False, phase_set: Optional[int] = None, allele_depths: Optional[Tuple[int, ...]] = None, total_depth: Optional[int] = None, genotype_quality: Optional[int] = None) dataclass

Bases: DataclassSerializable

One sample's genotype at one variant locus.

The alleles tuple encodes the observed alleles using VCF GT semantics: 0 is the reference allele, 1 is the first ALT listed on the VCF row, 2 is the second, and so on. None indicates a no-call on that haplotype.

For varcode's variant-level API, note that Variant.alt is a specific alt (a multi-allelic VCF row is split into one Variant per alt). When querying zygosity relative to a Variant, use the variant's alt_allele_index from the collection's metadata and add 1 to get the GT-encoded index, then call :meth:zygosity_for_alt or :meth:carries_alt.

is_called: bool property

True if at least one allele is non-None.

ploidy: int property

Number of alleles in the call (including missing).

from_sample_info(sample_info) classmethod

Build a Genotype from pyvcf's call.data._asdict() output.

Handles the keys varcode normally sees: GT, AD, DP, GQ, PS. Missing keys default to None.

Source code in varcode/genotype.py

@classmethod
def from_sample_info(cls, sample_info):
    """Build a Genotype from pyvcf's ``call.data._asdict()`` output.

    Handles the keys varcode normally sees: ``GT``, ``AD``, ``DP``,
    ``GQ``, ``PS``. Missing keys default to ``None``.
    """
    if sample_info is None:
        return cls(raw_gt="./.", alleles=(None, None), phased=False)
    gt_str = sample_info.get("GT")
    if gt_str is None:
        gt_str = "./."
    alleles, phased = parse_gt_string(gt_str)
    ad = sample_info.get("AD")
    return cls(
        raw_gt=gt_str,
        alleles=alleles,
        phased=phased,
        phase_set=sample_info.get("PS"),
        allele_depths=tuple(ad) if ad is not None else None,
        total_depth=sample_info.get("DP"),
        genotype_quality=sample_info.get("GQ"),
    )

carries_alt(alt_index: int) -> bool

True if this sample's genotype contains the given alt.

alt_index uses VCF GT encoding: 1 is the first alt on the row, 2 is the second, etc. (i.e. one more than alt_allele_index from the VariantCollection metadata).

Source code in varcode/genotype.py

def carries_alt(self, alt_index: int) -> bool:
    """True if this sample's genotype contains the given alt.

    ``alt_index`` uses VCF GT encoding: ``1`` is the first alt on
    the row, ``2`` is the second, etc. (i.e. one more than
    ``alt_allele_index`` from the VariantCollection metadata).
    """
    return any(
        a == alt_index
        for a in self.alleles
        if a is not None
    )

copies_of_alt(alt_index: int) -> int

Number of haplotypes carrying the given alt.

Source code in varcode/genotype.py

def copies_of_alt(self, alt_index: int) -> int:
    """Number of haplotypes carrying the given alt."""
    return sum(
        1 for a in self.alleles
        if a is not None and a == alt_index
    )

zygosity_for_alt(alt_index: int) -> Zygosity

Classify the sample's zygosity relative to one alt allele.

Multi-allelic aware: GT=1/2 queried for alt 1 returns HETEROZYGOUS (one copy of this alt, one of a different alt); queried for alt 3 it returns ABSENT.

Source code in varcode/genotype.py

def zygosity_for_alt(self, alt_index: int) -> Zygosity:
    """Classify the sample's zygosity relative to one alt allele.

    Multi-allelic aware: ``GT=1/2`` queried for alt ``1`` returns
    ``HETEROZYGOUS`` (one copy of this alt, one of a different
    alt); queried for alt ``3`` it returns ``ABSENT``.
    """
    called = [a for a in self.alleles if a is not None]
    if len(called) == 0:
        return Zygosity.MISSING
    n_copies = sum(1 for a in called if a == alt_index)
    if n_copies == 0:
        return Zygosity.ABSENT
    if n_copies == len(called):
        return Zygosity.HOMOZYGOUS
    return Zygosity.HETEROZYGOUS

depth_for_alt(alt_index: int) -> Optional[int]

Per-allele read depth for a given alt, from the AD field.

AD is indexed with ref at position 0 and alt #1 at position 1, etc., so alt_index should use GT encoding (1 = first alt).

Source code in varcode/genotype.py

def depth_for_alt(self, alt_index: int) -> Optional[int]:
    """Per-allele read depth for a given alt, from the ``AD`` field.

    ``AD`` is indexed with ref at position 0 and alt #1 at
    position 1, etc., so ``alt_index`` should use GT encoding
    (``1`` = first alt).
    """
    if self.allele_depths is None:
        return None
    if alt_index >= len(self.allele_depths):
        return None
    return self.allele_depths[alt_index]

varcode.Zygosity

Bases: Enum

Zygosity of a sample's genotype relative to a specific alt allele.

ABSENT is distinct from MISSING: ABSENT means the call exists but doesn't include the alt in question (e.g. the sample is ref-ref, or carries a different alt at a multi-allelic site). MISSING means the call itself is ./. or the sample wasn't called.

VariantCollection transforms

varcode.transforms.pair_breakends(vc)

Merge MATEID-paired BND rows into a single combined StructuralVariant per rearrangement event. (reduces)

For each pair of :class:~varcode.StructuralVariant rows where row A's MATEID references row B's VCF ID (and vice versa), emit one combined row carrying both endpoints. The combined variant's source_variants attribute holds the two originals.

Non-BND variants, single-row TRA, and single-ended BNDs (no MATEID) pass through unchanged with source_variants=().

Pairing rules:

• Primary key: MATEID field on each variant's info dict against the VCF row ID stored in the collection's source metadata.
• Alias: PARID (used by older GRIDSS) is treated as MATEID.
• Symmetric: row A's MATEID must equal row B's ID and row B's MATEID must equal row A's ID. Asymmetric references are warned and left unpaired.
• If a MATEID points to an ID not present in this collection (filtered out, chunked load), a warning is emitted and the variant passes through unpaired.
• If three or more rows share a MATEID group, the whole group is left unpaired with a warning (pairing is ambiguous).
• Already-paired input (source_variants non-empty) passes through unchanged — :func:pair_breakends is idempotent.

Metadata merge for paired rows:

• Genotype: both halves must agree on the per-sample GT. The combined variant inherits the shared genotype. Disagreement raises :class:ValueError.
• alt_assembly: if exactly one half carries an assembled sequence, the combined variant inherits it. If both differ, A's wins (deterministic via lex source-ID ordering) with a warning.
• filter: union of both halves' FILTER tokens; PASS is dropped if any non-PASS label is present (stricter wins).
• qual: minimum of the two halves' quality scores.
• Other INFO fields: prefer A's; the other half is reachable via combined.source_variants.

PARAMETER	DESCRIPTION
vc	Input collection. May contain a mix of structural and non-structural variants. TYPE: VariantCollection

RETURNS	DESCRIPTION
VariantCollection	A new collection. SV rows that were halves of paired BNDs are replaced by combined rows; everything else (including SNVs/indels/MNVs and unpaired SVs) passes through.

Source code in varcode/transforms/__init__.py

def pair_breakends(vc):
    """Merge MATEID-paired BND rows into a single combined
    ``StructuralVariant`` per rearrangement event. (reduces)

    For each pair of :class:`~varcode.StructuralVariant` rows where
    row A's ``MATEID`` references row B's VCF ID (and vice versa),
    emit one combined row carrying both endpoints. The combined
    variant's ``source_variants`` attribute holds the two originals.

    Non-BND variants, single-row TRA, and single-ended BNDs (no
    ``MATEID``) pass through unchanged with ``source_variants=()``.

    Pairing rules:

    * Primary key: ``MATEID`` field on each variant's ``info`` dict
      against the VCF row ID stored in the collection's source
      metadata.
    * Alias: ``PARID`` (used by older GRIDSS) is treated as ``MATEID``.
    * Symmetric: row A's ``MATEID`` must equal row B's ID and row B's
      ``MATEID`` must equal row A's ID. Asymmetric references are
      warned and left unpaired.
    * If a ``MATEID`` points to an ID not present in this collection
      (filtered out, chunked load), a warning is emitted and the
      variant passes through unpaired.
    * If three or more rows share a ``MATEID`` group, the whole group
      is left unpaired with a warning (pairing is ambiguous).
    * Already-paired input (``source_variants`` non-empty) passes
      through unchanged — :func:`pair_breakends` is idempotent.

    Metadata merge for paired rows:

    * Genotype: both halves must agree on the per-sample ``GT``. The
      combined variant inherits the shared genotype. Disagreement
      raises :class:`ValueError`.
    * ``alt_assembly``: if exactly one half carries an assembled
      sequence, the combined variant inherits it. If both differ,
      A's wins (deterministic via lex source-ID ordering) with a
      warning.
    * ``filter``: union of both halves' FILTER tokens; ``PASS`` is
      dropped if any non-PASS label is present (stricter wins).
    * ``qual``: minimum of the two halves' quality scores.
    * Other INFO fields: prefer A's; the other half is reachable
      via ``combined.source_variants``.

    Parameters
    ----------
    vc : VariantCollection
        Input collection. May contain a mix of structural and
        non-structural variants.

    Returns
    -------
    VariantCollection
        A new collection. SV rows that were halves of paired BNDs
        are replaced by combined rows; everything else (including
        SNVs/indels/MNVs and unpaired SVs) passes through.
    """
    variant_to_id, id_to_variants = _build_id_indices(vc)

    # For each BND variant, what does it claim as its mate?
    bnd_to_mateid = {}
    for variant in vc:
        if not _is_breakend(variant):
            continue
        if getattr(variant, "source_variants", ()):
            continue
        mate_id = _mate_reference(variant)
        if mate_id is None:
            continue
        bnd_to_mateid[variant] = mate_id

    # How many rows reference each ID as their mate? An ID with
    # in-degree > 1 indicates ambiguity — a clean pair requires the
    # MATEID pointer to be 1:1.
    mate_in_degree = {}
    for mate_id in bnd_to_mateid.values():
        mate_in_degree[mate_id] = mate_in_degree.get(mate_id, 0) + 1

    # Walk BND variants and resolve each to one of: paired (with a
    # specific mate variant), ambiguous (skip with warning naming the
    # connected component), missing-mate (warn), asymmetric (warn).
    replacement = {}     # original variant -> combined variant
    processed = set()    # variants whose pairing decision is final
    ambiguity_warned = set()   # connected-component ids already warned about

    def _component_ids(seed_id):
        """Connected component of mate references reachable from
        ``seed_id``. Used to produce one warning per ambiguous group
        rather than one per row."""
        visited = set()
        stack = [seed_id]
        while stack:
            node = stack.pop()
            if node in visited:
                continue
            visited.add(node)
            # Outgoing: variants whose ID is `node` -- find their mateids.
            for v in id_to_variants.get(node, ()):
                mid = bnd_to_mateid.get(v)
                if mid is not None and mid not in visited:
                    stack.append(mid)
            # Incoming: any mateid pointing AT `node`.
            for v, mid in bnd_to_mateid.items():
                if mid == node:
                    own = variant_to_id.get(v)
                    if own is not None and own not in visited:
                        stack.append(own)
        return visited

    for variant in vc:
        if not _is_breakend(variant):
            continue
        if variant in processed:
            continue
        if getattr(variant, "source_variants", ()):
            continue
        own_id = variant_to_id.get(variant)
        if own_id is None:
            continue
        mate_id = bnd_to_mateid.get(variant)
        if mate_id is None:
            continue
        if mate_id not in id_to_variants:
            warnings.warn(
                "pair_breakends: BND %r references MATEID/PARID %r "
                "which is not in this collection; left unpaired. "
                "This typically means the mate row was filter-dropped "
                "or split across chunked loads."
                % (own_id, mate_id),
                stacklevel=2)
            processed.add(variant)
            continue
        # Ambiguity: this variant or its claimed mate has incoming
        # mate-degree > 1, meaning at least one other row also points
        # at the same target. Whole connected component is unsafe.
        if (mate_in_degree.get(own_id, 0) > 1
                or mate_in_degree.get(mate_id, 0) > 1):
            component = _component_ids(own_id)
            component_key = frozenset(component)
            if component_key not in ambiguity_warned:
                ambiguity_warned.add(component_key)
                warnings.warn(
                    "pair_breakends: BND ID group %r has ambiguous "
                    "MATEID/PARID references (at least one ID is "
                    "referenced by more than one row); left unpaired. "
                    "Each rearrangement should produce exactly two BND "
                    "rows with 1:1 mate pointers."
                    % sorted(component),
                    stacklevel=2)
            # Mark every variant in the component as processed so we
            # don't re-warn from a different vantage point.
            for cid in component:
                for v in id_to_variants.get(cid, ()):
                    processed.add(v)
            continue
        # Identify the specific mate variant. With in-degree 1 there's
        # exactly one BND row at `mate_id`.
        mate_candidates = [
            v for v in id_to_variants[mate_id]
            if _is_breakend(v)
        ]
        if not mate_candidates:
            warnings.warn(
                "pair_breakends: BND %r references MATEID %r but that "
                "ID is held by a non-BND variant; left unpaired."
                % (own_id, mate_id),
                stacklevel=2)
            processed.add(variant)
            continue
        mate = mate_candidates[0]
        # Symmetric check: mate must point back at us.
        if bnd_to_mateid.get(mate) != own_id:
            warnings.warn(
                "pair_breakends: BND %r mate reference is asymmetric "
                "(this row -> %r but mate -> %r); left unpaired."
                % (own_id, mate_id, bnd_to_mateid.get(mate)),
                stacklevel=2)
            processed.add(variant)
            processed.add(mate)
            continue
        # Clean pair.
        a, b = sorted([variant, mate], key=lambda v: variant_to_id[v])
        a_id = variant_to_id[a]
        b_id = variant_to_id[b]
        combined = _build_combined(a, b, a_id, b_id)
        replacement[a] = combined
        replacement[b] = combined
        processed.add(a)
        processed.add(b)

    # Build output VC. Pass-through variants land first-occurrence; paired
    # variants are replaced by their combined variant the first time either
    # half is encountered, then the second half is skipped to preserve
    # ordering and cardinality.
    out_variants = []
    emitted_combined = set()
    for variant in vc:
        combined = replacement.get(variant)
        if combined is None:
            out_variants.append(variant)
            continue
        if id(combined) in emitted_combined:
            continue
        emitted_combined.add(id(combined))
        out_variants.append(combined)

    # Build output metadata. Pass-through variants reuse the existing
    # metadata entries by reference; combined variants get a freshly
    # merged entry written to every source path that contained either
    # half.
    out_metadata = {path: {} for path in vc.source_to_metadata_dict}
    for path, by_variant in vc.source_to_metadata_dict.items():
        out_by_variant = out_metadata[path]
        for variant, meta in by_variant.items():
            combined = replacement.get(variant)
            if combined is None:
                out_by_variant[variant] = meta
            elif combined not in out_by_variant:
                a, b = combined.source_variants
                a_meta = by_variant.get(a)
                b_meta = by_variant.get(b)
                a_id = variant_to_id[a]
                b_id = variant_to_id[b]
                out_by_variant[combined] = _merge_metadata(
                    a_meta, b_meta, a_id, b_id)

    return VariantCollection(
        variants=out_variants,
        sources=vc.sources,
        source_to_metadata_dict=out_metadata,
    )

varcode.transforms.left_align_indels(vc)

Shift indels to their canonical leftmost equivalent position. (preserves)

Indels in homopolymer or short-tandem-repeat regions can be represented at any of several equivalent positions — CTT->T inside a CT-repeat means the same biological event as CT->_ two positions to the left. Tools that compare variants by (contig, start, ref, alt) see those representations as distinct calls. Left-alignment normalizes to a single canonical representation per indel: the leftmost equivalent position.

The algorithm is the standard variant-normalization left-shift used by bcftools norm and GATK LeftAlignAndTrimVariants, applied as an opt-in VariantCollection -> VariantCollection transform rather than baked into VCF load.

Reference sequence is read via the genome the variants carry — no explicit reference parameter. Coverage tiers (see :mod:varcode.genome_sequence):

• Chromosome FASTA attached (via :class:varcode.Genome's fasta slot): indels everywhere shift to canonical positions, including in introns and intergenic regions.
• No FASTA (default pyensembl install): indels fully within an exon shift via the transcript cDNA fallback. Intronic and intergenic indels pass through unchanged. Indels that start exonic but would shift across an exon boundary stop at the boundary and carry info["left_align_partial"] = True.

PARAMETER	DESCRIPTION
vc	Input collection. May contain a mix of SNVs, MNVs, indels, complex variants, and SVs — only pure indels (length-different REF/ALT with one side empty after suffix trimming) are considered for shifting. Everything else passes through. TYPE: VariantCollection

RETURNS	DESCRIPTION
VariantCollection	A new collection with indels at their canonical leftmost positions. Variants that shifted carry source_variants=(original,); everything else passes through with source_variants=().
See	TYPE: doc:`/transforms` for the behavior table covering all
six (location × FASTA-attached) combinations and the metadata
fields the transform writes.

Examples:

>>> from varcode.transforms import left_align_indels
>>> normalized = left_align_indels(vc)

Source code in varcode/transforms/__init__.py

def left_align_indels(vc):
    """Shift indels to their canonical leftmost equivalent position. (preserves)

    Indels in homopolymer or short-tandem-repeat regions can be
    represented at any of several equivalent positions — ``CTT->T``
    inside a ``CT``-repeat means the same biological event as
    ``CT->_`` two positions to the left. Tools that compare variants
    by ``(contig, start, ref, alt)`` see those representations as
    distinct calls. Left-alignment normalizes to a single canonical
    representation per indel: the leftmost equivalent position.

    The algorithm is the standard variant-normalization left-shift
    used by ``bcftools norm`` and GATK ``LeftAlignAndTrimVariants``,
    applied as an opt-in ``VariantCollection -> VariantCollection``
    transform rather than baked into VCF load.

    Reference sequence is read via the genome the variants carry —
    no explicit ``reference`` parameter. Coverage tiers (see
    :mod:`varcode.genome_sequence`):

    * **Chromosome FASTA attached** (via :class:`varcode.Genome`'s
      ``fasta`` slot): indels everywhere shift to canonical
      positions, including in introns and intergenic regions.
    * **No FASTA** (default pyensembl install): indels fully within
      an exon shift via the transcript cDNA fallback. Intronic and
      intergenic indels pass through unchanged. Indels that *start*
      exonic but would shift across an exon boundary stop at the
      boundary and carry ``info["left_align_partial"] = True``.

    Parameters
    ----------
    vc : VariantCollection
        Input collection. May contain a mix of SNVs, MNVs, indels,
        complex variants, and SVs — only pure indels (length-different
        REF/ALT with one side empty after suffix trimming) are
        considered for shifting. Everything else passes through.

    Returns
    -------
    VariantCollection
        A new collection with indels at their canonical leftmost
        positions. Variants that shifted carry
        ``source_variants=(original,)``; everything else passes through
        with ``source_variants=()``.

    See :doc:`/transforms` for the behavior table covering all
    six (location × FASTA-attached) combinations and the metadata
    fields the transform writes.

    Examples
    --------

    >>> from varcode.transforms import left_align_indels
    >>> normalized = left_align_indels(vc)  # doctest: +SKIP
    """
    # original variant -> (shifted variant, bounded_by_coverage flag).
    # The flag rides with the shifted variant so we don't need a
    # second data structure keyed on id() — the lifecycle of the
    # flag and the lifecycle of the shifted variant are identical.
    replacement = {}

    for variant in vc:
        if not variant.is_indel:
            continue
        shifted, partial = _left_align_one(variant, _reference_range)
        if shifted is variant:
            continue
        replacement[variant] = (shifted, partial)

    if not replacement:
        return vc

    out_variants = [
        replacement[v][0] if v in replacement else v
        for v in vc
    ]

    out_metadata = {}
    for path, by_variant in vc.source_to_metadata_dict.items():
        out_by_variant = {}
        for variant, meta in by_variant.items():
            entry = replacement.get(variant)
            if entry is None:
                out_by_variant[variant] = meta
            else:
                shifted, partial = entry
                new_meta = dict(meta) if meta else {}
                info = dict(new_meta.get("info") or {})
                info["original_start"] = variant.start
                if partial:
                    info["left_align_partial"] = True
                new_meta["info"] = info
                out_by_variant[shifted] = new_meta
        out_metadata[path] = out_by_variant

    return VariantCollection(
        variants=out_variants,
        sources=vc.sources,
        source_to_metadata_dict=out_metadata,
    )

File loading

varcode.vcf.load_vcf(path, genome=None, reference_vcf_key='reference', only_passing=True, allow_extended_nucleotides=False, include_info=True, chunk_size=10 ** 5, max_variants=None, sort_key=variant_ascending_position_sort_key, distinct=True, normalize_contig_names=True, convert_ucsc_contig_names=True, parse_structural_variants=False, genome_fasta=None)

Load reference name and Variant objects from the given VCF filename.

Local files are parsed directly. HTTP/HTTPS URLs are downloaded to a temporary file and load_vcf recurses on the local copy; pandas doesn't reliably stream gzipped HTTP responses, so we materialize first.

PARAMETER	DESCRIPTION
path	Path to VCF (.vcf) or compressed VCF (.vcf.gz). TYPE: str
genome	Optionally pass in a PyEnsembl Genome object, name of reference, or PyEnsembl release version to specify the reference associated with a VCF (otherwise infer reference from VCF using reference_vcf_key) TYPE: pyensembl.Genome, reference name, Ensembl version int DEFAULT: pyensembl.Genome
reference_vcf_key	Name of metadata field which contains path to reference FASTA file (default = 'reference') TYPE: str DEFAULT: 'reference'
only_passing	If true, any entries whose FILTER field is not one of "." or "PASS" is dropped. TYPE: bool DEFAULT: True
allow_extended_nucleotides	Allow characters other that A,C,T,G in the ref and alt strings. TYPE: bool DEFAULT: False
include_info	Whether to parse the INFO and per-sample columns. If you don't need these, set to False for faster parsing. TYPE: bool DEFAULT: True
chunk_size	Number of records to load in memory at once. DEFAULT: 10 ** 5
max_variants	If specified, return only the first max_variants variants. TYPE: int DEFAULT: None
sort_key	Function which maps each element to a sorting criterion. Set to None to not to sort the variants. TYPE: fn DEFAULT: variant_ascending_position_sort_key
distinct	Don't keep repeated variants TYPE: bool DEFAULT: True
normalize_contig_names	By default contig names will be normalized by converting integers to strings (e.g. 1 -> "1"), and converting any letters after "chr" to uppercase (e.g. "chrx" -> "chrX"). If you don't want this behavior then pass normalize_contig_names=False. TYPE: bool DEFAULT: True
convert_ucsc_contig_names	Convert chromosome names from hg19 (e.g. "chr1") to equivalent names for GRCh37 (e.g. "1"). By default this is set to True. If None, it also evaluates to True if the genome of the VCF is a UCSC reference. TYPE: bool DEFAULT: True
genome_fasta	Optionally attach a chromosome FASTA to the resolved genome before parsing. Equivalent to wrapping with varcode.Genome(genome, fasta=genome_fasta) and passing the wrapper; see :class:varcode.Genome for the accepted types and verification behavior. Without this, features that need raw genomic bases (cryptic-exon scoring, indel left-alignment, sequence-aware splice prediction) fall back to transcript-cDNA coverage only. TYPE: str, path-like, or FASTA object DEFAULT: None

Source code in varcode/vcf.py

def load_vcf(
        path,
        genome=None,
        reference_vcf_key="reference",
        only_passing=True,
        allow_extended_nucleotides=False,
        include_info=True,
        chunk_size=10 ** 5,
        max_variants=None,
        sort_key=variant_ascending_position_sort_key,
        distinct=True,
        normalize_contig_names=True,
        convert_ucsc_contig_names=True,
        parse_structural_variants=False,
        genome_fasta=None):
    """
    Load reference name and Variant objects from the given VCF filename.

    Local files are parsed directly. HTTP/HTTPS URLs are downloaded to a
    temporary file and ``load_vcf`` recurses on the local copy; pandas
    doesn't reliably stream gzipped HTTP responses, so we materialize first.

    Parameters
    ----------

    path : str
        Path to VCF (*.vcf) or compressed VCF (*.vcf.gz).

    genome : {pyensembl.Genome, reference name, Ensembl version int}, optional
        Optionally pass in a PyEnsembl Genome object, name of reference, or
        PyEnsembl release version to specify the reference associated with a
        VCF (otherwise infer reference from VCF using reference_vcf_key)

    reference_vcf_key : str, optional
        Name of metadata field which contains path to reference FASTA
        file (default = 'reference')

    only_passing : bool, optional
        If true, any entries whose FILTER field is not one of "." or "PASS" is
        dropped.

    allow_extended_nucleotides : bool, default False
        Allow characters other that A,C,T,G in the ref and alt strings.

    include_info : bool, default True
        Whether to parse the INFO and per-sample columns. If you don't need
        these, set to False for faster parsing.

    chunk_size: int, optional
        Number of records to load in memory at once.

    max_variants : int, optional
        If specified, return only the first max_variants variants.

    sort_key : fn
        Function which maps each element to a sorting criterion.
        Set to None to not to sort the variants.

    distinct : bool, default True
        Don't keep repeated variants

    normalize_contig_names : bool, default True
        By default contig names will be normalized by converting integers
        to strings (e.g. 1 -> "1"), and converting any letters after "chr"
        to uppercase (e.g. "chrx" -> "chrX"). If you don't want
        this behavior then pass normalize_contig_names=False.

    convert_ucsc_contig_names : bool, default True
        Convert chromosome names from hg19 (e.g. "chr1") to equivalent names
        for GRCh37 (e.g. "1"). By default this is set to True. If None, it
        also evaluates to True if the genome of the VCF is a UCSC reference.

    genome_fasta : str, path-like, or FASTA object, optional
        Optionally attach a chromosome FASTA to the resolved genome
        before parsing. Equivalent to wrapping with
        ``varcode.Genome(genome, fasta=genome_fasta)`` and passing the
        wrapper; see :class:`varcode.Genome` for the accepted types
        and verification behavior. Without this, features that need
        raw genomic bases (cryptic-exon scoring, indel
        left-alignment, sequence-aware splice prediction) fall back
        to transcript-cDNA coverage only.
    """

    require_string(path, "Path or URL to VCF")
    parsed_path = parse_url_or_path(path)

    if parsed_path.scheme and parsed_path.scheme.lower() != "file":
        # pandas.read_table nominally supports HTTP, but it tends to crash on
        # large files and does not support gzip. Switching to the python-based
        # implementation of read_table (with engine="python") helps with some
        # issues but introduces a new set of problems (e.g. the dtype parameter
        # is not accepted). For these reasons, we're currently not attempting
        # to load VCFs over HTTP with pandas directly, and instead download it
        # to a temporary file and open that.

        (filename, headers) = urllib.request.urlretrieve(path)
        try:
            # The downloaded file has no file extension, which confuses pyvcf
            # for gziped files in Python 3. We rename it to have the correct
            # file extension.
            new_filename = "%s.%s" % (
                filename, parsed_path.path.split(".")[-1])
            os.rename(filename, new_filename)
            filename = new_filename
            return load_vcf(
                filename,
                genome=genome,
                reference_vcf_key=reference_vcf_key,
                only_passing=only_passing,
                allow_extended_nucleotides=allow_extended_nucleotides,
                include_info=include_info,
                genome_fasta=genome_fasta,
                chunk_size=chunk_size,
                max_variants=max_variants,
                sort_key=sort_key,
                distinct=distinct,
                normalize_contig_names=normalize_contig_names,
                convert_ucsc_contig_names=convert_ucsc_contig_names,
                parse_structural_variants=parse_structural_variants)
        finally:
            logger.info("Removing temporary file: %s", filename)
            os.unlink(filename)

    # Loading a local file.
    # The file will be opened twice: first to parse the header, then by
    # pandas to read the data. Normalize away any file:// scheme so the
    # opener (and pandas, below) get a plain filesystem path.
    if parsed_path.scheme and parsed_path.scheme.lower() == "file":
        path = parsed_path.path
    header = VCFHeader.from_path(path)

    ####
    # The following code looks a bit crazy because it's motivated by the
    # desired to preserve UCSC reference names even though the Variant
    # objects we're creating will convert them to EnsemblRelease genomes
    # with different reference names.
    #
    # For example, if a VCF is aligned against 'hg19' then we want to create a
    # variant which has 'hg19' as its genome argument, so that serialization
    # back to VCF will put the correct reference genome in the generated
    # header.
    if genome is None:
        if reference_vcf_key not in header.metadata:
            raise ValueError("Unable to infer reference genome for %s" % (path,))
        genome = header.metadata[reference_vcf_key]

    genome, genome_was_ucsc = infer_genome(genome)
    if genome_was_ucsc:
        genome = ensembl_to_ucsc_reference_names[genome.reference_name]

    if convert_ucsc_contig_names is None:
        convert_ucsc_contig_names = genome_was_ucsc

    if genome_fasta is not None:
        # Wrap the resolved genome in varcode.Genome so the FASTA is
        # part of the genome object rather than a side-channel
        # attachment. Late import keeps the cost on callers who opt in.
        from .genome import Genome as VarcodeGenome
        genome = VarcodeGenome(genome, fasta=genome_fasta)

    df_iterator = read_vcf_into_dataframe(
        path,
        include_info=include_info,
        sample_names=header.samples if include_info else None,
        chunk_size=chunk_size)

    if include_info:
        sample_info_parser = header.parse_samples
    else:
        sample_info_parser = None

    variant_kwargs = {
        'genome': genome,
        'allow_extended_nucleotides': allow_extended_nucleotides,
        'normalize_contig_names': normalize_contig_names,
        'convert_ucsc_contig_names': convert_ucsc_contig_names,
    }

    variant_collection_kwargs = {
        'sort_key': sort_key,
        'distinct': distinct
    }

    # TODO: drop chrMT variants from hg19 and warn user about it

    return dataframes_to_variant_collection(
        df_iterator,
        source_path=path,
        info_parser=header.parse_info if include_info else None,
        only_passing=only_passing,
        max_variants=max_variants,
        sample_names=header.samples if include_info else None,
        sample_info_parser=sample_info_parser,
        variant_kwargs=variant_kwargs,
        variant_collection_kwargs=variant_collection_kwargs,
        parse_structural_variants=parse_structural_variants)

varcode.load_maf(path, optional_cols=[], sort_key=variant_ascending_position_sort_key, distinct=True, raise_on_error=True, encoding=None, nrows=None)

Load reference name and Variant objects from MAF filename.

PARAMETER	DESCRIPTION
path	Path to MAF (*.maf). TYPE: str
optional_cols	A list of MAF columns to include as metadata if they are present in the MAF. Does not result in an error if those columns are not present. TYPE: list DEFAULT: []
sort_key	Function which maps each element to a sorting criterion. Set to None to not to sort the variants. TYPE: fn DEFAULT: variant_ascending_position_sort_key
distinct	Don't keep repeated variants TYPE: bool DEFAULT: True
raise_on_error	Raise an exception upon encountering an error or just log a warning. TYPE: bool DEFAULT: True
encoding	Encoding to use for UTF when reading MAF file. TYPE: str DEFAULT: None
nrows	Limit to number of rows loaded TYPE: int DEFAULT: None

Source code in varcode/maf.py

def load_maf(
        path,
        optional_cols=[],
        sort_key=variant_ascending_position_sort_key,
        distinct=True,
        raise_on_error=True,
        encoding=None,
        nrows=None):
    """
    Load reference name and Variant objects from MAF filename.

    Parameters
    ----------

    path : str
        Path to MAF (*.maf).

    optional_cols : list, optional
        A list of MAF columns to include as metadata if they are present in the MAF.
        Does not result in an error if those columns are not present.

    sort_key : fn
        Function which maps each element to a sorting criterion.
        Set to None to not to sort the variants.

    distinct : bool
        Don't keep repeated variants

    raise_on_error : bool
        Raise an exception upon encountering an error or just log a warning.

    encoding : str, optional
        Encoding to use for UTF when reading MAF file.

    nrows : int, optional
        Limit to number of rows loaded
    """
    # pylint: disable=no-member
    # pylint gets confused by read_csv inside load_maf_dataframe
    maf_df = load_maf_dataframe(
        path,
        nrows=nrows,
        raise_on_error=raise_on_error,
        encoding=encoding)

    if len(maf_df) == 0 and raise_on_error:
        raise ValueError("Empty MAF file %s" % path)

    ensembl_objects = {}
    variants = []
    metadata = {}
    for _, x in maf_df.iterrows():
        contig = x.Chromosome
        if isnull(contig):
            error_message = "Invalid contig name: %s" % (contig,)
            if raise_on_error:
                raise ValueError(error_message)
            else:
                logging.warn(error_message)
                continue

        start_pos = x.Start_Position
        ref = x.Reference_Allele

        # it's possible in a MAF file to have multiple Ensembl releases
        # mixed in a single MAF file (the genome assembly is
        # specified by the NCBI_Build column)
        ncbi_build = x.NCBI_Build
        if ncbi_build in ensembl_objects:
            genome = ensembl_objects[ncbi_build]
        else:
            if isinstance(ncbi_build, int):
                reference_name = "B%d" % ncbi_build
            else:
                reference_name = str(ncbi_build)
            genome, _ = infer_genome(reference_name)
            ensembl_objects[ncbi_build] = genome

        # have to try both Tumor_Seq_Allele1 and Tumor_Seq_Allele2
        # to figure out which is different from the reference allele
        if x.Tumor_Seq_Allele1 != ref:
            alt = x.Tumor_Seq_Allele1
        else:
            if x.Tumor_Seq_Allele2 == ref:
                error_message = (
                    "Both tumor alleles agree with reference %s: %s" % (
                        ref, x,))
                if raise_on_error:
                    raise ValueError(error_message)
                else:
                    logging.warn(error_message)
                    continue
            alt = x.Tumor_Seq_Allele2

        variant = Variant(
            contig,
            start_pos,
            str(ref),
            str(alt),
            genome)

        # keep metadata about the variant and its TCGA annotation
        metadata[variant] = {
            'Hugo_Symbol': x.Hugo_Symbol,
            'Center': x.Center,
            'Strand': x.Strand,
            'Variant_Classification': x.Variant_Classification,
            'Variant_Type': x.Variant_Type,
            'dbSNP_RS': x.dbSNP_RS,
            'dbSNP_Val_Status': x.dbSNP_Val_Status,
            'Tumor_Sample_Barcode': x.Tumor_Sample_Barcode,
            'Matched_Norm_Sample_Barcode': x.Matched_Norm_Sample_Barcode,
        }
        for optional_col in optional_cols:
            if optional_col in x:
                metadata[variant][optional_col] = x[optional_col]

        variants.append(variant)

    return VariantCollection(
        variants=variants,
        source_to_metadata_dict={path: metadata},
        sort_key=sort_key,
        distinct=distinct)

varcode.load_maf_dataframe(path, nrows=None, raise_on_error=True, encoding=None)

Load the guaranteed columns of a TCGA MAF file into a DataFrame

PARAMETER	DESCRIPTION
path	Path to MAF file TYPE: str
nrows	Optional limit to number of rows loaded TYPE: int DEFAULT: None
raise_on_error	Raise an exception upon encountering an error or log an error TYPE: bool DEFAULT: True
encoding	Encoding to use for UTF when reading MAF file. TYPE: str DEFAULT: None

Source code in varcode/maf.py

def load_maf_dataframe(path, nrows=None, raise_on_error=True, encoding=None):
    """
    Load the guaranteed columns of a TCGA MAF file into a DataFrame

    Parameters
    ----------
    path : str
        Path to MAF file

    nrows : int
        Optional limit to number of rows loaded

    raise_on_error : bool
        Raise an exception upon encountering an error or log an error

    encoding : str, optional
        Encoding to use for UTF when reading MAF file.
    """
    require_string(path, "Path to MAF")

    n_basic_columns = len(MAF_COLUMN_NAMES)

    # pylint: disable=no-member
    # pylint gets confused by read_csv
    df = pandas.read_csv(
        path,
        comment="#",
        sep="\t",
        low_memory=False,
        skip_blank_lines=True,
        header=0,
        nrows=nrows,
        encoding=encoding)

    if len(df.columns) < n_basic_columns:
        error_message = (
            "Too few columns in MAF file %s, expected %d but got  %d : %s" % (
                path, n_basic_columns, len(df.columns), df.columns))
        if raise_on_error:
            raise ValueError(error_message)
        else:
            logging.warn(error_message)

    # check each pair of expected/actual column names to make sure they match
    for expected, actual in zip(MAF_COLUMN_NAMES, df.columns):
        if expected != actual:
            # MAFs in the wild have capitalization differences in their
            # column names, normalize them to always use the names above
            if expected.lower() == actual.lower():
                # using DataFrame.rename in Python 2.7.x doesn't seem to
                # work for some files, possibly because Pandas treats
                # unicode vs. str columns as different?
                df[expected] = df[actual]
                del df[actual]
            else:
                error_message = (
                    "Expected column %s but got %s" % (expected, actual))
                if raise_on_error:
                    raise ValueError(error_message)
                else:
                    logging.warn(error_message)

    return df

Exceptions

varcode.ReferenceMismatchError(variant, transcript, expected_ref, observed_ref, transcript_offset=None, genome_start=None, genome_end=None)

Bases: ValueError

Raised when a variant's reported ref allele does not match the reference genome at the variant's position.

This most often means one of:

• The variant was called against a different reference build than the one being used for annotation (e.g. GRCh37 vs GRCh38).
• The variant's ref field was populated with the patient's germline allele rather than the canonical reference. VCF requires the ref field to match the reference genome; germline variants at the same position should be encoded as separate variants.
• Strand confusion: the variant is specified on the negative strand but varcode expects positive-strand coordinates.

Callers who would rather continue past this error can pass raise_on_error=False to :meth:Variant.effects to receive Failure effects instead.

Source code in varcode/errors.py

def __init__(self, variant, transcript, expected_ref, observed_ref,
             transcript_offset=None, genome_start=None, genome_end=None):
    self.variant = variant
    self.transcript = transcript
    self.expected_ref = expected_ref
    self.observed_ref = observed_ref
    self.transcript_offset = transcript_offset
    self.genome_start = genome_start
    self.genome_end = genome_end

    location = ""
    if transcript_offset is not None:
        location = " at transcript offset %d" % transcript_offset
    if genome_start is not None and genome_end is not None:
        location += " (chromosome positions %d:%d)" % (
            genome_start, genome_end)

    message = (
        "Reference allele mismatch for %s on %s%s: variant reports "
        "ref=%r but the reference genome has %r at this position.\n"
        "This usually means the variant was called against a "
        "different genome build, the ref field was filled in with "
        "the patient's germline allele rather than the reference, "
        "or the variant is on the wrong strand. Pass "
        "raise_on_error=False to .effects() to receive a Failure "
        "effect instead of raising." % (
            variant, transcript, location, observed_ref, expected_ref)
    )
    super().__init__(message)

varcode.SampleNotFoundError

Bases: KeyError

Raised when genotype info is requested for a sample that isn't present in the VariantCollection's source VCF(s).

varcode.GenomeBuildMismatchError(somatic_reference, germline_reference)

Bases: ValueError

Raised when a germline VCF and a somatic VCF were called against different reference genome builds (e.g. GRCh37 vs GRCh38). Effect coordinates from the two VCFs cannot be meaningfully composed.

Subclasses :class:ValueError so callers that already catch ValueError for ReferenceMismatchError continue to work. Set validate_reference=False on the call site if the user has explicitly lifted over one VCF into the other build and knows what they're doing.

Source code in varcode/germline.py

def __init__(self, somatic_reference, germline_reference):
    self.somatic_reference = somatic_reference
    self.germline_reference = germline_reference
    super().__init__(
        "Genome build mismatch: somatic VCF uses %r, germline "
        "context uses %r. Effect coordinates from the two cannot "
        "be composed without lift-over. If you've already lifted "
        "over one VCF into the other build, pass "
        "validate_reference=False to skip this check." % (
            somatic_reference, germline_reference))