生信（11）htslib处理bam或sam文件的简单示例

本文给出了一个示例，介绍如何使用htslib编写c程序来处理bam/sam文件。

（本文写于2020年初，随着将来htslib和samtools库的更新，本文部分内容可能会不适用，请读者注意官网的更新动态。）

我们通常会使用samtools软件来处理bam/sam文件。但有时候我们也需要对bam/sam文件进行一些个性化的处理，这个时候就需要自己编写程序。

从samtools的github官网上可以看到，原来的samtools已经被拆分成三个小项目，分别是htslib，samtools以及bcftools。其中htslib是一个处理高通量数据通用文件格式的库，是samtools软件和bcftools软件依赖的核心库。如果要用c/c++来操作bam/sam文件，一定要了解htslib。

最初samtools只有c语言版本，由于samtools软件的强大功能和广泛应用，许多编程语言都对samtools进行了封装。比如pysam就是python语言对samtools（当然也包括htslib库）的封装。pysam提供了一套完整的操作bam/sam的API以及对应的说明文档，上手比较容易。但是，htslib的官网上并没有像pysam那样完善详尽的说明文档，所以对于新手，往往不知道从哪里开始学习。

本文给出了一个示例，详细介绍基于htslib库编写c程序来操作bam/sam文件的完整步骤。按照这些步骤，大家可以编写并运行自己的第一个基于htslib的c/c++程序。与此同时，读完本文后大家也会了解可以通过学习htslib官网上的.h文件来熟悉这个库的API，通过学习.c文件来熟悉一些使用API的代码示例。具体可以分为以下几步：

• 下载安装htslib库
• 编写c程序
• 编译c程序
• 运行c程序

下载安装htslib库

方便起见，我们在家目录下新建一个项目文件夹，命名为demo。首先我们从官网上下载htslib库。在Linux中运行以下命令

mkdir ~/demo
cd ~/demo        # cd <path_to_project_dir>
git clone https://github.com/samtools/htslib.git    # download htslib

下载完成后在demo文件夹下会出现一个新文件夹，htslib。

接下来是安装htslib，按照官网上的安装步骤：（如果提示找不到autoheader或者autoconf，那么需要先安装autotools套件。如果安装过程中出现错误，最常见的就是缺少某些库，那么按照报错信息进行操作即可。）

cd ~/demo/htslib    # cd <path_to_htslib_project_dir>
autoheader   # If using configure, generate the header template...
autoconf     # ...and configure script (or use autoreconf to do both)
./configure  # Optional but recommended, for choosing extra functionality
make
make install # sudo make install

安装成功后htslib项目文件夹中应该会有以lib开头命名的文件。比如libhts.a文件和libhts.so文件。

编写c程序

安装好了htslib库之后，该如何用它操作bam/sam文件呢？

就像利用pysam提供的API去操作bam/sam文件一样，我们要利用htslib库，首先得熟悉这个库提供的API。由于htslib库没有提供详尽的API说明文档，所以我们只能去看源代码。其实，htslib这个库提供的API基本上都在一些头文件（.h文件）中。比如htslib库中的sam.h文件就包含了很多实用的API。很多头文件中的说明还是很详细的，多看看慢慢地就会熟悉了。

其次，可以通过阅读这个库中的一些.c文件来学习如何使用这些API。比如，htslib库中有一个test_view.c文件就给出了一个很好的使用API的示例（文末有链接）。由于test_view.c还是很长，笔者据此进行修改，写了一个更简单的示例，命名为samtest.c。这个程序的作用是从bam/sam文件中提取全部或者部分区域的比对结果（不包含头部信息）。代码如下：

#include <stdio.h>
#include <unistd.h>
#include <getopt.h>
#include <string.h>
#include "htslib/sam.h"

enum filetype {
	FBAM = 1,   // BAM file
	FSAM = 2,   // SAM file
};

int ftype;

