我有以下FASTA文件：
>header1
CGCTCTCTCCATCTCTCTACCCTCTCCCTCTCTCTCGGATAGCTAGCTCTTCTTCCTCCT
TCCTCCGTTTGGATCAGACGAGAGGGTATGTAGTGGTGCACCACGAGTTGGTGAAGC
>header2
GGT
>header3
TTATGAT
我想要的输出：
>header1
117
>header2
3
>header3
7
# 3 sequences, total length 127.
这是我的代码：
awk '/^>/ {print; next; } { seqlen = length($0); print seqlen}' file.fa
我用这段代码得到的输出是：
>header1
60
57
>header2
3
>header3
7
我需要一个小修改来处理多个序列行。
我还需要一种方法来获得总序列和总长度。任何建议都将受到欢迎......请用bash或awk表示。我知道在Perl / BioPerl中很容易做到这一点，实际上我有一个脚本可以用这些方式来做。
一个 awk / gawk 解决方案可以分为三个阶段：
每次 header 发现应该执行这些操作：
- 打印上一个seqlen 如果存在。
- 打印标签。
- 初始化 seqlen。
- 为了
sequence 我们只需要这些线条 积累总数。
- 终于来了
END 舞台我们打印 残余的seqlen。
评论代码：
awk '/^>/ { # header pattern detected
if (seqlen){
# print previous seqlen if exists
print seqlen
}
# pring the tag
print
# initialize sequence
seqlen = 0
# skip further processing
next
}
# accumulate sequence length
{
seqlen += length($0)
}
# remnant seqlen if exists
END{if(seqlen){print seqlen}}' file.fa
一个 oneliner：
awk '/^>/ {if (seqlen){print seqlen}; print ;seqlen=0;next; } { seqlen += length($0)}END{print seqlen}' file.fa
对于总数：
awk '/^>/ { if (seqlen) {
print seqlen
}
print
seqtotal+=seqlen
seqlen=0
seq+=1
next
}
{
seqlen += length($0)
}
END{print seqlen
print seq" sequences, total length " seqtotal+seqlen
}' file.fa
一个 awk / gawk 解决方案可以分为三个阶段：
每次 header 发现应该执行这些操作：
- 打印上一个seqlen 如果存在。
- 打印标签。
- 初始化 seqlen。
- 为了
sequence 我们只需要这些线条 积累总数。
- 终于来了
END 舞台我们打印 残余的seqlen。
评论代码：
awk '/^>/ { # header pattern detected
if (seqlen){
# print previous seqlen if exists
print seqlen
}
# pring the tag
print
# initialize sequence
seqlen = 0
# skip further processing
next
}
# accumulate sequence length
{
seqlen += length($0)
}
# remnant seqlen if exists
END{if(seqlen){print seqlen}}' file.fa
一个 oneliner：
awk '/^>/ {if (seqlen){print seqlen}; print ;seqlen=0;next; } { seqlen += length($0)}END{print seqlen}' file.fa
对于总数：
awk '/^>/ { if (seqlen) {
print seqlen
}
print
seqtotal+=seqlen
seqlen=0
seq+=1
next
}
{
seqlen += length($0)
}
END{print seqlen
print seq" sequences, total length " seqtotal+seqlen
}' file.fa
我想与klashxx的答案分享一些可能有用的调整。它的输出不同之处在于它在一行上打印序列id及其长度，它不再是单行，所以缺点是你必须将它保存为脚本文件。
它还根据空白区域从标题行解析序列id（chrM 在 >chrM gi|251831106|ref|NC_012920.1|）。然后，您可以通过设置变量来选择基于id的特定序列 target 像这样： $ awk -f seqlen.awk -v target=chrM seq.fa。
BEGIN {
OFS = "\t"; # tab-delimited output
}
# Use substr instead of regex to match a starting ">"
substr($0, 1, 1) == ">" {
if (seqlen) {
# Only print info for this sequence if no target was given
# or its id matches the target.
if (! target || id == target) {
print id, seqlen;
}
}
# Get sequence id:
# 1. Split header on whitespace (fields[1] is now ">id")
split($0, fields);
# 2. Get portion of first field after the starting ">"
id = substr(fields[1], 2);
seqlen = 0;
next;
}
{
seqlen = seqlen + length($0);
}
END {
if (! target || id == target) {
print id, seqlen;
}
}