双方向ベストヒットを取る
何エントリか前に紹介した遺伝子対応表を作る際に、一般的には双方向ベストヒットということをよくやります。これは2種類の生物種にコードされた遺伝子セット全体を比較する際、クエリとした生物種1とデータベースとした生物種2でBLASTをかけた後にクエリに生物種2、データベースに生物種1を割り当てて逆に計算し、その対応関係において双方向で相思相愛になるペアだけを選ぶというものです。詳しくは坊農の博士論文[PDF] 2.3オーソロガス遺伝子の推定(p.12-14)を見てください。それをやる具体的な手順を今回紹介します。
BLASTをかけて-outfmt 6で得られた結果に対して、cutコマンドを使って1カラム目と2カラム目だけを週出してuniqコマンドで重複を除きます。 [shell] cut -f1,2 kaiko_human.txt |uniq > kaiko-human1.txt [/shell] また逆のBLAST検索の結果も同様に処理します。 [shell] cut -f2,1 human_kaiko.txt |uniq > kaiko-human2.txt [/shell] で、得られた結果の2つのファイルの同じ部分が双方向ベストヒットのペアということになります。 それを抽出するにはちょっとトリッキーですが、以下のようにしています。 [shell] cat kaiko-human1.txt kaiko-human2.txt | sort > sort_only.txt cat kaiko-human1.txt kaiko-human2.txt | sort | uniq > sort_uniq.txt diff sort_only.txt sort_uniq.txt [/shell] 2つのファイルを結合(cat)してソート(sort)しただけのものと、さらに重複を除いた(uniq)ものの差分(diff)をとっています。両方のファイルに含まれていた(=双方向ベストヒット)のものは2回出てくるわけですから、これで差分となって現れるわけです。 diffの出力から双方向ベストヒットのペアだけを抽出したいのであれば以下の様なコマンドで処理するといいでしょう。 [shell] diff sort_only.txt sort_uniq.txt | grep ^< | awk '{ print $2 "t" $3 }' [/shell] 最後のところはもちろんcutコマンドでもできるのですが、今回はawkを使っております。