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舞台上我们打印了残余的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