blast only specific sequences from a large file

Question

[Kate]

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?

[5heikki]

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.

Answer 1

[Pierre Lindenbaum]

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

[Kate]

All right this works! On a test file, my previous method (blasting against all database) took 17 minutes on 24 cpu, yours takes half of the time (8 minutes). I thought it would go much much quicker, but well...

Answer 2

[mikhail.shugay]

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).

[Kate]

Thanks, that would be a good idea except that I wouldn't know how to start computing this... I want to get the normalized score out of the sequences to get an idea of the highest score they can generate.

[5heikki]

Pretty sure you can calculate the scores without doing a single blast utilizing only the information of sequence lengths, although there's of course database size parameter which will affect at least e-values, but maybe not bitscores (or at least raw scores)

[Kate]

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.

Answer 3

[mxs]

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

[Kate]

Thanks, this works too. As I only need the normalized score, running blastp without any other information works well. The speed is comparable to what I wrote on Pierre comment. thanks!

Answer 4

[Jean-Karim Heriche]

Use a script to wrap a call to blastp for each protein using the two sequences form of blast: blastp -query ... -subject ...

[Kate]

But that would mean creating 600,000 database files (makeblastdb) prior to running blast... It would be quite heavy.

[SES]

No, create one database with all your sequences, then blast the reads against that. This is very common. You would then filter out self hits and parse out whatever you are interested in from the report (it's not clear exactly what that is though).

[Kate]

That is what I did and I had to give up as each sequence is compared to hundreds of thousands which is a loss of time (it took 18 hours @ 24 cpu).

[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.