高通量数据BAM格式分析包HTSlib使用方法

首先安装依赖的包

wget https://tukaani.org/xz/xz-5.2.4.tar.gz
tar -zxvf xz-5.2.4.tar.gz
cd xz-5.2.4
./configure
make && make install

wget ftp://ftp.invisible-island.net/ncurses/ncurses-6.0.tar.gz
tar -zxvf ncurses-6.0.tar.gz
cd ncurses-6.0
./configure
make && make install

其次

htslib和samtools

wget https://github.com/samtools/htslib/releases/download/1.9/htslib-1.9.tar.bz2
tar -vxjf htslib-1.9.tar.bz2
cd htslib-1.9
make

wget https://github.com/samtools/samtools/releases/download/1.9/samtools-1.9.tar.bz2
tar -vxjf samtools-1.9.tar.bz2
cd samtools-1.9
mak

写好测试代码 parseBam.cpp ， 进行编译

#include <iostream>
#include <stdio.h>
#include <stdexcept>
#include <sam.h>
#include <hts.h>

using namespace std;

int usage() {
    cerr << "./parse_bam example.bam";
    return 0;
}

int parse_bam(int argc, char* argv[]) {
    if(argc < 2) {
        return usage();
    }

    string bam = string(argv[1]);
    string region_ = ".";
    if(argc > 2) {
        region_ = string(argv[2]);
    }
    if(!bam.empty()) {

        htsFile *in = hts_open(bam.c_str(), "r");
        bam_hdr_t *header = sam_hdr_read(in);
        hts_itr_t *iter = NULL;

        iter  = sam_itr_querys(idx, header, region_.c_str());

        bam1_t *aln = bam_init1();
        while(sam_itr_next(in, iter, aln) >= 0) {
            printf ("Read %s\tChr: %s" , bam_get_qname(aln) , header->target_name[aln->core.tid]);
            printf ("\tPos: %d" , aln->core.pos);
            string seq, qual;
            uint8_t *quali = bam_get_qual(aln);
            uint8_t *seqi = bam_get_seq(aln);
            for (int i = 0; i < aln->core.l_qseq; i++) {
                seq += seq_nt16_str[bam_seqi(seqi, i)];
                qual += 33 + quali[i];
            }

            printf("\tSeq: %s\tQual: %s" , seq.c_str(),  qual.c_str() );
            printf("\n") ;
        }

        hts_itr_destroy(iter);
        bam_destroy1(aln);
        bam_hdr_destroy(header);
        
        sam_close(in);
    }
    return 0;
}

int main(int argc, char* argv[]) {
    try {
        parse_bam(argc, argv);
    } catch (const runtime_error& e) {
        cerr << e.what();
    }
}

假设 htslib和samtools 都安装在root路径下，按照如下命令编译  

g++ parseBam.cpp -o parseBAM -I/root/htslib/htslib -L/root/samtools -lhts -std=c++11

生成可执行文件parseBam， 第一个参数是bam文件，第二个参数是chrN:start-end，  最后执行程序如下:

 parseBam example.bam chr9:79327223-79387223