18

我有一个数据，它总是以以下格式（称为 FASTQ）以四个为一组：

@SRR018006.2016 GA2:6:1:20:650 length=36
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNGN
+SRR018006.2016 GA2:6:1:20:650 length=36
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!+!
@SRR018006.19405469 GA2:6:100:1793:611 length=36
ACCCGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+SRR018006.19405469 GA2:6:100:1793:611 length=36
7);;).;);;/;*.2>/@@7;@77<..;)58)5/>/

是否有一种简单的 sed/awk/bash 方法可以将它们转换成这种格式（称为 FASTA）：

>SRR018006.2016 GA2:6:1:20:650 length=36
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNGN
>SRR018006.19405469 GA2:6:100:1793:611 length=36
ACCCGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

原则上，我们希望提取每个 4 块中的前两行并替换@为>.

4

13 回答 13

26

这是一个古老的问题，并且提供了许多不同的解决方案。由于接受的答案使用 sed，但有一个明显的问题（即当 @ 符号作为质量行的第一个字母出现时，它将用 > 替换 @），我觉得有必要提供一个基于 sed 的简单解决方案，它确实有效：

sed -n '1~4s/^@/>/p;2~4p'

唯一的假设是每次读取在 FASTQ 文件中恰好占据 4 行，但根据我的经验，这似乎很安全。

fastx 工具包中的 fastq_to_fasta 脚本也可以使用。（值得一提的是，您需要指定 -Q33 选项以适应现在常见的 Phred+33 质量编码。这很有趣，因为它无论如何都会丢弃质量数据！）

于 2012-04-28T00:19:07.207 回答
9

sed 没有死。如果我们在打高尔夫球：

sed '/^@/!d;s//>/;N'

或者，模仿Pierre 发布的http://www.ringtail.tsl.ac.uk/david-studholme/scripts/fastq2fasta.pl，它只打印第一行的第一个单词（id）并且（一些）错误处理：

#!/usr/bin/sed -f
# Read a total of four lines
$b error
N;$b error
N;$b error
N
# Parse the lines
/^@\(\([^ ]*\).*\)\(\n[ACGTN]*\)\n+\1\n.*$/{
  # Output id and sequence for FASTA format.
  s//>\2\3/
  b
}
:error
i\
Error parsing input:
q

似乎有很多现有的工具可以转换这些格式；您可能应该使用这些而不是此处发布的任何内容（包括上述内容）。

于 2009-10-09T14:45:31.543 回答
9

正如 Cock, et al (2009) NAR 中详述的那样，这些解决方案中的许多都是不正确的，因为“'@' 标记字符 (ASCII 64) 可能出现在质量字符串中的任何位置。这意味着任何解析器都不得处理以'@' 表示下一条记录的开始，没有额外检查质量字符串的长度到目前为止是否匹配序列的长度。”

有关详细信息，请参阅http://ukpmc.ac.uk/articlerender.cgi?accid=PMC2847217。

于 2010-05-09T20:06:48.470 回答
7

只需要awk，不需要其他工具

# awk '/^@SR/{gsub(/^@/,">",$1);print;getline;print}' file
>SRR018006.2016 GA2:6:1:20:650 length=36
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNGN
>SRR018006.19405469 GA2:6:100:1793:611 length=36
ACCCGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
于 2009-10-09T11:57:09.287 回答
4

请参阅http://www.ringtail.tsl.ac.uk/david-studholme/scripts/中的fastq2fasta.pl

于 2009-10-09T07:45:23.820 回答
3

我会写

awk '
    NR%4 == 1 {print ">" substr($0, 2)}
    NR%4 == 2 {print}
' fastq > fasta
于 2011-06-30T22:33:48.867 回答
2

这是我最快的，我把它放在我的 .bashrc 文件中：

alias fq2fa="awk '{print \">\" substr(\$0,2);getline;print;getline;getline}'"

它不会在以 @... 开头的罕见但并非不可能的质量行上失败，但在包装的 FASTQ 上确实会失败，如果这甚至是合法的（尽管它存在）。

于 2011-10-27T23:06:31.307 回答
2

您可能对 bioawk 感兴趣，它是 awk 的改编版本，用于处理 fasta 文件

bioawk -c fastx '{ print ">"$name ORS $seq }' file.fastq

注意： BioAwk基于Brian Kernighan 的 awk，该 awk 记录在Al Aho、Brian Kernighan 和 Peter Weinberger 的“The AWK Programming Language”中（Addison-Wesley，1988，ISBN 0-201-07981-X） 。我不确定这个版本是否与POSIX兼容。

于 2018-12-12T11:27:05.597 回答
1 
awk 'BEGIN{P=1}{if(P==1||P==2){gsub(/^[@]/,">");print}; if(P==4)P=0; P++}' data

>SRR018006.2016 GA2:6:1:20:650 length=36
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNGN
>SRR018006.19405469 GA2:6:100:1793:611 length=36
ACCCGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

以下

awk '{gsub(/^[@]/,">"); print}' data

其中 data 是您的数据文件。我收到了：

>SRR018006.2016 GA2:6:1:20:650 length=36
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNGN
+SRR018006.2016 GA2:6:1:20:650 length=36
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!+!
>SRR018006.19405469 GA2:6:100:1793:611 length=36
ACCCGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+SRR018006.19405469 GA2:6:100:1793:611 length=36
7);;).;);;/;*.2>/@@7;@77<..;)58)5/>/
于 2009-10-09T07:38:33.573 回答
1

这是我刚刚从 SO 中学到的问题的“跳过所有其他行”部分的解决方案：

while read line
do
    # print two lines
    echo "$line"
    read line_to_print
    echo "$line_to_print"

    # and skip two lines
    read line_to_skip
    read line_to_skip
done

如果需要做的只是将一个更改@为>，那么我认为

while read line
do
    echo "$line" | sed 's/@/>/'
    read line
    echo "$line"

    read line_to_skip
    read line_to_skip
done

将完成这项工作。

于 2009-10-09T07:39:06.273 回答
1

就像是：

awk 'BEGIN{a=0}{if(a==1){print;a=0}}/^@/{print;a=1}' myFastqFile | sed 's/^@/>/'

应该管用。

于 2009-10-09T07:39:40.920 回答
1

我认为，使用 gnu grep 可以这样做：

grep -A 1 "^@" t.txt | grep -v "^--" | sed -e "s/^@/\>/"
于 2009-10-09T07:40:23.797 回答
1

我知道我在未来，但为了谷歌员工的利益：

您可能想使用fastx 工具包中的 fastq_to_fasta。不过，它将保留@ 符号。它还将删除带有 Ns 的行，除非您告诉它不要这样做。

于 2011-06-30T19:21:45.140 回答