如何使用Bash脚本计算FASTA文件的GC-content？

GC-content的计算公式为：(G数量 + C数量) / (A数量 + T数量 + G数量 + C数量)。我之前尝试用wc命令直接计算，但没得到正确结果，后来查阅资料整理出了可行的Bash脚本方案：

filename=$@                     # 接收所有传入的文件名作为参数
for f in $filename              # 循环处理每个文件
do
    echo " $f 正在处理..."
    # 统计G和C的总数：过滤序列头，匹配所有大小写的G/C并逐行输出，统计行数
    gc=$( grep -v ">" "$f" | grep -io 'g\|c' | wc -l)
    # 统计有效碱基总数：过滤序列头，移除空白和换行符后统计字符数
    total=$( grep -v ">" "$f" | tr -d '\s\r\n' | wc -c)
    # 计算GC百分比，保留两位小数
    percent=$( echo "scale=2;100*$gc/$total" | bc -l)
    echo " $f 的GC含量为: "$percent"%"
    echo
done

关键步骤解释： #

• grep -v ">" "$f"：FASTA文件中以>开头的是序列注释行，用-v参数过滤掉这些行，只保留碱基序列内容
• grep -io 'g\|c'：-i忽略大小写，确保同时匹配大写的G/C和小写的g/c；-o让每个匹配的碱基单独占一行，这样wc -l就能直接统计G和C的总数量
• tr -d '\s\r\n'：移除序列中的空格、制表符、换行符，避免这些空白字符干扰碱基总数的统计
• bc -l：Bash本身不支持浮点运算，用bc工具来计算百分比，scale=2指定结果保留两位小数

我正在学习生物信息学。

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