生信（八）zlib库操作fq-gz文件

关键词：fq; gz; zlib

近期感谢yongzhe同学的需求，让我有机会能够用c来实操fq.gz的处理。

具体需求很简单：

输入一个index，将fq1和fq2（两个都是gz文件）中能够匹配该index的reads输出。输出文件也要是gz格式。

假设输入的index是ACCGAATG，那么下图中红色框中的字符串需要与输入index匹配才会将那条reads输出。

其实对fq.gz文件的处理，lh3（李恒）大牛已经写过一个头文件kseq.h（http://lh3lh3.users.sourceforge.net/kseq.shtml）。“拿来主义”当然没问题，但是对于新手而言，造轮子也是一个很好的进步途径。所以，我决定自己用c写。

C版本说明

fq.gz是gzip压缩过的文件，想要读取其内容需要用到zlib库（http://www.zlib.net/manual.html），我的程序中主要用到了gzopen、gzgetc、gzeof、gzerror以及gzclose这几个函数，分别是打开文件、读取一个字符、判断是否到文件末尾、获取出错信息、关闭文件的功能。

关于读取文件内容到内存中，不同于lh3使用了一个buffer数组，我只是简单使用了gzgetc函数。

另外，这次写getline类型的函数额外注意了EOF的处理。以前一直默认EOF之前一个字符应该是’\n’。这一次将EOF前一个字符不是’\n’的情况也考虑进去了。

最后，由于用到了zlib库，编译的时候要加上-lz参数。比如这样：

性能：

据yongzhe同学反映，c的版本所花时间大约是python版本的1/5。

具体的C代码

声明部分：

主程序：

从gz文件中读取一行（包括’\n’）

读取一个reads序列：

比较reads的index是否匹配输入index：

补充：仅处理一个fq的话

如果仅处理一个fq.gz文件，即仅打印fq1或fq2中匹配index的reads，可以这样做：

（假设要处理的gz文件是test.fq.gz，index序列是ACCGAATG）

使用grep –A命令：

zcat test.fq.gz | grep –A 3 ‘:ACCGAATG$’| gzip –c > out1.fq.gz

或者用sed命令

zcat test.fq.gz| sed –n ‘/:ACCGAATG$/{N;N;N;p}’ | gzip –c > out2.fq.gz

下面两种方式更准确：

用sed命令：

zcat test.fq.gz| sed –n ‘h;N;N;N;x;/:ACCGAATG$/{x;p}’ | gzip –c > out3.fq.gz

或者用awk命令。由于命令较长，所以写入了一个脚本文件，命名为index.awk。

zcat test.fq.gz| awk –f index.awk | gzip –c > out4.fq.gz

具体的awk命令如下：

四种方式的最终结果是一致的：