`vg haplotypes` crashed: top-level chain 7964 is a loop; haplotype sampling cannot be used with this graph

[GooLey1025]

error log from vg haplotype

$ vg haplotypes -v 2 -t 24 -H targets.hapl targets.gbz 
warning[vg haplotypes] Thread count (-t) is greater than the maximum number of threads available (16), capping to that value
Loading GBZ from targets.gbz
[vg haplotypes] Generating haplotype information
[vg haplotypes] Guessing that distance index is targets.dist
[vg haplotypes] Loading distance index from targets.dist
[vg haplotypes] Building minimizer index
Building MinimizerIndex with k = 29, w = 11, payload = none
1846289 keys (1782180 unique)
Minimizer occurrences: 1973556
Load factor: 0.44019
Construction time: 0.656898 seconds
[vg haplotypes] Built the minimizer index in 0.65691 seconds
[vg haplotypes] Guessing that r-index is targets.ri
Loading r-index from targets.ri
HaplotypePartitioner: 362227 fragments for 2978 haplotype sequences
Partitioning parameters:
- target length 10000 bp
- 32 jobs
Determining construction jobs
error[vg haplotypes] HaplotypePartitioner::partition_haplotypes(): top-level chain 7964 is a loop; haplotype sampling cannot be used with this graph
.This can sometimes be resolved by using the vg index -P option to specify a reference backbone when computing the distance index.

My gbz generated procedure:

vg find -x mc-cactus.xg -R target_regions.paths.bed -c 2000 -L \
| vg combine - > targets.vg

vg view -g targets.vg > targets.gfa

vg gbwt \
  -G targets.gfa \
  --max-node 1024 \
  --num-jobs 16 \
  -p \
  -g targets.gbz

vg index -t 24 -j targets.dist --no-nested-distance targets.gbz -P  FaHC_P8_Nipponbare_TEJ
vg gbwt -p -r targets.ri -Z targets.gbz
vg haplotypes -v 2 -t 24 -H targets.hapl targets.gbz 

What I am doing here is to extract the subgraph of MC-whole-genome-graph, then I can do giraffe mapping on the specific regions to save computation time.

I am not sure whehter my -P FaHC_P8_Nipponbare_TEJ provided is right?
I have generated the relevant backbone-reference paths for debugging:

$ vg paths -L -x targets.gbz | grep "FaHC_P8_Nip" | head -n 10
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr1[30749]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr1[163334]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr1[174338]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr1[178544]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr1[183832]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr1[214569]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr1[244971]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr1[299632]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr1[299783]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr1[300236] 
$ vg paths -L -x targets.gbz | grep "FaHC_P8_Nip" | tail -n 10
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr9[22744737]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr9[22765169]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr9[22814343]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr9[22828158]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr9[22834274]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr9[22855416]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr9[22862931]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr9[22888249]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr9[22901706]
FaHC_P8_Nipponbare_TEJ#0#P8.Nipponbare.TEJ.Chr9[22936364]

[jltsiren]

Was FaHC_P8_Nipponbare_TEJ also the sample you used as the reference when building the original graph with Minigraph–Cactus? Using any other sample as the backbone in vg index can produce unpredictable results. At least if your graph is the default/clipped graph rather than the full graph.

Are the intervals in the BED file in FaHC_P8_Nipponbare_TEJ or in another sample?

Haplotype sampling may not work if you have multiple graph components for a chromosome. At least you should rename the reference paths (before building the .hapl) file so that each component has a distinct contig/sequence/locus name. But you could also use Giraffe directly without haplotype sampling (and without building the .hapl file).

If you have whole-genome reads, you should always map them to a whole-genome graph. Otherwise you can get false positive mappings, and I'm not sure you will save any time by using a partial reference.

[GooLey1025]

THanks for your reply.

1

Was FaHC_P8_Nipponbare_TEJ also the sample you used as the reference when building the original graph with Minigraph–Cactus?

Yes, this is the same reference genome as the original graph with Minigraph-Cactus.

2

Are the intervals in the BED file in FaHC_P8_Nipponbare_TEJ or in another sample?

The intervals corrdinates is in FaHC_P8_Nipponbare_TEJ.

3

Haplotype sampling may not work if you have multiple graph components for a chromosome. At least you should rename the reference paths (before building the .hapl) file so that each component has a distinct contig/sequence/locus name. But you could also use Giraffe directly without haplotype sampling (and without building the .hapl file).

Yep, I would rather use haplotypes sampling than directly giraffe mapping, because encountering error when building index:

$ vg autoindex --workflow giraffe -g targets.gfa  -p targets
[vg autoindex] Executing command: vg autoindex --workflow giraffe -g targets2.gfa -p targets
[IndexRegistry] Checking for haplotype lines in GFA.
[vg autoindex] Guessing that targets.dist is Giraffe Distance Index
[IndexRegistry] Constructing a GBZ from GFA input.
[IndexRegistry] Constructing minimizer index and associated zipcodes.
[IndexRegistry]     using parameters -k 29 -w 11
error[VPKG::load_one]: Correct input type not found in targets.dist while loading bdsg::SnarlDistanceIndex

By the way, could you give me some command usage tips on renaming the reference paths? rename it in human-readable gfa file or how to rename it?

4

If you have whole-genome reads, you should always map them to a whole-genome graph. Otherwise you can get false positive mappings, and I'm not sure you will save any time by using a partial reference.

Yes, you are right. We are going to develop a fast-genotyping plugin. When mapping a partial reference, of course there are false positive mappings because of homologus sequences. First We are trying to map against paritial reference, and then extract mapped bam and convert it into fastq. Then use this mapped fastq and re-map against full-reference. Then we can delete the homologus sequences and reduce false positives.

[jltsiren]

If you add option --no-guessing to the vg autoindex command, it will rebuild the distance index instead of trying to reuse the existing one (which doesn't work for some reason).

[adamnovak]

If you want to rename a path, it probably is easiest to do it in the GFA, before you feed it to vg. We don't have a good tool for renaming paths.

[GooLey1025]

Thanks for vg team's help. I appreciate it much!