了解解开 fasta 文件的 awk 公式

Question

我刚刚找到了一个可用于解开 fasta 文件的公式。在给出公式之前，我需要解释一下什么是解包 fasta 文件。简而言之，fasta格式是这样的：

>name_of_sequence$
xxxxxxxxxxxxxxxxxxxxxx$
>name_of_sequence_2$
xxxxxxxxxxxxxxxxxxxxxx$
>name_of_sequence_3$
xxxxxxxxxxxxxxxxxxxxxx$

这将是一个普通的 fasta 文件，因为每个序列只有一行（xxxxxx ...）。美元符号是换行符。

但有时候，你会发现包裹fasta 文件是这样的：

>name_of_sequence$
xxxxxxxxx$
xxxxxxxxx$
xxxx$
>name_of_sequence_2$
xxxxxxxxx$
xxxxxxxxx$
xxxx$
>name_of_sequence_3$
xxxxxxxxx$
xxxxxxxxx$
xxxx$

在这里，您仍然只有三个序列，但每个序列都分为三个部分。解开 fasta 文件意味着将后一种格式转换为前一种格式（每个序列一行）。

为此，您需要从后一个文件中删除换行符，但不是全部。您需要在序列名称后面保留换行符（例如此处：>name_of_sequence$），并在序列末尾保留换行符（例如此处：xxxx$）。

看来这个公式是这样做的：

cat infasta | awk '/^>/{print s? s"\n"$0:$0;s="";next}{s=s sprintf("%s",$0)}END{if(s)print s}' > outfasta

我的问题是：有人可以向我解释它是如何工作的吗？

Answer 1

这是你的awk脚本：

/^>/ {
    print s ? s "\n" $0 : $0;
    s = "";
    next;
}

{
    s = s sprintf("%s", $0);
}

END {
    if (s)
      print s;
}

>仅当以 开头的行（即 fasta 标题行）触发第一个块。

在第一个块中，打印了一些内容。那东西是s ? s "\n" $0 : $0。这意味着“如果s非零（或未设置），则使用s并向其添加换行符，后跟整个当前行，否则仅使用整个当前行”。在此程序中，s将是属于最近处理的标题行的部分读取序列，当程序命中标题行时，此print语句将输出最后一个序列（现在已完成），如果有的话，后面跟着新发现的标题行位于新行上。

然后该块设置s为空字符串（我们还没有读取属于该标头的任何序列），然后跳到下一个输入行。

下一个块针对所有输入行执行（但不针对标题行，因为由于next上一个块中的 ，这些行将被跳过）。它只是将当前行附加到s. sprintf已使用，但我不太确定为什么（s = s $0也可能有效）。

最后一个块将在读取所有输入行后执行。它将打印属于最后一个标题行的序列（如果有）。

概括：

该awk脚本通过将所有单独的序列行保存在变量中来连接它们。当找到标题行时，它会在自己的一行上输出到目前为止读取的序列以及新标题。最后，输出属于最后一个标头的序列。

---

awk不在变量中存储序列的替代脚本（如果你的 fasta 文件中有非常大的基因组，可能会很有用）：

/^>/ {
    if (NR == 1) {
        print;  # 1st header line, just print it.
    } else {
        # Print a newline for the prev. sequence, then the header line on its own line.
        printf("\n%s\n", $0);
    }
    next; # Skip to next input line.
}

{
    printf("%s", $0); # Print sequence without newline.
}

END {
    printf("\n"); # Add final newline to output.
}

作为“一行”：

awk '/^>/{if(NR==1){print}else{printf("\n%s\n",$0)}next} {printf("%s",$0)} END{printf("\n")}' sequence.fasta

Answer 2

FWIW 提供了一个基于“sed”的解决方案来包装 fasta 文件。 sed 方法的基本流程是在找到序列名称行后，首先我们将该行单独显示在一行上，然后开始将序列行累积到模式空间本身中，并同时剥离换行符。当我们到达下一个序列名称行或 eof 时，此流程就会中断。

sed -e '
  /^>/{                  # caught sequence name line
     n                   # print seq name, next line into pattern space
     :loop
        N                # read next line into PS, if not print PS/quit
        /\n>/!s/\n//     # join successive sequences
     /\n/!bloop          # go back for more seq if new seq name not got yet
     P;D                 # print the current seq then delete it, branch to the top with PS having new seq name
  }
' your_fasta_file