基于指定PFAM ID提取蛋白序列并按家族分组的脚本实现求助

问题背景 #

已有存储PFAM家族ID的文件pfamacc（每行一个PFAM编号），需要从多个.faa.out文件中筛选对应PFAM的蛋白accession号，再从fasta序列文件中提取这些蛋白的序列，最终按PFAM家族分别保存到独立文件中。原脚本将所有结果输出到同一文件，现提供正确实现方案。

实现脚本 #

#!/bin/bash

# 1. 生成PFAM ID的正则匹配模式，用于筛选.out文件
PFAM_PATTERN=$(tr '\n' '|' < pfamacc | sed 's/|$//')

# 2. 从所有.out文件中提取目标PFAM对应的accession，建立PFAM与accession的映射
grep -E "$PFAM_PATTERN" *.faa.out | awk '{print $4, $1}' > pfam_to_acc.txt

# 3. 提前创建每个PFAM对应的输出文件，避免写入异常
while read pfam_id; do
    touch "${pfam_id}.fasta"
done < pfamacc

# 4. 根据映射关系，从fasta文件中提取序列并分组保存
awk '
BEGIN {
    # 加载PFAM与accession的映射表
    while ((getline line < "pfam_to_acc.txt") > 0) {
        split(line, arr, " ")
        acc_map[arr[2]] = arr[1]
    }
}
/^>/ {
    # 提取当前蛋白的accession号
    curr_acc = substr($0, 2)
    curr_pfam = acc_map[curr_acc]
}
curr_pfam != "" {
    # 将序列写入对应PFAM的文件
    print >> curr_pfam ".fasta"
}
' your_proteins.fasta

# 可选：清理临时映射文件
rm pfam_to_acc.txt

关键步骤说明 #

• 生成匹配模式：将pfamacc中的每行PFAM ID转换为grep支持的正则表达式（如PF12312|PF43555），实现一次性筛选所有目标PFAM。
• 建立映射关系：遍历所有.faa.out文件，筛选出包含目标PFAM的行，提取PFAM ID和对应的蛋白accession，保存到临时映射文件。
• 初始化输出文件：提前创建每个PFAM对应的空fasta文件，确保后续写入时不会因文件不存在报错。
• 提取并分组序列：用awk读取fasta文件，结合映射表判断当前蛋白所属的PFAM家族，将序列追加到对应家族的文件中。

注意事项 #

• 将脚本中的your_proteins.fasta替换为实际的fasta序列文件名；若有多个fasta文件，可改为*.fasta并循环处理。
• 若.faa.out文件不在当前目录，需将*.faa.out替换为具体路径（如/data/out_files/*.faa.out）。
• 确保使用GNU awk（大部分Linux系统默认自带），避免写入文件时出现兼容性问题。

内容的提问来源于stack exchange，提问作者aleksowski