如何用Bash统计多序列FASTA文件中每条序列的指定碱基出现次数

下面是几个无需额外依赖的实用方法，直接用基础命令就能实现需求：

方法一：用Awk（推荐，灵活易改） #

这个方法能自动识别序列名和对应序列，统计指定字符的出现次数，还能轻松切换要统计的目标字符：

脚本版（方便复用） #

创建一个名为count_char_per_seq.awk的文件，内容如下：

BEGIN{char="G"}  # 替换成你要统计的字符，比如"A"、"C"
/^>/ {
    if (seqname != "") {
        print seqname ": " count
    }
    seqname = substr($0,2)  # 提取>后面的序列名称
    count=0
    next
}
{
    count += gsub(char,"&")  # 统计当前行的目标字符数量并累加
}
END{print seqname ": " count}  # 输出最后一条序列的统计结果

运行命令：

awk -f count_char_per_seq.awk my_sequence.fasta

单行命令版（直接执行） #

如果不想保存脚本，直接用单行命令，把G换成你需要统计的字符即可：

awk 'BEGIN{char="G"} /^>/{if(seqname!="")print seqname": "count; seqname=substr($0,2); count=0; next} {count+=gsub(char,"&")} END{print seqname": "count}' my_sequence.fasta

方法二：Sed+Awk合并序列行后统计 #

先把每条序列的名称和对应序列合并成一行，再统计字符数：

sed '/^>/!{H;$!d};x;s/\n//g' my_sequence.fasta | awk -v char='G' '{n=split($0,a,">"); for(i=2;i<=n;i++){split(a[i],b," "); cnt=gsub(char,"&",b[2]); print b[1]": "cnt}}'

示例输出（针对你提供的FASTA文件，统计G的次数） #

sequence1: 5
sequence2: 6
sequence3: 4
sequence4: 5

内容的提问来源于stack exchange，提问作者Gerald Vasquez Aleman