生信（八）zlib库操作fq-gz文件

关键词：fq; gz; zlib

近期感谢yongzhe同学的需求，让我有机会能够用c来实操fq.gz的处理。
具体需求很简单：
输入一个index，将fq1和fq2（两个都是gz文件）中能够匹配该index的reads输出。输出文件也要是gz格式。
假设输入的index是ACCGAATG，那么下图中红色框中的字符串需要与输入index匹配才会将那条reads输出。

其实对fq.gz文件的处理，lh3（李恒）大牛已经写过一个头文件kseq.h（http://lh3lh3.users.sourceforge.net/kseq.shtml）。“拿来主义”当然没问题，但是对于新手而言，造轮子也是一个很好的进步途径。所以，我决定自己用c写。

C版本说明
fq.gz是gzip压缩过的文件，想要读取其内容需要用到zlib库（http://www.zlib.net/manual.html），我的程序中主要用到了gzopen、gzgetc、gzeof、gzerror以及gzclose这几个函数，分别是打开文件、读取一个字符、判断是否到文件末尾、获取出错信息、关闭文件的功能。

关于读取文件内容到内存中，不同于lh3使用了一个buffer数组，我只是简单使用了gzgetc函数。

另外，这次写getline类型的函数额外注意了EOF的处理。以前一直默认EOF之前一个字符应该是’\n’。这一次将EOF前一个字符不是’\n’的情况也考虑进去了。

最后，由于用到了zlib库，编译的时候要加上-lz参数。比如这样：

性能：
据yongzhe同学反映，c的版本所花时间大约是python版本的1/5。

具体的C代码
声明部分：

#include <stdio.h>
#include <zlib.h>
#include <stdlib.h>
#include <string.h>
#define RLINES 4        // set to be >= 1
#define MAXREADS 1000struct Reads;
typedef struct Reads *preads;
int getGzLine(gzFile gfile, char* line);
preads getReads(gzFile gfile, preads rd);
int cmpIndex(const char* reads, const int nl, const char* index, const int nx);struct Reads {int len;        // length of whole reads.int len1;       // length of 1st line of reads.char s[MAXREADS];
};

