gmap_build error

[StefanF]

Hey guys,

I am trying to build a custom reference (mm9+ERCC spike-in) to reproduce some single cell RNA-seq results. I appended the mm9 genome (as a fasta file) with the ERCC sequences (also as fasta) and run the following command:

gmap_build -d mm9_ERCC92 mm9_ERCC92.fa -k 15 -D ~/stefan/scLVM/

All goes well until:

Building suffix array
SACA_K called with n = 2725848238, K = 5, level 0
/usr/local/bin/gmapindex -d mm9_ERCC92 -F /stefan//scLVM//mm9_ERCC92 -D /stefan/scLVM//mm9_ERCC92 -S failed with return code 9 at /usr/local/bin/gmap_build line 375.

Any ideas what could be going on here?

[Michael.Ante]

Hi Stefan,

only a shot in the dark: did you check the fasta's window length? Try fastx-toolkit's fasta_formatter:

fasta_formatter -w 60 -i mm9_ERCC92.fa -o mm9_ERCC92_w60.fa
diff -s mm9_ERCC92.fa mm9_ERCC92_w60.fa

Cheers,
Michael

[StefanF]

Hi Michael,

I just compared the two files and they are indeed different.

[Michael.Ante]

Thinking about it twice, the two files can differ and still it will not work (if the window-length is constant but not 60 nts).

Just check the line length visually with the first ERCC-transcript (should be ERCC-00002):

grep -A1 -B1 ERCC-00002 mm9_ERCC92.fa

[StefanF]

Okay, the ERCC sequences are shorter than 60 nts. So gmap_build can't handle that? But it should work with the the corrected file mm9_ERCC92_w60.fa?

Edit: Alright, I get the same error with the formatted file... I really appreciate the help though. Maybe you or someone else has another idea?

Edit 2: The error seems to be in the suffix array construction. Is it really necessary to build that or does it only slow things down later?