<p>I am seemingly stuck with something that should be very simple and I hope I haven't overlooked something obvious. </p>

<p>Question: How can I make a valid Blast-database <em>with Taxids</em> from a NCBI query export? </p>

<p>What I have tried so far:</p>

<p>For a meta-genomics project I need a custom made blast database which I wish to generate from the result of the following N<a href='http://www.ncbi.nlm.nih.gov/nuccore'>CBI Nucleotide query</a>:</p>

<pre><code>Viruses[Organism] AND srcdb_refseq[PROP] NOT cellular organisms [ORGN]
</code></pre>

<p>The result is 3986 entries which I exported and saved (via 'Send to') in FASTA and ASN1 format.  (Both files are seemingly containing the right amount of entries) As this is a meta-genomics project I would love to have the taxon ids in the blastdb.</p>

<p>I was successful with making a valid blast database from the FASTA file using makeblastdb, but the FASTA header doesn't include taxids, hence I tried to make a blast database from the ASN1 export using the following command (it is not clear from <a href='http://www.ncbi.nlm.nih.gov/books/NBK1763/#CmdLineAppsManual.4_User_manual'>the documentation</a> which formats can be used to create the database):</p>

<pre><code>    $ makeblastdb -in AllViralDNARefSeq.asn1 -dbtype nucl  -out ViralASN1 -title  "All Viral RefSeq DNA from NCBI ASN1"

Building a new DB, current time: 12/20/2011 10:37:28
New DB name:   ViralASN1
New DB title:  All Viral RefSeq DNA from NCBI ASN1
Sequence type: Nucleotide
Keep Linkouts: T
Keep MBits: T
Maximum file size: 1073741824B
Adding sequences from FASTA; **added 10 sequences** in 0.00906897 seconds.
</code></pre>

<p>As you can see, this does not work as it adds only 10 sequences.</p>

<p>Any help to get in the Taxonids is appreciated it doesn't have to be elegant, I just need the database from that query. I am using Blast+ 2.2.25 </p>

I think this one requires some scripting.  Here are a few ideas.

First, I'd fetch viral Refseq sequences in Genbank format from the FTP site:

    wget ftp://ftp.ncbi.nih.gov/refseq/release/viral/viral.1.genomic.gbff.gz
    gunzip viral.1.genomic.gbff.gz
    grep -c "^LOCUS" viral.1.genomic.gbff # 3984

Then I'd parse the Genbank file to extract an ID, description and sequence for each entry. Since taxon IDs are contained in a field of the form `/db_xref="taxon:NNNN` where NNNN = taxon ID, they can be extracted and written into the header of a new fasta file.

Some quick and dirty Bioperl to illustrate:

    #!/usr/bin/perl -w
    
    use strict;
    use Bio::SeqIO;
    
    my $file  = "viral.1.genomic.gbff";
    my $seqio = Bio::SeqIO->new(-file => $file, -format => "genbank");
    my $fasta = Bio::SeqIO->new(-file => ">$file.fa", -format => "fasta");
    
    while(my $seq = $seqio->next_seq) {
        my $taxid = "";
        for my $feat($seq->get_SeqFeatures) {
        if($feat->has_tag("db_xref")) {
            for my $id($feat->get_tag_values("db_xref")) {
            if($id =~/taxon:(\d+)/) {
                $taxid = $1;
            }
            }
           }
        }
        my $fa = Bio::Seq->new(-id => $seq->id,
                               -desc => $taxid.", ".$seq->description,
                               -seq => $seq->seq);
        $fasta->write_seq($fa);
    }

This gives a fasta file *viral.1.genomic.gbff.fa*, with headers that look like:

    >NC_003038 176652, Invertebrate iridescent virus 6, complete genome.