vg autoindex on Minigraph output fails

[gamerino]

1. What were you trying to do?
Build vg giraffe (sr) indexes for a pangenome graph built with minigraph, combining GRCh38, CHM13, and GRCh37 using the first one as.

2. What did you want to happen?
Generate the indexes that vg giraffe requires for mapping short reads against my graph genome for then being able to surject the alignment to GRCh38 (and the same for each assembly)

3. What actually happened?

vg autoindex crashed while building the GBWT

[vg autoindex] Executing command: vg autoindex -w sr-giraffe --prefix vg_index/GRCh38-CHM13-GRCh37_v2_index -r fasta/GRCh38/GRCh38.fa --gfa mg_out/GRCh38-CHM13-GRCh37_v2.gfa -t 6 --target-mem 180G
[IndexRegistry] Checking for haplotype lines in GFA.
[IndexRegistry] Constructing VG graph from GFA input.
[IndexRegistry] Constructing XG graph from VG graph.
[IndexRegistry] Constructing a greedy path cover GBWT
terminate called after throwing an instance of 'std::runtime_error'
  what():  MetadataBuilder: Duplicate path for sample path_cover_0 contig 16 phase 0 count 0
━━━━━━━━━━━━━━━━━━━━
Crash report for vg v1.70.0 "Zebedassi"
Caught signal 6 raised at address 0x23d333c; tracing with backward-cpp
Stack trace (most recent call last):
#19   Object "/opt/conda/bin/vg", at 0x6824b4, in _start
#18   Object "/opt/conda/bin/vg", at 0x238f606, in __libc_start_main
#17   Object "/opt/conda/bin/vg", at 0x238dd69, in __libc_start_call_main
#16   Object "/opt/conda/bin/vg", at 0x6830cb, in vg_main(int, char**)
#15   Object "/opt/conda/bin/vg", at 0xfb4ddb, in vg::subcommand::Subcommand::operator()(int, char**) const
#14   Object "/opt/conda/bin/vg", at 0xded126, in main_autoindex(int, char**)
#13   Object "/opt/conda/bin/vg", at 0x12e20f3, in vg::IndexRegistry::make_indexes(std::vector<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&)
#12   Object "/opt/conda/bin/vg", at 0x12d21bc, in vg::IndexRegistry::execute_recipe(std::pair<std::set<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, std::allocator<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > > >, unsigned long> const&, vg::IndexingPlan const*, vg::AliasGraph&)
#11   Object "/opt/conda/bin/vg", at 0x12be9bd, in std::_Function_handler<std::vector<std::vector<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > > >, std::allocator<std::vector<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > > > > > (std::vector<vg::IndexFile const*, std::allocator<vg::IndexFile const*> > const&, vg::IndexingPlan const*, vg::AliasGraph&, std::set<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, std::allocator<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&), vg::VGIndexes::get_vg_index_registry()::{lambda(std::vector<vg::IndexFile const*, std::allocator<vg::IndexFile const*> > const&, vg::IndexingPlan const*, vg::AliasGraph&, std::set<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, std::allocator<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&)#35}>::_M_invoke(std::_Any_data const&, std::vector<vg::IndexFile const*, std::allocator<vg::IndexFile const*> > const&, vg::IndexingPlan const*&&, vg::AliasGraph&, std::set<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, std::allocator<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&)
#10   Object "/opt/conda/bin/vg", at 0x12be71f, in vg::VGIndexes::get_vg_index_registry()::{lambda(std::vector<vg::IndexFile const*, std::allocator<vg::IndexFile const*> > const&, vg::IndexingPlan const*, vg::AliasGraph&, std::set<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, std::allocator<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&)#35}::operator()(std::vector<vg::IndexFile const*, std::allocator<vg::IndexFile const*> > const&, vg::IndexingPlan const*, vg::AliasGraph&, std::set<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, std::allocator<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&) const [clone .constprop.0]
#9    Object "/opt/conda/bin/vg", at 0x1971225, in gbwtgraph::path_cover_gbwt(handlegraph::PathHandleGraph const&, gbwtgraph::PathCoverParameters const&, bool, std::function<bool (handlegraph::path_handle_t const&)> const*)
#8    Object "/opt/conda/bin/vg", at 0x1989e25, in gbwtgraph::MetadataBuilder::add_haplotype(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, unsigned long, unsigned long, unsigned long)
#7    Object "/opt/conda/bin/vg", at 0x5f4abc, in gbwtgraph::MetadataBuilder::add_path(handlegraph::PathSense, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, unsigned long, unsigned long, std::pair<unsigned long, unsigned long> const&, unsigned long) [clone .cold]
#6    Object "/opt/conda/bin/vg", at 0x22c7b08, in __cxa_throw
#5    Object "/opt/conda/bin/vg", at 0x22c79a6, in std::terminate()
#4    Object "/opt/conda/bin/vg", at 0x22c793b, in __cxxabiv1::__terminate(void (*)())
#3    Object "/opt/conda/bin/vg", at 0x63cfe7, in __gnu_cxx::__verbose_terminate_handler() [clone .cold]
#2    Object "/opt/conda/bin/vg", at 0x63f72f, in abort
#1    Object "/opt/conda/bin/vg", at 0x23a67e5, in raise
#0    Object "/opt/conda/bin/vg", at 0x23d333c, in __pthread_kill

Library locations:
ERROR: Signal 6 occurred. VG has crashed. Visit https://github.com/vgteam/vg/issues/new/choose to report a bug.
━━━━━━━━━━━━━━━━━━━━
Context dump:
        Thread 0: Starting 'autoindex' subcommand
