2
$\begingroup$

I use bedtools to convert genomic coordinates (in the form of bed files, chrX x1:x2) to genomic sequence content in the form of fasta files ('ACGT...G').

The bedtools format for extracting sequence (.fa file) from coordinates (.bed file) is:

bedtools getfasta -fi  input_sequences.fa -bed  genomic_coordinates.bed -fo output_sequences.fa

Very frequently I get a "Segmentation fault (cored dumped)" during the run.

Warning: the index file is older than the FASTA file.
Segmentation fault (core dumped)

The strange thing is that the error will go away, if for example I add additional data (chromosomes) to the input fasta file (input_sequences.fa), OR remove some chromosomes?? I am having trouble finding a logical path to explain this behavior so I can work around it.

There is no changing of fasta files while bedtools is working - in fact there is nothing in the primary directory other than the two files, the fasta file input_sequence.fa and the bed file genomic_coordinates.bed. There is 256GB of RAM memory on the machine, but that's not an issue, especially since the error goes away sometimes when additional chromosomes added to the fasta file.

Anyone experienced issues with bedtools getfasta erratic Segmentation faults? Known workarounds? I am using the current version of bedtools v2.30.0

What are some of the alternatives to "bedtools getfasta" for converting genomic coordinates to sequences?

Thanks...

$\endgroup$
3
  • 1
    $\begingroup$ This sounds like C thing. I don't do eukaryotes (mostly). I am aware of the strong reputation of bedtools, but thats exactly what C does if the data, code and RAM don't mix $\endgroup$
    – M__
    Oct 26, 2022 at 3:17
  • $\begingroup$ Thanks @M__ you are right, the error is likely caused by bedtools accessing (inside it's C/C++ code) un-allocated memory space. It's possible an array get's defined under certain conditions where the expected/allocated size ends up not being correct - often it's a programming oversight - what I am trying to understand is the condition that may be leading up to it, so I can avoid the condition: something related to makeup of the input fasta file as related to the coordinates in the bed file. $\endgroup$
    – Zebra Fish
    Oct 26, 2022 at 3:48
  • 1
    $\begingroup$ Have you tried taking a sample of your bed? Say, for example, splitting it into 3 parts and running the code again. Try to narrow down what is causing that error. Is the bed a big file? Maybe, you don't have enough RAM? $\endgroup$
    – user324810
    Oct 26, 2022 at 7:51

1 Answer 1

1
$\begingroup$

The error about old index files is insightful to me. Bedtools presumably assumes that the input files are static / unchanging. If bedtools is automatically indexing the files before extracting sequences (when index files don't already exist), and the input files keep changing, that would be a good explanation for why problems are happening.

Do these input fasta files have bedtools and/or samtools index files in the same directory (usually fileName.fa.fai). Are you re-indexing the fasta files after changing them (e.g. via samtools faidx fileName.fa)?

I don't mean changing during the running of bedtools, I mean changing between different bedtools runs. If you create a file input_sequences.fa, then the first time it is subset, it will create an index file input_sequences.fa.fai. The next time bedtools is run to subset an input_sequences.fa within the same directory, it will check to see if an index file already exists; if one exists, it won't be regenerated. If the file input_sequences.fa has changed in the interim, then the index file will be incorrect, and the results of subsetting will be unpredictable.

This is why Bedtools reports the error: it can cause strange things to happen with output (intermittent segmentation faults, for example).

In fact, bedtools is explicit about this index generation in its output:

$ ls
sampleTagsNew.bed  sampleTagsNew.fa

$ bedtools getfasta -fi sampleTagsNew.fa -bed sampleTagsNew.bed 
index file sampleTagsNew.fa.fai not found, generating...
>NewSampleTagX01:0-20
GTTGTCAAGATGCTACCGTT
>NewSampleTagX02:1-20
TTGTCAAGATGCTACCGTT

$ ls
sampleTagsNew.bed  sampleTagsNew.fa  sampleTagsNew.fa.fai

$ perl -i -pe 's/NewSampleTagX01/NewSampleTagX01_updated/' sampleTagsNew.fa

$ bedtools getfasta -fi sampleTagsNew.fa -bed sampleTagsNew.bed 
Warning: the index file is older than the FASTA file.
>NewSampleTagX01:0-20
updatedGTTGTCAAGATG
>NewSampleTagX02:1-20
TagX02GTTGTCAAGATG

$ rm sampleTagsNew.fa.fai

$ bedtools getfasta -fi sampleTagsNew.fa -bed sampleTagsNew.bed 
index file sampleTagsNew.fa.fai not found, generating...
WARNING. chromosome (NewSampleTagX01) was not found in the FASTA file. Skipping.
>NewSampleTagX02:1-20
TTGTCAAGATGCTACCGTT

In the example above, I carry out the following operations:

  1. List the files in the directory
  2. Run bedtools getfasta to subset a fasta file. Bedtools reports that an index file did not exist, so it is generating a new one.
  3. List the files in the directory. Note that there is a new file: sampleTagsNew.fa.fai, which is the fasta file index.
  4. Modify the sample tags fasta file [to change the name of a sequence]
  5. Run bedtools getfasta again with an identical command line. It produces the same warning you have seen, and produces unexpected output, i.e. incorrect results.
  6. Remove the fasta index file
  7. Run bedtools getfasta again with an identical command line. This time the output is correct, stating that the first sequence cannot be found in the input fasta file.
$\endgroup$
1
  • $\begingroup$ Thank you @gringer... update is SUPER helpful. Indeed the fasta file contents do change between subsequent bedtools calls, while the input fasta file name stays the same. There IS an *.fai file which I deleted and reran the bedtools call that was consistently causing segmentation fault. This time I got no segmentation fault!!! As you suspected index file *.fai seems to be the issue, but not all outdated *.fai files cause problems. I just tested that by mixing and matching all kinds of outdated *.fai files, no segmentation fault. For now I am happy to manually delete *.fai's. Thanks again!! $\endgroup$
    – Zebra Fish
    Oct 27, 2022 at 3:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.