Shell循环结合angsd/realSFS/AWK处理BAM文件执行逻辑咨询

这是一套二代测序群体遗传学分析流程，核心目的是批量计算BAM文件对应的样本全基因组杂合度：

• 逐行读取bams文件中存储的BAM文件路径
• 对每个BAM文件调用angsd生成位点等位基因频率似然（SAF）索引文件
• 调用realSFS基于SAF索引计算最大似然等位基因频谱
• 从频谱结果中提取总有效位点数、杂合位点数，计算杂合度，最终将所有样本的统计结果汇总到同一个结果文件

1. awk中-v file赋值的变量对应哪个文件 #

file变量的值来自$F，也就是执行basename $FF得到的带.bam后缀的原始BAM文件名（仅去掉了路径前缀），既不是去掉.bam后缀的前缀名，也不是saf.idx或ml文件。

注意：去掉.bam后缀的操作是赋值给F_PREFIX变量时单独做的，没有修改$F本身的值。

2. awk后传入的${F_PREFIX}.ml参数的作用 #

这个参数是awk的指定输入文件，awk会逐行读取这个文件的内容，对每一行执行你写的print计算逻辑。

你的初步理解存在几处偏差，且原写法存在会导致运行失败的逻辑bug：

• 管道完全失效：你写的realSFS ${F_PREFIX}.saf.idx >${F_PREFIX}.ml用>把realSFS的所有标准输出重定向到了ml文件，管道左侧没有任何内容传递给awk，|后面的awk命令根本接不到realSFS的输出
• 时序错误：重定向和管道连接的命令是并行启动的，你启动awk读ml文件的时候，realSFS可能还没开始往ml文件里写内容，大概率会读到空文件或者不完整的文件
• 逻辑认知偏差：awk不会读取BAM文件内容，它的输入只有你传入的ml文件；awk的计算结果是直接追加到goodbams.goodsites.het，不会写入ml文件，ml文件里存的始终是realSFS的原始输出。

如果你的需求是保留realSFS的原始ml输出，同时统计杂合度汇总到总文件，正确写法如下，同时修复了原for cat遍历文件容易因文件名带空格/特殊字符报错的问题：

#!/bin/bash
while IFS= read -r FF
do
    F=$(basename "$FF")
    F_PREFIX=${F/.bam/}
    # 等angsd运行完成后再跑realSFS
    angsd -i "$F" -anc "$GENOME_REF" $FILTERS -GL 1 -doSaf 1 -doCounts 1 -out "${F_PREFIX}" && \
    # 等realSFS写完ml文件后再执行awk计算
    realSFS "${F_PREFIX}.saf.idx" > "${F_PREFIX}.ml" && \
    # 如果需要file变量存不带.bam的样本名，把$F改成$F_PREFIX即可
    awk -v file="$F" '{sum=$1+$2+$3; print file"\t"sum"\t"$2/sum}' "${F_PREFIX}.ml" >> goodbams.goodsites.het
done < bams

补充说明：单样本的realSFS输出默认只有1行3列，三个值分别对应纯合参考、杂合、纯合替代的有效位点计数，你写的$1+$2+$3计算总位点数、$2/总和计算杂合度的逻辑是正确的。

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