Skip to content

Getting started with tomahawk

About

This is an introductory tutorial for using Tomahawk. It will cover:

  • Importing data into Tomahawk
  • Computing linkage-disequilibrium
  • Subsetting and filtering output data
  • Aggregating datasets

Usage instructions

The CLI of Tomahawk comprises of several distinct subroutines (listed below). Executing tomahawk gives a list of commands with brief descriptions and tomahawk <command> gives detailed details for that command.

All primary Tomahawk commands operate on the binary Tomahawk twk and Tomahawk output two file format. Interconversions between twk and vcf/bcf is supported through the commands import for vcf/bcf->twk and view for twk->vcf. Linkage disequilibrium data is written out in binary two format or human-readable ld format.

Command Description
aggregate data rasterization framework for TWO files
calc calculate linkage disequilibrium
scalc calculate linkage disequilibrium for a single site
concat concatenate TWO files from the same set of samples
import import VCF/VCF.gz/BCF to TWK
sort sort TWO file
view TWO->LD/TWO view, TWO subset and filter
haplotype extract per-sample haplotype strings in FASTA/binary format
relationship compute marker-based pair-wise sample relationship matrices
decay compute LD-decay over distance
prune perform graph-based LD-pruning of variant sites

Importing into Tomahawk

By design Tomahawk only operates on diploid and bi-allelic SNVs and as such filters out indels and complex variants. Tomahawk does not support mixed phasing of genotypes in the same variant (e.g. 0|0, 0/1). If mixed phasing is found for a record, all genotypes for that site are converted to unphased genotypes. This is a conscious design choice as this will internally invoke the correct algorithm to use for mixed-phase cases.

Importing standard files to Tomahawk involes using the import command. The following command imports a bcf file and outputs file.twk while filtering out variants with >20% missingness and sites that deviate from Hardy-Weinberg equilibrium with a probability < 0.001.

1
tomahawk import -i file.bcf -o file -m 0.2 -H 1e-3
1
2
$ md5sum 1kgp3_chr6.bcf
8c05554ebe1a51e99be2471c96fad4d9  1kgp3_chr6.bcf
1
2
$ bcftools view 1kgp3_chr6.bcf -HG | wc -l
5024119
1
$ time tomahawk import -i 1kgp3_chr6.bcf -o 1kgp3_chr6

Auto-completion of file extensions

You do not have to append the .twk suffix to the output name as Tomahawk will automatically add this if missing. In the example above "1kgp3_chr6" will be converted to "1kgp3_chr6.twk" automatically. This is true for most Tomahawk commands when using the CLI.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
Program:   tomahawk-264d039a-dirty (Tools for computing, querying and storing LD data)
Libraries: tomahawk-0.7.0; ZSTD-1.3.8; htslib 1.9
Contact: Marcus D. R. Klarqvist <mk819@cam.ac.uk>
Documentation: https://github.com/mklarqvist/tomahawk
License: MIT
----------
[2019-01-21 11:58:15,692][LOG] Calling import...
[2019-01-21 11:58:15,692][LOG][READER] Opening 1kgp3_chr6.bcf...
[2019-01-21 11:58:15,695][LOG][VCF] Constructing lookup table for 86 contigs...
[2019-01-21 11:58:15,695][LOG][VCF] Samples: 2,504...
[2019-01-21 11:58:15,695][LOG][WRITER] Opening 1kgp3_chr6.twk...
[2019-01-21 11:58:39,307][LOG] Duplicate site dropped: 6:18233985
[2019-01-21 11:59:26,749][LOG] Duplicate site dropped: 6:55137646
[2019-01-21 11:59:39,315][LOG] Duplicate site dropped: 6:67839893
[2019-01-21 11:59:47,176][LOG] Duplicate site dropped: 6:74373442
[2019-01-21 11:59:51,440][LOG] Duplicate site dropped: 6:77843171
[2019-01-21 12:00:41,784][LOG] Duplicate site dropped: 6:121316830
[2019-01-21 12:01:12,830][LOG] Duplicate site dropped: 6:148573620
[2019-01-21 12:01:16,557][LOG] Duplicate site dropped: 6:151397786
[2019-01-21 12:01:42,644][LOG] Wrote: 4,784,608 variants to 9,570 blocks...
[2019-01-21 12:01:42,644][LOG] Finished: 03m26,952s
[2019-01-21 12:01:42,644][LOG] Filtered out 239,511 sites (4.76722%):
[2019-01-21 12:01:42,645][LOG]    Invariant: 15,485 (0.308213%)
[2019-01-21 12:01:42,645][LOG]    Missing threshold: 0 (0%)
[2019-01-21 12:01:42,645][LOG]    Insufficient samples: 0 (0%)
[2019-01-21 12:01:42,645][LOG]    Mixed ploidy: 0 (0%)
[2019-01-21 12:01:42,645][LOG]    No genotypes: 0 (0%)
[2019-01-21 12:01:42,645][LOG]    No FORMAT: 0 (0%)
[2019-01-21 12:01:42,645][LOG]    Not biallelic: 26,277 (0.523017%)
[2019-01-21 12:01:42,645][LOG]    Not SNP: 194,566 (3.87264%)

Computing linkage-disequilibrium

All pairwise comparisons

In this first example, we will compute all-vs-all LD associations using the data we imported in the previous section. To limit compute, we restrict our attention to a 1/45 section of the data by passing the -c and -C job parameters. We will be using 8 threads (-t), but you may need to modify this to match the hardware available on your host machine. This job involves comparing a pair of variants >141 billion times and as such takes around 30 min to finish on most machines.

