7

I am trying to make local databases for the ncbi-blast+ package (version 2.60). I am doing so for 4 T-cell receptor genes. 3 of the 4 (TRAV, TRAJ, TRBV) have worked fine, but I am having problems with TRBJ. I say this because they should be the same format and if one works the rest works.

They take the format below:

>K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48 nt|3| | | |15 AA|15+0=15| | |
NTEAFFGQGTRLTVV

I used the following command (like I say, it was fine for the rest, oddly enough):

makeblastdb -in TRBJ_AA.fa -parse_seqids -dbtype prot

and I get the following error message:

BLAST options error: TRBJ_AA.fa does not match input format type, default input type is FASTA

One issue I checked was the length of the sequences being too short, but I believe 11 is the minimum length required (for inputs, don't know if that is the same for databases)- but would someone be able to confirm, that strikes me as the most likely candidate?

I'm happy to use some other metric, but I was hoping on keeping BLAST for consistency's sake for the time being.

Thank you.

edit: I know it's to late, but for reference I'm on Ubuntu 16.04.2 LTS x86_64 for all steps.

TW93
  • 449
  • 3
  • 11
  • 1
    What operating system are you running this on and what operating system was the input file created on? Also, please [edit] your question and show us the actual file you are using. – terdon Sep 21 '17 at 13:32

1 Answers1

8

This looks like a bug in makeblastdb. Removing the final | | | from your sequence's description makes it work:

>K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48 nt|3| | | |15 AA|15+0=15

So does removing everything after the 1st space:

>K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48

There is nothing in the FASTA definition that would explain this and, even more weirdly, the presence of the | | | alone isn't enough to trigger the issue. I played around with this a little and it looks like it doesn't like the combination of + or -along with empty | |, but there's no clear pattern I can see:

These fail:

>| | |15-| | |
>| | |15+| | |
>| | |15-1| | |
>| | |15+1| | |

These work:

>| |15+| |
>| | |151| | |

So yeah, sounds like a bug. I suggest you report it to the NCBI devs. In the meantime, a simple fix would be to just keep the first space-separated field of your fasta header. Pass your fasta file through this:

awk '{print $1}' fasta.fa > new.fa

That will print the 1st field of each line and, assuming you don't have spaces in your sequence, it should effectively sort this out for you.

terdon
  • 10,071
  • 5
  • 22
  • 48
  • 1
    I'm not convinced that this is a bug per se -- FASTA headers can use any combination of these specific identifiers. It's parsing it somehow (assuming it's figuring out which of the possible identifiers to use based the parsing the + and - as strings), but it must be getting confused by the addition of three empty columns. – OrdiNeu Sep 21 '17 at 14:34
  • 1
    @OrdiNeu what it is after > is left for the user to decide if NCBI decides to add these identifiers or more is their own decision but it should be able to work with other identifiers – llrs Sep 21 '17 at 14:46
  • 1
    @OrdiNeu the FASTA format itself allows anything you care to put there. The identifiers you linked to are just what NCBI tends to use in its FASTA files, but they have nothing to do with the FASTA format. If a program expects a specific format, then the program is buggy. The fact that there's no clear pattern and that the issue occurs both with and without the -parse_seqids flag is further indication of a bug. – terdon Sep 21 '17 at 14:50
  • 1
    Oddly I found on my system that whether or not it 'works' also depends on the sequence lines. For example, take header >| | |151| | | and sequence ATGAAGAG it fails, with ATGAAGAG\nATG it works (in both cases running makeblastdb -dbtype nucl -in my_fasta.fa) – Chris_Rands Sep 21 '17 at 15:10
  • 1
    @Chris_Rands curiouser and curiouser. But that makes it even more likely that this is a bug. – terdon Sep 21 '17 at 15:11
  • 1
    @terdon Yes, can you replicate what I described? – Chris_Rands Sep 21 '17 at 15:12
  • @terdon That's ahem... interesting! The 3 empty columns can be found in the other IMGT sequence headers, which were parsed successfuly- whereas the numerics with +/- are only found in the problem file which does place the blame firmly at their door. It also still doesn't work when leaving out the -parseseq_id flag. Will report it. Thanks for maintaining my sanity. – TW93 Sep 21 '17 at 15:18
  • 1
    @Chris_Rands yes and no. The one you say works, kinda works on my system but complains about BLAST Database creation error: Defline lacks a proper ID around line 1. Changing the defline to >a| | |151| | |, however, reproduces what you describe perfectly: sequence ATGAAGAG fails but sequence ATGAAGAG\nATG works. Weird. – terdon Sep 21 '17 at 15:26
  • @TW93 you're welcome. I know how you feel. I'd have wondered about my own sanity too had I encountered this in the wild. – terdon Sep 21 '17 at 15:27