reciprocal blast

[renesh]

How to do reciprocal blast for two datasets?

[tomc]

for two datasets A & B
where records a are in dataset A and records b are in dataset B

Blast query a against target dataset B to obtain best hit b'
Blast query b' against target dataset A to obtain best hit a'

if a == a' then b' is the reciprocal best hit of a (rbh)
otherwise a and b' are not reciprocal best hits

this is done for all a in A.
(and possible for any b in B not a rbh of something in A)

sorry if this too abstract. but the question is not so specific.

[kenietz]

Hi,
firstly i'm sorry for reviving the thread but i have a more specific question on the subject.

Further i will use the notation used by TOMC in previous post.

So i blast query a against set B and get the best hit b'. But then should i blast with b' or the full sequence of b'? I ask cos i found cases where if i blast with b' then a == a' but if blast with the full seq of b' then a != a' (they don't equal). Logically one should blast with the full seq of b' at least in my opinion but im not sure.

Any help is appreciated.
Thank you

[tomc]

Hi kenietz

Since each data set can be one or many sequences of any lenght,
the question remains less than fully specified,
but under the assumption that both datasets are composed of many shorter sequences each which represent some logical unit (transcript, illumima read) then yes
you blast the entire representation of b identified by the hit b' against dataset A.

Under different assumptions this could make no sense and you may need to decide what logical unit containing b' to blast against A

[clarissaboschi]

And how can I do this comparison?

I mean I already have the output results from blast (axb and bxa). So how can I compare them and get the matches between the 2 files?

Should I use unix commands or python scripts?

[tomc]

Q. Should I use unix commands or python scripts?

A. Yes.


A less misunderstandable answer is; it all depends on your current situation.

A ready made script in any language with extra bells and whistles is nice if it does what you want and runs in your environment; but they often do not age gracefully.

Writing your own script in any language is a great skill to have and may well be critical at some point so practice is always good. But learning takes time, mistakes will be made and an entire universe of details not directly related to the biological problem at hand enter consideration; time / space trade off algorithmic complexity yada ya.

In my opinion, by default, the shell is the option to beat. In other words it pays to demonstrate that assembling a shell script is insufficient before moving on to writing your own solution in a general purpose language.

In your case you say you have blast results from A-B and B-A.
I have no idea what that means.
Blast can be output in many different formats and any sort of filtering.

But lets say you do coerce each of the two datasets you have into
a list of pairs of sequence identifiers where the second in each pair
is the best hit of the first. i.e.


<File_A_B>
seqida1 seqidbx
seqida2 seqidby
seqida3 seqidbz
...

<File_B_A>
seqidbx seqida3
seqidby seqida2
seqidbz seqida1
...


the task of finding the reciprocal best hits boils down
to finding a pair of identifiers in the first file (a row)
exactly matching a reversed pair of identifiers in the second file.

e.g.
seqida2 seqidby == reverse(seqidby seqida2)


Again I have no idea what size your lists are but assuming they fit comfortably in your machine's memory. one approach may be to:

1) sort the first file because when a list is sorted you can be certain a match could not be behind where you are in the list

sort -u File_B_A > File_A_B_sorted

2) reverse the second list

awk '{print $2,$1}' File_B_A > File_B_A_rev

3) sort the second file

sort -u File_B_A_rev > File_B_A_rev_sorted

4) see what the files have in common

comm -12 File_A_B_sorted File_B_A_rev_sorted > File_AB_RBH

The list produces would be the set of RBHs according to your blast.

Hopefully that gives you a rough outline

edit: add in the `-u` (unique) on sort to remove duplicate from each list

[clarissaboschi]

Thanks tomc for your detailed explanation and time to answer it.

This was exactly what I did. I wrote a shell script in different steps. I was not able to write a python script because I am learning...
I found a few python scripts to do it but did not work for me.

But my shell script is working well. In the last step I got the duplicates of my combined file instead of the common ones.

[maubp]

Here's a BLAST RBH tool I wrote earlier in Python, which does consider duplicates and gives warnings about them:

https://github.com/peterjc/galaxy_blast/blob/master/tools/blast_rbh/blast_rbh.py

It has a Galaxy wrapper but you can ignore that, other than perhaps reading the help text included it it - which suggests thinking about setting minimum identity and minimum alignment lengths and reading this paper:

Punta and Ofran (2008) The Rough Guide to In Silico Function Prediction,
or How To Use Sequence and Structure Information To Predict Protein
Function. PLoS Comput Biol 4(10): e1000160.

http://dx.doi.org/10.1371/journal.pcbi.1000160