不管做生信算法还是分析数据,都会涉及到对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文件格式或压缩形式的文件,速度差别会更大。