Segmentation fault on getting split reads from long read file?

[cmdcolin]

Hi there, just following up on #42 I received this crash

 samtools view -h reads_lr_skbr3.fa_ngmlr-0.2.3_mapped.readsort.bam|/home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster --ignoreUnmated -s samp.split-sam -o /dev/null
samblaster: Version 0.1.24
samblaster: Inputting from stdin
samblaster: Opening /dev/null for write.
samblaster: Opening samp.split-sam for write.
samblaster: Loaded 86 header sequence entries.
Segmentation fault (core dumped)

The segmentation fault happens after the samp.split-sam file is 3.2Gb

I ran it again with valgrind and got this

ubuntu@ip-172-31-43-125:/data2$ samtools view -h reads_lr_skbr3.fa_ngmlr-0.2.3_mapped.readsort.bam| valgrind /home/linuxbrew/.linuxbrew/Cellar/samblaster/0
.1.24/bin/samblaster --ignoreUnmated -s samp.split-sam -o /dev/null
==6489== Memcheck, a memory error detector
==6489== Copyright (C) 2002-2017, and GNU GPL'd, by Julian Seward et al.
==6489== Using Valgrind-3.13.0 and LibVEX; rerun with -h for copyright info
==6489== Command: /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster --ignoreUnmated -s samp.split-sam -o /dev/null
==6489==
samblaster: Version 0.1.24
samblaster: Inputting from stdin
samblaster: Opening /dev/null for write.
samblaster: Opening samp.split-sam for write.
samblaster: Loaded 86 header sequence entries.
==6489== Invalid read of size 4
==6489==    at 0x40524D: hashTableInsert(hashTable*, unsigned long) (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==    by 0x4038A4: markDupsDiscordants(splitLine*, state_struct*) (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==    by 0x404216: processSAMBlock(splitLine*, state_struct*) (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==    by 0x401738: main (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==  Address 0x1005bec998 is 1,927,528 bytes inside a block of size 2,367,424 in arena "dinfo"
==6489==
==6489== Invalid read of size 4
==6489==    at 0x405250: hashTableInsert(hashTable*, unsigned long) (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==    by 0x4038A4: markDupsDiscordants(splitLine*, state_struct*) (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==    by 0x404216: processSAMBlock(splitLine*, state_struct*) (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==    by 0x401738: main (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==  Address 0x1005bec99c is 1,927,532 bytes inside a block of size 2,367,424 in arena "dinfo"
==6489==
==6489== Invalid read of size 8
==6489==    at 0x405272: hashTableInsert(hashTable*, unsigned long) (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==    by 0x4038A4: markDupsDiscordants(splitLine*, state_struct*) (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==    by 0x404216: processSAMBlock(splitLine*, state_struct*) (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==    by 0x401738: main (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==  Address 0x1005bec990 is 1,927,520 bytes inside a block of size 2,367,424 in arena "dinfo"
==6489==
==6489== Invalid read of size 8
==6489==    at 0x405318: hashTableInsert(hashTable*, unsigned long) (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==    by 0x4038A4: markDupsDiscordants(splitLine*, state_struct*) (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==    by 0x404216: processSAMBlock(splitLine*, state_struct*) (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==    by 0x401738: main (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==  Address 0x4466d9570f662424 is not stack'd, malloc'd or (recently) free'd
==6489==
==6489==
==6489== Process terminating with default action of signal 11 (SIGSEGV)
==6489==  General Protection Fault
==6489==    at 0x405318: hashTableInsert(hashTable*, unsigned long) (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==    by 0x4038A4: markDupsDiscordants(splitLine*, state_struct*) (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==    by 0x404216: processSAMBlock(splitLine*, state_struct*) (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==    by 0x401738: main (in /home/linuxbrew/.linuxbrew/Cellar/samblaster/0.1.24/bin/samblaster)
==6489==
==6489== HEAP SUMMARY:
==6489==     in use at exit: 14,692,245 bytes in 812 blocks
==6489==   total heap usage: 1,855 allocs, 1,043 frees, 28,039,720 bytes allocated
==6489==
==6489== LEAK SUMMARY:
==6489==    definitely lost: 0 bytes in 0 blocks
==6489==    indirectly lost: 0 bytes in 0 blocks
==6489==      possibly lost: 0 bytes in 0 blocks
==6489==    still reachable: 14,692,245 bytes in 812 blocks
==6489==         suppressed: 0 bytes in 0 blocks
==6489== Rerun with --leak-check=full to see details of leaked memory
==6489==
==6489== For counts of detected and suppressed errors, rerun with: -v
==6489== ERROR SUMMARY: 4 errors from 4 contexts (suppressed: 0 from 0)
Segmentation fault (core dumped)

Not very good with reading stack traces like this but thought it could help

[cmdcolin]

File https://s3.amazonaws.com/jbrowse.org/genomes/hg19/skbr3/reads_lr_skbr3.fa_ngmlr-0.2.3_mapped.readsort.bam (readsorted version of a file here http://schatz-lab.org/publications/SKBR3/)

[GregoryFaust]

I'll look into this as soon as the file finishes downloading. I also use valgrind in these situations, so we agree on use of tools. Having the file is even better!

[cmdcolin]

Thanks. I am currently doing some perl crimes to sidestep this...just get all duplicate read names. https://gist.github.com/cmdcolin/7d4ab1fe0c8dd2a5fa01b3e059063b17

[GregoryFaust]

I have found the problem. It is in fact the same issue as #26 but on steroids. Basically in samblaster 0.1.24, I fixed the problem in #26 by padding the "address space" for alignment start/end positions at the beginning/end of every contig in the reference genome the file was aligned against. This guards against reads that have a large clip and are aligned at the very beginning/end of a contig. In 0.1.24, I added enough padding for any reasonable Illumina pair-end read. However, in your data, the (at least first) read that caused the segfault was aligned near the start of chrom 1, and started with a clip of length 15,891! This resulted in a calculated alignment start position outside the actual genome, which in turn led to the segfault. In the upcoming release I will add a more general solution to this issue, as well as a guard that will fail more gracefully if this were to happen again.

In the mean time, and I should have suggested this in the first place, as you are not trying to mark duplicates, you should use the --acceptDupMarks (or just -a) option to the samblaster command to pull splitters. This bypasses the genome location mapping needed for finding duplicates and therefore will avoid the problematic code segment. Also, it will run faster and with less memory.

Thanks,
Greg

[cmdcolin]

Thanks for catching and debugging all this!

[cmdcolin]

I'll try re-running

[GregoryFaust]

Added the --maxReadLength parameter in samblaster release 0.1.25 to fix this problem. One can now specify as long a read length as necessary to avoid having an alignment whose calculated coordinate falls outside of the reference. In addition, I added an error check that will cause a fatal error (instead of a segfault) in such cases with output telling you the length of the read causing the problem. Finally, reads that are longer than --maxReadLength can have calculated coordinates that are within the reference as a whole, yet fall outside the reference sequence (chrom) in which they aligned. In such cases, there is a small chance that they could be marked as a duplicate of an alignment in an adjacent reference sequence (chrom). To help with such cases, samblaster now keeps track of the number of reads longer than --maxReadLength, and the maximum length of such reads, and outputs both at the end of the run. See the release notes for 0.1.24 for more details.

I could have just added a static humongous --maxReadLength internal to the program with no user visible parameter, but this has too downsides. First, someone will always end up with data longer than the fixed value. Second, having very long --maxReadLength with a reference genome that has a very high number of contigs (for example, a genome that is not fully assembled) would lead to unacceptably high memory consumption.