使用 SED/AWK 将 FASTQ 转换为 FASTA

StackOverflow https://stackoverflow.com/questions/1542306

  •  20-09-2019
  •  | 
  •  

题

我有一个数据，它总是以四个块的形式出现 采用以下格式（称为 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 块中的前两行 并替换 @ 和 >.

有帮助吗？

解决方案

这是一个老问题，也有过许多提供不同的解决方案。由于接受的答案使用SED，但有一个明显的问题（这是它会与>当@符号出现质量行的第一个字母代替@），我觉得有必要提供一个简单的sed式的解决方案，实际工作：

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

做的唯一的假设是，每一次读中占有FASTQ文件正好是4行，但似乎非常安全的，在我的经验。

在fastx工具包中的fastq_to_fasta脚本也有效。 （值得一提的是，你需要指定-Q33选项，以适应现在常见的PHRED + 33 QUAL编码。这很有趣，因为它反正扔掉质量数据！）

其他提示

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

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

或者，仿真 HTTP：//www.ringtail.tsl .ac.uk /大卫-studholme /脚本/ fastq2fasta.pl 张贴由Pierre，只从第一行打印的第一个字（的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）NAR，许多这些解决方案是自“的‘@’标记字符（ASCII 64）也可以在质量串的任何地方发生不正确。这意味着，任何解析器必须不治疗线开始的“@”为指示下一个记录开始时，无需额外检查质量字符串的长度迄今为止相匹配的序列的长度。“

请参阅 http://ukpmc.ac.uk/articlerender.cgi?accid=PMC2847217 ，获取的信息。

刚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

我会写

awk '
    NR%4 == 1 {print ">" substr($0, 2)}
    NR%4 == 2 {print}
' fastq > fasta

这是我得最快的，我把它贴在我的.bashrc文件：

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

它不会失败的一个以@ ...开始但并未能对包裹FASTQ不常见的，但不是不可能的音质线，如果这甚至法律（尽管它存在）。

下面是解决，我只是从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

将做的工作。

是这样的：

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

应该工作。

我认为，与GNU grep的，这可能与此来完成：

grep -A 1 "^@" t.txt | grep -v "^--" | sed -e "s/^@/\>/"

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

，其中数据是数据文件。 我已经收到：

>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/>/

我知道我的未来，但对于Google员工的利益：

您可能需要使用 fastq_to_fasta从fastx工具。这将保持@符号，虽然。它也将删除与NS线，除非你告诉它不要。

您可能对 bioawk 感兴趣，它是 awk 的改编版本，经过调整可以处理 fasta 文件

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

笔记： 生物Awk 是基于 Brian Kernighan 的 awk 这记录在 “AWK 编程语言”， 作者：Al Aho、Brian Kernighan 和 Peter Weinberger （Addison-Wesley，1988年，ISBN 0-201-07981-X） . 。我不确定这个版本是否兼容 POSIX.

许可以下： CC-BY-SA 和 归因
不隶属于 StackOverflow
scroll top