主程序：

 int main(int argc, char* argv[]) {gzFile infq1, infq2, outfq1, outfq2;preads reads1, reads2;const char* index;int m, err;if (argc < 6) {fprintf(stderr, "Usage: %s <index> <in-fq1> <in-fq2> <out-fq1> <out-fq2>\n", argv[0]);exit(1);}index = argv[1];m = strlen(index);if ((infq1 = gzopen(argv[2], "rb")) == NULL) {fprintf(stderr, "Error: cannot open %s\n", argv[2]);exit(2);}if ((infq2 = gzopen(argv[3], "rb")) == NULL) {fprintf(stderr, "Error: cannot open %s\n", argv[3]);exit(3);}if ((outfq1 = gzopen(argv[4], "wb")) == NULL) {fprintf(stderr, "Error: cannot open %s\n", argv[4]);exit(4);}if ((outfq2 = gzopen(argv[5], "wb")) == NULL) {fprintf(stderr, "Error: cannot open %s\n", argv[5]);exit(5);}if ((reads1 = (preads) malloc(sizeof(struct Reads))) == NULL ||(reads2 = (preads) malloc(sizeof(struct Reads))) == NULL) {fputs("Error: out of space!\n", stderr);exit(6);}while (getReads(infq1, reads1) != NULL && getReads(infq2, reads2) != NULL) {if (cmpIndex(reads1->s, reads1->len1 - 1, index, m)) {  // match Indexif (! gzwrite(outfq1, reads1->s, reads1->len)) {fprintf(stderr, "Error: %s\n", gzerror(outfq1, &err));exit(7);}if (! gzwrite(outfq2, reads2->s, reads2->len)) {fprintf(stderr, "Error: %s\n", gzerror(outfq2, &err));exit(8);}}}// further check lines-of-fq1 == lines-of-fq2free(reads1);free(reads2);gzclose(infq1);gzclose(infq2);gzclose(outfq1);gzclose(outfq2);
}

从gz文件中读取一行（包括’\n’）

int getGzLine(gzFile gfile, char* line) {
/* get one line with '\n'* return length of line. 0 if just a EOF in one line.*/char c;int n = 0;int err;while ((c = gzgetc(gfile)) != -1 && c != '\n')line[n++] = c;if (c == '\n') {line[n++] = '\n';line[n] = '\0';return n;} else if (gzeof(gfile)) {if (n == 0)return 0;else {line[n++] = '\n';       // add a '\n' to the end of file.line[n] = '\0';return n;}} else {        // Errorfprintf(stderr, "Error: %s\n", gzerror(gfile, &err));exit(101);}
}

读取一个reads序列：

 preads getReads(gzFile gfile, preads rd) {
/* return NULL if EOF*/int nline, nchar;char* p = rd->s;nline = nchar = 0;if ((nchar = getGzLine(gfile, p)) == 0) // EOFreturn NULL;rd->len1 = rd->len = nchar;p += nchar;while (++nline < RLINES && (nchar = getGzLine(gfile, p)) > 0) {rd->len += nchar;p += nchar;}if (nline < RLINES) {   // EOFfputs("Error: wrong fq\n", stderr);     // Need More Error Msgexit(201);}return rd;
}

比较reads的index是否匹配输入index：

 int cmpIndex(const char* reads, const int nl, const char* index, const int nx) {
// index located between the last ':' and the '\n'.int i, j;for (i = nl - 1, j = nx - 1; i >= 0 && j >= 0 && reads[i] == index[j]; i--, j--);if (j < 0 && i >= 0 && reads[i] == ':')return 1;elsereturn 0;
}

补充：仅处理一个fq的话
如果仅处理一个fq.gz文件，即仅打印fq1或fq2中匹配index的reads，可以这样做：
（假设要处理的gz文件是test.fq.gz，index序列是ACCGAATG）
使用grep –A命令：

zcat test.fq.gz | grep –A 3 ‘:ACCGAATG$’ | gzip –c > out1.fq.gz

或者用sed命令

zcat test.fq.gz | sed –n ‘/:ACCGAATG$/{N;N;N;p}’ | gzip –c > out2.fq.gz

下面两种方式更准确：
用sed命令：

zcat test.fq.gz | sed –n ‘h;N;N;N;x;/:ACCGAATG$/{x;p}’ | gzip –c > out3.fq.gz

或者用awk命令。由于命令较长，所以写入了一个脚本文件，命名为index.awk。

zcat test.fq.gz | awk –f index.awk | gzip –c > out4.fq.gz

具体的awk命令如下：

 #!/usr/bin/awk{f[0] = $0;for (i = 1; i < 4; i++) {getline;f[i] = $0;}if (f[0] ~ /:ACCGAATG$/) {for (i = 0; i < 4; i++)print f[i];}
}

四种方式的最终结果是一致的：

公众号：生信了

生信（八）zlib库操作fq-gz文件相关推荐

python使用 docx 库操作 docx 格式文件
docx 库文章结构: 一.docx 基本用,创建 docx 文件并添加数据二.深入理解文本格式(format),并设置所格式属性(attribute) 三.深入理解样式(styles),以及如何 ...
bam文件转fq.gz文件
bam转fq过程 1.fq文件格式 2.bam文件格式 3.转换思路 3.1 软件bedtools自带功能 3.2 自己写代码 3.3 代码示例 4.参考资料 1.fq文件格式 fastq格式是一 ...
Linux读取下机数据.fq.gz文件
#读取压缩文件 $ zcat R1.fq.gz 会唰唰唰全部读取,按Ctrl + C终止 #只查看部分 $ less R1.fq.gz 按↓键可往下翻页,按q退出
用 python 来操作 docx， xlsx 格式文件（二）（使用 docx 库操作 docx 格式文件
docx 库文章结构: 一.docx 基本用,创建 docx 文件并添加数据二.深入理解文本格式(format),并设置所格式属性(attribute) 三.深入理解样式(styles),以及如何 ...
python-docx库_用 python 来操作 docx， xlsx 格式文件（二）（使用 docx 库操作 docx 格式文件...
docx 库文章结构: 一.docx 基本用,创建 docx 文件并添加数据二.深入理解文本格式(format),并设置所格式属性(attribute) 三.深入理解样式(styles),以及如何 ...
[Zlib]_[初级]_[使用zlib库解压提取文件]
场景: 1. zlib库跨平台,Windows和MacOSX都可以使用,还支持64位编译,轻量级,没有不用的道理. 2. 处理.zip,apk,docx文件时,因为这类文件都是zip格式,使用zlib ...
使用zlib库解压*.zip文件
1.编辑zlib静态库 1)实现zip文件压缩和解压的zlib库部分代码(基于zlib-1.2.5) 这是zlib-1.2.5源码:http://download.csdn.net/download/ ...
送书 | 令附生信专用简明 Python 文字和视频教程
欢迎来到Python的世界,本教程将带你遨游Python,领悟Python的魅力.本教程专注于帮助初学者,尤其是生物信息分析人员快速学会Python的常用功能和使用方式,因此只精选了部分Python的 ...
生信人的自我修养：Linux 命令速查手册
标题:生信人的自我修养:Linux 命令速查手册目标:致力于为生信人打造一个完整的 Linux 命令速查手册作者:简佐义(jianzuoyi@qq.com) 版本:1.0 日期:2020-11-2 ...

生信（八）zlib库操作fq-gz文件

生信（八）zlib库操作fq-gz文件相关推荐

最新文章

热门文章