1
$ tomahawk calc -pi 1kgp3_chr6.twk -o 1kgp3_chr6_1_45 -C 1 -c 45 -t 8

Valid job balancing partitions

When computing genome-wide LD the balancing requires that number of sub-problems (-c) must be a member of the function: $$ \binom{c}{2} + c = \frac{c^2 + c}{2}, c > 0 $$ This function is equivalent to the upper-triangular of a square (c-by-c) matrix plus the diagonal. Read more about load partitioning in Tomahawk.

Here is the first 100 valid partition sizes:
1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 66, 78, 91, 105, 120, 136, 153, 171, 190, 210, 231, 253, 276, 300, 325, 351, 378, 406, 435, 465, 496, 528, 561, 595, 630, 666, 703, 741, 780, 820, 861, 903, 946, 990, 1035, 1081, 1128, 1176, 1225, 1275, 1326, 1378, 1431, 1485, 1540, 1596, 1653, 1711, 1770, 1830, 1891, 1953, 2016, 2080, 2145, 2211, 2278, 2346, 2415, 2485, 2556, 2628, 2701, 2775, 2850, 2926, 3003, 3081, 3160, 3240, 3321, 3403, 3486, 3570, 3655, 3741, 3828, 3916, 4005, 4095, 4186, 4278, 4371, 4465, 4560, 4656, 4753, 4851, 4950, 5050

Massive workload

By default, Tomahawk will attempt to compute the genome-wide linkage-disequilibrium for all the provided variants and samples. This generally require huge amount of pairwise variant comparisons. Depending on your interests, you may not want to undertake this large compute and need to parameterize for your particular use-case.

Different load balancing

As a consequence of the job balancing approach in Tomahawk, the diagonal jobs will always involve less variant comparisons (\( \binom{n}{2} \)) compared to off-diagonal jobs (\(n \times m\)), where \(n\) and \(m\) are the number of records in either block. Keep this in mind if you are choosing a partition size given the run-times of a given job. For example, examine the total run-time of the second job (-C 2) — the first off-diagonal job — as a proxy for the expected average run-time per job.

Output ordering

By design, Tomahawk computes pairwise associations using an out-of-order execution paradigm and further permutes the input ordering by using cache-blocking. Additionally, each thread/slave transiently store a number of computed associations in a private buffer that will be flushed upon reaching some frequency threshold. Because of these technicalities, Tomahawk will neither produce ordered output nor will the permutation of the output order be identical between runs on identical data. This has no practical consequence to most downstream applications in Tomahawk with the exception of query speeds. To address this limitation we have introduced a powerful sorting paradigm that will be discussed below.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
Program:   tomahawk-264d039a-dirty (Tools for computing, querying and storing LD data)
Libraries: tomahawk-0.7.0; ZSTD-1.3.8; htslib 1.9
Contact: Marcus D. R. Klarqvist <mk819@cam.ac.uk>
Documentation: https://github.com/mklarqvist/tomahawk
License: MIT
----------
[2019-01-21 12:13:50,210][LOG] Calling calc...
[2019-01-21 12:13:50,211][LOG][READER] Opening 1kgp3_chr6.twk...
[2019-01-21 12:13:50,212][LOG] Samples: 2,504...
[2019-01-21 12:13:50,212][LOG][BALANCING] Using ranges [0-1063,0-1063] in square mode...
[2019-01-21 12:13:50,212][LOG] Allocating 1,063 blocks...
[2019-01-21 12:13:50,213][LOG] Running in standard mode. Pre-computing data...
[2019-01-21 12:13:50,213][LOG][SIMD] Vectorized instructions available: SSE4...
[2019-01-21 12:13:50,213][LOG] Constructing list, vector, RLE...
[2019-01-21 12:13:50,213][LOG][THREAD] Unpacking using 7 threads: ....... Done! 01,291s
[2019-01-21 12:13:51,504][LOG] 531,500 variants from 1,063 blocks...
[2019-01-21 12:13:51,504][LOG][PARAMS] square=TRUE,window=FALSE,low_memory=FALSE,bitmaps=FALSE,single=FALSE,force_phased=TRUE,force_unphased=FALSE,compression_level=1,block_size=500,output_block_size=10000,l_surrounding=500000,minP=1.000000,minR2=0.100000,maxR2=100.000000,minDprime=0.000000,maxDprime=100.000000,n_chunks=45,c_chunk=0,n_threads=8,ldd_type=3,cycle_threshold=0
[2019-01-21 12:13:51,504][LOG] Performing: 141,245,859,250 variant comparisons...
[2019-01-21 12:13:51,504][LOG][WRITER] Opening 1kgp3_chr6_1_45.two...
[2019-01-21 12:13:51,505][LOG][THREAD] Spawning 8 threads: ........
[2019-01-21 12:13:51,533][PROGRESS] Time elapsed       Variants           Genotypes         Output  Progress    Est. Time left
[2019-01-21 12:14:21,533][PROGRESS]      30,000s  3,685,746,500   9,229,109,236,000        988,858   2.60945%   18m39s
[2019-01-21 12:14:51,533][PROGRESS]   01m00,000s  6,425,868,750  16,090,375,350,000      1,689,422   4.54942%   20m58s
[truncated]
[2019-01-21 12:39:51,546][PROGRESS]   26m00,013s   140,131,382,750 350,888,982,406,000     48,635,114    99.211%    12s
[2019-01-21 12:40:04,317][PROGRESS] Finished in 26m12,784s. Variants: 141,245,859,250, genotypes: 353,679,631,562,000, output: 49,870,388
[2019-01-21 12:40:04,317][PROGRESS] 89,806,242 variants/s and 224,874,830,855 genotypes/s
[2019-01-21 12:40:04,321][LOG][PROGRESS] All done...26m12,815s!