int sam_test_extract(int argc, char **argv, int optind, htsFile *in, htsFile *out) {
		sam_hdr_t *hdr;
		bam1_t *b;
		hts_idx_t *idx = NULL;
		hts_itr_t *iter = NULL;
		int ret;

		if ((hdr = sam_hdr_read(in)) == NULL) {
			fprintf(stderr, "[E::%s] couldn't read header for '%s'\n", __func__, argv[optind]);
			return -1;
		}
		if ((b = bam_init1()) == NULL) {
			fprintf(stderr, "[E::%s] Out of memory allocating BAM struct.\n", __func__);
			goto fail;
		}
		if (ftype == FBAM && optind + 2 <= argc) { // BAM input and has a region.
			if ((idx = sam_index_load(in, argv[optind])) == 0) {
				fprintf(stderr, "[E::%s] fail to load the index for '%s'\n", __func__, argv[optind]);
				goto fail;
			}
			if ((iter = sam_itr_querys(idx, hdr, argv[optind + 1])) == 0) {
				fprintf(stderr, "[E::%s] fail to parse region '%s'\n", __func__, argv[optind + 1]);
				goto fail;
			}
			while ((ret = sam_itr_next(in, iter, b)) >= 0) {
				if (sam_write1(out, hdr, b) < 0) {
					fprintf(stderr, "[E::%s] Error writing output.\n", __func__);
					goto fail;
				}
			}
			if (ret < -1) {
				fprintf(stderr, "[E::%s] Error reading input.\n", __func__);
				goto fail;
			}
			hts_itr_destroy(iter);
			iter = NULL;
			hts_idx_destroy(idx);
			idx = NULL;
		} else if (optind + 2 > argc) {
			while ((ret = sam_read1(in, hdr, b)) >= 0) {
				if (sam_write1(out, hdr, b) < 0) {
					fprintf(stderr, "[E::%s] Error writing alignments.\n", __func__);
					goto fail;
				}
			}
			if (ret < -1) {
				fprintf(stderr, "[E::%s] Error parsing input.\n", __func__);
				goto fail;
			}
		} else { // SAM input and has a region.
			fprintf(stderr, "[E::%s] couldn't extract alignments directly from raw sam file.\n", __func__);
			goto fail;
		}
		bam_destroy1(b);
		sam_hdr_destroy(hdr);
		return 0;
	fail:
		if (iter) sam_itr_destroy(iter);
		if (b) bam_destroy1(b);
		if (idx) hts_idx_destroy(idx);
		if (hdr) sam_hdr_destroy(hdr);
		return 1;
}

int main(int argc, char **argv) {
	htsFile *in, *out;
	int c, ret, exit_code;
	char moder[8];
	//char modew[800];
	char *outfn = "-";

	ftype = FSAM;
	exit_code = 0;
	strcpy(moder, "r");
	while ((c = getopt(argc, argv, "bo:")) >= 0) {
		switch (c) {
			case 'b': strcat(moder, "b"); ftype = FBAM; break;
			case 'o': outfn = optarg; break;
		}
	}
	if (optind + 1 > argc) {
		fprintf(stderr, "Usage: %s [-b] [-o out.sam] <in.bam>|<in.sam> [region]\n", argv[0]);
		fprintf(stderr, "Options:\n");
		fprintf(stderr, "\t-b:\tUse BAM as input if this option is set, otherwise use SAM as input.\n");
		fprintf(stderr, "\t-o:\tPath to the output file. Output to stdout if this option is not set.\n");
		return -1;
	}
	if ((in = hts_open(argv[optind], moder)) == NULL) {
		fprintf(stderr, "Error opening '%s'\n", argv[1]);
		return -3;
	}
	if ((out = hts_open(outfn, "w")) == NULL) {
		fprintf(stderr, "Error opening '%s'\n", argv[2]);
		return -3;
	}
	if ((ret = sam_test_extract(argc, argv, optind, in, out)) != 0) {
		fprintf(stderr, "Error extracting alignment from '%s'\n", argv[optind]);
		exit_code = -5;
	}
	if ((ret = hts_close(out)) < 0) {
		fprintf(stderr, "Error closing output.\n");
		exit_code = -3;
	}
	if ((ret = hts_close(in)) < 0) {
		fprintf(stderr, "Error closing input.\n");
		exit_code = -3;
	}
	return exit_code;
}

