Minimap2: samtools mpileup segfault

[fbemm]

Dear Heng,

I jump a bit ahead to see how me can quickly start working with minimap2. I mapped an assembly on a reference genome (-xasm5 -a), sorted and converted the output to a bam file. Then I passed the bam file to samtools mpileup (-v -u -B) to generate a VCF file. Works for the first to chromosomes in my reference but then fails on the third. The alignment that is causing the fail is 44kb long with about 440 CIGAR edits. Other, way longer alignments with more CIGAR edits have been processed on the first two chromosomes. I am not sure what is happening. Might be related to samtools/htslib (btw. I compiled the latest stable with your "Work-around for supporting CIGARs with >64k operations in BAM files") and I am not even sure the issues belongs here. I uploaded the alignment and bot qry and ref here https://owncloud.tuebingen.mpg.de/index.php/s/44pZgooEWvDgboD. Any idea how and where to start looking for the problem?

Thanks a lot especially for minimap2 :+1:
Felix

[lh3]

If samtools can convert SAM to BAM, it is more likely that samtools is buggy. Perhaps the ~20k insertion in your bam is killing mpileup. Generally, if you want to generate VCF from asm-to-asm alignment, I would recommend htsbox:

htsbox pileup -q10 -vcf ref.fa sorted.bam > var.vcf

It is tailored for such alignment. This command line works with your BAM.

[fbemm]

Thanks! htsbox works great. vg and vcflib don't like the output VCF at the moment. Looks like working on calls across or within N-containing regions especially in larger SVs is not yet possible. Either it is a pure size problem or the N's themselves. Pruning all N-containing variants from the VCF does not work atm.

[daviesrob]

I've opened issue samtools/samtools#712 for this.

[senaj]

Hi lh3,
I had a similar issue with samtools mpileup failing but htsbox was the work around. Thank you very much for recommending htsbox!

[lh3]

Samtools-1.7, which has just been released, should give you the correct pileup.