一句代码快速提取fastq文件中的fasta

NGS数据输出一般为fastq格式，包含序列的质量等信息，有时候我们只想提取fasta序列文件，可以通过多种NGS序列处理软件。此外，可以使用一行代码快速提取。

可以通过两种方式，第一种用sed，第二种用awd，个人觉得第1种sed 命令比较精妙，分享供大家参考：

复习一下典型的fastq文件格式

$ cat test.fq
@ST-E00126:128:HJFLHCCXX:2:1101:7405:1133
TTGCAAAAAATTTCTCTCATTCTGTAGGTTGCCTGTTCACTCTGATGATAGTTTGTTTTGG
+
FFKKKFKKFKF<KK<F,AFKKKKK7FFK77<FKK,<F7K,,7AF<FF7FKK7AA,7<FA,,

使用sed命令

读文献时看到这个命令sed '/^@/!d;s//>/;N'，顿时感到太精妙了，忍不住分享一下:

$ sed '/^@/!d;s//>/;N'  test.fq 
>ST-E00126:128:HJFLHCCXX:2:1101:7405:1133
TTGCAAAAAATTTCTCTCATTCTGTAGGTTGCCTGTTCACTCTGATGATAGTTTGTTTTGG

想要理解这个命令，首先需要理解sed 中N的用法：

N命令

N命令的介绍来自SED 简明教程，想详细了解sed的可以看看这篇文章

先来看N命令 —— 把下一行的内容纳入当成缓冲区做匹配。

演示文本为pets.text：

$ cat pets.txt
This is my cat
  my cat's name is betty
This is my dog
  my dog's name is frank
This is my fish
  my fish's name is george
This is my goat
  my goat's name is adam

下面的的示例会把原文本中的偶数行纳入奇数行匹配，而s只匹配并替换一次，所以，就成了下面的结果：

$ sed 'N;s/my/your/' pets.txt
This is your cat
  my cat's name is betty
This is your dog
  my dog's name is frank
This is your fish
  my fish's name is george
This is your goat
  my goat's name is adam

也就是说，原来的文件成了：

This is my cat\n  my cat's name is betty
This is my dog\n  my dog's name is frank
This is my fish\n  my fish's name is george
This is my goat\n  my goat's name is adam

这样一来，下面的例子你就明白了，

$ sed 'N;s/\n/,/' pets.txt
This is my cat,  my cat's name is betty
This is my dog,  my dog's name is frank
This is my fish,  my fish's name is george
This is my goat,  my goat's name is adam

这样，在命令sed '/^@/!d;s//>/;N' test.fq中，首先使/^@/’ 匹配开头为@的行，就是第一行，后面加了N就是说把开头@行和后面一行（序列行）算成1行，然后用/！反向匹配，即匹配开头是@行以外的行，d是删除，所以/^@/!d;N 这部分代码就是说把第一行第二行当成一行匹配，之外的行的全部删除，就是说原文件第3，4行被删除。然后s//>/将剩余的第一第二行中匹配的^@替换成">"。所以这一个命令使用了正则表达式（^）,sed 命令连用（；）反向匹配符号（！），还有sed 的N命令。

2. 第二个是使用awk命令，比较常规

$ awk '{if(FNR%4==1){sub("@",">",$0);print $0};if(FNR%4==2)print$0}' test.fq 
>ST-E00126:128:HJFLHCCXX:2:1101:7405:1133
TTGCAAAAAATTTCTCTCATTCTGTAGGTTGCCTGTTCACTCTGATGATAGTTTGTTTTGG