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

Question

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

>Nameofsequence
Sequencedata like: ATGCATCG
GCACGACTCGCTATATTATA
>Nameofsequence2
Sequencedata

在此处找到该行：

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

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

Ed Morton 编辑：这是上面的 awk 脚本，由 gawk -o-:

清晰地格式化

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

[=12=] !~ ">" {
    c += length([=12=])
}

END {
    print c
}

Answer 1

首先格式化命令：

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

代码将使用 c 作为字符 count.This 计数从值 0 开始，并且每次解析带有 > 的行时都会重置为 0。
当输入行没有 >.
时，输入行的长度被添加到 c c 的值必须在一个序列之后打印，因此当它发现一个新的 >（不在第一行）或当整个文件被解析时（块 END）。
正如您现在可能已经了解的那样：
breaking down the data section in substrings 是通过将输入行与 > 和
匹配 resetting the counts with each new sequence 是通过在 [=22=] ~ ">".

块中使用 c=0 完成的

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

Answer 2

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

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

或者你可以这样做：

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

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

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

How does this awk line that counts the number of nucleotides in a fasta file work?

unix

awk

bioinformatics