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:
- List the files in the directory
- Run
bedtools getfasta
to subset a fasta file. Bedtools reports that an index file did not exist, so it is generating a new one.
- List the files in the directory. Note that there is a new file:
sampleTagsNew.fa.fai
, which is the fasta file index.
- Modify the sample tags fasta file [to change the name of a sequence]
- 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.
- Remove the fasta index file
- 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.
C
thing. I don't do eukaryotes (mostly). I am aware of the strong reputation of bedtools, but thats exactly whatC
does if the data, code and RAM don't mix $\endgroup$