Found 1 threads with context.
━━━━━━━━━━━━━━━━━━━━
Please include this entire error log in your bug report!

4. If you got a line like Stack trace path: /somewhere/on/your/computer/stacktrace.txt, please copy-paste the contents of that file here:

Don't have it

5. What data and command can the vg dev team use to make the problem happen?

OUTDIR=vg_index

REFERENCE_FA=fasta/GRCh38/GRCh38.fa

PANGENOME_RGFA=mg_out/GRCh38-CHM13-GRCh37_v2.gfa

PANGENOME_IDX=GRCh38-CHM13-GRCh37_v2_index

vg autoindex -w sr-giraffe --prefix ${OUTDIR}/${PANGENOME_IDX} -r ${REFERENCE_FA} --gfa ${PANGENOME_RGFA} -t 6 --target-mem 180G"

GFA was built using the command line below. GRCh38.fa has chromosome names as 1, 2, ..., X, Y, while the other two assembly files have CHM13#1, CHM13#2, ... and GRCh37#1, GRCh37#2, ... respectively

"minigraph -c -x ggs --ggen -t 10 fasta/GRCh38/GRCh38.fa T2T-CHM13v2.0.fa GRCh37.fa > mg_out/GRCh38-CHM13-GRCh37_v2.gfa"

6. What does running vg version say?
1.7.0

[jltsiren]

Can you share the graph?

Also, unrelated to the issue at hand, it looks like vg autoindex does not handle minigraph graphs properly. When you have a GFA graph built by aligning sequences, it should convert the GFA directly to GBZ without passing through .vg and .xg. But somehow vg autoindex does not recognize the paths as haplotypes and assumes that the graph contains nodes not covered by any path.

You should be able to avoid the issue by building the GBZ directly and then giving it to vg autoindex:

vg gbwt -G graph.gfa --gbz-format -g graph.gbz
vg autoindex -w sr-giraffe --gbz graph.gbz

But even if that works, I would still like to find the cause for the original failure.

[gamerino]

Thanks!

I can share the file, which is the best way to do it?

I would also like to understand why it happened because when I run the same pipeline but using CHM13 as backbone for minigraph and then when building the giraffe index it worked without any error. I was also able to run vg giraffe using the generated index and then surject my alignments to the linear version of CHM12. However, when I tried with GRCh38 or with GRCh37 as backbone, I got the same error I shared in my first comment.

I tried running what you suggested to build the GBZ first but it failed

terminate called after throwing an instance of 'std::runtime_error'
  what():  No paths or walks in the GFA file
━━━━━━━━━━━━━━━━━━━━
Crash report for vg v1.70.0 "Zebedassi"
Caught signal 6 raised at address 0x23d333c; tracing with backward-cpp
Stack trace (most recent call last):
#15   Object "/opt/conda/bin/vg", at 0x6824b4, in _start
#14   Object "/opt/conda/bin/vg", at 0x238f606, in __libc_start_main
#13   Object "/opt/conda/bin/vg", at 0x238dd69, in __libc_start_call_main
#12   Object "/opt/conda/bin/vg", at 0x6830cb, in vg_main(int, char**)
#11   Object "/opt/conda/bin/vg", at 0xfb4ddb, in vg::subcommand::Subcommand::operator()(int, char**) const
#10   Object "/opt/conda/bin/vg", at 0xeae2dd, in main_gbwt(int, char**)
#9    Object "/opt/conda/bin/vg", at 0xeabbb6, in step_1_build_gbwts(vg::GBWTHandler&, GraphHandler&, GBWTConfig&)
#8    Object "/opt/conda/bin/vg", at 0x19151b1, in gbwtgraph::gfa_to_gbwt(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, gbwtgraph::GFAParsingParameters const&)
#7    Object "/opt/conda/bin/vg", at 0x5f0ab3, in gbwtgraph::check_gfa_file(gbwtgraph::GFAFile const&, gbwtgraph::GFAParsingParameters const&) [clone .cold]
#6    Object "/opt/conda/bin/vg", at 0x22c7b08, in __cxa_throw
#5    Object "/opt/conda/bin/vg", at 0x22c79a6, in std::terminate()
#4    Object "/opt/conda/bin/vg", at 0x22c793b, in __cxxabiv1::__terminate(void (*)())
#3    Object "/opt/conda/bin/vg", at 0x63cfe7, in __gnu_cxx::__verbose_terminate_handler() [clone .cold]
#2    Object "/opt/conda/bin/vg", at 0x63f72f, in abort
#1    Object "/opt/conda/bin/vg", at 0x23a67e5, in raise
#0    Object "/opt/conda/bin/vg", at 0x23d333c, in __pthread_kill

Library locations:
ERROR: Signal 6 occurred. VG has crashed. Visit https://github.com/vgteam/vg/issues/new/choose to report a bug.
━━━━━━━━━━━━━━━━━━━━
Context dump:
        Thread 0: Starting 'gbwt' subcommand
Found 1 threads with context.
━━━━━━━━━━━━━━━━━━━━
Please include this entire error log in your bug report!

Thanks!

[jltsiren]

Ok, I forgot that minigraph does not include the underlying sequences as paths vg can use. Then there is no choice but to generate arbitrary paths that cover the graph in order to use it in Giraffe. And because none of the underlying samples (CHM13, GRCh37, GRCh38) is included in the graph, you cannot surject the alignments to any reference.

You can share the graph in any way that works for you.