This run generated 935817820 bytes (935.8 MB) of output data in binary format using the default compression level (1):

1
2
$ ls -l 1kgp3_chr6_1_45.two
-rw-rw-r-- 1 mk819 mk819 935817820 Jan 21 12:40 1kgp3_chr6_1_45.two

If this data was stored in plain, human-readable, text format it would use over 45 GB:

1
2
$ tomahawk view -i 1kgp3_chr6_1_45.two | wc -c
4544702179

Sliding window

If you are working with a species with well-know LD structure (such as humans) you can reduce the computational cost by limiting the search-space to a fixed-sized sliding window (-w). In window mode you are free to choose any arbitrary sub-problem (-c) size.

Window properties

By default, Tomahawk computes the linkage disequilibrium associations for all pairs of variants within -w bases of the target SNV on either direction. In order to maximize computational throughput, we made the design decision to maintain blocks of data that overlap to the target interval. This has the consequence that adjacent data that may not directly overlap with the target region will still be computed. The technical reasoning for this is, in simplified terms, that the online repackaging of internal data blocks is more expensive than computing a small (relatively) number of non-desired (off-target) associations.

1
time tomahawk calc -pi 1kgp3_chr6.twk -o 1kgp3_chr6_4mb -w 4000000
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
Program:   tomahawk-264d039a-dirty (Tools for computing, querying and storing LD data)
Libraries: tomahawk-0.7.0; ZSTD-1.3.8; htslib 1.9
Contact: Marcus D. R. Klarqvist <mk819@cam.ac.uk>
Documentation: https://github.com/mklarqvist/tomahawk
License: MIT
----------
[2019-01-21 14:10:33,616][LOG] Calling calc...
[2019-01-21 14:10:33,616][LOG][READER] Opening 1kgp3_chr6.twk...
[2019-01-21 14:10:33,619][LOG] Samples: 2,504...
[2019-01-21 14:10:33,619][LOG][BALANCING] Using ranges [0-9570,0-9570] in window mode...
[2019-01-21 14:10:33,619][LOG] Allocating 9,570 blocks...
[2019-01-21 14:10:33,620][LOG] Running in standard mode. Pre-computing data...
[2019-01-21 14:10:33,620][LOG][SIMD] Vectorized instructions available: SSE4...
[2019-01-21 14:10:33,620][LOG] Constructing list, vector, RLE...
[2019-01-21 14:10:33,620][LOG][THREAD] Unpacking using 7 threads: ....... Done! 15,774s
[2019-01-21 14:10:49,394][LOG] 4,784,608 variants from 9,570 blocks...
[2019-01-21 14:10:49,394][LOG][PARAMS] square=TRUE,window=TRUE,low_memory=FALSE,bitmaps=FALSE,single=FALSE,force_phased=TRUE,force_unphased=FALSE,compression_level=1,block_size=500,output_block_size=10000,window_size=4000000,l_surrounding=500000,minP=1.000000,minR2=0.100000,maxR2=100.000000,minDprime=0.000000,maxDprime=100.000000,n_chunks=1,c_chunk=0,n_threads=8,ldd_type=3,cycle_threshold=0
[2019-01-21 14:10:49,394][LOG] Performing: 11,446,234,464,528 variant comparisons...
[2019-01-21 14:10:49,394][LOG][WRITER] Opening 1kgp3_chr6_4mb.two...
[2019-01-21 14:10:49,402][LOG][THREAD] Spawning 8 threads: ........
[2019-01-21 14:10:49,424][PROGRESS] Time elapsed       Variants           Genotypes         Output  Progress    Est. Time left
[2019-01-21 14:11:19,424][PROGRESS]      30,000s  2,311,240,500   5,787,346,212,000        853,368         0    0
[2019-01-21 14:11:49,424][PROGRESS]   01m00,000s  4,562,731,500  11,425,079,676,000      1,717,956         0    0
[truncated]
[2019-01-21 16:12:49,492][PROGRESS] 02h02m00,068s   527,889,726,500  1,321,835,875,156,000    470,354,846         0 0
[2019-01-21 16:13:17,236][PROGRESS] Finished in 02h02m27,812s. Variants: 529,807,522,528, genotypes: 1,326,638,036,410,112, output: 473,514,018
[2019-01-21 16:13:17,237][PROGRESS] 72,104,114 variants/s and 180,548,701,880 genotypes/s
[2019-01-21 16:13:17,257][LOG][PROGRESS] All done...02h02m27,854s!

Single interval vs its local neighbourhood

In many cases we are not interested in computing large-scale linkage-disequilibrium associations but have a limited genomics window of interest. If this is the case, you can use the specialization subroutine of calc called scalc that is parameterizable with a target region and a neighbourhood size in bases.

In this example, we will compute the regional associations of variants mapping to the interval chr6:10e6-10.1e6 and a 100kb neighbourhood.

1
tomahawk scalc -i 1kgp3_chr6.twk -o test -I 6:10e6-10.1e6 -w 100000

Neighbourhood interval

