Phasing¶
Two nearby variants, for example in the same codon, can have a different combined effect than either one alone, but only if they sit on the same copy of the chromosome (in cis). Variants on different copies (in trans) act separately. Phasing evidence tells Varcode which is the case.
Pass phase_resolver= when you have that evidence, from a phased VCF or RNA
reads. Varcode then predicts combined effects for variants in cis. To classify
somatic variants against the patient's own germline sequence, also pass
germline context.
Phased VCF¶
Load the jointly phased calls and select the sample whose phase should be used:
from varcode import load_vcf, VCFPhaseResolver
phased_variants = load_vcf("merged_phased.vcf", genome=81)
phaser = VCFPhaseResolver(phased_variants, sample="TUMOR")
effects = phased_variants.effects(phase_resolver=phaser)
The resolver uses phased GT and PS fields from the loaded collection.
Use one resolver per sample. Phase-set tags are local to their source VCF:
matching tags in independently produced tumor and normal files do not establish
relative phase. A joint call set or molecular evidence is needed.
Individual effects remain alongside joint results. The selected annotator owns
the joint prediction: the default returns HaplotypeEffect for supported point
edits, while annotator="transcript_model" returns classified phase/splice
candidates. Unsupported groups are retained as Unresolved with their
variants and phase_source; they are never silently dropped. The default
declines joint germline composition rather than ignoring the patient baseline.
See the annotator contract.
For germline-aware somatic annotation, the resolver's collection should contain
both sets of alleles, but call somatic_variants.effects(germline=...,
phase_resolver=phaser) on the somatic variants you want to classify.
RNA evidence¶
For direct read/fragment co-occurrence in an RNA-seq BAM:
from varcode import MolecularPhaseResolver, RNAReadPhasingSource
source = RNAReadPhasingSource("tumor.rna.bam")
phaser = MolecularPhaseResolver(source)
effects = variants.effects(phase_resolver=phaser)
Here variants is the loaded collection to annotate. Short- and long-read
data can both leave phase unknown; coverage and linked alleles determine whether
a pair is resolved. RNAReadPhasingSource sorts the fragments covering both
loci into four combinations: both alt alleles, either variant's alt allele
alone, or neither. A combination counts when it has at least min_alt_reads
fragments and is more than reads showing the wrong allele (phasing_error_rate,
1%) would leak into it, by a one-sided binomial test at
max_p_value_for_phasing (0.05). The pair is cis when both alts are seen
together and not each alone, and trans when each is seen alone and not
together; anything else, including all three, is unknown. A variant that arose
later on the other's copy, or next to a germline variant, is always with it,
while the earlier one also appears alone; that is cis, not trans. The rule,
from varcode.rna_read_phasing.four_gamete_phase, matches Isovar 1.39.7. A
source that only reports co-observed
partners establishes cis, and leaves every other pair unknown; sources that see
reference alleles too, such as RNAReadPhasingSource, report trans through
their own in_cis. Raw BAM phasing does not provide an assembled
MutantTranscript. Assembly-backed sources can provide one; see the
source protocols.
ReadPhaseResolver remains a compatibility name for MolecularPhaseResolver.
RNA phase from Isovar¶
Isovar's results can serve as the phasing source. Give run_isovar the matched
germline variants so they are recognized in the assembled RNA:
from isovar import IsovarReadPhasing, run_isovar
from varcode import MolecularPhaseResolver
results = run_isovar(
variants=somatic_variants,
alignment_file="tumor.rna.bam",
germline_variants=germline_ctx.variants,
)
phaser = MolecularPhaseResolver(IsovarReadPhasing(results))
effects = somatic_variants.effects(germline=germline_ctx, phase_resolver=phaser)
With Isovar 1.36 or later, IsovarReadPhasing.in_cis answers from fragments
that cover both variants. Since Isovar 1.39.7 it uses the same four-combination
rule as RNAReadPhasingSource above, with min_shared_fragments_for_phasing
fragments, phasing_error_rate and max_p_value_for_phasing. Earlier versions
called cis or trans by majority, which could call a nested subclone's pair
trans or a read-error pattern cis.
For a matched germline variant, Isovar 1.38 or later reads the germline site in fragments that carry the somatic alt allele: the germline alt allele there means cis, its reference allele trans. Fragments with the somatic reference allele are not counted, since the germline alt allele also comes from normal cells and homozygous sites. When those fragments do not decide, a germline edit in the somatic variant's assembled RNA is cis; when they say trans but the assembly has the edit, the answer is unknown. Isovar 1.36 and 1.37 only give the assembly-based cis.
Any pair Isovar cannot call stays unknown, and Varcode keeps every phase hypothesis for it.
Unknown phase¶
When a somatic and germline variant share a codon and relative phase is unknown,
the result can be a PhaseCandidateSet, with one classified effect per
haplotype hypothesis. Read .candidates to retain the alternatives. Resolved
phase can reduce the set to one effect, as in the example below.
Two variants in one codon¶
CFTR has a somatic at GRCh38 7:117531100 T→A. A neighbouring
germline at 7:117531101 T→C lands in the same codon.
from pyensembl import cached_release
from varcode import Variant, VariantCollection, GermlineContext
g = cached_release(81)
cftr = g.transcript_by_id("ENST00000003084")
somatic = Variant("7", 117_531_100, "T", "A", genome=g)
germline = Variant("7", 117_531_101, "T", "C", genome=g)
ctx = GermlineContext.from_variants([germline], reference_name="GRCh38")
Without germline context¶
eff = somatic.effect_on_transcript(cftr)
print(type(eff).__name__, eff.short_description)
# Substitution p.L159M
With germline context, phase unknown¶
eff = somatic.effect_on_transcript(cftr, germline=ctx)
print(type(eff).__name__, eff.short_description)
# PhaseCandidateSet ?p.L159M
for c in eff.candidates:
ev = c.evidence
print(f" haplotype={ev['haplotype']:<2} "
f"germline_in_cis={[v.short_description for v in ev['germline_variants']]} "
f"=> {c.effect.short_description}")
# haplotype=B germline_in_cis=[] => p.L159M
# haplotype=A germline_in_cis=['chr7 g.117531101T>C'] => p.S159T
The ? prefix on ?p.L159M flags the description as the
most-likely candidate of a PhaseCandidateSet. In application code,
check isinstance(eff, PhaseCandidateSet) and iterate
eff.candidates (a tuple of EffectCandidate objects carrying
per-hypothesis evidence keys); use eff.effects if only the
inner classified MutationEffects are needed.
With known phase¶
A real pipeline gets the cis/trans answer from a phased VCF or an RNA assembly. For this demo, two stub resolvers force each answer so you can see the candidate set collapse to a single effect:
class ForceCis:
source = "demo"
def in_cis(self, v1, v2, transcript=None): return True
class ForceTrans:
source = "demo"
def in_cis(self, v1, v2, transcript=None): return False
eff_cis = somatic.effect_on_transcript(
cftr, germline=ctx, phase_resolver=ForceCis())
print("cis ->", type(eff_cis).__name__, eff_cis.short_description)
# cis -> Substitution p.S159T
eff_trans = somatic.effect_on_transcript(
cftr, germline=ctx, phase_resolver=ForceTrans())
print("trans ->", type(eff_trans).__name__, eff_trans.short_description)
# trans -> Substitution p.L159M
For real evidence, use the phased VCF or RNA
resolvers above. The stubs implement only in_cis(...), which is all the
codon-collapse path consults. Richer resolvers implement more of the
duck-typed interface, including mutant_transcript and phased_partners; the
phasing API documents the built-in resolvers and
source protocols. No public PhaseResolver base class is required.
Combining evidence¶
Given a somatic VariantCollection, patient context, phase resolver, and RNA
resolver:
effects = somatic_variants.effects(
germline=germline_ctx,
phase_resolver=phaser,
rna_resolver=rna,
)
Order: germline modifies the transcript first, phase collapses the
candidate set next, and RNA refines multi-candidate effects. Splice
mechanism sets are reconciled to observed mechanisms; other
multi-outcome effects append observed-only candidates. Cross-axis key is
EffectCandidate.evidence["haplotype"], so an
RNA observation tagged with the same haplotype tag aligns with the
right germline-aware outcome.
Known deletion haplotypes in RNA alignments¶
Splice-aware aligners can encode a known deletion's RNA sequence with N
rather than D, or split the gap among mismatches and smaller deletions.
Opt in to a local sequence hypothesis instead of interpreting the gap as
proof of a DNA deletion:
from varcode import MolecularPhaseResolver, RNAReadPhasingSource
source = RNAReadPhasingSource("tumor.rna.bam", min_alt_reads=1)
source.register_haplotype([known_deletion, adjacent_snv_1, adjacent_snv_2])
resolver = MolecularPhaseResolver(source)
effects = variants.effects(phase_resolver=resolver)
The variants must share one genome dataset and contig. The genome must provide reference sequence across the interval, through a FASTA or annotated transcript. Every retained base between five-base reference flanks must match one alignment and pass the configured quality/read-edge filters. No mates or separate partial haplotypes are stitched together. Registration invalidates cached counts.
Registration does not assert cis phase: it supplies a candidate to test. For registered variants, only full-context matches count as alternate support; nonmatches remain unknown, not trans evidence. The current local mode supports nonoverlapping substitutions and deletions. Register competing combinations separately. Unregistered variants keep the ordinary CIGAR-based behavior.
For the audited GRCh38 MAP2 example, the combination of 209694769 C>A,
209694770 T>G, and the 28-base deletion at 209694773 gives the same local
sequence under 28D, 28N, and 22D2M6D alignments. An unrelated 19,449-base
exon skip lacks the required local anchors and is not supporting evidence.
Implementation limits¶
See germline limitations for the phase-hypothesis cap and normalization requirements, and the phasing API for custom sources.