Usage

Regardless of the data input format, Twitcher’s functionality is essentially the same:

• identify a mutation cluster
• prepare and extract the Query and Reference sequences surroundings the cluster
• re-align the two sequences allowing for template switches
• report the case if the TS-alignment scores better than the simple forward alignment

Depending on the data format, the exact steps required for this vary a bit.

Basic concepts

Twitcher accepts VCF and BAM (or CRAM) data. The two data types have different modes and the type has to be specified in the sub-command: twicher vcf or twitcher reads.

Twitcher expects that the data are from one individual and can extract the target individual from a multi-sample VCF file. The BAM reads have to be long enough to contain full mutations clusters. The program has been tested with PacBio HiFi reads and assembled haplotype data.

Mutation clusters

A TSM is a single mutational event but it typically coded in data as multiple base differences between the original sequence and the mutated copy. Here, we refer to these “multiple base differences” as “mutation clusters”.

Twitcher searches for mutation clusters using a “sliding window” and when triggered, extends and extracts the region fulfilling the cluster criteria. There are two modes for defining the number of differences: number of events (“legacy”), and edit-mass of events (“edit-mass”). While the former works well with individual SNPs, the latter is more appropriate (and recommended) for modern VCF data where variants are clustered into complex MNPs. See Concepts -> Clustering.

The definition of mutation cluster is arbitrary and the challenge of TSM analyses is that any mutation pattern can be explained as a combination of base change and indel mutations. There is (currently) no means for distinguishing the mutational mechanisms from the mutation footprint only.

Creation of alternative haplotype

A read in a BAM file is an observation with the full context. It is fairly straightforward to identify mutation clusters using the CIGAR string and then extract the locus with surrounding context from the read and the reference genome.

For a mutation cluster in a VCF file, Twitcher creates the mutated copy by placing the variant changes on the reference sequence and then adding the fixed-length flanking sequences on both sides of the cluster. It should be noted that the context is copied from the reference, and is guaranteed to correctly represent the study individual’s genome only for the span of “sliding window” up- and downstream of the mutation cluster. See Concepts -> Clustering and Parameters.

The mutation cluster created by a TSM is expected to be inherited as a whole from a parent to the offspring, and base changes belonging to the same cluster be in phase. Bases in a BAM read are naturally in phase, and ideally the variants in VCF data would also be fully phased. Often, that is not the case and VCF data from widely-used variant-calling software contain unphased mutation clusters. Fortunately, TSM analyses don’t require genome-wide phasing and even the short sequencing reads contain enough information to physical phase nearby sites. See Concepts -> Phasing.

To further complicate things, the VCF and BAM files do not always match. For example, GATK Haplotypecaller reassembles messy regions and calls the variants from the constructed haplotypes. The variant coordinates in the VCF are then based on haplotype alignments and may not have matching to the BAM data. Variant phasing error codes are explained in Concepts/Unresolvable clusters.

Alignment model

The tsalign algorithm is generic and, with an adjusted model, can be utilised for studying various types repeat events. The structure of the cost model is explained in Aligner Costs.

By default, Twitcher uses a simplified cost model where only reverse-complement copy events are allowed, copy events have to be between the two different sequences, and all flanking sequences are modelled with the same costs. (A credible TSM should show high similarity across the copy region and, by default, Twitcher models that with a separate model.) If one wants to adjust the cost model, the changes have to be done within a file describing the full model.

The tsalign algorithm allows for multiple switch-and-copy events. The copy events can be “to reference from reference”, “to reference from query”, “to query from query”, and “to query from reference”. If two different copy events have the same score, one of them is arbitrarily (but reproducibly) reported. By default, all events within one mutation cluster have to be of the same type.

It is thought that template switching could happen both “intra-strand”, i.e. replication jumping backwards in the newly copied sequence, or “inter-strand”, i.e. the replication jumping somewhere in the opposite strand. If the sequences were identical before the TSM, intra-strand events cannot be distinguished from inter-strand events. (However, inter-strand events include types of mutations not possible intra-strand. See Löytynoja and Goldman (2017) for details.)

