Clustering

VCF

Consecutive variants are grouped into a candidate cluster as long as the gap between them stays within --cluster-max-gap. How a group is then scored and accepted depends on --cluster-strategy (see Clustering strategies below). Under the default legacy strategy a group becomes a cluster when it has at least --cluster-min-records variants (default 2).

Only records that can carry a mutation take part in this: a record is skipped when it does not pass the input FILTER, when it lies fully outside --targets, when its genotype has no alternate allele on either haplotype, or when it has no alternate allele at all. Skipped records are holes, not cluster boundaries: they neither break a candidate cluster nor count towards --cluster-min-records, and the gap to the next variant is measured across them.

Since VCF clusters carry information about more than one haplotype, an alternate allele needs to be chosen when selecting records. This will always be the alternate allele with the highest edit mass. Later, records will be assigned to haplotypes, on which the clustering parameters will be re-evaluated.

BAM

When reading reads from a BAM/SAM/CRAM file, the CIGAR string and reference file are used to find clusters. Each contiguous run of mismatches, insertion, or deletion in the CIGAR is treated as one mutation event, scored exactly like the equivalent VCF record would be: a run of n mismatches counts as an n-base MNV, and an n-base insertion or deletion counts as the corresponding indel. The same clustering parameters and --cluster-strategy choice apply: --cluster-max-gap controls the maximum number of matching bases allowed between consecutive mutation events before the cluster is broken; and --cluster-min-records sets the minimum number of events required. Reference skips, soft clips, hard clips and padding always end a cluster.

Differences between VCF and BAM mode

Mass, span and gap are computed identically in both modes, but two things cannot be made to line up, because a read has no records:

  • Event granularity. A caller may emit three adjacent SNVs where the CIGAR shows a single 3-base mismatch run. The mass is the same either way, but that is three events against one, so --cluster-min-records bites differently. In particular, a read whose only variation is one mismatch block is a single event and can never reach the legacy minimum of two.
  • Allele choice. A read carries one sequence, so there is nothing to choose. For a multi-allelic VCF record one alternate allele has to be picked (see Clustering strategies below).

A cluster’s region, on the other hand, means the same thing in both modes: the reference bases the mutations actually edit. VCF puts an anchor base in front of an indel (A → AT at p records position p, but edits nothing there), and every shared leading base is stripped the same way, so the region starts at p + 1. That is exactly the base a CIGAR D operation covers, and a lone insertion, which no I operation covers at all, gets a zero-width region. This is the region that is realigned, that identifies a record, and that every log line and report names, so the two subcommands describe the same event with the same coordinates, and an alignment cached in --db by one of them can be reused by the other.

Two smaller sources of divergence remain: shared trailing bases are not stripped (AAT → AT gives a two-base region for a one-base deletion, where the CIGAR gives 1D), and --cluster-max-gap is measured from the record’s POS, so a VCF cluster can span one base more than the equivalent read cluster.

Clustering Parameters

  • --cluster-strategy (default: legacy): how candidate clusters are scored and validated; see Clustering strategies.
  • --cluster-max-gap (default: 20): maximum gap in bases between consecutive variants to still be grouped together.
  • --cluster-min-records: minimum number of mutation events required to form a cluster: VCF records carrying an alternate allele, or mutation runs in the CIGAR. Defaults to 2 under legacy and 1 under edit-mass.
  • --cluster-min-mass (default: 2.0): minimum total edit mass required for a cluster. Only valid with --cluster-strategy edit-mass.
  • --cluster-min-density (default: 0, hidden): minimum density score, i.e. the edited fraction of the cluster region. Disabled by default and not recommended: see work item 76. The descriptions below leave it out; when it is set, a cluster must additionally clear it.

Clustering strategies

--cluster-strategy selects between two ways of scoring a candidate cluster. They differ only in the per-event mass, the cluster span, and the validity gates.

Both read a variant as a pair of edited allele lengths: the reference allele and the chosen alternate allele with the bases they share at the front stripped off, giving lengths r and a. The shared bases are not edited, so they carry no mass; the VCF anchor base of an indel is one of them. An SNV is therefore (1, 1), an n-base MNV (n, n), an n-base insertion (A → AT…) (0, n), an n-base deletion (n, 0), and a length-changing MNV such as A → CT is (1, 2). No special case is needed to tell these apart: for an MNV there is nothing to strip, so the trimming is a no-op. For a multi-allelic record the alternate allele is the one local phasing assigned; where there is no such assignment, the alternate allele with the largest mass is used, and the same allele then also determines the span.

The reference extent both strategies build on is the number of reference bases the cluster covers, from the point where the first variant starts editing to the end of the last variant’s reference allele. Shared leading bases shift a variant’s start to the right, so an insertion recorded as A → AT at position p starts editing at p + 1. A single SNV therefore has an extent of 1, and two SNVs 10 bases apart an extent of 10. Matching bases between two variants are inside the extent; matching bases after the last one are not.

legacy (default): the original behaviour, preserved for backwards compatibility.

  • Complexity is the sum of per-variant complexity weights: 1.0 + log₂(max(r, a)). An SNV contributes 1.0, and an indel, MNV or delins editing n bases on its longer side 1.0 + log₂(n).
  • Span is the reference extent of the cluster.
  • A group is a cluster when it has at least --cluster-min-records variants (default 2).

edit-mass: makes single large indels and MNVs easier to target.

  • Edit mass is the sum of per-variant edit lengths, each simply max(r, a). An SNV contributes 1, an n-base MNV or n-base indel n, and a delins such as A → CT contributes 2.
  • Span is the symmetric extent of the cluster: the larger of its reference extent and its query extent. The query extent adjusts the reference extent by the net length change of the variants (insertions lengthen it, deletions shorten it), so an insertion and the mirror-image deletion get the same span. A pure insertion covers no reference base at all, so only this strategy gives it a meaningful span.
  • A group is a cluster when it has at least --cluster-min-records variants (default 1) and a total edit mass of at least --cluster-min-mass. Because the record minimum defaults to 1, a single sufficiently long indel or MNV forms a cluster on its own, while a lone SNV (mass 1) is filtered out by --cluster-min-mass.