使用 awk 计算每个物种 id 在多个 fasta 文件中出现的次数答案

【问题标题】：Using awk to count how many time each species id occurs in multi fasta files使用 awk 计算每个物种 id 在多个 fasta 文件中出现的次数
【发布时间】：2014-08-07 03:22:20
【问题描述】：

我搜索了这个主题，但没有找到。我有 5593 个多 fasta 文件，我需要计算每个物种 id 在每个文件中出现的次数。我只能识别每个物种的序列总数，但我无法识别输入文件。

输入

file1.fasta：

>hsa
ATCGATCGATCAGACTACG

>eco
ATCGATCGATCAGACTACG

file2.fasta：

>hsa
GATCGATCAGACTACGAAA

>hsa
GATCGATCACAGACTACGAAA

file3.fasta：

>hsa
CTAGACTAGATAGACACATAGAGA

>ecj
CTAGACTAGCTAGACCCATAGAGA

>mmu
CTAGACAAGATAGACACAAAGAGA

>eco
CTAGACTACATCGACACATAGAGA

预期输出

file1.fasta >hsa [count]
file1.fasta >eco [count]
file2.fasta >hsa [count]
file3.fasta >hsa [count]

file3.fasta >ecj [count]
file3.fasta >mmu [count]
file3.fasta >eco [count]

awk /^>.../ {print $1} *.* | sort | uniq -c | sort -nr

输出

[total counts]>hsa

[total counts]>eco

[total counts]>mmu

[total counts]>ecj

【问题讨论】：

标签： unix awk count fasta

【解决方案1】：

假设“物种”行以> 开头，并且示例输出中的方括号表达式是简单数字，则：

awk 'BEGIN { SUBSEP = " " }
     /^>/ { per_file[FILENAME,$1]++; total[$1]++ }
     END { for (k in per_file) print k, per_file[k]
           for (k in total)    print total[k], k
         }' *.fasta

您可能需要在awk 或之后的某处进行排序，因为无法保证for (index in array) 循环呈现的数据将按任何特定顺序排列。

如果没有BEGIN 块（或其他机制）设置SUBSEP，文件名之后和物种键之前会有一个\034 字符。通过将SUBSEP 设置为空白，文件名与物种键之间用空白分隔。

【讨论】：

【解决方案2】：

这是使用awk的一种方式：

awk '
BEGIN { SUBSEP = FS }
!(FILENAME in name) { fileorder[++num]=FILENAME; name[FILENAME]++ }
/^>/ { count[FILENAME,$1]++ }
END { 
  for (order=1; order<=num; order++) {
      for (total in count) {
          split(total, keys, SUBSEP)
          if (fileorder[order]==keys[1]) {
              print fileorder[order], keys[2], count[fileorder[order],keys[2]]
          }
      }
  }
}' file1 file2 file3

输出：

file1 >eco 1
file1 >hsa 1
file2 >hsa 2
file3 >eco 1
file3 >mmu 1
file3 >hsa 1
file3 >ecj 1

在BEGIN 块中，我们将SUBSEP 设置为FS。这允许复合键用空格分隔。
我们维护两个数组，name 和 fileorder。如果我们的文件不在name 数组中，我们将使用文件名填充fileorder 数组。
对于以 fasta 模式开头的行，我们填充另一个名为 count 的数组，其中文件名和模式作为复合键并计为值。
在END 块中，我们遍历fileorder 数组，然后是count 数组。我们使用split 函数拆分复合键，并检查文件名是否与我们的键的第一部分匹配。如果是，我们打印文件名、模式和计数，然后继续下一次迭代。

【讨论】：