Phasing and germline API

See Phasing and Germline-aware annotation for input requirements and examples.

Phasing

varcode.ReadPhasingSource

Bases: Protocol

Reports per-variant read-level evidence and co-occurring partners.

The minimum interface needed to establish cis from read-level data — co-observation on the same supporting reads, same long-read fragment, same assembled contig, etc. Implementations decide how strong the evidence is; consumers only see a boolean membership question.

Partners alone never establish trans: a variant missing from another's partners may simply be uncovered or outside the source's input. A source that observes one variant's alt allele together with the other's reference allele reports trans by also implementing in_cis(v1, v2, transcript=None), which returns True, False or None.

has_evidence(variant) -> bool

True if this source has any alt-supporting evidence for variant.

Source code in varcode/phasing.py

def has_evidence(self, variant) -> bool:
    """True if this source has any alt-supporting evidence for
    ``variant``."""
    ...

partners_in_cis(variant) -> Sequence

Variants observed in cis with variant — i.e. on the same supporting reads / fragment / contig. May include germline SNPs, nearby somatic variants, or variant itself (implementations pick their convention). Empty sequence when no evidence covers variant.

Source code in varcode/phasing.py

def partners_in_cis(self, variant) -> Sequence:
    """Variants observed in cis with ``variant`` — i.e. on the
    same supporting reads / fragment / contig. May include
    germline SNPs, nearby somatic variants, or ``variant`` itself
    (implementations pick their convention). Empty sequence when
    no evidence covers ``variant``."""
    ...

varcode.MutantTranscriptSource

Bases: Protocol

Reports an observed mutant transcript for (variant, transcript).

Independent of phasing. A source can implement just :class:ReadPhasingSource, just :class:MutantTranscriptSource, or both. Consumers iterating effects use this channel to substitute the observed mutant protein (from RNA assembly, long-read calling, etc.) for the reference-inferred one.

mutant_transcript(variant, transcript)

The :class:~varcode.MutantTranscript for (variant, transcript), or None when this source has no observed transcript for that pair.

Source code in varcode/phasing.py

def mutant_transcript(self, variant, transcript):
    """The :class:`~varcode.MutantTranscript` for
    ``(variant, transcript)``, or ``None`` when this source has
    no observed transcript for that pair."""
    ...

varcode.MolecularPhaseResolver(phasing_source: ReadPhasingSource)

