converting coordinates to sequences using bedtools getfasta - [segmentation fault] errors

Question

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...

Answer 1

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.