含grep和sed的Shell脚本提取配对低F_MISS个体出错排查

需求与问题现状 #

需求是从Relatedness_3rdDegree.txt读取个体配对，对比每个配对中两个个体在filename.imiss里的F_MISS值，保留数值更低的个体。但原脚本执行到第91行时出现(standard_in) 1: syntax error，且输出结果不符合预期。

相关文件示例 #

• Relatedness_3rdDegree.txt：

Individual1 Individual2
Individual5 Individual23
Individual50 Individual65

• filename.imiss：

INDV    N_DATA  N_GENOTYPES_FILTERED    N_MISS  F_MISS
Individual1 375029  0   782 0.00208517
Individual2 375029  0   341 0.000909263
Individual3 375029  0   341 0.000909263

原脚本 #

numlines=$(wc -l Relatedness_3rdDegree.txt|awk '{print $1}')

for line in `seq 1 $numlines`
do
ind1=$(sed -n "${line}p" Relatedness_3rdDegree.txt|awk '{print $1}')
ind2=$(sed -n "${line}p" Relatedness_3rdDegree.txt|awk '{print $2}')
miss1=$(grep $ind1 filename.imiss|awk '{print $5}')
miss2=$(grep $ind2 filename.imiss|awk '{print $5}')
if echo "$miss1 > $miss2" | bc -l | grep -q 1
then
echo $ind1 >> miss.txt
else
echo $ind2 >> miss.txt
fi
echo "$line / $numlines"
done

问题原因分析 #

1. grep部分匹配导致取值错误：直接用grep $ind1会匹配所有包含该字符串的行，比如Individual1会误匹配Individual10，导致拿到错误的F_MISS值，甚至多行结果，传给bc后触发语法错误。
2. 逻辑完全颠倒：原脚本中当miss1 > miss2时输出ind1，但我们需要保留F_MISS更低的个体，这时候应该输出ind2，逻辑搞反直接导致结果不符合预期。
3. 未处理个体不存在的情况：如果某个个体在filename.imiss中找不到，miss1或miss2会是空值，传给bc后就会报语法错误，这就是第91行报错的核心原因。
4. 循环方式低效且易出问题：用seq+sed逐行读取的方式不仅效率低，还容易在文件包含特殊字符时出错。

修复后的脚本 #

推荐用awk一次性处理，效率更高且避免上述所有问题：

awk '
# 第一次读取filename.imiss，存储个体和对应的F_MISS值
NR == FNR {
    if (NR > 1) {  # 跳过表头
        miss[$1] = $5
    }
    next
}
# 读取Relatedness_3rdDegree.txt的配对行
{
    ind1 = $1
    ind2 = $2
    # 检查个体是否存在
    if (!(ind1 in miss) || !(ind2 in miss)) {
        print "Warning: " ind1 " or " ind2 " not found in filename.imiss" > "/dev/stderr"
        next
    }
    # 对比F_MISS值，输出更低的那个
    if (miss[ind1] < miss[ind2]) {
        print ind1
    } else {
        print ind2
    }
}
' filename.imiss Relatedness_3rdDegree.txt > miss.txt

如果要修复原Shell脚本，可修改如下：

# 直接循环读取文件行，避免seq+sed的问题
while read -r ind1 ind2
do
    # 用grep -w匹配完整单词，-w确保只匹配整行开头的个体名，避免部分匹配
    miss1=$(grep -w "^$ind1" filename.imiss | awk '{print $5}')
    miss2=$(grep -w "^$ind2" filename.imiss | awk '{print $5}')
    
    # 检查是否获取到有效值
    if [[ -z $miss1 || -z $miss2 ]]; then
        echo "Line error: $ind1 or $ind2 not found" >&2
        continue
    fi
    
    # 修正逻辑：输出F_MISS更低的个体
    if (( $(echo "$miss1 < $miss2" | bc -l) )); then
        echo "$ind1" >> miss.txt
    else
        echo "$ind2" >> miss.txt
    fi
done < Relatedness_3rdDegree.txt

内容的提问来源于stack exchange，提问作者Gf.Ena