vg index forever at 100% cpu

[NickSto]

I'm trying to index a small graph, but the vg index command just keeps working at 100% cpu forever. Or, at least I assume it's forever, since construct generated the graph in a few seconds, but index has spent over an hour on it.

The command I'm using is:
$ vg index -k 11 -x nvc-filt.xg -g nvc-filt.gcsa nvc-filt.vg
This is with the v1.4.0-1347-g846bf0b static binary, but it's also happened with a non-static version run from docker. The .xg file is created instantly and then never changes size, and the .gcsa file never appears. I've also tried making a rocksdb index, but that also never completes.

I tried editing the original vcf to include only a few of the indels I'm interested in, and it successfully indexes that graph (nvc-curated.vg) in like a second:

$ grep -v '^##' nvc-curated.vcf | cut -f -7
#CHROM  POS ID  REF ALT QUAL    FILTER
SC8-ch  8268    .   TAGC    T   .   .
SC8-ch  8271    .   C   CACCCCCTCT,CT,CCCCTCT   .   .
SC8-ch  8289    .   TAGA    TA,T,TGA    .   .
SC8-ch  8290    .   AGAG    A,AAG   .   .

Maybe my graph is unusually complex. It's short, with the non-linear sequence only occurring over 1kb. But there are over 1,000 nodes in that space. And some areas are really complicated:

Supporting files:
asm.fa: reference sequence
nvc-filt.vg: problem graph
nvc-filt.vcf: VCF used to create it
nvc-curated.vg: subset which indexes fine
nvc-curated.vcf: VCF used to create it
nvc-filt.xg: incomplete xg file created early in the run of index

[ekg]

@NickSto I think that graph complexity is probably getting the indexer stuck.

To see where it's doing it, try running with progress reporting. -p I think. This will show at what stage the slowdown happens.

My guess is that there are way too many kmers in some of these regions.

A quick trick is to cut the graph at places where there are a series of bifurcations prior to generating the kmer index. The original graph is unchanged, as manipulation and local alignment don't necessarily suffer from the complexity as kmer enumeration. It's just that the kmers used to build up the GCSA index will need to be subset. The cutting is built into vg index:

vg index -k 11 --edge-max 3 -x nvc-filt.xg -g nvc-filt.gcsa nvc-filt.vg

This will not produce kmers that cross more than 3 bifurcations.

It is also possible to drop the reference path kmers back in. This means that we always have at least the reference sequence to go from. @adamnovak do we have an example of this from the wiki?

There is a technique which @jltsiren has implemented in GCSA2 that will allow us to reduce the complexity of the indexing in these regions provided we retain a set of reference paths (or haplotypes), but I have not matched up vg to that interface nor tested it. Doing so correctly is non trivial.

[jltsiren]

Extracting the 11-mers from that graph with vg kmers takes 35-40 CPU hours. The resulting .graph file contains only 2.8 million kmers, so there are probably many duplicates that get filtered out.

It's possible to build a 22-mer index from the kmers, but the construction requires over 60 gigabytes of memory. The resulting index takes around 33 MB, because most of the 22-mers are redundant. The extend() step produces 1.8 billion 22-mers, but they are merged into 6.5 million nodes in the pruned de Bruijn graph. The intermediate files produced by the second doubling step exceed the 500 GB size limit, so building a 44-mer index would be impossible with the hardware I can use.

Using the terminology from the original GCSA paper, the graph is highly nondeterministic. For example, node 721 has two predecessors (712 and 720) with label T. For every path starting with (712, 721), there is a path starting with (720, 721) with the same label. Hence index construction must generate all such paths up to the desired length (e.g. 44 or 88). This causes the exponential growth in the number of intermediate paths.

It may be possible to avoid the exponential growth by determinizing the graph before indexing. This would require three new components:

1. Transform the variation graph explicitly into a GCSA2 input graph with directed edges and single-character node labels.
2. Determinize the input graph and provide a mapping from the nodes of the determinized graph to (sets of) positions in the variation graph.
3. Extract the input kmers from the determinized input graph.

The mapping produced by component 2 would be used in the final phase GCSA2 construction to create an index that maps to the original variation graph.

[ekg]

I am not sure that determinizing will help. Maybe @jltsiren could advise.

My impression is that we must limit the input redundancy in some way.

On Sat, Oct 15, 2016, 19:31 Jouni Siren notifications@github.com wrote:

Extracting the 11-mers from that graph with vg kmers takes 35-40 CPU
hours. The resulting .graph file contains only 2.8 million kmers, so
there are probably many duplicates that get filtered out.

