SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
BLAST+ creating custom blast database and using blast+ filtering features deniz Bioinformatics 3 07-07-2019 08:04 AM
BLAST database error - when changing to new BLAST+ local program biobio Bioinformatics 4 06-15-2011 05:20 AM
What is reciprocal best hit? Dinesh Bioinformatics 1 02-08-2011 03:51 AM
Reciprocal Blast With Multiple Contigs Per Gene jwhittall Bioinformatics 0 07-21-2010 01:20 PM
Need Help Regarding Reciprocal best hit in BLAST Bharat Bioinformatics 1 03-11-2010 03:04 AM

Reply
 
Thread Tools
Old 04-26-2012, 02:39 PM   #1
renesh
Junior Member
 
Location: Louisiana

Join Date: Sep 2011
Posts: 9
Default reciprocal blast

How to do reciprocal blast for two datasets?
renesh is offline   Reply With Quote
Old 04-26-2012, 10:27 PM   #2
tomc
Member
 
Location: Oregon

Join Date: Feb 2011
Posts: 29
Default

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.

Last edited by tomc; 04-27-2012 at 02:23 PM. Reason: typo
tomc is offline   Reply With Quote
Old 01-02-2013, 06:49 PM   #3
kenietz
Member
 
Location: Singapore

Join Date: Nov 2011
Posts: 85
Default

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
kenietz is offline   Reply With Quote
Old 11-06-2014, 01:04 PM   #4
tomc
Member
 
Location: Oregon

Join Date: Feb 2011
Posts: 29
Default

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
tomc is offline   Reply With Quote
Old 02-08-2017, 10:57 AM   #5
clarissaboschi
Member
 
Location: US

Join Date: Apr 2010
Posts: 63
Default

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?
clarissaboschi is offline   Reply With Quote
Old 02-09-2017, 12:16 PM   #6
tomc
Member
 
Location: Oregon

Join Date: Feb 2011
Posts: 29
Default

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

Last edited by tomc; 02-09-2017 at 12:56 PM. Reason: splain stuff
tomc is offline   Reply With Quote
Old 02-09-2017, 12:47 PM   #7
clarissaboschi
Member
 
Location: US

Join Date: Apr 2010
Posts: 63
Default

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.
clarissaboschi is offline   Reply With Quote
Old 02-10-2017, 08:33 AM   #8
maubp
Peter (Biopython etc)
 
Location: Dundee, Scotland, UK

Join Date: Jul 2009
Posts: 1,543
Default

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_bl...h/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
maubp is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 06:38 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO