使用awk批量统计Fasta文件reads数量的技术求助

问题描述 #

需要批量处理大量FASTA文件，统计每个文件中以> NAME开头的reads数量，输出文件需包含文件名和对应的reads数（所有文件至少含1条reads）。

FASTA文件示例 #

File1.fasta

>sequence A
ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatattcatat
tctgttgccagaaaaaacacttttaggctatattagagccatcttctttgaagcgttgtc
>sequence B
ggtaagtgctctagtacaaacacccccaatattgtgatataattaaaattatattcatat
tctgttgccagattttacacttttaggctatattagagccatcttctttgaagcgttgtc
tatgcatcgatcgacgactg

File2.fasta

>sequence A
ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatattcatat
tctgttgccagaaaaaacacttttaggctatattagagccatcttctttgaagcgttgtc
>sequence B
ggtaagtgctctagtacaaacacccccaatattgtgatataattaaaattatattcatat
tctgttgccagattttacacttttaggctatattagagccatcttctttgaagcgttgtc
tatgcatcgatcgacgactg
>sequence C
ggtaagtgctctagtacaaacacccccaatattgtgatataattaaaattatattcatat
tctgttgccagattttacacttttaggctatattagagccatcttctttgaagcgttgtc
tatgcatcgatcgacgactg

无效脚本分析 #

脚本1 #

#!/bin/bash
for file1 in ~/Test/*.fasta
do
outfile1=${file1}_readcount.txt
awk -F ' ' -v out1=$outfile1 '{
if(NR==1) {lines[0]=lines[0] OFS FILENAME  > out1;
}
if(FNR==NR) {grep -c "^>" $file1 > out1;
}
}'  $file1
done

问题：无输出的核心原因是awk脚本内不能直接调用shell的grep命令，同时NR/FNR的判断逻辑混乱，没有实际统计作用。

脚本2 #

BEGIN   { OFS="	" } #output field delimiter is tab

FNR==1  { lines[0]=lines[0] OFS FILENAME } #append filename to header record

FNR==NR {grep -c "^>" FILENAME } # counts number of ">" at the beginning of lines

END     { for (i=0;i<=FNR;i++) #loop through the line numbers
              print lines[i]
        } #printing each line 
' *fasta > countreads.txt

问题：仅生成表头是因为awk内调用grep无效，且FNR==NR仅在处理第一个文件时生效，后续文件无统计逻辑，循环输出空数组元素导致大量空行。

正确解决方案 #

方案1：纯awk单命令（简洁高效） #

直接遍历所有FASTA文件，统计每个文件中>开头的行数，输出文件名与统计数：

awk '/^>/ {count++} ENDFILE {print FILENAME, count; count=0}' *.fasta > countreads.txt

如果需要用制表符分隔输出，添加BEGIN块设置分隔符：

BEGIN {OFS="\t"} /^>/ {count++} ENDFILE {print FILENAME, count; count=0} *.fasta > countreads.txt

逻辑说明：

• /^>/ {count++}：每匹配到以>开头的行，计数器加1
• ENDFILE：处理完单个文件后触发，打印当前文件名和计数，重置计数器

方案2：bash循环+awk（适合复杂路径场景） #

若FASTA文件分布在不同路径，或需要自定义文件名格式，用bash循环逐个处理：

#!/bin/bash
output="countreads.txt"
> "$output" # 清空已有输出文件
for file in ~/Test/*.fasta; do
    read_count=$(awk '/^>/ {c++} END {print c}' "$file")
    # 仅保留文件名（不含路径），需要完整路径则直接用"$file"
    filename=$(basename "$file")
    echo -e "$filename\t$read_count" >> "$output"
done

预期输出 #

File1.fasta	2
File2.fasta	3

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