The neighbourhood is computed as the region from [start of interval - neighbourhood, start of interval) and (end of interval + neighbourhood). If the start of the interval minus the neighbourhood falls outside the chromosome (position < 0) then this interval is truncated to 0 in the left end.

Scalability

This subroutine was designed to be used with a relatively short interval to maximize computability in this case. It is however possible to provide any valid interval as the target region to compute linkage-disequilibrium for. Providing large intervals will result in poor performance and potentially undesired output.

Concatenating multiple archives

On of the immediate downsides of partitioning compute into multiple non-overlapping sub-problems is that we will generate a large number of independent files that must be merged prior to downstream analysis. Fortunatly, this is a trivial operation in Tomahawk and involves the concat command.

First lets compute some data using the first 3 out of 990 partitions of the dataset used above.

1
for i in {1..3}; do time tomahawk calc -pi 1kgp3_chr6.twk -c 990 -C $i -o part$i\_3.two; done

Next, since we only have three files, we can concatenate (merge) these files together into a single archieve by passing each file name to the command:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
$ time tomahawk concat -i part1_3.two -i part2_3.two -i part3_3.two -o part3_concat
[2019-01-22 10:55:47,799][LOG] All files are compatible. Beginning merging...
[2019-01-22 10:55:47,800][LOG] Appending part1_3.two... 397.953344 Mb/98.987848 Mb
[2019-01-22 10:55:47,889][LOG] Appending part2_3.two... 359.538884 Mb/54.763322 Mb
[2019-01-22 10:55:47,938][LOG] Appending part3_3.two... 346.075500 Mb/52.072741 Mb
[2019-01-22 10:55:47,988][LOG] Finished. Added 3 files...
[2019-01-22 10:55:47,988][LOG] Total size: Uncompressed = 1.103568 Gb and compressed = 205.823911 Mb

real    0m0.226s
user    0m0.036s
sys     0m0.190s

Passing every single input file in the command line does not become feasble when we have many hundreds to thousands of files. To address this, it is possible to first store the target file paths in a text file and then pass that to Tomahawk. Using the data from above, we save these file paths to the file part_file_list.txt.

1
for i in {1..3}; do echo part$i\_3.two >> part_file_list.txt; done

Then we simply pass this list of file paths to Tomahawk using the -I argument:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
$ time tomahawk concat -I part_file_list.txt -o part3_concat 
adding file=part1_3.two
adding file=part2_3.two
adding file=part3_3.two
[2019-01-22 10:57:43,227][LOG] All files are compatible. Beginning merging...
[2019-01-22 10:57:43,262][LOG] Appending part1_3.two... 397.953344 Mb/98.987848 Mb
[2019-01-22 10:57:43,350][LOG] Appending part2_3.two... 359.538884 Mb/54.763322 Mb
[2019-01-22 10:57:43,400][LOG] Appending part3_3.two... 346.075500 Mb/52.072741 Mb
[2019-01-22 10:57:43,447][LOG] Finished. Added 3 files...
[2019-01-22 10:57:43,447][LOG] Total size: Uncompressed = 1.103568 Gb and compressed = 205.823911 Mb

real    0m0.477s
user    0m0.040s
sys     0m0.241s

Possible duplications

The concat subroutine will dedupe input file paths but will neiter check nor guarantee that the actual input data is not duplicated. Therefore, it is possible to get duplicates in your data if you are not careful. These duplicate entries could corrupt any downstream insights!

Sorting output files

All subroutines in Tomahawk will work on unsorted output files. However, sorted files are much faster to query and will result in considerable time savings. This is especially true for very large files. Sorting files in Tomahawk is trivial and involves calling the sort command with the required fields input (-i) and output (-o).

Disk usage

Tomahawk can easily sort huge .two files with the help of external merge sorting algorithms followed by a single k-way merge, memory assisted operation. The external merge step require additional disk space approximately equal to the size of the input dataset. Then, the merge operation will write a, generally, smaller output file. On average, we can approxiate that the sort operation require an additional two times the input file size on disk. After the sorting procedure is completed, the intermediate files are removed.

Temporary files

By default, Tomahawk will generate temporary files in the same directory as the input file. If this file path is undesired then you have to modify the appropriate parameters.

In this example we will sort the almost 500 million (473,514,826) records computed from the sliding window example above. This archive corresponds to >50 Gb of binary data.

