C++使用htslib库读入和写出bam文件

  有时候我们需要使用C++处理bam文件，比如取出read1或者read2等符合特定条件的序列，根据cigar值对序列指定位置的碱基进行统计或者对序列进行处理并输出等，这时我们可以使用htslib库。htslib可以用来处理SAM, BAM,CRAM 和VCF文件，是samtools、bcftools的核心库。

#include <stdio.h>

#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};

int main(int argc, char **argv)

{

    bam_hdr_t *header;

    bam1_t *aln = bam_init1();

    samFile *in = sam_open(argv[1], "r");

    htsFile *outR1 = hts_open(argv[2], "wb");

    header = sam_hdr_read(in);

    if (sam_hdr_write(outR1, header) < 0) {

        fprintf(stderr, "Error writing output.\n");

        exit(-1);

    }

    uint8_t *seq;

    int32_t lseq;

    uint32_t *cigar;

    char* qname;

    while (sam_read1(in, header, aln) >= 0) {

        if (bam_is_read1(aln)){

            sam_write1(outR1, header, aln);

        }

        else {

            seq = bam_get_seq(aln);

            lseq = aln->core.l_qseq;

            qname  = bam_get_qname(aln);

            printf("%s\n",qname);

            cigar = bam_get_cigar(aln);

            for(int i=0; i < aln->core.n_cigar;++i){

                int icigar = cigar[i];

                printf("%d%d\n",bam_cigar_op(icigar),bam_cigar_oplen(icigar));

            }

            for(int i=0; i < lseq;++i){

               printf("%c", Base[bam_seqi(seq, i)]);

            }

            printf("\n");

        }

    }

    sam_close(in);

    sam_close(outR1);

}

cigar值存储形式

32位int，通过bam_get_cigar获得地址，aln->core.n_cigar返回cigar operation的个数

• 低 4位存储cigar operation；通过函数bam_cigar_op()获得operation
  
  
• 高28位存储cigar值的长度；通过函数,bam_cigar_oplen获得
  

seq存储形式

8位int，4位存储一个碱基，1,2,4,8,15分别代表A、C、G、T、N，高四位存储坐标数较小的碱基，可通过bam_seqi(seq,i)遍历。

参考资料

htslib sam.h文件:https://github.com/samtools/htslib/blob/develop/htslib/sam.h

htslib sam文件格式说明:https://github.com/samtools/hts-specs/blob/master/SAMv1.pdf