cufflinks hates my long reads: segmentation faults on bwa output

[dagarfield]

Hi folks,

I've run into what appears to be a memory usage issue in cufflinks. I've been aligning 454 reads against a reference genome using Lastz and BWA. Cufflinks seems to run just fine (or at least runs) on the lastz output. However, I get a segmentation fault on the BWA output.

To see if the issue lies in the output format, I put both files into a similar, 11 column formatted SAM file.

Here are the entries for one particular 454 read in the Lastz output file.

FEKCLA404H56KG	0	Scaffold1	59767	255	29M1D50M178H	*	0	0	AGACGAAGGGCAACGATGCTTACAAGGTCACCGTCTTCACCCGAAGAAGTCGAAGAGACAGAGTCGAGATGCTGAAGGT	*
FEKCLA404H56KG	0	Scaffold1	60102	255	75H38M1I20M1I54M68H	*	0	0	AGGTGCAAATAAAGTGATAAAGAGGCAAAGAAGTACATCCAAAAACCCGAGGAGGGATCAAAAACGGATAAGAAACAGCGAAAAAGAAACAGCCTGATCATCCCTTCTGCAGAG	*

And here is the output for the same read in the BWA data formatted to look like Lastz (i.e. the first ten columns plus *)

FEKCLA404H56KG	0	Scaffold1	59767	27	29M1D50M178S	*	0	0	AGACGAAGGGCAACGATGCTTACAAGGTCACCGTCTTCACCCGAAGAAGTCGAAGAGACAGAGTCGAGATGCTGAAGGTGCAAATAAAGTGATAAAGAGGCAAAGAAGTACATCCAAAAACCCGAGGAGGGATCAAAAACGGATAAGAAACAGCGAAAAAGAAACAGCCTGATCATCCCTTCTGCAGAGAGTCTTAAAGCTGACCATCCCTCTGTAACCAGACTATCAAACGCCAGCGATGAATTAGCAGAGGATCT	*
FEKCLA404H56KG	0	Scaffold1	60102	28	75S38M1I20M1I54M68S	*	0	0	AGACGAAGGGCAACGATGCTTACAAGGTCACCGTCTTCACCCGAAGAAGTCGAAGAGACAGAGTCGAGATGCTGAAGGTGCAAATAAAGTGATAAAGAGGCAAAGAAGTACATCCAAAAACCCGAGGAGGGATCAAAAACGGATAAGAAACAGCGAAAAAGAAACAGCCTGATCATCCCTTCTGCAGAGAGTCTTAAAGCTGACCATCCCTCTGTAACCAGACTATCAAACGCCAGCGATGAATTAGCAGAGGATCT	*
FEKCLA404H56KG	0	Scaffold1	60677	22	187S63M7S	*	0	0	AGACGAAGGGCAACGATGCTTACAAGGTCACCGTCTTCACCCGAAGAAGTCGAAGAGACAGAGTCGAGATGCTGAAGGTGCAAATAAAGTGATAAAGAGGCAAAGAAGTACATCCAAAAACCCGAGGAGGGATCAAAAACGGATAAGAAACAGCGAAAAAGAAACAGCCTGATCATCCCTTCTGCAGAGAGTCTTAAAGCTGACCATCCCTCTGTAACCAGACTATCAAACGCCAGCGATGAATTAGCAGAGGATCT	*

Everything seems more or less the same EXCEPT for the length of the sequences -- In the BWA output, the whole read is written, while in the Lastz output only the section that aligns to the reference is written.

The result of running these files: A segmentation fault. However, if I take only a subset of the lines in the sam file things run just fine. Here is what you see when you run only the first 100,000 lines. For BWA

[bam_header_read] EOF marker is absent.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
File bwa_as_Lastz.sam doesn't appear to be a valid BAM file, trying SAM...
[17:21:25] Inspecting reads and determining fragment length distribution.
> Processed 11701 loci.                        [*************************] 100%
> Map Properties:
>	Total Map Mass: 99150.92
>	Read Type: 35bp single-end
>	Fragment Length Distribution: Gaussian (default)
>	              Estimated Mean: 203.48
>	           Estimated Std Dev: 75.29

The important thing to note: We have 35bp reads and the total Map Mass is less than the number of lines.

Now here's what you get with the Lastz output.

[bam_header_read] EOF marker is absent.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
File lastZ.sam doesn't appear to be a valid BAM file, trying SAM...
[17:23:35] Inspecting reads and determining fragment length distribution.
> Processed 3360 loci.                         [*************************] 100%
> Map Properties:
>	Total Map Mass: 100000.00
>	Read Type: 30bp single-end
>	Fragment Length Distribution: Gaussian (default)
>	              Estimated Mean: 203.00
>	           Estimated Std Dev: 75.72
[17:23:36] Assembling transcripts and estimating abundances.

Now the reads are 30bp and the Total Map Mass is exactly the same as the number of lines!?!

OK, one final bit of oddity. If I run all 13544847 lines of the Lastz output (that's 2.2G total file size), things seem to work. Here's the start of the process

[bam_header_read] EOF marker is absent.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
File allSeqs_noFlankingNs_to2.5_lastz_header_sorted.sam.sorted_noHead doesn't appear to be a valid BAM file, trying SAM...
[17:26:11] Inspecting reads and determining fragment length distribution.
> Processed 478607 loci.                       [*************************] 100%
> Map Properties:
>	Total Map Mass: 13544847.00
>	Read Type: 30bp single-end
>	Fragment Length Distribution: Gaussian (default)
>	              Estimated Mean: 203.00
>	           Estimated Std Dev: 75.72
[17:27:11] Assembling transcripts and estimating abundances.

But when I try to run the first 300,000 (of ~700,000) lines of the BWA output, things fail.

[bam_header_read] EOF marker is absent.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
File bwa_as_Lastz.sam doesn't appear to be a valid BAM file, trying SAM...
[17:29:06] Inspecting reads and determining fragment length distribution.
Segmentation fault

I'm deeply confused here on what to do, and I can't find anything on the web (as far as I can tell, the Lastz output shouldn't run at all missing the XS:A:-/+ field).

Thanks,

DG

[seb567]

If you use the GNU stack:

add the following option to the compilation flags:

-g

then, run the command leading to a segmentation fault in gdb.

[B]seb@godzilla:~$ cat sf.c [/B]
int main(){
	int*i=0;
	*i=0; /* will generate a segmentation fault */
	return 0;
}

[B]seb@godzilla:~$ gcc -g sf.c -o Program[/B]

[B]seb@godzilla:~$ gdb Program [/B]
GNU gdb (GDB) 7.1-ubuntu
Copyright (C) 2010 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.  Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-linux-gnu".
For bug reporting instructions, please see:
<http://www.gnu.org/software/gdb/bugs/>...
Reading symbols from /home/seb/Program...done.

[B](gdb) run[/B]
Starting program: /home/seb/Program 

Program received signal SIGSEGV, Segmentation fault.
0x00000000004004d4 in main () at sf.c:3
3		*i=0;

[B](gdb) bt[/B]
#0  0x00000000004004d4 in main () at sf.c:3

Then send the backtrace to the authors of the software so they can fix it.