不管做生信算法还是分析数据,都会涉及到对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 #include #include #include "kseq.h" // STEP 1: declare the type of file handler and the read() functionKSEQ_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 \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文件格式或压缩形式的文件,速度差别会更大。
参考网址:
https://www.cnblogs.com/xudongliang/p/5307462.html
http://lh3lh3.users.sourceforge.net/kseq.shtml
往期回顾:还没开学?经典生物信息书免费送!免费送!免费送!总价值4170元!
Tax4Fun: 16S 预测群落代谢功能实例——手把手示范3-差异分析
- Tax4Fun: 16S 预测群落代谢功能实例——手把手示范2
- Tax4Fun: 16S 预测群落代谢功能实例——手把手示范1
- 三代序列比对算法minimap2,就是快,快,快!
- 测序数据过滤神器Trimmomatic:去除接头和低质量序列
- 基于三代测序数据的结构变异检测,PBHoney方法解读
- usearch序列聚类地位不保?看看Linclust算法,快你没商量!
- BMC Bioinformatics期刊三审被拒,如此悲催,如何避免?
- NCBI数据下载速度50M/S以上?太快了吧
- rHAT,国内首个三代序列比对算法
- 扩增子OTU聚类软件:SeekDeep方法解读(NAR方法文章)
- mothur QIIME usearch,三足鼎立,谁主沉浮?
- 三代测序序列比对利器-BLASR,更小更快更方便
- 生信算法“八股文”,发表算法不再难!