kseq.h读取序列文件，快到想不到！

本文来自“生信算法”公众号。 

       不管做生信算法还是分析数据，都会涉及到对FASTA或FASTQ序列文件进行读取再处理，我们可以自己写脚本或程序进行读取，如果有现成快速的读取软件包或者库函数，那是最好不过了。今天给大家介绍的就是C语言编写的一个头文件，可以快速读取序列文件，尤其是对比较大的文件，速度是相当的快。

        之前读取序列文件主要用自己写的Python脚本统计这些信息， 但当序列文件非常大时，就变的很慢。通过网上不断搜索或看一些序列算法的源代码后，发现好多优秀的算法都用kseq.h来读取序列，简单好用，且读取速度快，下面我们看一下并比较一下和我自己写的Python脚本的读取速度。

        首先下载kseq.h头文件，这个不是C语言默认的头文件，需要自己下载放在文件夹下引用。如下图所示：

        下载头文件kseq.h请点击网站进行下载（复制后粘贴即可）：http://lh3lh3.users.sourceforge.net/kseq.shtml

        记得一定保存为kseq.h头文件，不能换成其他名字和文件类型。

        然后我们编写一个C程序如下，文件名为main.c，记得一定要包含kseq.h头文件（代码中第4行）：

#include <zlib.h>  
#include <stdio.h>
#include <string.h>  
#include "kseq.h"  
// STEP 1: declare the type of file handler and the read() function
KSEQ_INIT(gzFile, gzread)
  
int main(int argc, char *argv[])
{
    gzFile fp;
    kseq_t *seq;
    long seqs = 0;
    long bases = 0;
    int l;
    if (argc == 1) {
        fprintf(stderr, "Usage: %s <in.seq>\n", argv[0]);
        return 1;
    }
    fp = gzopen(argv[1], "r"); // STEP 2: open the file handler
    seq = kseq_init(fp); // STEP 3: initialize seq
    while ((l = kseq_read(seq)) >= 0) { // STEP 4: read sequence
        //printf("name: %s\n", seq->name.s);
        //if (seq->comment.l) printf("comment: %s\n", seq->comment.s);
        //printf("seq: %s\n", seq->seq.s);
        //if (seq->qual.l) printf("qual: %s\n", seq->qual.s);
        bases += strlen(seq->seq.s);
        seqs += 1;
    }
    //printf("return value: %d\n", l);
    printf("reads: %ld\n", seqs);
    printf("bases: %ld\n", bases);
    kseq_destroy(seq); // STEP 5: destroy seq
    gzclose(fp); // STEP 6: close the file handler
    return 0;
}

        然后编译：

gcc -o fastx_read_length -lz mai.c

        因为调用zlib读取压缩文件，所以编译时需要添加-lz 选项；编译后生成的可运行文件为：fastx_read_length，然后在当前文件夹下运行：

./fastx_read_length

        我们用两个例子看一下kseq.h读取文件的速度怎么样，和我自己编的Python脚本比较一下：

        读取111万条数据只用了1.05秒的时间。

而我自己编写的Python脚本用了3.9s。

我们再看看读取全长16S序列，319万条序列，平均长度1400bp左右，文件大小3.8G。

利用kseq.h读取大约需要14S的时间。

而我自己写的Python需要43S。

这只是FASTA文件格式，如果是FASTQ文件格式或压缩形式的文件，速度差别会更大。

参考网址：

1. https://www.cnblogs.com/xudongliang/p/5307462.html
2. http://lh3lh3.users.sourceforge.net/kseq.shtml