这条计算 fasta 文件中核苷酸数量的 awk 行是如何工作的？

Question

我目前正在学习使用 awk，并找到了一个我需要的 awk 命令，但不完全理解其中发生了什么。这行代码采用一个名为 fasta 的基因组文件，并返回每个序列的所有长度它。对于那些不熟悉 fasta 文件的人来说，它们是 txt 文件，可以包含多个称为 contigs 的基因序列。它遵循以下一般结构：

>Nameofsequence
Sequencedata like: ATGCATCG
GCACGACTCGCTATATTATA
>Nameofsequence2
Sequencedata

该行在这里找到：

cat file.fa | awk '$0 ~ ">" {if (NR > 1) {print c;} c=0;printf substr($0,2,100) "\t"; } $0 !~ ">" {c+=length($0);} END { print c; }'

我知道 cat 正在打开 fasta 文件，检查它是否是序列名称行，并在某个时候计算数据部分中的字符数。但我不明白它是如何分解子字符串中的数据部分的，也不明白它是如何用每个新序列重置计数的。

---

由 Ed Morton 编辑：这是上面由 gawk -o- 清晰格式化的 awk 脚本：

$0 ~ ">" {
    if (NR > 1) {
        print c
    }
    c = 0
    printf substr($0, 2, 100) "\t"
}

$0 !~ ">" {
    c += length($0)
}

END {
    print c
}

Answer 1

先格式化命令：

awk '
  $0 ~ ">" {
    if (NR > 1) {print c;}
    c=0;
    printf substr($0,2,100) "\t";
  }
  $0 !~ ">" {
    c+=length($0);
  }
  END { print c; }
  ' file.fa

代码将使用c 进行字符计数。此计数从值0 开始，每次解析带有> 的行时将重置为0。
当输入行没有> 时，输入行的长度将添加到c。
c 的值必须在一个序列之后打印，所以当它找到一个新的>（不在第一行）或解析完整文件时（用END 块）。
正如您现在可能已经理解的那样：
breaking down the data section in substrings 是通过将输入行与 > 匹配，
resetting the counts with each new sequence 是通过在块中使用 c=0 和 $0 ~ ">" 来完成的。

看Ed的评论：printf语句用错了。我不知道 %s 在 fasta 文件中出现的频率，但这并不重要：使用 %s 输入字符串。

Answer 2

@WalterA 已经通过解释脚本的作用回答了您的问题，但如果您有兴趣，这里有一个改进版本，包括几个小错误修复，供您使用 printf input 并在输入文件时打印空行是空的，改进了对相同条件的冗余测试两次并测试 > 并单独删除它而不是一次全部删除：

BEGIN { OFS="\t" }
sub(/^>/,"") {
    if (lgth) { print name, lgth }
    name = $0
    lgth = 0
    next
}
{ lgth += length($0) }
END {
    if (lgth) { print name, lgth }
}

或者你可以这样做：

BEGIN { OFS="\t" }
sub(/^>/,"") {
    if (seq != "") { print name, length(seq) }
    name = $0
    seq = ""
    next
}
{ seq = seq $0 }
END {
    if (seq != "") { print name, length(seq) }
}

但附加到变量很慢，因此为序列的每一行调用 length() 实际上可能更有效。