bash - FASTA文件的序列长度

Question

我有以下 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 中很容易做到，实际上，我有一个脚本可以用这些方式做到这一点。

Answer 1

awk/gawk 解决方案可以由三个阶段组成:

1. 每次找到 header 时，都应执行以下操作:

   • 如果存在，打印上一个seqlen。
   • 打印标签。
   • 初始化 seqlen。
2. 对于 sequence 行，我们只需要累加总数。
3. 最后，在 END 阶段，我们打印了 remnant 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