1
$ time tomahawk sort -i 1kgp3_chr6_4mb.two -o 1kgp3_chr6_4mb_sorted
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
Program:   tomahawk-264d039a-dirty (Tools for computing, querying and storing LD data)
Libraries: tomahawk-0.7.0; ZSTD-1.3.8; htslib 1.9
Contact: Marcus D. R. Klarqvist <mk819@cam.ac.uk>
Documentation: https://github.com/mklarqvist/tomahawk
License: MIT
----------
[2019-01-22 11:20:43,978][LOG] Calling sort...
[2019-01-22 11:20:43,993][LOG] Blocks: 47,358
[2019-01-22 11:20:43,993][LOG] Uncompressed size: 50.192950 Gb
[2019-01-22 11:20:43,993][LOG] Sorting 473,514,826 records...
[2019-01-22 11:20:43,993][LOG][THREAD] Data/thread: 6.274119 Gb
[2019-01-22 11:20:43,993][PROGRESS] Time elapsed       Variants  Progress   Est. Time left
[2019-01-22 11:20:43,994][LOG][THREAD] Slave-0: range=0->5919/47358 and name 1kgp3_chr6_4mb_sorted_fileSucRkC.two
[2019-01-22 11:20:43,994][LOG][THREAD] Slave-1: range=5919->11838/47358 and name 1kgp3_chr6_4mb_sorted_filewPI10T.two
[2019-01-22 11:20:43,994][LOG][THREAD] Slave-2: range=11838->17757/47358 and name 1kgp3_chr6_4mb_sorted_fileo1zcHb.two
[2019-01-22 11:20:43,994][LOG][THREAD] Slave-3: range=17757->23676/47358 and name 1kgp3_chr6_4mb_sorted_fileBvRnnt.two
[2019-01-22 11:20:43,994][LOG][THREAD] Slave-4: range=23676->29595/47358 and name 1kgp3_chr6_4mb_sorted_filewbBz3K.two
[2019-01-22 11:20:43,994][LOG][THREAD] Slave-5: range=29595->35514/47358 and name 1kgp3_chr6_4mb_sorted_fileCrNLJ2.two
[2019-01-22 11:20:43,994][LOG][THREAD] Slave-6: range=35514->41433/47358 and name 1kgp3_chr6_4mb_sorted_filerCfYpk.two
[2019-01-22 11:20:43,994][LOG][THREAD] Slave-7: range=41433->47358/47358 and name 1kgp3_chr6_4mb_sorted_filexx9a6B.two
[2019-01-22 11:21:13,994][PROGRESS]      30,000s     70,594,275   14.9086%  2m51s
[2019-01-22 11:21:43,994][PROGRESS]   01m00,000s    146,732,835    30.988%  2m13s
[2019-01-22 11:22:13,994][PROGRESS]   01m30,000s    245,535,675   51.8539%  1m23s
[2019-01-22 11:22:43,995][PROGRESS]   02m00,001s    323,064,235   68.2268%  55s
[2019-01-22 11:23:13,995][PROGRESS]   02m30,001s    421,592,790   89.0348%  18s
[2019-01-22 11:23:24,951][LOG][THREAD] Finished: 1kgp3_chr6_4mb_sorted_filewbBz3K.two with 23676-29595. Sorted n=59190000 variants with size=1.176156 Gb
[2019-01-22 11:23:25,114][LOG][THREAD] Finished: 1kgp3_chr6_4mb_sorted_fileBvRnnt.two with 17757-23676. Sorted n=59190000 variants with size=1.286938 Gb
[2019-01-22 11:23:25,385][LOG][THREAD] Finished: 1kgp3_chr6_4mb_sorted_filerCfYpk.two with 35514-41433. Sorted n=59190000 variants with size=1.219126 Gb
[2019-01-22 11:23:25,710][LOG][THREAD] Finished: 1kgp3_chr6_4mb_sorted_fileCrNLJ2.two with 29595-35514. Sorted n=59190000 variants with size=1.161081 Gb
[2019-01-22 11:23:26,384][LOG][THREAD] Finished: 1kgp3_chr6_4mb_sorted_filexx9a6B.two with 41433-47358. Sorted n=59184826 variants with size=1.215231 Gb
[2019-01-22 11:23:26,425][LOG][THREAD] Finished: 1kgp3_chr6_4mb_sorted_fileSucRkC.two with 0-5919. Sorted n=59190000 variants with size=1.242923 Gb
[2019-01-22 11:23:29,966][LOG][THREAD] Finished: 1kgp3_chr6_4mb_sorted_fileo1zcHb.two with 11838-17757. Sorted n=59190000 variants with size=1.686758 Gb
[2019-01-22 11:23:31,300][LOG][THREAD] Finished: 1kgp3_chr6_4mb_sorted_filewPI10T.two with 5919-11838. Sorted n=59190000 variants with size=2.013554 Gb
[2019-01-22 11:23:31,338][PROGRESS] 02m47,344s  473,514,826 (2,829,587 variants/s)
[2019-01-22 11:23:31,338][PROGRESS] Finished!
[2019-01-22 11:23:31,341][LOG] Spawning 112 queues with 2.380952 Mb each...
[2019-01-22 11:23:35,090][LOG][WRITER] Opening "1kgp3_chr6_4mb_sorted.two"...
[2019-01-22 11:23:35,091][PROGRESS] Time elapsed       Variants  Progress   Est. Time left
[2019-01-22 11:24:05,091][PROGRESS]      30,000s     39,757,023   8.39615%  5m27s
[2019-01-22 11:24:35,091][PROGRESS]   01m00,000s     76,990,000   16.2593%  5m9s
[2019-01-22 11:25:05,091][PROGRESS]   01m30,000s    107,840,000   22.7744%  5m5s
[2019-01-22 11:25:35,091][PROGRESS]   02m00,000s    133,420,000   28.1765%  5m5s
[2019-01-22 11:26:05,092][PROGRESS]   02m30,000s    170,215,235   35.9472%  4m27s
[2019-01-22 11:26:35,092][PROGRESS]   03m00,001s    208,070,000   43.9416%  3m49s
[2019-01-22 11:27:05,092][PROGRESS]   03m30,001s    247,010,000   52.1652%  3m12s
[2019-01-22 11:27:35,092][PROGRESS]   04m00,001s    287,690,000   60.7563%  2m35s
[2019-01-22 11:28:05,092][PROGRESS]   04m30,001s    325,480,000    68.737%  2m2s
[2019-01-22 11:28:35,093][PROGRESS]   05m00,001s    364,299,966   76.9353%  1m29s
[2019-01-22 11:29:05,093][PROGRESS]   05m30,002s    401,367,427   84.7634%  59s
[2019-01-22 11:29:35,093][PROGRESS]   06m00,002s    439,290,000   92.7722%  28s
[2019-01-22 11:30:02,180][PROGRESS] 06m27,089s  473,514,826 (1,223,269 variants/s)
[2019-01-22 11:30:02,180][PROGRESS] Finished!
[2019-01-22 11:30:02,194][LOG] Finished merging! Time: 06m27,103s
[2019-01-22 11:30:02,194][LOG] Deleting temp files...
[2019-01-22 11:30:02,194][LOG] Deleted 1kgp3_chr6_4mb_sorted_fileSucRkC.two
[2019-01-22 11:30:02,194][LOG] Deleted 1kgp3_chr6_4mb_sorted_filewPI10T.two
[2019-01-22 11:30:02,194][LOG] Deleted 1kgp3_chr6_4mb_sorted_fileo1zcHb.two
[2019-01-22 11:30:02,194][LOG] Deleted 1kgp3_chr6_4mb_sorted_fileBvRnnt.two
[2019-01-22 11:30:02,194][LOG] Deleted 1kgp3_chr6_4mb_sorted_filewbBz3K.two
[2019-01-22 11:30:02,194][LOG] Deleted 1kgp3_chr6_4mb_sorted_fileCrNLJ2.two
[2019-01-22 11:30:02,194][LOG] Deleted 1kgp3_chr6_4mb_sorted_filerCfYpk.two
[2019-01-22 11:30:02,194][LOG] Deleted 1kgp3_chr6_4mb_sorted_filexx9a6B.two
[2019-01-22 11:30:02,678][LOG] Finished!

