Aligner Costs

twitcher bundles two alignment algorithms: 1. tsalign, the A*-Aligner written by Sebastian Schmidt. This is the default and the most capable and customizable aligner. 2. A reimplementation of fpa, the four-point-aligner written originally by Ari Löytynoja. This aligner is faster, but less capable. Use the command line flag --fpa to switch to the four-point aligner. FPA has some limitations: - FPA cannot align indels in the inner part of a template switch. It will only find template switches that consist of only (mis-)matches in the inner sequence. - Because of the algorithm design, it cannot be gap-affine. This means that only the gap-extend costs are used, even when opening a gap. Consider adjusting these to a higher value to compensate for missing the gap-open costs. - FPA only aligns inter-strand template switches, i.e. the four geometries qrr, rqr, qrf and rqf (see below). The intra-strand base costs rrf, rrr, qqf and qqr are ignored, and with them the RRQQOffset table. - FPA does not support flanks, which means that all flank-related values will be ignored. - Furthermore, the following costs that could apply are not considered in the FPA algorithm: Length, LengthDifference, ForwardAntiDescendantGap, ReverseAntiDescendantGap.

Whenever a cost file configures something FPA disregards, `--fpa` emits a warning naming it (see [Warnings about settings](#warnings-about-settings)).

Both aligners can be customized with costs for each edit operation between the two sequences, including jump costs for template switches. When using custom costs, you need to point the tool to a custom .tsa file that contains all definitions. You need to provide the whole file. Providing only a section (e.g. only the # Base Cost section) will lead to an error.

Get the template file from here.

In this document, we will take a look at the sections in the cost file. All values are costs / penalties, i.e. if something has a higher value, the aligner will try to avoid it more.

Limits

The two values describe the length of the flanking regions that surround a template switch. Since the flanks are a bit buggy and not supported, keeping left_flank_length and right_flank_length at value 0 is highly recommended. This also renders the sections # Left Flank Edit Costs and # Right Flank Edit Costs obsolete.

Base Cost

Here, you define a flat cost for each template switch. There are eight values for the eight template switch geometries:

The first letter is r or q, meaning reference or query, respectively. This is the descendant (former: primary) sequence, i.e. the sequence where the points 1 and 4 of the template switch are located.
The second letter is, again, r or q, meaning reference or query, respectively. This is the ancestor (former: secondary) sequence, i.e. the sequence where the points 2 and 3 of the template switch are located. This sequence may or may not be the same as the descendant.
The third letter describes the direction; r for reverse (complement) and f for forwards. This refers to the direction in which the inner template switch sequence is aligned to the template.

So, for a template switch that occurred in the derived sequence and copied the reverse complement of the reference sequence, you’d need to adjust the value qrr_cost.

Jump Costs

This section lets you define finer constraints and preferences for the template switch alignment. Each table is structured as two lines which describe a step function. This is perhaps best described with an example:

ExampleValue
-inf -5 -2 3 6
  10  7  0 4 8

This translates to: - For input values (-inf, -5), the output will be 10 - For [-5, -2), the output is 7 - For [-2, 3), the output is 0 - For [3, 6), the output is 4 - For [6, inf), the output is 8

Or as a little graph:

10 ...-------
 9          |
 8          |                                -------...
 7          ----------                       |
 6                   |                       |
 5                   |                       |
 4                   |              ----------
 3                   |              |
 2                   |              |
 1                   |              |
 0                   ----------------
------------+--------+--------------+--------+---------
   ... -7 -6 -5 -4 -3 -2 -1  0  1  2  3  4  5  6  7 ...

Note the V-shape. This is required for these tables. So the following would be invalid, as it has more of an M-shape:

ExampleValue
-inf -5 -2 3 6
   5  0  3 0 5

Offset

The offset costs price the jump from point 1 to 2. This can be useful if e.g. you want to favour local jumps over ones that reach further up- or downstream. The wider the range of offsets with finite cost, the longer the runtime of the algorithm, as it considers more and further jumps. These costs are split into two cases:

RQQROffset

This defines the jump costs for inter-strand template switches.

RRQQOffset

This defines the jump costs for intra-strand template switches. For example, to exclude forward intra-strand template switches (which are not plausible from a biological point of view), you can set this table to something like:

RRQQOffset
-inf -100   1
 inf    0 inf

Length

This one is self-explanatory: The length of the inner part of the template switch can be taken into account, i.e. the length of the sequence from point 2 to 3.

Length Difference

You can assign a cost to the length difference between the inner template switch sequence and the “gap” in the forward alignment that it leaves. That is, the difference between the distance from point 2 to 3 and the distance from point 1 to 4.

Anti-Descendant Gap

The “Anti-Descendant-Gap” (former: “Anti-Primary-Gap”) is the gap just mentioned. It is the distance of points 1 to 4, covered by the inner alignment in the non-descendant sequence. So, in the earlier “QRR”-template switch example, this describes the length in the reference sequence that is no longer covered by a forward alignment, but instead by an alignment to the inner template switch part. Note that the Anti-Descendant-Gap can be negative, meaning point 4 is left of point 1, implying a repeat of some bases.

ForwardAntiDescendantGap

These costs are used when the template switch copies from a sequence in forward direction.

ReverseAntiDescendantGap

These costs are used when the template switch copies from a sequence’s reverse complement.

Edit costs

These may look more familiar to the trained bioinformatician: These costs simply describe substitution and gap costs for a gap-affine alignment. Each table includes costs for replacing one base with another, and costs for opening and extending a gap. The latter can be specified on a per-base basis.

Primary Edit Costs

The normal, forward aligning costs that apply when outside a template switch.

Secondary Forwards Edit Costs

The edit costs that apply when inside a forwards-copying template switch.

Secondary Reverse Edit Costs

The edit costs that apply when inside a reverse-complement-copying template switch.

Left/Right Flank Edit Costs

When the length of the left or right flank is >0, these costs are used in the flanks around a template switch, left or right, respectively. The idea was that this could be used to encourage a high sequence identity immediately up- and downstream of a template switch. However, it was not fully implemented and, as mentioned above, we recommend ignoring this setting.

Warnings about settings

Some combinations of command line arguments and cost tables are legal but contradictory: they price template switches the aligner can never reach, they enable a search that cannot succeed, or they configure something the selected aligner ignores. On startup, twitcher checks for these and emits a WARN for each. None of them abort the run, since an invocation may take hours, and a heuristic must not be able to kill it. To silence them, lower the log level.

The bundled default_costs.tsa is consistent with the default --padding 200, so a plain invocation emits none of these.

--min-length-strategy, --uncertainty-range-mode and --allow-mixed-descendants only configure the template switch search of the A* aligner. Combining them with --no-ts or --fpa is accepted rather than rejected, so that one set of flags can be reused across a batch of runs, but it warns. Passing the default value explicitly changes nothing and is not reported.

Warning	Trigger	Resolution
`RQQROffset` / `RRQQOffset` prices unreachable jumps	The offset table is finite beyond `±padding`. A jump can only land inside the padded window, so the rest of the table is dead weight.	Raise `--padding` to the table’s bound, or narrow the table.
`Length` permits switches longer than the available template	The last `Length` breakpoint exceeds `2 × padding`.	Lower the breakpoint, or raise `--padding`. An unbounded `Length` table is taken as deliberate and not reported.
No template switch can ever be found	Every base cost is infinite, the `Length` table is infinite everywhere, or every reachable offset table is.	Fix the cost file, or pass `--no-ts`, which skips the template switch machinery and is faster.
`preprocess-filter` only prunes reverse template switches	A forward geometry (`rrf`, `rqf`, `qrf`, `qqf`) has a finite base cost while the minimum length strategy is `preprocess-filter`, which is implemented for reverse switches only and therefore prunes nothing.	`--min-length-strategy lookahead`.
A template switch flag has no effect	`--no-ts` or `--fpa` is combined with a non-default `--min-length-strategy`, `--uncertainty-range-mode` or `--allow-mixed-descendants`.	Drop either the mode or the other flag. Harmless if the flags are shared across a batch of runs.
`--fpa` ignores parts of the cost model	The cost file configures something the four-point aligner disregards (see the limitations at the top of this document).	Use the A* aligner, or accept that those tables have no effect.
Flanks are not supported	`left_flank_length` or `right_flank_length` is non-zero. Besides being incomplete, this disables extending an alignment beyond its range at equal cost.	Set both to `0`.
`--cluster-min-density` can never reject a cluster	The threshold is at or below `1 / (--cluster-max-gap + 1)`, which is the lowest density any cluster can have.	Raise the density, or widen `--cluster-max-gap`.
`--cluster-min-density` above 1 under `legacy`	A run of adjacent SNVs has a density of exactly `1.0`; only insertions, which cover next to no reference bases, get above it. So nearly nothing is clustered.	Lower the density, or use `--cluster-strategy edit-mass`.
`--cluster-min-records 1` under `legacy`	Every single variant becomes a cluster and is realigned.	Raise the minimum, or use `--cluster-strategy edit-mass`, where a one-record minimum is the intended default.