It is important to mention that this content is provided by Ragnar Groot Koerkamp, author of curiouscoding.nl.
Otherwise, bad things will happen.
Anyone should feel free to reach out to Ragnar directly for any clarifications.
Also, vote left and green! Save the planet!
With sassy, we can quickly do approximate string matching (ASM) or fuzzy
searching of small patterns such as 23bp barcodes or 23bp Crispr guides againsta fasta
file.
With multi-threading, searching a human genome takes as little as 30ms per
pattern on my laptop. And even less with the work-in-progress GPU version.
It would be great to also search in HPRCv2, but searching 466 full copies of a
human genome would be unnecessarily slow as they are highly redundant.
So here we go about deduping HPRCv2 into a spectrum preserving string set
(SPSS) that contains each k-mer at least once. We choose \(k=40\) to accommodate
searches of patterns of length around \(m=32\) while allowing up to \(k-m=8\) errors.
Here are my notes of developing pandedup, a tool for deduplicating (human) pangenomes.
As input we take the 3.3GB AGC file that we parse using the (vibe coded) ragc
Rust port.
The implementation uses multiple decoder instances to parse multiple samples in
parallel and sends them to a single thread (using a mpsc) that does the actual
deduplication.
We can use the simd-minimizers crate to extract syncmers, which (in my
notation) are substrings (windows) of length \(\ell = k+w-1\) that have their minimizer
(smallest k-mer) at the start or end. It is guaranteed that consecutive syncmers
overlap by at least \(k\) (although \(k-1\) is sufficient for our application).
We then incrementally build a hashset containing the 128-bit hashes of all syncmers seen so far.
When parsing a new contig, we take the union of all new syncmers and output the
corresponding substrings to the output.
For \(w=100\) we get:
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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
| process sample 0: 8.747084871s
new bp: 2934.064 Mbp (10620.9 bp/range)
new taken: 58.283 M (211.0 /range)
syncmers taken: 58.283 M (89.6%)
num_ranges: 0.276 M
output_bp: 2.934 Gbp (94.1%)
process sample 1: 11.235665811s
new bp: 541.665 Mbp (255.9 bp/range)
new taken: 8.061 M (3.8 /range)
syncmers taken: 66.344 M (51.2%)
num_ranges: 2.393 M
output_bp: 3.476 Gbp (56.3%)
process sample 2: 10.015981828s
new bp: 291.877 Mbp (238.7 bp/range)
new taken: 4.251 M (3.5 /range)
syncmers taken: 70.595 M (36.6%)
num_ranges: 3.616 M
output_bp: 3.768 Gbp (41.0%)
...
process sample 10: 9.581835449s
new bp: 87.308 Mbp (215.8 bp/range)
new taken: 1.230 M (3.0 /range)
syncmers taken: 91.651 M (13.4%)
num_ranges: 10.097 M
output_bp: 5.238 Gbp (15.7%)
...
process sample 25: 10.337954934s
new bp: 35.677 Mbp (202.2 bp/range)
new taken: 0.492 M (2.8 /range)
syncmers taken: 107.409 M (6.6%)
num_ranges: 15.382 M
output_bp: 6.362 Gbp (8.1%)
...
process sample 100: 11.310624909s
new bp: 20.608 Mbp (199.8 bp/range)
new taken: 0.282 M (2.7 /range)
syncmers taken: 148.540 M
num_ranges: 29.598 M
output_bp: 9.318 Gbp (3.1%)
...
process sample 200: 11.294771273s
new bp: 26.053 Mbp (203.0 bp/range)
new taken: 0.360 M (2.8 /range)
syncmers taken: 179.428 M
num_ranges: 40.543 M
output_bp: 11.552 Gbp (1.9%)
...
process sample 300: 10.509615192s
new bp: 11.874 Mbp (196.9 bp/range)
new taken: 0.162 M (2.7 /range)
syncmers taken: 201.636 M
num_ranges: 48.568 M
output_bp: 13.166 Gbp (1.5%)
...
process sample 400: 11.201425417s
new bp: 13.529 Mbp (198.9 bp/range)
new taken: 0.186 M (2.7 /range)
syncmers taken: 220.593 M
num_ranges: 55.347 M
output_bp: 14.536 Gbp (1.2%)
...
process sample 464: 12.979684545s
new bp: 7.983 Mbp (199.6 bp/range)
new taken: 0.109 M (2.7 /range)
syncmers taken: 231.548 M
num_ranges: 59.269 M
output_bp: 15.329 Gbp (1.1%)
process sample 465: 13.323026054s
new bp: 13.655 Mbp (197.8 bp/range)
new taken: 0.188 M (2.7 /range)
syncmers taken: 231.736 M
num_ranges: 59.338 M
output_bp: 15.343 Gbp (1.1%)
|
This finishes in 1:25 hours using 5h of CPU time.
The main bottleneck is that ragc does not support multithreaded decompression,
and so I have to create one decompressor instance per thread and these take a
lot of memory, so I can only have 4 at a time.
Observe:
- The first sample adds 2.9 Gbp of new sequence, and 90% of the syncmers is taken.
- Following samples quickly add less and less new content: around
30 Mbp/sample at sample 25, 20 Mbp/sample at sample 100, and just 10 Mbp/sample at sample 250. - Each newly added range has average length around 200bp, which is \(2w\) or \(w +
2.5k\). TODO: figure out which is mathematically correct.
- Most likely, each added range contains a single point mutation or structural variation.
It’s a bit long to add 200 bases if all we need is \(2k-1=79\) bases to
cover it.
- Each range is the span of 2.8 syncmers on average.
- After 466 samples, we have 231 million syncmers in total in our hashset.
- After 466 samples, we have 15.3 Gbp of data accumulated over 59M contigs, which is equivalent
to around 5 human genomes.
- If we bit-pack the data, it would be 3.8 GB, which is only slightly larger
than the 3.3 GB input file.
What’s somewhat annoying is that each mutation gives multiple new syncmers that
together span an unnecessarily large interval. Let’s try to fix that.
2 Syncmer-parse rather than syncmers (whoops)
Link to heading
Instead of fixed-length syncmers, we can use the minimizer parse: compute all
\((k, w)\) minimizers, and then extract the substrings spanned by each pair of
consecutive minimizers. These form strict substrings of the syncmers and should
be less prone to nearby mutations.
This results in:
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
| process sample 0: 7.513588348s
new bp: 2926.985 Mbp (9043.9 bp/range)
new taken: 57.840 M (178.7 /range)
syncmers taken: 57.840 M (88.9%)
num_ranges: 0.324 M
output_bp: 2.927 Gbp (93.9%)
process sample 1: 9.204434037s
new bp: 511.038 Mbp (236.8 bp/range)
new taken: 7.639 M (3.5 /range)
syncmers taken: 65.480 M (50.5%)
num_ranges: 2.481 M
output_bp: 3.438 Gbp (55.7%)
process sample 2: 7.965897865s
new bp: 418.844 Mbp (222.6 bp/range)
new taken: 6.141 M (3.3 /range)
syncmers taken: 71.621 M (37.1%)
num_ranges: 4.363 M
output_bp: 3.857 Gbp (42.3%)
...
process sample 10: 8.435886889s
new bp: 179.178 Mbp (213.1 bp/range)
new taken: 2.589 M (3.1 /range)
syncmers taken: 89.146 M (13.1%)
num_ranges: 10.340 M
output_bp: 5.083 Gbp (15.3%)
...
process sample 25: 8.458732794s
new bp: 32.617 Mbp (182.8 bp/range)
new taken: 0.446 M (2.5 /range)
syncmers taken: 103.603 M (6.3%)
num_ranges: 15.712 M
output_bp: 6.119 Gbp (7.8%)
...
process sample 100: 8.336943788s
new bp: 18.526 Mbp (178.4 bp/range)
new taken: 0.247 M (2.4 /range)
syncmers taken: 140.545 M (2.2%)
num_ranges: 30.125 M
output_bp: 8.816 Gbp (2.9%)
...
process sample 300: 8.166004583s
new bp: 10.600 Mbp (175.1 bp/range)
new taken: 0.140 M (2.3 /range)
syncmers taken: 187.128 M (1.0%)
num_ranges: 49.271 M
output_bp: 12.285 Gbp (1.4%)
|
Compared to the syncmer parse:
- Ranges (mostly for point mutations) are 10% shorter.
- We get 10% less new output bases.
- We sample fewer minimizer phrases than syncmers.
(Not sure anymore exactly which parameters this has.)
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
| process sample 0: 7.682570099s
new bp: 2912.950 Mbp (6754.1 bp/range)
new taken: 56.843 M (131.8 /range)
syncmers taken: 56.843 M (88.2%)
num_ranges: 0.431 M
output_bp: 2.913 Gbp (93.4%)
process sample 1: 9.894358837s
new bp: 434.693 Mbp (187.0 bp/range)
new taken: 6.091 M (2.6 /range)
syncmers taken: 62.934 M (49.0%)
num_ranges: 2.755 M
output_bp: 3.348 Gbp (54.3%)
process sample 2: 8.430858872s
new bp: 350.879 Mbp (177.4 bp/range)
new taken: 4.805 M (2.4 /range)
syncmers taken: 67.739 M (35.4%)
num_ranges: 4.733 M
output_bp: 3.699 Gbp (40.6%)
...
process sample 10: 8.516717746s
new bp: 65.566 Mbp (159.8 bp/range)
new taken: 0.850 M (2.1 /range)
syncmers taken: 81.128 M (12.0%)
num_ranges: 10.792 M
output_bp: 4.707 Gbp (14.1%)
...
process sample 100: 10.845207061s
new bp: 14.189 Mbp (144.6 bp/range)
new taken: 0.175 M (1.8 /range)
syncmers taken: 118.237 M (1.8%)
num_ranges: 30.154 M
output_bp: 7.550 Gbp (2.5%)
...
process sample 200: 10.901156258s
new bp: 18.143 Mbp (146.5 bp/range)
new taken: 0.222 M (1.8 /range)
syncmers taken: 137.706 M (1.1%)
num_ranges: 40.685 M
output_bp: 9.110 Gbp (1.5%)
|
Before, we were using \(w=100\), \(k=40\) minimizers, but in fact it is sufficient
to use shorter minimizers with \(k’=8\) to extract ‘key positions’ at most 100
apart, and then simply extract sufficiently overlapping subsequences starting at
these positions and ending 40 beyond the next sampled position.
For random minimizers, \(k’=8\) is large enough to not loose any density, while
being more stable under changes.
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
| process sample 0: 9.144635842s
new bp: 2909.696 Mbp (6250.1 bp/range)
new taken: 58.868 M (126.4 /range)
syncmers taken: 58.868 M (87.9%)
num_ranges: 0.466 M
output_bp: 2.910 Gbp (93.3%)
process sample 1: 8.989641454s
new bp: 419.542 Mbp (177.9 bp/range)
new taken: 5.700 M (2.4 /range)
syncmers taken: 64.568 M (48.4%)
num_ranges: 2.824 M
output_bp: 3.329 Gbp (54.0%)
process sample 2: 8.699237103s
new bp: 338.204 Mbp (169.4 bp/range)
new taken: 4.485 M (2.2 /range)
syncmers taken: 69.053 M (34.8%)
num_ranges: 4.820 M
output_bp: 3.667 Gbp (40.2%)
process sample 3: 9.111089474s
new bp: 159.663 Mbp (161.9 bp/range)
new taken: 2.078 M (2.1 /range)
syncmers taken: 71.132 M (27.0%)
num_ranges: 5.806 M
output_bp: 3.827 Gbp (31.5%)
...
process sample 10: 9.04328886s
new bp: 143.675 Mbp (168.2 bp/range)
new taken: 1.911 M (2.2 /range)
syncmers taken: 81.645 M (11.6%)
num_ranges: 10.895 M
output_bp: 4.640 Gbp (13.9%)
...
process sample 465: 10.648589872s
new bp: 9.454 Mbp (143.5 bp/range)
new phrases: 0.116 M (1.8 /range)
total phrases: 170.376 M (0.6%)
num_ranges: 58.554 M
output_bp: 11.732 Gbp (0.8%)
|
So 15.3 to 11.7; pretty good gain :)
In 4344s overall on 3 threads.
Turns out AGC is returning u8 vectors of 0123, not ACTG, and so the
minimizer/PFP computation was only using 1 of the 2 bits (because going from
ACTG to 0123 normally shifts right by 1).
Fixing that gives:
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
| process sample 0: 9.24019449s
new bp: 2908.803 Mbp (6423.3 bp/range)
new phrases: 57.031 M (125.9 /range)
total phrases: 57.031 M (86.4%)
num_ranges: 0.453 M
output_bp: 2.909 Gbp (93.3%)
process sample 1: 11.256846815s
new bp: 420.293 Mbp (178.5 bp/range)
new phrases: 5.640 M (2.4 /range)
total phrases: 62.671 M (47.9%)
num_ranges: 2.807 M
output_bp: 3.329 Gbp (54.0%)
process sample 2: 9.406656924s
new bp: 221.560 Mbp (169.2 bp/range)
new phrases: 2.883 M (2.2 /range)
total phrases: 65.554 M (33.7%)
num_ranges: 4.117 M
output_bp: 3.551 Gbp (38.6%)
process sample 3: 9.715269299s
new bp: 277.431 Mbp (166.4 bp/range)
new phrases: 3.604 M (2.2 /range)
total phrases: 69.158 M (26.8%)
num_ranges: 5.784 M
output_bp: 3.828 Gbp (31.5%)
...
process sample 10: 8.702274402s
new bp: 142.963 Mbp (168.0 bp/range)
new phrases: 1.888 M (2.2 /range)
total phrases: 79.541 M (11.6%)
num_ranges: 10.853 M
output_bp: 4.640 Gbp (13.9%)
|
4.2 With read-write lock of hashset
Link to heading
- Same output, but only 500ms instead of 10s lock on the hashmap.
Let’s also try the approach of sampling all kmers with fractional hash value at
most \(2/(w+1)\) (to have the same overall density).
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
| process sample 0: 10.844735133s
new bp: 2927.617 Mbp (6567.0 bp/range)
new taken: 72.765 M (163.2 /range)
syncmers taken: 72.765 M (91.7%)
num_ranges: 0.446 M
output_bp: 2.928 Gbp (93.9%)
process sample 1: 10.579291664s
new bp: 550.039 Mbp (202.1 bp/range)
new taken: 6.931 M (2.5 /range)
syncmers taken: 79.696 M (51.6%)
num_ranges: 3.168 M
output_bp: 3.478 Gbp (57.4%)
process sample 2: 11.18861265s
new bp: 301.227 Mbp (202.4 bp/range)
new taken: 3.574 M (2.4 /range)
syncmers taken: 83.269 M (35.8%)
num_ranges: 4.656 M
output_bp: 3.779 Gbp (41.5%)
process sample 3: 12.66143639s
new bp: 189.336 Mbp (196.2 bp/range)
new taken: 2.146 M (2.2 /range)
syncmers taken: 85.415 M (27.5%)
num_ranges: 5.622 M
output_bp: 3.968 Gbp (32.7%)
...
process sample 10: 11.550021873s
new bp: 78.163 Mbp (195.8 bp/range)
new taken: 0.816 M (2.0 /range)
syncmers taken: 96.193 M (11.3%)
num_ranges: 10.578 M
output_bp: 4.950 Gbp (14.9%)
|
Overall, this is worse than the minimizer variant.
6 Minimizers, with smaller \(w=50\), \(k=40\), \(k’=8\)
Link to heading
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
| process sample 0: 15.900242517s
new bp: 2887.140 Mbp (4294.3 bp/range)
new phrases: 113.153 M (168.3 /range)
total phrases: 113.153 M (88.6%)
num_ranges: 0.672 M
output_bp: 2.887 Gbp (92.6%)
process sample 1: 18.765048945s
new bp: 335.610 Mbp (133.2 bp/range)
new phrases: 8.261 M (3.3 /range)
total phrases: 121.414 M (47.9%)
num_ranges: 3.191 M
output_bp: 3.223 Gbp (52.2%)
process sample 2: 15.640636325s
new bp: 265.662 Mbp (127.1 bp/range)
new phrases: 6.365 M (3.0 /range)
total phrases: 127.779 M (34.0%)
num_ranges: 5.282 M
output_bp: 3.488 Gbp (38.3%)
push sample 8 (3.0 Gbp): read: 23.327796953s minis: 5.11133427s
process sample 3: 16.567902729s
new bp: 123.190 Mbp (122.6 bp/range)
new phrases: 2.898 M (2.9 /range)
total phrases: 130.677 M (26.1%)
num_ranges: 6.287 M
output_bp: 3.612 Gbp (29.7%)
...
process sample 10: 16.139335034s
new bp: 111.328 Mbp (129.2 bp/range)
new phrases: 2.709 M (3.1 /range)
total phrases: 145.220 M (10.7%)
num_ranges: 11.407 M
output_bp: 4.233 Gbp (12.7%)
|
So maybe 10-20% smaller overall result, but 2x more phrases and thus 2x larger hashset.
7 \(k=64\), \(w=100\), \(k’=8\)
Link to heading
Gives a 15.0 Gbp file, up from 11.8 Gbp for \(k=40\).
process sample 465: 135.466742ms
new bp: 13.638 Mbp (192.3 bp/range)
new phrases: 0.162 M (2.3 /range)
total phrases: 203.864 M (30.3%)
num_ranges: 59.379 M
output_bp: 14.860 Gbp (1.1%)
604s walltime on 64 threads (backus)
With canonical (in 582s):
process sample 465: 135.578722ms
new bp: 14.003 Mbp (194.4 bp/range)
new phrases: 0.163 M (2.3 /range)
total phrases: 202.118 M (30.6%)
num_ranges: 59.692 M
output_bp: 15.063 Gbp (1.1%)
With non-canonical anti-lex kmers (in 603s):
1
2
3
4
5
6
| process sample 465: 129.883241ms
new bp: 13.647 Mbp (191.0 bp/range)
new phrases: 0.156 M (2.2 /range)
total phrases: 196.372 M (30.3%)
num_ranges: 59.486 M
output_bp: 14.848 Gbp (1.1%)
|
After speading up the ragc reader: 438s.
After a smaller initial batch to reduce 64x 2s wait: 383s
We can shard the hashtable into 256 parts and process the phrases in each sample
part-by-part.
If we randomize the order of the parts in each sample, we get much larger output
(398s).
1
2
3
4
5
6
| process sample 465: 12.971602ms
new bp: 10.967 Mbp (212.6 bp/range)
new phrases: 0.133 M (2.6 /range; 100.0 %)
total phrases: 196.372 M (100.0%)
num_ranges: 107.403 M
output_bp: 18.175 Gbp (1.3%)
|
Back to linear order in 362s:
1
2
3
4
5
6
7
| push sample 465 (3.0 Gbp): read: 11.05s minis: 4.82s phrases: 5.97s sort: 1.45s lookups: 9.07s sort: 0.13s
process sample 465: 15.361753ms
new bp: 13.014 Mbp (189.5 bp/range)
new phrases: 0.148 M (2.1 /range; 100.0 %)
total phrases: 196.372 M (100.0%)
num_ranges: 59.484 M
output_bp: 14.848 Gbp (1.1%)
|
128 threads instead of 64, to avoid low utilization (of only 20 cores at times):
380s; 750GB RAM.
With RwLock instead of Mutex: 366s also. Slightly larger/more ranges because of
unordered insertions
1
2
3
4
5
6
| process sample 465: 14.683933ms
new bp: 11.650 Mbp (190.2 bp/range)
new phrases: 0.132 M (2.2 /range; 100.0 %)
total phrases: 196.372 M (100.0%)
num_ranges: 61.175 M
output_bp: 14.968 Gbp (1.1%)
|
7.2 Per-contig processing and per-thread output
Link to heading
Instead of sending all output to a single thread, we can output directly from
the worker thread. We can also process only single contigs at a time to
significantly reduce memory usage: 471s.
Note that we also process the reference first on a single thread.
1
2
3
4
5
6
| push sample 463 (3.1 Gbp 85 ctg): read: 6.58s minis: 4.49s phrases: 1.58s sort: 0.27s lookups: 23.59s sort: 0.13s lock: 0.26s output: 0.15s
new bp: 11.004 Mbp (212.7 bp/contig)
new phrases: 0.133 M (2.6 /contig; 0.2%)
unique phrases: 196.372 M (0.7%)
num_contigs: 60.215 M
output_bp: 14.899 Gbp (1.1%)
|
7.3 Pre-filtering against reference:
Link to heading
We store the reference and the position of each phrase in it.
When parsing a new sequence, if the phrase already occurs in the reference, we
then greedily extend from there and skip any phrases contained in this region.
307s
1
2
3
4
5
6
| push sample 465 (3.0 Gbp 79 ctg): read: 7.23s minis: 4.79s phrases: 6.70s sort: 0.09s lookups: 1.84s sort: 0.04s lock: 1.22s output: 0.14s
new bp: 13.141 Mbp (189.5 bp/contig)
new phrases: 0.149 M (2.1 /contig; 0.2%; 68.0% filtered)
unique phrases: 196.372 M (0.7%)
num_contigs: 60.225 M
output_bp: 14.899 Gbp (1.1%)
|
8 \(k=64\), \(w=25\), \(k’=8\)
Link to heading
Gives a 11.5 Gbp file, or 3.1 GB after zst-compression.
96 threads, anti-lex: 1000s; 950GB RAM
1
2
3
4
5
6
| process sample 465: 13.606797ms
new bp: 8.943 Mbp (138.3 bp/range)
new phrases: 0.347 M (5.4 /range; 100.0 %)
total phrases: 552.892 M (100.0%)
num_ranges: 57.789 M
output_bp: 11.288 Gbp (0.8%)
|
Latest version, in 386s and <100GB RAM:
1
2
3
4
5
6
| push sample 465 (3.0 Gbp 79 ctg): read: 7.37s minis: 5.48s phrases: 11.59s sort: 0.21s lookups: 4.03s sort: 0.10s lock: 0.22s output: 0.15s
new bp: 9.042 Mbp (138.3 bp/contig)
new phrases: 0.351 M (5.4 /contig; 0.2%; 81.9% filtered away)
unique phrases: 552.892 M (0.5%)
num_contigs: 58.045 M
output_bp: 11.312 Gbp (0.8%)
|
With \(w=10\), in 587s and 96GB RAM:
1
2
3
4
5
6
| push sample 465 (3.0 Gbp 79 ctg): read: 7.45s minis: 5.70s phrases: 17.49s sort: 0.32s lookups: 7.64s sort: 0.20s lock: 0.33s output: 0.22s
new bp: 8.250 Mbp (128.1 bp/contig)
new phrases: 0.706 M (11.0 /contig; 0.1%; 86.6% filtered away)
unique phrases: 1179.914 M (0.5%)
num_contigs: 57.766 M
output_bp: 10.625 Gbp (0.8%)
|
- 400 s
- <100 GB RAM
- 1.9 GB .zst output
1
2
3
4
5
6
| push sample 465 (3.0 Gbp 79 ctg): read: 4.87s minis: 4.08s phrases: 17.66s sort: 0.28s lookups: 6.14s sort: 0.09s lock: 0.16s output: 0.15s
new bp: 3.652 Mbp (67.0 bp/contig)
new phrases: 0.317 M (5.8 /contig; 0.1%; 85.7% filtered away)
unique phrases: 777.793 M (0.3%)
num_contigs: 54.224 M
output_bp: 6.390 Gbp (0.5%)
|
10.1 Using it in deacon: bug in ambiguous bases
Link to heading
Bede used hprcv2.k64.zst in deacon. For k=31 and w=21, it gives
353,721,077 minimizers, whereas running on the full hprcv2 .agc directly
gives 353,690,959 minimizers, which is 30118 smaller, and so I’m
over-reporting kmers. Hrmmm…
Added a bunch of tests to pandedup, but those all seem to work fine.
It turns out I assumed that AGC never returns ambiguous bases, but that’s not
true. Thus, doing b"ACGT"[b%4] quietly collapses ambiguous bases.
Fixing this should resolve the additional kmers.
*A second bug c
First, we do the rough dedup using pandedup:
1
2
3
4
5
6
7
8
9
10
| > cargo run -r -- -k 64 -w 100 hprcv2.agc
...
push sample 465 (3.0 Gbp 79 ctg): read: 0.00s minis: 4.95s phrases: 5.08s sort: 0.08s lookups: 1.05s sort: 0.04s lock: 0.00s output: 0.10s
new bp: 13.039 Mbp (189.5 bp/contig)
new phrases: 0.148 M (2.1 /contig; 0.2%; 68.0% filtered away)
unique phrases: 196.372 M (0.7%)
num_contigs: 59.592 M
output_bp: 14.856 Gbp (1.1%)
|
Takes 1700s on my laptop and outputs 14.9Gbp. On a server, it runs in 230s
with 64 threads.
ggcat on the output takes forever, see https://github.com/algbio/ggcat/issues/76
(now fixed; it had a bug for even \(k\)).
1
| ggcat build --greedy-matchtigs -o hprcv2.dedup.ggcat.fa --kmer-length 64 -m 512 -j 128 hprcv2.dedup.fa
|
Instead, we can use Cuttlefish to build unitigs in 20min. Some small remarks here https://github.com/COMBINE-lab/cuttlefish/issues/70.
1
| cuttlefish build --ref --seq hprcv2.dedup.fa --kmer-len 63 --threads 64 --work-dir work --output hprcv2.unitigs.fa -m 512
|
- Sadly, this increases the file size from 15GB to 18GB.
- Note that Cuttlefish only takes odd k, so we have \(k=63\) now.
Now, let’s shrink things using greedy matchtigs (note that matchtigs uses
--greedytigs for greedy matchtigs): 1h, small issue here https://github.com/algbio/matchtigs/issues/8:
1
| matchtigs -k 63 -t 64 --fa-in hprcv2.unitigs.fa --greedytigs-fa-out hprcv2.greedytigs.fa
|
This outputs 8.7GB (down from 10.6GB with \(w=25\)), or 2.38GB after zstd compression!
Download the result here:
https://ragnargrootkoerkamp.nl/upload/hprcv2-k63-greedytigs.fa.zst
GGCAT implements fout types of tig output:
(Refer to this post summarizing types of tigs.)
Unitigs: the default.
Simplitigs: just concatenate unitigs, reimplemented and very fast.
Eulertigs: optimally concatenate unitigs, reimplemented and very fast.
Greedy matchtigs: the smallest and directly forwards to the matchtigs library.
- These also contract kmers at a distance, causing duplicate kmers.
It turns out that ggcat had an issue for even k that has now been fixed. For
\(k=63\), it outputs 8.9GB (slightly larger only because of file headers) and runs in 50min with:
1
2
3
4
5
6
7
8
| phase: reads bucketing => 12.59s
phase: kmers merge => 19.18s
phase: unitigs joining => 38.37s
phase: maximal unitigs links building [step 1] => 1.82s
phase: maximal unitigs links building [step 2] => 5.21s
phase: maximal unitigs links building [step 3] => 124.00s
phase: greedy matchtigs building [step1] => 2682.47s
phase: greedy matchtigs building [step2] => 160.27s
|
Note that the unitig part takes only 200s compared to the 20
minutes used by cuttlefish. The matchtigs part takes around 47min, which is also
~10min faster than running matchtigs directly. I assume that’s because ggcat has
more optimized I/O, and that either way the matchtigs part is largely single-threaded.
10.9GB in 160s
1
2
3
4
5
6
7
| TOTAL TIME: 160.60s
Max virtual fs usage: 9.48 GiB
Max disk usage: 0.00 octets
Final stats:
phase: reads bucketing => 12.47s
phase: kmers merge => 12.98s
phase: unitigs joining => 135.15s
|
10.7GB in 162s
- This crashed for
k=64, so this is for k=63 instead.
1
2
3
4
5
6
7
8
9
10
| TOTAL TIME: 162.43s
Max virtual fs usage: 9.46 GiB
Max disk usage: 0.00 octets
Final stats:
phase: reads bucketing => 12.40s
phase: kmers merge => 12.76s
phase: unitigs joining => 41.64s
phase: eulertigs building part 1 => 35.81ms
phase: eulertigs building part 2 => 9.89ms
phase: eulertigs building part 3 => 95.58s
|
TODO 10.3.4 Masked superstrings
Link to heading
There was a second bug: ggcat would throw away contigs shorter than k, but we
explicitly need those if we use the index for a smaller final value of k.
Thus, deacon would give 353,690,956 minimizers rather than 353,690,959.
Explicitly outputting short kmers to a separate file, indexing that
with deacon, and merging the two indices recovers the last 3 missing minimizers
from the index obtained when running on the .agc directly.
See examples/filter_short_contigs.rs in pandedup, and this (slop!) PR/FR to add
--preserve-short-contigs to ggcat.
11.2 Rerun by calling ggcat via the API
Link to heading
Instead of writing intermediate files to disk, it’s easier and faster to just call the
ggcat API directly.
This takes as input the hprcv2.agc to pandedup, and then repeatedly calls
ggcat on this output in different modes.
I had codex create a ggcat branch that uses the GGCAT API to output all
possible types of output:
1
| cargo run -r -- hprcv2/hprcv2.agc -k 63 -w 100 -j 64 -o hprcv2/out/hprcv2.dedup.fa.zst --ggcat
|
- pandedup output: 3.5 min, 15GB
- unitigs: 69s, 21GB (incl ~3GB record headers)
- simplitigs: 60s, 10.746GB (~1GB record headers)
- eulertigs: 67s, 10.745GB (~1GB record headers)
- (I would have expected a larger difference here.)
- greedy matchtigs: 5200s, 8.9GB (of which ~0.2GB headers), 2.6GB zstd compressed
- takes 75min of Dijkstra. Earlier this was much faster; maybe some badly
balanced buckets? Most of the time, only 2 threads are active.
- 4min to find Eulerian cycle.
- masked superstring: 5h, 7.8GB (no headers), 2.5GB zstd compressed
- Builds directly from eulertigs:
1
| kmercamel compute -k 63 -o hprcv2.masked.msfa -S hprcv2-k63-eulertigs.fa
|
Note that ggcat uses a fast reimplementation of the simplitigs and eulertigs
methods, but falls back to the matchtigs library for greedy matchtigs.
Somewhat surprisingly, after plain zstd compression, the masked superstring is
only 60MB smaller, which is possibly less than the overhead of numbered records
in the greedy matchtigs.
stats: (note that ggcat outputs include record lengths in the header, at ~9B/record)
| file | time (min) | tigs | size (GB) | zstd compressed size (GB) |
|---|
| hprcv2.pandedup.fa | 3.5 | 60M | 14.9 | 3.9 |
| hprcv2.ggcat-unitigs.fa | 1 | 185M | 21.2 | |
| hprcv2.ggcat-simplitigs.fa | 1 | 59M | 10.7 | |
| hprcv2.ggcat-eulertigs.fa | 1 | 60M | 10.7 | |
| hprcv2.ggcat-greedy-matchtigs.fa | 90 | 27M | 8.9 | 2.6 |
| hprcv2.kmercamel-masked.fa | 500 | 1 | 7.8 | 2.5 |
12.1 Faster deacon index build
Link to heading
- Make the output deterministic by sorting minimizer values by target
bucket and value before writing them.
- Add a
deacon index reorder <path> command to re-sort existing index files - Parallellize the build in two steps:
- let each thread collect its own minimizer values and then single-threadedly
collect those.
- Secondly, re-shard by hash-table order in parallel: each thread builds a
Vec<Vec<T>> with one inner Vec for each target shard, and then appends
each inner list to the global Vec<Mutex<Vec<T>>>. Then, each global Vec
is sorted and deduplicated.
1
| cargo run -r -- index build -t 12 $d/chm13v2.0.fa -o index-new
|
Before: 23.9s (66s CPU)
1
2
3
4
5
6
7
8
9
10
| 19:42:19 INFO Deacon v0.18.0; mode: build; input: single; options: threads=12
19:42:19 INFO Building index (k=31, w=15)
Processed 7 sequences (1093031482bp), 116696808 minimizers
Processed 12 sequences (2100591497bp), 223674031 minimizers
Processed 25 sequences (3117292070bp), 316436671 minimizers
19:42:35 INFO Collecting minimizers into a single set
19:42:35 INFO Indexed 316436671 minimizers from 25 record(s) (3117292070bp)
19:42:35 INFO Sorting minimizers..
19:42:40 INFO Writing minimizers..
19:42:43 INFO Completed build in 23.90s
|
Collect minimizers per thread: 42.9s
1
2
3
4
5
6
7
8
9
10
11
| 19:37:33 INFO Deacon v0.18.0; mode: build; input: single; options: threads=12
19:37:33 INFO Building index (k=31, w=15)
19:37:39 INFO Processed 7 sequences (1032781085bp)
19:37:41 INFO Processed 15 sequences (2048986257bp)
19:37:45 INFO Processed 25 sequences (3117292070bp)
19:37:46 INFO Collecting minimizers into a single set
19:38:09 INFO Indexed 316436671 minimizers from 25 record(s) (3117292070bp)
19:38:09 INFO Sorting minimizers..
19:38:12 INFO Checking that values are sorted by hash.
19:38:13 INFO Writing minimizers..
19:38:16 INFO Completed build in 42.92s
|
Multi-threaded re-sharding at the end: 19.2s (140s CPU time)
1
2
3
4
5
6
7
8
9
10
11
| 19:34:17 INFO Deacon v0.18.0; mode: build; input: single; options: threads=12
19:34:17 INFO Building index (k=31, w=15)
19:34:23 INFO Processed 7 sequences (1032781085bp)
19:34:25 INFO Processed 15 sequences (2207406112bp)
19:34:29 INFO Processed 25 sequences (3117292070bp)
19:34:30 INFO Re-sharding and sorting minimizers
19:34:30 INFO Re-shard minimizers
19:34:31 INFO Sort shards
19:34:33 INFO Indexed 316436671 minimizers from 25 record(s) (3117292070bp)
19:34:33 INFO Writing minimizers..
19:34:36 INFO Completed build in 19.22s
|
To check, all have the same sha256:
1
2
| > sha256sum index-new
bdbfcd394ab8e89e11cdb37e997a582197abebca3a6dcf53c4ba77cd7de5e162 index-new
|
Vibe-coded branch: 9.1s (47s CPU)
(Note that this one has a different hash… apparently because it sorts the
shards by bucket value, but not the elements within the shards.)
Time to figure out what makes the big difference here..
12.1.2 Learning from the vibe-coded version
Link to heading
The main difference is that we spend a lot of time inserting minimizers into
hashmaps. Pushing them onto a list and regularly sorting+deduplicating those
is more efficient. In particular, this version has a length-1024 thread-local
buffers for every shard that are flushed to the global shards once full. Then,
global shards are sorted and deduplicated every time they double in size (and
have at least 4M elements).
Collect-to-vec:
So, let’s first collect onto a local Vec<T> that is pushed to a global Vec<T> per shard that we deduplicate
once, at the end. 20s
This is already much better: 11s
1
2
3
4
5
6
7
8
| 20:29:02 INFO Processed 7 sequences (1232028657bp)
20:29:02 INFO Processed 12 sequences (2100591497bp)
20:29:03 INFO Processed 25 sequences (3117292070bp)
20:29:03 INFO Re-shard minimizers
20:29:06 INFO Sort shards
20:29:08 INFO Indexed 316436671 minimizers from 25 record(s) (3117292070bp)
20:29:08 INFO Writing minimizers..
20:29:11 INFO Completed build in 10.85s
|
voracious_radix_sort:
Now let’s use voracious_radix_sort::RadixSort::voracious_sort instead of
Vec::sort_unstable. 10.9s
It turns out the sorting to deduplicate minimizers isn’t so much a bottleneck
since at least on this dataset there is not too much duplication.
Faster sort-by-hash. Currently, I do shard.sort_unstable_by_key(|x| (bucket(x), x)) to deterministically sort the hashing. Let’s instead to sort
only by bucket(x) first, and then do a second pass to sort values with the
same bucket. 11.5s.
Let’s just ignore these last 2 ideas.
1
2
| cargo build -r
time cargo run -r -- index build -t 12 ./hprcv2-k63-greedytigs.fa.zst -o index
|
Baseline: 82s (CPU 147s)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
| 21:17:53 INFO Building index (k=31, w=15)
Processed 4956378 sequences (1000001644bp), 76332117 minimizers
Processed 9840366 sequences (2000024274bp), 140810438 minimizers
Processed 14452692 sequences (3000017524bp), 198972527 minimizers
Processed 18651198 sequences (4000001222bp), 252520490 minimizers
Processed 22183944 sequences (5000017030bp), 303053561 minimizers
Processed 24701748 sequences (6000046266bp), 352937691 minimizers
Processed 26150688 sequences (7000025386bp), 404442635 minimizers
Processed 26897844 sequences (8000047407bp), 457857180 minimizers
21:19:01 INFO Collecting minimizers into a single set
21:19:01 INFO Indexed 480863905 minimizers from 27092669 record(s) (8424467217bp)
21:19:02 INFO Sorting minimizers..
21:19:08 INFO Writing minimizers..
21:19:15 INFO Completed build in 82.56s
|
New version: 35s (CPU 157s)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
| 21:16:10 INFO Deacon v0.18.0; mode: build; input: single; options: threads=12
21:16:10 INFO Building index (k=31, w=15)
21:16:13 INFO Processed 4956492 sequences (1000022657bp)
21:16:15 INFO Processed 9840252 sequences (2000003572bp)
21:16:18 INFO Processed 14452692 sequences (3000011065bp)
21:16:21 INFO Processed 18651312 sequences (4000026167bp)
21:16:23 INFO Processed 22183944 sequences (5000019134bp)
21:16:26 INFO Processed 24701748 sequences (6000046266bp)
21:16:28 INFO Processed 26150688 sequences (7000025386bp)
21:16:31 INFO Processed 26897844 sequences (8000024736bp)
21:16:33 INFO Dedup shards
21:16:36 INFO Re-shard minimizers
21:16:38 INFO Sort shards
21:16:41 INFO Indexed 480863905 minimizers from 27092669 record(s) (8424467217bp)
21:16:41 INFO Writing minimizers..
21:16:45 INFO Completed build in 35.31s
|
LLM version: 35.6s (CPU 132s)
1
2
3
| ...
Indexed 480863905 minimizers from 27092669 record(s) (8424467217bp)
Completed build in 35.65s
|
Currently, we only dedup everything at the very end, but this is tricky when
there are many duplicate minimizers. This has two problems: identical minimizers might
occur very frequently within a thread, and the same minimizer can occur in all
threads. Thus, we need constant global deduplication to prevent storing many
minimizers once per thread. Indeed this is what the LLM does. It uses a
“deduplicating vector” that does sort+dedup every time it doubles in size.
We also reduce the size of per-thread buffers to at most 1024 elements per shard
and actively push to the global buffers. The total time remains higher than the
35s we had originally:
DeduppingVec: 40s (CPU 180s)
With .fa instead of .fa.zst input, this reduces to 28s!
If we further increase the threshold where dedupping starts to happen to 4M
elements per shard, then the time reduces to 20s (153s CPU) and it nearly matches the LLM
code again :)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
| > time ./target/release/deacon index build -t 12 ./hprcv2-k63-greedytigs.fa -o index
22:27:00 INFO Building index (k=31, w=15)
22:27:01 INFO Processed 4956492 sequences (1000020972bp)
22:27:02 INFO Processed 9840252 sequences (2000003572bp)
22:27:03 INFO Processed 14452692 sequences (3000005896bp)
22:27:04 INFO Processed 18651312 sequences (4000026167bp)
22:27:05 INFO Processed 22183944 sequences (5000006629bp)
22:27:05 INFO Processed 24701748 sequences (6000046266bp)
22:27:06 INFO Processed 26150688 sequences (7000000451bp)
22:27:07 INFO Processed 26897844 sequences (8000047407bp)
22:27:07 INFO Dedup shards
22:27:11 INFO Re-shard minimizers
22:27:13 INFO Sort shards
22:27:16 INFO Indexed 480863905 minimizers from 27092669 record(s) (8424467217bp)
22:27:16 INFO Writing minimizers..
22:27:20 INFO Completed build in 20.55s
|
12.1.5 Prevent spamming the terminal
Link to heading
Also added a quick commit to redirect output from stdout to /dev/null when
no explicit output is given and the output is a terminal rather than a pipe/redirect.
Decompressing hprcv2-k63-greedytigs.fa.zst from 2.6GB to 8.7GB takes 15s on a
single thread,
while the entire deacon pipeline takes 10s when using .fa input.
There happens to be a (vibe coded?) PR from this week adding multithreading to
zstd: https://github.com/facebook/zstd/pull/4796.
Why not gzip? Because gzip hprcv2-k63-greedytigs.fa takes over 20 minutes.
1
| time ./target/release/deacon index build -t 12 ../pandedup/hprcv2-k63-greedytigs.<ext>
|
(excluding writing a file to disk now)
| extension | size (GB) | time (s) | CPU (s) |
|---|
| .fa | 8.7 | 18 | 157 |
| .fa.zst | 2.9 | 33 | 156 |
| .cbq | 2.0 | 19 | 164 |
Do note that .cbq is stored without record-id headers here.
Cbq is as small as zstd, but as fast as plain fasta!
- Various ggcat, kmercamel, and matchtigs explorations; see below.
13.1 Pandedup-GGCAT integration
Link to heading
- Pandedup now preserves short kmers by manually appending them to the ggcat output file
14.1.1 What does it do and why is it slow?
Link to heading
- Single-threaded currently.
- Iterates over overlap from \(k-1\) down to 0.
- for each, builds a hashset of all prefixes and rc suffixes, and tries to match
them to suffixes and rc prefixes.
- does some tricky stuff with a
next array to handle the fact that multiple
prefixes/rc suffixes can have the same value, and so we need the equivalent of
a C++ unordered_multimap, but khash doesn’t support this directly. - half the time is in hash table inserts/lookups
- half the time is in encoding/decoding the (rc) kmer at the start/end of each input sequence
- There is a
memory_reduction constant of 16 so that the hashmap is only built
for 1/16th of the data at a time. This saves 16x max memory, but makes it (up
to) 16x slower also. In practice my laptop has plenty RAM to just have this
constant at 1 instead.
- Implemented in
pandedup/examples/mss.rs (for masked super string). - Instead of hashing, we can collect all \(4n\) terminal kmers (head, rc tail) and
(tail, rc head) in two vectors and sort those using
voracious_radix_sort
with multiple threads. - Then, do a mergesort, and link matching heads and tails. Fill the
next: Option<(usize, u8)> and prev: Option<(usize, u8)> of the two nodes
accordingly, containing the adjacent node ID and the length of the overlap.
- on laptop; up to 12 threads
- k=63
| method | input | records | input Gbp | output Gbp | time |
|---|
| kmercamel | eulertigs | 27.0M | 9.58 | 7.8 | 5 h |
| fwd | eulertigs | 58.3M | 9.58 | 8.11 | 237 s |
| fwd+rc | eulertigs | 58.3M | 9.58 | 7.88 | 361 s |
| fwd | greedy matchtigs | 27.0M | 8.42 | 8.08 | 117 s |
| fwd+rc | greedy matchtigs | 27.0M | 8.42 | 8.06 | 191 s |
Unsurprisingly, passing greedy matchings is ~2x faster as they have ~2x fewer
input records. But the output is not as small, because the matchtigs introduce
additional bases that we cannot (or at least do not) get rid of anymore.
14.1.4 Reusing the sorted tail array
Link to heading
We can generate and sort the tails only once and then just shift-right
on-demand to get the exact key. We can then either compactify the vec by
removing used tails, or just repeatedly iterate the full array. This reduces
the time for fwd+rc greedy matchtigs to 160s.
14.1.5 Reusing the sorted head array
Link to heading
Between consecutive overlaps, the head array loses the 2 most significant bits
of all values, so it becomes 4 sorted runs. We can merge-sort those instead of
doing a full radix sort. This increases the time to 220s, but note that now
the bulk of the work is single-threaded instead.
35% of the running time is now spend in heads.retain() to filter out
already-used tig ends. Skipping that check gives back 150s.
We can use multithreading for the re-sorting of the heads array: we can split
each of the 4 increasing slices into \(t\) sub-slices (one per thread) and merge
those in parallel. 120s.
Currently, the heads and tails array slowly fill up with kmers corresponding
to already-used ends. During the merge, we can directly compactify them and only
write elements that are still unlinked at the time.
These elements could still be linked via the reverse-complement match. In that
case we simply discard them in the next iteration.
Note that we can’t greedily check for every element whether it’s still active,
as that involves a useless cache-miss for most elements. Thus, we only do this
once an element has a match to another end: there, we anyway already have a
cache-miss, and so we only pay the additional cache-miss for those elements that
are indeed redundant now.
90s
We can split the heads and tails into one chunk per thread and merge those
in parallel. We shard the ends array into a list of \(16\cdot t\) chunks each locked by
their own mutex so we can concurrently write to it.
Entries of heads and tails are set to usize::MAX as soon as the
corresponding end point has been linked, and these are then filtered out in a
final single-threaded retain after the merge. Note that this also uses a
bitmask to indicate used links. 50s
We can also do a multi-threaded retain: first retain sub-slices in parallel,
and then merge those into a single array. 60s
Dropping the used bitmask from this: 55s
Still it’s not great: the main thread spends half its time compacting the
per-thread compacted data. Let’s first count how much space each thread will
need for its compacted output, and then write things directly to the right
location. 60s
Better: we can only call retain when over half of the entries is dead: 30s
14.1.10 Reducing further single-threaded bottlenecks
Link to heading
Now, half the time of the main thread is spent determining where the heads array
has boundaries between the 4 character ranges. Instead of a linear scan, let’s
do binary searches. 27s (Turns out the main thread really isn’t the bottleneck
anymore anyway.)
Now, just reading all sequences and writing an array of their start/end kmer
values is ~30% of the time. Let’s also multithread that: 25s
Note: I’m now testing the eulertigs file with roughly 2x the contigs instead of the greedytigs.
- Base: 30GB
- Drop input data as soon as we read their head and starts: 24GB
- Single scratch buffer for head and tail: 23GB
- Drop the 4
u128 values from the ends array; that now only stores the links
from the head and tail of each node: 18.6GB - Store input data in a
PackedSeqVec together with ranges pointing into it: 12.8GB- Now the peak is reached when the scratch buffer is allocated for
retain
- We can shrink our types a bit and use
u32 indices, to get 11.9GB - The heads and tails are still
(u128, u32, bool). Aligning them to u32
reduces memory to 9.2GB and also reduces time from 54s to 44s. - 1-byte alignment further reduces memory to 8.3GB
- intermezzo: Making the algorithm generic over
u128 and u64 for \(k\leq 32\), we reduce
memory further to 6.3GB, and time for \(k=31\) from 26s to 21s.
- Let’s pack the
index: u32 and reverse: bool together in a single u32,
and also use packed/unaligned Link types: 7.8GB (compared to 8.3 above for
k=64) and 42s.
| type | k | input GB | contigs | walltime | CPU time | output Gbp |
|---|
| eulertigs | 17 | 6.9 | | 8 | 188 | 1.96 |
| eulertigs | 31 | 6.4 | | 19 | 113 | 4.39 |
| eulertigs | 63 | 11 | | 44 | 305 | 7.83 |
| greedytigs | 63 | 8.7 | | 25 | 152 | 8.06 |
1
| ggcat eulertigs: 655.37s wall, 7589.32s CPU, 1158.0% average CPU, 6904923140 bytes (hprcv2.dedup.ggcat-eulertigs.fa)
|
- do k=64 and k=32 instead of 63 and 31?
- test for N/ambiguous chars in input and panic on them
- end-to-end benchmarks
- construct actual output sequence with masking => cleanup/check the slop
- add MaskedPackedSeq and use it in Deacon
- compare compression schemes
- parallel cycle reconstruction
- test paraseq for IO
- add
pandedup mss -j
Eulertigs output size: 7833783806
Greedytigs output size: ?? 067
- I should try inputting the bcaml graph rather than just the fasta: https://github.com/algbio/matchtigs/issues/8
- Currently loading input (and presumably reconstruction the graph) takes 40min on 1 thread
- Dijkstra takes 16min on many threads
- Instead of dijkstra with a binary heap, a bucket queue of size k should be enough
- let’s try to parallellize things; the initial graph construction also seems single-threaded
15.1.1 Building the unitig graph
Link to heading
Building unitigs as input:
1
2
3
4
5
6
| > cargo run -r --example ggcat hprcv2.dedup.fa.zst -k 63 --mode none
11:32:18 INFO Read hprcv2.dedup.fa.zst: 3941236558 bytes, 59591914 records, 14855772704 bases, 249.3 bases/record
11:32:20 INFO ggcat input: 59584907 in-memory contigs, 14855772704 bases (from hprcv2.dedup.fa.zst)
11:38:03 INFO Wrote hprcv2.dedup-k63-none.fa: 21181450134 bytes, records unknown, bases unknown, bases/record unknown
11:38:03 INFO ggcat: 381.23s wall [2054.63s CPU; 539.0%], input 3941236558 bytes, output 21181450134 bytes (hprcv2-unitigs.fa)
|
- 3.9Gbp in => 21.2 Gbp out over 186M contigs of average length 114bp, but
median length probably close to 64
Time for some profiling:
1
| cargo build -r && samply record ./target/release/matchtigs -k 63 -t 6 --fa-in $pd/hprcv2-k63-unitigs.fa --greedytigs-fa-out $pd/hprcv2-k63-greedytigs-2.fa
|
A lot of time is spent on parsing the input and building the graph. As
recommended by the docs, let’s instead include graph information in the fasta
headers in BCALM format. We use the GGCAT API, but this is the -e option in
the CLI.
Building unitigs with links:
1
2
3
4
5
6
| > cargo run -r -- ggcat hprcv2.dedup.fa.zst -k 63 --mode unitig-links
14:17:59 INFO Read hprcv2.dedup.fa.zst: 3941236558 bytes, 59591914 records, 14855772704 bases, 249.3 bases/record
14:18:00 INFO ggcat input: 59584907 in-memory contigs, 14855772704 bases (from hprcv2.dedup.fa.zst)
14:26:11 INFO Wrote hprcv2.dedup-k63-unitig-links.fa: 29386798229 bytes, records unknown, bases unknown, bases/record unknown
14:26:11 INFO ggcat: 529.31s wall [2598.50s CPU; 490.9%], input 3941236558 bytes, output 29386798229 bytes (hprcv2.dedup-k63-unitig-links.fa)
|
unitig graph:
- bases: 21.2Gbp
- unitigs: 186M
- edges: 267M
We can try running matchtigs again, but it crashes after ~12 mins with
SIGKILL, so I assume it went out of memory:
1
| matchtigs -k 63 -t 6 --bcalm-in $pd/hprcv2-k63-unitig-links.fa --greedytigs-fa-out $pd/hprcv2-k63-greedytigs-2.fa
|
Investigating the unitigs, we see this:
1
2
3
4
5
6
7
8
9
10
11
12
| >10097112 LN:i:70 L:+:10097111:+ L:-:10097113:- L:-:10097114:-
AAGAAAGAAAGAAAGAAAGAAAGAAATAATAAAAAGAAAGAAAAAGAAAGAAAGAAAGAAAAAGAAAGAA
>10097113 LN:i:67 L:+:10097112:+ L:-:10097117:-
GAAAGAAGAAAGAAAGAAAGAAAGAAAGAAATAATAAAAAGAAAGAAAAAGAAAGAAAGAAAGAAAA
>10097114 LN:i:67 L:+:10097112:+ L:-:135476963:-
GAAGAAAGAAAGAAAGAAAGAAAGAAAGAAATAATAAAAAGAAAGAAAAAGAAAGAAAGAAAGAAAA
>10097115 LN:i:67 L:+:137890227:+ L:-:10097117:-
GAAAGAAGAAAGAAAGAAAGAAAGAAAGAAATAATAAAAAGAAAGAAAAAGAAAGAAAGAAAAAGAA
>10097116 LN:i:67 L:+:137890227:+ L:-:135476963:-
GAAGAAAGAAAGAAAGAAAGAAAGAAAGAAATAATAAAAAGAAAGAAAAAGAAAGAAAGAAAAAGAA
>10097117 LN:i:76 L:+:10097113:+ L:+:10097115:+ L:-:135471706:-
AAGAAAAAGAAAAAGAAAGAAGAAAGAAAGAAAGAAAGAAAGAAATAATAAAAAGAAAGAAAAAGAAAGAAAGAAA
|
In particular, nodes 10097113 to 10097116 all have 1 incoming (L:-:...) and 1
outgoing (L:+:...) edge. Still, they cannot be merged because supposedly there
are other edges joining in those nodes.
15.1.4 Graph representation: Unitig-centric
Link to heading
The first option is to number the unitigs, as in the snippet above, and for each
unitig store the IDs (and strand) of al connecting unitigs.
When a \(k-1\)-mer occurs at the start of 4 unitigs and end of 4 unitigs, that
means we get a total of 16 edges connecting them, so that we store 32 IDs.
15.1.5 Graph representation: overlap-centric
Link to heading
Alternatively, we could ID the \(k-1\)-mers at the start and end of every
unitig. I estimate their count is roughly 1/2 to 2/3 of the number of unitigs,
corresponding respectively to average degree 4 and 3.
In this case, each unitig becomes an edge between two \(k-1\)-mers and thus
corresponds to storing 2 \(k-1\)-mer IDs for the forward and reverse direction.
I suspect this this turns out more compact overall.
Most unitigs are only slightly longer than \(k\), and we could store not the
target ID, but the bases that are added to the kmer.
Together with a PHF mapping kmers to IDs, that could maybe give a compact representation?
Here’s the prompt:
In src/mss.rs there is an algorithm for creating a superstring (spectrum preserving string set) of a fasta file input. Follow that overall style.
We’ll add a new algorithm in src/matchtigs.rs.
The input is again a PackedSeqVec and a list of ranges. This time, we assume the input is unitigs (whereas before simplitigs also worked). First we must recover the graph. For each unitig, add (head, id) and (rc tail, id) to a list of heads, and (tail, id) and (rc head, id) to a list of tails. Then, sort both lists. Then, build the graph: for each range of heads and tails with the same kmer, add an edge from each head id to each tail id. Do this by simply pushing (from id, to id) on a new list. Then sort that list, and build a graph out of it. The representation is one list of edges where the parity of the to_id encodes from which side we enter the unitid id. Then we keep a separate list with for every unitig two entries: the u32 index where its outgoing edges start.
Once the graph representation is done, we start the algorithm: we again have a list of endpoints/tails that will become the “receiving ends” from the left, and we have a list of tails from where we grow sequences to the right. This time, we store the starting tail id, and the current tail id (which is initially is the same).
Then do a first “extend edge” step, where for each unitig tail, we replace it using the graph by all outgoing edges, specifically the tail id corresponding to each of those unitigs. Collect them in a vec and then sort the vec. Then we do a linear pass to link paths to remaining active receiving tails. We have a separate Vec<Links> to store active links: whenever a new match is found, check if the two corresponding unitig-ends are still active; if not, then mark these nodes as u32::MAX. Otherwise, link the two unitig ends (the starting node of the path to the receiving unitig).
Then filter out all done unitigs, and repeat: extend all active paths by 1 edge, sort tail ids, find matches, and shrink the 2 vectors.
All edges also have a weight though, corresponding to the number of 1-character extensions to go from the end of the current unitig to the end of the next unitig.
One thing is tricky though: we don’t want to include the full cost of the last unitig. So I think we should first run the matching stage, and then update the cost of all edges. We want to use a bucket queue: we keep a separate vec for the edges at each distance up to k. In each iteration, we only work on the first non-empty vec
- packed structs
- drop intermediate path parent pointers, since just final links is enough
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
| > cargo run -r -- matchtigs hprcv2-k63-unitigs.fa.zst -k 63
02:09:28 INFO Reading unitigs..
02:10:21 INFO Read hprcv2-k63-unitigs.fa.zst: 4914797023 bytes, 186730681 records, 17488538356 bases, 93.7 bases/record
02:10:21 INFO Building graph of 186723674 unitigs..
02:10:33 INFO Sorting heads
02:11:00 INFO Sorting tails
02:11:28 INFO Matching unitig ends
02:11:34 INFO Sorting 534216886 edges
02:12:00 INFO Building CSR adjacency
02:12:00 INFO Filling offsets: 1424 MB
02:12:02 INFO Mapping pairs
02:12:03 INFO Allocating 373447348 links: 1780 MB
02:12:03 INFO Allocating 373447348 receiving ids: 1424 MB
Pushing to initial bucket
02:12:04 INFO Distance 0: 373447348 active paths, 373447348 receiving ends
02:12:52 INFO Distance 1: 22900904 active paths, 118578822 receiving ends
02:12:56 INFO Distance 2: 22658491 active paths, 112868300 receiving ends
02:16:14 INFO Distance 16: 317083268 active paths, 84779338 receiving ends
|
using ~32GB at this point
This crashes with out-of-memory eventually in a separate run:
1
2
| 02:38:42 INFO After distance 25, bucket queue sizes: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 887547819, 362200595, 162538168, 124284752, 86530464, 69234377, 58137007, 50224293, 43420318, 40512191, 37494701, 35084655, 32624626, 31644801, 33588580, 33041442, 32333127, 30225164, 30976702, 31286851, 32149128, 31869221, 33070349, 32540403, 35147195, 31435462, 31503491, 28257588, 29517232, 27807232, 27708123, 23382365, 20760981, 17088280, 18330846, 16236818, 0]
02:38:42 INFO Distance 26: 887547819 active paths, 75642574 receiving ends
|
at this point there are 2.6G queue elements in total, for 186M unitigs…
- voracious radix sort brings the running time to get to distance 20 excluding input reading from 8.5min to 3.5min.
I think: for each distance, iterate over all active unitigs, and for each do a
full DFS/BFS backwards search to find a match. Kinda weird but it does seem to
work, although slowly. Memory usage is around 18GB.
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
| 13:24:06 INFO Reading unitigs..
13:24:58 INFO Read hprcv2-k63-unitigs.fa.zst: 4914797023 bytes, 186730681 records, 17488538356 bases, 93.7 bases/record
13:24:58 INFO Building graph of 186723674 unitigs..
13:25:09 INFO Sorting heads
13:25:15 INFO Sorting tails
13:25:21 INFO Matching unitig ends
13:25:27 INFO Sorting 534216886 edges
13:25:29 INFO Building CSR adjacency
13:25:29 INFO Filling offsets: 1424 MB
13:25:31 INFO Mapping pairs
13:25:33 INFO Distance 0: 373447348 receiving ends
13:29:03 INFO After distance 0: 127434263 links
13:29:03 INFO Distance 1: 7997244 receiving ends
13:29:17 INFO After distance 1: 130289575 links
13:29:17 INFO Distance 2: 7046483 receiving ends
13:48:16 INFO Distance 60: 1422199 receiving ends
13:48:32 INFO After distance 60: 157940371 links
13:48:32 INFO Distance 61: 1425403 receiving ends
13:48:50 INFO After distance 61: 158111508 links
13:48:50 INFO Distance 62: 1423215 receiving ends
13:49:07 INFO After distance 62: 158281419 links
13:49:08 INFO Reconstruct output
13:51:30 INFO Wrote hprcv2-k63-unitigs.fa.matchtigs.msfa: 8257821647 bytes, 1 records, 8257821616 bases, 8257821616.0 bases/record
13:51:30 INFO matchtigs: 1643.62s wall [1783.16s CPU; 108.5%], output 8257821616 bases (hprcv2-k63-unitigs.fa.matchtigs.msfa)
|
Output:
- 8.26Gbp or 8.3GB; actually an improvement over the 8.7GB greedytigs!
- 1 record
- 1640s walltime (1780s CPU).
Again, for comparison (with tool):
- unitigs (ggcat): 17.5Gbp
- eulertigs (ggcat): 10.7Gbp
- greedy matchtigs (matchtigs): 8.71Gbp
- greedy matchtigs (us): 8.26Gbp << new here
- masked superstring (us): 7.83Gbp
- Rough current algorithm:
- iterate distance upwards
- for each dist, do multiple passes
- within a pass, iterate over all nodes and find those at the current distance
(the implicit frontier)
- for each that has an outgoing edge to a target that wasn’t used yet, find a
“receiver” (weird/opposite naming?) ie the source of a path of the given
length that ends here in the right distance, by doing a DFS over the
spanning tree.
Output (aborted):
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
| cargo run -r -- matchtigs hprcv2-k63-unitigs.fa.zst -k 63
13:53:40 INFO Reading unitigs..
13:54:33 INFO Read hprcv2-k63-unitigs.fa.zst: 4914797023 bytes, 186730681 records, 17488538356 bases, 93.7 bases/record
13:54:34 INFO Building graph of 186723674 unitigs..
13:54:46 INFO Sorting heads
13:54:52 INFO Sorting tails
13:54:57 INFO Matching unitig ends
13:55:03 INFO Sorting 534216886 edges
13:55:06 INFO Building CSR adjacency
13:55:06 INFO Filling offsets: 1424 MB
13:55:08 INFO Mapping pairs
13:56:21 INFO Distance 0, pass 1: 373447348 frontier nodes, 119690289 links
13:58:31 INFO Distance 0, pass 2: 134066770 frontier nodes, 6513605 links
13:58:58 INFO Distance 0, pass 3: 121039560 frontier nodes, 949707 links
13:59:11 INFO Distance 0, pass 4: 119140146 frontier nodes, 280662 links
13:59:20 INFO Distance 0, pass 5: 118578822 frontier nodes, 0 links
13:59:20 INFO After distance 0: 127434263 links
13:59:22 INFO Distance 1, pass 1: 18940594 frontier nodes, 2798292 links
13:59:37 INFO Distance 1, pass 2: 13912819 frontier nodes, 56205 links
13:59:39 INFO Distance 1, pass 3: 13809173 frontier nodes, 764 links
13:59:40 INFO Distance 1, pass 4: 13807697 frontier nodes, 0 links
13:59:40 INFO After distance 1: 130289524 links
13:59:42 INFO Distance 2, pass 1: 12593328 frontier nodes, 2099130 links
13:59:57 INFO Distance 2, pass 2: 8862016 frontier nodes, 65518 links
13:59:59 INFO Distance 2, pass 3: 8725709 frontier nodes, 4509 links
14:00:00 INFO Distance 2, pass 4: 8711163 frontier nodes, 327 links
14:00:01 INFO Distance 2, pass 5: 8709765 frontier nodes, 24 links
14:00:02 INFO Distance 2, pass 6: 8709644 frontier nodes, 1 links
14:00:03 INFO Distance 2, pass 7: 8709636 frontier nodes, 0 links
14:00:03 INFO After distance 2: 132459033 links
14:09:22 INFO After distance 16: 144867407 links
^C
|
Currently quite some time is wasted in repeated passes that iterate over all
nodes again and again, even though only few changed. Instead we can just update
the graph directly after finding a new link and do a local re-check.
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
| > cargo run -r -- matchtigs hprcv2-k63-unitigs.fa.zst -k 63 1230.515s 14:14
14:55:24 INFO Reading unitigs..
14:56:16 INFO Read hprcv2-k63-unitigs.fa.zst: 4914797023 bytes, 186730681 records, 17488538356 bases, 93.7 bases/record
14:56:18 INFO Building graph of 186723674 unitigs..
14:56:29 INFO Sorting heads
14:56:35 INFO Sorting tails
14:56:41 INFO Matching unitig ends
14:56:47 INFO Sorting 534216886 edges
14:56:49 INFO Building CSR adjacency
14:56:49 INFO Filling offsets: 1424 MB
14:56:51 INFO Mapping pairs
14:56:54 INFO Matching 186723674 unitigs, 373447348 ends, k 63: two ends per unitig; overlap = k minus one minus distance; examined = receiver checks; links made = connections; receiving ends remain = unlinked ends scheduled below k; free ends = initial ends minus twice cumulative links
14:56:54 INFO Before distance 0: 0 receiving ends examined, 0 links made, 373447348 receiving ends remain, 373447348 free ends, 17488538356 estimated bases
14:56:54 INFO Node visits count repeated inspections: direct_overlap scans graph neighbors at the initial distance; nearest_receiver includes its reverse search; connect touches two link slots; remove_source includes tree, boundary, and repair scans
14:56:54 INFO Estimated bases starts as the sum of unitig lengths and drops by every link overlap; cycle closures can make it lower than the final linear output
14:59:27 INFO After distance 0: 246013085 receiving ends examined, 127434263 links made, 118578822 receiving ends remain, 118578822 free ends, 9587614050 estimated bases
14:59:27 INFO After distance 0 node visits: direct_overlap 257853327 (total 257853327), nearest_receiver 0 (total 0), connect 254868526 (total 254868526), remove_source 9061551581 (total 9061551581)
14:59:52 INFO After distance 1: 115723510 receiving ends examined, 2855312 links made, 77167699 receiving ends remain, 112868198 free ends, 9413440018 estimated bases
14:59:52 INFO After distance 1 node visits: direct_overlap 0 (total 257853327), nearest_receiver 117871611 (total 117871611), connect 5710624 (total 260579150), remove_source 1493213306 (total 10554764887)
15:00:09 INFO After distance 2: 3028201 receiving ends examined, 2169644 links made, 72820395 receiving ends remain, 108528910 free ends, 9283261378 estimated bases
15:00:09 INFO After distance 2 node visits: direct_overlap 0 (total 257853327), nearest_receiver 3356136 (total 121227747), connect 4339288 (total 264918438), remove_source 1659115818 (total 12213880705)
15:28:46 INFO After distance 62: 1158960 receiving ends examined, 169648 links made, 0 receiving ends remain, 56887566 free ends, 8255130536 estimated bases
15:28:46 INFO After distance 62 node visits: direct_overlap 0 (total 257853327), nearest_receiver 1186320 (total 230476753), connect 339296 (total 316559782), remove_source 1729083327 (total 172785885593)
15:28:47 INFO Final node visits (including closure links): direct_overlap 257853327, nearest_receiver 230476753, connect 373447348, remove_source 172785885593
15:28:47 INFO Reconstruct output
15:31:10 INFO Length summary: 8255130536 estimated bases after links, 8257861763 output bases
15:31:14 INFO Wrote hprcv2-k63-unitigs.fa.matchtigs.msfa: 8257861794 bytes, 1 records, 8257861763 bases, 8257861763.0 bases/record
15:31:14 INFO matchtigs: 2149.84s wall [2284.33s CPU; 106.3%], output 8257861763 bases (hprcv2-k63-unitigs.fa.matchtigs.msfa)
|
New output: 2149s (2284s CPU), 8.26Gbp, ie same as before.
- I’m a bit surprised this is actually slower than before.
Every time a link is made, we have to repair a part of the graph around it.
There were 2 DFSes and one loop before; instead let’s just do a single DFS.
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
| 15:35:31 INFO Reading unitigs..
15:36:22 INFO Read hprcv2-k63-unitigs.fa.zst: 4914797023 bytes, 186730681 records, 17488538356 bases, 93.7 bases/record
15:36:24 INFO Building graph of 187M unitigs..
15:36:35 INFO Sorting heads
15:36:41 INFO Sorting tails
15:36:47 INFO Matching unitig ends
15:36:52 INFO Sorting 534M edges
15:36:55 INFO Building CSR adjacency
15:36:55 INFO Filling offsets: 1.49GB
15:36:56 INFO Mapping pairs
15:36:59 INFO Matching 187M unitigs, 373M ends, k 63: two ends per unitig; overlap = k minus one minus distance; examined = receiver checks; links made = connections; receiving ends remain = unlinked ends scheduled below k; free ends = initial ends minus twice cumulative links
15:36:59 INFO Before distance 0: 0 receiving ends examined, 0 links made, 373M receiving ends remain, 373M free ends, 17.5G estimated bases
15:36:59 INFO Node visits count repeated inspections: direct_overlap scans graph neighbors at the initial distance; nearest_receiver includes its reverse search; connect touches two link slots; remove_source includes tree, boundary, and repair scans
15:36:59 INFO Estimated bases starts as the sum of unitig lengths and drops by every link overlap; cycle closures can make it lower than the final linear output
15:39:01 INFO After distance 0: 246M receiving ends examined, 127M links made, 119M receiving ends remain, 119M free ends, 9.59G estimated bases
15:39:01 INFO After distance 0 node visits: direct_overlap 258M (total 258M), nearest_receiver 0 (total 0), connect 255M (total 255M), remove_source 3.50G (total 3.50G)
15:39:20 INFO After distance 1: 116M receiving ends examined, 2.86M links made, 77.2M receiving ends remain, 113M free ends, 9.41G estimated bases
15:39:20 INFO After distance 1 node visits: direct_overlap 0 (total 258M), nearest_receiver 118M (total 118M), connect 5.71M (total 261M), remove_source 305M (total 3.80G)
15:39:32 INFO After distance 2: 3.03M receiving ends examined, 2.17M links made, 72.8M receiving ends remain, 109M free ends, 9.28G estimated bases
15:39:32 INFO After distance 2 node visits: direct_overlap 0 (total 258M), nearest_receiver 3.35M (total 121M), connect 4.34M (total 265M), remove_source 311M (total 4.11G)
15:57:27 INFO After distance 62: 1.16M receiving ends examined, 170k links made, 0 receiving ends remain, 56.9M free ends, 8.26G estimated bases
15:57:27 INFO After distance 62 node visits: direct_overlap 0 (total 258M), nearest_receiver 1.19M (total 230M), connect 339k (total 317M), remove_source 233M (total 26.0G)
15:57:28 INFO Final node visits (including closure links): direct_overlap 258M, nearest_receiver 230M, connect 373M, remove_source 26.0G
15:57:28 INFO Reconstruct output
15:59:49 INFO Length summary: 8.26G estimated bases after links, 8.26G output bases
15:59:54 INFO Wrote hprcv2-k63-unitigs.fa.matchtigs.msfa: 8257965107 bytes, 1 records, 8257965076 bases, 8257965076.0 bases/record
15:59:54 INFO matchtigs: 1462.87s wall [1467.26s CPU; 100.3%; MaxRSS 26.18 GiB], output 8.26G bases (hprcv2-k63-unitigs.fa.matchtigs.msfa)
|
- 1462s (1467s CPU); 26GB memory
The DFS can simply infer the owner trees on-the-fly; let’s save some memory and
not store these things explicitly.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
| 16:27:38 INFO Reading unitigs..
16:28:29 INFO Read hprcv2-k63-unitigs.fa.zst: 4914797023 bytes, 186730681 records, 17488538356 bases, 93.7 bases/record
16:28:31 INFO Building graph of 187M unitigs..
16:28:42 INFO Sorting heads
16:28:49 INFO Sorting tails
16:28:55 INFO Matching unitig ends
16:29:01 INFO Sorting 534M edges
16:29:04 INFO Building CSR adjacency
16:29:04 INFO Filling offsets: 1.49GB
16:29:05 INFO Mapping pairs
16:29:07 INFO Matching 187M unitigs, 373M ends, k 63: two ends per unitig; overlap = k minus one minus distance; examined = receiver checks; links made = connections; receiving ends remain = unlinked ends scheduled below k; free ends = initial ends minus twice cumulative links
16:29:07 INFO Before distance 0: 0 receiving ends examined, 0 links made, 373M receiving ends remain, 373M free ends, 17.5G estimated bases
16:29:07 INFO Node visits count repeated inspections: direct_overlap scans graph neighbors at the initial distance; nearest_receiver includes its reverse search; connect touches two link slots; remove_source includes tree, boundary, and repair scans
16:30:51 INFO After distance 0: 246M receiving ends examined, 127M links made, 119M receiving ends remain, 119M free ends, 9.59G estimated bases
16:30:51 INFO After distance 0 node visits: direct_overlap 258M (total 258M), nearest_receiver 0 (total 0), connect 255M (total 255M), remove_source 3.57G (total 3.57G)
16:32:04 INFO After distance 1: 116M receiving ends examined, 2.86M links made, 77.2M receiving ends remain, 113M free ends, 9.41G estimated bases
16:32:04 INFO After distance 1 node visits: direct_overlap 0 (total 258M), nearest_receiver 752M (total 752M), connect 5.71M (total 261M), remove_source 328M (total 3.90G)
16:32:16 INFO After distance 2: 3.03M receiving ends examined, 2.17M links made, 72.8M receiving ends remain, 109M free ends, 9.28G estimated bases
16:32:16 INFO After distance 2 node visits: direct_overlap 0 (total 258M), nearest_receiver 29.7M (total 782M), connect 4.34M (total 265M), remove_source 353M (total 4.25G)
16:35:44 INFO After distance 16: 2.02M receiving ends examined, 533k links made, 46.7M receiving ends remain, 83.7M free ends, 8.61G estimated bases
16:35:44 INFO After distance 16 node visits: direct_overlap 0 (total 258M), nearest_receiver 52.2M (total 1.54G), connect 1.07M (total 290M), remove_source 383M (total 9.07G)
^C
|
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
| 16:36:53 INFO Reading unitigs..
16:37:44 INFO Read hprcv2-k63-unitigs.fa.zst: 4914797023 bytes, 186730681 records, 17488538356 bases, 93.7 bases/record
16:37:45 INFO Building graph of 187M unitigs..
16:37:56 INFO Sorting heads
16:38:03 INFO Sorting tails
16:38:09 INFO Matching unitig ends
16:38:15 INFO Sorting 534M edges
16:38:18 INFO Building CSR adjacency
16:38:18 INFO Filling offsets: 1.49GB
16:38:20 INFO Mapping pairs
16:38:21 INFO Matching 187M unitigs, 373M ends, k 63: two ends per unitig; overlap = k minus one minus distance; examined = receiver checks; links made = connections; receiving ends remain = unlinked ends scheduled below k; free ends = initial ends minus twice cumulative links
16:38:21 INFO Before distance 0: 0 receiving ends examined, 0 links made, 373M receiving ends remain, 373M free ends, 17.5G estimated bases
16:38:21 INFO Node visits count repeated inspections: direct_overlap scans graph neighbors at the initial distance; nearest_receiver includes its reverse search; connect touches two link slots; remove_source scans changed nodes and edges used to recompute distances
16:39:41 INFO After distance 0: 246M receiving ends examined, 127M links made, 119M receiving ends remain, 119M free ends, 9.59G estimated bases
16:39:41 INFO After distance 0 node visits: direct_overlap 258M (total 258M), nearest_receiver 0 (total 0), connect 255M (total 255M), remove_source 2.09G (total 2.09G)
16:40:51 INFO After distance 1: 116M receiving ends examined, 2.86M links made, 77.2M receiving ends remain, 113M free ends, 9.41G estimated bases
16:40:51 INFO After distance 1 node visits: direct_overlap 0 (total 258M), nearest_receiver 752M (total 752M), connect 5.71M (total 261M), remove_source 178M (total 2.27G)
16:41:00 INFO After distance 2: 3.03M receiving ends examined, 2.17M links made, 72.8M receiving ends remain, 109M free ends, 9.28G estimated bases
16:41:00 INFO After distance 2 node visits: direct_overlap 0 (total 258M), nearest_receiver 29.7M (total 782M), connect 4.34M (total 265M), remove_source 168M (total 2.44G)
16:48:22 INFO After distance 62: 1.16M receiving ends examined, 170k links made, 0 receiving ends remain, 56.9M free ends, 8.26G estimated bases
16:48:22 INFO After distance 62 node visits: direct_overlap 0 (total 258M), nearest_receiver 4.26M (total 3.56G), connect 339k (total 317M), remove_source 57.2M (total 7.79G)
16:48:23 INFO Final node visits (including closure links): direct_overlap 258M, nearest_receiver 3.56G, connect 373M, remove_source 7.79G
16:48:23 INFO Reconstruct output
16:50:44 INFO Length summary: 8.26G estimated bases after links, 8.26G output bases
16:50:49 INFO Wrote hprcv2-k63-unitigs.fa.matchtigs.msfa: 8257965169 bytes, 1 records, 8257965138 bases, 8257965138.0 bases/record
16:50:49 INFO matchtigs: 836.30s wall [973.04s CPU; 116.4%; MaxRSS 26.18 GiB], output 8.26G bases (hprcv2-k63-unitigs.fa.matchtigs.msfa)
|
Now we have 836s! (and still 26GB mem)
16.1.9 Lazy remove_source: 2-3x slower
Link to heading
Instead of eagerly updating the tree on node-removal, we can lazily update nodes
when we have a candidate match at the right distance.
This turns out much slower. I’m not exactly sure why, but probably because
(somehow?) nodes get relazed multiple times each time a node disappears.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
| 17:09:59 INFO Reading unitigs..
17:10:51 INFO Read hprcv2-k63-unitigs.fa.zst: 4914797023 bytes, 186730681 records, 17488538356 bases, 93.7 bases/record
17:10:52 INFO Building graph of 187M unitigs..
17:11:04 INFO Sorting heads
17:11:10 INFO Sorting tails
17:11:16 INFO Matching unitig ends
17:11:22 INFO Sorting 534M edges
17:11:24 INFO Building CSR adjacency
17:11:24 INFO Filling offsets: 1.49GB
17:11:26 INFO Mapping pairs
17:11:27 INFO Matching 187M unitigs, 373M ends, k 63: two ends per unitig; overlap = k minus one minus distance; examined = receiver checks; links made = connections; receiving ends remain = unlinked ends scheduled below k; free ends = initial ends minus twice cumulative links
17:11:27 INFO Before distance 0: 0 receiving ends examined, 0 links made, 373M receiving ends remain, 373M free ends, 17.5G estimated bases
17:11:27 INFO Node visits count repeated inspections: direct_overlap scans graph neighbors at the initial distance; nearest_receiver includes its bounded DFS; connect touches two link slots
17:11:38 INFO After distance 0: 246M receiving ends examined, 127M links made, 119M receiving ends remain, 119M free ends, 9.59G estimated bases
17:11:38 INFO After distance 0 node visits: direct_overlap 258M (total 258M), nearest_receiver 0 (total 0), connect 255M (total 255M)
17:12:00 INFO After distance 1: 116M receiving ends examined, 2.86M links made, 113M receiving ends remain, 113M free ends, 9.41G estimated bases
17:12:00 INFO After distance 1 node visits: direct_overlap 0 (total 258M), nearest_receiver 291M (total 291M), connect 5.71M (total 261M)
17:12:23 INFO After distance 2: 111M receiving ends examined, 2.17M links made, 109M receiving ends remain, 109M free ends, 9.28G estimated bases
17:12:23 INFO After distance 2 node visits: direct_overlap 0 (total 258M), nearest_receiver 302M (total 592M), connect 4.34M (total 265M)
17:17:42 INFO After distance 13: 87.7M receiving ends examined, 620k links made, 87.0M receiving ends remain, 87.1M free ends, 8.69G estimated bases
17:17:42 INFO After distance 13 node visits: direct_overlap 0 (total 258M), nearest_receiver 327M (total 4.08G), connect 1.24M (total 286M)
<aborted>
|
16.1.10 Not quite meet-in-the-middle
Link to heading
Let’s alternate relaxing from the start and from the end. That way, exploration
trees should be much smaller.
Here, GPT did a bunch of stuff that is not exactly what I intended I think, but
it seems to work.
My idea (that’s not really what is happening now; see below): every time we relax an edge and detect that we already have a distance
to the other end, push the edge on a bucket queue for that total distance.
Then we only need to relax edges from nodes at distance at most \(k/2\) or so, and
process the ‘bridges’/meeting points up to distance \(k\).
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
| 18:21:08 INFO Reading unitigs..
18:21:59 INFO Read hprcv2-k63-unitigs.fa.zst: 4914797023 bytes, 186730681 records, 17488538356 bases, 93.7 bases/record
18:22:01 INFO Building graph of 187M unitigs..
18:22:12 INFO Sorting heads
18:22:18 INFO Sorting tails
18:22:24 INFO Matching unitig ends
18:22:30 INFO Sorting 534M edges
18:22:33 INFO Building CSR adjacency
18:22:33 INFO Filling offsets: 1.49GB
18:22:34 INFO Mapping pairs
18:22:36 INFO Matching 187M unitigs, 373M ends, k 63: two ends per unitig; overlap = k minus one minus distance; examined = bridge entries checked; links made = connections; receiving ends remain and free ends = initial ends minus twice cumulative links
18:22:36 INFO Meeting search: source and receiver fronts each store distances below 31; a path through one unitig has cost left distance + unitig extensions + right distance
18:22:36 INFO Before distance 0: 0 receiving ends examined, 0 links made, 373M receiving ends remain, 373M free ends, 17.5G estimated bases
18:22:36 INFO Node visits count repeated inspections: direct_overlap scans graph neighbors at distance zero; meet checks and schedules central unitig traversals and traces endpoints; connect touches two link slots; remove_source repairs both bounded fronts
18:23:45 INFO After distance 0: 246M receiving ends examined, 127M links made, 119M receiving ends remain, 119M free ends, 9.59G estimated bases
18:23:45 INFO After distance 0 node visits: direct_overlap 258M (total 258M), meet 0 (total 0), connect 255M (total 255M), remove_source 1.95G (total 1.95G)
18:24:08 INFO After distance 1: 8.58M bridge entries examined, 2.86M links made, 113M receiving ends remain, 113M free ends, 9.41G estimated bases
18:24:08 INFO After distance 1 node visits: direct_overlap 0 (total 258M), meet 77.7M (total 1.33G), connect 5.71M (total 261M), remove_source 161M (total 2.11G)
18:24:14 INFO After distance 2: 10.2M bridge entries examined, 2.17M links made, 109M receiving ends remain, 109M free ends, 9.28G estimated bases
18:24:14 INFO After distance 2 node visits: direct_overlap 0 (total 258M), meet 83.6M (total 1.42G), connect 4.34M (total 265M), remove_source 150M (total 2.26G)
18:27:17 INFO After distance 61: 1.54M bridge entries examined, 172k links made, 57.2M receiving ends remain, 57.2M free ends, 8.26G estimated bases
18:27:17 INFO After distance 61 node visits: direct_overlap 0 (total 258M), meet 9.78M (total 4.01G), connect 343k (total 316M), remove_source 12.4M (total 4.96G)
18:27:19 INFO Final node visits (including closure links): direct_overlap 258M, meet 4.01G, connect 373M, remove_source 4.96G
18:27:19 INFO Reconstruct output
18:29:39 INFO Length summary: 8.26G estimated bases after links, 8.26G output bases
18:29:44 INFO Wrote hprcv2-k63-unitigs.fa.matchtigs.msfa: 8258031942 bytes, 1 records, 8258031911 bases, 8258031911.0 bases/record
18:29:44 INFO matchtigs: 515.75s wall [659.01s CPU; 127.8%; MaxRSS 26.19 GiB], output 8.26G bases (hprcv2-k63-unitigs.fa.matchtigs.msfa)
|
At this point, our main algorithm is only like half the total time!
- 60s reading input
- 100s prep; building graph; linking $(k-1)$-overlaps
- 200s: Our actual algorithm
- 140s writing output
Let’s try to understand the details:
- First link adjacent unitigs
- Then compute distances from free ends to remaining unitigs ends
- Then push all ‘bridge’ unitigs on a bucket queue corresponding to the length
of the path they complete.
- Process the lowest-cost unitigs, drop their endpoints, and update
corresponding distances.
- If the cost of a unitig changed, re-push it at the new distance so we try
again later.
Some small optimizations:
- To compute the weight of a bridge (the middle edge), we can just take
dist[head]+dist[tail]-cost(unitig). (Minus, instead of plus, because both
paths to the middle already include the cost of the unitig iteslf.) The only
times we cannot do this is when one of the costs is 0 (because it is a
starting point) or infinity. - There was quite a lot of code to avoid self-loops where a unitig-end connects
back onto itself from the other direction. This stuff should be quite rare, so
let’s just allow that and detect this during string reconstruction, to save
250 lines of code.
- The initialization of all distances was weird: they started all at 0 (because
everything is an open end) and then in the first iteration many of them were
removed and trees updated. Now we only write those elements that remain after
the first pass and do a multi-source dijkstra from those, updating everything
at once.
PackedSeq::unpack_into() nearly doubles the speed of string reconstruction, but it
still pretty slow at 80s (from 140s).
At this point, the total time is 420s.
16.1.11 Multithreading (failed)
Link to heading
I tried a few ideas, but none really work yet:
- We could reconstruct in parallel, but most time ends in 1 thread (that has 1
long cycles?)
- We could check for good bridges in parallel, but a lot of the work is in
dynamically updating the graph.
16.1.12 Reducing bridges (failed)
Link to heading
- Currently we always have most unitigs as a bridge in a queue somewhere. In
practice, when there is a path from A to B, we only need to consider the
middle edge.
- We then need to dynamically push new bridges whenever
remove_source leads to
new ‘middle’ edges appearing. - Below, it seems that we do indeed check only half the number of bridges of
before, but on the other hand, we do more work in remove_source, and this
appears not to be worth it.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
| 22:44:42 (00:00:00) INFO Reading unitigs..
22:45:36 (00:00:53) INFO Read hprcv2-k63-unitigs.fa.zst: 4914797023 bytes, 186730681 records, 17488538356 bases, 93.7 bases/record
22:45:38 (00:00:55) INFO Building graph of 187M unitigs..
22:45:48 (00:01:06) INFO Sorting heads
22:45:55 (00:01:12) INFO Sorting tails
22:46:01 (00:01:18) INFO Matching unitig ends
22:46:06 (00:01:23) INFO Sorting 534M edges
22:46:09 (00:01:26) INFO Building CSR adjacency
22:46:09 (00:01:26) INFO Filling offsets: 1.49GB
22:46:10 (00:01:27) INFO Mapping pairs
22:46:11 (00:01:28) INFO Matching 187M unitigs, 373M ends, k 63: two ends per unitig; overlap = k minus one minus distance; examined = bridge entries checked; links made = connections; a self-link occupies one end, other links occupy two
22:46:11 (00:01:28) INFO Meeting search: source and receiver fronts each store distances below 31; a path through one unitig has cost left distance + unitig extensions + right distance
22:46:11 (00:01:28) INFO Before distance 0: 0 receiving ends examined, 0 links made, 373M receiving ends remain, 373M free ends, 17.5G estimated bases
22:46:11 (00:01:28) INFO Node visits count repeated inspections: direct_overlap scans graph neighbors at distance zero; build_distances initializes bounded shortest paths; meet checks and schedules central unitig traversals and traces endpoints; connect touches link slots; remove_source repairs both bounded fronts
22:46:57 (00:02:14) INFO After distance 0: 246M receiving ends examined, 127M links made, 119M receiving ends remain, 119M free ends, 9.59G estimated bases
22:46:57 (00:02:14) INFO After distance 0 node visits: direct_overlap 258M (total 258M), build_distances 1.35G, meet 0 (total 0), connect 255M (total 255M), remove_source 0 (total 0)
22:47:20 (00:02:38) INFO After distance 1: 8.58M bridge entries examined, 2.86M links made, 113M receiving ends remain, 113M free ends, 9.41G estimated bases
22:47:20 (00:02:38) INFO After distance 1 node visits: direct_overlap 0 (total 258M), meet 29.9M (total 707M), connect 5.71M (total 261M), remove_source 295M (total 295M)
22:47:31 (00:02:48) INFO After distance 2: 7.79M bridge entries examined, 2.17M links made, 109M receiving ends remain, 109M free ends, 9.28G estimated bases
22:47:31 (00:02:48) INFO After distance 2 node visits: direct_overlap 0 (total 258M), meet 25.9M (total 733M), connect 4.34M (total 265M), remove_source 273M (total 568M)
22:50:12 (00:05:30) INFO After distance 32: 2.17M bridge entries examined, 320k links made, 70.8M receiving ends remain, 70.8M free ends, 8.37G estimated bases
22:50:12 (00:05:30) INFO After distance 32 node visits: direct_overlap 0 (total 258M), meet 12.8M (total 1.20G), connect 640k (total 303M), remove_source 90.2M (total 4.75G)
^C
|
17.1 Continuing matchtigs again
Link to heading
17.1.1 Meet-in-the-middle for real
Link to heading
Let’s try mitm again:
- Dijkstra from both sides
- Relax all edges at distance d
- We iterate over all nodes; not sure if we should make a queue
- Push new bridges to a queue for their cost.
- find & use all bridges with total length 2d-1 and 2d or so.
- We iterate over all edges; not sure if we should make a queue
- Note that a bridge is only considered if it is the ‘middle’ edge: the cost at
both ends must be at most half the total length of the path.
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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
| 10:25:49 (00:00:00) INFO Reading unitigs..
10:26:41 (00:00:52) INFO Read hprcv2-k63-unitigs.fa.zst: 4914797023 bytes, 186730681 records, 17488538356 bases, 93.7 bases/record
10:26:43 (00:00:54) INFO Building graph of 187M unitigs..
10:26:53 (00:01:04) INFO Sorting heads
10:26:59 (00:01:10) INFO Sorting tails
10:27:04 (00:01:15) INFO Matching unitig ends
10:27:11 (00:01:21) INFO Sorting 534M edges
10:27:13 (00:01:24) INFO Building CSR adjacency
10:27:13 (00:01:24) INFO Filling offsets: 1.49GB
10:27:15 (00:01:26) INFO Mapping pairs
10:27:16 (00:01:27) INFO Matching 187M unitigs, 373M ends, k 63: two ends per unitig; overlap = k minus one minus distance; examined = bridge entries checked; links made = connections; a self-link occupies one end, other links occupy two
10:27:16 (00:01:27) INFO Meeting search: grow source and receiver distances below 31 one layer at a time; a path through one unitig has cost left distance + unitig extensions + right distance
10:27:16 (00:01:27) INFO Before distance 0: 0 receiving ends examined, 0 links made, 373M receiving ends remain, 373M free ends, 17.5G estimated bases
10:27:16 (00:01:27) INFO Node visits count repeated inspections: direct_overlap scans graph neighbors at distance zero; build_distances scans and relaxes distance layers; meet scans unitig traversals and traces endpoints; connect touches link slots; remove_source repairs settled distances
10:27:33 (00:01:44) INFO Half-distance 0: 119M nodes settled
10:27:33 (00:01:44) INFO After distance 0: 246M receiving ends examined, 127M links made, 119M receiving ends remain, 119M free ends, 9.59G estimated bases
10:27:33 (00:01:44) INFO After distance 0 node visits: direct_overlap 258M (total 258M), build_distances 987M, meet 0 (total 0), connect 255M (total 255M), remove_source 0 (total 0)
10:27:36 (00:01:47) INFO Half-distance 1: 18.9M nodes settled
10:27:40 (00:01:51) INFO After distance 1: 2.86M bridges found, 2.91M entries examined, 2.86M links made, 113M receiving ends remain, 113M free ends, 9.41G estimated bases
10:27:40 (00:01:51) INFO After distance 1 node visits: direct_overlap 0 (total 258M), build_distances 485M (total 1.47G), meet 279M (total 279M), connect 5.71M (total 261M), remove_source 62.9M (total 62.9M)
10:27:45 (00:01:56) INFO After distance 2: 2.81M bridges found, 2.86M entries examined, 2.17M links made, 109M receiving ends remain, 109M free ends, 9.28G estimated bases
10:27:45 (00:01:56) INFO After distance 2 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 1.47G), meet 302M (total 581M), connect 4.34M (total 265M), remove_source 27.7M (total 90.5M)
10:27:47 (00:01:58) INFO Half-distance 2: 8.71M nodes settled
10:27:51 (00:02:02) INFO After distance 3: 2.68M bridges found, 2.71M entries examined, 1.73M links made, 105M receiving ends remain, 105M free ends, 9.18G estimated bases
10:27:51 (00:02:02) INFO After distance 3 node visits: direct_overlap 0 (total 258M), build_distances 444M (total 1.92G), meet 312M (total 893M), connect 3.45M (total 268M), remove_source 28.1M (total 119M)
10:27:56 (00:02:07) INFO After distance 4: 1.77M bridges found, 1.79M entries examined, 1.40M links made, 102M receiving ends remain, 102M free ends, 9.10G estimated bases
10:27:56 (00:02:07) INFO After distance 4 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 1.92G), meet 316M (total 1.21G), connect 2.80M (total 271M), remove_source 21.1M (total 140M)
10:27:57 (00:02:08) INFO Half-distance 3: 5.38M nodes settled
10:28:02 (00:02:13) INFO After distance 5: 1.72M bridges found, 1.74M entries examined, 1.18M links made, 99.9M receiving ends remain, 99.9M free ends, 9.03G estimated bases
10:28:02 (00:02:13) INFO After distance 5 node visits: direct_overlap 0 (total 258M), build_distances 423M (total 2.34G), meet 324M (total 1.53G), connect 2.37M (total 274M), remove_source 20.7M (total 160M)
10:28:07 (00:02:18) INFO After distance 6: 1.30M bridges found, 1.31M entries examined, 1.06M links made, 97.8M receiving ends remain, 97.8M free ends, 8.97G estimated bases
10:28:07 (00:02:18) INFO After distance 6 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 2.34G), meet 328M (total 1.86G), connect 2.12M (total 276M), remove_source 17.5M (total 178M)
10:28:08 (00:02:19) INFO Half-distance 4: 5.04M nodes settled
10:28:14 (00:02:25) INFO After distance 7: 1.38M bridges found, 1.39M entries examined, 975k links made, 95.8M receiving ends remain, 95.8M free ends, 8.92G estimated bases
10:28:14 (00:02:25) INFO After distance 7 node visits: direct_overlap 0 (total 258M), build_distances 417M (total 2.76G), meet 336M (total 2.20G), connect 1.95M (total 278M), remove_source 18.7M (total 197M)
10:28:19 (00:02:30) INFO After distance 8: 1.07M bridges found, 1.08M entries examined, 877k links made, 94.1M receiving ends remain, 94.1M free ends, 8.87G estimated bases
10:28:19 (00:02:30) INFO After distance 8 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 2.76G), meet 339M (total 2.54G), connect 1.75M (total 279M), remove_source 16.2M (total 213M)
10:28:20 (00:02:31) INFO Half-distance 5: 4.22M nodes settled
10:28:25 (00:02:36) INFO After distance 9: 1.13M bridges found, 1.14M entries examined, 809k links made, 92.5M receiving ends remain, 92.5M free ends, 8.83G estimated bases
10:28:25 (00:02:36) INFO After distance 9 node visits: direct_overlap 0 (total 258M), build_distances 411M (total 3.17G), meet 345M (total 2.88G), connect 1.62M (total 281M), remove_source 16.7M (total 229M)
10:28:31 (00:02:42) INFO After distance 10: 873k bridges found, 882k entries examined, 727k links made, 91.0M receiving ends remain, 91.0M free ends, 8.79G estimated bases
10:28:31 (00:02:42) INFO After distance 10 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 3.17G), meet 347M (total 3.23G), connect 1.45M (total 282M), remove_source 14.2M (total 244M)
10:28:32 (00:02:43) INFO Half-distance 6: 3.70M nodes settled
10:28:38 (00:02:49) INFO After distance 11: 944k bridges found, 952k entries examined, 692k links made, 89.6M receiving ends remain, 89.6M free ends, 8.76G estimated bases
10:28:38 (00:02:49) INFO After distance 11 node visits: direct_overlap 0 (total 258M), build_distances 407M (total 3.57G), meet 354M (total 3.58G), connect 1.38M (total 284M), remove_source 15.0M (total 259M)
10:28:43 (00:02:54) INFO After distance 12: 786k bridges found, 795k entries examined, 660k links made, 88.3M receiving ends remain, 88.3M free ends, 8.72G estimated bases
10:28:43 (00:02:54) INFO After distance 12 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 3.57G), meet 357M (total 3.94G), connect 1.32M (total 285M), remove_source 13.9M (total 272M)
10:28:45 (00:02:56) INFO Half-distance 7: 3.29M nodes settled
10:28:50 (00:03:01) INFO After distance 13: 838k bridges found, 846k entries examined, 620k links made, 87.1M receiving ends remain, 87.1M free ends, 8.69G estimated bases
10:28:50 (00:03:01) INFO After distance 13 node visits: direct_overlap 0 (total 258M), build_distances 404M (total 3.98G), meet 362M (total 4.30G), connect 1.24M (total 286M), remove_source 14.5M (total 287M)
10:28:56 (00:03:07) INFO After distance 14: 694k bridges found, 702k entries examined, 586k links made, 85.9M receiving ends remain, 85.9M free ends, 8.66G estimated bases
10:28:56 (00:03:07) INFO After distance 14 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 3.98G), meet 364M (total 4.66G), connect 1.17M (total 288M), remove_source 13.2M (total 300M)
10:28:58 (00:03:09) INFO Half-distance 8: 3.15M nodes settled
10:29:04 (00:03:15) INFO After distance 15: 740k bridges found, 747k entries examined, 557k links made, 84.8M receiving ends remain, 84.8M free ends, 8.64G estimated bases
10:29:04 (00:03:15) INFO After distance 15 node visits: direct_overlap 0 (total 258M), build_distances 402M (total 4.38G), meet 369M (total 5.03G), connect 1.11M (total 289M), remove_source 13.6M (total 314M)
10:29:10 (00:03:21) INFO After distance 16: 628k bridges found, 635k entries examined, 533k links made, 83.7M receiving ends remain, 83.7M free ends, 8.61G estimated bases
10:29:10 (00:03:21) INFO After distance 16 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 4.38G), meet 371M (total 5.41G), connect 1.07M (total 290M), remove_source 12.7M (total 326M)
10:29:11 (00:03:22) INFO Half-distance 9: 2.93M nodes settled
10:29:17 (00:03:28) INFO After distance 17: 691k bridges found, 698k entries examined, 519k links made, 82.7M receiving ends remain, 82.7M free ends, 8.59G estimated bases
10:29:17 (00:03:28) INFO After distance 17 node visits: direct_overlap 0 (total 258M), build_distances 401M (total 4.78G), meet 376M (total 5.78G), connect 1.04M (total 291M), remove_source 13.7M (total 340M)
10:29:24 (00:03:35) INFO After distance 18: 586k bridges found, 593k entries examined, 497k links made, 81.7M receiving ends remain, 81.7M free ends, 8.57G estimated bases
10:29:24 (00:03:35) INFO After distance 18 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 4.78G), meet 378M (total 6.16G), connect 994k (total 292M), remove_source 12.8M (total 353M)
10:29:25 (00:03:36) INFO Half-distance 10: 2.85M nodes settled
10:29:31 (00:03:42) INFO After distance 19: 632k bridges found, 639k entries examined, 476k links made, 80.7M receiving ends remain, 80.7M free ends, 8.55G estimated bases
10:29:31 (00:03:42) INFO After distance 19 node visits: direct_overlap 0 (total 258M), build_distances 400M (total 5.18G), meet 383M (total 6.54G), connect 953k (total 293M), remove_source 13.4M (total 366M)
10:29:38 (00:03:49) INFO After distance 20: 534k bridges found, 541k entries examined, 453k links made, 79.8M receiving ends remain, 79.8M free ends, 8.53G estimated bases
10:29:38 (00:03:49) INFO After distance 20 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 5.18G), meet 385M (total 6.93G), connect 906k (total 294M), remove_source 12.4M (total 379M)
10:29:39 (00:03:50) INFO Half-distance 11: 2.64M nodes settled
10:29:46 (00:03:57) INFO After distance 21: 591k bridges found, 599k entries examined, 446k links made, 78.9M receiving ends remain, 78.9M free ends, 8.51G estimated bases
10:29:46 (00:03:57) INFO After distance 21 node visits: direct_overlap 0 (total 258M), build_distances 398M (total 5.58G), meet 389M (total 7.32G), connect 893k (total 295M), remove_source 13.3M (total 392M)
10:29:53 (00:04:04) INFO After distance 22: 513k bridges found, 521k entries examined, 434k links made, 78.1M receiving ends remain, 78.1M free ends, 8.49G estimated bases
10:29:53 (00:04:04) INFO After distance 22 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 5.58G), meet 391M (total 7.71G), connect 868k (total 295M), remove_source 13.0M (total 405M)
10:29:54 (00:04:05) INFO Half-distance 12: 2.58M nodes settled
10:30:01 (00:04:12) INFO After distance 23: 555k bridges found, 563k entries examined, 420k links made, 77.2M receiving ends remain, 77.2M free ends, 8.48G estimated bases
10:30:01 (00:04:12) INFO After distance 23 node visits: direct_overlap 0 (total 258M), build_distances 398M (total 5.98G), meet 395M (total 8.10G), connect 840k (total 296M), remove_source 13.6M (total 419M)
10:30:07 (00:04:18) INFO After distance 24: 474k bridges found, 481k entries examined, 401k links made, 76.4M receiving ends remain, 76.4M free ends, 8.46G estimated bases
10:30:07 (00:04:18) INFO After distance 24 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 5.98G), meet 396M (total 8.50G), connect 803k (total 297M), remove_source 12.7M (total 431M)
10:30:09 (00:04:20) INFO Half-distance 13: 2.41M nodes settled
10:30:16 (00:04:27) INFO After distance 25: 509k bridges found, 516k entries examined, 388k links made, 75.6M receiving ends remain, 75.6M free ends, 8.45G estimated bases
10:30:16 (00:04:27) INFO After distance 25 node visits: direct_overlap 0 (total 258M), build_distances 397M (total 6.37G), meet 400M (total 8.90G), connect 777k (total 298M), remove_source 13.2M (total 445M)
10:30:23 (00:04:34) INFO After distance 26: 445k bridges found, 451k entries examined, 379k links made, 74.9M receiving ends remain, 74.9M free ends, 8.43G estimated bases
10:30:23 (00:04:34) INFO After distance 26 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 6.37G), meet 401M (total 9.30G), connect 757k (total 299M), remove_source 12.6M (total 457M)
10:30:24 (00:04:35) INFO Half-distance 14: 2.35M nodes settled
10:30:31 (00:04:42) INFO After distance 27: 490k bridges found, 498k entries examined, 374k links made, 74.1M receiving ends remain, 74.1M free ends, 8.42G estimated bases
10:30:31 (00:04:42) INFO After distance 27 node visits: direct_overlap 0 (total 258M), build_distances 396M (total 6.77G), meet 405M (total 9.70G), connect 748k (total 299M), remove_source 13.8M (total 471M)
10:30:38 (00:04:49) INFO After distance 28: 426k bridges found, 433k entries examined, 361k links made, 73.4M receiving ends remain, 73.4M free ends, 8.41G estimated bases
10:30:38 (00:04:49) INFO After distance 28 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 6.77G), meet 407M (total 10.1G), connect 721k (total 300M), remove_source 13.2M (total 484M)
10:30:39 (00:04:50) INFO Half-distance 15: 2.31M nodes settled
10:30:47 (00:04:58) INFO After distance 29: 458k bridges found, 464k entries examined, 347k links made, 72.7M receiving ends remain, 72.7M free ends, 8.40G estimated bases
10:30:47 (00:04:58) INFO After distance 29 node visits: direct_overlap 0 (total 258M), build_distances 396M (total 7.17G), meet 410M (total 10.5G), connect 695k (total 301M), remove_source 13.4M (total 498M)
10:30:54 (00:05:05) INFO After distance 30: 398k bridges found, 405k entries examined, 336k links made, 72.0M receiving ends remain, 72.0M free ends, 8.39G estimated bases
10:30:54 (00:05:05) INFO After distance 30 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 7.17G), meet 411M (total 10.9G), connect 672k (total 301M), remove_source 12.9M (total 511M)
10:30:55 (00:05:06) INFO Half-distance 16: 2.24M nodes settled
10:31:03 (00:05:14) INFO After distance 31: 435k bridges found, 441k entries examined, 328k links made, 71.4M receiving ends remain, 71.4M free ends, 8.38G estimated bases
10:31:03 (00:05:14) INFO After distance 31 node visits: direct_overlap 0 (total 258M), build_distances 395M (total 7.56G), meet 414M (total 11.3G), connect 656k (total 302M), remove_source 13.1M (total 524M)
10:31:10 (00:05:21) INFO After distance 32: 384k bridges found, 391k entries examined, 322k links made, 70.7M receiving ends remain, 70.7M free ends, 8.37G estimated bases
10:31:10 (00:05:21) INFO After distance 32 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 7.56G), meet 416M (total 11.8G), connect 644k (total 303M), remove_source 13.1M (total 537M)
^C
|
This initial implementation has faster initialization but then ends up slower.
Mostly, for every bridge it needs to recompute the distances to both endpoints.
Instead, we will redo the layout so that computed distances do not yet include
the cost of crossing the unitig. Then, computing the total cost of a unitig is trivial.
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
| 10:43:31 (00:00:00) INFO Reading unitigs..
Pending lost events: 232
Pending lost events: 65
10:44:22 (00:00:51) INFO Read hprcv2-k63-unitigs.fa.zst: 4914797023 bytes, 186730681 records, 17488538356 bases, 93.7 bases/record
10:44:23 (00:00:52) INFO Building graph of 187M unitigs..
10:44:35 (00:01:03) INFO Sorting heads
10:44:42 (00:01:10) INFO Sorting tails
10:44:48 (00:01:16) INFO Matching unitig ends
10:44:54 (00:01:22) INFO Sorting 534M edges
10:44:56 (00:01:25) INFO Building CSR adjacency
10:44:56 (00:01:25) INFO Filling offsets: 1.49GB
10:44:58 (00:01:26) INFO Mapping pairs
10:44:59 (00:01:28) INFO Matching 187M unitigs, 373M ends, k 63: two ends per unitig; overlap = k minus one minus distance; examined = bridge entries checked; links made = connections; a self-link occupies one end, other links occupy two
10:44:59 (00:01:28) INFO Meeting search: grow source and receiver distances below 31 one layer at a time; a path through one unitig has cost left distance + unitig extensions + right distance
10:44:59 (00:01:28) INFO Before distance 0: 0 receiving ends examined, 0 links made, 373M receiving ends remain, 373M free ends, 17.5G estimated bases
10:44:59 (00:01:28) INFO Node visits count repeated inspections: direct_overlap scans graph neighbors at distance zero; build_distances scans and relaxes distance layers; meet scans unitig traversals and traces endpoints; connect touches link slots; remove_source repairs settled distances
10:45:13 (00:01:42) INFO Half-distance 0: 117M nodes settled
10:45:13 (00:01:42) INFO After distance 0: 246M receiving ends examined, 127M links made, 119M receiving ends remain, 119M free ends, 9.59G estimated bases
10:45:13 (00:01:42) INFO After distance 0 node visits: direct_overlap 258M (total 258M), build_distances 1.08G, meet 0 (total 0), connect 255M (total 255M), remove_source 0 (total 0)
10:45:16 (00:01:45) INFO Half-distance 1: 18.0M nodes settled
10:45:20 (00:01:49) INFO After distance 1: 2.86M bridges found, 2.91M entries examined, 2.86M links made, 113M receiving ends remain, 113M free ends, 9.41G estimated bases
10:45:20 (00:01:49) INFO After distance 1 node visits: direct_overlap 0 (total 258M), build_distances 460M (total 1.54G), meet 211M (total 211M), connect 5.71M (total 261M), remove_source 109M (total 109M)
10:45:24 (00:01:52) INFO After distance 2: 2.81M bridges found, 2.86M entries examined, 2.17M links made, 109M receiving ends remain, 109M free ends, 9.28G estimated bases
10:45:24 (00:01:52) INFO After distance 2 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 1.54G), meet 208M (total 419M), connect 4.34M (total 265M), remove_source 55.2M (total 164M)
10:45:26 (00:01:55) INFO Half-distance 2: 6.93M nodes settled
10:48:34 (00:05:03) INFO Half-distance 31: 0 nodes settled
10:48:36 (00:05:05) INFO After distance 61: 185k bridges found, 188k entries examined, 172k links made, 57.2M receiving ends remain, 57.2M free ends, 8.26G estimated bases
10:48:36 (00:05:05) INFO After distance 61 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 13.3G), meet 191M (total 11.9G), connect 343k (total 316M), remove_source 18.0M (total 1.50G)
10:48:37 (00:05:06) INFO Final node visits (including closure links): direct_overlap 258M, build_distances 13.3G, meet 11.9G, connect 373M, remove_source 1.50G
10:48:37 (00:05:06) INFO Reconstruct output
10:49:51 (00:06:19) INFO Length summary: 8.26G estimated bases after links, 8.26G output bases
10:49:56 (00:06:24) INFO Wrote hprcv2-k63-unitigs-greedytigs.msfa: 8257928676 bytes, 1 records, 8257928645 bases, 8257928645.0 bases/record
10:49:56 (00:06:24) INFO matchtigs: 384.81s wall [535.37s CPU; 139.1%; MaxRSS 26.18 GiB], output 8.26G bases (hprcv2-k63-unitigs-greedytigs.msfa)
|
Great! It’s faster now 384s (535s CPU)
One weird thing is that in relax_at, which relaxes edges from all nodes at
distance \(d\), it first checks if the distance is still up-to-date. Instead, we
should just make sure all finite distances are always correct directly in
remove_source so the check here can be skipped.
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
| 11:01:55 (00:00:00) INFO Reading unitigs..
Pending lost events: 293
11:02:46 (00:00:50) INFO Read hprcv2-k63-unitigs.fa.zst: 4914797023 bytes, 186730681 records, 17488538356 bases, 93.7 bases/record
11:02:47 (00:00:51) INFO Building graph of 187M unitigs..
11:02:58 (00:01:02) INFO Sorting heads
11:03:04 (00:01:08) INFO Sorting tails
11:03:10 (00:01:14) INFO Matching unitig ends
11:03:16 (00:01:20) INFO Sorting 534M edges
11:03:18 (00:01:22) INFO Building CSR adjacency
11:03:18 (00:01:22) INFO Filling offsets: 1.49GB
11:03:20 (00:01:24) INFO Mapping pairs
11:03:22 (00:01:26) INFO Matching 187M unitigs, 373M ends, k 63: two ends per unitig; overlap = k minus one minus distance; examined = bridge entries checked; links made = connections; a self-link occupies one end, other links occupy two
11:03:22 (00:01:26) INFO Meeting search: grow source and receiver distances below 31 one layer at a time; a path through one unitig has cost left distance + unitig extensions + right distance
11:03:22 (00:01:26) INFO Before distance 0: 0 receiving ends examined, 0 links made, 373M receiving ends remain, 373M free ends, 17.5G estimated bases
11:03:22 (00:01:26) INFO Node visits count repeated inspections: direct_overlap scans graph neighbors at distance zero; build_distances scans and relaxes distance layers; meet scans unitig traversals and traces endpoints; connect touches link slots; remove_source repairs settled distances
11:03:36 (00:01:40) INFO Half-distance 0: 117M nodes settled
11:03:36 (00:01:40) INFO After distance 0: 246M receiving ends examined, 127M links made, 119M receiving ends remain, 119M free ends, 9.59G estimated bases
11:03:36 (00:01:40) INFO After distance 0 node visits: direct_overlap 258M (total 258M), build_distances 1.08G, meet 0 (total 0), connect 255M (total 255M), remove_source 0 (total 0)
11:03:37 (00:01:41) INFO Half-distance 1: 18.0M nodes settled
11:03:43 (00:01:47) INFO After distance 1: 2.86M bridges found, 2.91M entries examined, 2.86M links made, 113M receiving ends remain, 113M free ends, 9.41G estimated bases
11:03:43 (00:01:47) INFO After distance 1 node visits: direct_overlap 0 (total 258M), build_distances 412M (total 1.50G), meet 211M (total 211M), connect 5.71M (total 261M), remove_source 152M (total 152M)
11:03:47 (00:01:51) INFO After distance 2: 2.81M bridges found, 2.86M entries examined, 2.17M links made, 109M receiving ends remain, 109M free ends, 9.28G estimated bases
11:03:47 (00:01:51) INFO After distance 2 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 1.50G), meet 208M (total 419M), connect 4.34M (total 265M), remove_source 103M (total 255M)
11:03:48 (00:01:52) INFO Half-distance 2: 6.93M nodes settled
11:06:34 (00:04:38) INFO Half-distance 31: 0 nodes settled
11:06:36 (00:04:40) INFO After distance 61: 185k bridges found, 188k entries examined, 172k links made, 57.2M receiving ends remain, 57.2M free ends, 8.26G estimated bases
11:06:36 (00:04:40) INFO After distance 61 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 12.5G), meet 191M (total 11.9G), connect 343k (total 316M), remove_source 18.0M (total 2.14G)
11:06:37 (00:04:41) INFO Final node visits (including closure links): direct_overlap 258M, build_distances 12.5G, meet 11.9G, connect 373M, remove_source 2.14G
11:06:37 (00:04:41) INFO Reconstruct output
11:07:49 (00:05:53) INFO Length summary: 8.26G estimated bases after links, 8.26G output bases
11:07:53 (00:05:57) INFO Wrote hprcv2-k63-unitigs-greedytigs.msfa: 8257928676 bytes, 1 records, 8257928645 bases, 8257928645.0 bases/record
11:07:53 (00:05:57) INFO matchtigs: 357.53s wall [503.61s CPU; 140.9%; MaxRSS 26.18 GiB], output 8.26G bases (hprcv2-k63-unitigs-greedytigs.msfa)
|
- 357s! (503s CPU) Another half-minute speedup, where the first 90s and last
60s are IO.
We can get some more perf by scanning the list of all unitigs for bridges not
for every distance, but only every other distance. This again saves 35s, for
321s total!
Currently we scan all bridges at every distances. Better is to keep a bucket
queue storing ’live’ bridges: 306s
- We can first determine all cycles in the graph, and then in parallel write
those out.
- We have to use
with_max_len(1) to ensure a single thread does not grab all
adjacent heavy paths. We further avoid this by splitting cycles every 32MB. - To shrink output a bit, short cycles are split at a min-overlap point.
- Now it’s 20s to find the cycles and 10s to concatenate them over multiple threads.
Result: 265s
- At the end, the algorithm was joining unmatched unitig ends into cycles, but that’s not
allowed, as it creates new kmers. Instead these should remain as loose ends.
Time to use heaptrack. We start out at 26.2GB:
- Most of it is input + ranges + temporary 128bit head and tail of every node.
- Reduce by only storing the smaller of fwd and rc link?
- Instead, it seems the peak is reached when constructing the final output, as
this needs many things in memory:
- the bitpacking input sequence
- the ranges of it
- the links, for the overlaps
- the paths of unitig IDs
- the work-in-progress output
It does help slightly to eagerly drop data structures as soon as we don’t need
them anymore: 17.86GB MaxRSS during the algorithm, and ~21GB at the end.
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
| 15:51:27 (00:00:00) INFO Reading unitigs..
15:51:50 (00:00:22) INFO Seq: 4.37G GB
15:51:50 (00:00:22) INFO Ranges: 2.99G GB
15:51:50 (00:00:22) INFO Read hprcv2-k63-unitigs.fa: 21181450134 bytes, 186730681 records, 17488538356 bases, 93.7 bases/record
15:51:52 (00:00:24) INFO Building graph of 187M unitigs..
15:51:55 (00:00:27) INFO Sorting heads (3.81G B)
15:51:58 (00:00:30) INFO Sorting tails (3.66G B)
15:52:00 (00:00:33) INFO Matching unitig ends
15:52:05 (00:00:37) INFO Sorting 534M edges
15:52:08 (00:00:40) INFO Building CSR adjacency
15:52:08 (00:00:40) INFO Filling offsets: 1.49GB
15:52:09 (00:00:42) INFO Mapping pairs
15:52:10 (00:00:43) INFO Matching 187M unitigs, 373M ends, k 63: two ends per unitig; overlap = k minus one minus distance; examined = bridge entries checked; links made = connections; a self-link occupies one end, other links occupy two
15:52:10 (00:00:43) INFO Meeting search: grow source and receiver distances below 31 one layer at a time; a path through one unitig has cost left distance + unitig extensions + right distance
15:52:10 (00:00:43) INFO Before distance 0: 0 receiving ends examined, 0 links made, 373M receiving ends remain, 373M free ends, 17.5G estimated bases
15:52:10 (00:00:43) INFO Node visits count repeated inspections: direct_overlap scans graph neighbors at distance zero; build_distances scans and relaxes distance layers; meet checks and queues bridges and traces endpoints; connect touches link slots; remove_source repairs distance labels
15:52:29 (00:01:02) INFO Half-distance 0: 117M nodes settled
15:52:29 (00:01:02) INFO After distance 0: 246M receiving ends examined, 127M links made, 119M receiving ends remain, 119M free ends, 9.59G estimated bases
15:52:29 (00:01:02) INFO After distance 0 node visits: direct_overlap 258M (total 258M), build_distances 1.08G, meet 117M (total 117M), connect 255M (total 255M), remove_source 0 (total 0)
15:52:31 (00:01:04) INFO Half-distance 1: 18.0M nodes settled
15:52:37 (00:01:09) INFO After distance 1: 5.73M bridge entries queued, 5.78M entries examined, 2.86M links made, 113M receiving ends remain, 113M free ends, 9.41G estimated bases
15:52:37 (00:01:09) INFO After distance 1 node visits: direct_overlap 0 (total 258M), build_distances 412M (total 1.50G), meet 71.0M (total 188M), connect 5.71M (total 261M), remove_source 152M (total 152M)
15:52:41 (00:01:13) INFO After distance 2: 6.37M bridge entries queued, 6.42M entries examined, 2.17M links made, 109M receiving ends remain, 109M free ends, 9.28G estimated bases
15:52:41 (00:01:13) INFO After distance 2 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 1.50G), meet 41.0M (total 229M), connect 4.34M (total 265M), remove_source 103M (total 255M)
15:52:42 (00:01:15) INFO Half-distance 2: 6.93M nodes settled
15:54:46 (00:03:19) INFO After distance 60: 850k bridge entries queued, 853k entries examined, 175k links made, 57.6M receiving ends remain, 57.6M free ends, 8.26G estimated bases
15:54:46 (00:03:19) INFO After distance 60 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 12.5G), meet 7.82M (total 1.09G), connect 351k (total 316M), remove_source 19.4M (total 2.12G)
15:54:46 (00:03:19) INFO Half-distance 31: 0 nodes settled
15:54:47 (00:03:20) INFO After distance 61: 811k bridge entries queued, 814k entries examined, 172k links made, 57.2M receiving ends remain, 57.2M free ends, 8.26G estimated bases
15:54:47 (00:03:20) INFO After distance 61 node visits: direct_overlap 0 (total 258M), build_distances 0 (total 12.5G), meet 7.48M (total 1.09G), connect 344k (total 316M), remove_source 18.2M (total 2.14G)
15:54:47 (00:03:20) INFO Final node visits: direct_overlap 258M, build_distances 12.5G, meet 1.09G, connect 316M, remove_source 2.14G
15:54:48 (00:03:20) INFO Dropping graph
15:54:48 (00:03:20) INFO Reconstruct output
15:55:22 (00:03:55) INFO Reconstructing 28.7M parts in parallel
15:55:37 (00:04:09) INFO Dropping links
15:55:37 (00:04:10) INFO Length summary: 8.26G estimated bases after links, 8.26G output bases in 28.7M records
15:55:49 (00:04:21) INFO Wrote hprcv2-k63-greedytigs.fa: 8819238227 bytes, 28654356 records, 8257262210 bases, 288.2 bases/record
15:55:49 (00:04:21) INFO matchtigs: 261.64s wall [487.26s CPU; 186.2%; MaxRSS 21.37 GiB], output 8.26G bases in 28.7M records (hprcv2-k63-greedytigs.fa)
|
- One additional thing we can do is including overlap lengths in the paths we
construct, so that we can
drop(links) before constructing the string output: 19.5GB - While collecting edges of the graph, we find matching elements in
heads and
tails via a merge sort. Instead, we can iterate backwards and free memory in
1MB blocks as we process them. 14.7GB during the main loop (but still 19.5
at the end) - We can replace
Vec<Range<usize>> by simply Vec<usize> since ranges have
matching ends and starts. 13.4GB during the main loop (18.1GB MaxRSS now
at the end!) - Lastly, we could first build the bitpacked representation of the output, then
drop everything else, and then directly stream the unpacked representation to
disk, but this won’t give too much since we’re now already close to the 17.5GB that is
used during the bulk of the algorithm.
- Instead, we can also flush everything to disk directly: each thread can build
a small buffer and then take the global file writer lock. => 13.4GB overall
MaxRSS :)