Read Phasing (for VCF data)

Within a cluster the phase of unphased heterozygous variants (e.g. 0/1) is ambiguous: it is not clear which haplotype each mutation sits on. When that happens, the cluster cannot be realigned and is not further considered. If you have the sequencing reads the variant calls were derived from, Twitcher can try to resolve this ambiguity by local read-based phasing, determining which mutations co-occur on the same reads, which then allows realignments.

Provide an indexed BAM file:

twitcher vcf input.vcf reference.fa --phasing-bam reads.bam -o output.vcf

Logic

1. Phasing runs only for clusters that contain at least one unphased heterozygous record.

2. Also fetched is the reference sequence around the cluster, which is what each record’s alleles are compared against.

3. All reads spanning the cluster are fetched. For each record, Twitcher works out a comparison span — a stretch of reference — and what sequence each of the record’s alleles predicts a read will spell across it. A read is then matched by taking its own sequence over that same span (deleted positions drop out, inserted bases are spliced in) and looking for the allele that predicts exactly that. A read matching no allele is uninformative there; it is never assigned the reference allele by elimination.

   How wide the span is depends on the record:

   • Substitutions (reference and alternate of equal length) claim each base at a fixed position, so the span shrinks to just the bases that actually differ. An SNV’s span is one base wide, which keeps a neighbouring variant’s insertion from leaking into the comparison.
   • Insertions and deletions can be written at more than one reference position and still spell the same haplotype — inside a homopolymer or an STR, the aligner’s placement in a read and the VCF’s left-aligned placement routinely disagree. Their span is therefore widened across the whole repeat the indel could have been placed in, plus a margin base on each side. Any equivalent placement then yields the same sequence, so a shifted indel is recognised instead of being misread as the reference allele. The search for equivalent placements stops after 100 bases in either direction; an indel in a longer repeat cannot be pinned down that way, so such a record contributes no read evidence at all rather than risk a false reference call.
   • Complex alleles (both sides non-empty and of different lengths) are not widened across repeats, but get the same margin base on each side, so bases the alternate allele inserts at either edge of the replaced stretch fall inside the comparison.

   The cost of comparing whole spans is that a read must cover the span and match it exactly: a read stopping inside a long deletion’s span no longer confirms the reference allele there, and neither does one carrying an unrelated variant or a sequencing error inside the span.

4. From these reads it builds a co-occurrence matrix recording, for every pair of variants, how many reads place them on the same haplotype versus on different haplotypes. A read that is reference at both positions counts too: at two heterozygous sites it says the two reference alleles are in cis, hence so are the two alternate alleles, which is exactly as informative as a read carrying both mutations.

5. Variants that are already phased in the input act as anchors: their haplotype is known, so each unphased variant is assigned to a haplotype by majority vote of the reads connecting it to the anchors. When there are no anchors, the variants are instead partitioned greedily into two haplotypes from their co-occurrence alone.

6. A variant is only phased when the read evidence is strong enough (see the parameters below). Variants that stay unresolved remain unphased and continue to block their cluster, which is then output unchanged.

Records phased this way are written with an invented phase set. They carry the TWITCHERPHASE INFO flag, so they can be told apart from phasing inherited from the input VCF (see Genotype changes).

Per-haplotype realignment

Once a multi-record cluster is unambiguous, it is decomposed into per-haplotype sub-clusters that are each realigned on their own. Every sub-cluster has to stand on its own: --cluster-max-gap and the validity gates of the active strategy are applied again to it, and sub-clusters that no longer pass are dropped. Removing the other haplotype’s variants can open a gap or thin out the density, so a split can leave a haplotype with nothing to realign. A phased heterozygous variant contributes to whichever haplotype carries its alternate allele. A homozygous-alternate variant (e.g. 1/1) sits on both haplotypes, so how it is handled depends on the rest of the cluster:

• No phased heterozygous variants: the two haplotypes are identical across the cluster, so each homozygous-alternate stretch is realigned once and reported as homozygous (1/1).
• Phased heterozygous variants present: the haplotypes genuinely differ, so every homozygous-alternate variant is included on each haplotype’s reconstruction. This keeps each haplotype’s realignment input complete (a mutation that is physically present on a haplotype is never omitted) and prevents an unrelated homozygous variant from opening a spurious gap that would split a real sub-cluster.

Phasing Parameters

• --phasing-bam <FILE>: the indexed BAM/CRAM file with the reads. Without it, no local phasing is performed.
• --phasing-min-reads (default: 3): minimum number of reads that must support the winning haplotype assignment. Guards against assignments made on too little coverage.
• --phasing-confidence (default: 0.8): minimum fraction of informative reads that must agree on the assignment, in [0.5, 1.0]. A read is informative for a variant when it covers the position and carries a recognisable allele. Higher values reject conflicting or noisy evidence; lower values phase more aggressively.
• --phasing-min-mapq (default: 20): minimum read mapping quality (MAPQ). Reads below this are skipped entirely, so mis-mapped reads cannot contribute spurious allele calls.
• --phasing-min-base-qual (default: 20): minimum base quality (Phred) at the variant position. A read whose base at the locus is below this is treated as uninformative there, so likely sequencing errors do not become votes. Deletions (which have no read base) are unaffected.

If the anchors within a cluster belong to more than one phase set, the cluster is considered too ambiguous and phasing is skipped.

Unresolvable clusters

A cluster that cannot be decomposed into haplotypes is not realigned; its records are passed through to the output VCF unchanged. To find out which clusters these were, and why, pass --output-unresolvable <FILE>:

twitcher vcf input.vcf reference.fa --phasing-bam reads.bam --output-unresolvable unresolvable.tsv -o output.vcf

The file is a headerless TSV with one line per unresolvable cluster and three columns:

#	Column	Content
1	Region	The cluster’s span as contig:start-end, 1-based and inclusive.
2	Phasing status	What local phasing did with the cluster (see below).
3	Reason	Why the cluster is unresolvable: MoreThanOneUnphasedHet (more than one heterozygous record stayed unphased, so the haplotype assignment is ambiguous) or MissingAllele (a record’s genotype does not name an alternate allele that can be placed on a haplotype).

The phasing status is one of:

• trivial: the cluster needed no phasing (nothing was ambiguous to begin with).
• no reads: no --phasing-bam was given, so no phasing was attempted.
• error(region resolution): no region could be derived from the cluster’s records, so no reads could be fetched.
• error(multiple phasesets), error(bam error), error(other): phasing ran but failed; the parenthesised text names the cause.
• phased[D,D,…]: phasing ran, with one decision D per record it attempted, in cluster order:
  • OK_DIRECT / OK_FLIP: resolved; the record keeps the slot order, or has it swapped.
  • NO_EVIDENCE: no informative read covers the record.
  • MIN_READS: too few reads support the winning assignment (--phasing-min-reads).
  • MIN_CONFIDENCE: the reads disagree too much (--phasing-confidence).
  An empty list (phased[]) means every record was already phased in the input.

Example:

chr1:1000523-1000547    phased[OK_DIRECT,NO_EVIDENCE]   MoreThanOneUnphasedHet
chr7:88214-88231    no reads    MoreThanOneUnphasedHet
chr7:4410902-4410915    trivial MissingAllele