htslib安装和使用

htslib安装和使用

本文作者：Sunny-King

发布时间：2022-11-30 22:09:47 星期三

本文链接：https://www.cnblogs.com/Sunny-King/p/Bioinformatics-htslib.html

最近工作中遇到了基于htslib实现所需功能，发现官方的文档描述很模糊，只能去看代码及注释。并且网上给出的例子大多模糊不清，并且有Bug，不能准确编译执行。记录下自己安装和编译过程，以及基于网上例子的注释。

一、下载安装

1、从GitHub下载源文件

pwd
/home/test/
git clone https://github.com/samtools/htslib.git
cd htslib

2、安装

# 指定安装路径
./configure --prefix=/home/test/htslib
make
make install

会在/home/test/htslib路径下生成新的文件夹{share、lib、include、bin}等。

二、添加环境变量

1、添加动态链接库

需要将对应路径添加到环境中，或者拷贝至动态链接库的路径下。

export LD_LIBRARY_PATH=/home/test/htslib:$LD_LIBRARY_PATH

2、编译执行

cd /home/test/
g++ samtest.cpp -o samtest -Ihtslib -Lhtslib -lhts

./samtest test.bam out.bam > test.log

该段代码的功能是读取原始bam文件，将reads1输出到输出的bam中，将reads2的qname、cigar、seq输出到终端

三、具体代码和注释

C++代码及注释

#include <iostream>
#include <stdlib.h>
#include "htslib/sam.h"

using namespace std;

#define bam_is_read1(b) (((b)->core.flag&BAM_FREAD1) != 0)

uint8_t Base[16] = {0,65,67,0,71,0,0,0,84,0,0,0,0,0,0,78}; //1,2,4,8,15分别代表A、C、G、T、N，

int main(int argc, char **argv)
{
    bam_hdr_t *header;
    bam1_t *aln = bam_init1();

    samFile *in = sam_open(argv[1], "r");  //读方式打开bam文件
    htsFile *outR1 = hts_open(argv[2], "wb"); //写方式打开二进制文件
    header = sam_hdr_read(in); //获得bam的header
    //header写入新文件
    if (sam_hdr_write(outR1, header) < 0) {
        cerr << "Error writing output." << endl;
        exit(-1);
    }

    uint8_t *seq;
    int32_t lseq;
    uint32_t *cigar;
    char* qname;
    while (sam_read1(in, header, aln) >= 0) {
        /*
        >= 0    成功读取一条新记录
        = -1    读取到结束标记
        < -1    报错
        */
        //read1写入文件
        if (bam_is_read1(aln)){
            int sam_return = sam_write1(outR1, header, aln); //写入文件
        }
        //read2输出到终端
        else {
            seq = bam_get_seq(aln); //获取序列
            lseq = aln->core.l_qseq; 
            qname  = bam_get_qname(aln); //获取qname
            cout << qname << endl; //输出qname
            cigar = bam_get_cigar(aln); //获取cigar
            //逐段输出CIGAR标签
            cout << dec << aln->core.n_cigar << endl;
            for(int i=0; i < aln->core.n_cigar;++i){
                int icigar = cigar[i];
                cout << dec << icigar << "\t";
                cout << dec << bam_cigar_op(icigar) << "\t" << dec << bam_cigar_oplen(icigar) << endl;
            }
            for(int i=0; i < lseq;++i){
                cout << Base[bam_seqi(seq, i)];
            }
            cout << endl;
        }
    }
    //回收资源
    bam_destroy1(aln);
    bam_hdr_destroy(header);
    //关闭文件
    sam_close(in);
    sam_close(outR1);
}