基于BLAST结果筛选FASTA文件：保留对应最常见基因的物种序列

步骤1：获取对应最常见基因的物种列表 #

你之前的问题出在shell变量传递到awk的方式错误——单引号包裹的awk代码里，$var1不会被shell解析。下面是两种可靠的解决方式：

方式一：分步实现（先拿基因名，再筛物种） #

先通过你的脚本获取最常见基因名并赋值给shell变量：

var1=$(nawk -f don.awk blastfile)

再用-v参数将变量传递给awk，用精确匹配替代正则（避免基因名含特殊字符时出错）：

awk -F'\t' -v target="$var1" '$2 == target {print $1}' blastfile > species_list.txt

执行后species_list.txt里就存着所有需要保留的物种名称。

方式二：单awk脚本直接输出目标物种 #

无需单独生成基因名，一次遍历统计+二次遍历筛选：

awk -F'\t' '
# 第一次遍历：统计基因出现次数，记录物种-基因映射
NR == FNR {
    count[$2]++
    if (count[$2] > max_count) {
        max_count = count[$2]
        most_common = $2
    }
    sp_gene[$1] = $2
    next
}
# 第二次遍历：输出对应最常见基因的物种
FNR == NR {
    if (sp_gene[$1] == most_common) print $1
}' blastfile blastfile > species_list.txt

步骤2：用物种列表过滤FASTA文件 #

假设你的FASTA文件标题行以>开头，且后面直接跟物种名，用以下awk脚本过滤最可靠（支持多行序列）：

awk '
BEGIN {
    # 读取物种列表到数组
    while ((getline < "species_list.txt") > 0) {
        keep[$1] = 1
    }
}
# 处理标题行，标记是否保留当前序列
/^>/ {
    current_sp = substr($0, 2) # 去掉开头的>
    is_keep = (current_sp in keep)
}
# 标记为保留的序列，输出所有行
is_keep {print}' input.fasta > filtered.fasta

如果物种名无特殊字符，也可以用更简洁的grep组合（仅适合单行序列的FASTA）：

grep -Ff <(sed 's/^/>/' species_list.txt) -A 1 input.fasta > filtered.fasta

一步完成（无中间文件） #

把所有逻辑合并到一个awk脚本，直接从blast文件到过滤后的FASTA：

awk -F'\t' '
# 处理blast文件：统计最常见基因，保存物种-基因映射
ARGIND == 1 {
    count[$2]++
    if (count[$2] > max_count) {
        max_count = count[$2]
        most_common = $2
    }
    sp_map[$1] = $2
    next
}
# 处理FASTA文件：筛选目标物种的序列
ARGIND == 2 {
    if (/^>/) {
        sp = substr($0, 2)
        keep = (sp_map[sp] == most_common)
    }
    if (keep) print
}' blastfile input.fasta > filtered.fasta

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