I have 600,000 + unique protein sequences that I want to blast against themselves. For example, I want to run sequence 1 against sequence 1, sequence 2 against sequence 2, ... sequence 3453 against sequence 3453 etc.
I have local blast installed on my computer and blastdb is prepared for the protein sequences. Problem is, each sequence has to be blast against the whole database when I will only retain the hit from a unique sequence against itself. Needless to say, I need a much more efficient method! Any idea?
4 answers
curl -ks "ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz" | gunzip -c |\
awk '/^>/{printf("\n%s\n",$0);next;} {printf("%s",$0);} END{printf("\n");}' | grep -v '^$' |\
while read T
do
echo "$T" > prot.fa
read S
echo "$S" >> prot.fa
blastp -query prot.fa -subject prot.fa
done
•
written
by
Pierre Lindenbaum
Just iterate over your sequences and sum the diagonal element of substitution matrix (http://homepages.ulb.ac.be/~dgonze/TEACHING/stat_scores.pdf) for each amino acid, apply some basic math and you'll get the result (I suspect you are looking for a way to compute E-values).
Thanks, that is my goal, so doing a blast is not something I look for as I already know my sequence are a perfect match and I don't need the evalue. I manage to find the different scoring matrices from NCBI, and I made a script to calculate raw score as suggested Mikhail. Now yes, I'll put in the length parameter to normalize it and I just have to find out how I can get the K and lambda parameters right.
Hi,
If I understood the problem correctly what I would do is to first create the database and blast one sequence against it using the default output format option. In that output file (if I remember correctly) there should be information about the effective database length (EFSIZE). Next, I would extract that number and then create a script like this one:
perl -e 'for(1..600000){system("blastp -query in.fa.$_ -subject in.fa.$_ -dbsize EFSIZE -out bout.$_ -evalue ....")}'
where:
in.fa.1 - single sequence
Of course this means that the initial database should be divided accordingly.
I must admit that I haven't done any real blasting in a while so there is a chance I might be missing something here.
cheers
mxs
Use a script to wrap a call to blastp for each protein using the two sequences form of blast: blastp -query ... -subject ...
•
written
by
Jean-Karim Heriche
For the record, this way of using blast doesn't require a database. See Pierre Lindenbaum's post for a development of what I meant.
Log in to answer this question.
-
Comparing alignments of differing length - what's the right test?
written by ngarber •
-
Using 'blastn' to do a remote search on NCBI BLAST against the nt collection database
written by maciwuk •
-
Blastp command line to get only full matches of query sequence
written by rlee03 •
-
Obtaining BLAST results with only one subject per query and only one query per subject
written by tal.shalev •
-
Is it possible that sequence have blast hit in arabidopsis database and no blast hit for viridiplan…
written by tcf.hcdg
-
problem running Blast
written by nagem7 •
-
local or remote blast BLASTp against pdb database (formatdb)
written by money11gupta •
-
How To Blast Against A Local Database
written by mealbahrani •
-
Blast Multifasta Against Ncbi Nr Database, And Download All The Results
written by cristianrohr768 •
-
Blastp Against Human Proteome
written by Sophie •
I'm curios of what kind of results you're expecting, as clearly the best hit in every case will have 100% similarity covering the entire query sequence length.