Segfault in VCF loading
Closed this issue · 5 comments
Hi,
I was trying to use merfin with a deepvariant vcf to explore polishing a hifi assembly. When loading the variants, there was a segfault. I've posted the stacktrace below, where it seems the issue is a gt variable is empty and causes the segfault.
From what I could tell, it happened immediately, where ii was on its first iteration in L225 of vcf.c. It is probably related, but it appears that most records were rejected, as below
Loaded 43 records with 43 unique contig(s) from 0 contig IDs.
while excluding 460383 invalid records
The first few lines of the variants were
NC_006853.1_RagTag 2226 . GCC G 1.9 RefCall . GT:GQ:DP:AD:VAF:PL ./.:4:11:8,2:0.181818:0,2,30
NC_006853.1_RagTag 2468 . CG C 0 RefCall . GT:GQ:DP:AD:VAF:PL 0/0:32:13:10,3:0.230769:0,31,47
NC_006853.1_RagTag 18569 . GC G 35.7 PASS . GT:GQ:DP:AD:VAF:PL 0/1:7:104:51,39:0.375:34,0,5
stacktrace
0x0000000000414880 in gtAllele::parseGT (this=0x3e71a50, gt=0x0,
ref=0x693080 "\005\331\030\203\276N\251\263\252\373\214?\002@\266J\310+\036\023e\243\070\235$\260\256z\225p5ܐ\235\023o\016\022\032Ā\214\274\241\204",
alts=0x676b60, alleles=0x3e71a70, isValid=@0x6c33f0: true)
at merfin/vcf.C:59
#2 0x0000000000414817 in gtAllele::gtAllele (this=0x3e71a50, record=0x6c3380)
at merfin/vcf.C:43
#3 0x000000000041618d in posGT::posGT (this=0x3d2d000, record=0x6c3380)
at merfin/vcf.H:138
#4 0x000000000041560c in vcfFile::loadFile (this=0x664650,
inName=0x7fffffffcb15 "../../DV_test/output/asm.output.vcf.gz")
at merfin/vcf.C:228
#5 0x00000000004150b6 in vcfFile::vcfFile (this=0x664650,
inName=0x7fffffffcb15 "../../DV_test/output/asm.output.vcf.gz")
at merfin/vcf.C:154
#6 0x0000000000408b05 in main (argc=18, argv=0x7fffffffc418)
at merfin/merfin.C:985
Oddly, after uncompressing the vcf file and re-running, I know get a std::out_of_range error in the varMer.C file. It seems like an issue if there is nothing in kstrs, so calling kstrs.at(0) throws the out of range.
This attempt loaded the variants seemingly correctly with
Loaded 7017436 records with 683 unique contig(s) from 712 contig IDs.
while excluding 0 invalid records`
I also tried recompressing the vcf file with bcftools and gzip, and each time it gave different values for the contigs with most records invalid. Not sure what is causing that.
I cut the stacktrace ignoring the std lib allocator fails, but the relevant part is
Generating fasta index.
Processing 'NC_006853.1_RagTag'
terminate called after throwing an instance of 'std::out_of_range'
...
0x000000000041bf5b in varMer::getAvgAbsK (this=0x1a8a2ed80, idx=0)
at merfin/varMer.C:399
#10 0x00000000004068bc in varMers (
seqName=0x7fffffffcab6 "asm.scaffolds.fasta", sfile=0x1a8a02b10,
vfile=0x664650, rlookup=0x65fdd0, alookup=0x65ffb0,
out=0x7fffffffcb46 "test_merfin.fa", comb=15, nosplit=false,
copyKmerDict=..., bykstar=false, threads=8) at merfin/merfin.C:617
#11 0x0000000000408b98 in main (argc=18, argv=0x7fffffffc438)
at merfin/merfin.C:988
Hello @ASLeonard, thanks for tracking this down.
Could you try with an uncompressed vcf, filtering out RefCalls and try again?
Filtering out refcalls fixes the problem, it successfully generates the debug and polish.vcf files. Thanks for that tip.
I tried again by compressing the filtered vcf, and ran into the same weird issues of near random contig numbers and mostly invalid records.
I just tried out the new commits, and it now works correctly reading compressed vcf. I also tried using the original unfiltered vcf (containing RefCalls), and this also worked fine. The updated meryl version is great as well.
It is out of scope for this issue, but the varMer section still doesn't appear to be threading correctly. I added back the private(varMerId) to see if that was causing the serialisation, but it didn't fix it. I added an extra call here to see which thread was being used, and it does use different threads, but only in serial.
Processing 'ptg000283l' on thread 6
Processing 'ptg000604l' on thread 12
Processing 'ptg000502l' on thread 10
Processing 'ptg000457l' on thread 9
Processing 'ptg000606l' on thread 12
I'll hold off on opening a new issue for this as it looks like a work in progress already.
Thanks!