Hi Biostar,
I am working on a microbial gene annotation project and I am interested in taking a large number of sequences (say 20,000) and blasting them against the NR database. However, even with a local copy of nr and a decent computer this will take forever (days, anyway). My sequences, however, are not random so I wanted to make a subset of NR that will allow me to perform a much faster BLAST query.
My sequences are 454 reads of PCR amplicons made with highly degenerate primers targeting a family of proteins so what I would like to do is as follows:
So before reinventing the wheel, I was wondering if there is a trivially easy way to do this that someone knows about, or has done. In the meantime I will try writing a biopython script to run the blast and to use the 'gi' and 'sseqids' of the hits to make a new fasta file.
thanks, zach cp
•
written
by
Zach Powers
No need for a python script. The command line is your friend.
Assuming your initial sequences are in known_sequences.fa, first do the blastx against a local copy of nr (add your favourite blastx parameters) using the blast+ toolkit:
blastx -query known_sequences.fa -db nr -outfmt 6 >known_sequences.blastx.nr.hits.txt
Then, extract the GIs from the second column of the hits file:
cut -f2 known_sequences.blastx.nr.hits.txt | cut -f2 -d "|" >known_sequences.blastx.nr.gids.txt
Then make a blastdb using blastdb_aliastool:
blastdb_aliastool -dbtype prot -gilist known_sequences.blastx.nr.gids.txt -db nr -out nr_subset
nice Sujal, Thanks! I will ty this. One downside ( i will check) of using the GIs is that many of the blast queries hit back large sequences (whole genomes, ESTs) but I think the speed gains will still be awesome ---- time to give it a go.
I think you should consider ALchEmiXt's broader suggestions though - a taxon restricted database is more general than just picking the ones your sequences hit.
Sujai just a tip, to extract the GI's instead of using two cut commands you could simply use awk:
"awk -F"|" '{print $2}' known_sequences.blastx.nr.hits.txt > known_sequences.blastx.nr.gids.txt
Learning just the basics about awk or sed can really save quite some time on the long run.
Agreed, I was just being lazy ;-)
If you are interested in just a specific taxon you could limit your BLAST searches by taxonid. Furthermore the blast+ utilities will allow you to generate a subset multiFasta or DB from a list of gi's. More specifically the BLASTDBCMD command.
So there are several ways to go with this;
In a completely different application, I need to limit my BLAST results to a particular taxon---say, danio rerio. How do I do this from the command line for the c++ version of blastn?
Oh, found it.
Search the Entrez Protein database (or Nucleotide, etc) with query: "txid7742[ORGN]" or "danio rerio[orgn]
Select "Send to File" and choose format "GI list"
Then, for blastn, us -gilist <file>.
I figured this out from the following thread.
http://www.biostars.org/post/show/6528/vertebrate-subset-nr-database-build-my-own/
Log in to answer this question.
In step 2 "take those sequences", you mean take the hits to those sequences? Otherwise the procedure doesn't make sense.
@neilfws - yes thats what I mean. corrected. thanks