It's possible to build a 22-mer index from the kmers, but the construction
requires over 60 gigabytes of memory. The resulting index takes around 33
MB, because most of the 22-mers are redundant. The extend() step produces
1.8 billion 22-mers, but they are merged into 6.5 million nodes in the
pruned de Bruijn graph. The intermediate files produced by the second
doubling step exceed the 500 GB size limit, so building a 44-mer index
would be impossible with the hardware I can use.

Using the terminology from the original GCSA paper, the graph is highly
nondeterministic. For example, node 721 has two predecessors (712 and 720)
with label T. For every path starting with (712, 721), there is a path
starting with (720, 721) with the same label. Hence index construction must
generate all such paths up to the desired length (e.g. 44 or 88). This
causes the exponential growth in the number of intermediate paths.

It may be possible to avoid the exponential growth by determinizing the
graph before indexing. This would require three new components:

1. Transform the variation graph explicitly into a GCSA2 input graph
   with directed edges and single-character node labels.
2. Determinize the input graph and provide a mapping from the nodes of
   the determinized graph to (sets of) positions in the variation graph.
3. Extract the input kmers from the determinized input graph.

The mapping produced by component 2 would be used in the final phase GCSA2
construction to create an index that maps to the original variation graph.

—
You are receiving this because you commented.

Reply to this email directly, view it on GitHub
#511 (comment), or mute
the thread
https://github.com/notifications/unsubscribe-auth/AAI4EbJQtQkwmowKY4xB_-Jiyn-2nAdIks5q0Q3-gaJpZM4KXmo-
.

[jltsiren]

Here we have 2.8 billion 22-mers after the first extend() step. If we do a prune() step, merging a lexicographic range of 22-mers if they all start from the same position, we end up with 576 million paths. If we do a merge() step instead, merging the range if every 22-mer starts from the same set of positions, as well as all identical 22-mers regardless of starting positions, we get just 6.5 million path nodes.

The 90-fold reduction in the number of paths comes from two sources:

1. Merging ranges of 22-mers, where each distinct 22-mer in the range starts from the same set (of size >1) of positions. This can be significant, if we have nondeterministic edges in complex regions.
2. Merging identical 22-mers from different parts of the graph. This can be significant, if we have repetitive regions with enough variants.

If reason 1 is dominant, determinization can help greatly, as it should reduce the number 22-mers created by extend() by a similar factor. Simplifying the graph helps with both 1 and 2, but we lose potentially interesting paths in the process. (Simplifying also makes the final index smaller.)

[NickSto]

Wow, alright, thank you for confirming my suspicions.

So I tried --edge-max 3, but that didn't do the trick. I turned on progress reporting, though, and the transcripts are here and here for the two regions described below.

It sounds like @jltsiren has already gone way past this, but just FYI, I've narrowed down the problem region(s) and observed more of the behavior of index.

It turns out there's at least two parts of the graph that create this problem, and it's not in the area I pictured. It was able to index that area in only 15 seconds.

I've narrowed it to two regions, each containing about 400 nodes and covering 250bp of the reference: coords 7752-8000 and 8001-8256. And I've now let each run until they fail, and they're different types of failures:

7752-8000 (4 hr runtime):

terminate called after throwing an instance of 'std::bad_alloc'
  what():  std::bad_alloc
Aborted (core dumped)

8001-8256 (1 hr runtime):

PathGraphBuilder::write(): Size limit exceeded, construction aborted

They ran until their temp files filled my entire hard drive (about 9GB of space). The bulk of the space was taken by the .gcsa_* files.

The regions:
nvc-debug.7752-8000.vg: graph for 7752-8000
nvc-debug.7752-8000.vcf: VCF used to create it
nvc-debug.8001-8256.vg: graph for 8001-8256
nvc-debug.8001-8256.vcf: VCF used to create it
asm.fa: the reference sequence

[jltsiren]

I have an idea that can be implemented fairly quickly:

1. We can already extract the input graph with vg kmers -gB -k 1.
2. The determinization code from the old GCSA should work with minimal changes, as long the sink node is reachable from every node in the variation graph.
3. Because we extract the input graph with vg kmers, we can easily produce 16-mers for index construction with four steps of prefix-doubling.

[NickSto]

I haven't had much time for this issue recently, but I just wanted to add a quick note that I've narrowed the problem areas further. I don't have the graph files, but here are the outer bounds of the regions:

7795-7800
7815-7825
7840-7855
7990-8070
7940-7960

I've managed to easily index a graph that excludes VCF variants from those regions, so whatever the problems are, they're confined to those areas.