real    9m18.793s
user    23m31.390s
sys     0m53.197s

Faster queries

Success: You can now use the newly created sorted archieve exactly as unsorted files but with faster query times.

Format descriptions

Memory layout (technical)

A frequent design choice when designing a file format is how to align memory of records: either as lists of records (array-of-struct) versus column-orientated layouts (struct-of-array). The pivoted layout (struct-of-array) of twk1_two_t structs in Tomahawk will result considerable savings in disk space usage at the expense of querying speed of the resulting data. We decided to build the two format around the classical array-of-struct memory layout as we want to maximize the computability of the data. This is especially true in our application as we will never support individual columnar slicing and subset operations.

LD format

Tomahawk can output binary two data in the human-readable ld format by invoking the view command. The general schema for ld is below:

Column Description
FLAG Bit-packed boolean flags (see below)
CHROM_A Chromosome for marker A
POS_A Position for marker A
CHROM_B Chromosome for marker B
POS_B Position for marker B
REF_REF Inner product of (0,0) haplotypes
REF_ALT Inner product of (0,1) haplotypes
ALT_REF Inner product of (1,0) haplotypes
ALT_ALT Inner proruct of (1,1) haplotypes
D Coefficient of linkage disequilibrium
DPrime Normalized coefficient of linkage disequilibrium (scaled to [-1,1])
R Pearson correlation coefficient
R2 Squared pearson correlation coefficient
P Fisher's exact test P-value of the 2x2 haplotype contigency table
ChiSqModel Chi-squared critical value of the 3x3 unphased table of the selected cubic root (α, β, or δ)
ChiSqTable Chi-squared critical value of table (useful comparator when P = 0)

The 2x2 contingency table, or matrix, for the Fisher's exact test (P) for haplotypes look like this:

REF-A REF-B
REF-B A B
ALT-B C D

The 3x3 contigency table, or matrix, for the Chi-squared test for the unphased model looks like this:

0/0 0/1 1/1
0/0 A B C
0/1 D E F
1/1 G H J

The two FLAG values are bit-packed booleans in a single integer field and describe a variety of states a pair of markers can be in.

Bit position Numeric value One-hot Description
1 1 0000000000000001 Used phased math.
2 2 0000000000000010 Acceptor and donor variants are on the same contig.
3 4 0000000000000100 Acceptor and donor variants are far apart on the same contig.
4 8 0000000000001000 The output contingency matrix has at least one empty cell (referred to as complete).
5 16 0000000000010000 Output correlation coefficient is perfect (1.0).
6 32 0000000000100000 Output solution is one of >1 possible solutions. This only occurs for unphased pairs.
7 64 0000000001000000 Output data was generated in 'fast mode'.
8 128 0000000010000000 Output data is estimated from a subsampling of the total pool of genotypes.
9 256 0000000100000000 Donor vector has missing value(s).
10 512 0000001000000000 Acceptor vector has missing value(s).
11 1024 0000010000000000 Donor vector has low allele count (<5).
12 2048 0000100000000000 Acceptor vector has low allele count (<5).
13 4096 0001000000000000 Acceptor vector has a HWE-P value < 1e-4.
14 8192 0010000000000000 Donor vector has a HWE-P value < 1e-4.

Viewing and manipulating output data

It is possible to filter two output data by: 1) either start or end contig e.g. chr1, 2) position in that contig e.g. chr1:10e6-20e6; 3) have a particular contig mapping e.g. chr1,chr2; 4) interval mapping in both contigs e.g. chr1:10e3-10e6,chr2:0-10e6

1
tomahawk view -i 1kgp3_chr6_1_45.two -I 6:10e3-10e6,6:0-10e6 -H | head -n 6

Viewing all data