The tsalign algorithm allows both “to-self-from-self” and “to-self-from-other” copy events, but this is meant to model the evolutionary time rather that the intra- and inter-strand process. If the reference and query sequences have diverged, the source sequence at the time of the TSM may resemble more the current “self” that the current “other”, and thus “self-to-self” provides a solution with a better cost. However, for reproducibility of TSM events with an identical mutation pattern, the Twitcher default cost model only allows to-reference-from-query and to-query-from-reference copy events, and the default switch costs are:

Abbreviation	Meaning	Cost
rrf	Ref-from-Ref-Forward	inf
rqf	Ref-from-Qry-Forward	inf
qrf	Qry-from-Ref-Forward	inf
qqf	Qry-from-Qry-Forward	inf
rrr	Ref-from-Ref-Reverse	inf
rqr	Ref-from-Qry-Reverse	2
qrr	Qry-from-Ref-Reverse	2
qqr	Qry-from-Qry-Reverse	inf

It is noteworthy that, e.g. in an analysis of two modern humans, one should not expect the study individual to be the mutated copy. It is possible (and very common) that the reference genome is the derived form and the study individual carries the ancestral form of the sequence.

If one wants to study the derived mutations in the query (e.g. a patient with a disease), one should set rqr_cost = inf and thus only allow copy events to the query from the reference.

Although not part of the cost model, there is a parameter that is relevant here and may affect the output. To speed up the alignment, the tsalign algorithm restricts the switch out and switch in events to the proximity of the mutation cluster. The width of those regions is defined with --range-extension <BASES> and is by default 5 bp. See Parameters -> CLI Reference.

Regions and outputs

Similarly to BCFtools, Twitcher allows defining target regions for the analysis. The regions are inclusive, but the site filtering is done first and mutation clusters are not expanded beyond the region boundaries. To reduce the boundary effect, one should avoid defining many short regions and rather perform the analysis over large regions, and remove unwanted regions in the postprocessing stage.

Twitcher outputs VCF and CSV data, explained in Output and Output -> CSV.

With VCF output, mutation clusters inferred to be caused by TSM events are merged into complex MNVs; variants deemed not to belong to the TSMs are separated and reported independently. Instead of a full, re-annotated VCF file, one may output the variant clusters with an inferred TSM event only with or without the context. With BAM data, only CSV output is provided. The BAM data can be highly redundant and Twitcher outputs each copy of the inferred TSM as an independent row.

For its simplicity and parsing efficiency, the CSV output was chosen as the primary format for downstream processing and visualisation. When analysing VCF data, it is recommended that also the CSV output is specified.

twitcher viz provides an easy-to-use interface to visualisation functionalities of tsalign. The visualisation works via a CSV file and generates SVG graphic files. See Parameters -> CLI-reference.

Twitcher Workbench and other helper tools

In genome-scale analyses, Twitcher should be used with Twitcher Workbench, a graphical interface for post-processing and visualisation of Twitcher output data. The package includes a command-line tool for the post-processing steps and data filtering.

The CSV/TSV toolkit is a command-line tool for the processing of CSV-formatted data. Some of the examples expect it to be installed on the system.

Simple analysis

The VCF file below describes a mutation cluster.

##fileformat=VCFv4.2
##reference=file://hg38.no_alt.fa.gz
##contig=<ID=chr11,length=135086622,md5=1511375dc2dd1b633af8cf439ae90cec>
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  NA24385
chr11   133464064   .   CCC C   .   .   .   GT  0|1
chr11   133464068   .   A   C   .   .   .   GT  0|1
chr11   133464070   .   A   T   .   .   .   GT  0|1
chr11   133464071   .   G   T   .   .   .   GT  0|1
chr11   133464072   .   A   C   .   .   .   GT  0|1

Writing the VCF file as temp.vcf and having a copy of the reference genome, we can run Twitcher with the command:

twitcher vcf --reference hg38.no_alt.fa.gz -o temp.out.vcf.gz --csv temp.out.csv temp.vcf

This produces two output files, one in compressed VCF format (temp.out.vcf.gz) and the other in CSV format (temp.out.csv).

We can visualise the inferred TSM event using the CSV file:

twitcher viz --reference hg38.no_alt.fa.gz -o svgs/ -z 10 -a -e full temp.out.csv

This produces an SVG file:

Twitcher visualisation of a TSM event.

The output reveals that the haplotype carried by the sample individual is, in fact, the ancestral state and the reference is the derived form.

TSIGAR format

