c - 如何使用 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