1
tomahawk view -i 1kgp3_chr6_1_45.two -H | head -n 6
FLAG CHROM_A POS_A CHROM_B POS_B REF_REF REF_ALT ALT_REF ALT_ALT D DPrime R R2 P ChiSqModel ChiSqTable
2059 6 89572 6 214654 4999 6 0 3 0.000597965 1 0.577004 0.332934 4.01511e-09 1667.33 0
2059 6 89573 6 214654 4999 6 0 3 0.000597965 1 0.577004 0.332934 4.01511e-09 1667.33 0
2059 6 122855 6 214654 4988 17 0 3 0.000596649 1 0.38664 0.149491 5.44908e-08 748.648 0
3 6 143500 6 212570 4421 387 82 118 0.0195352 0.54402 0.331324 0.109775 1.53587e-69 549.755 0
3 6 143500 6 213499 4419 387 84 118 0.0194949 0.537523 0.329068 0.108286 7.57426e-69 542.296 0

Slicing a range

1
tomahawk view -i 1kgp3_chr6_1_45.two -H -I 6:5e6-6e6 | head -n 6
FLAG CHROM_A POS_A CHROM_B POS_B REF_REF REF_ALT ALT_REF ALT_ALT D DPrime R R2 P ChiSqModel ChiSqTable
2063 6 5022893 6 73938 5000 7 0 1 0.000199362 1 0.353306 0.124825 0.00159744 625.125 0
2063 6 5020564 6 89339 5000 7 0 1 0.000199362 1 0.353306 0.124825 0.00159744 625.125 0
7 6 5018459 6 156100 1150 1283 388 2187 0.0804323 0.509361 0.348864 0.121706 1.17649e-138 609.504 0
7 6 5018691 6 156100 861 811 677 2659 0.0693918 0.339199 0.31898 0.101748 1.10017e-109 509.554 0
7 6 5018910 6 156100 861 811 677 2659 0.0693918 0.339199 0.31898 0.101748 1.10017e-109 509.554 0

Slicing matches in B string

1
tomahawk view -i 1kgp3_chr6_1_45.two -H -I 6:5e6-6e6,6:5e6-6e6 | head -n 6
FLAG CHROM_A POS_A CHROM_B POS_B REF_REF REF_ALT ALT_REF ALT_ALT D DPrime R R2 P ChiSqModel ChiSqTable
3 6 5000012 6 5073477 4988 9 5 6 0.0011915 0.544089 0.465743 0.216917 1.05037e-13 1086.32 0
1027 6 5000160 6 5072168 5000 2 4 2 0.000398404 0.4994 0.407677 0.166201 7.1708e-06 832.333 0
1027 6 5000160 6 5078340 5000 2 4 2 0.000398404 0.4994 0.407677 0.166201 7.1708e-06 832.333 0
3 6 5000482 6 5079621 4993 6 2 7 0.0013931 0.777199 0.64641 0.417846 3.9492e-18 2092.57 0
3 6 5000482 6 5082227 4994 6 1 7 0.00139362 0.874675 0.685808 0.470333 8.78523e-19 2355.43 0

As alluded to in the sorting section, sorted files have generally much faster query times. We can demonstrate this difference by using the sorted and unsorted data from the sliding window example above.

1
2
3
4
5
$ time tomahawk view -i 1kgp3_chr6_4mb.two -H -I 6:1e6-5e6,6:1e6-5e6 > /dev/null

real    3m1.029s
user    2m41.929s
sys     0m5.682s

Sorted

1
2
3
4
5
$ time tomahawk view -i 1kgp3_chr6_4mb_sorted.two -H -I 6:1e6-5e6,6:1e6-5e6 > /dev/null

real    0m38.167s
user    0m37.579s
sys     0m0.192s

Viewing ld data from the binary two file format and filtering out lines with a Fisher's exact test P-value < 1e-4, minor haplotype frequency < 5 and have both markers on the same contig (bit 2)

1
tomahawk view -i file.two -P 1e-4 -a 5 -f 2

Example output

FLAG CHROM_A POS_A CHROM_B POS_B REF_REF REF_ALT ALT_REF ALT_ALT D Dprime R R2 P ChiSqModel ChiSqTable
15 20 1314874 20 2000219 5002 5 0 1 0.00019944127 1 0.4080444 0.16650023 0.0011980831 0 833.83313
15 20 1315271 20 1992301 5005 2 0 1 0.00019956089 1 0.57723492 0.33320019 0.00059904153 0 1668.6665
15 20 1315527 20 1991024 5004 0 3 1 0.00019952102 1 0.49985018 0.2498502 0.00079872204 0 1251.2498
15 20 1315763 20 1982489 5006 0 1 1 0.00019960076 1 0.70703614 0.49990013 0.00039936102 0 2503.4999
15 20 1315807 20 1982446 5004 3 0 1 0.00019952102 1 0.49985018 0.2498502 0.00079872204 0 1251.2498

Example unphased output. Notice that estimated haplotype counts are now floating values and the bit-flag 1 is not set.