Phase resolver backed by a :class:ReadPhasingSource (#269, #259).

Two variants are cis if the source reports them as co-observed. That's direct molecular evidence — not a probabilistic call. Not being co-observed is not evidence of trans, so without the source's own in_cis the answer is then None.

Usage::

# Any object satisfying ReadPhasingSource works. Common
# implementation: an Isovar adapter shipped by openvax/isovar.
resolver = MolecularPhaseResolver(source)
effects = variants.effects(phase_resolver=resolver)

If source also satisfies :class:MutantTranscriptSource, the resolver routes :meth:mutant_transcript calls to it — so any effect whose (variant, transcript) is covered by the source gets its :attr:~MutationEffect.mutant_transcript populated with the observed mutant transcript.

Sources may also expose their own in_cis(v1, v2, transcript=None) method. When present, the resolver delegates to that method instead of reducing through :meth:ReadPhasingSource.partners_in_cis. This is how direct-read sources report trans from co-covering fragments and None for pairs without enough of them.

Source code in varcode/phasing.py

def __init__(self, phasing_source: ReadPhasingSource):
    self.phasing_source = phasing_source
    self.phase_source = getattr(
        phasing_source, "source", self.phase_source)

has_evidence(variant) -> bool

Convenience passthrough to the wrapped source.

Source code in varcode/phasing.py

def has_evidence(self, variant) -> bool:
    """Convenience passthrough to the wrapped source."""
    return self.phasing_source.has_evidence(variant)

mutant_transcript(variant, transcript)

Return the observed :class:MutantTranscript for (variant, transcript), or None when the wrapped source doesn't implement :class:MutantTranscriptSource or has no transcript for that pair.

Source code in varcode/phasing.py

def mutant_transcript(self, variant, transcript):
    """Return the observed :class:`MutantTranscript` for
    ``(variant, transcript)``, or ``None`` when the wrapped source
    doesn't implement :class:`MutantTranscriptSource` or has no
    transcript for that pair."""
    get = getattr(self.phasing_source, "mutant_transcript", None)
    if get is None:
        return None
    return get(variant, transcript)

in_cis(v1, v2, transcript=None) -> Optional[bool]

Return the source's own in_cis answer when it has one.

Otherwise True when the source reports v1 and v2 as co-observed and None when it doesn't: a variant missing from the other's partners may be uncovered, so partners alone never establish trans.

transcript is accepted for interface symmetry with :class:VCFPhaseResolver.in_cis and passed to the source's own in_cis; the partners fallback is isoform-agnostic.

Source code in varcode/phasing.py

def in_cis(self, v1, v2, transcript=None) -> Optional[bool]:
    """Return the source's own ``in_cis`` answer when it has one.

    Otherwise ``True`` when the source reports ``v1`` and ``v2`` as
    co-observed and ``None`` when it doesn't: a variant missing from
    the other's partners may be uncovered, so partners alone never
    establish trans.

    ``transcript`` is accepted for interface symmetry with
    :class:`VCFPhaseResolver.in_cis` and passed to the source's own
    ``in_cis``; the partners fallback is isoform-agnostic.
    """
    source = self.phasing_source
    source_in_cis = getattr(source, "in_cis", None)
    if source_in_cis is not None:
        return source_in_cis(v1, v2, transcript=transcript)
    if ((source.has_evidence(v1) and v2 in source.partners_in_cis(v1))
            or (source.has_evidence(v2) and v1 in source.partners_in_cis(v2))):
        return True
    return None

phased_partners(variant, transcript=None) -> Sequence

Variants co-observed with variant — i.e. the cis set. Empty when the source has no evidence for variant.

Source code in varcode/phasing.py

def phased_partners(self, variant, transcript=None) -> Sequence:
    """Variants co-observed with ``variant`` — i.e. the cis set.
    Empty when the source has no evidence for ``variant``."""
    if not self.phasing_source.has_evidence(variant):
        return ()
    return tuple(self.phasing_source.partners_in_cis(variant))

varcode.ReadPhaseResolver(phasing_source: ReadPhasingSource)

Bases: MolecularPhaseResolver

Compatibility name for :class:MolecularPhaseResolver.

The old name is kept because it shipped in varcode 5.0.0. New code should prefer :class:MolecularPhaseResolver, which better describes sources backed by reads, paired fragments, assembled contigs, or other direct molecular evidence.

Source code in varcode/phasing.py

def __init__(self, phasing_source: ReadPhasingSource):
    self.phasing_source = phasing_source
    self.phase_source = getattr(
        phasing_source, "source", self.phase_source)

varcode.RNAReadPhasingSource(bam_path: str, *, variants=None, min_mapping_quality: int = 20, min_base_quality: int = 20, min_alt_reads: int = 2, phasing_error_rate: float = 0.01, max_p_value_for_phasing: float = 0.05, max_distance_from_read_edge: Optional[int] = 5, require_proper_pair: bool = True, skip_duplicates: bool = True, skip_secondary: bool = True, skip_supplementary: bool = True)

BAM-backed phasing source for RNA read co-occurrence.

This is the lightweight alternative to an Isovar-style assembly source. It reads quality-filtered alignments from an RNA-seq BAM and answers whether variants are observed on the same read or paired-end fragment. It does not assemble contigs and does not provide mutant_transcript; callers that need observed mutant proteins should use an assembly-backed source.

Usage::

source = RNAReadPhasingSource("tumor.rna.bam")
resolver = MolecularPhaseResolver(source)
effects = variants.effects(phase_resolver=resolver)

PARAMETER	DESCRIPTION
bam_path	Coordinate-sorted, indexed BAM path. TYPE: str
variants	Optional universe used by :meth:partners_in_cis. Variants seen through :meth:has_evidence or :meth:in_cis are registered automatically, so this is mainly a convenience for callers that query MolecularPhaseResolver.phased_partners directly. TYPE: sequence DEFAULT: None
min_mapping_quality	Minimum MAPQ for reads contributing evidence. TYPE: int DEFAULT: 20
min_base_quality	Minimum base quality for SNV/MNV and insertion allele calls. TYPE: int DEFAULT: 20
min_alt_reads	Minimum alt-supporting reads/fragments required for has_evidence, and fragments an allele combination needs to count toward a cis/trans call. TYPE: int DEFAULT: 2
phasing_error_rate	Chance that a read shows the wrong allele at a variant's locus, from sequencing error, RNA editing or mismapping. Raise it for error-prone reads, such as ONT's. TYPE: float DEFAULT: 0.01
max_p_value_for_phasing	One-sided binomial p-value at or below which :meth:in_cis counts an allele combination as more than such errors. TYPE: float DEFAULT: 0.05
max_distance_from_read_edge	Discard allele calls whose queried bases are closer than this many bases to either read edge. Set to None to disable. TYPE: int DEFAULT: 5
require_proper_pair	For paired reads, discard fragments not marked proper pair. Unpaired reads are still accepted. TYPE: bool DEFAULT: True
skip_duplicates	Standard SAM flag filters. TYPE: bool DEFAULT: True
skip_secondary	Standard SAM flag filters. TYPE: bool DEFAULT: True
skip_supplementary	Standard SAM flag filters. TYPE: bool DEFAULT: True

Source code in varcode/rna_read_phasing.py

def __init__(
        self,
        bam_path: str,
        *,
        variants=None,
        min_mapping_quality: int = 20,
        min_base_quality: int = 20,
        min_alt_reads: int = 2,
        phasing_error_rate: float = 0.01,
        max_p_value_for_phasing: float = 0.05,
        max_distance_from_read_edge: Optional[int] = 5,
        require_proper_pair: bool = True,
        skip_duplicates: bool = True,
        skip_secondary: bool = True,
        skip_supplementary: bool = True):
    try:
        import pysam
    except ImportError as e:
        raise ImportError(
            "RNAReadPhasingSource requires pysam. Install with "
            "`pip install varcode[rna]`.") from e
    for name, value in (("phasing_error_rate", phasing_error_rate),
                        ("max_p_value_for_phasing", max_p_value_for_phasing)):
        if not 0 < value < 1:
            raise ValueError("%s must be between 0 and 1, not %r" % (name, value))
    self._pysam = pysam
    self._bam = pysam.AlignmentFile(bam_path, "rb")
    self.min_mapping_quality = min_mapping_quality
    self.min_base_quality = min_base_quality
    self.min_alt_reads = min_alt_reads
    self.phasing_error_rate = phasing_error_rate
    self.max_p_value_for_phasing = max_p_value_for_phasing
    self.max_distance_from_read_edge = max_distance_from_read_edge
    self.require_proper_pair = require_proper_pair
    self.skip_duplicates = skip_duplicates
    self.skip_secondary = skip_secondary
    self.skip_supplementary = skip_supplementary
    self._known_variants = []
    self._known_variant_keys = set()
    self._support_cache = {}
    self._phase_cache = {}
    self._contig_cache = {}
    self._haplotypes = []
    if variants is not None:
        self.register_variants(variants)

close()

Close the underlying BAM handle.

Source code in varcode/rna_read_phasing.py

def close(self):
    """Close the underlying BAM handle."""
    self._bam.close()

register_variants(variants)

Register variants used by :meth:partners_in_cis.

This is optional for ordinary MolecularPhaseResolver.in_cis use, where both queried variants are registered automatically.

Source code in varcode/rna_read_phasing.py

def register_variants(self, variants):
    """Register variants used by :meth:`partners_in_cis`.

    This is optional for ordinary ``MolecularPhaseResolver.in_cis`` use,
    where both queried variants are registered automatically.
    """
    for variant in variants:
        self._register_variant(variant)

register_haplotype(variants, *, flanking_bases=5)

Test a known local allele combination by its anchored RNA sequence.

PARAMETER	DESCRIPTION
variants	Nonoverlapping substitutions/deletions on one genome and contig. All alternate alleles form one hypothesis, not an assumed phase. Register competing hypotheses separately when necessary. TYPE: sequence of Variant
flanking_bases	Unchanged reference bases required on each side (default five). TYPE: int DEFAULT: 5

Notes

A match requires the entire sequence, including both flanks, on ONE alignment, with the configured base-quality and read-edge filters. Equivalent D/N/split-gap encodings can then support the same known sequence. This does not establish a genomic deletion from an RNA skip. A fully covered high-quality mate with a different sequence vetoes that fragment's support. Missing coverage and low-quality context do not. Reference is fetched from the variants' genome; unavailable sequence, overlapping edits and reference mismatches raise ValueError. No BAM bases are corrected, no missing bases filled, and no mates assembled. Registered variants use these full-context matches instead of individual CIGAR calls. Nonmatches remain unknown, not reference/trans evidence.

Source code in varcode/rna_read_phasing.py

def register_haplotype(self, variants, *, flanking_bases=5):
    """Test a known local allele combination by its anchored RNA sequence.

    Parameters
    ----------
    variants : sequence of Variant
        Nonoverlapping substitutions/deletions on one genome and contig.
        All alternate alleles form one hypothesis, not an assumed phase.
        Register competing hypotheses separately when necessary.
    flanking_bases : int
        Unchanged reference bases required on each side (default five).

    Notes
    -----
    A match requires the entire sequence, including both flanks, on ONE
    alignment, with the configured base-quality and read-edge filters.
    Equivalent D/N/split-gap encodings can then support the same known
    sequence. This does not establish a genomic deletion from an RNA skip.
    A fully covered high-quality mate with a different sequence vetoes that
    fragment's support. Missing coverage and low-quality context do not.
    Reference is fetched from the variants' genome; unavailable sequence,
    overlapping edits and reference mismatches raise ValueError. No BAM
    bases are corrected, no missing bases filled, and no mates assembled.
    Registered variants use these full-context matches instead of individual
    CIGAR calls. Nonmatches remain unknown, not reference/trans evidence.
    """
    from .genome_sequence import reference_range

    variants = tuple(sorted(variants, key=lambda v: v.start))
    if not variants or isinstance(flanking_bases, bool) or not isinstance(flanking_bases, int) or flanking_bases < 1:
        raise ValueError("A nonempty haplotype and positive flanking_bases are required")
    first = variants[0]
    if any(v.contig != first.contig or v.genome is not first.genome for v in variants):
        raise ValueError("Haplotype variants must share a genome dataset and contig")
    if any(not v.ref or len(v.alt) > len(v.ref) for v in variants):
        raise ValueError("Anchored haplotypes currently support substitutions and deletions")
    if any(a.end >= b.start for a, b in zip(variants, variants[1:])):
        raise ValueError("Haplotype edits must not overlap")
    left, right = first.start - flanking_bases, variants[-1].end + flanking_bases
    if left < 1:
        raise ValueError("Haplotype lacks a left genomic anchor")
    sequence = reference_range(first.genome, first.contig, left, right).upper()
    if len(sequence) != right - left + 1 or set(sequence) - set("ACGT"):
        raise ValueError("Unambiguous reference sequence is required across the haplotype")
    for variant in reversed(variants):
        start = variant.start - left
        end = start + len(variant.ref)
        if sequence[start:end] != variant.ref.upper():
            raise ValueError("Haplotype reference allele disagrees with its genome")
        sequence = sequence[:start] + variant.alt.upper() + sequence[end:]
    if set(sequence) - set("ACGT"):
        raise ValueError("Unambiguous alternate sequence is required")
    keys = frozenset(self._variant_key(v) for v in variants)
    record = (keys, left - 1, right - 1, sequence)
    if record not in self._haplotypes:
        self._haplotypes.append(record)
        self.register_variants(variants)
        self._support_cache.clear()
        self._phase_cache.clear()

supports_variant(variant) -> Optional[int]

Count quality-filtered reads/fragments supporting variant.alt.

Returns None when the variant's contig is absent from the BAM. Results are cached until a new haplotype is registered.

Source code in varcode/rna_read_phasing.py

def supports_variant(self, variant) -> Optional[int]:
    """Count quality-filtered reads/fragments supporting ``variant.alt``.

    Returns ``None`` when the variant's contig is absent from the BAM.
    Results are cached until a new haplotype is registered.
    """
    self._register_variant(variant)
    variant_key = self._variant_key(variant)
    if variant_key in self._support_cache:
        return self._support_cache[variant_key]
    reads = self._fetch_reads_for_variant(variant)
    if reads is None:
        self._support_cache[variant_key] = None
        return None
    calls_by_fragment = defaultdict(set)
    for read in reads:
        call = self._read_allele(read, variant)
        if call is not None:
            fragment_key = (read.get_tag("RG") if read.has_tag("RG") else "", read.query_name)
            calls_by_fragment[fragment_key].add(call)
    count = sum(calls == {"alt"} for calls in calls_by_fragment.values())
    self._support_cache[variant_key] = count
    return count

has_evidence(variant) -> bool

True if the BAM has enough alt-supporting reads/fragments.

Source code in varcode/rna_read_phasing.py

def has_evidence(self, variant) -> bool:
    """True if the BAM has enough alt-supporting reads/fragments."""
    count = self.supports_variant(variant)
    return count is not None and count >= self.min_alt_reads

in_cis(v1, v2, transcript=None) -> Optional[bool]

Return cis/trans from RNA read or fragment co-occurrence.

Decided by :func:four_gamete_phase from the fragments covering both loci: True when both alt alleles are seen together beyond read errors and not each alone, False when each is seen alone and not together, None otherwise.

Source code in varcode/rna_read_phasing.py

def in_cis(self, v1, v2, transcript=None) -> Optional[bool]:
    """Return cis/trans from RNA read or fragment co-occurrence.

    Decided by :func:`four_gamete_phase` from the fragments covering
    both loci: ``True`` when both alt alleles are seen together beyond
    read errors and not each alone, ``False`` when each is seen alone
    and not together, ``None`` otherwise.
    """
    return four_gamete_phase(
        *self._phase_counts(v1, v2), min_fragments=self.min_alt_reads,
        error_rate=self.phasing_error_rate,
        max_p_value=self.max_p_value_for_phasing)

partners_in_cis(variant) -> Sequence

Known registered variants observed in cis with variant.

Source code in varcode/rna_read_phasing.py

def partners_in_cis(self, variant) -> Sequence:
    """Known registered variants observed in cis with ``variant``."""
    self._register_variant(variant)
    partners = []
    for other in self._known_variants:
        if other == variant:
            continue
        if self.in_cis(variant, other) is True:
            partners.append(other)
    return tuple(partners)

varcode.VCFPhaseResolver(variant_collection, sample)

Phase resolver backed by VCF GT + PS FORMAT fields.

Reads the phase data that varcode's VCF loader already parses into :class:~varcode.Genotype (via #267): whether the GT delimiter was | (phased) or / (unphased), the PS phase-set identifier, and the per-haplotype allele indices in :attr:Genotype.alleles.

Two variants are cis when they sit in the same phase set on the same haplotype slot, trans when they sit in the same phase set on different slots, and the resolver returns None ("no evidence") for variants that aren't both phased, don't share a phase set, or lack called alleles.

Compatible with any tool that writes standard-shaped VCF: WhatsHap, HapCUT2, DeepVariant, GATK HaplotypeCaller, long-read callers (PEPPER-DeepVariant, Clair3), population phasers (SHAPEIT5, Eagle2). varcode doesn't care which one wrote the file — it only reads GT and PS.

Multi-allelic sites are handled: varcode splits those rows into one :class:~varcode.Variant per ALT, each with an alt_allele_index preserved on the :class:~varcode.VariantCollection metadata. The resolver maps each variant to its GT-encoded index and asks "which haplotype slot carries this specific alt?".

Single-sample by construction. Phase is per-sample; multi-sample VCFs need one resolver per sample.

Currently supplies the cis/trans query but does not attach a :class:~varcode.MutantTranscript — DNA phasing alone doesn't produce an assembled contig. The natural next step is a HaplotypeEffect / multi-variant apply_variants_to_transcript helper that, when two or more cis variants overlap the same transcript, builds a single joint :class:MutantTranscript applying all edits at once. That's a separate PR — this resolver already has the inputs it needs (in_cis) to drive the grouping.

Source code in varcode/phasing.py

def __init__(self, variant_collection, sample):
    self._collection = variant_collection
    self._sample = sample

in_cis(v1, v2, transcript=None) -> Optional[bool]

Return True if v1 and v2 are on the same haplotype in the same phase set, False if they're on different haplotypes in the same phase set, None when the phase relationship can't be determined (unphased GT, different phase sets, uncalled alleles).

transcript is accepted for interface symmetry with :class:MolecularPhaseResolver.in_cis but isn't consulted — DNA-level phase is isoform-agnostic.

Source code in varcode/phasing.py

def in_cis(self, v1, v2, transcript=None) -> Optional[bool]:
    """Return ``True`` if ``v1`` and ``v2`` are on the same
    haplotype in the same phase set, ``False`` if they're on
    different haplotypes in the same phase set, ``None`` when
    the phase relationship can't be determined (unphased GT,
    different phase sets, uncalled alleles).

    ``transcript`` is accepted for interface symmetry with
    :class:`MolecularPhaseResolver.in_cis` but isn't consulted —
    DNA-level phase is isoform-agnostic.
    """
    g1 = self._genotype(v1)
    g2 = self._genotype(v2)
    if g1 is None or g2 is None:
        return None
    alt1 = self._collection._alt_index_for(v1)
    alt2 = self._collection._alt_index_for(v2)
    # Homozygous-alt on either side makes phase deterministic:
    # every haplotype carries this alt, so it's cis with anything
    # the OTHER variant sits on. Shortcut before requiring
    # phased=True.
    hom1 = self._is_homozygous_alt(g1, alt1)
    hom2 = self._is_homozygous_alt(g2, alt2)
    if hom1 and hom2:
        return True
    if hom1:
        # v1 is on every haplotype → cis with v2 iff v2 is called.
        return self._haplotype_slot(g2, v2) is not None
    if hom2:
        return self._haplotype_slot(g1, v1) is not None
    # Otherwise both must be phased and in the same phase set.
    if not g1.phased or not g2.phased:
        return None
    if g1.phase_set != g2.phase_set or g1.phase_set is None:
        return None
    s1 = self._haplotype_slot(g1, v1)
    s2 = self._haplotype_slot(g2, v2)
    if s1 is None or s2 is None:
        return None
    return s1 == s2

phased_partners(variant, transcript=None)

Variants in the collection that are cis with variant under this resolver — i.e. sit in the same phase set on the same haplotype slot. Empty when variant isn't phased or has no called alt in the sample.

Source code in varcode/phasing.py

def phased_partners(self, variant, transcript=None):
    """Variants in the collection that are cis with ``variant``
    under this resolver — i.e. sit in the same phase set on the
    same haplotype slot. Empty when ``variant`` isn't phased or
    has no called alt in the sample.
    """
    partners = []
    for other in self._collection:
        if other == variant:
            continue
        if self.in_cis(variant, other) is True:
            partners.append(other)
    return tuple(partners)

varcode.apply_phase_resolver_to_effects(effects, phase_resolver)

Post-process an :class:EffectCollection (or any iterable of :class:MutationEffect) to attach observed :class:MutantTranscript objects when the resolver has evidence.

Mutates each effect in place by setting effect.mutant_transcript. Effects whose transcript isn't resolvable or whose (variant, transcript) has no observed transcript are left untouched — so this is safe to call on a mixed collection where only some variants have RNA evidence.

Source code in varcode/phasing.py

def apply_phase_resolver_to_effects(effects, phase_resolver):
    """Post-process an :class:`EffectCollection` (or any iterable of
    :class:`MutationEffect`) to attach observed
    :class:`MutantTranscript` objects when the resolver has evidence.

    Mutates each effect in place by setting
    ``effect.mutant_transcript``. Effects whose transcript isn't
    resolvable or whose ``(variant, transcript)`` has no observed
    transcript are left untouched — so this is safe to call on a mixed
    collection where only some variants have RNA evidence.
    """
    if phase_resolver is None:
        return effects
    if not hasattr(phase_resolver, "mutant_transcript"):
        return effects
    for e in effects:
        transcript = getattr(e, "transcript", None)
        variant = getattr(e, "variant", None)
        if variant is None or transcript is None:
            continue
        mt = phase_resolver.mutant_transcript(variant, transcript)
        if mt is not None:
            # Intentional mutation: the effect's mutant_transcript
            # slot was either None (point variants, cryptic stubs,
            # etc.) or populated from DNA-only inference. An observed
            # mutant transcript is higher-confidence evidence, so it
            # wins.
            e.mutant_transcript = mt
    return effects

varcode.query_in_cis(phase_resolver, v1, v2, transcript=None) -> Optional[bool]

Ask phase_resolver whether v1 and v2 are in cis.

Returns True (cis), False (trans) or None (unknown). NumPy booleans count as booleans. No resolver, a non-boolean answer, or a resolver error gives None; errors are logged as warnings so a broken resolver is visible without failing annotation. transcript is passed only when given, so two-argument resolvers keep working.

Source code in varcode/phasing.py

def query_in_cis(phase_resolver, v1, v2, transcript=None) -> Optional[bool]:
    """Ask ``phase_resolver`` whether ``v1`` and ``v2`` are in cis.

    Returns ``True`` (cis), ``False`` (trans) or ``None`` (unknown). NumPy
    booleans count as booleans. No resolver, a non-boolean answer, or a
    resolver error gives ``None``; errors are logged as warnings so a
    broken resolver is visible without failing annotation. ``transcript``
    is passed only when given, so two-argument resolvers keep working.
    """
    in_cis = getattr(phase_resolver, "in_cis", None)
    if in_cis is None:
        return None
    try:
        answer = (in_cis(v1, v2) if transcript is None
                  else in_cis(v1, v2, transcript=transcript))
    except Exception as error:
        logger.warning(
            "Phase resolver %s failed for %s and %s: %s",
            type(phase_resolver).__name__, v1, v2, error)
        return None
    return bool(answer) if answer in (True, False) else None

Germline-aware annotation

varcode.GermlineContext(variants: 'VariantCollection', completeness: Completeness = Completeness.COMPLETE, reference_name: Optional[str] = None, metadata: Mapping[str, Any] = dict(), max_phase_hypotheses: int = DEFAULT_MAX_PHASE_HYPOTHESES) dataclass

The patient's germline, packaged with completeness metadata and reference-build info for cross-VCF validation.

Construct via the from_* classmethods rather than instantiating directly; the constructors apply the input-shape-specific validation each route needs.

ATTRIBUTE	DESCRIPTION
variants	The germline variants as a :class:~varcode.VariantCollection. Empty for :meth:empty contexts. TYPE: 'VariantCollection'
completeness	How to interpret absence-of-a-call (see :class:Completeness). TYPE: Completeness
reference_name	The genome reference these variants were called against — "GRCh37", "GRCh38", "hg19", etc. None when the context is empty or the reference is unknown. Used by :meth:validate_against to fail fast on cross-VCF build mismatches. TYPE: Optional[str]
metadata	Open-ended dict for caller-supplied annotations (source caller name, sample identifier, normalization tool, etc.). Not interpreted by varcode; rides along for downstream consumers and serialization. TYPE: Mapping[str, Any]
max_phase_hypotheses	Most cis/trans phase hypotheses to enumerate for one somatic variant (default 8, i.e. up to three unphased germline variants nearby). Beyond it, effects are :class:~varcode.HypothesisLimit. Change it with dataclasses.replace(ctx, max_phase_hypotheses=16). TYPE: int

Examples:

Route 1 — full germline call set::

ctx = GermlineContext.from_germline_vcf("normal.vcf")

Route 2 — multi-sample VCF, extract a column. The user must declare completeness explicitly because absence-from-a-multi- sample column rarely means ref/ref::

ctx = GermlineContext.from_multi_sample_vcf(
    "merged.vcf", sample="NORMAL", completeness=Completeness.SPARSE)

Direct construction (tests, custom pipelines)::

ctx = GermlineContext.from_variants(
    germline_variants, completeness=Completeness.COMPLETE,
    reference_name="GRCh38")

Explicit empty context — opt-in to reference-relative fallback::

ctx = GermlineContext.empty()

from_germline_vcf(path: str, *, completeness: Completeness = Completeness.COMPLETE, metadata: Optional[Mapping[str, Any]] = None, **load_vcf_kwargs) -> 'GermlineContext' classmethod

Load a full germline VCF into a context.

load_vcf_kwargs are passed through to :func:varcode.load_vcf — for example genome= or only_passing=False. The returned context defaults to Completeness.COMPLETE; pass completeness= only if the VCF is something other than a real germline call set.

Source code in varcode/germline.py

@classmethod
def from_germline_vcf(
        cls,
        path: str,
        *,
        completeness: Completeness = Completeness.COMPLETE,
        metadata: Optional[Mapping[str, Any]] = None,
        **load_vcf_kwargs) -> "GermlineContext":
    """Load a full germline VCF into a context.

    ``load_vcf_kwargs`` are passed through to
    :func:`varcode.load_vcf` — for example ``genome=`` or
    ``only_passing=False``. The returned context defaults to
    ``Completeness.COMPLETE``; pass ``completeness=`` only if the
    VCF is something other than a real germline call set.
    """
    # Lazy import to avoid pulling vcf.py into the import graph
    # for callers who don't need it.
    from .vcf import load_vcf
    vc = load_vcf(path, **load_vcf_kwargs)
    if len(vc) == 0:
        warnings.warn(
            "Loaded germline VCF %r contains zero variants. This "
            "is almost always a wrong-file error; effect "
            "prediction will silently fall through to "
            "reference-relative on every somatic variant." % path)
    return cls(
        variants=vc,
        completeness=completeness,
        reference_name=cls._reference_name_of(vc),
        metadata=dict(metadata or {}),
    )

from_multi_sample_vcf(path: str, sample: str, *, completeness: Completeness, metadata: Optional[Mapping[str, Any]] = None, **load_vcf_kwargs) -> 'GermlineContext' classmethod

Load a multi-sample VCF and extract one sample's calls as the germline.

completeness is required (no default) — multi-sample VCFs from somatic callers (Mutect2's NORMAL column, e.g.) are almost always sparse, but pure-germline multi-sample VCFs (1000G, gnomAD batch genotyping) are complete. Forcing the caller to declare prevents subtle correctness bugs from treating a sparse column as if absence implied ref/ref.

The sample is filtered post-load. If you need per-sample zygosity information, pass include_info=True (the default) and consult vc.metadata[variant]["sample_info"][sample] downstream.

Source code in varcode/germline.py

@classmethod
def from_multi_sample_vcf(
        cls,
        path: str,
        sample: str,
        *,
        completeness: Completeness,
        metadata: Optional[Mapping[str, Any]] = None,
        **load_vcf_kwargs) -> "GermlineContext":
    """Load a multi-sample VCF and extract one sample's calls as
    the germline.

    ``completeness`` is required (no default) — multi-sample VCFs
    from somatic callers (Mutect2's ``NORMAL`` column, e.g.) are
    almost always sparse, but pure-germline multi-sample VCFs
    (1000G, gnomAD batch genotyping) are complete. Forcing the
    caller to declare prevents subtle correctness bugs from
    treating a sparse column as if absence implied ref/ref.

    The sample is filtered post-load. If you need per-sample
    zygosity information, pass ``include_info=True`` (the default)
    and consult ``vc.metadata[variant]["sample_info"][sample]``
    downstream.
    """
    from .vcf import load_vcf
    vc = load_vcf(path, **load_vcf_kwargs)
    if sample not in cls._samples_of(vc):
        from .errors import SampleNotFoundError
        raise SampleNotFoundError(
            "Sample %r not present in %s. Available samples: %s" % (
                sample, path, sorted(cls._samples_of(vc))))
    # Filter to variants where the named sample has a non-ref call.
    # The base VariantCollection isn't intrinsically sample-aware;
    # this is a conservative subset that the caller can refine.
    filtered = cls._variants_called_in_sample(vc, sample)
    meta = dict(metadata or {})
    meta.setdefault("source_path", path)
    meta.setdefault("sample", sample)
    return cls(
        variants=filtered,
        completeness=completeness,
        reference_name=cls._reference_name_of(vc),
        metadata=meta,
    )

from_variants(variants, *, completeness: Completeness = Completeness.COMPLETE, reference_name: Optional[str] = None, metadata: Optional[Mapping[str, Any]] = None) -> 'GermlineContext' classmethod

Construct from an already-built :class:VariantCollection or any iterable of :class:Variant objects.

Useful for tests, hand-built pipelines, and downstream tools that already have variants in memory and don't need to re-parse a VCF. reference_name should be passed explicitly when not carried by the variants themselves; otherwise cross-VCF validation will be a no-op.

Source code in varcode/germline.py

@classmethod
def from_variants(
        cls,
        variants,
        *,
        completeness: Completeness = Completeness.COMPLETE,
        reference_name: Optional[str] = None,
        metadata: Optional[Mapping[str, Any]] = None) -> "GermlineContext":
    """Construct from an already-built :class:`VariantCollection`
    or any iterable of :class:`Variant` objects.

    Useful for tests, hand-built pipelines, and downstream tools
    that already have variants in memory and don't need to re-parse
    a VCF. ``reference_name`` should be passed explicitly when not
    carried by the variants themselves; otherwise cross-VCF
    validation will be a no-op.
    """
    from .variant_collection import VariantCollection
    if isinstance(variants, VariantCollection):
        vc = variants
    else:
        vc = VariantCollection(list(variants))
    if reference_name is None:
        reference_name = cls._reference_name_of(vc)
    return cls(
        variants=vc,
        completeness=completeness,
        reference_name=reference_name,
        metadata=dict(metadata or {}),
    )

empty() -> 'GermlineContext' classmethod

Explicit no-germline context. Use this in pipelines where germline= is structurally required but the caller has no germline data — it documents intent better than passing None, and downstream code can rely on the kwarg always being a :class:GermlineContext.

Effect prediction with an empty context falls through to reference-relative annotation (no patient transcript construction), with no warnings — the caller has explicitly opted in to the fallback.

Source code in varcode/germline.py

@classmethod
def empty(cls) -> "GermlineContext":
    """Explicit no-germline context. Use this in pipelines where
    ``germline=`` is structurally required but the caller has no
    germline data — it documents intent better than passing
    ``None``, and downstream code can rely on the kwarg always
    being a :class:`GermlineContext`.

    Effect prediction with an empty context falls through to
    reference-relative annotation (no patient transcript
    construction), with no warnings — the caller has explicitly
    opted in to the fallback.
    """
    from .variant_collection import VariantCollection
    return cls(
        variants=VariantCollection([]),
        completeness=Completeness.EMPTY,
        reference_name=None,
        metadata={},
    )

homozygous(variants) -> Tuple

The variants this context's sample carries on both haplotypes.

Needs genotypes for one sample: metadata["sample"] or the only sample in :attr:variants. Without them nothing is reported.

Source code in varcode/germline.py

def homozygous(self, variants) -> Tuple:
    """The ``variants`` this context's sample carries on both haplotypes.

    Needs genotypes for one sample: ``metadata["sample"]`` or the only
    sample in :attr:`variants`. Without them nothing is reported.
    """
    from .genotype import Zygosity
    sample = self.metadata.get("sample")
    if sample is None:
        samples = self.variants.samples if self.variants is not None else ()
        if len(samples) != 1:
            return ()
        sample = samples[0]
    return tuple(
        v for v in variants
        if self.variants.zygosity(v, sample) is Zygosity.HOMOZYGOUS)

__bool__() -> bool

Truthy when there's something to apply. EMPTY contexts are falsy so if germline_context: reads idiomatically.

Source code in varcode/germline.py

def __bool__(self) -> bool:
    """Truthy when there's something to apply. ``EMPTY`` contexts
    are falsy so ``if germline_context:`` reads idiomatically."""
    return self.completeness is not Completeness.EMPTY and len(self.variants) > 0

validate_against(somatic, *, validate_reference: bool = True) -> None

Cross-validate this context with a somatic :class:VariantCollection. Hard error on reference-build mismatch unless validate_reference=False; warn on suspicious shapes (empty germline, sparse coverage with no overlap with somatic, etc.).

Called automatically by :meth:Variant.effects / :meth:VariantCollection.effects when a context is supplied; callers running validation manually can do so up front to fail fast before annotation.

Source code in varcode/germline.py

def validate_against(
        self,
        somatic,
        *,
        validate_reference: bool = True) -> None:
    """Cross-validate this context with a somatic
    :class:`VariantCollection`. Hard error on reference-build
    mismatch unless ``validate_reference=False``; warn on
    suspicious shapes (empty germline, sparse coverage with no
    overlap with somatic, etc.).

    Called automatically by :meth:`Variant.effects` /
    :meth:`VariantCollection.effects` when a context is supplied;
    callers running validation manually can do so up front to fail
    fast before annotation.
    """
    if not isinstance(self, GermlineContext):
        raise TypeError(
            "validate_against expects self to be a GermlineContext")
    somatic_ref = self._reference_name_of(somatic)
    if validate_reference and self.reference_name and somatic_ref:
        if self.reference_name != somatic_ref:
            raise GenomeBuildMismatchError(
                somatic_reference=somatic_ref,
                germline_reference=self.reference_name)
    if (self.completeness is not Completeness.EMPTY
            and len(self.variants) == 0):
        warnings.warn(
            "GermlineContext is non-empty by completeness flag (%s) "
            "but holds zero variants — likely a wrong-file or "
            "filter-too-aggressive error." % self.completeness.value)

variants_in_window(contig: str, start: int, end: int) -> Tuple

Germline variants overlapping [start, end] on contig (inclusive on both ends).

Used by the window-based lookup machinery (slice 2 of #268). Lazy interval index is built on first call and cached on the instance — subsequent calls are O(log N) per contig.

Returns a tuple (immutable) so callers can safely cache the result without worrying about the underlying index mutating.

Source code in varcode/germline.py

def variants_in_window(
        self,
        contig: str,
        start: int,
        end: int) -> Tuple:
    """Germline variants overlapping ``[start, end]`` on
    ``contig`` (inclusive on both ends).

    Used by the window-based lookup machinery (slice 2 of #268).
    Lazy interval index is built on first call and cached on the
    instance — subsequent calls are O(log N) per contig.

    Returns a tuple (immutable) so callers can safely cache the
    result without worrying about the underlying index mutating.
    """
    index = self._index()
    starts, variants = index.get(contig, ((), ()))
    if not starts:
        return ()
    # Variants are indexed by start. Find candidates whose
    # start <= end, then filter by their actual end span. Most
    # germline variants are SNVs / short indels so the candidate
    # window is small.
    cutoff = bisect.bisect_right(starts, end)
    result = []
    for i in range(cutoff):
        v = variants[i]
        v_end = getattr(v, "end", None) or v.start
        if v_end >= start:
            result.append(v)
    return tuple(result)

varcode.Completeness

Bases: Enum

How exhaustive the germline call set is — the load-bearing flag that pins what absence of a call at a position means.

The same data structure ("a list of germline variants") can come from very different pipelines, and downstream effect prediction cannot make the right call without knowing which:

• If a position is absent from a real germline VCF emitted by a germline caller that examined the entire normal BAM, the patient is ref/ref there. Effect prediction proceeds reference-relative at that codon.
• If a position is absent from the NORMAL column of a somatic-caller VCF, it likely means the somatic caller didn't emit a row — not that the position is ref/ref. The patient's germline state at that codon is unknown. The honest output is a possibility set including "unknown germline."
• If a position is absent from a panel-of-normals filter list, it definitely doesn't imply ref/ref — the file only lists curated hotspots.

Mis-treating "absent" as "ref/ref" silently produces wrong germline-aware effects on somatic variants in long stretches of the genome the somatic caller never touched. The flag exists so that mistake fails loud (or at least produces an honest possibility set) instead of silently corrupting clinical annotation.

Values

+-------------------+-----------------------------------+--------------------------+ | Value | Typical pipeline of origin | Absence at a position | +===================+===================================+==========================+ | :attr:COMPLETE | Germline caller (DeepVariant, | ⇒ ref/ref | | | HaplotypeCaller, Strelka2 | | | | germline) on the normal BAM | | +-------------------+-----------------------------------+--------------------------+ | :attr:SPARSE | NORMAL column of a somatic | ⇒ unknown (probably | | | tumor-vs-normal VCF (Mutect2, | ref/ref but not | | | Strelka2 somatic, VarScan2 | queried). Honest output| | | somatic) | is a possibility set. | +-------------------+-----------------------------------+--------------------------+ | :attr:HOTSPOTS_ | Panel-of-normals filter list, | ⇒ definitely unknown. | | ONLY | ClinVar pathogenic list, single- | Strictly weaker | | | hotspot allowlists | evidence than SPARSE. | +-------------------+-----------------------------------+--------------------------+ | :attr:EMPTY | "I have no germline data" — | n/a (no germline-aware | | | explicit fallback, used so users | logic runs; equivalent to| | | opt into reference-relative | not passing germline= at | | | annotation rather than getting it | all) | | | by accident from a missing kwarg | | +-------------------+-----------------------------------+--------------------------+

What downstream slices do with this

Slice 3 of #268 wires germline= through annotator dispatch. When a somatic variant lands in a transcript window that has no germline calls, the annotator reads this flag to decide between:

• COMPLETE → patient is ref/ref in this window; emit a single reference-relative effect.
• SPARSE / HOTSPOTS_ONLY → patient's germline is unknown in this window; emit a possibility set including the reference-relative effect plus "germline-unknown" outcomes so the user sees the uncertainty.
• EMPTY → no germline-aware logic; reference-relative.

Constructors and defaults

:meth:GermlineContext.from_germline_vcf defaults to COMPLETE because that's almost always what a real germline VCF is.

:meth:GermlineContext.from_multi_sample_vcf requires the caller to declare completeness explicitly (no default) — a multi-sample VCF could be either, and silently defaulting either direction is a correctness bug waiting to happen.

:meth:GermlineContext.empty always sets EMPTY.

varcode.predict_germline_aware_effect(somatic_variant, transcript, germline_ctx: GermlineContext, annotator, phase_resolver=None, window_fn=default_germline_window, max_hypotheses: Optional[int] = None)

Predict the effect of somatic_variant on transcript against the patient's germline-applied transcript.

Single entry point for germline-aware effect prediction. :func:varcode.effects.predict_variant_effects calls this whenever a non-empty :class:GermlineContext is supplied; otherwise it bypasses the germline path entirely and the existing annotator dispatch produces today's reference-relative output unchanged.

Behaviour by case:

• Somatic allele matches a germline allele — a GermlineAlleleOverlap, without inferring zygosity, LOH, or a new allele sequence.
• Phase cannot matter (no germline in the somatic's window, or a noncoding or incompletely annotated transcript) — delegate to annotator. With no germline nearby, SPARSE / HOTSPOTS_ONLY contexts mark the result effect.germline_unknown = True.
• Phase known (resolver answers, homozygous germline, or a hemizygous contig) — classify against the single patient haplotype via :func:_classify_against_patient_baseline.
• Phase unknown — classify each hypothesis and wrap them in a :class:~varcode.effects.effect_classes.PhaseCandidateSet.
• Too many hypotheses — a :class:~varcode.HypothesisLimit recording the phase partition and the reference-relative effect.

max_hypotheses overrides germline_ctx.max_phase_hypotheses. window_fn is the pluggable window selector — defaults to :func:default_germline_window (codon-level, with splice-signal expansion when the somatic is splice-adjacent). Callers that need different windows pass their own.

Source code in varcode/germline.py

def predict_germline_aware_effect(
        somatic_variant,
        transcript,
        germline_ctx: GermlineContext,
        annotator,
        phase_resolver=None,
        window_fn=default_germline_window,
        max_hypotheses: Optional[int] = None):
    """Predict the effect of ``somatic_variant`` on ``transcript``
    against the patient's germline-applied transcript.

    Single entry point for germline-aware effect prediction.
    :func:`varcode.effects.predict_variant_effects` calls this whenever
    a non-empty :class:`GermlineContext` is supplied; otherwise it
    bypasses the germline path entirely and the existing annotator
    dispatch produces today's reference-relative output unchanged.

    Behaviour by case:

    * **Somatic allele matches a germline allele** — a
      ``GermlineAlleleOverlap``, without inferring zygosity, LOH, or a
      new allele sequence.
    * **Phase cannot matter** (no germline in the somatic's window, or a
      noncoding or incompletely annotated transcript) — delegate to
      ``annotator``. With no germline nearby, SPARSE / HOTSPOTS_ONLY
      contexts mark the result ``effect.germline_unknown = True``.
    * **Phase known** (resolver answers, homozygous germline, or a
      hemizygous contig) — classify
      against the single patient haplotype via
      :func:`_classify_against_patient_baseline`.
    * **Phase unknown** — classify each hypothesis and wrap them in a
      :class:`~varcode.effects.effect_classes.PhaseCandidateSet`.
    * **Too many hypotheses** — a :class:`~varcode.HypothesisLimit`
      recording the phase partition and the reference-relative effect.

    ``max_hypotheses`` overrides ``germline_ctx.max_phase_hypotheses``.
    ``window_fn`` is the pluggable window selector — defaults to
    :func:`default_germline_window` (codon-level, with splice-signal
    expansion when the somatic is splice-adjacent). Callers that
    need different windows pass their own.
    """
    from pyensembl import Transcript
    from .effects.effect_classes import (
        GermlineAlleleOverlap,
        HypothesisLimit,
        PhaseCandidateSet,
    )

    limit = _check_max_hypotheses(
        germline_ctx.max_phase_hypotheses if max_hypotheses is None
        else max_hypotheses)
    if not isinstance(transcript, Transcript):
        # Mirrors annotator entry: SVs and other non-Transcript
        # consumers don't go through germline-aware prediction.
        return annotator.annotate_on_transcript(somatic_variant, transcript)

    contig, start, end = window_fn(somatic_variant, transcript)
    germline_in_window = germline_ctx.variants_in_window(contig, start, end)
    if detect_germline_overlap(somatic_variant, germline_in_window):
        return GermlineAlleleOverlap(somatic_variant, transcript)

    if not (germline_in_window and transcript.is_protein_coding
            and transcript.complete):
        effect = annotator.annotate_on_transcript(somatic_variant, transcript)
        if effect is NotImplemented:
            return effect
        if not germline_in_window and germline_ctx.completeness in (
                Completeness.SPARSE, Completeness.HOTSPOTS_ONLY):
            effect.germline_unknown = True
        return effect

    phase = partition_germline_by_phase(
        somatic_variant, germline_in_window, phase_resolver,
        homozygous=germline_ctx.homozygous(germline_in_window))
    try:
        hypotheses = phase.hypotheses(limit)
    except HypothesisLimitError:
        from .effects import predict_variant_effect_on_transcript
        return HypothesisLimit(
            somatic_variant, transcript, limit, phase,
            predict_variant_effect_on_transcript(
                somatic_variant, transcript, annotator=annotator))

    if len(hypotheses) == 1:
        effect = _classify_against_patient_baseline(
            somatic_variant, transcript, hypotheses[0])
        # Stash the phase metadata on the effect so consumers /
        # serializers can recover it. Single-hypothesis effects don't
        # need a possibility set, but the evidence is still useful.
        effect.germline_phase_state = hypotheses[0].phase_state
        effect.germline_variants_in_window = tuple(germline_in_window)
        return effect

    return PhaseCandidateSet(
        variant=somatic_variant,
        transcript=transcript,
        candidates=tuple(
            _classify_against_patient_baseline(somatic_variant, transcript, h)
            for h in hypotheses),
        hypotheses=hypotheses,
        germline_variants=tuple(germline_in_window))

varcode.apply_germline_to_transcript(transcript, germline_ctx, somatic_variant=None)

Apply germline variants from germline_ctx to transcript, returning the patient's baseline :class:MutantTranscript.

Lower-level entry point for callers that want the patient transcript directly without going through full effect prediction. Used internally by :func:predict_germline_aware_effect; exposed publicly for downstream tools (Isovar, Exacto) that want to compute a custom analysis on the patient haplotype.

Behaviour:

• If germline_ctx is empty, returns None.
• If somatic_variant is provided, restricts germline to the somatic's window (per :func:default_germline_window); else applies all germline variants overlapping any exon of the transcript.
• If germline edits conflict (overlapping cDNA ranges) or land outside the CDS, returns None and the caller falls back.

The returned object is the same shape that :func:varcode.mutant_transcript.apply_variants_to_transcript produces: a :class:MutantTranscript carrying the germline edits with mutant_protein_sequence populated when the edits land after the CDS start.

Source code in varcode/germline.py

def apply_germline_to_transcript(transcript, germline_ctx, somatic_variant=None):
    """Apply germline variants from ``germline_ctx`` to ``transcript``,
    returning the patient's baseline :class:`MutantTranscript`.

    Lower-level entry point for callers that want the patient
    transcript directly without going through full effect prediction.
    Used internally by :func:`predict_germline_aware_effect`; exposed
    publicly for downstream tools (Isovar, Exacto) that want to
    compute a custom analysis on the patient haplotype.

    Behaviour:

    * If ``germline_ctx`` is empty, returns ``None``.
    * If ``somatic_variant`` is provided, restricts germline to the
      somatic's window (per :func:`default_germline_window`); else
      applies all germline variants overlapping any exon of the
      transcript.
    * If germline edits conflict (overlapping cDNA ranges) or land
      outside the CDS, returns ``None`` and the caller falls back.

    The returned object is the same shape that
    :func:`varcode.mutant_transcript.apply_variants_to_transcript`
    produces: a :class:`MutantTranscript` carrying the germline
    edits with ``mutant_protein_sequence`` populated when the edits
    land after the CDS start.
    """
    if not germline_ctx:
        return None
    from .mutant_transcript import apply_variants_to_transcript
    if somatic_variant is not None:
        contig, start, end = default_germline_window(
            somatic_variant, transcript)
        germline_in_window = germline_ctx.variants_in_window(
            contig, start, end)
    else:
        # Whole-transcript window: cover every exon. Useful for
        # callers building a single patient transcript independent of
        # any specific somatic — e.g., to translate the patient
        # protein for cohort-level analysis.
        germline_in_window = []
        try:
            for exon in transcript.exons:
                germline_in_window.extend(
                    germline_ctx.variants_in_window(
                        exon.contig, exon.start, exon.end))
        except Exception:
            return None
    if not germline_in_window:
        return None
    return apply_variants_to_transcript(
        list(germline_in_window), transcript)

varcode.partition_germline_by_phase(somatic_variant, germline_variants, phase_resolver=None, homozygous=(), cis_variants=()) -> PhasePartition

Split germline_variants by their phase relative to somatic_variant.

Variants in homozygous sit on both haplotypes and are always cis, as is every germline variant on a hemizygous contig (chrM, chrY; chrX is treated as diploid, since Varcode has no sex information). The rest are placed by :func:~varcode.phasing.query_in_cis answers between every pair of variants, so phase can follow through another germline variant. Variants whose phase is known only relative to each other become linked groups. An answer that contradicts earlier ones is ignored with a warning.

cis_variants supplies additional variants known to be in cis with the somatic variant. They anchor the phase graph before resolver answers, so germline phase can follow through any member of that group. Homozygous alleles do not link the two haplotypes through this graph.

Source code in varcode/germline.py

def partition_germline_by_phase(
        somatic_variant, germline_variants, phase_resolver=None,
        homozygous=(), cis_variants=()) -> PhasePartition:
    """Split ``germline_variants`` by their phase relative to ``somatic_variant``.

    Variants in ``homozygous`` sit on both haplotypes and are always cis,
    as is every germline variant on a hemizygous contig (chrM, chrY; chrX
    is treated as diploid, since Varcode has no sex information). The
    rest are placed by :func:`~varcode.phasing.query_in_cis` answers
    between every pair of variants, so phase can follow through another
    germline variant. Variants whose phase is known only relative to each
    other become ``linked`` groups. An answer that contradicts earlier
    ones is ignored with a warning.

    ``cis_variants`` supplies additional variants known to be in cis with
    the somatic variant. They anchor the phase graph before resolver answers,
    so germline phase can follow through any member of that group.
    Homozygous alleles do not link the two haplotypes through this graph.
    """
    germline = tuple(germline_variants)
    homozygous = tuple(homozygous)
    if str(somatic_variant.contig).lstrip("chr").upper() in ("M", "MT", "Y"):
        return PhasePartition(germline=germline, cis=germline, implicit=True)
    rest = [g for g in germline if g not in homozygous]
    # Union-find over the somatic variant (node 0), ``rest``, and any
    # additional cis-group members. Parity is 0 for the same haplotype
    # as the parent node and 1 for the other one.
    nodes = [somatic_variant] + rest
    cis_variants = tuple(cis_variants)
    for variant in cis_variants:
        if not detect_germline_overlap(variant, nodes + list(homozygous)):
            nodes.append(variant)
    parent = list(range(len(nodes)))
    parity = [0] * len(nodes)
    for i, variant in enumerate(nodes):
        if detect_germline_overlap(variant, cis_variants):
            parent[i] = 0

    def find(i):
        if parent[i] != i:
            root, root_parity = find(parent[i])
            parent[i], parity[i] = root, parity[i] ^ root_parity
        return parent[i], parity[i]

    for i in range(len(nodes)):
        for j in range(i + 1, len(nodes)):
            answer = query_in_cis(phase_resolver, nodes[i], nodes[j])
            if answer is None:
                continue
            (root_i, parity_i), (root_j, parity_j) = find(i), find(j)
            relation = 0 if answer else 1
            if root_i != root_j:
                parent[root_i] = root_j
                parity[root_i] = parity_i ^ parity_j ^ relation
            elif parity_i ^ parity_j != relation:
                logger.warning(
                    "Ignoring phase answer for %s and %s that contradicts "
                    "earlier answers", nodes[i], nodes[j])
    somatic_root, somatic_parity = find(0)
    cis, trans, groups = [], [], {}
    for i, g in enumerate(rest, 1):
        root, node_parity = find(i)
        if root != somatic_root:
            groups.setdefault(root, []).append((g, node_parity))
        elif node_parity == somatic_parity:
            cis.append(g)
        else:
            trans.append(g)
    return PhasePartition(
        germline=germline,
        cis=tuple(g for g in germline if g in homozygous or g in cis),
        trans=tuple(trans),
        linked=tuple(
            (tuple(g for g, p in members if p == members[0][1]),
             tuple(g for g, p in members if p != members[0][1]))
            for members in groups.values() if len(members) > 1))

varcode.PhasePartition(germline: Tuple = (), cis: Tuple = (), trans: Tuple = (), linked: Tuple = (), implicit: bool = False) dataclass

Bases: DataclassSerializable

Germline variants near a somatic variant, split by what is known about their phase relative to it.

germline holds every nearby germline variant in input order; cis and trans are those whose phase relative to the somatic variant is known, and :attr:unphased is the rest. linked lists groups of unphased variants whose phase relative to each other is known, as (one_haplotype, other_haplotype) pairs with the group's first variant on one_haplotype; each group changes sides as a unit. implicit is true on hemizygous contigs (chrM, chrY), where every germline variant is cis without evidence. Build one with :func:partition_germline_by_phase.

blocks: Tuple property

Units of unknown phase as (one_haplotype, other_haplotype) pairs: each linked group, then each other unphased variant.

hypothesis_count: int property

One hypothesis per assignment of each block to a haplotype.

hypotheses(max_hypotheses: int = DEFAULT_MAX_PHASE_HYPOTHESES) -> Tuple[PhaseHypothesis, ...]

Enumerate phase hypotheses, keeping known phase fixed in each.

Known phase gives one hypothesis (phase_state "phased", or "implicit" on hemizygous contigs). Otherwise there is one "unknown" hypothesis per assignment of the :attr:blocks. Haplotype labels depend only on the cis set: "A" when every germline variant is cis, "B" when none is, otherwise "A_mixed_<mask>" over :attr:germline.

RAISES	DESCRIPTION
HypothesisLimitError	When :attr:hypothesis_count exceeds max_hypotheses.

Source code in varcode/germline.py

def hypotheses(
        self, max_hypotheses: int = DEFAULT_MAX_PHASE_HYPOTHESES,
) -> Tuple[PhaseHypothesis, ...]:
    """Enumerate phase hypotheses, keeping known phase fixed in each.

    Known phase gives one hypothesis (``phase_state`` ``"phased"``,
    or ``"implicit"`` on hemizygous contigs). Otherwise there is one
    ``"unknown"`` hypothesis per assignment of the :attr:`blocks`.
    Haplotype labels depend only on the cis set: ``"A"`` when every
    germline variant is cis, ``"B"`` when none is, otherwise
    ``"A_mixed_<mask>"`` over :attr:`germline`.

    Raises
    ------
    HypothesisLimitError
        When :attr:`hypothesis_count` exceeds ``max_hypotheses``.
    """
    max_hypotheses = _check_max_hypotheses(max_hypotheses)
    blocks = self.blocks
    if not blocks:
        return (self._hypothesis(
            self.cis, "implicit" if self.implicit else "phased"),)
    if self.hypothesis_count > max_hypotheses:
        raise HypothesisLimitError(max_hypotheses)
    return tuple(
        self._hypothesis(self.cis + tuple(
            v for i, (one, other) in enumerate(blocks)
            for v in (one if mask >> i & 1 else other)), "unknown")
        for mask in range(self.hypothesis_count))

varcode.enumerate_phase_hypotheses(somatic_variant, germline_in_window, phase_resolver=None, max_hypotheses: int = DEFAULT_MAX_PHASE_HYPOTHESES) -> Tuple[PhaseHypothesis, ...]

Enumerate the phase configurations of somatic_variant relative to germline_in_window.

Shorthand for :func:partition_germline_by_phase followed by :meth:PhasePartition.hypotheses. Raises :class:~varcode.HypothesisLimitError when the unphased variants need more than max_hypotheses hypotheses.

Source code in varcode/germline.py

def enumerate_phase_hypotheses(
        somatic_variant,
        germline_in_window,
        phase_resolver=None,
        max_hypotheses: int = DEFAULT_MAX_PHASE_HYPOTHESES,
) -> Tuple[PhaseHypothesis, ...]:
    """Enumerate the phase configurations of ``somatic_variant``
    relative to ``germline_in_window``.

    Shorthand for :func:`partition_germline_by_phase` followed by
    :meth:`PhasePartition.hypotheses`. Raises
    :class:`~varcode.HypothesisLimitError` when the unphased variants
    need more than ``max_hypotheses`` hypotheses.
    """
    return partition_germline_by_phase(
        somatic_variant, germline_in_window, phase_resolver
    ).hypotheses(max_hypotheses)

varcode.detect_loh(somatic_variant, germline_in_window) -> Optional[bool]

Deprecated: these allele-only inputs cannot establish LOH.

Returns None (not assessed), never a positive or negative LOH call. Use :func:detect_germline_overlap for the previous overlap query. LOH requires an independently supported tumor/normal allelic-state assessment; this annotation API does not perform that analysis.

Source code in varcode/germline.py

def detect_loh(somatic_variant, germline_in_window) -> Optional[bool]:
    """Deprecated: these allele-only inputs cannot establish LOH.

    Returns ``None`` (not assessed), never a positive or negative LOH call.
    Use :func:`detect_germline_overlap` for the previous overlap query. LOH
    requires an independently supported tumor/normal allelic-state assessment;
    this annotation API does not perform that analysis.
    """
    warnings.warn(
        "detect_loh cannot infer LOH from allele overlap; returning None. "
        "Use detect_germline_overlap for allele matching.",
        DeprecationWarning, stacklevel=2)
    return None

varcode.detect_germline_overlap(somatic_variant, germline_in_window) -> bool

Whether the same allele occurs in the supplied germline variants.

Match reference assembly, contig, normalized position, REF, and ALT. This says nothing about genotype, allele balance, copy number, or LOH.

Source code in varcode/germline.py

def detect_germline_overlap(somatic_variant, germline_in_window) -> bool:
    """Whether the same allele occurs in the supplied germline variants.

    Match reference assembly, contig, normalized position, REF, and ALT. This
    says nothing about genotype, allele balance, copy number, or LOH.
    """
    for g in germline_in_window:
        if getattr(g, "is_structural", False) or getattr(somatic_variant, "is_structural", False):
            if g == somatic_variant:
                return True
            continue
        if (getattr(g, "reference_name", None) == getattr(somatic_variant, "reference_name", None)
                and g.contig == somatic_variant.contig
                and g.trimmed_base1_start == somatic_variant.trimmed_base1_start
                and g.trimmed_ref == somatic_variant.trimmed_ref
                and g.trimmed_alt == somatic_variant.trimmed_alt):
            return True
    return False

varcode.GermlineAlleleOverlap(variant, transcript)

Bases: TranscriptMutationEffect

This reported allele is already present in the supplied germline.

No new allele sequence is established. This does not assess allele dosage, loss of the other allele, expression, or the allele's pathogenicity.

Source code in varcode/effects/effect_classes.py

def __init__(self, variant, transcript):
    Intragenic.__init__(self, variant, gene=transcript.gene)
    self.transcript = transcript

varcode.default_germline_window(somatic_variant, transcript) -> Tuple[str, int, int]

Default window for looking up germline variants relevant to a somatic variant on a transcript.

Returns (contig, start, end) covering the codon containing the somatic variant — three reference bases on each side of somatic_variant.start. This is the window from #268's table for in-exon coding variants.

Larger windows (splice signal region for splice-adjacent variants, same exon for frameshift candidates) are useful refinements but don't change the API. Callers that need them pass a custom window_fn to :func:predict_germline_aware_effect.

Splice-adjacent: when the somatic is within 6bp of an exon-intron boundary, expand to a 12bp window centered on the boundary so germline edits to the donor / acceptor signal show up in the lookup. This catches the "germline broke the splice site" case without forcing the caller to wire up a separate window function.

Source code in varcode/germline.py

def default_germline_window(somatic_variant, transcript) -> Tuple[str, int, int]:
    """Default window for looking up germline variants relevant to a
    somatic variant on a transcript.

    Returns ``(contig, start, end)`` covering the codon containing the
    somatic variant — three reference bases on each side of
    ``somatic_variant.start``. This is the window from #268's table
    for in-exon coding variants.

    Larger windows (splice signal region for splice-adjacent variants,
    same exon for frameshift candidates) are useful refinements but
    don't change the API. Callers that need them pass a custom
    ``window_fn`` to :func:`predict_germline_aware_effect`.

    Splice-adjacent: when the somatic is within 6bp of an exon-intron
    boundary, expand to a 12bp window centered on the boundary so
    germline edits to the donor / acceptor signal show up in the
    lookup. This catches the "germline broke the splice site" case
    without forcing the caller to wire up a separate window function.
    """
    contig = somatic_variant.contig
    pos = somatic_variant.start
    # Default: the codon containing the somatic. Three bases on either
    # side — slightly wider than strictly necessary so overlapping
    # frame-aware codon membership is conservative.
    start = pos - 3
    end = (getattr(somatic_variant, "end", None) or pos) + 3
    # Splice-adjacent expansion: if any of this transcript's exon
    # boundaries is within 6bp of the somatic, widen to capture the
    # canonical splice signal region (MAG | GURAGU and YAG | R, ~12bp).
    try:
        for exon in transcript.exons:
            for boundary in (exon.start, exon.end):
                if abs(boundary - pos) <= 6:
                    start = min(start, boundary - 6)
                    end = max(end, boundary + 6)
    except Exception:
        # Hand-built transcripts in tests may not have exons; the
        # default codon window is a fine fallback.
        pass
    return contig, max(1, start), end