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

david393

已于 2022-10-02 15:23:20 修改

阅读量743

点赞数

分类专栏：三代序列比对文章标签： c++ 开发语言算法

于 2022-09-12 21:16:31 首次发布

原文链接：https://mp.weixin.qq.com/s/2cOTKsE_2rPNeMGuPtAefA

版权

三代序列比对专栏收录该内容

10 篇文章 5 订阅

订阅专栏

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

不管做生信算法还是分析数据，都会涉及到对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;
}

然后编译：