The template switch process is described using the TSIGAR format (see Output -> TSIGAR), remotely similar to CIGAR widely used in genome analyses.

To understand the template switch process and the data format, we use csvtk and extract the fields fw_cigar and ts_cigar:

$ csvtk cut -f fw_cigar,ts_cigar temp.out.csv 
fw_cigar,ts_cigar
5=2D1=1X1=3X5=,"3=[TSRQR:[0,2]:[-3,0]:-7:13=:11]2="

fw_cigar is the reconstructed forward alignment for the sequences for the mutation cluster; the amount of sequence flanking the cluster is defined with --range-extension <BASES> and, here, it is 5 bp (5=2D1=1X1=3X5=). The internal part describes the variants from the VCF: a deletion of two, one match, one mismatch, one match, three mismatches.

ts_cigar is the inferred template switch alignment. The external part, 3= and 2=, are the parts of the context not included in the template switch solution, and the part in the square brackets is the actual jump process. TSRQR means that it is a TS-event where the Reference copies from Query in Reverse-complement orientation. The last three values -7:13=:11 describe the jumps and the copy event: the first jump, from A1 to A2, is -7 bases; that is followed by thirteen identities, =13, from A2 to A3; finally, the end point of the return jump, from A3 to A4, is 11 bases from A1.

The values in square brackets describe the uncertainty intervals of the switch positions. In the figure above, these intervals are shown in light and bright blue letters. [0,2] means that the A1->A2 jump could not have occurred any earlier (further left) in the forward part (0), but it could have occurred two bases later (2) than shown in the solution. (The bright blue letters CT could be moved where the light blue letters TC are now.) Similarly, [-3,0] means that the A3->A4jump have occurred three bases earlier (-3) than shown, but not any later (0).

For the tsalign algorithm, any solution within those intervals scores equally well and one is arbitrarily reported. One should note tha the visualisation option used above (-e full) draws the uncertainty regions in both places and some characters are thus duplicated. See Parameters -> CLI-reference.

Complex Twitcher analyses

VCF data

A more comprehensive Twitcher command for CSV data could look like this:

twitcher vcf \
  --reference $ref --phasing-bam $bam $vcf \
  --regions $chr --db known_clusters.$chr.db \
  --output $out.$chr.vcf.gz --csv $out.$chr.csv \
  --stat-output $out.$chr.stats.txt \
  --output-failed-alignments $out.$chr.failed.txt \
  --output-unresolvable $out.$chr.unresolvable.txt \
  --output-alignment-stats $out.$chr.alignment-stats.txt \
  --cluster-strategy edit-mass --only-realigned \
  --memory-limit 4gb --threads 10

Here, the first line defines the reference genome ($ref), the BAM file used for variant phasing ($bam; see Concepts -> Read phasing) and the input VCF file ($vcf). The second line restricts the analysis to a specific chromosome (allowing for parallelisation and more specific outputs) and defines a database file for known clusters. This DB works as a cache and allows for resolving known cases, either in the same individual or in other individuals, without the alignment step. See Concepts -> Database.

The following lines define various kinds of additional outputs: stats is a simple summary, failed and unresolvable list regions that could not analysed e.g. due to out-of-memory errors or VCF being inconsistent (and thus non-phasable) with the given BAM file. alignment-stats gives highly detailed information of each cluster. See Output.

The second last line defines the cluster detection strategy and limits the VCF output to realigned regions only. The last line defines the memory for each thread and the number of threads.

BAM data

A command for PacBio HiFi BAM data could look like this:

twitcher reads \
  --reference $ref $bam \
  --regions $chr --db known_clusters.$chr.db \
  --output $out.$chr.csv \
  --stat-output $out.$chr.stats.txt \
  --output-failed-alignments $out.$chr.failed.txt \
  --cluster-strategy edit-mass \
  --memory-limit 4gb --threads 10

The main difference is that the output can be in CSV format only and each read with a copy of (germline) mutation is reported as an independent line in the CSV file. (This highly redundant CSV output can be deduplicated with Twitcher Workbench.) Even when no DB file is specified, Twitcher uses within-run cache of known clusters, massively speeding up analyses of BAM data.

It is notable that the reproducability of the results and the DB cache depend on the amount of memory provided. The DB records the amount of resources provided and a cluster with OOM error is not retried unless more memory is provided.

Twitcher visualisation of a TSM event.