FLAG CHROM_A POS_A CHROM_B POS_B REF_REF REF_ALT ALT_REF ALT_ALT D Dprime R R2 P ChiSqModel ChiSqTable
46 20 1882564 20 1306588 5003.9996 2.0003999 1.0003999 0.99960008 0.00019936141 0.49950022 0.40779945 0.1663004 0.0011978438 0.0011996795 832.83241
46 20 1895185 20 1306588 5004.9998 1.0001999 1.0001999 0.99980011 0.0001994811 0.49970025 0.49970022 0.24970032 0.00079864228 0.00069960005 1250.4992
46 20 1306588 20 1901581 5003.9996 2.0003999 1.0003999 0.99960008 0.00019936141 0.49950022 0.40779945 0.1663004 0.0011978438 0.0011996795 832.83241
46 20 1901581 20 1306588 5003.9996 2.0003999 1.0003999 0.99960008 0.00019936141 0.49950022 0.40779945 0.1663004 0.0011978438 0.0011996795 832.83241
46 20 1306649 20 1885268 5006 1 1.2656777e-08 0.99999999 0.00019960076 1 0.70703614 0.49990013 0.00039936102 0.00039969285 2503.4999

Example output for forced phased math in fast mode. Note that only the REF_REF count is available, the fast math bit-flag is set, and all P and Chi-squared CV values are 0.

FLAG CHROM_A POS_A CHROM_B POS_B REF_REF REF_ALT ALT_REF ALT_ALT D Dprime R R2 P ChiSqModel ChiSqTable
67 20 1345922 20 1363176 4928 0 0 0 0.011788686 0.98334378 0.86251211 0.74392718 0 0 0
67 20 1345922 20 1367160 4933 0 0 0 0.011800847 0.98336071 0.89164203 0.79502547 0 0 0
67 20 1345958 20 1348347 4944 0 0 0 0.0092644105 0.97890127 0.85316211 0.72788554 0 0 0
67 20 1345958 20 1354524 4938 0 0 0 0.0092493389 0.86871892 0.80354655 0.64568704 0 0 0
75 20 1345958 20 1356626 4945 0 0 0 0.0033518653 0.99999994 0.51706308 0.26735422 0 0 0

Aggregating and visualizing datasets

Tomahawk generally output many millions to many hundreds of millions to billions of output linkage disequilibrium (LD) associations generated from many millions of input SNVs. It is technically very challenging to visualize such large datasets. Read more about aggregation in Tomahawk.

Aggregation Description
R Pearson correlation coefficient
R2 Squared pearson correlation coefficient
D Coefficient of linkage disequilibrium
Dprime Scaled coefficient of linkage disequilibrium
Dp Alias for dprime
P Fisher's exact test P-value of the 2x2 haplotype contigency table
Hets Number of (0,1) or (1,0) associations
Alts Number of (1,1) associations
Het Alias for hets
Alt Alias for alts
Reduction Description
Mean Mean number of aggregate
Max Largest number in aggregate bin
Min Smallest number in aggregate bin
Count Total number of records in a bin
N Alias for count
Total Sum total of aggregated number in a bin

In this example, we will aggregate the genome-wide data generated from the sliding example above by aggregating on R2 values and reducing into count in a (4000,4000)-bin space. If a bin have less than 50 observations (-c) we will drop all the value and report 0 for that bin.

1
twk aggregate -i 1kgp3_chr6_4mb_sorted.two -x 4000 -y 4000 -f r2 -r count -c 50 -t 4 -O b -o 1kgp3_chr6_4mb_aggregate.twa
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
Program:   tomahawk-7f8eef9b-dirty (Tools for computing, querying and storing LD data)
Libraries: tomahawk-0.7.0; ZSTD-1.3.8; htslib 1.9
Contact: Marcus D. R. Klarqvist <mk819@cam.ac.uk>
Documentation: https://github.com/mklarqvist/tomahawk
License: MIT
----------
[2019-01-25 14:55:43,299][LOG] Calling aggregate...
[2019-01-25 14:55:43,299][LOG] Performing 2-pass over data...
[2019-01-25 14:55:43,299][LOG] ===== First pass (peeking at landscape) =====
[2019-01-25 14:55:43,299][LOG] Blocks: 47,352
[2019-01-25 14:55:43,299][LOG] Uncompressed size: 50.192950 Gb
[2019-01-25 14:55:43,299][LOG][THREAD] Data/thread: 12.548238 Gb
[2019-01-25 14:55:43,299][PROGRESS] Time elapsed       Variants  Progress   Est. Time left
[2019-01-25 14:56:13,299][PROGRESS]      30,000s    344,930,000   72.8446%  11s
[2019-01-25 14:56:27,979][PROGRESS] 44,679s 473,514,826 (10,597,935 variants/s)
[2019-01-25 14:56:27,979][PROGRESS] Finished!
[2019-01-25 14:56:27,979][LOG] ===== Second pass (building matrix) =====
[2019-01-25 14:56:27,979][LOG] Aggregating 473,514,826 records...
[2019-01-25 14:56:27,979][LOG][THREAD] Allocating: 2.560000 Gb for matrices...
[2019-01-25 14:56:27,979][PROGRESS] Time elapsed       Variants  Progress   Est. Time left
[2019-01-25 14:56:57,979][PROGRESS]      30,000s    348,560,000   73.6112%  10s
[2019-01-25 14:57:12,499][PROGRESS] 44,520s 473,514,826 (10,635,926 variants/s)
[2019-01-25 14:57:12,499][PROGRESS] Finished!
[2019-01-25 14:57:13,713][LOG] Aggregated 473,514,826 records in 16,000,000 bins.
[2019-01-25 14:57:13,713][LOG] Finished.

You can now use the output binary format .twa in your downstream analysis. Alternatively, it is possible to output a human-readable (x,y)-matrix by setting the -O parameter to u. This tab-delimited matrix can now be loaded in any programming language and used as input for graphical visualizations or for analysis. It is easiest to use these files directly in rtomahawk, the R-bindings for tomahawk.