2種の遺伝子配列同士の類似配列検索

2種の生物の遺伝子やアミノ酸配列があったとして，両者間に似た遺伝子が存在するかどうかを検索する方法．

計算用サーバーに移動し，

tmux a

でいずれかのノードに入る．tmuxウィンドウの変更はctrl+b n，tmuxを抜けるときはctrl+b d．

fastaが置いてあるディレクトリに移動し，2種どちらかのfastaファイルをサーバー上でプライベートデータベース化する．

makeblastdb -in 目的のfastaファイル名 -dbtype nucl（塩基配列）か prot（アミノ酸配列） -out ******-database（作成DB名）

比較対象の配列をNCBIからとってくる場合は，

wget https://ftp.ncbi.nlm.nih.gov/genomes/GCF_014107475.1_ASM1410747v1_protein.faa（目的配列の圧縮ファイルURL）

ゲノムが決定されている生物は，個別ページの一番上に四角が表示される．この枠内にあるDownload sequences in FASTA format for genome, protein のプロテインからアミノ酸配列DLのURLを抽出（ゲノムではないのはblast検索のため）．geneがある場合は，これを選択．

各クエリ配列が作成したデータベース内のリードに何度ヒットしたかという情報を，全クエリ配列と共に一覧化する．

tblastn -db ******-database（作成したDB名） -query /home/*/*/*/GCF_014107475.1_ASM1410747v1_protein.faa（クエリ配列が違うディレクトリの場合はパスを指定） -outfmt 7 -evalue 1e-15 -out /home/*/*/*/*/*****_vs_*******_tblastn_1e-15.txt（アウトプットするファイル名，必ずblast形式とe-valueを併載）

e-valueは，塩基同士なら1e-10，塩基対アミノ酸だと1e-15など．とりあえず1e-5とかで緩めにして，後からフィルタリングしてもよい．

-outfmt オプションの7はタブ区切りのみのシンプルなフォーマット (tabular with comment lines)なので，プログラムでの処理がしやすくなる．一方，一対一の配列比較などでアラインメントをしっかり見たいときはデフォルトの0のペアワイズで見るとよい．

クエリ配列の中に遺伝子が何個入っているかを数える．

grep ">" /home/*/*/*/*/GCF_014107475.1_ASM1410747v1_protein.faa | wc -l

grep コマンドによってダブルクオーテーションに囲まれた行を書き出し（-v オプションを使うと，反対に該当しない行を書き出せる），wc -lで行数を数え上げることができる．

fastaなので行頭の>の数を数えれば，それが遺伝子の数と一致する．今回は1003遺伝子と算出．

次に，作成した2種の遺伝子配列同士の一致を調べた一覧の中で，DBの遺伝子がヒットしなかったクエリ配列の数を検索する．

grep "# 0 hits found" /home/*/*/*/*/*****_vs_Wolbachia-pipientis_tblastn_1e-15.txt | wc -l

DBにヒットしなかったクエリ配列は# 0 hits foundという文字列を含むので，それを数えるとよい。今回は993遺伝子がDBにヒットしなかったものとして算出．

1003-993=10遺伝子が2種間で共通する遺伝子として1e-15のe-valueでヒットしたことになる．