其中：
htsFile 结构：储存了sam文件的信息；
sam_hdr_t 结构：储存了sam文件头部的信息；
bam1_t 结构：储存了一条比对结果的信息；
hts_idx_t结构：储存了index文件的信息；
hts_itr_t结构：是一个迭代器，每次返回一条比对结果；

这些结构都有对应的函数对它们进行处理。比如sam_hdr_read(htsFile*)可以读取sam文件的头部信息。其它函数的作用可以参考官网上的头文件。

编译c程序

为什么要将编译单独写一小节呢？我们平时编译c程序时，都很简单，比如：

gcc program.c -o program

这是因为我们平时写的c程序都很简单，用到的库基本都是标准库。但是当用到像htslib这样的第三方库时，编译就会相对复杂一点。以上面的samtest.c为例，理论上编译的命令是：

gcc samtest.c -o samtest -I<htslib_dir> -L<htslib_dir> -l<htslib_name>

如果以我们上面的demo文件夹来说，命令就是

cp samtest.c ~/demo   # if necessary
cd ~/demo
gcc samtest.c -o samtest -Ihtslib -Lhtslib -lhts

编译成功后demo文件夹里就会有一个samtest程序了。

上面-I<dir>选项表示到除了标准头文件目录之外的某个目录下去寻找头文件；-L<dir>选项表示到除了标准库文件目录之外的某个目录下去寻找库文件；-l<libname>选项用来指定库文件。

更多关于编译的知识可以参考《GCC编译器30分钟入门教程》（文末有链接）

运行程序

编译成功后，并不意味着就可以立即运行程序了。因为编译的时候用到了动态链接库，那么就需要确保运行程序的时候系统能找到这些动态链接库，否则程序就运行不了。

在本文的例子中，samtest程序需要加载libhts.so.3这个动态链接库。如果不能正确加载该动态链接库，程序会报错，像这样：

不同的系统具有不同的加载链接库的方法。在Linux系统中，我们可以将动态链接库复制到标准库文件目录下（例如/usr/lib或者/usr/local/lib），或者设置一个合适的环境变量，例如LD_LIBRARY_PATH。当然还有一些其它方法（参考文末链接文章）。

此次我们选择设置环境变量来确保samtest程序能正确加载要用到的动态链接库。在Linux shell中运行：

# export LD_LIBRARY_PATH=<htslib_dir>:<samtools_dir>:$LD_LIBRARY_PATH
export LD_LIBRARY_PATH=~/demo/htslib:~/demo/samtools:$LD_LIBRARY_PATH

再运行samtest程序就成功了：

当然，上面的export命令只能对本次shell会话有效，要想每次开机（登入shell）都能自动运行export命令，可以将其写入~/.bashrc文件中：

# echo 'export LD_LIBRARY_PATH=<htslib_dir>:<samtools_dir>:$LD_LIBRARY_PATH' >> ~/.bashrc
echo 'export LD_LIBRARY_PATH=~/demo/htslib:~/demo/samtools:$LD_LIBRARY_PATH' >> ~/.bashrc

备注：如上文所说，将动态链接库复制到标准库文件目录也可以让程序正确加载动态链接库。对于本文中的例子，在Linux shell中运行：

sudo cp ~/demo/htslib/libhts.so ~/demo/htslib/libhts.so.3 /usr/local/lib

参考

• htslib库源码地址：https://github.com/samtools/htslib
• samtools库源码地址：https://github.com/samtools/samtools
• 使用htslib API的一个示例程序：https://github.com/samtools/htslib/blob/develop/test/test_view.c
• 《GCC编译器30分钟入门教程》：http://c.biancheng.net/gcc/
• 《共享库01-静态库和共享库的使用》：https://blog.csdn.net/Lincoln_cz/article/details/81982908