关键词: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代码
声明部分:
#include <stdio.h>
#include <zlib.h>
#include <stdlib.h>
#include <string.h>
#define RLINES 4 // set to be >= 1
#define MAXREADS 1000
struct Reads;
typedef struct Reads *preads;
int getGzLine(gzFile gfile, char* line);
preads getReads(gzFile gfile, preads rd);
int cmpIndex(const char* reads, const int nl, const char* index, const int nx);
struct Reads {
int len;