数据结构与算法

C语言实战课-让程序能够处理gz文件

之前的程序不能够处理压缩文件，而事实上为了节约空间，基本上fastq都会压缩成gz格式，因此这一课就是让程序能够支持压缩文件。

这里有两种思路，一种是利用管道，将之前的压缩文件通过zcat程序读取然后传递给我们的程序，另一种则是在程序中调用zlib库，让程序能够直接处理gz文件.

背景知识

针对思路1: 系统在每个C语言运行的时候都至少会提供三个流，标准输入(stdin)，标准输出(stdout)和标准错误(stderr). 在之前的程序我们就用到了stdout，用于将结果输出到屏幕上, 即fprintf(stdout, "%s", line);。 同样，我们可以修改fgets(line, MAX_LINE_LENGTH, fi) 中的fi为stdin使得程序能够接受管道传递的数据。

针对思路2: 利用zlib读取gz文件并不复杂，只需要将原来的函数名前面加上或者改成gz。毕竟优秀的代码应该符合人的直觉，下面就是我们将要用到的几个函数

  • 声明水管: gzFile
  • 连接水管: gzopen(const char *path, const char *mode)
  • 读取一行: gzgets(gzFile file, char *buf, int len)
  • 输出到gz文件中: gzprintf(gzFile file, const char *format, ...)
  • 关闭水管: gzclose(gzFile file)

我们以一个非常简单的代码进行展示，它会从gz文件中按行读取然后输出

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#include <stdio.h>
#include <zlib.h>

int main(int argc, char *argv[])
{

	gzFile fi;
	fi = gzopen(argv[1], "r");
	char string[1000];

	while ( gzgets(fi, string, 1000) != NULL) {
		printf("%s", string);
	}
	return 0;
}

它不是只能处理压缩文件，对于没有压缩的文件也是可以识别的，因此不需要写专门的逻辑判断语句，根据文件名来决定是否使用zlib相关的函数。

gzopen can be used to read a file which is not in gzip format; in this case gzread will directly read from the file without decompression. When reading, this will be detected automatically by looking for the magic two-byte gzip header

由于zlib不是标准库，我们在编译程序的时候需要额外加上-lz参数，用于链接动态库。

动手搞事情

有了以上的背景知识后，我们就可以来修改我们之前的代码了.

zlib在大部分服务器上都属于默认安装，因此只需要多加一行#include <zlib.h>. 但是还是有一小部分服务器居然没有安装zlib，因此我们需要先下载编译（root权限可以直接用apt/yum进行安装）。

1
2
3
wget https://www.zlib.net/zlib-1.2.11.tar.gz
tar xf zlib-1.2.11.tar.gz
cd zlib-1.2.11 && ./configure && make -j 8

修改后的fq2fa.c的代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "zlib.h"

#define MAX_LINE_LENGTH 1000

int main(int argc, char *argv[]){

	gzFile fi;
	gzFile fo;
	if (argc < 2 ){
		return -1;
	}
	fi = gzopen( argv[1], "rb");
	if ( argc == 3 ){
	    fo = gzopen (argv[2], "wb");
	}

	char line[MAX_LINE_LENGTH];
	int count = 0;
	while ( gzgets(fi,line, MAX_LINE_LENGTH) != NULL) {
		if ( count == 0) {
			line[0] = '>';
		}
		if ( count < 2){
			if ( argc == 2) {
				fprintf(stdout, "%s", line);
			}
			if ( argc ==3 ){
			    gzprintf(fo, "%s", line);
			}
		}
		if ( ++count > 3 ){
			count = 0;
		}

	}
	gzclose(fi);
	if (argc == 3) gzclose(fo);
	return 0;
}

编译方法为gcc -lz -Lzlib -o fq2fa fq2fa.c。关于代码的讲解和zlib的知识点可以看视频介绍, https://www.bilibili.com/video/av84768292/。

代码目前存在的问题是行缓冲是1000字符，对于三代测序数据，按照宣传可以有几M，一般也是几百k，因此当前程序处理三代程序肯定不行。一个方法是将MAX_LINE_LENGTH设置为更大的值，比如说设置为10M，或者是将其设置为一个参数，可以手动指定。

不过下一次，将会学习一个新的轮子klib，用里面的函数来解决这个问题。