计数 fasta 文件中的读取次数（循环）

Question

我想使用 awk 计算 fasta 文件的读取次数。 fasta 文件中的读取以“> NAME”开头，后跟一行和“DNA 代码”。就我而言，fasta 文件至少包含 1 次读取（因此不为空）。我还有很多 fasta 文件，我想循环播放它。因此，我还想将文件名复制到输出中，并将读取次数粘贴到输出文件中的文件名旁边。

Fasta 文件示例（文件 1.fasta）

>sequence A
ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatattcatat
tctgttgccagaaaaaacacttttaggctatattagagccatcttctttgaagcgttgtc
>sequence B
ggtaagtgctctagtacaaacacccccaatattgtgatataattaaaattatattcatat
tctgttgccagattttacacttttaggctatattagagccatcttctttgaagcgttgtc
tatgcatcgatcgacgactg

Fasta 文件示例（文件 2.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

它没有给出错误，但也没有输出

脚本2

awk '
BEGIN   { OFS="\t" } #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

这里我只得到文件中的标题和数千个空行。

我想得到的预期输出

File1   2
File2   3

Answer 1

如果您打算将每个

/^>/

计数为一个序列，然后按文件名汇总序列计数，您可以在 awk 中执行此操作：

awk 'FNR==1{if (name) print name, cnt; name=FILENAME; cnt=0}
/^>/{cnt++}
END{print name, cnt}' *.fasta

这也有效，但文件名的顺序可能与 awk 读取它们的顺序不同：

awk '/^>/{file_cnt[FILENAME]++}
END{for (fn in file_cnt) print fn, file_cnt[fn]}' *.fasta

您也可以直接使用

grep

来统计匹配次数：

grep -c "^>" *.fasta

总之，如果您想自定义计数或打印的内容，请使用

awk

。如果您只想计算正则表达式匹配的总数并按文件名进行汇总，请使用

grep

和 glob。

Answer 2

使用 GNU 实用程序，您可以执行以下操作：

grep -cZ '^>' ~/Test/*.fasta | tr '\0' '\t'

• -Z

  在文件名后面打印空值而不是冒号

• tr

  将其转换为选项卡

如果您知道文件名不包含冒号，则这应该适用于任何版本：

grep -c '^>' ~/Test/*.fasta | tr : '\t'