vgteam
ERROR: Signal 6 occurred. VG has crashed - vg safari crashes with short paired-end reads
1. What were you trying to do?
I am trying to align short aDNA reads against a Minigraph-Cactus pangenome made of multiple chr-level haplotype-resolved genomes using vg safari. I didn't merge overlapping reads. The median read length is 90 bp.
I have two batches of reads:
- with fragment length distribution and stdv (estimated by vg safari, e.g. 161.12 +/- 25.6574 or 206.627 +/- 30.5905) that have a final Fragment distance limit > read distance limit and vg safari doesn't fall back into single-end mapping
- with fragment length distribution and stdv (estimated by vg safari, e.g 92.3316 +/- 18.5327) that have a final Fragment distance limit < read distance limit and makes vg safari falling back into single-end mapping.
Any help in setting the right parameters to be able to align reads with these fragment length distributions are appreciated.
2. What did you want to happen?
I wanted to align the reads without falling back into single-end mapping and without crashing.
3. What actually happened?
For the first batch, I didn't specify any fragment length distribution or any other parameter. For most of the fastqs it worked, but for some VG crashed:
Loading Distance Index v2
Initializing MinimizerMapper
Loading and initialization: 43.7847 seconds
Mapping reads to "-" (GAM)
--hit-cap 10
--hard-hit-cap 50000
--score-fraction 0.9
--max-min 500
--num-bp-per-min 1000
--exclude-overlapping-min
--max-extensions 800
--max-alignments 8
--cluster-score 50
--pad-cluster-score 0
--cluster-coverage 0.3
--extension-score 1
--extension-set 20
--max-multimaps 1
--distance-limit 200
--paired-distance-limit 2
--rescue-subgraph-size 4
--rescue-seed-limit 100
--rescue-attempts 15
--rescue-algorithm dozeu
Using fragment length estimate: 206.193 +/- 29.5261
vg: src/dozeu_interface.cpp:126: vg::DozeuInterface::graph_pos_s vg::DozeuInterface::calculate_max_position(const vg::DozeuInterface::OrderedGraph&, const vg::DozeuInterface::graph_pos_s&, size_t, bool, const std::vector<const dz_forefront_s*>&): Assertion `forefronts.at(max_node_index)->mcap != nullptr' failed.
ERROR: Signal 6 occurred. VG has crashed. Visit https://github.com/vgteam/vg/issues/new/choose to report a bug.
Stack trace path: /tmp/vg_crash_w344Gk/stacktrace.txt
Please include the stack trace file in your bug report!
For the second batch, the same crash happened if I specify a lower --distance-limit 150 or --paired-distance-limit 3/4 to avoid vg falling back on single-end mapping:
Loading Distance Index v2
Initializing MinimizerMapper
Loading and initialization: 23.0865 seconds
Mapping reads to "-" (GAM)
--hit-cap 10
--hard-hit-cap 50000
--score-fraction 0.9
--max-min 500
--num-bp-per-min 1000
--exclude-overlapping-min
--max-extensions 800
--max-alignments 8
--cluster-score 50
--pad-cluster-score 0
--cluster-coverage 0.3
--extension-score 1
--extension-set 20
--max-multimaps 1
--distance-limit 150
--paired-distance-limit 2
--rescue-subgraph-size 4
--rescue-seed-limit 100
--rescue-attempts 15
--rescue-algorithm dozeu
Using fragment length estimate: 108.274 +/- 25.8251
vg: src/dozeu_interface.cpp:126: vg::DozeuInterface::graph_pos_s vg::DozeuInterface::calculate_max_position(const vg::DozeuInterface::OrderedGraph&, const vg::DozeuInterface::graph_pos_s&, size_t, bool, const std::vector<const dz_forefront_s*>&): Assertion `forefronts.at(max_node_index)->mcap != nullptr' failed.
ERROR: Signal 6 occurred. VG has crashed. Visit https://github.com/vgteam/vg/issues/new/choose to report a bug.
Stack trace path: /tmp/vg_crash_aN9npU/stacktrace.txt
Please include the stack trace file in your bug report!
5. What data and command can the vg dev team use to make the problem happen?
first batch:
$vg safari -t 32 -p \
--gbz-name graph.d4.gbz \
-d k"$k".w"$w".graph.d4.dist \
-m k"$k".w"$w".safari.min \
-q k"$k".w"$w".safari.ry \
-f "$READS"_R1.fastq.gz \
-f "$READS"_R2.fastq.gz \
--deam-3p "$READS"_BWAaln_sort.bam.3p.prof \
--deam-5p "$READS"_BWAaln_sort.bam.5p.prof \
> "$READS"_safari_k"$k".w"$w".gam
second batch:
$vg safari -t 32 -p \
--gbz-name graph.d4.gbz \
-d k"$k".w"$w".graph.d4.dist \
-m k"$k".w"$w".safari.min \
-q k"$k".w"$w".safari.ry \
-f "$READS"_R1.fastq.gz \
-f "$READS"_R2.fastq.gz \
--deam-3p "$READS"_BWAaln_sort.bam.3p.prof \
--deam-5p "$READS"_BWAaln_sort.bam.5p.prof \
--distance-limit 150 \
> "$READS"_safari_k"$k".w"$w".gam
6. What does running vg version say?
vg version 79cc189
Compiled with g++ (Ubuntu 9.4.0-1ubuntu1~20.04.2) 9.4.0 on Linux
Linked against libstd++ 20210601
Built by jdru@node07
Many thanks!!
Update: I re ran the commands with vg giraffe just to check if the problem was the vg frozen version used for vg safari:
$vg giraffe -t $cpus -p \
--gbz-name 5.1_MC18.d4.gbz \
-d 5.1_MC18.d4.dist \
-m 5.1_MC18.d4.k"$k".w"$w".min \
-f "$READS"_R1.fastq.gz \
-f "$READS"_R2.fastq.gz \
--distance-limit 100 \ #needed for the second batch samples
> "$READS"_giraffe_k"$k".w"$w".gam
I got the same crash but with more details:
Guessing that 5.1_MC18/5.1_MC18.d4.xg is XG
Preparing Indexes
Loading Minimizer Index
Loading GBZ
Loading Distance Index v2
Paging in Distance Index v2
Initializing MinimizerMapper
Loading and initialization: 19.8915 seconds
Of which Distance Index v2 paging: 0.453282 seconds
Mapping reads to "-" (GAM)
--watchdog-timeout 10
--match 1
--mismatch 4
--gap-open 6
--gap-extend 1
--full-l-bonus 5
--max-multimaps 1
--hit-cap 10
--hard-hit-cap 500
--score-fraction 0.9
--max-min 500
--num-bp-per-min 1000
--distance-limit 100
--max-extensions 800
--max-alignments 8
--cluster-score 50
--pad-cluster-score 20
--cluster-coverage 0.3
--extension-score 1
--extension-set 20
--rescue-attempts 15
--max-fragment-length 2000
--paired-distance-limit 2
--rescue-subgraph-size 4
--rescue-seed-limit 100
--chaining-cluster-distance 100
--precluster-connection-coverage-threshold 0.3
--min-precluster-connections 10
--max-precluster-connections 50
--max-lookback-bases 100
--min-lookback-items 1
--lookback-item-hard-cap 15
--chain-score-threshold 100
--min-chains 1
--chain-min-score 100
--max-chain-connection 100
--max-tail-length 100
--max-dp-cells 16777216
--rescue-algorithm dozeu
Using fragment length estimate: 105.157 +/- 29.3463
vg: src/dozeu_interface.cpp:126: vg::DozeuInterface::graph_pos_s vg::DozeuInterface::calculate_max_position(const vg::DozeuInterface::OrderedGraph&, const vg::DozeuInterface::graph_pos_s&, size_t, bool, const std::vector<const dz_forefront_s*>&): Assertion `forefronts.at(max_node_index)->mcap != nullptr' failed.
━━━━━━━━━━━━━━━━━━━━
Crash report for vg v1.61.0 "Plodio"
Stack trace (most recent call last) in thread 2920938:
#20 Object "", at 0xffffffffffffffff, in
#19 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0x21cd5bf, in __clone3
#18 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0x2126aaa, in start_thread
#17 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0x20c933d, in gomp_thread_start
#16 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0x20cbc87, in gomp_team_barrier_wait_end
#15 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0x20c358a, in gomp_barrier_handle_tasks
#14 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0xed419a, in unsigned long vg::io::paired_for_each_parallel_after_wait<vg::Alignment>(std::function<bool (vg::Alignment&, vg::Alignment&)>, std::function<void (vg::Alignment&, vg::Alignment&)>, std::function<bool ()>, unsigned long) [clone ._omp_fn.1]
#13 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0xd37031, in std::_Function_handler<void (vg::Alignment&, vg::Alignment&), main_giraffe(int, char**)::{lambda()#1}::operator()() const::{lambda(vg::Alignment&, vg::Alignment&)#6}>::_M_invoke(std::_Any_data const&, vg::Alignment&, vg::Alignment&)
#12 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0x11a0516, in vg::MinimizerMapper::map_paired(vg::Alignment&, vg::Alignment&, std::vector<std::pair<vg::Alignment, vg::Alignment>, std::allocator<std::pair<vg::Alignment, vg::Alignment> > >&)
#11 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0x119f1a1, in vg::MinimizerMapper::map_paired(vg::Alignment&, vg::Alignment&)
#10 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0x119237a, in void vg::MinimizerMapper::process_until_threshold_c<double>(unsigned long, std::function<double (unsigned long)> const&, std::function<bool (unsigned long, unsigned long)> const&, double, unsigned long, unsigned long, vg::LazyRNG&, std::function<bool (unsigned long)> const&, std::function<void (unsigned long)> const&, std::function<void (unsigned long)> const&) const [clone .constprop.0]
#9 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0x11a322b, in vg::MinimizerMapper::map_paired(vg::Alignment&, vg::Alignment&)::{lambda(unsigned long)#16}::operator()(unsigned long) const
#8 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0x11a20aa, in vg::MinimizerMapper::attempt_rescue(vg::Alignment const&, vg::Alignment&, vg::VectorView<vg::MinimizerMapper::Minimizer> const&, bool)
#7 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0xec1282, in vg::Aligner::align_xdrop(vg::Alignment&, handlegraph::HandleGraph const&, std::vector<handlegraph::handle_t, std::allocator<handlegraph::handle_t> > const&, std::vector<vg::MaximalExactMatch, std::allocator<vg::MaximalExactMatch> > const&, bool, unsigned short) const
#6 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0xfe2931, in vg::DozeuInterface::align(vg::Alignment&, handlegraph::HandleGraph const&, std::vector<handlegraph::handle_t, std::allocator<handlegraph::handle_t> > const&, std::vector<vg::MaximalExactMatch, std::allocator<vg::MaximalExactMatch> > const&, bool, signed char, unsigned short)
#5 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0xfe040f, in vg::DozeuInterface::calculate_max_position(vg::DozeuInterface::OrderedGraph const&, vg::DozeuInterface::graph_pos_s const&, unsigned long, bool, std::vector<dz_forefront_s const*, std::allocator<dz_forefront_s const*> > const&)
#4 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0x20f52a5, in __assert_fail
#3 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0x5f26db, in __assert_fail_base.cold
#2 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0x5f27b3, in abort
#1 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0x20fb8b5, in raise
#0 Object "/lustre/fs5/vgl/scratch/ssecomandi/BIN/cactus-bin-v2.9.3/bin/vg", at 0x212840c, in __pthread_kill
ERROR: Signal 6 occurred. VG has crashed. Visit https://github.com/vgteam/vg/issues/new/choose to report a bug.
Please include this entire error log in your bug report!
This probably shouldn't happen. We might still not really be sure what we mean by "fragment distance". We might in some places think of it as distance between the two reads (in which case it can in theory be very small or even negative when reads overlap), and in other places think of it as length of the molecule containing the two reads (in which case it can never be smaller than the longest read).
Somehow we're pulling out a rescue subgraph for some read that's not acceptable to Dozeu when we go to align the other read to it, and makes Dozeu crash.
Did the crash log include any lines about the "context" that the crash happened in? Those should include the names of the read pairs that were in the process of being mapped when the program crashed. You can fetch out those reads by name (zcat whatever_R1.fastq.gz | grep -A3 whatever_read_name > whatever_subset_r1.fastq and similarly for R2), and run Giraffe again on the subset files single-threaded (-t1) and isolate the failing read pair.
If you can provide the read pair and the graph and the k/w settings, we should be abler to replicate and fix the problem.
Hi Adam, I am sorry for the late reply. I'll replicate the error and extract the "problematic" reads early next week. Meanwhile, how could I share the graph and the reads with you? I would prefer via email, I can send you a GoogleDrive link or similar. The .gbz file is 3 Gb.
Many thanks!
Hi Adam, unfortunately there is no such information in the crash file (like which reads caused the crash). I already have the -p flag in vg giraffe/safari.
I am also trying to map paired-end reads with another species and another graph (for a course I'm preparing) and it's falling back to single read mapping. I generated a small pangenome with MC for chr22 of two bird individuals and I wanted to map 7x modern short reads with vg giraffe (v1.67.0 "Vetria", from conda):
vg giraffe -t 32 -p \
--gbz-name 2_bTaeGut_pangenome/bTaeGut_pangenome.d1.gbz \
-m 2_bTaeGut_pangenome/bTaeGut_pangenome.d1.min \
-d 2_bTaeGut_pangenome/bTaeGut_pangenome.d1.dist \
-f 4_short_reads/SRR16569049_1.fastq.gz \
-f 4_short_reads/SRR16569049_2.fastq.gz \
-N wildtype10 \
> 5.1_vg_giraffe/SRR16569049_vg_giraffe.gam
Guessing that 2_bTaeGut_pangenome/bTaeGut_pangenome.d1.xg is XG
Rebuilding minimizer index to include zipcodes
Preparing Indexes
[IndexRegistry]: Constructing minimizer index and associated zipcodes.
use parameters -k 29 -w 11
Loading Minimizer Index
Loading Zipcodes
Loading GBZ
Loading Distance Index v2
Paging in Distance Index v2
Initializing MinimizerMapper
Loading and initialization: 2.82146 seconds
Of which Distance Index v2 paging: 0.0087179 seconds
Mapping reads to "-" (GAM)
--watchdog-timeout 10
--batch-size 512
--match 1
--mismatch 4
--gap-open 6
--gap-extend 1
--full-l-bonus 5
--max-multimaps 1
--hit-cap 10
--hard-hit-cap 500
--score-fraction 0.9
--max-min 500
--min-coverage-flank 250
--num-bp-per-min 1000
--downsample-window-length 18446744073709551615
--downsample-window-count 0
--distance-limit 200
--max-extensions 800
--max-alignments 8
--cluster-score 50
--pad-cluster-score 20
--cluster-coverage 0.3
--max-extension-mismatches 4
--extension-score 1
--extension-set 20
--rescue-attempts 15
--max-fragment-length 2000
--paired-distance-limit 2
--rescue-subgraph-size 4
--rescue-seed-limit 100
--mapq-score-window 0
--mapq-score-scale 1
--zipcode-tree-score-threshold 50
--pad-zipcode-tree-score-threshold 20
--zipcode-tree-coverage-threshold 0.3
--zipcode-tree-scale 2
--min-to-fragment 4
--max-to-fragment 10
--max-direct-chain 0
--gapless-extension-limit 0
--fragment-max-graph-lookback-bases 300
--fragment-max-graph-lookback-bases-per-base 0.03
--fragment-max-read-lookback-bases 18446744073709551615
--fragment-max-read-lookback-bases-per-base 1
--max-fragments 18446744073709551615
--fragment-max-indel-bases 2000
--fragment-max-indel-bases-per-base 0.2
--fragment-gap-scale 1
--fragment-points-per-possible-match 0
--fragment-score-fraction 0.1
--fragment-max-min-score 1.79769e+308
--fragment-min-score 60
--fragment-set-score-threshold 0
--min-chaining-problems 1
--max-chaining-problems 2147483647
--max-graph-lookback-bases 3000
--max-graph-lookback-bases-per-base 0.3
--max-read-lookback-bases 18446744073709551615
--max-read-lookback-bases-per-base 1
--max-indel-bases 2000
--max-indel-bases-per-base 0.2
--item-bonus 0
--item-scale 1
--gap-scale 1
--points-per-possible-match 0
--chain-score-threshold 100
--min-chains 4
--min-chain-score-per-base 0.01
--max-min-chain-score 200
--max-skipped-bases 0
--max-chains-per-tree 1
--max-chain-connection 100
--max-tail-length 100
--max-dp-cells 18446744073709551615
--max-tail-gap 18446744073709551615
--max-middle-gap 18446744073709551615
--max-tail-dp-length 30000
--max-middle-dp-length 2147483647
--wfa-max-mismatches 2
--wfa-max-mismatches-per-base 0.1
--wfa-max-max-mismatches 20
--wfa-distance 10
--wfa-distance-per-base 0.1
--wfa-max-distance 200
--min-unique-node-fraction 0
--rescue-algorithm dozeu
warning[vg::giraffe]: Encountered 100000 ambiguously-paired reads before finding enough
unambiguously-paired reads to learn fragment length distribution. Are you sure
your reads are paired and your graph is not a hairball?
warning[vg::giraffe]: Finalizing fragment length distribution before reaching maximum sample size
mapped 299 reads single ended with 100000 pairs of reads left unmapped
mean: 0, stdev: 1
warning[vg::giraffe]: Cannot cluster reads with a fragment distance smaller than read distance
Fragment length distribution: mean=0, stdev=1
Fragment distance limit: 2, read distance limit: 200
warning[vg::giraffe]: Falling back on single-end mapping
Using fragment length estimate: 0 +/- 1
Mapped 47931638 reads across 32 threads in 135.615 seconds with 12.6935 additional single-threaded seconds.
Mapping speed: 11012.7 reads per second per thread
Used 3866.35 CPU-seconds (including output).
Achieved 12397.1 reads per CPU-second (including output)
Used 18920245329548 CPU instructions (not including output).
Mapping slowness: 0.394734 M instructions per read at 4893.56 M mapping instructions per inclusive CPU-second
Memory footprint: 1.05947 GB
What should I do in these situations?
If you can't get vg to log the name of the read pair that is crashing it, and you also can't transfer the whole read file, you could try cutting the file in half repeatedly to binary search for the offending read (or at least to get a small file that triggers the problem).
If vg giraffe is not able to learn your fragment length distribution from the data, you can provide the fragment length distribution information yourself with the --fragment-mean and --fragment-stdev parameters. They don't need to be exactly right, but the real fragment lengths need to be plausible under the normal distribution they describe.
(You also need to use these parameters when binary searching for an offending read that causes a crash, because as you get down to smaller and smaller files, vg giraffe will stop having enough data to learn the fragment length distribution, and will stop crashing at rescue because it will start falling back to single-end mapping. Your example crashing run said it got 105.157 +/- 29.3463 for the fragment length distribution, so you would pass --fragment-mean 105.157 --fragment-stdev 29.3463).