关键词: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命令如下:
四种方式的最终结果是一致的: