如何使用kseq.h解析FASTA文件

Question

我已经知道这个来自Heng Li的图书馆已经有一段时间了,但直到现在我还没有尝试过使用它,主要是因为到目前为止python已经足够快了.

以下是标题的链接:http://lh3lh3.users.sourceforge.net/kseq.shtml

当我尝试使用以下内容来解析fasta文件时,它返回-1作为序列行的长度.我查看了Li的代码,这似乎主要是为FASTQ解析而设计的,但他确实在他的网页上说它也支持FASTA格式.

这是我的代码:

#include <stdio.h>
#include <stdlib.h>
#include "kseq.h"  
// STEP 1: declare the type of file handler and the read() function  
KSEQ_INIT(FILE*, read)


int main(int argc, char** argv) {
    FILE* fp = fopen(argv[1], "r"); // STEP 2: open the file handler
    kseq_t *seq = kseq_init(fp); // STEP 3: initialize seq 

    int l;

    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);  
    }

    printf("return value: %d\n", l);  
    kseq_destroy(seq); // STEP 5: destroy seq
    fclose(fp);

    return (0);
}

我一直用来测试的FASTA是Hg19 GRCH37 ChrY.fa文件,可从包括Broad Institute在内的多个来源获得.

任何帮助,将不胜感激.

Answer 1

首先你应该检查fopen()的返回值:

FILE* fp = fopen(argv[1], "r"); // STEP 2: open the file handler
if(fp == 0) {
    perror("fopen");
    exit(1);
}

其次,我查看了头文件,我认为kseg_init采用的是fd而不是FILE*.您可以使用fileno()从FILE*获取fd.

kseq_t *seq = kseq_init(fp); // STEP 3: initialize seq 

应该:

kseq_t *seq = kseq_init(fileno(fp)); // STEP 3: initialize seq