研究メモ〜バイオインフォっぽいこと〜

こんにちは、nkjmyです。

平成最後の記事を書いたきり、新しい環境で研究していたら、

すっかり令和初めての記事が8月に・・・。

徐々に新天地での研究も進み、

色んな微生物を対象に実験をできています。

培養楽しい。

 さて、8月末から学会2連ちゃんで（国内学会は発表せんけど）、

その準備やら論文でバタバタしているわけですが、

先日、twitterでアドバイスをもらった研究関連の小技？をメモがてら残しておきます。

私のようにバイオインフォっぽいことをやりたくても、

ウェットごりごりで、やれることと言ったらfastaファイルをオンライン上のBLASTに投げるくらいです！

みたいなパターンってあると思うんですよね。

ようやくMacのターミナルでなんか色々いじれるようになってきた

レベルの小技ですが、私と同じような悩みを解決できるようにと

やることは

Local BLASTで大量のfastaファイルをデータベースとして検索する

です

(Local BLASTそのものの導入はここで解説するまでもなく、色んなわかりやすいサイトや統合TVさんがありますので、そちらを参照に。https://togotv.dbcls.jp/20170606.html)

通常、オンラインでは、

type materialの指定や、uncultureなものの除去、

特定の分類群内での検索あるいは

fasta vs fasta　のような2sequencesの比較が可能です

でも、

自分の(集めてきた)ゲノムA, B, C, D, ・・・に対して、

何か塩基配列でもアミノ酸配列でもmulti fastaファイルをqueryとして投げて検索したい！

みたいなことがあり得ますよね。

subject側のゲノムが少なければ、

-makeblastdbでそれぞれのファイルをデータベース化して、

-db "db1 db2"みたいな感じで複数個指定するなり、

-subject file.fna にすると直接データベース化なしで検索するなり

で解決可能なのかなと。

では、ゲノムA, B, C, ,,,,,AA, AB, AC,,,,,BB, BC,,,,,と50個、100個、1000個あったらどうか。

今回私の案件では、GTDBのrepresentative genomeを約24000個持ってきて(そこから一部を使う)というものだったので、

-makebalstdbにしても-subjectにしてもを数千回レベルでやるのはさすがに面倒・・・

じゃあ、

-makeblastdb file1.fna file2.fna

-subject file1.fna file2.fna  みたいに複数指定できるかというと、エラーが返ってくる

ツイッターできいてみたら、やはり、そのままでは複数指定することはできないとのことだったが、

「catで元のゲノム配列を繋ぐのがいいのでは」というアドバイスをもらった。

cat

についてはターミナルのコマンドの基本の1つ（ということも学んだ）

cat fileA fileB > out_file

みたいな感じで、

テキストファイルを「開いてコピペして別のファイルにペースト」みたいなことをしなくていい

それをさらに発展させた形で別の方にも

「-subjectでもできますよ」と教わった。

コマンドは以下

cat file1 file2 file3 file4 file5 | blastn -query query.fasta -subject /dev/stdin (追加のオプションがあればこの後に)

という感じ

ミソは

cat　〜〜〜  |　←これパイプで繋いで、前の処理をこの後のコマンドに引き渡すというものらしい

つまり、

cat ~~でたくさんファイルを繋いで、それを次の処理(blastnなりblastpなりtblastnなり)に渡して、

-queryは何か投げかけたい配列なので、まぁそのままの使い方

そして次のミソが、

-subject /dev/stdin

詳しいことは自分自身ちゃんと理解できていませんが、「標準入力」というものを指定していて、そこからデータを受け取っていると。

つまり、catで繋いだものが

（仮想的な状態で）/dev/stdin　に存在していて、

それを-subjectの指定先にしているので、複数指定には該当せず、うまく走る

という形です。（一応うまくいってるけど、どこまで正確な理解かは・・・うん・・・）

これで、cat でつなぐ時は、フォルダのファイルをドラッグなり何なりでたくさん洗濯しておいて、

それをズズズッとターミナルに引っ張っていけば大量のファイル名が入力されていくので問題ない。

|　

これを忘れず。

色々いじって見た結果、

2600個くらいまでのゲノムなら繋いで渡せる。

（サイズが色々ばらけてるせいかな）

それ以上になると、「多すぎる」ってターミナルさんに怒られる

2500個のゲノムで

cat ~~~ | tblastn -query fileA.faa -subject /dev/stdin -outfmt 6 -evalue 1E-10 -out result_file

みたいな感じでやったら(queryはアミノ酸配列20個くらい)、

1回の解析で15分とか20分くらいやったかな？

ちなみに、メモリがごっそり持っていかれる

32GBの搭載で通常、20〜22GBくらいが空いている状態でしたが、

上のを走らせると、時折、2GBくらいまで埋まった。

ちなみに、結果を色々見ている中で他に使ったコマンドは、

uniq aaaa1 aaaa2

あいうえお

あいうえお

かきくけこ

かきくけこ

さしすせそ

↓

あいうえお

かきくけこ

さしすせそ

みたいに同じ文字列の行を消して(aaaa1)、新しいテキストファイルを作れる(aaaa2)

あるいは、

grep -c "あいう" aaaa1

で、aaaa1内にある「あいう」の数を数えられる

（このままだと、あいうえお　の前半部分を認識して2と出るはずだが）

これでとある遺伝子のhitの数とかも数えられる

（かつて発現解析いじってた時にも使ってた）

少しレベルが上がった気がしました、まる。

終わり