3

我从 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 在内的多个来源获得。

任何帮助，将不胜感激。

4

2 回答 2

7

首先你应该检查 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
于 2013-10-15T20:56:58.427 回答
2

这是为我工作的完整代码

#include <zlib.h>
#include <stdio.h>
#include "kseq.h"
KSEQ_INIT(int, read)

int main(int argc, char **argv)
{
        FILE* fp;
        kseq_t *seq;
        int n = 0, slen = 0, qlen = 0;
        fp = fopen(argv[1], "r");
        seq = kseq_init(fileno(fp));
        while (kseq_read(seq) >= 0)
                ++n ;//slen += seq->seq.l, qlen += seq->qual.l;
        printf("%d\t%d\t%d\n", n, slen, qlen);
        kseq_destroy(seq);
        fclose(fp);
        return 0;
}
于 2018-01-26T03:53:15.077 回答