Fastaファイルからデータベースにヒットした配列のみを抽出する #BLAST #DIAMOND #mmseqs2

手持ちのfastaファイルをデータベースに照合し、ヒットした配列のみを抽出する方法です。おすすめはDIAMONDです。

Blastで検索し、抽出する場合

# protein.faaはfasta形式のタンパク質データベース
makeblastdb -in protein.faa -dbtype prot -hash_index
# 塩基配列の場合
makeblastdb -in nucleotide.fasta -dbtype nucl -hash_index


# blast検索
blastp -query A.fasta -db protein.faa -evalue 1e-5 -outfmt 7 -num_alignments 1 > blast_out.txt

# Top HItを抽出
grep -v "#" blast_out.txt > blast_out.tsv

# IDリスト作成
awk -F "\t" '{print $1}' blast_out.tsv > label

# 冗長なラベルを処理
uniq label > label_2

# label_2をもとにA.fastaから配列を抽出
# -f: 抽出したいIDリストを指定
# -o: 出力ファイル名を指定
seqkit grep -f label_2 A.fasta -o A_Blast.fasta

DIAMONDで検索し、抽出する場合

# protein.faaはfasta形式のタンパク質データベース
diamond makedb --in protein.faa -d diamond_db

# DIAMONDでblastp検索
  diamond blastp \
    --query A.fasta \
    --db diamond_db.dmnd \
    --out diamond_out.tsv \
    --outfmt 6 \
    --evalue 1e-5 \
    --threads 30 \
    --max-target-seqs 1 \
    --query-cover 50 \
    --subject-cover 50

# outputファイルからIDを抽出
awk '{print $1}' diamond_out.tsv > label

# 冗長なラベルを処理
uniq label > label_2

# label_2をもとにA.fastaから配列を抽出
# -f: 抽出したいIDリストを指定
# -o: 出力ファイル名を指定
seqkit grep -f label_2 A.fasta -o A_diamond.fasta

mmseqs2で検索し、抽出する場合

# protein.faaはfasta形式のタンパク質データベース
mmseqs createdb protein.faa mmseqs_db

# mmseqs2で検索
mmseqs easy-search A.fasta mmseqs_db mmseqs_out.txt tmp_mmseqs --threads 20 --format-mode 4 --max-seqs 1
  
# outputファイルからIDを抽出
awk '{print $1}' mmseqs_out.txt > label

# 冗長なラベルを処理
uniq label > label_2

# label_2をもとにA.fastaから配列を抽出
# -f: 抽出したいIDリストを指定
# -o: 出力ファイル名を指定
seqkit grep -f label_2 A.fasta -o A_mmseqs.fasta