蛋白配列や核酸配列の冗長性を除きたいときには、blastclustを使うとよいらしい。

blastclust -i infile -o outfile -p F -L .9 -b T -S 95

上の例だと、infileの配列を、配列の90%以上の領域で95%以上の同一性を持つものをクラスタリングして出力します。

でも、実際に検索したり、配列いじくったりする側からすると、最初に配列リストありきみたいな状況は稀で、大概次のような処理フローになると思うから、コレだけだとイマイチ使えん。

1. NR,NT,SNP,ESTなんかにblast検索
2. ヒットしたIDから配列リストをゲット
3. マルチプルアライメント or クラスタリング

blastallでfastalistを吐くか、id-listを吐けば　blastall | fastacmd が出来るから、

blastall -p blastn -d nt -i XXX.seq -e 1e-20 | blastclust -o outfile -p F -L .9 -b T -S 95

こんな感じでワンライナーになって使いやすいかな。いまのとこ、同じことやる自前のスクリプトのほうが取り回しがらくだ。

ん？

1をタブ区切りで出力させて、fastacmdで2を実行するコマンド作っておけばワンライナーへの道は開けるか？

やっとPubmedがRSS配信するようになったのね。

というわけで、PubmedのRSSをbloglinesかfeed on feedsでチェックすればよいんだけど、既読と未読の管理をしないと効率よく捌けないので、なにか作らないといけないな。

Catalyst使うときが近づいているのか？

なんて。

やっぱ、「暗黙知の共有」が一番おもろいよね～と思うことしばしば。

ってことで、ごく限られたチームの中で hubmedとか使えないかナァと思ったり。 これとblosxomを組み合わせて適当にプラグイン入れれば何かがokな感じもするので、やってみようかなぁと思ったりするわけだヨ。

blastからfastaを出力するプログラム欲しいってことで自分で書いてみた。

最初はGetoptモジュール使ってちゃんとやろうと思ったんだが、blastallってオプション多すぎていちいち面倒なんで、ARGV適当につなげて、オプションのチェックはblastallにまかせた。

普通にblastallのオプションが使えて

# bl2fasta -p blastp -i top2.seq -d pdbaa >gi|4929914|pdb|1BJT| Topoisomerase Ii Residues 409 - 1201 ENALKKSDGTRKSRITNYPKLEDANKAGTKEGYKCTLVLTEGDSALSLAVAGLAVVGRDYYGCYPLRGKMLNVREASADQ ILKNAEIQAIKKIMGLQHRKKYEDTKSLRYGHLMIMTDQDHDGSHIKGLIINFLESSFLGLLDIQGFLLEFITPIIKVSI TKPTKNTIAFYNMPDYEKWREEESHKFTWKQKYYKGLGTSLAQEVREYFSNLDRHLKIFHSLQGNDKDYIDLAFSKKKAD DRKEWLRQYEPGTVLDPTLKEIPISDFINKELILFSLADNIRSIPNVLDGFKPGQRKVLYGCFKKNLKSELKVAQLAPYV SECTAYHHGEQSLAQTIIGLAQNFVGSNNIYLLLPNGAFGTRATGGKDAAAARYIYTELNKLTRKIFHPADDPLYKYIQE DEKTVEPEWYLPILPMILVNGAEGIGTGWSTYIPPFNPLEIIKNIRHLMNDEELEQMHPWFRGWTGTIEEIEPLRYRMYG RIEQIGDNVLEITELPARTWTSTIKEYLLLGLSGNDKIKPWIKDMEEQHDDNIKFIITLSPEEMAKTRKIGFYERFKLIS PISLMNMVAFDPHGKIKKYNSVNEILSEFYYVRLEYYQKRKDHMSERLQWEVEKYSFQVKFIKMIIEKELTVTNKPRNAI IQELENLGFPRFNKEGKPYYGSPNDEIAEQINDVKGATSDEEDEESSHEDTENVINGPEELYGTYEYLLGMRIWSLTKER YQKLLKQKQEKETELENLLKLSAKDIWNTDLKAFEVGYQEFLQRDAEARGGNVPNKGSKTKGKGKRKLVDDED >gi|3891992|pdb|1AB4| 59kda Fragment Of Gyrase A From E. Coli VGRALPDVRDGLKPVHRRVLYAMNVLGNDWNKAYKKSARVVGDVIGKYHPHGDSAVYDTIVRMAQPFSLRYMLVDGQGNF GSIDGDSAAAMRYTEIRLAKIAHELMADLEKETVDFVDNYDGTEKIPDVMPTKIPNLLVNGSSGIAVGMATNIPPHNLTE VINGCLAYIDDEDISIEGLMEHIPGPDFPTAAIINGRRGIEEAYRTGRGKVYIRARAEVEVETIIVHEIPYQVNKARLIE KIAELVKEKRVEGISALRDESDKDGMRIVIEGEVVLNNLYSQTQLQVSFGINMVALHHGQPKIMNLKDIIAAFVRHRREV VTRRTIFELRKARDRAHILEALAVALANIDPIIELIRHAPTPAEAKTALVANPWQLGNVAAMLEDAARPEWLEPEFGVRD GLYYLTEQQAQAILDLRLQKLTGLEHEKLLDEYKELLDQIAELLRILGSADRLMEVIREELELVREQFGDKRRTEIT >gi|66361362|pdb|1Z84|B Chain B, X-Ray Structure Of Galt-Like Protein From Arabidopsis Thaliana At5g18200

または、パイプでつなげて

# bl2fasta -p blastp -i top2.seq -d pdbaa | blastclust May 30, 2005 8:59 PM Start clustering of 6 queries 52696130 66361362 4929914 3891992 18159024 15826535

みたいにも使える。

使い道？

• ホモロジーモデリングの際blast->clustalwをかける際に無駄に重複する配列を除く
• SNPのデータベースに当てたとき、冗長な配列を除く

とか、お手軽に違いない。まぁ、30分で書いたからバグだらけと思うが、バグったらゴメンヨ。